Evolution of cold streams in hot gaseous halos

In the prevailing model of galaxy formation and evolution, the process of gas accretion onto central galaxies undergoes a transition from cold-dominated to hot-dominated modes. This shift occurs when the mass of the parent dark matter halos exceeds a critical threshold known as $M_{shock}$. Moreover, cold gas usually flows onto central galaxies through filamentary structures, currently referred to as cold streams. However, the evolution of cold streams in halos with masses around $M_{shock}$, particularly how they are disrupted, remains unclear. To address this issue, we conduct a set of idealised hydrodynamic simulations. Our simulations show that (1) for a gas metallicity $Z=0.001-0.1Z_{\odot}$, cold stream with an inflow rate $\sim 3\, \rm{M_{\odot}}/yr$ per each can persist and effectively transport cold and cool gas to the central region ($<0.2$ virial radius) in halos with mass $10^{12}\, \rm{M_{\odot}}$, but is disrupted at a radius around $0.2$ virial radius due to compression heating for halos with mass $3 \times 10^{12}\, \rm{M_{\odot}}$. (2) At $z\sim 2$, the maximum halo mass that capable of hosting and sustaining cold streams $M_{stream}$ is between $1\times 10^{12} \rm{M_{\odot}}$ and $1.5\times 10^{12}\rm{M_{\odot}}$ for gas metallicity $Z=0.001Z_{\odot}$, while for a higher gas metallicity $Z=0.1Z_{\odot}$, this value increases to $\sim 1.5\times 10^{12}\rm{M_{\odot}}$. (3) The evolution and ultimate fate of cold streams are determined primarily by the rivalry between radiative cooling and compression. Stronger heating due to compression in halos more massive than $M_{stream}$ can surpass cooling and heat the gas in cold streams to the hot ($\geq 10^6\,$ K) phase.


INTRODUCTION
The acquisition of gas is a pivotal facet of galaxy formation and evolution, a process whose understanding has improved significantly over the past five decades.This knowledge has unfolded in tandem with the advancements in ΛCDM cosmology.Once dark matter halos are formed, they accrete gas from their cosmic surroundings.Following the epoch of reionization, completed at redshift  ∼ 5.5, the gas that was initially acquired typically has a temperature of approximately 10 4 K, the typical temperature of the intergalactic medium heated by the cosmic ultraviolet background (UVB).In the conventional scenario, the gas accreting onto halos will be heated to the virial temperature by shocks produced by the gravitational collapsing, then subsequently loses pressure support because of radiative cooling, and eventually settles into a rotation-supported disc, where stars were born (e.g., Rees & Ostriker 1977;Silk 1977;White & Rees 1978;Fall & Efstathiou 1980;Bertschinger 1985;White & Frenk 1991).This picture was partly supported by hydrodynamical simulations conducted in the early 2000s, which had demonstrated that a substantial fraction of cold gas was indeed shock heated following its accretion (e.g., Benson et al. 2001;Helly et al. 2003).
However, in a theoretical study on the collapse of protogalactic cloud, Binney (1977) concluded that specific conditions related to ★ E-mail: zhuwshan5@mail.sysu.edu.cn(WSZ) the infalling velocity and the surface density of the descending gas are required to maintain adiabatic collapse shocks, which make the post-shock gas achieve virial temperature.Since the early 1990s, numerous simulations of galaxy formation with gas cooling have suggested that a substantial portion of the gas in galaxies has never been heated above 10 5 K (e.g.Katz & Gunn 1991;Kay et al. 2000;Fardal et al. 2001).Later, Birnboim & Dekel (2003) used an analytical model and one-dimensional simulations, revealing that virial shocks fail to develop in halos with masses below a critical threshold, approximately 10 12 M ⊙ at a redshift of  = 2 and are usually called as M shock , due to radiative cooling.Consequently, gas entering these less massive halos remains cold (10 4 − 10 5 K) along their journeys to the central region, usually called the cold accretion mode, and can drive a high star formation rate at high redshifts.Meanwhile, gas falling into halos more massive than the critical mass will be shock-heated to near the virial temperature.This scenario has been confirmed and further improved by several theoretical and simulation works in the past two decades.The transitions from cold-mode dominated accretion in low-mass halos and galaxies to hot-mode dominated accretion in massive systems have been explored in great detail (e.g.Kereš et al. 2005;Dekel & Birnboim 2006;Ocvirk et al. 2008;Dekel et al. 2009;Kereš et al. 2009;van de Voort et al. 2011;Nelson et al. 2013;Wright et al. 2021).The critical transition halo mass, M shock , exhibits modest variation in different studies, spanning a range from 10 11.4 to 10 12.0  ⊙ at a redshift of  ∼ 2. Remarkably, simulations have illustrated that cold accretion predominantly occurs along filamentary structures and are often referred to as 'cold streams' (e.g., Kereš et al. 2005;Dekel & Birnboim 2006;Dekel et al. 2009;Nelson et al. 2013).Furthermore, these cold streams can penetrate through shocks within halos with masses above M shock at  ≳ 1 − 2, and the maximum halo mass that can host cold streams, denoted as M stream , increases with redshift (Dekel & Birnboim 2006).
Since the late 2000s, great efforts have been made to detect cold streams.Theoretical studies based on simulations suggest that identifying observationally cold streams through metal absorption lines is relatively difficult while the approach of Lyman-alpha absorption and emission is more promising (e.g.Dĳkstra & Loeb 2009;Goerdt et al. 2010;Faucher-Giguère & Kereš 2011;Fumagalli et al. 2011;Kimm et al. 2011;Goerdt et al. 2012).However, despite several observed features that could potentially be attributed to cold accretion flows, there is still a lack of conclusive observational evidence for the existence of cold streams (e.g., Crighton et al. 2013;Martin et al. 2014;Daddi et al. 2021;Fu et al. 2021).More recently, Daddi et al. (2022b) reports the finding of observational evidence for the transition from cold to hot accretion modes by measuring the ratio of Lyman-alpha luminosity from groups and clusters and the expected baryonic accretion rate at  ∼ 2. Nevertheless, interpreting the observed features faces notable challenges and uncertainties.For instance, the Lyman-alpha emission may be partly powered by outflows and photoionization stemming from star formation and active galactic nuclei.On the other hand, the kinetics, clumpiness, and disrupted radius of cold streams remain somewhat unclear within the theoretical framework.
To consolidate the scenario of cold accretion through streams and to have a more robust interpretation of current and upcoming observations, it is imperative to dive into a deeper understanding of the evolution and fate of cold streams in massive halos.More specifically, several important open questions need to be addressed.How these cold streams are disrupted and lead to the closure of the cold accretion channel when the halo mass becomes comparable to M stream ?How long can these cold streams survive against various perturbations and instabilities, including tidal forces due to interaction, Kelvin-Helmholtz, Rayleigh-Taylor and thermal instabilities (e.g.Wang et al. 2014;Mandelker et al. 2016;Berlok & Pfrommer 2019;Mandelker et al. 2019Mandelker et al. , 2020;;Ledos et al. 2023)?What impact does feedback from central galaxies and outflows have on the cold streams?
In this work, we utilize idealized three-dimensional hydrodynamic simulations to investigate the evolution of cold streams in isolated intermediate mass halos, aiming to gain insight into the mechanisms and processes of cold stream disruption in these halos.The layout of this paper is as follows: Section 2 describes the setup of our models for cold streams, hot gas halos, the simulations, and the typical timescales of relevant physical mechanisms.Detailed results from the simulations are presented in Section 3. In Section 4, we discuss the impact of the conduction process and the caveats of our study, as well as some factors that can be improved in future work.Finally, we summarise our findings in Section 5.

Models of hot gaseous halos and cold streams
We adopt a simplified model of cold streams in the hot gaseous halo, which allows us to focus on the evolution of cold streams from the outside of the halo to the boundary of the central galaxy.We firstly set up a hot gaseous halo that is assumed to be in hydrostatic equilibrium under the gravitational potential of a dark matter halo with a given mass around the critical mass M shock , which indicates the gas accretion onto central galaxies transits from cold-mode dominated to hot-mode dominated at redshift  = 2.The properties of the dark matter halo are assumed to be non-evolving during the simulation.Subsequently, we introduce a continuous injection of cold gas into the hot gaseous halo from a location beyond the virial radius.In the gravitational well, the cold gas flows towards the central region of the gaseous halo and experiences radiative cooling and compression along its streamline.Figure 1 provides a schematic diagram of our model.We track the evolution of cold streams in halos of different masses.The specifics of our models are outlined in the following subsections.
Before delving into the physics of our numerical models, it is essential to emphasise that our simplified model intentionally excludes the central galaxy and the associated physics of its interstellar medium (ISM), stellar component, and central supermassive black hole.Our primary objective is to concentrate on the evolutionary dynamics of the cold stream originating from external regions beyond the halo and extending downward to the boundary of the central galaxy.The latter has a radius around 0.10-0.15times the virial radius of the halo.Consequently, we acknowledge that the results of our study possess limited predictive capacity in the inner region, given the absence of a detailed and realistic model in that particular domain.

