Dissecting Cosmological Filaments at High Redshifts: Emergence of Spaghetti-type Flow Inside DM Haloes

We use high-resolution zoom-in simulations to study the fueling of the central galaxies by gas accretion from cosmological filaments at high redshifts, z>=2. Their parent haloes with similar DM masses of log(M_vir/M})~11.65, have been chosen at z=6, 4, and 2, in high/low overdensity environments, with the goal of comparing evolution within similar M at different z, under dual action of cosmological accretion and galactic outflows -- forming the circumgalactic medium (CGM). We focus on the filamentary and diffuse gas accretion within few virial radii, R_vir, down to the central galaxy. Using a hybrid d-web/entropy method we have mapped the gaseous filaments, and invoking particle kinematics allowed us to separate inflows from outflows, thus resolving thermodynamic and kinematic signatures of the CGM. We find that (1) The CGM is multiphase and not in thermodynamic or dynamic equilibrium; (2) accretion rates via individual filaments display a lower accretion rate and densities at lower redshifts. The inflow velocities along the filaments decrease with redshift, z~ 6-2, from 200-30 kms^-1 by a factor of 2; (3) Temperature within the filaments increases inside R_vir, faster at lower redshifts, in tandem with decrease in the accretion rate; (4) The filaments show a complex structure along their spines: a core radial flow surrounded by a lower density envelope. The core exhibits an elevated density and lower temperature, with no obvious metallicity gradient in the filament cross sections. It also tends to separate the filament into different infall velocity regions and density cores, thus producing a spaghetti-type flow; (6) Inside the inner ~ 30\,h^-1 kpc, the filaments develop the Kelvin-Helmholtz instability which ablates and dissolves them, and triggers turbulence along the filament spine; (7) Finally, the galactic outflows affect mostly the inner ~ 0.5R_vir~ 100 h^-1 kpc of the CGM.


INTRODUCTION
Galaxies reside inside dark matter (DM) haloes which extend to about ten times the galaxy size.The DM haloes lie within the cosmological network formed by filaments, walls and their intersections, which, together with voids, form the so-called large-scale structure of the universe extensively studied over the last 3-4 decades (e.g., Bardeen et al. 1986;Geller & Huchra 1989;Bond et al. 1996;Springel et al. 2006).Deep within these parent haloes, galactic morphology in contemporary universe has been analyzed for the last century with a great success, solidified by Hubble (1936, see reviews by Binney & Tremaine (2008); Mo et al. (2010); Kormendy (2013); Shlosman (2013)).
However, the gray zone of contemporary cosmology lies between these two extremes -how does the large-scale structure -the cosmic web, connect to the galactic structure?How do properties of the infalling material change within the halo, and to what extent it preserves its identity on the way to the central galaxy?How does it feed the galaxy growth and affect their morphology.And finally, how does this process evolve with the cosmological time?
In this work, we attempt to answer more modest questions.We ★ E-mail: dbi224@g.uky.edufocus on comparing the evolution of baryonic filaments and diffuse accretion in the immediate vicinity of galaxies and their host DM haloes, focusing on the kinematic and thermodynamic properties of the gas.To distinguish the baryonic flows from the extensions of DM filaments inside the haloes, we denote them as streamers.
For this purpose we use our high-resolution cosmological zoom-in simulations presented in Bi et al. (2022a,b).Two primary modes of galaxy growth exist: mergers with other galaxies and smooth accretion of matter (e.g., Rees 1977;Fall & Efstathiou 1980).High-redshift galaxies exhibit very high star formation rates for extended periods of time which cannot be supported by galaxy merger events only (e.g., Chapman, et al. 2004;Keres et al. 2005;Genzel et al. 2006;Dekel & Birnboim 2006;Bi et al. 2022a).Smooth gas accretion plays an important and even dominant role in the growth of these galaxies, affecting their morphology.In the past decade or so, numerical modeling confirmed that gas accretion can dominate over galaxy mergers at high redshifts, becoming an important component of the hierarchical scenario of galaxy growth (e.g., Keres et al. 2005;Dekel et al. 2009;Devriendt et al. 2010;Romano-Díaz et al. 2011, 2014).
It is widely accepted that the two complementary mechanisms for gas accretion can operate: the so-called cold and hot accretion modes (e.g., Keres et al. 2005;Dekel & Birnboim 2006).Dekel & Birnboim (2006) have argued that in the hot accretion, the incoming gas, which is supersonic and not virialized, will trigger a shock positioned around the halo virial radius,  vir .In more recent numerical simulations, the position of this virial shock has been found to lie outside the virial radius (e.g., Nelson et al. 2016), and we discuss this point in section 4. The postshocked gas, which is heated up to the virial temperature,  vir , forms a quasistatic envelope, starts to cool and is accreted radially thereafter as it can potentially cool down over less than the Hubble time and hence, contribute as well to the galaxy growth (e.g., Birnboim & Dekel 2003).On the other hand, the cold gas penetrates deep inside the halo and streams towards the central galaxy, without being shocked and heated up, allowing for an efficient delivery of a potentially starforming fuel to the central galaxy disk (Brooks et al. 2009).This can even happen when the hot shock is present (e.g., Dekel & Birnboim 2006;Dekel et al. 2009;Agertz et al. 2009).Actually, recent theoretical work has shown that the cold baryonic accretion rather than the hot accretion rate dominates below a certain critical DM halo mass (e.g., Birnboim & Dekel 2003;Keres et al. 2005;Dekel & Birnboim 2006;Ocvirk et al. 2008;Dekel et al. 2009;Keres et al. 2009) and the cold accretion rate decays steeply down at lower redshifts (e.g., Keres et al. 2005;Dekel & Birnboim 2006;Romano-Díaz et al. 2017).Therefore, it is especially important to investigate the details of the cold accretion flows in the early universe.
Several numerical works focused on the cold flows and the baryon cycle of accretion versus outflow in galaxies.Some of them favored low numerical resolution in order to follow an statistical approach to galaxy growth from gas accretion (e.g., Keres et al. 2005;Ocvirk et al. 2008;Dekel et al. 2009), while others used "zoom-in" simulations to address the process on galactic scales at specific redshifts (e.g., Faucher-Giguere et al. 2011;Fumagalli et al. 2011;Kimm et al. 2011;Stewart et al. 2011;Goerdt et al. 2012;Shen et al. 2013;Romano-Díaz et al. 2017).
In this work, we focus on the baryonic accretion onto the central galaxy within its parent DM halo at the Cosmic Dawn, ending at three selected redshifts,  f = 6, 4, and 2. Using high-resolution zoom-in simulations (Bi et al. 2022a), we follow the galaxy evolution within similar mass DM haloes, log  vir ∼ 11.75 ± 0.05 M ⊙ , which are only a factor of 2 below those haloes which are the most efficient in producing stars (e.g., Behroozi et al. 2013).
The chosen mass range ensures that the selected halos can create favorable conditions for both the hot and cold accretion flows, and are expected to form sufficiently massive galaxies to be properly resolved numerically in our simulations.In addition, these halos are expected to contain  * galaxies at their respective  f .Hence, they can be compared with present day  * galaxies.Our final redshifts,  f , encompass the end of the reionization epoch,  > ∼ 6, and the subsequent time period of ∼ 2.5 Gyr, when the star formation in the universe peaks.Overall, galaxy evolution for  ∼ 9 − 2 is analyzed.
Our main goal is to study accretion of the cosmological gas from the large scale, i.e., from ∼ 4 vir , down to the central galaxies, where the gas settles in the circumgalactic space.We define this gas as the circumgalactic medium (CGM, e.g., Tomlinson et al. 2017, for a recent review).This basically corresponds to the gaseous component of the DM halos (e.g., Sadoun et al. 2019), but its radial extension will be more precisely defined in sections 3 and 4. and we analyze its specific properties, i.e., the kinematics, chemical composition and thermodynamics.Within this spatial range, the gas transform from the extragalactic medium to the starforming ISM.This gas is expected to supply the new material for star formation from the IGM, satellite galaxies, the hot and cold accretion phases, and from the recycled material injected by the disk stars and active galactic nuclei (AGN) -as a result of the galactic feedback processes.This region also has the potential of hiding the missing baryons in the universe.Moreover, the CGM dynamic and thermodynamic states can be altered by the stellar and AGN feedback mechanisms which in this case will be responsible for the quenching of star formation in galaxies (e.g., Davies et al. 2020).
Major efforts have been aimed at understanding the feedback of winds from massive stars, supernovae (SN), and AGN on galaxy evolution.However, specific details of this feedback and its implementation and fine tuning to observational and computational models are far from being understood (e.g., Yajima et al. 2015;Sadoun et al. 2016).In particular, the spatial extent of this feedback is still unclear -is it limited to the inner DM haloes, to the halo virial radius, or extends beyond it?
This paper has been organized as follows.Section 2 describes the numerical issues and the methods used.Section 3 presents our results.This is followed by the discussion section and the summary.