Hot Gaseous Halos
In our study, the background gravitational potential is determined by a dark matter halo that follows the Navarro-Frenk-White (NFW) Table 1.Parameters for the hot gaseous halos and associated dark matter halos with different masses.The gaseous halos are assumed to be isothermal with a temperature  ℎ .The scale radius   , concentration c, and scale density   are derived using the python package COLOSSUS with Planck18 cosmological parameters and ishiyama20 model at z=2.The gas density at the halo center   is derived according to the cosmic baryon fraction (see Section 2.1 for more details).profile (Navarro et al. 1996), where  ℎ and   are the characteristic overdensity parameter and the critical density of the universe, respectively;   denotes the scale density, and   is the scale radius.Consequently, the background gravitational potential can be expressed as, In this work, we consider dark matter halos with masses of 1.0, 1.5, and 3.0 × 10 12 M ⊙ at redshift  = 2.The values of   ,   and c are calculated using the Python package COLOSSUS (Diemer 2018) with the Planck18 cosmological parameters (Planck Collaboration et al. 2020), and the Ishiyama20 model (Ishiyama et al. 2021).Under the static potential of the dark matter halos, the density profile of the associated gaseous halo under the hydrostatic equilibrium condition would be (Mo et al. 2010) where   = (0) is the gas density at the halo center, and b is given by where   is the mass of a proton,   is the Boltzmann constant, and  ℎ is the temperature of halo gas.We determine the value of   by assuming that the ratio of the gas mass to the dark matter mass within the virial radius is equal to the cosmological mean Ω  /(Ω  − Ω  ).The temperature of hot halo gas  ℎ is set to around the virial temperature.The parameters relevant to the gaseous halo and their specific values are listed in Table.1.To reduce the possible boundary effect in our simulations, we will extend the gaseous halos to 1.5 times the virial radius   , and the side length of the simulation box will be set to twice this enlarged scale.

Cold streams
To mimic the cold streams, cold gas is injected continuously into the hot gaseous halo through a circular cross-section at the boundary region of the simulation box.Following the consideration in Wang et al. (2014), we set the radius of the cross-section to be 15 kpc, which will have a covering factor close to previous studies (e.g., Goerdt et al. 2012;Fumagalli et al. 2013).Under the influence of the gravitational field, the cold gas develops a conical stream entering the halo's central region.The injected cold gas is endowed with an initial velocity   of 100 km/s along the x-direction at the boundary of the simulation box, which is close to the typical inflow velocity of cold flow at  = 1.5   for halos with  ℎ ∼ 10 12 M ⊙ , that is, ∼ 0.7 times of   = √︁   ℎ /  , at redshift  = 2 in cosmological simulation (Goerdt & Ceverino 2015).The temperature of the cold streams is 3.3 × 10 4 K, which is moderately higher than the mean temperature of the IGM at  = 2 that is heated by the cosmic UV background (see, e.g.Gaikwad et al. 2021).However, it is comparable to the typical temperature of cold streams in cosmological hydrodynamic simulations (e.g., Nelson et al. 2013).The mean density of cold gas is set to 3.0 × 10 −27 g • cm −3 , which is an order of magnitude lower than that in Wang et al. (2014).The inflow rate of a cold stream is 30 M ⊙ /yr in Wang et al. (2014).An injection rate of 3-5 M ⊙ /yr should be more reasonable, as the cold accretion rate to halo is about 10-20 M ⊙ /yr at  ∼ 2 for a dark matter halo with mass  ℎ ∼ 10 12 M ⊙ (e.g.van de Voort et al. 2011;Zhu et al. 2022), and usually there are 3-4 streams in a single halo.For simplicity, the injection rate of each cold stream is fixed and does not depend on the halo mass.We will discuss the potential limitation of this assumption in Section 4.
Note that previous simulations show that the gas in the cold streams is inhomogenous and partly clumpy (e.g., Dekel & Birnboim 2006;Ocvirk et al. 2008;Nelson et al. 2013).Moreover, this inhomogeneity rises outside of the halos.Actually, the density of the intergalactic medium (IGM) indeed fluctuates to a certain extent.In the mildly non-linear regime, the probability distribution function (PDF) of the density of the IGM can be modelled by a log-normal random field (Bi & Davidsen 1997).In the highly non-linear regime, the IGM exhibits turbulent features (Zhu et al. 2013;Zhu & Feng 2017), and shows an extended density distribution.Taking into account these factors, we further assume that the gas density of the injected cold gas follows a log-normal distribution of supersonic turbulence (e.g., Konstandin et al. (2016)).Specifically, the mean of the logarithmic density is denoted as s = () and is given by : in which   ,   are the mean values of the natural logarithm of the mass density and the mean of the original mass density, respectively.The standard deviation of   is related to the Mach number as follows, where   can be interpreted as the ratio of the compressive Mach number to the total root-mean-square Mach number of the flow, i.e.
The value of   depends on the driving force of turbulence.For purely solenoidal forcing,   = 1/3, and for purely compressive forcing,   = 1.In this work, we adopt   = 0.5, since both solenoidal and compressive modes are presented in the turbulent IGM (Zhu et al. 2013).

Simulations
To investigate the evolution of cold streams in halos, we run hydrodynamic simulations using the ATHENA++ code (Stone et al. 2020).The simulations are performed in a 3D Cartesian geometry with a uniform cubic grid.The simulation box has a side length of 3  , and the halo center is located at the origin (0,0,0).The cold stream is injected at a specific position (-1.5 , 0.3  , 0), taking into account that cold streams could have a non-negligible angular moment.
Several relevant physical processes may affect the evolution of cold streams, including cooling, conduction, and the strength of the gravitational potential.These factors were considered in our simulations to identify the mechanisms that lead to the disruption of cold streams.In addition, the impacts of the cold streams' density perturbation and the resolution of the simulation were also explored.We have run a set of simulations to investigate the roles of each factor by turning on/off different processes or modifying relevant parameters in different simulations.Hereafter, the simulations are named 'X(X.Y)e12-ĲK-abc(-HZ)', where 'X(X.Y)e12' represents the mass of host dark matter halo, 'ĲK' is the number of grids along each box side, the first character in 'abc' indicates whether there is a perturbation in the density of cold streams or not, denoted as 'p' or 'n' respectively.The second (third) character in 'abc' stands for whether cooling (conduction) is turned on or off, denoted by 'c', or 'n' correspondingly.'HZ' indicate a gas metallicity of  = 0.1 ⊙ , while other simulations adopt a metallicity of  = 0.001 ⊙ .The full list of our simulations can be found in Table 3.More details about the simulations are described below.

Radiative Cooling and thermal conduction:
In the simulations, the effects of radiative cooling and thermal conduction have been considered, which allows us to analyse how they may affect the fate of cold streams.We have implemented these two effects by solving the following hydrodynamic equations with the Athena++ code (Stone et al. 2020): where  and v stand for the density and velocity respectively,  represents the energy density,  = ( − 1) denotes the pressure,  signifies the internal energy density, and  = 5/3; Λ  =  2 Λ   stands for the net cooling rate per unit volume, where  is the number density and Λ   is the normalised cooling function and is obtained from the cooling table given by Sutherland & Dopita (1993) in this work;q is the conduction flux.Cooling and thermal conduction are handled by updating the energy density term  after solving hydrodynamic equations without source terms at each time step.That is, Λ  and ∇ • q are treated as source terms in the equation of energy density.
The cooling rate and conduction flux in each cell at each time step are calculated using the methods described below.
The cooling rate of the gas in a grid cell depends on its temperature, hydrogen number density, and metallicity.In this work, we assume that the cold streams and the gaseous halo have the same metallicity for simplicity.Yet, to explore the impact of metallicity, we adopt two different metallicities within our simulations, i.e.,  = 0.001 ⊙ , and  = 0.1 ⊙ , respectively.The former value is adopted because the metallicity of the diffuse IGM at  ∼ 2 is estimated to be 10 −2 − 10 −3  ⊙ (e.g., Aguirre et al. 2008;Wiersma et al. 2009;D'Odorico et al. 2016).The latter value is somewhat lower than the metallicity of the ISM at ∼ 2, i.e., 0.16 − 1.0 ⊙ , (e.g., Erb et al. 2010;Steidel et al. 2014;Torrey et al. 2019;Sanders et al. 2023).Thus, the metallicity of gas in the cold stream at  ∼ 2 is very likely between  = 0.001 ⊙ and  = 0.1 ⊙ .For each simulation with a given metallicity, we use the cooling table of the corresponding metallicity provided by Sutherland & Dopita (1993) to find the normalised cooling rate, Λ   , of gas in different cells by interpolation, where the temperature of the gas is inferred from the internal energy and density.The normalised cooling rate, Λ   , as a function of gas temperature for these two metallicities, is shown in Fig. 2. The vertical dotted lines mark the peak of the cooling rates at around 10 5 K.If the gas temperature approaches the corresponding value for a particular metallicity, the cooling rate is boosted, which could cause the gas to cool down rapidly if the increment of pressure due to compressional heating is insufficient.
The cooling function employed in our study takes into account the cooling of gas with a temperature higher than 10 4 K.That is, cooling below temperature 10 4 K is totally neglected.We expect that such an implementation has a minor effect on our main results about the evolution of the cold stream from somewhere outside of the halo inward to the boundary of the central galaxy because of the following considerations.First, the normalised cooling rate reduces significantly from around 10 −22.5 to 10 −24.5 erg s −1 cm 3 for a gas number density of 1/cm 3 , when the temperature drops from slightly above 10 4 K to several 10 3 K (e.g.Maio et al. 2007;Smith et al. 2017).Second, molecular cooling is the primary cooling mechanism of gas below 10 4 K.It usually requires a total hydrogen number density of   > 10/cm 3 to form molecules effectively (e.g., Glover et al. 2010;Wakelam et al. 2017), which is much higher than the typical number density of gas in the cold stream, i.e., ∼ 10 −3 − 10 −1 /cm 3 .
On the other hand, radiative heating from the cosmic UVB and central galaxy is not included in our simulation.The initial temperature adopted for gas in the cold stream has accounted for the heating of IGM by the UVB.Moreover, self-shielding would partially attenuate the heating of UVB.In addition, neglecting the cooling below 10 4 K can partially compensate for the absence of UVB.Yet, the omission of radiative heating from the central galaxy, including UV and X-rays from accretion to the SMBH, would lead to over-cooling for the hot gas in the central region.Mandelker et al. (2020) shut off cooling for gas with T >∼ 2 × 10 5 K to prevent the hot gas from over-cooling after long periods.To compensate for the probable overcooling of the hot gas in the central region in our simulations, we manually reduce the cooling rate, Λ net , of the gas at each simulation step such that the cooling time is not shorter than 1Gyr for gas cells with a temperature of T > 5 × 10 5 K and a number density of  > 1 × 10 −2 cm −3 .In practice, this remedy works only for gas within the inner 0.2   of halos in our simulations.
In addition to radiative cooling, thermal conduction, i.e., the transfer of energy due to temperature gradients, has also been found to play an important role in the survival of cold clouds within ambient gas (e.g., Armillotta et al. 2016;Armillotta et al. 2017).Thermal conduction may affect the cooling of hot gas and the diffusion of cold streams, especially in the diffuse and intermediate stream regimes (Ledos et al. 2023).The second term on the RHS of equation 8 represents the thermal conduction process.According to the Spitzer formula (Spitzer.L 1962), the heat conduction flux is where ▽ is the temperature gradient, and the heat conduction coefficient is: where Ψ is the Coulomb logarithm: Following the consideration and treatment in Armillotta et al. (2017), the heat conduction flux can be modified to: where  indicates the absolute ratio of the classical heat flux to the saturated heat flux.In this work, we use  = 49, because we find that such a value can reproduce the results well in Armillotta et al. (2017).
Before moving on to more details of the simulations, we first estimate the relative strength between cooling and thermal conduction.In this regard, comparing the Field length (Begelman & McKee 1990;Armillotta et al. 2016), defined by the maximum length scale on which heat energy transport is effective, with the stream radius, can serve as a valuable indicator.The Field length read as, where  = 0.1 1+ represents the heat conduction coefficient introduced above,   ,   is the number density and temperature of the cold stream, respectively, and  ℎ is the hot halo temperature.The function Λ(  ) represents the cooling rate with temperature   .In our models, the Field length is 10 (20) pc for halo mass 1.0× 10 12 M ⊙ (3.0 × 10 12 M ⊙ ) and metallicity Z=0.001 ⊙ , which is much smaller than the typical radius of a stream, i.e., 2-15 kpc.Even if the value of  is reduced to 2-5, the Field length is still far below the radius of the stream.Moreover, as we will see, the typical cooling timescale is much shorter than the thermal conduction timescale.On the basis of these results, we expect that thermal conduction is unlikely to affect our primary conclusions regarding the evolution of cold streams in the presence of radiative cooling.
Nevertheless, to provide a comprehensive evaluation, we will inspect the effect of thermal conduction on the evolution of cold streams by turning the thermal conduction on and off in different simulations.

Kelvin Helmholtz Instability (KHI) and comparison of multi-time scales:
Some previous studies suggest that Kelvin-Helmholtz instability(KHI) can play a significant role in the evolution of the stream (e.g., Mandelker et al. 2016Mandelker et al. , 2019Mandelker et al. , 2020)).For a better understanding of our simulations, it is beneficial to evaluate the probable overall effects of KHI.If radiative cooling is included, a cold stream with a radius greater than the critical radius   (Aung et al. 2019;Mandelker et al. 2020) would grow in mass, instead of being disrupted by the KHI.The critical stream radius   is given by Mandelker et al. (2020), where  0.1 = /0.1 = 0.21 × (0.8 (−M 2  ) + 0.2)/0.1,M  =   /(  +  ℎ ) is the total Mach number,   and  ℎ denote the sound speed of the stream and the background halo gas, respectively. 100 = /100, where  =   / ℎ represents the ratio of the density of gas in the stream to the density of the background halo gas.M ℎ =   / ℎ represents the Mach number of background halo gas.Furthermore,  4 =   /(10 4 K) indicates the temperature of the stream in units of 10 4 K. Additionally,  ,0.01 =   /0.01 cm −3 and Λ ,−22.5 = Λ   /10 −22.5 erg s −1 cm 3 .The subscript 'mix' indicates the turbulent mixing layer (Begelman & Fabian 1990) formed by the turbulent mixing driven by KHI at the interface of two-phase gas: cold, dense gas, and hot, diffuse gas.Recent studies (Gronke & Oh 2018;Mandelker et al. 2020;Chen et al. 2023) have revealed that cooling in these layers can significantly impact the fate of the cold streams.In our model, the critical stream radius  , is typically 1.54 (7.69 kpc) at r=1.0   for a halo mass of 1.0 × 10 12 M ⊙ (3.0 × 10 12 M ⊙ ) and a metallicity Z=0.001 ⊙ .These results are detailed in Table .4and are smaller than the original radius of the cold stream.As shown in Mandelker et al. (2020),  , would decrease with r in a trend similar to   .Therefore, the condition  , <   would hold throughout the evolution of a cold stream.According to Mandelker et al. (2020), it is reasonable to expect that the KHI would not dissolve the stream in our model.
Alternatively, the collective impact of various mechanisms influencing the evolution of a cold stream such as Kelvin-Helmholtz Instability (KHI), cooling, thermal conduction, and compressional heating can be estimated by considering their respective typical timescales.For a 3-dimensional cylindrical cold stream, Mandelker et al. (2016) and Mandelker et al. (2019) show that the characteristic time for the growth of perturbations with wavelengths  ∼ 1 − 10   due to KHI in the linear regime,  KH , is related to the sound crossing timescale,   , as: It suggests a  KH of around 0.9-1.8Gyr for the cold stream in our model, i.e., about 1/4 to 1/2 of the simulation time we set.
The thermal conduction timescale for the stream can be derived from the analysis in Ledos et al. (2023), It gives a thermal conduction timescale of up to 10 4 Gyr for our model, which is significantly larger than other timescales.
The cooling timescale reads as, We have evaluated both the cooling time of the gas in the hot halo and in the mixing layer at one and a half virial radius for both halos with mass 1.0 × 10 12 M ⊙ and 3.0 × 10 12 M ⊙ , respectively.Both two gas metallicity Z=0.001 ⊙ and Z=0.1 ⊙ are considered.As shown in Mandelker et al. 2020, the mean number density and temperature of the mixing layer to be    ≃ (  ×  ℎ ) 1/2 and    ≃ (  ×  ℎ ) 1/2 respectively.The values of estimated cooling timescales are listed in Table 4.For the hot halo gas, the cooling time is much longer than the Hubble time at  =   , and drops to several Gyr at  = 0.5   .In contrast, the cooling of gas in the mixing layer is much more effective.The corresponding cooling time scale is about 1-3 Gyr at  =   for Z=0.001 ⊙ , and decreases to less than 0.3 Gyr for Z=0.1 ⊙ .The cooling time of gas in the mixing layer decreases sharply as r decreases.Dekel & Birnboim (2006) demonstrated that if the cooling rate exceeds the compression rate required to restore pressure in the postshock gas, the post-shock gas will become unstable against radial gravitational contraction and cannot sustain the shock.Conversely, when the cooling rate is slower, the pressure acquired through compression can counterbalance the loss from radiative cooling, allowing the post-shock gas to remain stable against global gravitational collapse and thus support the shock.Specifically, the timescale of compressional heating is given by where Γ ≡ 3+2  (3−4) .The compressional heating timescales in all of our simulations are approximately 0.5 Gyr, which is shorter than the cooling timescale of hot halo gas in most regions.For the gas in the mixing layer, the cooling timescales are larger or comparable to the compression timescale for gas metallicity Z=0.001 ⊙ at  >   but are shorter for Z=0.1 ⊙ .
Given the timescales estimated above, compressional heating is expected to exert a more significant influence than Kelvin-Helmholtz instability (KHI) on the evolution of the cold stream in all the cases considered and dominate the evolution for the case Z=0.001 ⊙ and  ℎ = 3.0 × 10 12 M ⊙ .However, cooling is expected to become more important if the metallicity is increased or the halo is less massive and will be a dominant player in case Z =0.1 ⊙ and  ℎ = 1.0× 10 12 M ⊙ .Thermal conduction would play a limited role in our models.