Simulation Suite
We analyze the zoom-in simulations suite presented in Bi et al. (2022a) which used the hybrid -body/hydro code gizmo (Hopkins 2017) with the Lagrangian meshless finite mass (MFM) hydro solver.The full details of the simulations are given in Bi et al. (2022a).Here we only provide the highlights.
Uni-grid, DM-only initial conditions (IC) were generated at redshift  = 99 by means of the music code (Hahn & Abel 2011) within a box of 50ℎ −1 Mpc, which were evolved until redshift 2. From this parent simulation haloes of the same mass range were chosen at three different target redshifts,  f = 6, 4, and 2, to be re-simulated at much higher resolution.The selected DM haloes have been chosen to reside at low and high overdensities, δ ∼ 1 and δ ∼ 3, respectively (Bi et al. 2022a).
The individual zoom-in ICs, also generated with music, were composed of five nested levels of refinement on top of the base grid, i.e., from 2 7 to 2 12 .Their DM-only versions were first evolved in order to check for (and to avoid) contamination from massive, lowerresolution particles at the highest resolution-level volume.After this, baryons were included at the highest level of refinement in the reconstruction of their respective ICs.All the details of the models are displayed in the Table 1, as well as in Tables 1 and 2 of Bi et al. (2022a).The total number of selected haloes were 6, and in combination with the two different galactic winds feedbacks (see section 2.3) the final simulation suite was composed of 12 haloes.
Within this setup, the effective number of particles (DM and baryons) in our simulations is 2 × 4096 3 , leading to a mass resolution per particle of 3.5 × 10 4 M ⊙ for gas and (eventually) stars, and 2.3 × 10 5 M ⊙ for DM.The minimal adaptive gravitational softening (in comoving units) for gas was 74 pc, and for stars and DM 74 pc and 118 pc.This means that at the final redshifts,  f = 6, 4, and 2, the softening for stars in physical coordinates is 10.5 pc, 14.7 pc, and 24.6 pc, respectively.
For better angular momentum conservation and in order to resolve the Kelvin-Helmholtz instability which is expected to develop when the cosmological filaments penetrate the DM halo, the MFM hydro solver was employed instead of the "traditional" SPH solver with an adaptive gravitational softening for the gas.The hybrid multiphase model was invoked for the ISM and star formation (Springel & Hernquist 2003).In this, starforming particles contain the cold phase that forms stars, and the hot phase that results from the SN II heating.Metal enrichment is included: the metallicity increase in the starforming gas scales with the fraction of gas in the cold phase, the fraction of stars that turns into SN, and the metal yield per SN.A total of 11 metal species were followed in both gas and stars, including H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe.Metal diffusion is not implemented explicitly, but metals can be transported by winds (see section 2.3).The density threshold for star formation (SF) was set to  SF crit = 4 cm −3 .Furthermore, our simulations include the redshift-dependent cosmic UV background, which has been turned on at  = 11.7 (Faucher-Giguere et al. 2009).

Identification and Properties of DM haloes and Galaxies
DM haloes and their properties were identified by the group finder rockstar (Behroozi et al. 2013), with a Friends-of-Friends linking length of  = 0.28.The halo virial radius and the virial mass,  vir and  vir , have been defined by  200 and  200 (e.g., Navarro et al. 1996). 200 is the radius within which the mean interior density is 200 times the critical density of the universe at that time.
Galaxies have been identified by the group-finding algorithm hop (Eisenstein & Hut 1998), using the outer boundary threshold of baryonic density of 10 −2  SF crit , which ensured that both the host starforming gas and the lower density non-starforming gas are roughly bound to the galaxy (Romano-Díaz et al. 2014).This assures that identified galaxies are not imposed with a particular geometry.Note that with this definition, all the galaxies appear to be generally smaller than galaxies defined with 0.1 vir , which is typically used in the literature (e.g., Scannapieko et al. 2012;Marinacci et al. 2014).

Galactic Wind Models and Supernovae Feedback
We also have made used of the two wind models of Bi et al. (2022a) -the Constant Wind (Springel & Hernquist 2003) and the Variable Wind (Oppenheimer & Dave 2006) (hereafter CW and VW, respectively).Both wind models implement the concept of "the decoupled particle wind," when a wind particle decouples from the ambient gas particles, no longer interacts hydrodynamically, and move ballistically.The time period of decoupling stage depends on the shortest of either 10 6 yrs, or when the particle moves to a region (but still within the galaxy) where the background gas density is lower by a factor of 10 compared to the critical density of star formation.
Both CW and VW orientations have been assumed isotropic.For the CW model, the wind velocity is  w = 484 km s −1 , and the mass loading factor,  w ≡  w /  SF = 2 (Sadoun et al. 2016).Here  w is the mass loss rate by the wind, and  SF is the star formation rate (SFR) measured in  ⊙ yr −1 .For the VW model, the wind velocity scales with the physical escape velocity of the host halo.In this case, the mass loading factor has been calculated assuming the total wind energy, and it is given by the energy-driven and the momentumdriven winds.In general, by calculating the total kinetic energy of all the wind particles, the VW constitutes a stronger feedback, about 8 times stronger on the average than the CW.Yet, at each moment of evolution this ratio can fluctuate in both directions.
In addition to galactic winds, this modeling includes the SN feedback (Springel & Hernquist 2003).The SN deposition of energy in the gas is fixed at 4 × 10 48 erg M −1 ⊙ .The rate of the SN depends on the initial stellar mass function (IMF), which is taken as the Chabrier IMF (Chabrier 2003).

Cosmic Web Decomposition
The cosmic web can be divided into its various components, i.e., voids, sheets/walls, filaments, and clusters/knots (e.g., De Lapparent et al. 1986;Colless et al. 2003;Tegmark et al. 2004;Mehmet et al. 2014).A number of different methods have been attempted in the literature in both numerical simulations and observational data to separate the web into components (Aragón-Calvo et al. 2007a;Hahn et al. 2007a;Forero-Romero et al. 2009;Bond et al. 2010;Sousbie 2011;Hoffman et al. 2012;Cautun et al. 2013).We have employed a hybrid scheme, involving the d-web, which is based on the orbital stability using the local tidal tensor usually calculated for the DM distribution (Hahn et al. 2007b), as well as two additional methods discussed below.Because we are mainly interested in the evolution of the baryonic component within DM haloes, this method has been applied to scales smaller than few× vir and to the baryonic distribution, rather than to the DM one.In this way, we have been capable to follow the baryonic streamers, i.e., cold streams, inside the virial radius (e.g., Dekel et al. 2009).
We briefly describe the d-web method.A particle i moving in a peculiar gravitational potential,   , can be described by the following equation of motion in comoving coordinates, where the dots represent derivatives with respect to time.After linearizing the equation of motion, the system can be re-written as where  ĳ represents (in our case) the tidal tensor of baryonic gravitational potential within the DM haloes.By calculating the Hessian 1 matrix, we can define different matrix elements by an analogy with the Zel'dovich approximation (Zel'dovich 1970), which depend on the basis of the eigenvalues λ 1 , λ 2 , λ 3 of the tidal tensor, • voids: However, the particular details of this classification method, i.e., extension, thickness of the identified structures, depends on the length scale of the potential field's grid resolution.
As our goal is to extend the filament identification to the baryonic component inside the virial radius, we complement the d-web method outside the virial radius with an entropy-based method inside the such radius, as described below.
We define the entropy as  = / virial radius obtained from the d-web method.This cutoff is used inside the virial radius as well.We have tested this method, and it shows a smooth transitions across the virial radius.Furthermore, we extended the baryonic d-web determination inside  vir , down to 50ℎ −1 kpc, in tandem with the entropy method, for comparison.In addition, inside the virial radius, we apply the baryon kinematics method to distinguish the accretion flow from the outflow, in the region where both coexist.
We have tested the above hybrid method and found that it successfully extrapolates the d-web/entropy elements down to ∼ 10ℎ −1 kpc, where the boundary of the galaxy lies.

RESULTS
We present our results based on our full simulation suite, starting with the detection of cold baryonic streamers and diffuse accretion.Furthermore, we emphasize differences between the two galactic wind feedback models used in our simulations under identical initial conditions, the CW and VW galactic outflow models, and the environmental effects.