Tracer of cold streams
To trace the cold gas injected from the boundary, we use the passive scalar module in Athena++.Based on the tracer field, we define a quantity called '   ', which is the ratio of the density of gas that was injected as a cold stream, denoted as    , to the total gas density   , in each grid cell.More specifically,    at a grid (i,j k) reads as, where , ,  are the ranks of the grid along three dimensions.If    equals one at a particular grid, this grid cell is entirely occupied by gas injected as part of a cold stream at the boundary.When    equals 0, the cold stream has not yet reached the corresponding grid cell.

Resolution
To inspect the impact of resolution, we have run two simulations that track the evolution of a cold stream in a halo with mass 10 12 M ⊙ , using two different resolutions, 512 3 and 256 3 grids, respectively.Except for the resolution, all the other settings are the same in these two simulations, '1e12-256-pcn' and '1e12-512-pcn'.A gas metallicity of  = 0.001 ⊙ is adopted.We find that the overall results are similar between these two resolutions.However, the accumulation of cold gas in the inner region ( < 0.2  ) starts earlier in the higher resolution run (see Appendix B for more details).However, the central galaxy is omitted from our current model.A realistic model of the central galaxy is required to obtain more reliable results on gas evolution in the central region, in addition to higher resolution.
As claimed in Nelson et al. (2020), a resolution finer than Δ < 100 pc is necessary to fully resolve the small-scale structure of cold gas in the CGM.Consequently, a simulation with a box size of 300 kpc requires nearly 3000 grids in each dimension for a uniform grid, which is extremely costly for a 3D hydrodynamic simulation.Adaptive mesh refinement would be helpful in resolving this issue.However, designing the refinement strategy for the evolution of the Table 4.This table displays cooling timescales and typical critical radius of gas in the hot gas and in the mixing layer.The first and second columns represent the gas metallicity and the halo mass.The third indicates the distance from the center of the halo.The fourth, fifth, and sixth columns are the temperature, number density, and cooling timescales of hot halo gas, respectively.The seventh to ninth columns are the same as the fourth to sixth but for the gas in the mixing layer.The last column is the typical critical radius.cold stream in the CGM is not straightforward.On the other hand, the timescales described in the previous section indicate that the dominant process in our models would be cooling and compressional heating.Consequently, to optimize computational resources, we have performed additional simulations utilizing a default configuration of 256 3 fixed grids, i.e., reaching a resolution around 1.2 kpc.The default resolution we adopted may have somewhat underestimated the impact of instabilities on scales of ∼ 100 pc; however, it is expected to adequately capture the overall feature of stream evolution at radial distances  > 0.2  , concluded from a comparative analysis between the simulations '1e12-256-pcn' and '1e12-512-pcn'.

Boundary condition
In our model, the gravitational field is assumed to be static.We found that the outflow boundary conditions built into Athena++ are not suitable for our simulation, as they can cause instability in the ghost zone at the boundaries.To address this issue, we have implemented user-defined boundary conditions.Specifically, we design the boundary condition in those ghost zones adjacent to the stream to allow cold streams to continuously flow into the halo, while the properties of the gas in the other ghost zones are set to fix to the profile of the gaseous halo as described in the section 2.1.Under this setting, the inflow and outflow at the boundary adjacent to hot halo gas is somewhat suppressed.Nevertheless, the hot halo can sustain a steady state throughout the simulation.

Simulations without inflow cold streams
The density profile of the gaseous halo given by equation 3 is derived without considering cooling.Before we explore the evolution of the cold stream in the halo, it is necessary to check the evolution of the hot gaseous halo under the influence of cooling and the gravitational potential.To this end, we have performed four test simulations without injected cold streams for halo mass 1.0 and 3.0 × 10 12 M ⊙ , with metallicity 0.001 ⊙ and 0.1 ⊙ respectively.The temperature and density profiles as a function of radius at time 0, 2.0, and 4.0 Gyr for the simulations with a halo mass  ℎ = 1.0 × 10 12 M ⊙ are depicted in Figure 3 for both metallicities.After 4 Gyr, the temperature and number density of the halo have experienced notable changes.The change in the inner region is more violent, showing an evident temperature decrease and an increment of the density gradient.On the other hand, the gas temperature is moderately enhanced while the density decreases in the outer region.The results of the other two simulations with  ℎ = 3.0 × 10 12 M ⊙ are similar (Figures are not shown, for the sake of conciseness).These changes can be explained by the combined effect of cooling and compressional heating.The cooling time of the hot halo gas, roughly on the order of the Hubble time around the virial radius, results in relatively subtle changes in the outer region.As the radius decreases, this timescale decreases to 1Gyr at around 0.2   and further reduces below 1 Gyr in the innermost region.However, our recipe that avoids over-cooling ensures that the cooling time for gas with  > 5 × 10 5 K and number density  > 1 × 10 −2 cm −3 is not less than 1 Gyr.In all, the cooling time scale of the central region is shorter than the simulation run time of 4 Gyr.This is why the temperature dips and the number density is enhanced in the inner region, as shown Figure 3.Meanwhile, the gas in the outer region slowly moves inward under the gravitational field as the temperature drops in the inner region.Consequently, the temperature (density) of gas in the outer region gradually increases (decreases).Moreover, higher metallicity leads to more significant changes because of more efficient cooling.However, it is important to note that at t=4.0 Gyr, the average temperature of the gas remains above ∼ 10 6 K throughout, from the halo center to the region beyond the virial radius.The analysis based on the four test simulations affirms that if any cold gas appears within the halo's central region in the simulations incorporating injected cold streams, cold streams should be the game changer.

Gas phases
As the cold gas descends and interacts with the hot halo gas, undergoing processes such as compression, radiative cooling, and turbulent mixing, the gaseous medium manifests into multi-phases.To describe the evolution of different phases in the following context, we classify gaseous mediums into three different categories according to their temperature: (i) Cold and cool: gas with a temperature T < 10 5 K. (ii) Warm: gas with 10 5 ≤ T < 10 6 K. (iii) Hot: gas with T ≥ 10 6 K.

RESULTS
In this section, we present our simulation results starting from the ones with a gas metallicity  = 0.1 ⊙ , where the cooling of gas in the cold stream is more efficient and hence has a higher chance to survive in the hot gaseous halos.Then, we inspect the impact of gas metallicity using other simulations.

Overall evolution
We first inspect the overall evolution of inhomogeneous cold streams in hot gaseous halos.Figure 4 and Figure 5 show the results with halo mass 1×10 12 M ⊙ and 3×10 12 M ⊙ respectively, adopting a metallicity of  = 0.1 ⊙ .The first and second rows in each figure show the density and temperature of the gas in a slice of depth around 1 kpc.The third and fourth rows illustrate the tracer field of gas initially injected at the boundary and the velocity fields, respectively.
In the first ∼ 1.2 Gyr, cold streams fell toward the central region due to gravitation in both simulations.At  ∼ 1.2 Gyr, a tiny fraction of cold gas appears in the central region, i.e.,  < 0.2   (∼ 24 kpc), of the halo with  ℎ = 10 12 M ⊙ .Figure 4 only shows the evolution of the gas in a slice.The time since the existence of cold gas in the central region may be earlier in other slices.Starting from  = 2.4 Gyr, cold gas accumulates gradually and occupies a considerable fraction of the volume in the central region, i.e.,  < 0.2   , of halo  ℎ = 10 12 M ⊙ .Note that we neglected the central galaxies in our simplified model, which would somehow limit the application of our simulation in the very central region (for instance, the inner 10-20 kpc).
The evolution of the tracer field indicates that a small fraction of gas injected at the boundary had reached the central region no later than  = 1.2 Gyr, which was synchronous with the emergence of cold gas in the central region.The velocity field shows that the cold stream has a minor impact on the overall gas distribution, mainly limited to the stream itself and the central region after  ∼ 2.4 Gyr.In the central region ( < 0.2   ), the cold gas has a notable rotational velocity, which should result from the initial angular momentum.
Figure 5 shows the gas evolution in a more massive halo with  ℎ = 3 × 10 12 M ⊙ .On the contrary, the central region remains warm and hot till the end of the simulation.Although there is some residual cold gas around 0.2   after  = 2.4 Gyr, the cold gas is quickly heated again.Meanwhile, the tracer field suggests that some of the gas initially injected at the boundary had also arrived at the central region of  ≲ 0.2   since  = 1.2 Gyr.However, their presence in the central region of the halo does not lead to the accumulation of cold gas therein.On the other hand, the cold stream introduces moderate changes in the gas velocity, mainly around  ∼ 0.2   .
In both simulations, the features of KHI were not observed.The radius of the simulated cold stream  , can be roughly estimated by the fraction of the volume occupied by the cold gas in a particular shell.We find that the radii of the cold stream in our simulation,  , , are usually significantly greater than the critical radius  , .For instance, in the 1e12-256-pcn-HZ simulation,  , ∼ 40 and 500  , at  = 1.0 and 0.5   , respectively.Therefore, the role of KHI should have been suppressed by radiative cooling, which is consistent with the estimates in Section 2.2.However, it is worth noting that the limited resolution in our simulations could, to some extent, impede the growth of KHI.Furthermore, the development of KHI might have been attenuated due to the high Mach number and the considerable density contrast.
To track the evolution of cold streams, we have measured the averaged number density, radial velocity, and pressure of gas in three thermal phases as functions of distance from the halo center.The radial profiles of gas properties in halos with  ℎ = 1 × 10 12 M ⊙ and 3 × 10 12 M ⊙ are shown in Figure 6 and 7 respectively.
The density profiles in Figure 6 indicate that cold and cool gas appears in the central region ( < 0.2  ) of the halo  ℎ = 1 × 10 12 M ⊙ since  ∼ 1.0 Gyr, and accumulates gradually after that.During their path toward the center of the halo, the inward radial velocity of the cold gas increases from around 200 km/s at  ∼   to ∼ 400 km/s at  ∼ 0.1   , and then decreases sharply to ∼ 0 km/s at a smaller radius.The acceleration and deceleration are driven by the gravitation and pressure gradient, respectively.Given the temperature and velocity of gas, the cold and cool gas is falling inward with a supersonic speed (the Mach number of the gas in the cold stream   =   /  varies from 12.3 to 24.6, and the background Mach number   =   / ℎ varies from 1.23 to 2.47) in the range 0.1  <  <   .Meanwhile, as can be seen from the velocity and pressure plots in the second and third rows, a tiny fraction of the gas had turned to the warm phase, likely from the cold and hot phases.On average, 10% of the hot gas turns into the warm phase.The warm gas moves inward at a speed moderately slower than that of the cold gas while it is in thermal equilibrium with the hot phase.
Figure 7 indicates a different history for the halo with  ℎ = 3 × 10 12 M ⊙ .Cold and warm gases can not exist in the central 20 kpc, i.e., 0.1   , till the end of the simulation, despite their radial velocities in the outer region being faster than in  ℎ = 10 12 M ⊙ .Therefore, the cold and warm gas at r>20 kpc should have been heated to the 'hot' phase.As we have discussed in the previous context, the third row in Figure 5 demonstrates that some of the gas that was initially injected at the boundary as the cold stream had arrived at the central region earlier than  = 1.2 Gyr.However, their presence does not lead to the emergence of cold and cool phases in the central region of r<20 kpc, i.e.,  ≲ 0.1  .
The bottom row of Figure 6 and 7 indicate that the thermal pressure of the cold stream is lower than that of the hot halo gas.This discrepancy can be traced back to the initial conditions for the cold streams and the hot halo at the boundary.The parameters settings for the hot halo are detailed in 2.1.1,and the cold streams are listed in 2.1.2.The adopted properties of the cold streams and hot halo gas are similar to previous cosmological hydrodynamical simulations.It is important to note that the thermal pressures of the two phases are not expected to be identical, given their distinct motions.The hot phase remains nearly static, while the cold phase descends inward with a velocity exceeding 100 km/s.In addition, the imbalance in thermal pressure between two phases should have shaped the cone-like ge-ometry of a cold stream, as illustrated in previous works (e.g., Dekel et al. 2009;Nelson et al. 2013).

Competition between cooling and compressional heating
In the case of spherical accretion, Birnboim & Dekel (2003) had demonstrated that the rivalry between radiative cooling Λ  and the compress of the gas determines the stability of the accretion shock.According to the estimation on the typical time scales of multiphysical processes given in Section 2.2, the competition between radiative cooling Λ  and compressional heating caused by gas pressure should be the key to determining the thermal evolution of gas in the cold stream.The compressional heating rate can be given by (Birnboim & Dekel 2003) Before we present the detailed results on cooling and compression heating, we should remind the readers again that this work focusses on the evolution of cold streams in the CGM environment.In this work, the model in the central 0.2   region is oversimplified.If there are other heating sources in the central region, such as feedback from star formation and supermassive black holes, they should be added to the RHS of the above heating equation.The impact of central feedback on the cold stream is beyond the scope of this work.
We conduct a quantitative study on the gas cooling and compressional heating rates in our simulations.Moreover, we divide the gas into two categories, one with tracer density greater than 0, i.e.,    > 0, and the other is tracer-free.The former category contains gas injected at the boundary, while the latter is the origin of hot halo gas that has not yet mixed with the injected gas.From the snapshots of the simulation, we use post-simulation analysis to measure the mean cooling rate, compressional heating rate, the ratio of cooling/compression, and the average temperature of the two categories of gas as functions of the radius from the halo centre.The corresponding values are volume-weighted means that take gas cells within different radius bins, with an interval of  = 0.01   .The results of the run with halo mass  ℎ = 1 × 10 12 M ⊙ and gas metallicity  = 0.1 ⊙ are shown in Figure 8. From the left to the right columns, the results at times 0.8, 1.0, 1.2, and 1.4 Gyr are illustrated by turns.The lines in the first three rows are colour-coded with the average temperature.The solid triangle (circle) indicates the gas with    > 0 (   = 0).Because not all gas cells experience compression, there are some negative values for the compressional heating rate and the cooling/compression ratio.We mark these data points with downward triangles and pentagrams, showing their absolute values.
In the first 0.8 Gyr, the gas in the cold stream is still in its path, falling toward the centre.At the same radius, the mean cooling rate of gas with tracer density    > 0 is larger than their counterparts with    = 0 by about 2-3 orders of magnitude.The probable reason is as follows.The gas with    > 0 is the mixture of the injected cold gas and the origin hot halo gas and, therefore, has a temperature between 10 4 K (cold stream) and ∼ 10 6 K (hot halo).When the temperature of mixed gas approaches the plateau in the cooling curve shown in Figure 2, i.e., 10 4.9 − 10 5.4 K for metallicity  = 0.1  ⊙ , the cooling rate would increase dramatically.The mean cooling/compression ratio of the mixed gas ranges from 1 to 100, which is also ∼ 100 times higher with respect to the gas with    = 0 at the same radius.In the central region where the gas in the cold stream has not reached out (there is no sign of tracer at  < 0.2  ), the mean cooling/compression ratio is around 1, indicating that the mean cooling rate is either marginally balanced with or is moderately higher than the mean compressional heating rate for the hot gas therein.This is consistent with the results of relevant timescales presented in Section 2.2.
At  ∼ 1.0 Gyr, the frontier of the cold stream has reached the very central region of the halo.Since then, the cooling of the mixed gas in the central region has been boosted and its rate has grown considerably higher than the compressional heating rate.The corresponding cooling/compression ratio can be as large as 100 at  ≳ 1.0 Gyr.This remarkable feature mainly results from the temperature of the mixed gas in the central region becoming lower than 10 6 K and climbing to the peak of the cooling curve.Later, the mixed gas within  < 20 kpc experiences runaway cooling at a cooling/compression ratio greater than 100 since the temperature of some mixed gas is close to the plateau of cooling rate, i.e., 10 4.9−5.4K.This result also agrees with the analysis of cooling and compressional heating timescales introduced in Section 2.2.Consequently, the mean temperature of mixed gas can be lower than 10 5 K in the central region since  = 1.4 Gyr.
The cooling and compressional heating rates of the gas in the halo  ℎ = 3 × 10 12 M ⊙ are shown in Figure 9.In the first 1.0 Gyr, the cooling efficiency of gas is similar to that of the halo  ℎ = 10 12 M ⊙ .At  ∼ 1.0 Gyr, cold gas injected at the boundary had arrived at the inner region of the halo, as indicated by the tracer field.In contrast to the case in the halo of  ℎ = 10 12 M ⊙ , the increment in the cooling rate of mixed gas was moderate and comparable with the rise of compressional heating rate.Therefore, the cooling/compression ratio of the mixed gas fluctuates around 1.0, which is lower than that in the outside region.Consequently, the gas in the central region of the halo remains warm and hot till the end of the simulation.The difference between the two halos mentioned above is that a deeper potential well hosts a gaseous halo with higher temperatures, leading to a larger falling velocity and stronger compressional heating.Thus, the temperature of the mixed gas is far away from the plateau in the cooling curve for a halo with  ℎ = 3 × 10 12 M ⊙ , resulting in a significantly longer cooling time.