Identifying the filamentary streamers
By using our hybrid d-web/entropy method (section 2.4), we have been able to separate the large scale flows (outside  vir ) down to the central galaxy.
Figure 1 displays the mapped individual streamers within a box of side 1, 200ℎ −1 kpc, as computed with the d-web algorithm.Only CW feedback models are shown because at these scales there is no observable difference between both feedback schemes.All models are displayed at their final redshifts.
We follow the filamentary streamers and diffuse accretion flows from ∼ 3 vir down to the inner ∼ 10 ℎ −1 kpc from the halo centers.Our general implementation allow us to overcome the difficulty in following the accretion streams deep into the DM haloes and connecting them to galaxies.Furtheremore, the hybrid d-web/entropy method supplemented by the kinematics inside  vir allows us to separate the inflow from outflow.The streamers have been concatenated to the innermost flow using the baryons kinematic (see Appending sec:append).We discuss the details of gas motions in the innermost regions and their penetration into the central galaxies in section 3.5.
We observe that, overall, the DM haloes are typically associated with two or three main cosmological filaments at  f , as shown in Figure 1.

Radial properties of filamentary streamers
Figure 2 shows the whole-sky Hierarchical Equal Area isoLatitude Pixelisation (HEALPIX) projection maps of spherical shells of 10ℎ −1 kpc thickness at 0.1 vir , 0.5 vir and  vir .It reveals the streamers extending down from  vir to the haloes' central regions.Also shown are the walls (seen as the threads connecting the streamers), and diffuse accretion flow.The color represents radial influx of gas mass per solid angle.Sheets existing between the streamers at larger radii dissolve on the galactic scale, i.e., 0.1 vir , which are especially detectable at  f = 6 and  f = 4.
The color palette in Figure 2 reveals the difference in the accretion rates between the filaments and walls.The former channel up to 100  ⊙ yr −1 rad −2 , while the latter less than 40  ⊙ yr −1 rad −2 .Furthermore, notice how the smooth accretion (not from filaments The color coding represents: knots in blue, filaments, (i.e., streamers) red, walls green and voids in black.From left to right columns, haloes at their final redshifts  f = 6, 4, and 2 are shown.All the models correspond to the CW feedback as they are similar to the VW models at these scales.The upper panels display haloes in the high overdensity environments and the lower panels show haloes in low overdensity environments, (see Table 1).The arrowhead at the Z2L snapshot points to the filament #1 which is analyzed in detail in the subsequent sections.and, or walls) comes from all directions, although at a much lower rates (few × ⊙ yr −1 rad −2 ). Figure 3 reveals that the streamers differ in their density distribution at their final redshifts  f .At lower redshifts, the gas density profiles are lower than those at higher , with the decline becoming more prominent between  f = 4 and 2, as this represents the longest time interval between different  f .At  f = 6, their density ranges in the interval of  ∼ 10 −23 − 10 −25 g cm −3 over the distance of 400ℎ −1 kpc.At  f = 4, this density is about a factor of 3 lower.At  f = 2, the density is even lower and ranges between ∼ 10 −25 − 10 −27 g cm −3 .These results are indicative of a two-fold effect: the decrease in the gas mass flux in the streamers and the overall decrease of density in the universe with a decreasing redshift.Notice, however, that no trend is observed between the high and low overdensities models in this Figure .Figure 4 displays the radial distributions of temperature in the filamentary accreting gas (the diffuse accretion flow properties are discussed in the next section).We distinguish a different behavior inside  vir (∼ 200 ℎ −1 kpc) and outside it, with a transition region between 0.7 vir −1.3 vir .Some differences can also be seen between the CW and VW models.For  f = 6 models, the temperature is rising relatively mildly inside the virial radius.At  f = 4, this rise is much steeper.Streamers at  f = 2 are hotter everywhere, and still display an additional rise towards the center.In the same way, the central temperature increases towards lower  f .We also observe a larger amount of hot gas inside  vir at lower redshifts and for lower overdensity models.At higher redshifts, the gas accumulates mostly at  ∼ 10 4 K.The virial temperature of high-redshift haloes is higher than for the low redshift ones, which is   expected as these haloes are more centrally concentrated.Moreover, the flat part of the -curves broadens to ∼ few × 10 5 K at lower redshifts.
Figure 5 exhibits a diverse radial profile behavior of the gas metallicity.While the filamentary streamers fit in the range of (/ ⊙ ) ∼ 10 −4 − 1, the radial gradients sometimes exhibit opposing trends, i.e., the metallicity can increase or decrease with radius.This can be understood when the metal production occurs in small haloes embedded within the streamers at different radii.Note that the surrounding massive haloes have been removed, but the small ones and substructures persist in these Figures.Hence, the metal contamination contributes due to the outflows and the SN feedback.We also observe that the VW models appear to be less metal rich compared to the CW models, especially in the low overdensity models.This is the result of the VW models being more gas rich due to a reduced star formation everywhere (Bi et al. 2022a).
Figure 6 provides the radial profiles of the inflow velocities along the filamentary streamers.The obvious trend which shows up is that these velocities increase at smaller radii but decline overall at lower redshift, i.e., become less negative inside  vir .Around 50 kpc/h from the center they decrease sharply with decreasing  f .This effect is a direct consequence to DM halos of the same mass being less centrally concentrated, as they form in a lower density universe.
Also, inside ∼ 100 kpc/h, the difference between the CW and VW becomes large, which is probably related to stronger shocks in VW models due to the larger fraction of gas in the VW models.At  f = 6, the inflow velocity in CW models ranges within  r ∼ 100 − 300 km s −1 .At  f = 4, this range becomes smaller,  r ∼ 50−200 km s −1 .While at  f = 2 it decreases to  r ∼ 70−150 km s −1 .
The mass accretion rate of low-temperature gas onto the halo provides an important reservoir of gas that can join the galaxy and contributes to its SFR. Figure 7 shows the radial profiles of the accretion rate from baryonic filaments, , for all models at  f = 6, 4, and 2. The individual streamer contributions have been added up.
The  appears flat.
The observed variations in the mass accretion rates in Figure 7 are limited within a factor of 2 at smaller radii for  f = 6 galaxies and for high overdensity models at  f = 4, and are associated to the presence of substructures.At  f = 2, accretion rate in filaments is either flat or decreases at small radii for the high density haloes, with a slight increase for the low-density ones.
As expected, the accretion rates along the streamers decrease with decreasing  f .For  f = 6 it lies within  ∼ 50 − 100 M ⊙ y −1 .For  f = 4, it ranges within ∼ 10 − 100 M ⊙ y −1 , and for  f = 2, it ranges within ∼ 5 − 10 M ⊙ y −1 .The spikes visible in this Figure are the result of substructure and should be ignored.
This decrease in the accretion rates towards smaller redshifts appears in tandem with the decrease in the density within the filaments and the associated inflow velocities discussed above.We also observe this trend in the accretion rates of the diffuse gas, shown in the next section.However, we do not see any systematic differences between the accretion rates in the CW and VW models.This issue will be addressed in section 4.
In order to visualize the accretion flow on scales smaller than 50ℎ −1 kpc, Figure 8 shows whole-sky HEALPix maps on scales of 5, 10, 30, 40, and 50ℎ −1 kpc, using shells of 3ℎ −1 kpc thickness and colored by the gas surface density.In this Figure, we can detect the filamentary and diffuse accretion flows, but the walls that have been observed in Figure 2 on larger scales, have been washed out here.The inner most panels show the gas distribution within the galaxies (aligned along the equatorial planes) and its immediate surroundings.The presence and strength of the streamers decrease at smaller radius until  ∼ 10ℎ −1 kpc (typical radii of our galaxies) and also as a function of redshift.This situation is also independent of the wind mechanism employed.At this region, some streamers connect directly with the galaxy, but some others might miss it becoming outflows.Such interplay can be noticed by the intricate density pattern at the 10ℎ −1 kpc shell for all models, different and almost uncorrelated with respect to their external shells.
Figure 8 also shows that some of the material in the filamentary accretion misses the central galaxy, it is subsequently diverted and becomes a filamentary outflow.This is detected in the analyzed velocity field not shown here, and is typical for all the models.This outflow does not escape from the halo, reaching < ∼ 150ℎ −1 kpc.The gas in these outflows increases its entropy and is converted into a diffuse virialized gas.In addition, this Figure includes contribution from galactic outflows, i.e., CWs and VWs.

Radial properties of diffuse accretion flow
Next, we turn to the radial properties of the diffuse accretion flow in the vicinity and inside the DM haloes, < ∼ 2 vir .We follow their density, temperature, radial velocity, metallicity and diffuse accretion rate.
Figure 9 shows the density profile of the diffuse accreting gas.Its radial dependence is similar to that of the filamentary gas, but shifted down by a factor of ∼ 5 − 20.These profiles appear to be flat between (1 − 2) vir and  ∼  −1 inside the virial radius.The explanation to this behavior is similar to that in the filamentary accretion mentioned in section 3.2 -decrease in the gas supply as well as the expansion of the universe, both contribute to this effect.We observe no differences between the density profiles based on environment neither with respect to wind model.
The temperature in the diffuse accretion gas (Figure 10) differs from that in the filaments (Figure 4).The lion share of the gas mass in the filaments has a temperature of ∼ 10 4 K, except for  f = 2.The diffuse gas however, has such a low temperature only outside the virial radius.At smaller radii, the majority of the gas lies at  ∼ few × 10 5 − 10 7 K.The cold diffuse flow is ceased to exist inside  vir .At  f = 2, the amount of hot gas outside  vir is also increased.Some differences among the models at   regarding their local overdensities can be noticed, in particular the amount of diffuse gas around and beyond  vir , it is much larger the amount of diffuse gas for the high-density environment haloes (as expected) with respect to their low-density counterparts.In summary, the thermodynamic state of the gas in the filamentary and diffuse accretion differ, and this can have implications for the gas kinematics as well, due to the buildup of thermal pressure gradients.
Figure 11 exhibits the mass-weighted metallicity radial profiles of the diffuse accreting gas.The profiles are relatively flat with significant variations related to the embedded substructures and small neighboring haloes.Simulations with the VW feedback exhibit a lower metallicity than those with CW feedback, especially outside the virial radius.This is mainly due to the overall lower SFR in these models when compared to the CW ones.Furthermore, this is quite similar when compared with the filamentary accretion flows (Figure 5).Both distributions display a substantial fraction of low metallicity gas, / ⊙ < ∼ 10 −3 and even pristine gas outside  vir .The radial inflow velocity profiles of the diffuse accreting gas, while following the same trend as the filamentary accretion in general, differ from the gas flow in the filaments in detail (Figure 12).Namely, the velocities of the diffuse gas appear smaller at lower  f , but their radial velocity gradient is not as large as in the filaments shown in Figure 6.Overall, as in the filamentary accretion, the dif-     The VW (CW) models are given by the solid red (blue) lines.For both inflow and outflow rates, we have used the radial velocities exceeding the local dispersion velocities in the gas.
fuse gas accretes in more concentrated halos at smaller redhift.But the differences in the thermodynamic state between the diffuse and filamentary gas lead to variations in the innermost inflow velocity profiles.
Next, we calculate the mass accretion rates, , of the diffuse gas (Figure 13), and separate it from the outflows triggered by the feedback from the central galaxies.To distinguish the inflow and outflow from the virialized gas inside the DM haloes, we have counted only the inflow and outflow velocities which exceed the local gas velocity dispersions.The resulting mass accretion rates for the diffuse gas decrease with approaching the central galaxies.They also decline with decrease of  f from  ∼ 10 − 40  ⊙ yr −1 at  f = 6 to ∼ 5 − 10  ⊙ yr −1 at  f = 2.We do not see differences between the CW and VW models, except at the  f = 2 low overdensity Z2L model.
The outflows are clearly focused on the central galaxies and become more extended, i.e., reaching larger radii, at lower  f .At  f = 6, the outflows reach ∼  vir , at  f = 4, the VW outflow is reaching ∼ 1.5 vir , while at  f = 2, in the low-density environment haloes the CW reaches ∼ 1.3 vir while VW extends to ∼ 1.8 vir .In general, the outflows seen at larger radii belong to small haloes located well beyond  vir .
The outflow radial profiles decrease with radius.At the maximum, the outflows are ∼ 5 − 10  ⊙ yr −1 at all  f = 6, but decline below 0.1  ⊙ yr −1 at the maximal radii mentioned above.But in all cases, we find that the outflow extends to the virial radius, at least.In addition, some of the inflow misses the central galaxy and is diverted out entirely because of gravitational effects.It interacts with the inflow even outside  vir as we discuss in section 4. Based on this result, we call this region, which extends from the HOP-defined galaxy to the baryonic backsplash radius as the CGM.Thus, the CGM is affected both by outflows and inflows, and, therefore, has a complicated kinematics and other thermodynamic characteristics.

Dissecting the streamers: emerging spaghetti-type flow
While analyzing the radial profiles of various thermodynamic properties of streamers at specific times is necessary for obtaining the global picture, following the time-dependent processes is required as well (e.g., Arzoumanian et al. 2011).As a next step, we investigate the radial motions within the streamers in the range of ∼ 50 − 350 ℎ −1 kpc from the central galaxy, and determine their structure at different radii by dissecting them into slices of 10 ℎ −1 kpc radial thickness.This is complemented by the next section (section 3.5), where we focus on the inner 50 ℎ −1 kpc, where the inflow gradually dissolves, changing its kinematics before penetrating the central galaxy.
Analysis of the streamers radial profiles and their cross section structures at different radii is important for our understanding of their evolution within the DM haloes.It also reveals the fate of these streamers at radii comparable to galaxy sizes.As a prototype, we use the streamer #1 of the Z2L halo identified in Figure 1, and analyze it down to ∼ 10 ℎ −1 kpc.For this purpose, we select the gas particles outside  vir at  ∼ 2.4 and trace them inwards until  f = 2.This streamer falls into the central galaxy being inclined with respect to the disk plane by an angle of ∼ 30 • .
Properties such as the volume density, temperature, metallicity and the radial velocity of this gas show an interesting structure as seeing in Figure 14 (upper plot).Most prominent is the inner, central region of the streamer, i.e., the radial spine, at almost all radii, where one can define the core of the streamer for each of these variables.We also observe that the above physical properties correlate along the streamer.For example, the volume density and the temperature show a tight correlation.However, we do not see a gradient of metallicity between the core and the envelope of the streamer.
To verify this, we have dissected this streamer perpendicularly to its spine.Figure 15 provides a more revealing view to study the gradients between the core and the envelope of the streamer, and their radial dependence.For the volume density, the cross section at the innermost radius  ∼ 50 ℎ −1 kpc (in the midst of the inner CGM) displays a high density core of  ∼ 10 −25.5 g cm −3 , while somewhat less denser at 150 (in the midst of the outer CGM), 200 (the virial radius), 250, and 350 ℎ −1 kpc.The envelope density of the streamer drops by more than one order of magnitude compared to that at the inner boundary of  ∼ 50 ℎ −1 kpc.The core density has a weaker dependence on radius than the envelope.
Small scale structures and substructures are typically embedded in the streamer.For example, the slice situated at  ∼ 350 ℎ −1 kpc of Figure 15 (i.e., the lowest frame) displays a structure embedded in the spine of the streamer.This trapping of the structure in the core of the streamer explains partly its very high density.The temperature minimum typically lies within the streamer's core, this is also the case for this snapshot due to the cold gas contained within the structure.
Most of the slices exhibit almost pristine metallicity, (/ ⊙ ) ∼ 10 −4 .But high-metallicity gas particles of (/ ⊙ ) ∼ 10 −3 − 10 −1 are found scattered across the cross-section without a strong correlation with the density or temperature.This is in agreement with the metallicity radial profile in Figure 14.At the same time, the region surrounding the structure always has a high metallicity, and this is confirmed by the slice at  ∼ 350 ℎ −1 kpc.
We have measured the radial infall velocity at the innermost radius of  = 50 ℎ −1 kpc, and it has reached  r ∼ 250 − 200 km s −1 for the core gas.But the infall velocity for the envelope drops to as low as  r ∼ 25 − 50 km s −1 .At larger radii, the infall velocity gradually decreases from the core to the envelope from  r ∼ 200 km s −1 to below  r ∼ 50 km s −1 .Such velocity gradients will induce shear between the core and the envelope, and trigger turbulence.Furthermore, the high density zones, i.e., the cores, tend to separate into different infall velocity regions with a spaghetti-like morphology.Both the developed turbulence and the spaghetti-type flow are characteristics of dissolution of the filamentary flow and its gradual virialization.
The lower plot of Figure 14 displays the 0.1 vir thick slice of the same halo, Z2L, with the superposed velocity field.The kinematics of filamentary streamers is delineated here, supplemented by increasingly virialized gas around the central region.The major streamers can be seen partly hitting the central galaxy region and partly missing it converting the inflow into outflow.The emergence of the central Figure 14.Upper plot: the gas particles in the filamentary streamer #1 of the halo Z2L at  = 2.1 (see Figure 1), colored with density (top left), temperature (top right), metallicity (bottom left) and the infall velocity (bottom right).Lower plot: the CGM temperature in a slice with 0.1 vir thickness, normalized by  vir from 0.1 (galaxy's radius, inner circle) to  vir for the halo Z2L.The velocity field is traced by the blue streamlines.turbulent region can be observed as well.Clearly, the gas withing the halo is not in hydrostatic or thermal equilibrium.

Joining the central galaxy: dissolution of the accretion streamers
We proceed to address the evolution of the streamers in the innermost halo region of  ∼ 10 − 50 ℎ −1 kpc, where it is more difficult to separate the streamers from the smooth accretion and where the kinematics becomes more complex -thus characterizing the virialization process of the streamers.This region is the most interesting one because it borders with the HOP-defined central galaxy -we call it the inner CGM in contrast to the overall CGM which extends to the baryonic backsplash radius.
As stated earlier, in order to analyze the streamer's kinematic and thermodynamic properties, we identify their particles at larger radii and follow them individually to the center.For the sake of not  overestimating the inflow rate within this region, we only consider particles belonging to the streamer or to the diffuse flow if their radial velocity exceeds the local dispersion velocity in the gas.
Figure 16 displays the fate of the streamer #1 within the halo Z2L discussed earlier (see also Figure 1), using two scales: a small scale of 30 ℎ −1 kpc (top panels) and a large scale representing the virial radius, ∼ 200 ℎ −1 kpc (bottom panels).On both scales the streamer and the flow streamlines are presented in two perpendicular projections with respect to the central galactic disk, face-on (left column) and edge-on (right column).We noted before that this streamer approaches the galaxy with an angle of ∼ 30 • to the galactic equatorial plane.In fact, in our models, the streamers appear to be inclined to the galaxy planes at various angles, as has been shown in previous  simulations (e.g., Heller et al. 2007;Shlosman 2013;Krolewski et al. 2019).
In the lower panel of Figure 16, the gaseous streamer can be seen maintaining the temperature around  ∼ 10 4 − 10 5 K when falling through the virial radius, and it is heated up dramatically to  ∼ 10 6 K at around  ∼ 50 ℎ −1 kpc.The streamer has a thickness, i.e., diameter, of ∼ 120 ℎ −1 kpc.This is much wider than the galaxy scale of ∼ 20 ℎ −1 kpc.Most of the gas falls towards the center being aligned with the direction of the streamer until  ∼ 30 ℎ −1 kpc, where the gas streamlines become distorted and more turbulent (see also upper panels).Some of the infalling gas which missed the galaxy is converted into an outflow here, which extends to the splash radius -the first apocentric passage, which in most of the cases lies outside  vir .Although originally defined with respect to DM (e.g., Diemer & Kravtsov 2014;Adhikari et al. 2014), it has also a baryonic feature formed by the caustic generated by the piling up of the accreted material near apocenters.
Analysis of a turbulent motion is inherently difficult.We follow the method outlined in Choi et al. (2013), which relies on the vorticity and its cross product with the velocity field, ì  = ▽ × ì .We have pixelized the filament image and plotted the vorticity field in Figure 17 at the galaxy (upper frames) and the halo scales (lower frames) as specified above.On large scales and outside the streamer, the flow appears to be irrotational.Inside the streamer, the vorticity increases closer to the galaxy.This is partly due to the shear discussed above and the resulting Kelvin-Helmholtz instability which triggers ablation of the streamer gas inside  ∼ 30 ℎ −1 kpc.Ablated gas contributes to the turbulent media surrounding the galaxy.Understanding and quantifying this ablation process is of a prime importance for the evolution of streamers at smaller radii.
Figure 18 displays the growth rate of galaxies (dashed lines) in our sample overplotted on their SFRs (solid lines).For  f = 6 models, the growth rate exceeds the corresponding SFR by about an order of magnitude.For  f = 4 the difference between the growth rate and the SFR here decreases towards  f , independently of the environment and feedback.The decrease comes clearly because of the reduced growth rate, which depends on the net accretion rate.

DISCUSSION
We used high-resolution zoom-in cosmological simulations to trace the structure and evolution of the cosmological gas streamers penetrating DM haloes until the galaxy boundaries at selected final redshifts of  f = 6, 4, and 2. All the haloes have been chosen to have similar masses of log  vir /M ⊙ ∼ 11.75±0.05at their final redshifts, and evolving in high or low density environments, i.e., in different overdensities.Furthermore, the resulting central galaxies have been subjected to different types of feedback.We compared the thermody- namic and kinematic properties of streamers and diffuse accretion at these redshifts, starting outside the virial radii and down to the central galaxy regions.We analyze the dissolution process of streamers due to their interaction with the galactic environment extending to the halo virial radius and forming the CGM -the gaseous counterpart of DM halo.
We start by summarizing our results and analyze them.Specifically, we find that, • Using a hybrid d-web/entropy method applied to the gaseous filaments, i.e., streamers, allows us to map these streamers down to ∼ 10 ℎ −1 kpc, i.e., down to the central galaxy scales.Applying this method in tandem with the gas kinematics provides an efficient way to separate the inflow from outflow and from the virialized gas within the parent DM haloes.
• Accretion rates decrease with decreasing final redshifts, whether they proceed via streamers or via diffuse accretion -this is a direct consequence of the reduction of available gas and expansion of the universe.The typical density of both types of accretion declines with redshift as well, i.e., from  ∼ 10 −24 − 10 −25 g cm −3 at  f = 6 to ∼ 3×10 −26 −10 −27 g cm −3 at  f = 2.We do not observe a systematic difference in the accretion rates in CW and VW models, not in the filamentary or diffuse modes.This result needs to be explained, and we address it in this section.
The temperature inside the streamers increases inside the virial radii, and faster at lower redshifts.However, when mass-weighted, it shows that majority of gas remains at low temperature, ∼ 10 4 K, at  f = 6 and 4, and only at  f = 2 it spreads more equally in mass for  ∼ 10 4 − 10 6 K.The temperature of the diffuse mode of accretion is higher than in the streamers inside the haloes, and while most of the gas is at ∼ 10 4 K outside the virial radius, it forms a bi-modal distribution of  ∼ 10 4.5−5 K and ∼ 10 5.5−7 K.In the presence of a wide range of temperatures and velocities there, a wide range of ionization is expected as well.The maximal radial velocities within the streamers decrease with a decreasing redshift, from ∼ 200 − 300 km s −1 at  f = 6, down to ∼ 100 − 150 km s −1 at  f = 2.The infall velocities for the diffuse gas are smaller compared to the streamers, due to the difference in the thermodynamic state of the diffuse and filamentary gas, as we have discussed in section 3.3.
• Outside ∼ 2 vir , both the streamers and the diffuse gas accretion have metallicity in the range of / ⊙ ∼ 10 −4 − 10 0 , while inside the haloes, the metallicity distributions are flat.The VW models appear less metal rich than CW models, due to the larger gas fraction in VW galaxies and lower SFRs.The radial inflow velocities inside the central ∼ 100ℎ −1 kpc are typically lower for VW models than CW, possibly due to the stronger shocks in the former.
• Within the DM haloes, the streamers exhibit a core-envelope structure: the core radial flow surrounded by a lower density envelope.The core possesses an elevated density and lower temperature, a similar metallicity to that in the envelope, but a higher inflow velocity.This velocity gradient between the core and its envelope induces shear.At smaller radii, inside ∼ 50ℎ −1 kpc, the core-envelope structure is fading away and the filamentary flow splits and can be described as a spaghetti-type flow.Within the central ∼ 30 ℎ −1 kpc of the streamers, the shear triggers the Kelvin-Helmholtz instability and turbulence which contribute to the ablation process of the streamers and dissolves them.We have quantified this turbulence and mapped it in two projections.
• Galactic outflows in all our models reach the characteristic distance of ∼ 0.5 vir ∼ 100ℎ −1 kpc from the central galaxies.The maximal outflow rates are located close to those galaxies, ∼ few× ⊙ yr −1 , which is less than 10% of the inflow rate at  f = 6 and 4.But for  f = 2 models, the accretion rates have declined and the maxima of the galactic outflows have increased.Hence they differ by factor of ∼ 2 only.Moreover, the outflow rates decline substantially with the distance to the galaxy, to < ∼ 0.1  ⊙ yr −1 at ∼  vir .• The entire DM halo region has a complex kinematics involving filamentary and diffuse accretion, and galactic outflows.Some of the filamentary accretion is deflected and extends back to the baryonic backsplash radius, ∼  vir − 2 vir .Based on the modeled properties of such a multiphase CGM, one must conclude that it cannot be in thermal equilibrium, and its dynamic and thermodynamic state are strongly time-dependent.
• We do not find any significant dependence on the environment in the filamentary and diffuse accretion in our models.This is probably related to a relatively small difference between the overdensities involved, between  ∼ 1.3 and  ∼ 3 (Table 1).
In the cold accretion model framework, the filamentary accreting gas can be kept at low temperature within the DM haloes, without being shocked around the virial radius, in constrast with the diffuse accretion in more massive haloes.We have confirmed that the streamer gas temperature remains low, ∼ 10 4 K, as it crosses the virial radius, which is substantially lower than the halo virial temperature of ∼ 4 − 10 × 10 6 K, depending on  f .Inside the virial radius, some of the gas temperature starts to increase, reaching ∼ 10 5.5−6 K.The accreting gas is heated up by the galactic winds, is compressed and heated by the decay of the turbulent motions and small-scale shocks.However, the amount of the hot gas is not large, the majority of the filamentary flow remains cold, as can be seen in the Figure 4.Only at  f = 2, the cold stream appears to be disrupted inside  vir , and the filamentary gas is distributed evenly in the range of 10 4 − 10 6 K.
The properties of diffuse gas differ substantially from that of the filamentary one.Outside  vir , most of the gas remains around 10 4 K, but inside the halo, the gas heats up and its mass is distributed between 10 5 K and up to the virial temperature.
Figures 7 and 13 reveal that about equal mass accretion rates of filamentary and diffuse accretion cross the virial radii.But deeper in the haloes, the diffuse accretion rates start to dominate, as the filamentary flow dissolves gradually.Therefore, the hot gas starts to dominate in mass in the inner halo.
Fielding et al. ( 2020) have compared in detail the CGM properties in various numerical simulations at  = 0, invoking the pressureentropy plane.Following their work, we have attempted to characterize the properties of the gas filling up the CGM in our high resolution simulations into the pressure-entropy plane, separating the inner and outer CGM of three ZLCW haloes (i.e., Z2LCW, Z4LCW and Z6LCW), at  f = 6, 4 and 2 (Figure 20).This gas excludes the ISM of the central galaxy determined by the HOP.Moreover, we limit our analysis to CW models residing in the low overdensity regions.
We find that both CGM regions differ substantially in their properties.Generally, we distinguish three different regions which account for most of the gas mass in Figure 20.First, the top region, which extends from / vir ∼ 10 1 to ∼ 10 −2 .Second, the middle region of an isothermal gas at ∼ 10 4 K at / vir ∼ 10 −2 − 10 −4.5 .And third, the horizontally extended and, therefore, radially extended, low-entropy gas across a range in density at / vir ∼ 10 −4.5 − 10 −6 .
We can reconstruct the locations and masses of this gas in this halo.The total mass of the CGM is ∼ 10 10  ⊙ at  f = 6, ∼ 3 × 10 9  ⊙ at  f = 4, and ∼ 2 × 10 9  ⊙ at  f = 2.The first region defined above contains ∼ 3 × 10 9  ⊙ (inner CGM) and ∼ 4 × 10 8  ⊙ (outer CGM) of the virialized hot diffuse gas at  f = 6.This amount decreases by a factor of 3 (inner CGM) and stays the same (outer CGM) at  f = 4.It stays unchanged further on at  f = 2.
The second region has ∼ 4 × 10 9  ⊙ (inner CGM) and ∼ 2 × 10 9  ⊙ (outer CGM) of the nearly isothermal gas at  f = 6.This gas is in thermal equilibrium with the cosmological UV background.Its amount decreases by a factor of 4 (inner CGM) and stays the same (outer CGM) at  f = 4.It further decreases by a factor of 2 (inner CGM) and by a factor of 10 (outer CGM) at  f = 2.
The third region includes the starforming gas locked in the substructures.Its left boundary is limited by our critical density for star formation, 4 cm −3 .It extends horizontally to high pressure, and this extension depends on the star formation recipe -it can be seen in the log  − log  diagrams as well.Note also that because this gas is located in substructures, it should be normalized by the virial pressure of each object, but we ignore this change in normalization.The starforming gas mass has ∼ 5 × 10 8  ⊙ (inner CGM) and ∼ 4 × 10 7  ⊙ (outer CGM) at  f = 6.This amount decreases by a factor of 10 (inner CGM) and by a factor of 100 (outer CGM) at  f = 4.It further decreases by a factor of 100 (inner CGM) and by another factor of 10 (outer CGM) at  f = 2.
When the starforming gas exists, it has the highest pressure of all three components discussed above.The high entropy gas which is always present in simulations extends from low to intermediate pressure, while the isothermal gas exhibits the lowest pressure and intermediate entropy.
Overall, we observe similar structural details in pressure-entropy plane, as well as substantial differences between the gas properties in the inner and outer CGM, as well as redshift evolution in similar DM haloes.In the inner CGM, the amount of the starforming gas decreases abruptly with  f , and this gas is basically absent at  f = 2.The amount of 10 4 K gas decreases as well with decreasing redshift, but not so dramatically.In the outer CGM, the starforming gas disappears already by  f = 4.And the amount of 10 4 K gas decreases sharply as well, becoming negligible at  f = 2.In summary, the CGM analysis exposes its multiphase character and being away from thermal and dynamic equilibrium.
Note that masses quoted above represent the CGM in specific and similar haloes used in our simulations.The median properties of the CGM, such as mass fractions, have been found to differ among numerical simulations due to the galactic feedback (e.g., Davies et al. 2020).Such differences in the CGM become explicit in our nearly identical haloes at three different redshifts due to the changing galactic feedback.
Our results presented in sections 3.2 and 3.3 show that filamentary and diffuse accretion rates appear similar in CW and VW models outside the radius of ∼ 50ℎ −1 kpc.This appears puzzling as the VW feedback is about 8 times stronger than the CW one (e.g., section 2.3).However, one must remember that the galactic wind feedback is supplemented by the SN feedback which is strongly related to the star formation rate.
In order to demonstrate the relative effects of the wind and SN feedback, we have calculated the instantaneous energies of the wind and the SN energy depositions at various times.For example, at  f = 2, for the galaxies Z2LCW and Z2LVW, the ratio of the kinetic energies of the VW to CW wind feedback is ∼ 14.Hence, the wind energy VW is substantially higher than in the CW at this time.Next, we have estimated the energy contribution from the SN during the preceding 1 Gyr in each model, and found than the ratio for VW to CW energies is ∼ 0.38, much smaller for the VW model.This is understandable, as the efficiency of the SF in VW models is substantially lower than in CW, resulting in a much larger fraction of gas in galaxies (Bi et al. 2022a).Combining the contribution of the wind and the SN rate in both models, we find that the total energy feedback in VW model is only larger by a factor of ∼ 2.6.Hence, we conclude that the absence of a real difference between the mass accretion rates in VW and CW models is due to the relatively small net difference in the applied feedback.
One should be careful not to make far reaching statements about the state of the accreting gas across the virial radius based on numerical simulations.For example, the AGN feedback, which is absent in the current simulations, has been shown to have a strong effect on the state of the halo gas and hence its accretion rate onto the central galaxy (e.g., Woods et al. 2014;Übler et al. 2014;Nelson et al. 2015).We, therefore, limit our conclusions about the state of the CGM which follow from the types of feedback used here.
While the inflow velocity profiles along the streamers increase at smaller radii, within ∼ 100ℎ −1 kpc the trend is reversed.The maxima of this velocity decrease with lower  f (Figure 6).This decrease appears to correlate with the halo escape speed which decreases with time because the halo concentration of similar haloes at different redshifts decreases with lower  f .The velocity profiles of diffuse accreting gas are flatter (Figure 12), but generally follow the same trend as the filamentary accretion.Slicing the streamers perpendicular to their spines (e.g., Figure 15) allows us to define the high density core and low density envelope region(s) in each slice, e.g., Figures 14 and 15.The high and low density regions in the streamers display velocity gradients, i.e., the core region separates higher from lower infall velocities.Such gradients are expected to contribute to the ablation process of the core region in streamers by the Kelvin-Helmholtz instability.Indeed, at smaller radii, we observe the splitting in the streamer's core, and the streamer starts to resemble a spaghetti-type flow.
Along with the gas temperature and density radial profiles, we do find substantial gradients for the metallicity on the halo scale, i.e., between the inner and outer CGM.But on smaller scales, larger variations of metallicity prevail, forming pockets of high metallicty.When a (sub)structure is embedded in a streamer or in diffuse accretion, its metallicity diffuses out either by stellar feedback or other processes, forming these pockets.
As one of the important issues in understanding the evolution of the filamentary inflow inside the virial radius, we have measured the ratio of the cold-to-hot gas mass,  ch =  cold / hot , at  vir , 0.5 vir and 0.1 vir (not shown here).We define the cold gas with  < ∼ 3 × 10 4 K, and hot gas with the temperature above this threshold.We find that  ch declines monotonically with lower  f from ∼ 2 − 3 at  f = 6 to ∼ 0.25 − 0.5 at  f = 2 by an order of magnitude.The sharper decline happens inside 0.5 vir .Not accidentally this decline coincides with the region where the galactic outflows are strong, i.e., inside 100ℎ −1 kpc.At all radii,  ch falls below unity at  f = 2, so, a sharp decrease at lower redshift haloes.Note that our haloes have the same mass by design at all these final redshifts.
Finally, we analyse the structure of the CGM which populates the entire DM haloes outside the central galaxies.As an example, we present the evolution of the CGM around the Z6LCW galaxy -its mass growth, and temperature and metallicity evolution (Figure 19).We have divided the CGM which occupies the parent DM halo, into inner and outer CGM.Based on the ratio of cold-to-hot gas,  ch , discussed above, the inner CGM can be defined inside the inner 100ℎ −1 kpc, and the outer CGM region between this radius and the virial radius.The inner CGM is also affected by the galactic feedback stronger than the outer CGM.For models presented here, we find that the CGM is not in equilibrium, nor kinematically, thermodynamically, chemically or temporally.It is constantly perturbed by the influx mass, momentum and energy across  vir and from the central galaxy, and forms a multi-phase structure.This conclusion agrees with observations of low-redshift galaxies (e.g., Tomlinson et al. 2017, and refs. therein).
Furthermore, the filamentary flow which does not impact the central galaxy expands to the baryonic backsplash radius, which appears to lie at ∼ 1 − 2 vir , smaller than the DM backsplash radius, and, therefore, can provide additional perturbations to the CGM, especially for  > ∼ 1, when the filamentary accretion is more significant (e.g., Figure 14, lower panel).The shock associated with the splash radius is visible in our analysis but is not shown here.This result agrees well with previous determinations of the shock position (e.g., Nelson et al. 2016).
Recent work has raised an interesting issue that the results of zoom-in simulations can be affected to some extent by the 'butterfly effect'2 by introducing stochasticity based on the usage of random number generators implemented in various recipes, such as star formation routines, etc. (e.g., Genel et al. 2019;Ramesh & Nelson 2023).To what extent this affects our results?
We have encountered this effect as well as the 'timing effect' (Springel et al. 2008) in another set of the zoom-in numerical simulations (Goddard et al., in prep.).We found that the evolution initially separates, but then runs in parallel, with various events happening with a small time delay.However, we did not find that the evolution diverges as in chaotic systems.Hence, in the current set of simulations, even if both effects are present, they have only a small amplitude and do not grow with time.For example, the actual differences between the CW and VW models appear not to be affected by them.

Streamer properties in the innermost halo
Turning our attention to the innermost haloes, inside ∼ 30 ℎ −1 kpc, we observe that the geometrical width of the streamers exceeds substantially the central galaxy cross sections.In all our models, the streamers are inclined to the equatorial plane of the central galaxy.The inclination angle is never around 90 • or 0 • .For example, Figure 16, shows the inclination angle being ∼ 30 • for the streamer #1 (see also Figure 1 for the identification of this filament).
Inside this region, the streamers have been observed to dissolve gradually by the ablation process (e.g., Figure 16).This process is characterized by an increase in the turbulent layer associated with the streamer and decreasing radii.In order to quantify the turbulence, we have calculated the associated vorticity (see section 3.5) and have pixelized it (Figure 17).It shows that the dissolution is rapid and the turbulence, which is minimal at higher radii, is spreading sideways.A related question is to what extent this turbulent flow in the innermost halo region contributes to the turbulence in the underlying galactic disk (e.g., Choi et al. 2013) We can recognize few types of interactions between the baryonic filamentary inflow and the galactic disk.First, some of the inflow impacts the galaxy and has both prograde or retrograde motion with respect to the rotation of the stellar/gaseous galactic disk.We expect, that the prograde encounter between the filamentary inflow and the galactic disk is smoother than the retrograde one.Second, and more generally, most of the gas inflow misses the direct encounter with the disk and creates a 'tail' flow around the disk, which appears to be turbulent.This part of the inflow which has missed the disk contributes to buildup of the extended and highly turbulent region around the disk, and its kinematics differs substantially from that of the filamentary or diffuse accretion at larger radii.
We observe that all the halos in our simulations are connected typically by 2-3 filaments, and only in one case by 4 filaments.This is in agreement with other simulations (e.g.IllustrisTNG simulation, Nelson et al. (2018)).The reason for this is still not well understood.
Comparison between our profiles of thermodynamic variables provide a complementary information to those presented in Fielding et al. (2020).Direct comparison for our median profiles agree with the latter.But our hex-binned distributions of the temperatures in filamentary and diffuse accretion separate their gas mass as a function of temperature.For example, the increasing mass accumulation around  ∼ 10 4 K in the inner halo is evident in the filamentary gas (Figure 4), but not in the diffuse accretion (Figure 10).
The difference in growth rate with redshift in our galaxies (Figure 18) is also mainly due to the time employed by the galaxies to reach their final masses, i.e., for a given mass range, galaxies at high- grow faster than galaxies at lower-.

From the virial radius to the central galaxy
Finally, we take a look at the big picture of accretion and outflows inside the parent haloes of central galaxies and its immediate environment, i.e., within ∼ 2 vir .While the number of models presented here is limited compared to the cosmological simulations in large computational boxes, the advantage of our zoom-in simulations is the higher resolution which allows to obtain a detailed picture of kinematic and thermodynamic properties of these flows.
We find that the gaseous component of DM haloes in our models has a complex structure, both thermodynamically and kinematically.Filamentary and diffuse inflows differ in density, temperature, metallicity, and in the associated velocity field.Moreover, the observed outflows have multiple sources of origin.First come the galactic winds of various strength.Second contribution comes from filamentary accretion which missed the central galaxy and has been diverted radially out.The cross sections of the streamers are typically much wider than those of the galaxies.Much of the flow along the streamers, therefore, misses the galaxy, but turns around and can even interact with itself.
The outer boundary of this interaction defines the backsplash radius (e.g., Adhikari et al. 2014).The backsplash radius of the gas is expected to be smaller than the one for the DM due to the shocks and associated dissipation.By using the logarithmic derivative of the density with respect to the radius, we have checked for the position of the backsplash radii for our central galaxies.For  f = 6 and 4, we find a good corresponce with  200 , which lies close to  vir .But for  f = 2, we find that the backsplash radius lies further out of  vir , on the average around ∼ 1.7 vir .We conclude that the shock associated with the baryonic backsplash radius in our simulations lies at ∼ 1 − 2 vir .
We find that the properties of the filamentary accretion vary with distance to the central galaxy.We detect the ablation process of the filamentary inflow, associated with Kelvin-Helmholtz instability and generation of turbulence around the dissolving filaments, leading to the spaghetti-type flow.Thus a filament separates in a number of individual streamers which continue to be ablated and mixing with the CGM.Hence, we associate the CGM with the entire volume of the DM halo surrounding the central galaxy and define it as the transition region between the galactic ISM and the virial radius.This confirms that the CGM is highly inhomogeneous, i.e., multiphase, in its thermodynamic and kinematic properties.
We have compared our analysis of the filamentary accretion with that performed by Ramsoy et al. (2021) at high redshifts using the adaptive mesh refinement (AMR) code ramses, with a comparable resolution to ours.Both approaches have detected a developed turbulence measured by vorticity in the filaments.And in both cases, this vorticity has delineated the boundary of a filament, allowing to determine its radius, i.e., its cross section.While Ramsoy et al. have focused on the typical filament structure, we have addressed the dissolution process of the filament.The core-envelope structure of a typical filament in our simulations have been found to separate in a spaghetti-type flow (section 3.4).The Kelvin-Helmholtz instability further acted to ablate the streamers, ultimately mixing it with the diffuse environment.
We have also calculated the fraction of volume occupied by the walls shown in Figure 1, and found it is roughly independent of the final redshifts,  f .By definition, the sheets represent a onedimensional stable manifold, while a filamentary streamer represents a two-dimensional stable manifold.Note, that Hahn et al. (2007a) has measured evolution of this fraction and found that it increases untill  ∼ 2 and decreases thereafter.However, there is no contradiction with our results -our redshift dependence is not an evolutionary one, but compares the environment of similar mass haloes at different final redshifts.When we switch to evolutionary redshift for each halo separately, we confirm the Hahn et al. result.All the streamers exhibit metallicity in the range of / ⊙ ∼ 10 −3 − 10 −1 , with a radial gradient between the outer and inner halo.Inside  vir , the streamers have a mass-weighted metallicity / ⊙ > ∼ few × 10 −2 − few × 10 −1 , and agrees well with observations in the optical-X-ray bands (e.g., Werk et al. 2014;Lehner et al. 2019).For the diffuse accretion, the metallicity is smaller by a factor of a few.

CONCLUSIONS
Based on a set of high-resolution zoom-in cosmological simulations, we investigate the baryonic filaments (streamers) and the diffuse accretion flows which channel the gas across the virial radii and down to the galactic regions, under a dual action of accretion and galactic outflows.We choose a set of DM haloes with similar masses of log ( vir /M ⊙ ) ∼ 11.65 at three representative redshifts,  f = 6, 4 and 2 from Bi et al. (2022a).These haloes have been evolved in relatively low and high overdensities compared to the average density in the universe, and using two different stellar feedback.
Using a hybrid d-web/entropy method supplemented by kinematics in the innermost regions, we have mapped the filamentary streamers and seprated them from galactic outflows and diffuse accretion flows.This allowed us to analyze the dynamic and thermodynamic properties of the CGM in nearly identical haloes at different redshifts.We find that the CGM is highly inhomogeneous and multiphase, and not in thermodynamic or dynamic equilibrium.
We find that accretion rates in filamentary streamers exhibit decrease with decreasing redshifts, and the inflow velocities along these filaments decrease by a factor of ∼ 2 with lower , again in similar haloes.The temperature inside the CGM increases at smaller radii, as well as with decreasing redshit.
The filamentary streamers display a core-envelope structure inside the virial radius -a higher density, lower temperature core surrounded by a lower density, higher temperature envelope.The filament separates into spaghetti-type flow.Inside the inner ∼ 30 ℎ −1 kpc, we show that the filaments develop the Kelvin-Helmholtz instability, which triggers turbulence, ablates and dissolve them.
We find that the galactic outflows in tandem with diverted accretion flow affect the accretion flow, mostly within the inner CGM of ∼ 0.5 vir ∼ 100ℎ −1 kpc.Finally, we find that the thermodynamic properties of the CGM gas can be separated into three phases on the pressure-entropy plane (first used by Fielding et al. (2020)).These regions include the high entropy, low-density hot gas, isothermal gas at ∼ 10 4 K in equilibrium with the UV background, and the low entropy starforming gas.

Figure 1 .
Figure 1.Structure classification of the DM-haloes environment using the dweb method applied to baryons in boxes of length 1, 200 h −1 kpc, (∼ ±3 vir ).The color coding represents: knots in blue, filaments, (i.e., streamers) red, walls green and voids in black.From left to right columns, haloes at their final redshifts  f = 6, 4, and 2 are shown.All the models correspond to the CW feedback as they are similar to the VW models at these scales.The upper panels display haloes in the high overdensity environments and the lower panels show haloes in low overdensity environments, (see Table1).The arrowhead at the Z2L snapshot points to the filament #1 which is analyzed in detail in the subsequent sections.

Figure 2 .
Figure 2. Inflowing streamers, diffuse accretion, and walls (seen as threads connecting the streamers) at three representative radii in a thick shell of 10ℎ −1 kpc, at 0.1 vir , 0.5 vir and 1 vir for haloes at  f = 6,  f = 4 and  f = 2. Shown as whole-sky HEALPix projection maps (see the text).The color represents radial mass influx of gas per solid angle.Thin walls are seen between the streamers at larger radii until dissolved on the galaxy scale of ∼ 0.1 vir , especially at  f = 6 and  f = 4.

Figure 3 .
Figure3.Volume density radial profiles of gas representing the sum of all the streamer spines at  f = 6, 4, and 2, from left to right columns, respectively.The streamers have been divided into shells of 1 ℎ −1 kpc thickness.The corresponding DM haloes can be identified from Figure1.The VW feedback models are given by the red lines and CW models by the blue lines, which represent the median.The color shadows show the 20-80 percentile scale.The upper panels display haloes in the high overdensity environment and the lower panels show haloes in the low overdensity environment (see Table1for further details).

Figure 4 .
Figure 4.The hex-binned gas temperature radial profiles of the filamentary streamers within 2 vir at  f = 6, 4 and 2, with the color palette representing the gas mass in individual pixels.Note that the virial temperatures of these DM haloes is  ∼ 1.3 × 10 6 (1 + ) K.

Figure 5 .
Figure 5.As in Figure 3, but for the mass weighted average gas metallicity radial profiles of filamentary streamers.

Figure 9 .
Figure 9. Radial density profiles of diffuse gas.The VW (CW) models are given by the red (blue) lines.The color shadows represent the 20-80 percentile distributions.

Figure 10 .
Figure10.Hex-binned radial temperature profiles of a diffuse gas inside 2 vir with the color palette representing the gas mass in individual pixels at  f = 6, 4 and 2. Note that the virial temperatures of these DM haloes is  ∼ 1.3 × 10 6 (1 + ) K.

Figure 11 .
Figure11.Mass-weighted gas metallicity radial profiles of diffuse gas.The VW models are given by the red lines, and CW models by the blue lines.

Figure 12 .
Figure 12.Radial velocity profiles of diffuse gas.The VW (CW) models are given by the red (blue) lines.The color shadows show the 20-80 percentile scale.

Figure 13 .
Figure13.Mass accretion rate profiles of diffuse gas (solid lines) overplotted onto the outflows (dashed lines) triggered by feedback from the central galaxies.The VW (CW) models are given by the solid red (blue) lines.For both inflow and outflow rates, we have used the radial velocities exceeding the local dispersion velocities in the gas.

Figure 15 .
Figure 15.Representative slices of the gas in the filament #1 of halo Z2L (see Figure 1) at  f = 2. Shown the gas density (top) and metallicity (bottom).From left to right the slices represent the cuts at distances  = 50 (the inner CGM), 150 (the outer CGM), 200, 250, and 350 ℎ −1 kpc.Each slice has a 10 ℎ −1 kpc radial thickness.Each frame has 120 ℎ −1 kpc length on the side.

Figure 16 .
Figure16.Face-on (left) and edge-on (right) maps of the gas stream in filament #1 of the halo Z2L at  f = 2 (see Figure1), colored with the gas temperature.The upper panels represent a region of the 60 ℎ −1 kpc side, and the lower panels of 400 ℎ −1 kpc side.The dashed circles represent the virial radii.The central black contours in the upper boxes display the shape of the central galaxy, i.e., face-on and edge-on.

Figure 17 .
Figure 17.Face-on (left) and edge-on (right) vorticity maps of the gas particles in the filament #1 of the halo Z2L at  f = 2 (see Figure 1).The upper panels shown in the 60 ℎ −1 kpc boxes and the lower panel shown in the 400 ℎ −1 kpc boxes.The dashed circles represent the virial radii.The central white contours in the upper frames display the shape of the central galaxy, i.e., face-on (left) and edge-on (right).

Figure 18 .
Figure18.Evolution of the SFR (solid lines) and galaxy growth rates (dashed lines) in all models of our simulation suite.The CW and VW models are represented by blue and red lines, respectively.The galaxy growth rate (gas + stars) corresponds to the baryonic accretion rate minus the outflow rates.The galaxies residing in a high density environment are shown in the upper panels and those in low density environment lie in the lower panels.

Figure 19 .
Figure 19.Evolution of the CGM of the Z6LCW halo model.Left: the gas mass in the inner (red line) and outer (blue line) CGM, which envelops the central galaxy and extends to ∼  vir .The inner CGM is defined within 0.5 vir , and the outer CGM lies in the range of 0.5 vir −  vir ; Center: the median temperature evolution of the inner (red) and outer (blue) CGM.The shaddows represent the 20-80 percentiles.Right: average metallicity evolution of the inner (red) and outer (blue) CGM.

Figure 20 .
Figure 20.Pressure-entropy phase diagram in 0.1 vir thick shells in the inner (centered on 0.25 vir ) and outer (centered on 0.75 vir ) CGM of three haloes LCW at  f = 6, 4 and 2. The dashed lines show lines of a constant temperature (10 3 − 10 7 K) and constant density (10 −3 − 10 1 cm −3 ).Note that the color palette is given in  ⊙ .

Table 1 .
Table of the halo properties in the DM only simulations.All values are given at the final redshifts f = 6, 4, and 2. The columns correspond to (from left to right): the final redshift  f ; the model number (see definition in section 2.2; the virial mass of DM halo  vir at  f ; the halo virial radius (in comoving coordinates)  vir ; δ -the local overdensity.
Table 1 for further details).