Supply of cold and cool gas
Cold streams play an essential role in galaxy formation and evolution because they serve as the channel to provide cold and cool gas for the central galaxy in halos with a mass below and moderately higher than 10 12 M ⊙ at  ≳ 1.The critical factor is the mass rate of the cold and cool gas at which the central region of the halo can obtain.We probe the quantitative change of the gas mass in three phases, i.e., cold, warm, and hot, at three radii, i.e., r = 1.0, 0.5, and 0.2   respectively.We have measured the changes in two different ways.One is the averaged change of gas mass within a certain radius between two successive snapshots and is denoted as: where  (, ) is the total mass of gas in a particular phase within the radius  at the simulation time ,  is the time interval between snapshots.The other method describes the instant inflow rate at the spheres with three given radii, i.e., The evolution of ¯  (, ) and  (, ) in the two simulations with gas metallicity  = 0.1 ⊙ , and halo mass  ℎ = 10 12  ⊙ and 3×10 12  ⊙ are shown in Figure 10 and 11 respectively.In each figure, the upper row indicates ¯  (, ) at r = 1.0   (red line), r = 0.5   (green line) and r = 0.2   (blue line).Meanwhile, the bottom row presents the instant inflow rate  (, ).From the left to right columns, the plots show the measurements of cold and cool, warm, hot, and all the gas, respectively.
For the halo with  ℎ = 10 12 M ⊙ , the instant inflow rate suggests that the cold and cool gases begin to flow into  = 1.0, 0.5, and 0.2   at the time 0.6, 0.8 and 1.0 Gyr, respectively.From  = 1.4 Gyr to the end of the simulation, the inflow rate of cold and cool gas remains around 3.0 M ⊙ /yr at  = 1.0   , and ∼ 4 M ⊙ /yr at  = 0.5   .At  = 0.2   , the accretion rate of cold and cool gas increases from 0 at  = 0.8 Gyr to 50 M ⊙ /yr at  = 1.0 Gyr, and then decreases sharply to 5.0 M ⊙ /yr at  = 1.4 Gyr.This violent and short change is likely triggered by the arrival of the frontier of the cold stream at  = 0.2   , combining the results presented in Figures 4, 6 and 8. Between  = 1.4 and  = 4.0 Gyr,  (, ) of cold and cool gas at 0.2   stabilise at ∼ 5 M ⊙ /yr.Note that the injection rate of the cold gas at the boundary is ∼ 3 ⊙ /, which is equal to the instant accretion rate of cold and cool gas at  = 1.0   , but lower than the accretion rate at  = 0.5 and 0.2   by ∼ 30 − 50%.Cooling should be responsible for the increment from  = 1.0   to  = 0.5   and  = 0.2   .
On the other hand, the mass of cold and cool gas barely grows in the halos ( <   ) before  = 1.0 Gyr, as shown by the top-left panel in Figure 10.From  = 1.2 Gyr to  = 2.2 Gyr, the total mass of cold and cool gas in the central region, i.e.,  ≤ 0.2   increases at a rate greater than 50 M ⊙ /yr.The dominant cause should be the runaway cooling of the hot gas after mixing with the cold gas, as demonstrated by the plots in the second and third columns.After  = 2.2 Gyr, cold and cool gases accumulate steadily within  < 0.2   , at a rate approximately equal to the instant inflow rate at  = 0.2   .The mass of warm gas within three radii changes dramatically between  = 1.0 Gyr and  = 2.0 Gyr.The probable reason is that the warm phase is the temporary transition phase between cold and hot and cannot exist stably.During the epoch  = 1.0 and  = 2.0 Gyr, a considerable amount of hot gas first cools to the warm phase, resulting in a peak increment of warm gas, and then the warm gas is further cooled to the cold gas, resulting in a dip of warm gas at 1.8 Gyr.Between  = 1.0 Gyr and  = 2.0 Gyr, the mass of hot gas within the halo experiences a violent decline.In addition to the runaway cooling in the very central region, there are several other driving factors.Any imbalance between the gravitation and pressure gradient would lead to the inflow of hot halo gas.Furthermore, part of the gas in the cold stream was heated up to the hot phase.
We then extended our investigation to the case where  ℎ = 3 × 10 12 M ⊙ .The bottom left panel in Figure 11 shows that cold and cool gas flows steadily into spheres with radii of  = 1.0 and 0.5   at rates of 3 M ⊙ /yr and 4 M ⊙ /yr from approximately  ∼ 0.4 Gyr to the end of the simulation.However,  (, ) of cold and cool gas at the radius of  = 0.2  is positive only in the time range from t=1.6 Gyr and t=2.8 Gyr.Subsequently, this inflow rate decreases to nearly zero after t=2.8 Gyr.It is worth noting that the accumulation of cold and cool gas within  = 0.2   occurs at a considerably slower pace compared to the case with  ℎ = 10 12 M ⊙ , which is also indicated in the first row of Figure 7. Furthermore, no cold and cool gas is found within  ≤ 20 kpc at  = 4 Gyr.Along their streamlines towards the centre of the halo, the gas injected from the boundary experiences a more significant heating due to compression than those in halos with  ℎ = 10 12 M ⊙ .Interestingly, the variation in the warm gas throughout the simulation is rather gentle.The hot gas exhibits a notable inflow into the central region ( ≤ 0.2  ) during the initial 2 Gyr, gradually slowing down thereafter.This inflow can probably be attributed to the imbalance between the gas pressure and the gravitational potential.
Briefly, we find that cold streams are highly effective in transporting cold gas to the central region of a halo ( < 0.2,   ) with a mass of  ℎ = 10 12 M ⊙ .The delivery rate is equal to and sometimes exceeds the injection rate at  >   .Moreover, the cold stream can lead to a runaway cooling of the hot gas within 0.2,   of the halo with  ℎ = 10 12 M ⊙ .However, in the case of a halo with a larger mass,  ℎ = 3×10 12 M ⊙ , the cold stream is substantially disrupted by the heating due to compression, occurring at a radius slightly larger than 0.2   .Consequently, it fails to supply cold and cool gas into the central region of the halo.3 and Section 2.2.The left and right panels present ¯  ( < 0.2  ,  ) and  ( = 0.2  ,  ) respectively.See Equations 21 and 22 for definition.Note that a lower limit of 0.01 is set for dm/dt to enable full coverage in the left panel.Meanwhile, to prevent overlap, the lines have been slightly shifted upward in both panels.The magnitude of the shift for each line can be determined by the results at  < 0.5 Gyr, where the actual value is 0 for all simulations.The flat lines in both panels indicate no cold gas flow into or within the central 0.2   region.A value exceeding the value at  < 0.5 Gyr signifies cold gas inflow or generation of cold gas through cooling within the central 0.2   region at that particular time.

The transition halo mass and impact of metallicity
The results presented in the last sections show that the fate of cold streams in our simulations is mainly determined by the rivalry of cooling and compression of falling gas.The latter is governed by the halo mass, whereas the former is closely related to the metallicity of gas.This result is generally consistent with the bimodal accretion scenario proposed in previous theoretical and cosmological simulation studies (e.g., Birnboim & Dekel 2003), where spherical accretion was studied.Previous studies have concluded that there is a maximum halo mass,    , below which a cold stream can exist.To explore the transition halos mass that the cold stream no longer survives in the central region of halos for the metallicity Z=0.1  ⊙ , we have also performed a simulation with the same metallicity but for a halo mass 1.5 × 10 12 M ⊙ .For the sake of clarity, the plots of the slices and profiles are not shown.Figure 12 shows the accretion rate of cold and cool gas at 0.2   for different simulations.Note that, to avoid visual overlap, the lines have been slightly shifted upward in both panels of Figure 12.The magnitude of the offset for each line can be determined by the results at  < 0.5 Gyr, where the actual value is 0 for all simulations.In four simulations, there is no cold gas supply to the central region until the end.In the other five simulations, the cold gas accretion rate shows significant fluctuations.In addition, the limited number of stored snapshots have caused occasional spikes.For the simulation with  ℎ = 1.5 × 10 12 M ⊙ and Z=0.1  ⊙ (1.5e12-256-pcn-Hz), the instant accretion rate remains 0 till  ≈ 0.8 Gyr, then increases rapidly to ∼ 15M ⊙ /yr at  ≈ 1.0 Gyr, but then decreases to ∼ 5M ⊙ /yr afterward.However, the total mass of cold and cool gas in the central region only gains between 0.8-1.4Gyr and after 3.4 Gyr.That is, the cold gas supply to the central region is not continuous.Therefore, we speculate that the transition mass    is very likely close to 1.5 × 10 12 M ⊙ for Z=0.1  ⊙ .
The simulations introduced in the previous sections adopt a gas metallicity Z=0.1  ⊙ .To investigate the effect of metallicity, we have performed another set of simulations for halos with masses of 1.0, 1.5, and 3.0 × 10 12 M ⊙ , but with a lower metallicity of Z=0.001  ⊙ .Plots of radial profiles, cooling and compression, and accretion rates for the halo mass 1.0 ×10 12 M ⊙ can be found in the appendix C. The accretion rate of cold and cool at 0.2   for these three simulations is also shown in Figure 12.For a halo mass 1.5 × 10 12 M ⊙ , the cold stream cannot bring cold and cool gas into the central region with Z=0.001  ⊙ .Nevertheless, a cold stream can deliver cold and cool gas into the central region of the halo with  ℎ = 1.0 × 10 12 M ⊙ , although the accumulation of cold gas only starts at  ∼ 3.4 Gyr.In contrast to the high-metallicity case, runaway cooling is not observed.We conclude that the transition halo mass    for Z=0.001  ⊙ is between 1.0 × 10 12 M ⊙ and 1.5 × 10 12 M ⊙ , somewhat lower than the case with Z=0.1  ⊙ .
Obviously, gas metallicity also plays a nonnegligible role in the evolution of cold streams.Higher metallicity results in a higher cooling rate and shorter cooling timescale, as indicated in Table 4.The strength of compressional heating remains constant for a given halo mass, but the gas will cool more rapidly with increased metallicity.Consequently, the cold stream can penetrate the hot gaseous halo, causing the transition mass to rise with the gas metallicity.Our study assumes a uniform metallicity for the injected cold gas and the surrounding hot gaseous halo.However, it is important to note that the IGM, CGM, and ISM can exhibit varying metallicities.For a more refined modeling approach, it may be necessary to consider gas metallicity as a variable and track its evolution.
The transition mass found here for both metallicities is in general in agreement with the theoretical expectation (e.g., Dekel & Birnboim 2006) and recent results obtained by Daddi et al. (2022a,b), i.e. 2 × 10 12 M ⊙ .

DISCUSSION
In addition to the simulations presented in the last section, several simulations have been carried out to explore the effects of thermal conduction and resolution.We discuss the impacts of these factors and the limitations of our study in this section.

Conduction
Thermal conduction has been found to play an important role in the evolution of cold gas clumps within a hot gas environment (e.g., Armillotta et al. 2017).Thermal conduction may impede the development of hydrodynamical instabilities at the interface of the stream and CGM (Ledos et al. 2023), keeping the stream compact and therefore more difficult to destroy and thus helping the streams penetrate the central region of halo in the cold phase.One of our simulations, '1e12-256-pcc', has turned the thermal conduction on, in contrast to the simulation '1e12-256-pcn'.Figure 13 compares these two simulations on the evolution of the mean density, radial velocity, and pressure of three gas phases.Thermal conduction enables the stream to flow deeper into the central region of halo (reach out r= 0.1  ) with a somewhat higher speed at the end of the simulation.But, as we can see, including thermal conduction has induced a minor impact on the overall evolution of cold streams.In addition, a slight change can be found in the supply of cold gas to the central region, as shown in Figure 12.The main reason is that cooling is more efficient than thermal conduction, as manifested by the corresponding time scales shown in Section 2.2.
The role of thermal conduction may probably have been underestimated as a result of the limited spatial resolution.In Armillotta et al. (2017), the resolution could reach scales of around one pc.As mentioned in 2.2.4, a resolution of Δ < 100 pc is needed to fully resolve the small-scale properties of cold gas.In contrast, the resolution in this work is only around one kpc.Therefore, fine structures on scales smaller than one kpc are inevitably smoothed.Consequently, the effect of thermal conduction is likely damped somewhat in the simulation '1e12-256-pcc', as well as the KHI.Increasing the simulation resolution to pc scale, possibly by adaptively refined grids in the future study, would provide a more accurate estimation of the effects of thermal conduction and KHI.However, such an endeavor requires substantial computational resources.Nevertheless, we anticipate that the main results will remain consistent, grounded in considering the typical timescales of relevant mechanisms in our models, as detailed in Section 2.2.

Caveats
In this paper, we have implemented a simplified model to simulate the cold stream in a hot gaseous halo.Recognizing that several factors not included in our simulation may affect our results is necessary.First, our model does not account for the presence of central galaxies and their associated stars and interstellar medium, which will affect gas cooling in the central region.Moreover, feedback from central galaxies may have non-negligible effects on the fate of cold streams.For instance, Nelson et al. (2015) and Correa et al. (2018) showed that feedback can strongly suppress the inflow rate of the gas at all redshifts, regardless of the temperature history of newly acquired gas.Second, halo gas is usually in a multi-phase state due to inflow, outflow, and thermal instability.Here, we have assumed that the halo gas is isothermal and in hydrostatic equilibrium.A more realistic model of the CGM, combined with improved resolution, would provide a more vivid picture of the evolution of cold stream (e.g.Fielding et al. 2020;Stern et al. 2021;Faucher-Giguère & Oh 2023).Thirdly, throughout the simulation, the host dark matter halo and associated gravitational potential well are static.Considering the growth of dark matter mass as a function of time and merger, the strength of compressional heating will increase with time, potentially shortening the lifetime of cold streams.Furthermore, it is essential to recognize that our simulations exclude certain influential factors, such as magnetic fields and extrinsic turbulence, which could more or less influence the evolution of cold streams.The recent work by Ledos et al. (2023) indicates that magnetic fields can potentially suppress the Kelvin-Helmholtz instability and facilitate cold streams to fuel galaxies.Therefore, for a more comprehensive investigation in the future, including magnetic fields in MHD simulations could provide valuable insights.
On the other hand, we assume a constant accretion rate of 3 M ⊙ /yr for each cold stream onto halos with a mass of approximately 10 12 M ⊙ .This value is inferred from the results in cosmological simulations (e.g., Ocvirk et al. 2008;van de Voort et al. 2011;Zhu et al. 2022).However, we did not incorporate the dependence of the cold accretion rate on the halo mass in this study.Some literature suggests a moderate increase in the cold accretion rate for halos with masses ranging from 10 12 M ⊙ to 3 × 10 12 M ⊙ (e.g., van de Voort et al. 2011;Zhu et al. 2022).Cold stream with a higher accretion rate is expected to lift    .However, the increase in cold accretion to more massive halos may be partly caused by more cold streams.A comprehensive study on the inflow rate of a single cold stream in halos with varying masses could offer a more accurate estimation of    .Finally, it is noted that due to limitations in the version of Athena++ employed in this study, some modules conflict with the self-gravity of gas.The latest version of the Athena++ code has resolved this issue, allowing for the inclusion of the self-gravity of gas in future exploration.

CONCLUSIONS
In this paper, we have performed a set of 3D hydrodynamic simulations to track the evolution of cold streams in hot gaseous halos supported by the gravitational potential of host dark matter halos with a mass between 10 12 and 3.0 × 10 12 M ⊙ at  = 2.The general evolution and the cooling and compressional heating rate of cold streams in the simulations have been probed.The accretion rate of cold and cool gas at various radii has been measured to estimate the survivability of cold streams quantitatively.We have explored the impact of dark matter halo mass, gas metallicity, simulation resolution, inhomogeneity, and thermal conduction on the fate of cold streams.Our findings are summarized as follows: (i) For a gas metallicity Z=0.1  ⊙ , cold streams with mass injection rate 3 M ⊙ /yr per each can flow into the central region, i.e.,  < 0.2   , of hot gaseous dark matter halo with mass  ℎ = 10 12 M ⊙ at  = 2.In addition, cold streams induce runaway cooling in the central region.On the contrary, the cold streams could be substantially disrupted by the heating due to compression at a radius around 0.2   in a more massive halo with  ℎ = 1.5 and 3.0 × 10 12 M ⊙ .
(ii) The rivalry between cooling and compressional heating is the decisive factor governing the fate of the cold stream in the hot gaseous halo, which is similar to the case of spherical accretion.In the case of a  ℎ = 1 × 10 12 M ⊙ halo with a gas metallicity Z=0.1  ⊙ , the cooling/compression ratio of a mixed gas of the cold stream and the hot halo can reach 1000, which means a pronounced dominance of cooling.As a result, the cold stream can be sustained to the central region of the halo.However, for a  ℎ = 3 × 10 12 M ⊙ halo, the cooling rate of mixed gas either marginally balances with or falls below the compressional heating rate.Consequently, the cold accretion via stream into the central region of the halo is suppressed.
(iii) The transition halo mass,    , which determines whether cold streams can persist or not, increases with increasing gas metallicity.Higher metallicity increases the cooling rate of gas, thereby enhancing the likelihood of cold stream survival.For Z=0.001  ⊙ , the transition mass is in the range of 1.0 and 1.5 × 10 12 M ⊙ .If the metallicity increases to Z=0.1  ⊙ , the transition mass increases to approximately 1.5 × 10 12 M ⊙ .These values are generally consistent with previous theoretical studies of spherical accretion and recent observations (e.g., Daddi et al. 2022b).
(iv) In halos less massive than the transition mass,    , cold streams can efficiently deliver cold and cool gas to the central region through two mechanisms.The first is the gas within the cold streams that initially entered from outside the halo.The second mechanism arises from the cooling of hot gas within the central region when it mixes with the injected cold gas.In some extreme cases, with the runaway cooling, the latter can dominate over the former.For instance, in the simulation with  ℎ = 1 × 10 12 M ⊙ , the mass growth of cold and cool gas within  = 0.2  can be much higher than the inflow rate at the sphere of  = 0.2  during some stages.This suggests that the cold and cool gas within  = 0.2 vir is mainly sourced from the gas already within that same radius rather than from gas external to this region.Nevertheless, the second mechanism warrants further validation through a more realistic model of the central region of the halo, particularly including a central galaxy and its feedback.
In summary, our investigation provides insights into the mechanisms and processes that disrupt cold streams in intermediate-mass halos.

Figure 1 .
Figure 1.Schematic diagram of the model for cold stream and hot gaseous halo, showing a slice view of the central xy plane along the z direction.The simulation box has a cubic size of (3  ) 3 .Cold gas is injected into the box from the left boundary with an initial velocity along the x-axis, then will flow toward the halo center under the influence of static gravitational potential.

Figure 2 .
Figure 2. Normalised cooling rate as a function of the gas temperature for metallicity Z=0.1 ⊙ (red line) and Z=0.001 ⊙ (blue line).The vertical dotted lines indicate the local maximum around 10 5 K for both metallicity, i.e.,  = 10 4.9 K for  = 0.001 ⊙ and  = 10 5.35 K for  = 0.1 ⊙ respectively.

Figure 4 .
Figure 4. Slices through the central xy plane along z direction of the simulation 1e12-256-pcn-HZ, i.e., with  ℎ = 10 12 M ⊙ , inhomogeneous cold streams, cooling and  = 0.1 ⊙ , but without conduction.The first and second rows show the number density and temperature, respectively.The third row is the tracer weight    (mentioned in Section 2.1.2),and the last row is the velocity field.The three columns from left to right show results at time 1.2Gyr, 2.4Gyr, and 3.6Gyr, respectively.

Figure 6 .
Figure 6.Rows from top to bottom show the averaged number density, radial velocity, and pressure as functions of distance from the halo center, respectively, for the simulation 1e12-256-pcn-HZ.The results at times 1.0, 2.0, 3.0, and 4.0 Gyr are presented from the left column to the right column, respectively.Red, orange and blue lines indicate hot (T≥ 10 6 K), warm (10 5 ≤  < 10 6 K), and cold and cool gas ( < 10 5 K).Vertical dotted lines indicate r= 0.2 (green), 0.5 (orange), 1.0 (blue)   .The black dotted line indicates the simulation results without inflow cold stream.

Figure 8 .
Figure 8.The first, second, third, and fourth row show the mean cooling, compression rate, ratio of cooling/compression, and the average temperature as functions of the radius from halo center, respectively, at time 0.8, 1.0, 1.2, 1.4 Gyr (from left to right) for the simulation 1e12-256-pcn-HZ, in which  ℎ = 10 12  ⊙ .The data points in the first three rows are color-coded by the average temperature.The triangles (circles) represent the average values of gas cells (do not) containing tracers of injected cold streams.In the second row, the downward triangles (stars) indicate that the gas cells (do not) contain a tracer of injected cold streams and have a negative compression rate.Vertical dotted lines indicate the radius of r= 0.2 (green), 0.5 (orange), 1.0 (blue)   .

Figure 10 .
Figure 10.Gas accretion rate for the simulation 1e12-256-pcn-HZ, in which  ℎ = 10 12  ⊙ .The top row indicates ¯  ( ,  ) within r=1.0  (red), r=0.5  (green) and r=0.2  (blue).The bottom row presents the instant accretion rate  ( ,  ) at r=1.0  , r=0.5  and r=0.2  .See equation 21 and 22 for the definition of accretion rate.From the left column to the right, the plots show the cold and cool, warm, hot, and all the gas, respectively.

Figure 12 .
Figure12.The accretion rate of cold and cool gas at 0.2   measured in nine simulations.The setting of different simulations can be found in Table3and Section 2.2.The left and right panels present ¯  ( < 0.2  ,  ) and  ( = 0.2  ,  ) respectively.See Equations 21 and 22 for definition.Note that a lower limit of 0.01 is set for dm/dt to enable full coverage in the left panel.Meanwhile, to prevent overlap, the lines have been slightly shifted upward in both panels.The magnitude of the shift for each line can be determined by the results at  < 0.5 Gyr, where the actual value is 0 for all simulations.The flat lines in both panels indicate no cold gas flow into or within the central 0.2   region.A value exceeding the value at  < 0.5 Gyr signifies cold gas inflow or generation of cold gas through cooling within the central 0.2   region at that particular time.

Figure 13 .
Figure13.Rows from top to bottom show the averaged number density, radial velocity, and pressure of three gas phases as functions of distance from the halo center, respectively, for the simulation 1e12-256-pcn.The results at times 1.0, 2.0, 3.0, and 4.0Gyr are presented from the left column to the right column, respectively.The red, orange, and blue lines indicate hot (T≥ 10 6 K), warm (10 5 ≤  < 10 6 K), and cold and cool gas ( < 10 5 K).Solid (dotted) lines represent simulations with (without) thermal conduction.The vertical dotted lines indicate a radius of r= 0.2 (green), 0.5 (orange), 1.0 (blue)   respectively.The black dotted line indicates the simulation results without inflow cold stream.

Figure B1 .
Figure B1.The same as Figure 6, but for the simulation 1e12-256-pcn and 1e12-512-pcn, in which  ℎ = 1 × 10 12 M ⊙ .The solid line is the result of 256 resolution, the colored dashed line is the result of 512 resolution, and the black dashed line is the control result of no cold inflow.

Figure C1 .
Figure C1.Rows from the top to bottom show the averaged number density, radial velocity, and pressure, as functions of distance from the halo center respectively for the simulation 1e12-256-pcn, in which  ℎ = 10 12 M ⊙ and  = 0.001 ⊙ .Results at times 1.0, 2.0, 3.0, and 4.0 are presented from the left to the right columns, respectively.The red, orange, and blue lines indicate hot (T≥ 10 6 K), warm (10 5 ≤  < 10 6 K), and cold and cool gas ( < 10 5 K).The black dashed lines are the simulation results with the same conditions but without a cold stream.Vertical dotted lines indicate r= 0.2 (green), 0.5 (orange), 1.0 (blue)   .

Figure C2 .
Figure C2.Gas accretion rate for the simulation 1e12-256-pcn, in which  ℎ = 10 12 M ⊙ and  = 0.001 ⊙ .The upper row indicates ¯  ( ,  ) within r=1.0  (red), r=0.5  (green) and r=0.2  (blue).The bottom row presents the instant accretion rate  ( ,  ) at r=1.0  , r=0.5  and r=0.2  .See Equations 21 and 22 for the definition of accretion rate.From the left to the right column, the plots show the cold and cool, warm, hot, and all the gas, respectively.

Figure C3 .
Figure C3.The first, second, third, and fourth column show the mean cooling, compression rate, ratio of cooling/compression, and the average temperature as functions of the radius from halo center, respectively, at time 0.8, 1.0, 1.2, 1.4 Gyr (from top to bottom) for the simulation 1e12-256-pcn, in which  ℎ = 10 12 M ⊙ and  = 0.001 ⊙ .The data points in the first three columns are color-coded by the average temperature.The triangles (circles) represent the average values of gas cells (do not) containing tracers of injected cold streams.In the second column, the downward triangles (stars) indicate that the gas cells (do not) contain a tracer of injected cold streams and have a negative compression rate.The vertical dotted lines indicate r= 0.2 (green), 0.5 (orange), 1.0 (blue)   .

Table 2 .
Parameters for the cold stream: It has a log-normal density distribution with mean value   , temperature   , and an initial velocity of   =100 km/s.The radius of the cross-section at the initial position is 15 kpc, (Gyr)  ,(kpc)