The Structure and Dynamics of Massive High-𝑧 Cosmic-Web Filaments: Three Radial Zones in Filament Cross-Sections

We analyse the internal structure and dynamics of cosmic-web filaments that connect massive high-𝑧 haloes. Our analysis is based on a high-resolution AREPO cosmological simulation zooming-in on a volume encompassing three Mpc-scale filaments feeding three massive haloes of ∼ 10 12 M ⊙ at 𝑧 ∼ 4, embedded in a large-scale sheet. Each filament is surrounded by a cylindrical accretion shock of radius 𝑟 shock ∼ 50 kpc. The post-shock gas is in virial equilibrium with the potential well set by an isothermal dark-matter filament. The filament line-mass is ∼ 9 × 10 8 M ⊙ kpc − 1 , the gas fraction within 𝑟 shock is the universal baryon fraction, and the virial temperature is ∼ 7 × 10 5 K. These all match expectations from analytical models for filament properties as a function of halo-mass and redshift. The filament cross-section has three radial zones. In the outer “thermal” ( T ) zone, 𝑟 ≥ 0 . 65 𝑟 shock , inward gravity and ram-pressure forces are over-balanced by outwards thermal pressure forces, decelerating the inflowing gas expanding the shock outward. In the intermediate “vortex” ( V ) zone, 0 . 25 ≤ 𝑟 / 𝑟 shock ≤ 0 . 65, the velocity field is dominated by a quadrupolar vortex structure due to offset inflow along the sheet through the post-shock gas. The outwards force is dominated by centrifugal forces associated with these vortices, with additional contributions from global rotation and thermal pressure. The shear and turbulent forces associated with the vortices act inward. The inner “stream” ( S ) zone, 𝑟 < 0 . 25 𝑟 shock , is a dense isothermal core, 𝑇 ∼ 3 × 10 4 K and 𝑛 H ∼ 0 . 01 cm − 3 , defining the cold streams that feed galaxies. The core is formed by an isobaric cooling flow and is associated with a decrease in outwards forces, though it exhibits both inflows and outflows.


INTRODUCTION
In the last two decades, our understanding of galaxy formation in a cosmological context has undergone a paradigm shift.Galaxies, as we now know, do not form in isolation, but rather are connected through the "cosmic-web", which dominates the matter distribution on Mpc scales in both dark matter (DM) and baryons.This intricate structure has been predicted theoretically (Zel'dovich 1970;Bardeen et al. 1986;Bond et al. 1996;Springel et al. 2005), and can be seen in the distribution of galaxies in observational surveys such as 2dFGRS (Colless et al. 2001), SDSS (Tegmark et al. 2004), and the 2MASS redshift survey (Huchra et al. 2005).Perhaps the most prominent ★ E-mail: yul232@ucsd.edufeature of the cosmic-web is a network of intergalactic filaments, which contain roughly half of the mass in the Universe (Zel'dovich 1970;Bardeen et al. 1986;Cautun et al. 2014;Eckert et al. 2015;Libeskind et al. 2018).
The filaments connect the most massive galaxies at any given epoch, which are located at the nodes of the cosmic-web, while more typical galaxies are located along the filaments (see Fig. 1 below).At high redshift,  ∼ (2 − 6) near the peak of cosmic star-and galaxyformation, these filaments manifest as streams of cold, dense gas ( ∼ 10 4 K,  ∼ 10 −2 cm −3 ), predicted to be the main mode of accretion onto high- galaxies, feeding them directly from the cosmic-web.In massive haloes with virial mass  v > ∼ 10 12 M ⊙ above the critical mass for the formation of a stable virial shock (Rees & Ostriker 1977; White & Rees 1978;Birnboim & Dekel 2003;Fielding et al. 2017;Stern et al. 2021b), the streams are predicted to penetrate the hot circumgalactic medium (CGM), maintaining temperatures of  s > ∼ 10 4 K and reaching the central galaxies in roughly a virial crossing time (Dekel & Birnboim 2006;Dekel et al. 2009a).This gas supply sustains the large observed star formation rates (SFRs) of ∼ 100 M ⊙ yr −1 without galaxy mergers (e.g., Förster Schreiber et al. 2006, 2009;Genzel et al. 2006Genzel et al. , 2008;;Elmegreen et al. 2007;Shapiro et al. 2008;Stark et al. 2008;Wisnioski et al. 2015).
Despite their prominence in the large-scale structure of the cosmicweb, and their clear importance for galaxy formation, the structure and properties of filaments are incredibly poorly constrained.Even a question as fundamental as what prevents the filaments from collapsing under their own self-gravity remains open.Previous work has suggested that this could be due to thermal pressure gradients (Harford & Hamilton 2011;Klar & Mücket 2012;Ramsøy et al. 2021), coherent filament rotation (Birnboim et al. 2016;Mandelker et al. 2018;Wang et al. 2021;Xia et al. 2021), helical vortices (Codis et al. 2012;Laigle et al. 2015), turbulence (Mandelker et al. 2018) or magnetic pressure.However, all of these works relied either on highly simplified model assumptions or on low resolution simulations that did not resolve the internal structure of intergalactic filaments.Works that have measured filament properties in large-volume cosmological simulations (e.g., Cautun et al. 2014;Laigle et al. 2015;Libeskind et al. 2018;Ganeshaiah Veena et al. 2018, 2019, 2021;Uhlemann et al. 2020;Song et al. 2021;Galárraga-Espinosa et al. 2021, 2022) commonly involve DM-only N-body simulations or low-resolution hydrodynamic simulations that do not resolve the inner structure of filaments or their support against gravity.While cosmological 'zoom-in' simulations achieve much higher resolution in individual galaxies, and recently also in the CGM (van de Voort et al. 2019;Hummels et al. 2019;Peeples et al. 2019;Suresh et al. 2019;Bennett & Sĳacki 2020), cosmic-web filaments in the IGM typically lie outside the high-resolution region and are still poorly resolved.Inferring detailed filament properties from observations is extremely challenging, particularly in the IGM, owing to their low density, narrow size, and uncertain chemical and ionisation compositions.The properties of filaments as a function of halo mass, redshift, and environment thus remain poorly constrained.
A recent effort to study the properties of intergalactic filaments at high- using cosmological zoom-in simulations with relatively high resolution was made by Ramsøy et al. (2021).These authors used a cosmological zoom-in simulation of a Milky Way (MW)-mass halo,  v ∼ 5 × 10 11 M ⊙ at  = 0, taken from the NUT suite (Powell et al. 2011) and run with the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002).They studied the properties of the main filament that feeds the central halo in the redshift range  ∼ 3.5 − 8, focusing primarily on  ∼ 4 when the filament extended < ∼ 200 kpc on either side of the halo and did not seem to connect to any other massive object.The spatial resolution of the analysed filament was ∼ 1.2 kpc, with small patches around haloes embedded in the filament reaching a resolution of ∼ 0.6 kpc.These authors found the filament to be well described by a model of an isothermal, self-gravitating, infinite cylinder (Ostriker 1964), embedded inside an isothermal selfgravitating sheet that dominated the mass distribution at  > ∼ (15 − 20) kpc from the filament axis.This was found to be true for both the gas and the dark matter, whose density profiles had very similar shapes and widths.An accretion shock was identified around the filament by a sharp increase in gas temperature.However, the postshock gas was found to cool rapidly, and the shock width was limited to a few kpc.The azimuthally-averaged radial profile of the gas temperature increased by <  ∼ 60% between the central value and the peak of the shock.While significant vorticity was identified within the filament, consistent with previous results (Laigle et al. 2015), the filament was found to be primarily supported by thermal pressure.
The analysis and results of Ramsøy et al. (2021) are novel and extremely detailed.However, since their work focused on a single filament feeding an isolated high- progenitor of a MW-mass halo, it is unclear whether their results describe filaments feeding massive,  v > ∼ 10 12 M ⊙ , haloes as cold streams at high-, as these are much more massive than the studied system.Additional studies of filaments feeding different mass haloes, in different environments, and at different redshifts are needed to draw more general conclusions about the filament population.Finally, while the resolution of ∼ 1 kpc in intergalactic filaments is quite good and better than most state-of-the-art simulations, it is unclear whether this is enough to resolve the detailed internal dynamics of filaments given the filament radius of ∼ 15 kpc, and in particular the support due to turbulence and/or vorticity (for discussions in the context of support of the CGM see, e.g., Bennett & Sĳacki 2020;Lochhaas et al. 2023).
In this work, we seek to study filaments around more massive haloes, in less isolated regions, and at higher resolution.We use a novel suite of cosmological simulations first introduced in Mandelker et al. (2019b) and Mandelker et al. (2021), which zoom in on a large patch of the IGM in between two massive haloes, which at  ∼ 2 are  v ∼ 5 × 10 12 M ⊙ each and connected by a ∼ Mpc-scale cosmic filament.We hereafter refer to this suite of simulations as IPMSim, where "IPM" refers to the "intra-pancake medium" of multiphase gas within cosmological sheets (or "Zel'dovich pancakes," see Mandelker et al. 2021;Pasha et al. 2023).The large filament at  ∼ 2 was formed by the merger of several smaller filaments at  ∼ 2.5.Our analysis focuses on  ∼ 4, when the system contained several welldefined filaments that connect three haloes with  v ∼ 10 12 M ⊙ each (see Fig. 1 below).This simulation employs the resolution characteristic of a standard "zoom-in" simulation of a single halo within the entire IGM region between the two protoclusters.In the highest resolution version of the simulation used in this work (see Section 2.1), the typical cell size in the IGM at  ∼ 4 is ∼ 1 kpc near the midplane of the sheet and ∼ 300 pc near the cores of the filaments.A detailed convergence study presented in Mandelker et al. (2021) suggests that this appears sufficient to resolve the formation of multiphase gas in turbulent IPM.However, the morphology and distribution of the cold component are likely only marginally resolved.As shown in Section 3.1 below, the characteristic scale of the accretion shock that surrounds our filaments is  shock ∼ 50 kpc, while the radius of the cold stream is  stream ∼ 0.25 shock ∼ 12 kpc.Thus, we have ∼ 80 cells across the diameter of the cold stream, and many more across the diameter of the entire filament, sufficient to resolve the driving scale of turbulence and a partial inertial range, thus resolving most of the support due to turbulence and vorticity, if present.This work is organised as follows: In Section 2, we introduce the simulation, our method for selecting filament slices for analysis, and our division of the filament cross-section into three distinct radial zones.In Section 3, we analyse the radial structure of gas and dark matter in filaments, including their thermal structure.In Section 4, we discuss the thermal stability of gaseous filaments by comparing their cooling and free-fall time scales.In Section 5, we address whether the gaseous filaments are in virial equilibrium within the potential well set by the dark matter filament.In Section 6, we decompose the forces that support gaseous filaments against their self-gravity and use a simple model to explain their origin.We discuss the implications of our results and compare them to previous work in Section 7. Finally, in Section 8, we conclude.Throughout, we assume a flat ΛCDM cosmology with Ω m = 1 − Ω Λ = 0.3089, Ω b = 0.0486, ℎ = 0.6774,  8 = 0.8159, and   = 0.9667 (Planck Collaboration et al. 2016).

Simulation Method
We use the highest resolution version of IPMSim introduced in Mandelker et al. (2019b) and Mandelker et al. (2021).We briefly describe key aspects of the simulation below and refer the reader to those papers for additional details.The simulation was performed with the quasi-Lagrangian moving-mesh code AREPO (Springel 2010).The goal of the simulation was to zoom-in not just on a single halo, as is commonly done, but rather on two protoclusters and the cosmic-web elements that connect them.To select our target region, we first considered the 200 most massive haloes in the  ∼ 2.3 snapshot of the Illustris TNG1001 magnetohydrodynamic cosmological simulation (Springel et al. 2018;Nelson et al. 2018;Pillepich et al. 2018b).These have virial masses2  v in the range ∼ (1 − 40) × 10 12 M ⊙ .We then selected from all these halo pairs with a comoving separation in the range (2.5 − 4.0)ℎ −1 Mpc ∼ (3.7 − 5.9) Mpc.There are 48 such halo pairs, each of which is either connected by a cosmic-web filament or else embedded in the same cosmic sheet.We randomly selected one pair consisting of two haloes with  v ∼ 5 × 10 12 M ⊙ each, separated by a proper distance of  ∼ 1.2 Mpc at  ∼ 2.3.By  = 0, the two haloes evolve into mid-size groups with  v ∼ (1.6 − 1.9) × 10 13 M ⊙ separated by ∼ 2.7 Mpc, such that their comoving distance has decreased by < ∼ 30%.The zoom-in region is the union of a cylinder with length  and radius  ref = 1.5 ×  v,max ∼ 240 kpc that extends between the two halo centers, and two spheres of radius  ref centred on either halo, where  v,max is the larger of the two virial radii at  = 2.3.We trace all dark matter particles within this volume back to the initial conditions of the simulation at  = 127, refine the corresponding Lagrangian region to higher resolution, and rerun the simulation to a final redshift,  fin .In this paper, we use the highest resolution version of this simulation, dubbed ZF4.0 in Mandelker et al. (2021), which has a dark matter particle mass of  dm = 8.2 × 10 4 M ⊙ and a Plummer-equivalent gravitational softening of  dm = 250 pc comoving.Gas cells are refined so that their mass is within a factor of 2 of the target mass,  gas = 1.5 × 10 4 M ⊙ .Gravitational softening for gas cells is twice the cell size, down to a minimal gravitational softening  gas = 0.5 dm = 125 pc.The simulation was run until  fin ∼ 2.9, though our analysis here focuses on  ∼ 4.
The simulations were performed with the same physics model used in the TNG simulations, described in detail in Weinberger et al. (2017) and Pillepich et al. (2018a).This includes radiative cooling down to  = 10 4 K from Hydrogen and Helium (Katz et al. 1996), metal-line cooling (Wiersma et al. 2009), radiative heating from a spatially constant and redshift-dependent ionising ultraviolet background (UVB) (Faucher-Giguère et al. 2009), partial self-shielding of dense gas from the UVB (Rahmati et al. 2013), and additional heating from the radiation field of active galactic nuclei (AGN) within 3 v of haloes containing actively accreting supermassive black holes (Vogelsberger et al. 2013).Star-formation occurs stochastically in gas with densities greater than  thresh = 0.13 cm −3 , which is placed on an artificial equation of state meant to mimic the unresolved multiphase ISM (Springel & Hernquist 2003).We include feedback from both supernova and AGN following the Illustris-TNG model.
To identify dark matter (sub)haloes in the simulation, we first apply a Friends-of-Friends (FoF) algorithm with a linking length  = 0.2 to dark matter particles.We then assign gas and stars to FoF groups based on their nearest-neighbour dark matter particle.Finally, we apply SUBFIND (Springel et al. 2001;Dolag et al. 2009) to the total mass distribution in each FoF group.The most massive SUBFIND object in each FoF group is identified as the central halo, whose virial radius  v is defined using the Bryan & Norman (1998) spherical overdensity criterion, and the virial mass  v is the total mass of dark matter, gas, and stars within  v .

Filament Selection
We focus our analysis on  ∼ 4. At this time, our system consists of a well-defined cosmic sheet containing several prominent co-planar filaments, with end-points at massive haloes with  v ∼ 10 12 M ⊙ , and along which nearly all haloes with  v > 10 9 M ⊙ are located.The sheet is formed by the collision of two smaller sheets at  > ∼ 5, while most of the filaments merge by  ∼ 2.5 (Mandelker et al. 2019b(Mandelker et al. , 2021)).Figure 1 shows the Hydrogen column density,  H , of the sheet at  ∼ 4 in face-on (left) and edge-on (right) projections 3 .In order to minimise contamination of the filament properties by haloes, we avoid analysing filament regions that contain haloes with  v ≥ 10 10 M ⊙ , marked by black circles in Fig. 1.This is ∼ 1% of the mass of each of the three haloes at the nodes of the cosmic-web in this region, located at the bottom-left, bottom-right, and top-right of the left-hand panel in Fig. 1.We selected ten such slices from three separate filaments, marked with white rectangles and numbered from 1 − 10 in Fig. 1.Each slice is 30 kpc thick along the filament axis, slightly larger than the typical characteristic scale radius of the filaments,  0 , as described in Section 3.1 below.However, we note that using a slice thickness of (20 − 40) kpc does not change our key results.We note that most of our slices contain 1-2 haloes with masses 10 9 <  v /M ⊙ < 10 10 , some of which appear assoiated with local distortions in the density/temperature structure of the filament.However, these distortions average out when stacking the slices, and have no impact on our results.Likewise, all of our slices contain several haloes with  v < 10 8 M ⊙ , though none of these appear to have any impact on the filament structure.
The orientation of the filament axis in each slice is determined by eye by approximating a straight line to the ridge line of maximal  H in each slice, as seen in the face-on projection through the sheet (left-hand panel of Fig. 1).Given the relatively small length of our slices compared to the total length of the filament, small changes to the axial directions in each slice do not change any of our results.We hereafter define the -axis to refer to the local filament axis in each slice, while the (, ) plane represents the filament cross-section with the -axis always aligned within the sheet.
To find the filament centre in each slice, we project the gas density in each slice along its axial direction and initially place the filament centre at the point of maximum gas density.We then refine this using the "shrinking-cylinders" method.We begin by calculating the centre of mass of a cylinder of gas with radius  0 = 50 kpc and length  = 30 kpc about our initial centre, i.e., spanning the entire thickness of the slice, and update the centre to this position.We repeat this process five times, where the radius of the cylinder at each step is half that of the previous step, leading to a final radius of 50/2 5 ∼ 1.6 kpc, while keeping  unchanged.Changing the number of iterations from 4 − 6 has no effect on our results.We experimented with alternative definitions of the filament centre, such as using shrinking-circles of gas density projected along the length of the slice rather than shrinking cylinders in 3D, using spheres rather than cylinders/circles, and using total (gas plus dark matter) mass rather than gas mass.All these definitions yield centres within ∼ (1 − 4) kpc of each other for nearly all cases, and the differences have no bearing on our results.

Three Radial Zones
As will be made clear and expanded upon throughout Sections 3-6 below, we find that the filament cross section can be broadly divided into three radial zones, with a fourth zone outside the stream boundary.Since the details of each zone are motivated at different points in the paper, yet we refer to each throughout the paper and mark them on most figures, we introduce them here for the benefit of the reader.A schematic diagram showing this proposed structure is presented in Fig. 2.
∼  shock : The postshock region, where the gas is hot (at the virial temperature, see Section 5) and the dynamics are dominated by thermal and ram pressure.
(ii) Intermediate vortex (V) zone, 0.25 shock < ∼  < ∼ 0.65 shock : The gas is cooling and the dynamics are dominated by a quadrupolar vortex structure and increasing centrifugal forces (see Section 6).
(iii) Inner stream (S) zone, 0 < ∼  < ∼ 0.25 shock : A dense, isothermal core in the inner filament, representing the cold stream that penetrates the CGM of massive haloes.
In addition, we define a fourth zone outside the filament boundary, the Pre-shock (PS) zone at  > ∼  shock .This is the region outside the accretion shock surrounding the filament, where cold and diffuse gas is free-streaming towards the filament.While  shock is well-defined (Section 3.1), the boundaries between the T and V zones and between the V and S zones are not as sharp.We crudely assign a width of ∼ 0.05 shock to each of them, / shock ∼ (0.6−0.65) and (0.2−0.25) respectively, marked in the figures below.This roughly corresponds to the width of transition regions of various zone characteristics in stacked data, as detailed throughout the paper.However, we caution that the precise width of each zone likely depends on filament properties and is not universal.A schematic diagram showing the structure of the filament cross section, marking the key physical characteristics of the three zones described in Section 2.3 and detailed throughout the paper.The stream is surrounded by a cylindrical accretion shock, at a normalised radius of R = 1, which expands outwards as the filament grows.The shock-heated gas in the outer Thermal (T) zone, at R > ∼ 0.65, is at the virial temperature of the underlying dark matter filament, with low density and inwardsradial velocity driven by an isobaric cooling flow.In this region, the thermal pressure gradients balance both gravity and the ram pressure from the inflowing gas.In the intermediate Vortex (V) zone, at 0.65 > ∼ R > ∼ 0.25, the dynamics is dominated by a quadrupolar vortex structure which is induced by shearing motions between the sheet that feeds the filament and the post-shock filament gas.These vortices induce a complex combination of inflows, outflows, and rotation in this zone.Gravity here is balanced by a combination of thermal pressure gradients, centrifugal forces, and non-thermal random motions.The inner stream (S) zone, at R < ∼ 0.25, is a dense, cold, isothermal core which represents the cold streams that feed massive galaxies.While the stream boundary expands as the filament grows, the gas in this region is also characterised by both ithats and outflows due to non-radial inflow along the sheet.

FILAMENT RADIAL STRUCTURE
In this and the following sections, we analyse the radial structure of various filament properties.These include the gas density, temperature, and thermal pressure (Section 3.1), the baryon fraction (Section 3.2), the dark matter density and radial velocity dispersion (Section 3.3), the mass-per-unit-length of gas and the total mass (Section 3.4), the ratio of cooling time to free-fall time (Section 4), the virial parameter (Section 5), and the internal kinematics and dynamics (Section 6).We compute all profiles as a function of the cylindrical radius, i.e., the distance from the local filament axis as defined in Section 2.2, with each radial bin representing a cylindrical shell extending the thickness of the slice,  = 30 kpc.

Thermal Properties and Characteristic Radii
In Fig. 3, we address the distribution of gas density (top), temperature (centre), and thermal pressure (bottom), stacked among all ten slices.We present stacked projection maps on the left, and stacked radial profiles on the right.The projection maps represent weighted-averages along the filament axis (the  direction) across the slice thickness of  = 30 kpc.The average density and pressure are weighted by volume, while the temperature is weighted by mass.Prior to stacking, we align the -axes such that they lie in the sheet containing the filaments (see Fig. 1), while the -axes remain perpendicular to the sheet.We then normalize each property by its value at the filament centre, and take the average in log-space among all ten slices.The cross-sections are roughly circular with high-density, low-temperature cores surrounded by a low-density, high-temperature medium, and roughly constant pressure throughout.The vertical high-density region extending from the filament in the -direction breaking the circular symmetry represents the sheet.There are well-defined accretion shocks around both the filament and the sheet, visible in both the temperature and pressure maps.The filamentary accretion shocks identified here have typical Mach numbers of order M > ∼ 10 (Appendix A).While the accretion-shock around the sheet is discussed in Mandelker et al. (2019bMandelker et al. ( , 2021)); Pasha et al. (2023), those works did not explicitly identify filamentary accretion shocks as seen here.Such filamentary accretion shocks have been both predicted theoretically (e.g., Klar & Mücket 2012;Birnboim et al. 2016) and seen in other cosmological simulations (Ramsøy et al. 2021).While the sheet clearly breaks the circular symmetry of the filament cross-section, for most of our analysis (except Fig. 9 below) we ignore this complication and treat the filament as circular.A more accurate analysis distinguishing between the on-sheet and off-sheet properties and dynamics is deferred to future work.
When generating the profiles in the right-hand column, we first use 150 linearly spaced bins out to a radius of  = 100 kpc for each individual slice, namely a bin size of Δ = 100 kpc/150 ≈ 0.67 kpc.While this is larger than the typical cell size near the filament centre, it is slightly smaller than the typical cell size at large distances from the filament axis (see Fig. C1 below).However, due to the unstructured nature of the grid, all of our bins are sufficiently populated for good statistics, with ∼ (4000 − 6000) cells per bin in the S and V zones, and ∼ (2000 − 4000) cells per bin in the T Zone.The density and pressure (temperature) profiles represent the volumeweighted (mass-weighted) average of all gas cells in each cylindrical shell.In each slice, the temperature profile reaches a well-defined maximum, which we define as the shock-radius,  shock .Once  shock is identified, we recompute each profile using 150 linearly spaced bins from / shock = 0.1 − 2.0 and normalise them by their central values.We then compute the mean (red dots) and standard deviation (cyan shaded regions) in log-space among all ten slices in each bin to produce the stacked profiles shown in Fig. 3. Hereafter, all stacked profiles and projection maps presented used binonproportional to  shock of each slice.The location of the shock seen on the temperature and pressure maps is indeed at  ≃  shock , marked by a dashed-white circle in each projection.
At small radii, the filaments are characterised by a core of constant density and temperature extending to roughly ∼ (0.2 − 0.25)  shock , namely the S zone defined in Section 2.3 4 .Outside of this radius, the temperature begins increasing towards a maximum at the shock radius, marked in each profile panel with a vertical dashed blue line.On average, the peak temperature at the shock is ∼ (10 − 20) times larger than the core temperature.This is consistent with the shock Mach number, which is measured to be M ∼ (10 − 20) Figure 3. Thermal structure of filaments, stacked among the ten slices.We show projected maps integrated along the filament axis (left) and radial profiles (right) of gas density (top), temperature (middle), and thermal pressure (bottom) of the text for details regarding the averaging and stacking procedures.The temperature profiles reach a well-defined maximum that defines a shock radius,  shock , marked by white circle in the projection maps and vertical blue lines in the radial profiles.The temperature declines towards smaller radii, reaching an isothermal core in the S Zone,  < ∼ 0.25  shock , where the gas is dense and cold.In the V zone, / shock ∼ (0.25 − 0.65), the pressure is roughly constant while the temperature increases by a factor of ∼ 3.In the T zone, / shock ∼ (0.65 − 1.0), the pressure drops by a factor of ∼ 2, the temperature increases by a factor of ∼ 5, and the density develops a slope of  ∝  −3 .Beyond the shock radius, in the PS zone, the density becomes dominated by the sheet and develops a shallower profile.The boundaries between the S and V zones and between the V and T zones are marked by concentric circles in the projection maps and by shaded regions of width 0.05  shock in the radial profiles.At  <  shock , the density profile is well fit by an infinite self-gravitating isothermal cylinder in hydrostatic equilibrium, shown by the black dashed line, with a scale radius of  0 ∼ 0.5  shock marked by a yellow star in the density profile.However, this model is a very poor fit to the pressure profile, where strong gradients are only present in the T zone, while the S and V zones are nearly isobaric.for all slices (see Appendix A).Within the V zone,  ∼ (0.25 − 0.65) shock , the pressure remains roughly constant, decreasing by only ∼ 10% from its central value.The temperature increases by a factor of ∼ 3 from its central value, often considered the threshold for 'cool' gas in studies of multiphase gas mixing (e.g., Scannapieco & Brüggen 2015;Gronke & Oh 2018;Mandelker et al. 2020a), though the density declines by a larger factor in this region.The cold and dense filament core in the S zone thus seems to be in pressure equilibrium with the hot and diffuse mixed gas in the V zone.In the T zone,  ∼ (0.65 − 1.0) shock , the pressure decreases by another factor of < ∼ 2, the temperature rapidly increases by a factor of ∼ 5, and the density reaches a slope of roughly  ∝  −3 at  < ∼  shock .As will be discussed further below, thermal pressure gradients are indeed not the dominant force in the S and V zones, becoming important only in the T zone.Outside the shock, in the PS zone, the pressure and temperature rapidly decline, while the slope of the density profile becomes noticeably shallower.In this region, the density is dominated by the underlying sheet, rather than by the filament itself (see also Ramsøy et al. 2021).The boundaries between the four radial zones are marked by vertical lines/shaded regions in the radial profiles and by concentric circles in the projections.
We show the values of central density,  0 , central gas temperature,  0 , and shock radius,  shock , for our ten slices in Fig. 4. Most of the slices have central densities  0 ∼ (10 −26.5 − 10 −25.5 ) cm −3 , corresponding to a Hydrogen number density of  H ∼ 0.01 cm −3 and to a baryonic overdensity of < ∼ 100.Their central temperatures are  0 ∼ 2 × 10 4 K, placing the filament core in approximate thermal equilibrium with the UVB (e.g., Mandelker et al. 2020a).The shock radii are typically  shock ∼ (35 − 50) kpc, with an average value of ∼ 45 kpc among the ten slices.The typical core size is thus  core ∼ 0.25  shock ∼ 12 kpc.These values of  0 ,  0 , and  core are in good agreement with those predicted for cold streams feeding massive high- galaxies from the cosmic-web (Dekel et al. 2009a;Mandelker et al. 2018Mandelker et al. , 2020b)).It thus seems that the cold streams predicted to penetrate the CGM of massive haloes and feed their central galaxies can be associated with the isothermal cores of cosmic-web filaments, motivating our definition of the S zone.
Slices 05 and 10 seem to have much hotter cores, with  0 > ∼ 10 5 K.Both of these slices lie close to the intersection between two filaments (see Fig. 1 left), leading to a highly disturbed morphology where the cold region of the filament is off-centre.However, even in these cases, the minimum temperature in the filament is < ∼ 3 × 10 4 K. Similarly, slices 01 and 02 have unusually extended shocks, with  shock ∼ (70 − 80) kpc.Slice 01 is adjacent to a massive, ∼ 10 12 M ⊙ halo, which likely influences both the extent and the Mach number of the shock, the latter of which is > ∼ 50% larger than the rest of the sample (see Appendix A).The large shock radius in slice 02 seems to be due to curvature in the filament in this region, making it difficult to track the filament spine.However, these slice-to-slice variations do not influence our stacked results.
To further characterise the filament radial structure, we fit the stacked density profile shown in the upper-right hand panel of Fig. 3 to the density profile of an infinitely long, self-gravitating, isothermal cylinder in hydrostatic equilibrium (Ostriker 1964), with where  0 is the central density,  is the gravitational constant,  B is the Boltzmann constant,  p is the proton mass,  is the mean molecular weight, and  0 is the isothermal temperature.While the filament is clearly not isothermal within  shock , several previous works have attempted to model intergalactic filaments as such (e.g., Harford & Hamilton 2011;Mandelker et al. 2018;Ramsøy et al. 2021).It is, therefore, interesting to see how well this model can describe the actual density profile.This is shown by the dashed black line in the upper-right-hand panel of Fig. 3.As we can see, it actually fits the density profile quite well, except at large radii,  > ∼  shock , where the density becomes dominated by the sheet rather than the filament.Similar results were found in Ramsøy et al. (2021), who attempted to fit the density profile with a combination of an isothermal cylinder and an isothermal sheet, though we make no such attempt here.
When fitting our stacked density profile to eq. ( 1), the only free parameter is  0 / shock , since the density has been normalized to its central value so effectively  0 = 1 5 .From the fit, we derive a characteristic scale radius of  0 ∼ 0.5  shock , marked by a yellow star on the density profile.This is roughly twice as large as the radius of the cold, dense, isothermal core that represents the cold stream.It is slightly smaller than the outer radius of the V zone, though there is no clear relation between these two radii.
We can use the isothermal fit to the density profile to find the corresponding temperature and pressure profiles.The isothermal temperature is given by and the isothermal pressure profile is given by Plugging the average value of  0 ∼ 10 −26 cm −3 and the derived value of  0 ∼ 0.5 shock ∼ 22 kpc into eq.( 3) returns a value of  0 ∼ 3.6 × 10 4 K, which is roughly the average value of  isothermal measured among our ten slices (see Fig. 4).The density profile thus resembles that of an isothermal filament at the core temperature.We discuss this further in Section 7.1.The model pressure profile, unsurprisingly, drastically underpredicts the pressure outside the core region (S zone), where the temperature increases and deviates strongly from isothermality.This suggests that thermal pressure gradients are not the primary force supporting the filaments against gravity, at least not in the S and V zones where the pressure is nearly constant both on and off the sheet.From the projection map we see that in the T zone, the off-sheet pressure declines more rapidly, leading to stronger gradients in this region.This will be discussed further in Section 6.2.

Baryon Fraction
In Fig. 5 we show the baryon fraction as a function of radius, stacked among our ten slices as described in Section 3.1.The baryon fraction in each radial bin,  b (), is given by the ratio of the baryonic mass (gas and stars) to the total mass (gas, stars, and dark matter) enclosed within the radius .In practice, the stellar mass is negligible, contributing ∼ 0.3% of the total baryon mass, and  b can be thought of as the gas fraction.In the PS zone,  > ∼  shock ,  b () converges to the universal baryon fraction,  b ∼ 0.16, with very little scatter among the different slices.Within the filament core, the S zone,  b , saturates at roughly 0.5 on average, with a large scatter of ∼ (0.3 − 0.7).The gas is thus more centrally concentrated than dark matter due to efficient cooling, resulting in a comparable contribution of gas and dark matter in the core.This is similar to galaxies in dark matter haloes, for example, in the Milky Way the baryon-to-total mass ratio within the solar circle is ∼ 0.5.

Dark Matter Properties
Following the discussion in Section 3.1 and Section 3.2, we now examine the radial structure of the dark matter in our filament slices.
In the left panel of Fig. 6, we show a projected map of the dark matter density, oriented and stacked among all ten slices in an identical way to the projection map of gas density shown in Fig. 3. Similar to gas, the dark matter filament has a roughly circular cross-section if ignoring the contribution of the overdense pancake along the  direction.However, unlike the gas density which is very centrally concentrated, the dark matter density appears approximately constant within  shock , marked with a white circle.The dark matter distribution is also much clumpier than the gas, due to the presence of low mass haloes with  v < 10 10 M ⊙ (recall that higher mass haloes have already been removed when selecting the slices).Removing haloes with lower and lower masses prior to projecting the density results in smoother and smoother maps, until the clumpiness disappears once haloes with  v > ∼ 10 8 M ⊙ are removed.This clumpiness thus does not seem to be due to internal fragmentation within the filament, but rather to the cosmological halo mass function.
In the middle panel of Fig. 6, we present the radial profile of dark matter density, normalised by its central value and stacked among our ten filament slices, using the same procedure described in Section 3.1 for the gas profiles.We fit the dark matter density profile to the same parametric isothermal model used for the gas (eq.1).The fit is shown in the middle panel of Fig. 6 with a black-dashed line and is a good representation of the data.From Fig. 5 we know that the central dark matter density is comparable to the central gas density, namely an overdensity of ∼ 30 with respect to the mean matter density6 at  ∼ 4. From the fit, we extract the scale radius of the DM density  0,DM ≃ 1.17  shock , which is marked in Fig. 6 along with the scale radius of the gas density profile,  0 < ∼ 0.5  shock .The radius of the DM density core is  core, DM ∼ 0.8 shock , approximately 4 times larger than the gas density core.This is consistent with the baryon fraction profile (Fig. 5), which showed that the gas was much more centrally concentrated than dark matter.This is unlike the case studied in Ramsøy et al. (2021), where gas and dark matter were found to have similar concentrations and density profiles.We address this in Section 7.1.
In the right-hand panel of Fig. 6, we present the radial profile of dark matter effective temperature, normalised by its central value and stacked among our ten filament slices as described above.This is derived from the radial velocity dispersion of dark matter, which we compute in the standard way as , where   is the radial velocity with respect to the filament axis and ⟨. . .⟩ represents a mass-weighted average7 over all dark matter particles in a given  3.The shock radius,  shock , identified based on the gas temperature profiles, is marked by a white circle.As in Fig. 3, we also mark the boundaries between the four radial zones with concentric circles.Middle and right panels: Structure of the dark matter filaments, stacked among our ten slices.We show the radial profiles of dark matter density (middle) and effective temperature (right, computed using the dark matter radial velocity dispersion following eq.5).These were normalised, stacked, and averaged in the same way as the gas profiles shown in Fig. 3, with red points (cyan shaded region) representing the mean (standard deviation) in log space.The -axes in these two panels span the same range as the upper and middle right-hand panels of Fig. 3.The black dashed line in the density profile shows a fit to the same parametric isothermal model used for the gas (eq.1), which matches the data well.The yellow and green stars mark the scale radii for the gas and dark matter density profiles, respectively.The latter is slightly larger than the shock radius,  shock , marked with a vertical blue line.The vertical shaded regions mark the boundaries between the S and V zones and the V and T zones.The effective dark matter temperature increases by a factor of ∼ 2 from the filament core to a peak near  0,DM , and is thus more isothermal than the gas, whose temperature increases by a factor of ∼ 20 over the same radial range.cylindrical shell at radius .The effective dark matter temperature8 is then defined via From the profile, we see that  eff increases by a factor of ∼ 2 from the filament core to a peak just outside  shock , near  0,DM .This is much less than the factor ∼ 20 increase in gas temperature, showing that the isothermal approximation is more valid for dark matter than for gas.The average central value of  eff in the filament core is  eff,C ∼ 9.6 × 10 5 K. Using the average central dark matter density, our fit to  0,DM , and eq.( 3), we calculate the effective "isothermal" temperature of the dark matter,  iso,DM = 1.2 × 10 6 K.This is very similar to the measured central dark matter temperature, supporting the validity of the isothermal model for the dark matter.We note that dark matter haloes are also often modelled as isothermal spheres, though as singular isothermal spheres with a central density cusp,  ∝  −2 , rather than a central core as we find for filaments.

Mass per unit length (line-mass)
In Fig. 7 we show radial profiles of the mass per-unit-length (hereafter line-mass), Λ(< ), normalised and stacked among our ten filament slices.For a given slice, Λ(< ) is simply defined as the total mass (of gas and/or dark matter, see below) interior to radius , divided by the thickness of the slice,  = 30 kpc.Motivated by our fits of the gas and dark matter density to isothermal profiles (Section 3.1, Section 3.3), prior to stacking we normalise the line-mass in each slice by the maximal possible value of a self-gravitating isothermal cylinder in hydrostatic equilibrium (Ostriker 1964), where  2  =  B /( p ) is the isothermal sound speed squared, and  is the effective isothermal temperature given by the fit of the model using the average central density, the characteristic scale radius,  0 , and eq.( 3).This is similar to the central temperature in the filament core.Filaments with Λ ∼ Λ max are unstable to gravitational fragmentation and radial collapse (Mandelker et al. 2018;Aung et al. 2019), while much lower mass filaments are likely gravitationally stable.
In Fig. 7, we show separately the normalised line-mass profiles of gas alone (red) and total mass (blue, primarily gas and dark matter).Note that the total line-mass is obviously larger than the line-mass of gas alone, However, the gas has a much lower effective temperature than the dark matter, ∼ 3.6 × 10 4 K and ∼ 1.2 × 10 6 K respectively, while the effective temperature derived from the density profile of total mass is ∼ 8.6 × 10 5 K.This leads to a much smaller value of Λ max for the gas alone compared to the total mass -2.3 × 10 8 M ⊙ kpc −1 and 5.5 × 10 9 M ⊙ kpc −1 respectively.As a result, the normalised line-mass of gas alone is larger than the normalised total line-mass.
For the gas, Λ/Λ max ∼ 0.2 near the outer boundary of the S zone at  ∼  core ∼ 0.25  shock , Λ/Λ max ∼ 0.7 near  shock , and it approaches unity at 2  shock .This is consistent with the results of Mandelker et al. (2018) (see their equation 40), who predicted that cold streams feeding massive high- galaxies should have Λ ∼ Λ max , which could lead to gravitational fragmentation and star-formation in filaments outside of galaxies.However, we stress that this estimate ignores non-thermal support such as turbulence, which is important Profiles of the filament line-mass, Λ, accounting for only gas (red points) and total mass (mostly gas and dark matter, blue points), averaged among our ten filament slices.Each profile has been normalised by the corresponding, temperature-dependent, maximal line-mass for hydrostatic equilibrium, Λ max (eq.6).Cyan and orange shaded regions represent 1 −  scatter of Λ gas and Λ tot among our ten filament slices, respectively.Note that since all slices have the same length (Section 2.2), this result is also the same as the total line-mass of a single filament consisting of all our ten slices connected end-to-end.While the total line-mass is clearly larger than the line-mass of gas only, the higher temperature of the dark matter results in a higher value of Λ max , leading to a lower normalised line-mass.While the normalised linemass for gas approaches unity near the filament edge, the total line-mass is well below the critical value.
at small radii (see Section 6.2 below), and the non-isothermality of the gas which is important at large radii (see Fig. 3).The net effect on the gravitational fragmentation of star formation in filaments is thus unclear and is left for future work.When considering the total filament line-mass we have Λ/Λ max ∼ 0.02 and 0.15 within  core and  shock respectively, suggesting that gravitational fragmentation of the dark matter in the filament is unlikely, owing to its much larger effective temperature.

THERMAL EQUILIBRIUM
In analogy with the gaseous haloes around massive galaxies, groups, and clusters, we estimate the thermal stability of the shock-heated gas surrounding the isothermal core by examining the ratio of the gas cooling time,  cool , to the free-fall time,  ff .
The gas cooling time,  cool , is given by where  is the gas internal energy per unit mass,  is the gas density,  H is the hydrogen number density, and L = C −H is the net cooling rate, i.e., cooling minus heating, per unit density squared.All these properties are stored for each cell in the simulation, allowing us to evaluate  cool for each cell.We then compute the average cooling time as a function of radius,  cool (), by taking the mass-weighted average of eq. ( 7) among all cells in each radial bin with temperatures  >  0 = 3.6 × 10 4 K, the temperature of the cold isothermal core.
The gas with  < ∼  0 is near thermal equilibrium with UVB and often undergoes net heating, while virtually all gas with  >  0 is undergoing net cooling.
For an infinite cylinder, the free-fall time is given by where  mean (< ) = Λ(< )/( 2 ) is the mean density interior to .
For isobaric cooling,  ∝  ∼ const, so  cool ∝  −2 , while  ff ∝  −1/2 .Therefore, the ratio  cool / ff is smaller at smaller radii closer to the filament axis.When  cool / ff < ∼ 1, the system cannot maintain the shock-heated component, and the bulk of the gas cools and falls to the centre (e.g., Rees & Ostriker 1977;White & Rees 1978;Birnboim & Dekel 2003;Birnboim et al. 2016;Fielding et al. 2017;Stern et al. 2020aStern et al. , 2021a)).The outer radius where  cool / ff = 1 and the volume-filling medium transitions from hot to cool corresponds to the 'sonic radius' (Bertschinger 1989;Stern et al. 2020bStern et al. , 2023)).Even if  cool / ff > 1, the hot gas may be unstable to local thermal instabilities where perturbations create dense clouds that cool and condense out of the hot medium and "rain down" to the centre.In studies of the hot CGM around massive galaxies or the ICM in galaxy clusters, this process is often referred to as precipitation and occurs when  cool / ff < ∼ 10 (e.g., McCourt et al. 2012;Sharma et al. 2012;Gaspari et al. 2012;Voit & Donahue 2015;Voit et al. 2015a,b).The analogous threshold for cylindrical geometry may be different due to the different radial dependences of stabilising buoyancy forces.A detailed study of this is beyond the scope of this paper, and instead we adopt the threshold of  cool / ff ∼ 10 as representative of the threshold for condensation to occur.
Figure 8 shows the ratio  cool ()/ ff () as a function of / shock , stacked among the ten slices by taking the average and standard deviation.At  ∼  shock ,  cool / ff ∼ 3, which implies that while the shock may be stable to monolithic cooling it is likely unstable to precipitation and condensation.Indeed, we see several cases of cold clouds encompassed by hot gas in the T and V zones in individual slices.A more detailed analysis of these clouds is beyond the scope of the current study, and is left for future work.The cooling radius, where  cool / ff ∼ 1, occurs near the outer edge of the S zone, consistent with this being the boundary of the dense isothermal core.These results suggest that the cold stream is built by cooling of the post-shock filament gas, and that this cooling occurs roughly isobarically (see Fig. 3).

VIRIAL EQUILIBRIUM
In this section, we wish to ascertain to what extent the gas is in virial equilibrium within the gravitational potential of the filament.While dark matter haloes are expected to be bound and virialized in three dimensions, intergalactic filaments will at most be bound and virialized in the two dimensions perpendicular to the filament axes, while remaining unbound along their axis.Previous studies have speculated that intergalactic filaments may result from cylindrical collapse ending in virial equilibrium per-unit-length (Fillmore & Goldreich 1984;Mandelker et al. 2018).In this section, we test this hypothesis in our simulations.∼ 3 ff implying that the shocked medium is likely unstable to condensation and precipitation.In the S Zone,  cool < ∼  ff consistent with the formation of a cold isothermal core at these radii.

The Virial Theorem per-unit-length
An infinite, self-gravitating, collisionless and isolated cylinder in virial equilibrium (per-unit-length) obeys (Mandelker et al. 2018 where K is the kinetic energy per-unit-length and Λ is the linemass.If the motions are confined to the plane perpendicular to the filament axis, this becomes the two-dimensional version of the virial theorem.In a confined, gaseous system such as ours, we must account for additional terms related to thermal pressure, surface pressure, and magnetic pressure. As we show in Appendix B (Fig. B1), magnetic support is negligible in the filament gas.Therefore, we use the virial theorem for an unmagnetized gaseous system.For such a system in time-independent equilibrium, the virial theorem states (Krumholz 2015, Chapter 6.1) where T is the total thermal plus kinetic energy in the system, T  is the confining pressure on the surface, including both thermal and ram pressure, and W is the total gravitational energy of the system.In integral form, these terms can be expressed as In eqs.( 11)-( 13),  is the gas density,  is its velocity,  is its thermal pressure, Φ is the gravitational potential9 , and ← →  is the fluid pressure tensor given by where    is the Kronecker delta.We define the virial parameter Note that this differs from the standard virial parameter for unmagnetized systems only in the inclusion of the surface pressure term, which is important for non-isolated systems10 .When  vir > 1, the system expands due to kinetic plus thermal energy, while for  vir < 1 the system collapses due to the combined effects of surface pressure and gravity.
When evaluating eqs.( 11)-( 13) for our filament slices, we recall that we are only interested in virial equilibrium per-unit-length, in the two dimensions perpendicular to the filament axis.We thus make an approximation of cylindrical symmetry and treat the filament as an infinite cylinder.We are here neglecting the non-axisymmetric nature of the sheet, an approximation which will be validated by the results below.Therefore, we approximate the gravitational acceleration as with Λ tot (< ) being the total line-mass interior to radius .Similarly, the gas density and line-mass are related by We, therefore, approximate the gravitational potential energy as where  = 30 kpc is the filament slice thickness.We note that this factor  is only included for consistency with the formalism presented in Krumholz (2015), and in practice, we normalise all the energy terms in eqs.( 11)-( 13) by  to obtain energy per-unit-length before computing  vir .We numerically evaluate the integral in eq. ( 18) by performing the following double-sum: where the sum < ′ in the square brackets is taken over all particles (dark matter and stellar) and gas cells interior to radius  ′ .The sum  ′ is taken over all gas cells within the cylindrical shell of thickness Δ ′ at  ′ .
We numerically evaluate the volumetric kinetic and thermal terms (eq.11) as where  = √︃  2  +  2  +  2  is the total velocity.Note that we include  2  in T despite our assumption of filaments as infinite cylinders when computing W. This is because even when evaluating the energy per-unit-length, such that there is no net velocity or acceleration along the  direction, we still expect the velocity dispersion to be three-dimensional locally, such that  2  remains an important contribution to the kinetic energy density of the system.We comment in Section 5.3 on the effect of removing  2  from the kinetic term.The surface term, eq. ( 12), comprises three surfaces for each radius , one in the cylindrical shell, and two others at either end of the cylindrical slice.However, the latter two surfaces are irrelevant for evaluating the virial equilibrium per-unit-length.Therefore, we do not consider them here.We comment in Section 5.3 on the effect of including them in the calculation.The surface term is thus given by The second term, namely the integral of     , is not compatible with our assumption of an infinite cylinder, nor with the evaluation of T  () per-unit-length.For an infinite cylinder in equilibrium perunit-length,  z would be constant (or zero).Likewise, eqs.( 16)-( 18) are only valid for  which does not depend on , and in such cases  r should also be independent of .This term thus reduces to the integral of d over a symmetric interval and is, therefore, 0. In light of this, we neglect this term moving forward.In practice, it turns out that this term is negligible anyway, at most a few percent of the integral of  2  at all radii  > ∼ 0.15  shock , because     is indeed roughly independent of , consistent with our assumed symmetry.
In practice, when numerically evaluating T  we perform the following sum over gas cells, where  is the volume of each gas cell and Δ is the width of the bin.Therefore, we have approximated d ≃ /Δ.We have verified that this approximation is roughly independent of our choice of Δ, as long as this is within a factor of ∼ (2 − 3) of the typical cell size in the simulation.

Rapidly Flowing Gas on and off the Sheet
Finally, before we calculate the corresponding profiles, we wish to remove gas that is rapidly accreting onto or outflowing from the filaments.This gas is not expected to be in equilibrium within the filament, but rather to be unbound or on its first infall, and therefore will bias our results (see Lochhaas et al. 2021, for a similar discussion in the context of DM haloes).We do this by applying a cut in the radial velocity of gas cells included in the energy terms described above.At each radius , we calculate the corresponding "virial" velocity (from eq. 9) and only consider gas cells where |  | <  vir ().We apply this threshold to all terms appearing above, except for Λ tot which sets the gravitational potential.In practice, this removes almost exclusively gas which is rapidly inflowing along the sheet onto the filament.
To illustrate this, we present in the left-hand panel of Fig. 9 radial profiles of the mass-weighted average radial velocity, stacked among our ten filament slices as in previous figures.In the right-hand panel, we show the radial mass flux, , stacked among the ten slices.In both panels, we show separately the profiles for material within and outside the sheet.We crudely define the sheet region for all the slices as anything within || < 10 kpc ∼ 0.22  shock in the frame of the projection maps in Figs. 3 and 6.In the left-hand panel, we also show the profile of ± vir .It is evident that the radial velocities towards the filament are quite different within and outside the sheet, and that our velocity threshold of |  | <  vir () for the virial parameter predominantly removes rapidly inflowing gas (with   < 0) along the sheet in the T and V zones.The mean inflow velocity along the sheet exceeds  vir at  > ∼ 0.3  shock , while outside the sheet it only exceeds  vir at  > 1.2 shock .
Besides their use for our calculation of the virial parameter, it is worth discussing the structure of the inflow velocity and mass accretion profiles in detail.In the PS zone at  < ∼ 2 shock , and throughout the T zone, the off-sheet gas is decelerating as a result of the strong thermal pressure gradients generated by the shock (see Section 6.2).At  > 2 shock , the off-sheet gas is either accelerating or maintaining constant velocity.However, note that the gas does not completely stall at  shock , as might have been expected for a strong accretion shock, because the cooling time is only slightly longer than the free-fall time at  < ∼  shock (Fig. 8) 11 .This generates a cooling flow and allows the gas to continue flowing towards the centre, albeit at a reduced velocity.This is supported by the fact the mass-flux of off-sheet material remains constant throughout the T zone, suggesting that the material does not "pile-up" behind the shock.On the other hand, the on-sheet gas maintains a roughly constant inflow velocity throughout the T zone, as it interacts only weakly with the shock.This is similar to how the inflow velocity of cold streams feeding massive galaxies at high- through their shock-heated haloes are seen in simulations to maintain roughly constant inflow velocities throughout the hot CGM (Dekel et al. 2009a;Goerdt & Ceverino 2015).
Within the V zone, at 0.2 < ∼ / shock < ∼ 0.6, the off-sheet gas is outflowing with an average outflow velocity smaller than  vir .This seems to be caused by a strong quadrupolar vortex structure that develops in this region, as discussed in Section 6.3.These vortices cause the gas outside the sheet to swirl around and can lead to outwards radial motions.In this same region, the gas within the sheet decelerates due to a combination of shear against the shock-heated gas and interaction with the vortices, both of which drain momentum from the inflowing gas (Mandelker et al. 2016;Padnos et al. 2018;Mandelker et al. 2019a), and also due to a non-zero impact parameter of the sheet with respect to the filament, which decreases the radial component of the velocity closer to the filament centre.In general, mass accretion into the S zone is dominated by flow outside the sheet in the PS zone, and by flow along the sheet in the T and V zones.
Within the S zone, the off-sheet component is not well-defined due to our crude definition of the sheet as everything within || < 10 kpc ∼ 0.22  shock .The sheet material continues inflowing towards the filament axis, slowly decelerating from ∼ 0.1  vir ( shock ) in the outer S zone towards smaller radii.

Virial Equilibrium and the Virial Radius
Using eqs.( 15), ( 18)-(20), and ( 22), we evaluate  vir for each slice as a function of radius for gas with |  | <  vir ().We discuss the sensitivity of our results to this threshold below.We then stack all the ten slices by weighting  vir () by |W ()|.This is equivalent to treating all ten slices as one long cylinder when computing the virial Negative values refer to inflow towards the filament axis while positive values refer to outflowing gas.We distinguish between gas outside the sheet (blue lines) and within the sheet (red lines), where the sheet region is crudely defined as |  | < 10 kpc ∼ 0.22  shock (see Fig. 3).The dashed-orange line in each panel shows the results for all gases regardless of their position relative to the sheet.The value of   in each bin is the mass-weighted average radial velocity, and the value of  in each bin is   /Δ, where Δ is the width of the bin.Prior to stacking, the   profiles for each slice are normalised by the corresponding virial velocity at  shock ,  vir ( shock ) = (Λ shock ) 1/2 (eq.( 23)), while the profiles of  were normalised by Λ (  shock )  vir ( shock ).In the left-hand panel, the green lines show the profiles of ± vir ( )/ vir ( shock ).The sheet is inflowing at all radii, with |  | >  vir ( ) outside the S zone and with |  | <  vir ( ) inside the S zone.On the other hand, the off-sheet gas is inflowing in the T zone and outflowing in the V zone, but always with |  | <  vir ( ).
parameter per-unit-length, in accordance with our modelling of the filaments as infinite cylinders.
We present the stacked profile of  vir () in Fig. 10 in black, along with the ratios of the four energies in the numerator of  vir (volumetric and surface kinetic and thermal energies) to the gravitational energy.Due to the inclusion of the surface terms, the virial parameter should be well-defined at each radius rather than just at the outer boundary of the system.Examining the profile of  vir () is useful to locate the radius where  vir ∼ 1, which may be considered the "virial radius" of the system.The profile monotonically decreases from  vir ∼ 10 in the S zone to  vir ∼ 1 in the T and outer V zones, 0.5 < ∼ / shock < ∼ 1.It then increases towards larger radii as the surface thermal term sharply declines outside the shock.We may thus associate a virial radius of  vir ∼  shock with our filament sample.
We further note that the volumetric and surface thermal terms are both ∼ 1 in the T zone, where the gas temperature is roughly equal to the effective temperature of dark matter defined by its velocity dispersion (Section 3.3).This implies that dark matter is also expected to be in approximate virial equilibrium (per-unit-length) within  shock .This is consistent with previous estimates of the virial radii of filaments feeding massive haloes at high- (Mandelker et al. 2018, equation 14).For a filament feeding a 10 12 M ⊙ halo at  ∼ 4, these authors predict  vir ∼ 55 kpc, slightly larger than our average  shock ∼ 45 kpc and within the range of  shock values among our ten slices (Fig. 4).Note that this is different from the case of galaxy clusters, where the shock radius is found to be larger than the virial radius (Lau et al. 2015;Zinger et al. 2018;Aung et al. 2021).
Removing  2  from T and T  reduces  vir by ∼ 30%, such that  vir ∼ 0.6 near  shock .This is consistent with the approximate configuration between the three dimensions of kinetic energy.Our inclusion of  2  in the calculation is thus justified because if small-scale motion along this dimension were restricted, the same kinetic energy would likely be split between the remaining two dimensions.On the other hand, including the axial surface terms in T  results in  vir < 0 at all radii.The exact numerical value of  vir , in this case, is of little consequence, but suffice to say that the system is clearly not virialized or even bound along its axis, as expected.Finally, regarding our velocity cut, excluding the gas with |  | > 2 vir does not affect our results, showing that the exclusion is dominated by a very rapidly flowing gas.Changing our threshold in   to only exclude inflowing gas with   < − vir while including rapidly outflowing gas retains  vir ∼ 1 in the T and outer V zones, showing that our estimation of a virial radius at  vir ∼  shock is robust.At smaller radii, where the system is still far from equilibrium based on Fig. 10,  vir becomes negative due to a much larger surface kinetic term.Including a rapidly flowing gas also results in negative  vir everywhere, consistent with our interpretation that the gas currently inflowing towards the filament along the sheet is not yet in equilibrium.
In summary, the gas within  shock , and by extension the dark matter as well, is in virial equilibrium within the potential well set by the dark matter filament, allowing us to define a virial radius  vir ∼  shock .Within the T zone, the gas is approximately at the virial temperature, and additional kinetic motions are largely offset by surface pressure.However, in the V and S zones, these kinetic motions grow stronger with respect to the gravitational potential, and the surface pressure terms cannot keep up.This drives the system out of equilibrium in these regions, consistent with the outflow velocities seen in Fig. 9.  15) to the gravitational energy, | W ( ) |. Specifically, these are the volumetric kinetic and thermal terms, T dyn (green solid line) and T therm (solid magenta line), and the surface kinetic and thermal terms, T dyn, (green dashed line) and T therm, (solid magenta line).We see that  vir , T therm /| W |, and T therm, /| W | are all of order unity in the T and outer V zones, where the gas temperature is hot and comparable to the effective temperature of the dark matter.Therefore, the filaments are in approximate virial equilibrium per unit-length with  vir ∼  shock .

DYNAMICAL EQUILIBRIUM
In this section, we study the dynamical stability of gaseous filaments.
In particular, we seek to address what the dominant force is that supports the filament against gravitational collapse towards its axis.We begin in Section 6.1 by describing the mathematical formalism of our force decomposition and present results from our simulation in Section 6.2.In Section 6.3, we interpret the simulation results using an analytic toy model for the internal filament dynamics.In Appendix Section C, we discuss the robustness of the results presented in Section 6.2 to our numerical method.

Basic Framework
As we show in Appendix B (Fig. B1), magnetic fields are not dynamically important in the filaments, with typical plasma  values of  =  thermal / magnetic > 10 5 .Therefore, the equation of motion for gas in our simulation can be approximated by the Euler equation: We follow Lau et al. (2013) in using eq.( 24) to determine the relative importance of the different terms in maintaining dynamical equilibrium and supporting the system against gravitational collapse.Specifically, we use a modified version of the "summation method" detailed in that paper.While that work focused on gas in massive galaxy clusters assuming spherical symmetry, we generalise their method to study gas dynamics in intergalactic filaments assuming cylindrical symmetry.We describe the method below, both for general considerations and for formulas specific to cylindrical geometry.

Primary Decomposition
Gauss's law relates the total mass in any arbitrary volume to the gravitational potential at the boundary of that volume, where  tot is the total mass enclosed by an imaginary closed surface  over which the surface integral is taken.Combining eq. ( 24) and eq.( 25), we have the following.
Each of the terms in the integral on the right-hand-side of eq. ( 26) represents an acceleration.When combined, these terms must balance gravitational acceleration.From the left-hand side of eq. ( 26), we learn that these terms can be considered as "mass terms", contributing to the total mass interior of the surface .Thus, we define where is the acceleration term representing the temporal change of the velocity field normal to the bounding surface,12 is the thermal term representing the support of the filament against gravity by thermal pressure gradients, and is the inertial term.As we demonstrate below and expand upon in Section 6.3, this term can be decomposed into several fictitious forces resulting from the fluid motion, each of which can either help support the filament against gravity or else work together with gravity to induce collapse.
In order to estimate the overall dynamical equilibrium of the filament and the relative contribution of each term to the support of the filament against gravity, we analyse eqs.( 28)-( 30) as a function of the radius of each filament slice and compare their sum (eq.27) to the true mass enclosed within the radius  in the simulation.For each radius , we assume a Gaussian surface that is a cylinder of radius  and length  = 30 kpc, the thickness of each filament slice.We comment on the impact of  on our results in Section 6.2.Unlike the discussion of virial equilibrium presented in Section 5, here we are interested in the dynamical state in three dimensions rather than per-unit-length.Furthermore, Gauss's theorem, in this case, requires a closed surface around the three-dimensional volume of interest.This results in three surface integrals for each radius, one on the cylindrical shell itself, hereafter  0 with d 0 = d 0 r, and two others on either end of the cylinder, hereafter  ± , with d ± = d ± ± ẑ.Then, we have the following: and where R stands for radial and A stands for axial.Each of the integrals in eqs.( 28)-( 30) is decomposed into two parts, Using this notation, we can explicitly write out all the terms in eqs.( 28)-( 30) as follows: Note that all the integrals in eqs.( 34)-( 39) are functions of radius, , which we have highlighted on the left-hand-side of each of these equations.Eqs. ( 34), (36), and (38) represent radial accelerations toward or away from the filament axis, while eqs.( 35), ( 37) and ( 39) represent axial accelerations along the filament axis.

Further decomposing the Inertial Terms
The last term in the integrand of eq. ( 38) represents the centrifugal acceleration and is hereafter referred to as the rotation term, When evaluating eq. ( 40), we further decompose the azimuthal velocity at radius  into mean and residual components.The mean rotation is given by (41) Note that eq. ( 41) represents the volume-weighted average azimuthal velocity at radius , rather than the mass-or density-weighted average.This is necessary for consistency with eq. ( 40).The residual azimuthal velocity in each cell is then given by  ,res (, , ) =   (, , ) −  ,mean ().
With these definitions, it is straightforward to show that where we have further decomposed the rotation term into a mean and a residual component, resulting from the mean and residual rotation velocity, respectively.
The first three terms in the integrand of eq. ( 38) can be interpreted in two ways.The first, which we will demonstrate in detail in Section 6.3, is the additional fictitious forces resulting from the fluid motion.The second contributes to ram pressure, turbulent pressure, and shear forces acting on the fluid.The ram (turbulent) pressure in the radial direction is the flux of radial momentum associated with mean (residual) radial motions, while radial shear forces represent flux of radial momentum caused by non-radial motions.These three terms in eq. ( 38) represent the radial momentum advected into a volume element by radial, azimuthal, or axial motions, respectively.If one defines mean and residual radial/axial motions analogously to eqs. ( 41)-( 42), such that the mean is only a function of , then it is straightforward to show that these terms can be written as follows: We expect the first terms in eqs.( 46) and ( 47) are negligible compared to the second, because  ,mean () and  ,mean () can be taken outside the integral and the average of the derivatives of   ,res should be small.Analogous divisions13 into the ram, turbulent and shear forces in the axial direction can be applied to the three terms in eq. ( 39).
Finally, we also divide the gravitational term (eq.27) into its radial and axial components: and These represent the radial and axial gravitational fields in the filament, respectively.Note that since the gravitational field is  = −∇Φ,  grav,r and  grav,a are defined to be positive when the radial and axial gravitational fields are directed inwardstowards the filament centre.However, the thermal and inertial terms in eqs.( 36)-( 39) and eqs.( 43)-( 47) are defined as positive when the respective forces are directed outwards, away from the filament centre.The equation of force equilibrium within the filament in the radial directions is with an analogous equation for the axial direction.

Numerical Calculation
Using our simulations, we numerically evaluate eqs.( 34)-( 49) as a function of radius, , for each filament slice.Full details of how this is done, including convergence and robustness tests, are provided in the Appendix Section C. Here, we briefly review the key points of our fiducial method.For each slice, we begin by depositing all fluid properties, including the gravitational potential, onto a uniform cylindrical grid extending from (0.05 − 1.1) shock with a typical cell-size of Δ ∼ 300 pc, comparable to the size of the smallest gas cells in the simulated filaments.We assign to each grid-cell the fluid properties of its nearest-neighbour gas cell in the simulation.
Despite not being inherently conservative, this method conserves mass, momentum, and energy to better than < ∼ 2% both globally and locally.All numerical derivatives and integrals are evaluated on these grids.Before computing any partial derivatives, we smooth the relevant quantity with a one-dimensional Gaussian in the direction along which we are evaluating the derivative, with a width of  = 0.5 kpc.Evaluating both sides of eq. ( 25), we find that our numerical estimate of  grav deviates from  tot by at most < ∼ 5% at each radius.After computing the radial profiles of each mass term for each filament slice, we normalise these by  tot (< ) and stack the ten filament slices by taking the average of each normalised mass term in each radial bin.While it is the ratio of each mass term to  grav () that tells us the contribution of the respective force to the filament support at that radius, we normalise the mass terms by  tot (< ), which is a smooth, monotonically increasing function within < ∼ 5% of  grav (), which exhibits small fluctuations.Hereafter we use the notation  ... / tot ,  ... / grav , and  ... / grav interchangeably and think of these as ratios of forces/accelerations rather than ratios of masses.
Finally, we note that we did not directly compute the two acceleration terms in eqs.( 34)-( 35), representing a temporal change in the velocity field in the filament.Direct evaluation of these terms requires the use of multiple snapshots.As detailed in Appendix Section C, this is difficult due both to the time between adjacent snapshots being comparable to both the cooling time and the eddy crossing time of filament gas, and to complications in identifying the same volume elements in each adjacent timestep while the entire filament moves through the sheet.Therefore, we use eq.( 50) and the corresponding equation for the axial terms to infer the acceleration terms from the difference between the gravitational, thermal and inertial terms.

Results
In what follows, we present the results of our analysis of the dynamical state of filaments.We begin in Section 6.2.1 by analysing where the velocity field in the filaments is in steady state, with  accel ∝ / ≃ 0, and where it is changing over time.We then analyse the various components of the thermal and inertial terms in Section 6.2.2, to better characterise the force balance and dynamical state in the filaments.Throughout, we emphasise key trends and physical processes in the three radial zones.

Overall Force Balance
In Fig. 11, we assess to what extent thermal and inertial forces balance gravity within the filament, as a function of / shock .The dashed dark-blue (hereafter navy) and pink lines show the stacked profiles of  grav,r / tot and  grav,a / tot , respectively, representing the normalized radial and axial gravitational fields.Notice that the radial term is positive, while the axial term is negative.As discussed following eqs.( 48)-( 49), this implies that the radial gravitational field is directed inwards towards the filament axis, as expected, while the axial gravitational field is directed outwards.The latter is a manifestation of the fact that filaments are unbound along their axis (see Section 5), and are being stretched apart by large-scale tidal forces.Note that the influence of the background Hubble expansion on this result is negligible.For reference, the dashed magenta line and the cyan-shaded region show the mean and 1- standard deviation among our ten slices of the radial profile of  grav / tot = ( grav,r +  grav,a )/ tot .This ratio is ≃ 1 everywhere, as Gauss's Theorem tells us, and must be (eq.25).Although our numerical method introduces an error ∼ 5% in / shock < ∼ 0.3 and > ∼ 0.8 (see also Fig. C3), this is negligible compared to the variations in  grav,r and  grav,a .
The solid navy and pink lines show the stacked profiles of ( therm,r +  inertial,r )/ tot and ( therm,a +  inertial,a )/ tot , respectively (eqs.36-39).In a steady state where thermal and inertial forces balance gravity and the filament velocity field is constant in time, these must equal  grav,r / tot (eq.50) and  grav,a / tot , respectively.If they are larger (smaller) than the corresponding terms  grav,r and  grav,a at some radius , then the velocity field at that radius accelerates outwards (inwards) in the corresponding direction.
In the radial direction, the filament exhibits a dynamical steadystate with no temporal acceleration throughout the inner halves of both the V zone and the T zone, 0.25 < ∼ / shock < ∼ 0.4 and 0.6 < ∼ / shock < ∼ 0.8.However, each zone contains regions with net temporal acceleration where the velocity field is not in a steady state.
In the T zone, at  > ∼ 0.8 shock , there is an outwards radial acceleration which appears to be due to the outwards expansion of the shock.While a full analysis of filament properties as a function of time is beyond the scope of this paper, we repeated the analysis presented in Section 3.1 to identify the shock radius in our ten filament slices at  ∼ 4.05 and  ∼ 3.82, roughly 55 Myr before and after our fiducial snapshot at  ∼ 3.93.The typical  shock increases from ∼ 36 kpc to ∼ 51 kpc during this interval, corresponding to a shock velocity of  shock ∼ 130 km s −1 , while the sound speed in the post-shock gas is  s ∼ 85 km s −1 .The outer T zone near  shock has not reached equilibrium after being hit by the shock, which explains the positive acceleration in this region.
In the V zone, there is an outwards radial acceleration at 0.4 < ∼ / shock < ∼ 0.6.This is the region where the radial velocity of the sheet material begins to decrease and the off-sheet material begins to move outwards (Fig. 9).As we show in Section 6.3 below, this is also the location of peak vorticity in the filament, where a quadrupolar vortex structure induced by the non-radial flow of the sheet towards the filament centre, dominates the dynamics.
In the S zone, at  < ∼ 0.25 shock , there is an inwards acceleration.This corresponds to the outer boundary of the isothermal core Force equilibrium in the filaments.We compare the sum of thermal and inertial forces (solid lines) with the gravitational force (dashed lines), in both the radial and axial directions, as a function of / shock .The magenta dashed line and cyan shaded region show  grav ( )/ tot (<  ), which is ≃ 1 everywhere with a < ∼ 5% error at / shock < ∼ 0.3 and > ∼ 0.8 due to our numerical method.The navy and pink dashed lines show the ratios  grav,r ( )/ tot (<  ) and  grav,a ( )/ tot (<  ), respectively.The former is positive, indicating a radial gravitational field directed inwards, while the latter is negative, indicating an axial gravitational field directed outwards since the filament is unbound along its axis.Solid navy and pink lines show (  therm,r + inertial,r )/ tot (<  ) and (  therm,a + inertial,a )/ tot (<  ), respectively (eqs.36-39).Radii, where the navy or the pink solid lines are equal to the corresponding dashed lines, indicate regions in a dynamical steady state without temporal acceleration in the corresponding direction (eqs.(34) and/or (35) are equal to 0).If, however, the solid lines are greater (smaller) than the corresponding dashed lines, this indicates a negative (positive) temporal acceleration at this radius.Over a fairly large radial range in the V and T zones, / shock ∼ (0.25 − 0.40) and (0.60 − 0.80), there is a dynamical steady-state in the radial direction.  / > 0 at / shock ∼ (0.40 − 0.60) and (0.80 − 1.0), while   / < 0 in the S zone at / shock < ∼ 0.25.In the axial direction, there is no dynamical steady state at / shock > ∼ 0.2.
(Fig. 3), where the thermal pressure decreases locally due to strong cooling (Fig. 8).However, we do not interpret this as evidence that the cold-stream is collapsing with no support against gravity.Rather, when examining the radial velocity profiles (as in Fig. 9) at  ∼ 4.05 and  ∼ 3.82, roughly 55 Myr before and after our fiducial snapshot at  ∼ 3.93, we see that the mild inwards velocity along the sheet in the S zone increases slightly with amplitude towards lower redshift, from < ∼ 0.05 vir ( shock ) to > ∼ 0.1 vir ( shock ).This may be due to additional gas accumulation near the filament axis, due to both inflows along the sheet and the off-sheet cooling flow.This results in a slightly larger inwards gravitational acceleration in the central regions, which leads to slightly more rapid inflow along the sheet.However, since the radial velocities in this region are always small anyway, we defer a more detailed study of the temporal evolution of the central streams to future work.
Finally, we note that in the axial direction, there is no steadystate at  > ∼ 0.2  shock , as the filament is continually stretched by tidal forces and the material accelerates towards the massive haloes on either end.Note that the acceleration term is positive in the S zone and negative outside it.This suggests that the velocity of the cold streams towards the haloes increases with time, whereas that of the hot filament gas decreases with time.We discuss this further below.

Detailed Force Decomposition
In Fig. 12, we show the individual profiles of all the forces that went into the curves shown in Fig. 11, to assess their individual contributions to the dynamical state of the filament.On the left, the solid lines show  therm,r / grav (red),  inertial,r / grav (green),  therm,a / grav (orange), and  inertial,a / grav (yellow).The dashed lines are as in Fig. 11, showing  grav,r / grav (navy),  grav,a / grav (pink), and  grav / tot (magenta).First, we note that the four nongravitational forces obtain both positive and negative values.As discussed in the following eqs.( 48)-( 49), a positive value implies an acceleration outwards that counteracts the collapse, while a negative value implies an acceleration inwards.In the centre and right panels, we show the profiles of all the individual components that comprise the radial and axial inertial terms, respectively (eqs.38-39).We begin with the radial direction, discussing the left and centre panels simultaneously, and focusing on the different radial zones.
In the T zone, at  > ∼ 0.6 shock ,  therm,r is directed outwards and is stronger than the gravitational force inwards, monotonically increasing from ∼ | grav,r | at  ∼ 0.6  shock to ∼ 4| grav,r | at  ∼  shock .This excess in the force of the external thermal pressure is partially balanced by a strong inertial force in the inside,  inertial,r , which increases in amplitude from ∼ (0.3 − 2)| grav,r | in the same radial range.The inwards inertial force in this region is dominated by the −    / term (centre panel, red lines), which increases in amplitude from ∼ (1.5 − 3)| grav,r |.As discussed in Section 6.1, this term contributes to the radial components of both the ram pressure force (eq.44, dot-dashed red line) and the turbulent pressure force (eq.45, solid red line).The ram pressure force can be visualised by examining the radial velocity profiles in Fig. 9.While the velocity in the sheet is roughly constant in the T zone, the inwards velocity outside the sheet decreases from ∼  vir at  shock to ∼ 0 at 0.6  shock , resulting in −    / < 0. In the outer T zone, 0.8 < ∼ / shock where there is a positive radial acceleration term (Fig. 11), the ram pressure force dominates over the turbulent pressure force, which actually becomes negligibly small at  > ∼  shock .In the inner T zone, 0.6 < ∼ / shock < ∼ 0.8 where the radial velocity field is in steady-state (Fig. 11), the ram pressure and turbulent pressure forces are comparable.The two shear forces,  and  (eqs.46-47, orange and blue lines in the centre panel, respectively), are dominated by residual motions (solid lines), while the shear generated by the mean flow (dot-dashed lines) is negligible.Together with the centrifugal forces (purple dot-dashed and solid lines for mean and residual rotation, respectively), these four terms combine to produce a net outwards force that is smaller than both the thermal and ram pressure forces throughout this region.In summary, strong thermal pressure gradients in the T zone support the filament against both gravity and ram pressure from the inflowing gas.At  < ∼ 0.8  shock , the gas is roughly in a steady state, while at  > ∼ 0.8  shock the velocity field accelerates outwards in response to the shock.
The V zone, 0.25 < ∼ / shock < ∼ 0.6, is the only zone where both  therm,r and  inertial,r are directed outward.In the outer V zone, 0.4 < ∼ / shock < ∼ 0.6 where there is a positive radial acceleration (Fig. 11), the thermal pressure is dominant and still stronger than the radial gravitational force and is comparable in amplitude to its value in the inner T zone.In the inner V zone, 0.25 < ∼ / shock < ∼ 0.4 where We compare the radial thermal and inertial forces (solid red and green lines) with the radial gravitational force (navy dashed line, as in Fig. 11), and the axial thermal and inertial forces (solid orange and yellow lines) with the axial gravitational force (pink dashed line).Centre panel: We show the four constituent forces comprising the radial inertial force (eq.38), namely the centrifugal force (purple lines), the ram/turbulent pressure forces (red lines) and the shear forces (orange and blue lines for azimuthal and axial shear, respectively).For each force, we further distinguish between their mean (dot-dashed lines) and residual (solid lines) components.The legend for this panel is in the box located below all of the panels.

Right panel:
We show the three constituent forces that comprise the axial inertial force (eq.39).These are the ram/turbulent pressure forces (blue lines) and the two shear forces (red and yellow lines for radial and azimuthal shear, respectively).Based on all three panels, we conclude: In the T Zone, the radial thermal force is directed outwards and is stronger than the gravitational force inwards, while the inertial forces are directed inwards and are dominated by the ram pressure, which mainly offsets the excess thermal pressure.In the V Zone, the radial thermal force is directed outwards and decreases from ∼ 1.5| grav,r | to ∼ 0.6| grav,r | as  decreases.Radial inertial forces are directed outwards with a roughly constant amplitude of ∼ 0.6| grav,r |, and are dominated by centrifugal forces.The mean rotation has a roughly constant amplitude of ∼ | grav,r |.The residual rotation is much larger, but is largely offset by ram/turbulent pressure and azimuthal shear.In the S Zone, the radial thermal and inertial forces become negative at  > ∼ 0.15  shock and  < ∼ 0.15  shock , respectively, resulting in a net inwards acceleration in this region.The axial gravitational force is directed outwards, representing the tidal stretching of the filament along its axis.It is roughly balanced by axial ram pressure forces at all radii, while the axial thermal pressure force is relatively small, directed outwards in the S zone and inwards the T and V zones.
the radial velocity field is in steady-state (Fig. 11), the thermal pressure term decreases and becomes comparable to the inertial term, which maintains a value of < ∼ 0.6| grav,r | throughout the entire V zone.The change in sign of the inertial term, from negative in the T zone to positive in the outer V zone, is mainly due to a decrease in the force of ram pressure inwards and an increase in the centrifugal force outward.The mean rotation term increases from ∼ 0 at  ∼  shock to a constant ∼ | grav,r | throughout the S and inner V zones, at  < ∼ 0.4 shock .The residual rotation term is ∼ 1.5| grav,r | at  ∼  shock , and increases monotonically toward smaller .Throughout the V and S zones, at  < 0.6 shock , this is mainly offset by the inwards   shear force, −(  /)  /, with a smaller contribution from the ram/turbulent pressure.In Section 6.3, we will show that the behaviour of this trio of forces is generic, resulting from the quadrupolar vortex structure that dominates the filament dynamics throughout the V zone.To summarise the situation in the outer V zone, at 0.4  shock < ∼  < ∼ 0.6  shock , strong thermal pressure forces combined with increasing centrifugal forces and decreasing ram pressure forces lead to a net outwards acceleration of the velocity field.The lack of a steady state here may be linked to an evolving vorticity structure fuelled by ongoing accretion onto the inner filament (see Section 6.3 below).
In the inner V zone, at 0.25 shock < ∼  < ∼ 0.4 shock where the filament velocity field is in steady state with no acceleration (Fig. 11), the thermal force continues to be directed outwards but is now weaker than the gravitational force inwards, ∼ (0.3 − 0.6)| grav,r |.The inertial term continues to be directed outward, maintaining a constant amplitude of ∼ 0.6| grav,r |.The components of the inertial term behave similarly here to the outer V zone.The outwards inertial force is dominated by the centrifugal force, with the mean rotation term roughly balancing gravity with an amplitude of ∼ | grav,r |.The residual rotation term continues to increase towards smaller , with amplitudes in the range ∼ (4 − 7)| grav,r |, though this is again compensated for by shear and turbulent pressure (see Section 6.3 below).This trio yields a slight net inwards acceleration that lowers the total inertial term compared to the mean rotation.To summarise the situation in the inner V zone, at 0.25  shock < ∼  < ∼ 0.4shock, the filament remains in force equilibrium with roughly equal contributions from thermal pressure forces and centrifugal forces in the support against gravity.
In the S zone, at  < ∼ 0.25 shock , there is a net inwards acceleration term (Fig. 11).This is mainly driven by the thermal pressure force becoming negative at  ∼ (0.15−0.25) shock , near the outer boundary of the isothermal core, where the cooling is maximal.The inertial term becomes negative at  < ∼ 0.15 shock due to a slight drop in centrifugal and residual shear forces (dot-dashed purple and solid blue lines).This suggests that the dynamics of the isothermal cores of high- filaments, representing the cold streams feeding massive high- galaxies (Section 3.1), are far from a relaxed state of equilibrium, but are dominated by strong cooling flows and complex motions and accelerations of infalling and strongly rotating gas.
We broadly summarise the balance of the radial force throughout the filament, highlighting the key physical characteristics of each of the three radial zones.In the T (thermal) zone, thermal pressure forces dominate, roughly balancing gravity and ram pressure forces.In the V (vortex) zone, the thermal pressure forces become comparable to the centrifugal forces induced by mean rotation, and these combine to roughly balance gravity.At the same time, the dominant forces are actually residual rotation, shear, and turbulent pressure, which are induced by the vortex structure of the filament.Individually, each of these is stronger than gravity in absolute value, though they roughly cancel out, leaving a relatively small combined force, which is weaker than the mean centrifugal term.In the S (stream) zone, a drop in both centrifugal and thermal pressure forces due to the development of a strong cooling flow results in a net inwards acceleration of the velocity field.
In the axial direction, the situation is much simpler.The axial inertial force is comparable in amplitude and opposite in sign to the axial gravitational force, which both exhibit a very weak radial dependence. inertial,a is dominated by axial ram pressure forces encapsulated in the term −    /, while the other two terms, representing turbulent forces, are very small.The axial thermal pressure force is comparatively small in amplitude, directed outwards in the S zone and inwards at larger radii.This is consistent with a picture where the cold stream gas penetrates the hot CGM and free-falls towards the massive central galaxies on either end, while the hot filament gas does not penetrate the virial shock and instead builds up in the CGM and IGM around the haloes, providing additional pressure confinement that slightly weakens the tidal stretching of the filament in the axial direction.We will discuss this further in Section 7.2.

A Toy Model for the Filament Velocity Field
In Section 6.2, we saw that, in addition to thermal pressure gradients, the filament inertia plays a major role in its overall equilibrium.In particular, the residual rotation provides a very strong outwards force that significantly exceeds the inwards gravitational force at small radii, and is largely balanced by the shear and ram/turbulent forces, (  /)  / and     /.In this section, we wish to understand the origin of these terms and in particular how they represent fictitious forces acting on the gas along its streamlines, rather than our description of them as ram/turbulent/shear forces in Section 6.2.We begin in Section 6.3.1 by discussing the very simple case of an elliptical orbit in a Keplerian potential as an instructive example.In Section 6.3.2, we discuss the vorticity structure of the filaments in our simulation, highlighting a characteristic quadrupolar structure.Finally, in Section 6.3.3,we discuss how such dynamics can give rise to the inertia forces we find.

The Keplerian Orbit
Consider an elliptical Keplerian orbit of a test particle around a central mass .The equation of motion in the radial direction, with respect to the central mass, is where  is the distance from the central mass and  is the azimuthal angle.The azimuthal and radial velocities are   =  , and   = .
The first term on the left-hand side of eq. ( 51) is thus  2  / while the second is   /.The latter term can be expressed in two possible ways, depending on whether the orbit is parametrised with respect to the radius, , or the azimuthal angle14 , .
If we parametrize the orbit in terms of  we have Inserting this into eq.( 51) yields showing that in the frame of the particle, gravity is indeed balanced by the two relevant components of the inertial term in eq. ( 38), which act as fictitious forces along the orbit.We can explicitly evaluate eq. ( 53) for the Keplerian orbit by recalling that both the specific angular momentum,  =   , and the specific energy, E = − /(2) with  the semi-major axis, are conserved.Therefore, and One can easily verify that eqs.( 54)-( 55) obey eq. ( 53).
If we parametrize the orbit in terms of , eqs.( 52)-( 55) become where  is the eccentricity of the orbit, which is a conserved quantity.Again, one can easily verify that eqs.( 58)-( 59) obey eq. ( 57) and that gravity is balanced by the two relevant components of the inertia term in eq. ( 38), which act as fictitious forces along the orbit.
This simple example shows that in a steady flow around a central mass without thermal pressure, the three major terms appearing in eq. ( 38) act as fictitious forces that balance gravity.

The Velocity Structure of Filaments
In order to model the inertia forces seen in Fig. 12 using an analysis similar to that presented in Section 6.3.1, we must first characterise the gas velocity field within the filaments.In Fig. 13, we show the velocity field in slice 2 (see Fig. 1), projected along the filament cross-section and overlaid on a map of the average density along the filament's axis, in the same orientation as the stacked projections shown in Fig. 3.The arrows represent the 2D gas velocity field perpendicular to the filament's axis,   x +   ŷ, mass-weightedaveraged along the axis.Several features of the velocity field are apparent.Outside the shock radius (dashed-white circle), gas flows towards the filament from both inside and outside the sheet (aligned parallel to the -axis).The cylindrical shock is evident in the flow pattern of off-sheet material whose velocity suddenly changes orientation and develops vorticity, as expected from a curved shock.The sheet gas, on the other hand, maintains its approximate trajectory due to its higher density, which destabilises the shock there (see Dekel & Birnboim 2006, for an analogous discussion of cold streams in hot haloes).Inside  shock , the velocity is broadly characterised by four vortices centred near the transition between the T and V zones at  ∼ 0.6  shock , in a quadrupolar pattern with the upper-right and lower-left quadrants spinning counter-clockwise (positive vorticity), and the upper-left and lower-right quadrants spinning clockwise (negative vorticity).Such a quadrupolar vorticity structure has been seen in previous studies of intergalactic filaments using lower-resolution cosmological simulations (e.g., Pichon & Bernardeau 1999;Codis et al. 2012Codis et al. , 2015b;;Laigle et al. 2015;Xia et al. 2021;Ramsøy et al. 2021), and is thought to arise from asymmetric inflow into filaments from sheets with non-zero impact parameters.While this is the only source of vorticity for dark matter, an additional source of vorticity for the gas is shear due to the inflow along the sheet interacting with the hot gas in the T and V zones, combined with the vorticity creation due to the curved shocks.We defer a more detailed study of the origin of this vorticity structure in our simulations to future work.
The quadrupolar vortex structure is a generic feature of filament dynamics in our simulation, not unique to slice 2. In Fig. 14 (left), we show a stacked projection map of   , the vorticity component parallel to the filament axis.In each slice, we compute the mass-weighted average value of   along the line-of-sight, and weigh this by  circ ( shock )/ shock prior to stacking, where In the right-hand panel of Fig. 14, we address the fraction of vorticity aligned with the filament axis.For each slice, we compute the mass-weighted average of both  2  and  2 tot = ( 2  +  2  +  2  ) at each radius , and take the square-root of their ratio.We then present the mean and standard deviation of this ratio among the ten slices.If the vorticity were randomly oriented, this ratio would be ∼ 0.5, which is the average of cos() uniformly distributed between 0 and 1.Instead, the typical ratio is > ∼ 0.6 within the shock and increases slightly toward smaller radii, displaying a preference for the vorticity to align with the filament axis.This is broadly consistent with the results of Laigle et al. (2015) who found a preference for the filament vorticity to align with its axis, with an excess probability of ∼ 20% having an angle less than 60 • , corresponding to a cosine greater than 0.5.

Ideal Vortex Model
While the gas dynamics in filaments is inherently three-dimensional and complex, the quadrupolar vortex structure seen in Figs.13-14 appears to be the most robust feature common to all filament slices.Therefore, we examine to what extent such a velocity field alone can reproduce the robust features of the radial inertial forces seen in Fig. 12 (centre panel).To this end, we construct a simple twodimensional toy model for the filament velocity field, characterised by four ideal vortices in a quadrupolar configuration, such that the vorticity is positive (negative) in the first and third (second and fourth) quadrants, following Figs.13-14.This model is shown in Fig. 15.In the frame of each vortex (hereafter the primed frames), its velocity contribution is given by where the plus sign applies to the first and third quadrants and the minus sign applies to the second and fourth quadrants. sets the slope of the rotation curve with respect to the vortex centres while  is a normalisation constant. ′ is the distance to the centre of each vortex, which is assumed to be at (, ) = (±, ±).
The vorticity corresponding to such a velocity field, in the vortex frame, is: For any  ≠ −1, the vorticity is non-zero 16 and finite at all  ′ > 0.
Transforming the velocity field to the filament frame (hereafter the unprimed frame) yields where in every plus/minus combination, the upper symbol applies to the first and third quadrants, and the lower symbol applies to the second and fourth quadrants.The transformation of  ′ and  ′ The stacked profile of the contribution to the filament vorticity oriented along its axis.We calculate the mass-weighted average profiles of  2  and  2 tot for each slice and take the square root of their ratio.We then compute the average and standard deviation of this ratio among our ten slices.Filaments prefer vorticity to align with their axis, with a typical ratio of > ∼ 0.6, compared to 0.5, which is characteristic of a random distribution.A visual inspection of Fig. 13 indicates that, the rotation velocity increases away from the vortex centres.Therefore, we take  > 0. We further require that the ram/turbulent pressure term, −    / be negative, based on the middle panel of Fig. 12, which turns out to require  < 1.We use  = 0.5 as our fiducial value, but note that we obtain similar trends with any value of  ∼ (0.1 − 0.9).The value of the normalisation constant  is arbitrary, and hereafter we take  = 1.
Given parameters , , and , we compute the velocity field and the resulting fictitious forces using the same method we used in Section 6.1.We set up a cylindrical grid with linear spacing in the radial and azimuthal directions and evaluate the velocity in each bin.Note that since our assumed velocity field is explicitly twodimensional, all axial derivatives vanish, i.e., / = 0. Furthermore, in our highly idealised model, both the mean azimuthal and radial velocities are zero,  ,mean =   ,mean = 0.This leaves only the three components in the radial inertial term, which depend only on the residual velocities, namely the turbulent pressure (eq.45), residual azimuthal shear (eq.46), and residual centrifugal force (eq.43) (red, orange, and purple solid lines in the centre panel of Fig. 12).We evaluate these following our methodology described in Appendix Section C. Finally, we normalise all forces by (/)  , with  = 1.48.This accounts for the radial dependence of the gravitational force,  15).Since our model is completely symmetric, the mean parts of   and   are, by definition, 0, and therefore only the residual parts contribute.Each force term of our model has been scaled by an additional   to account for the radial dependence of  g , the denominator in the ratios presented in Fig. 12.Our simple model reproduces many of the trends seen in the simulation in the V zone, at 0.25 < ∼ / shock < ∼ 0.6 (Fig. 12, middle panel, solid lines), where the ram pressure force has decreased and vortices dominate the dynamics of the filament (see the text for details).
g , the denominator in all the profiles shown in Fig. 12, which is given by the total mass interior to radius  (blue line in Fig. 7).
The resulting force terms are plotted in Fig. 16.These reproduce many of the trends seen the corresponding curves in the centre panel of Fig. 12.At small radii ( < ∼ 0.5), the residual rotation is mostly balanced by azimuthal shear, both of which appear to diverge as  → 0 as  − .The turbulent pressure also contributes to the balance of the rotation term, but saturates at a finite value as  → 0.Even the sign-flip in this term at  < ∼ 0.6 , near the singularity present in our model at the location of the vortices, resembles the sharp decrease in the amplitude of the corresponding term in Fig. 12 near the location of the vortices at  ∼ 0.6  shock .The corresponding sign flip of the shear force at  =  is not seen in the data.However, in the simulations, the inertial forces in the T zone at  > 0.6  shock are dominated by the ram pressure of the inflowing gas, which is not included in our ideal vortex model.Therefore, it is expected that it will only represent the data at  <  = 0.6 .Furthermore, unlike our simple two-dimensional toy model, the actual filament dynamics is complex and three-dimensional, with vorticity that is not even truly aligned with the filament axis (Fig. 14).Despite these caveats, our model confirms that the structure of the quadrupolar vorticity in intergalactic filaments produces fictitious force terms that resemble the inertial terms dominating the dynamics of the filament at  < ∼ 0.6  shock .More detailed theoretical studies of cosmic velocity fields in filaments can be found in, e.g., Pichon & Bernardeau (1999)

Filament Thermal Structure and Size
We compare our results with those of Ramsøy et al. (2021), who used an AMR simulation to study the properties of an intergalactic filament that feeds the progenitor of a MW-mass halo ( v ∼ 10 11.5 M ⊙ at  = 0).While they examined the filament in the redshift interval  ∼ (3.5 − 8), their main analysis focused on  ∼ 4 as ours, and many of their results are similar to ours.They too find that the filament is surrounded by a cylindrical accretion shock and embedded in an intergalactic sheet with a planar accretion shock.The filament density in their simulation is well fit by the profile of an infinite, selfgravitating, isothermal cylinder, out to the radius where the sheet begins dominating the mass distribution.This was true for both gas and dark matter.Finally, similar to our results, their filaments were characterised by a quadrupolar vorticity structure with a weak radial dependence inside the filament shock radius.Despite these similarities, there are several important differences between our results.In their simulation, the density profiles of gas and dark matter were extremely similar.In particular, the characteristic scale radii,  0 , for the gas and dark matter were within a few tens of percent of each other, while these differed by a factor of > ∼ 3 in our simulations.Similarly, the effective dark matter temperature in Ramsøy et al. (2021) was only a factor of < ∼ 2 higher than the central gas temperature (see their figure 7), while it was ∼ 30 times higher in our simulations (Section 3.3).Furthermore, the shock surrounding the filament in their simulation was weaker and narrower than in ours.They measured an azimuthally-averaged gas-temperature increase of < ∼ 60% between the filament core and the shock, while the temperature in our simulation increases by a factor of ∼ (20 − 30).Likewise, they found the shock to be nearly isothermal, with the gas temperature reaching its core value ∼ 3 kpc from the shock front, while in our case the gas remains hot outside the cooling radius of ∼ 0.5  shock ∼ 25 kpc (Section 4) and only reached the core temperature ∼ (30 − 40) kpc from the shock front.Finally, while they find that thermal pressure alone can account for most of the filament support against gravity, we find this to be true only in the T zone near the shock radius, while the V and S zones are dominated by vortical, turbulent, and rotational motions driving the filament far from equilibrium.
The primary reasons for these differences are likely the different environments and line-masses of the filaments studied in this work compared to that studied in Ramsøy et al. (2021).While these authors studied a single, relatively low-mass and isolated filament feeding a single MW-progenitor halo,  v ∼ 10 10.5 M ⊙ at  ∼ 4, we analyse ten slices of three filaments in a crowded region between three haloes with  v ∼ 10 12 M ⊙ at  ∼ 4. In general, more massive filaments are predicted to have more prominent accretion shocks, similar to the transition from cold to hot CGM for haloes (Birnboim et al. 2016).More quantitatively, the shock Mach number is predicted to scale as M ∝  r ∝  vir ∝ Λ 0.5 (eq.23, see also Appendix A).Based on Eq. ( 10) from Mandelker et al. (2018), the filament line-mass scales with the halo mass as Λ ∝  0.77 v , implying that our filament is < ∼ 15 times more massive (per-unit-length) than theirs.For a filament in virial equilibrium per-unit-length (see Section 5), the effective DM temperature is expected to scale as  ∝ Λ, implying that the DM temperature in our filaments should be < ∼ 15 times larger than in Ramsøy et al. (2021).On the other hand, the central gas temperature is ∼ (2 − 3) × 10 4 K in both cases, set by thermal equilibrium with UVB, as the central gas can efficiently cool.We thus expect the ratio of DM to gas temperature to be ∼ 15 times larger in our simulations compared to Ramsøy et al. (2021), consistent with the simulation results cited above.This demonstrates how filament properties are sensitive to the masses of the haloes they feed, and validates several scaling relations.
Another important difference between our studies is the resolution in the filaments.In the simulations analysed by Ramsøy et al. (2021), the filaments were resolved by ∼ 1.2 kpc with small patches surrounding haloes embedded in the filament resolved with ∼ 600 pc.In our simulations, the typical cell size is ∼ (300 − 500) pc in the S zone and (0.5 − 1.0) kpc in the T zone (Fig. C1).Combined with the fact that our filaments are ∼ 3 times wider than the filament studied in Ramsøy et al. (2021), we have nearly an order of magnitude more cells across the filament diameter in our simulation, allowing us to better resolve the turbulent cascade and reduce artificial dissipation of turbulent energy.Finally, it is possible that some of the differences are driven by differences in the hydro-solver (Eulerian AMR versus moving-mesh) and the different sub-grid models implemented in the simulation.More work with larger samples of filaments in different simulations using different codes will be required to break these degeneracies.
Another interesting question is why the gas density profile is wellfit by an isothermal model with the correct core temperature, despite the temperature increasing by a factor of ∼ 20 from the core to the shock.If we assume that the gas initially traced the DM and was thus roughly isothermal at the post-shock temperature, then the gas density is initially cored with a comparable radius to the DM core, as in Ramsøy et al. (2021).This constant-density and constant-temperature core then begins to cool monolithically toward a temperature of ∼ 3×10 4 , set by thermal equilibrium with the UVB.If this cooling is isobaric, as suggested by the gas pressure profile (Fig. 3), then the density in the core increases by a factor of ∼ (15−20) as the gas cools.From the mass conservation in the core, we thus expect the core radius to shrink by a factor of ∼ (15 − 20) 1/2 ∼ 4, consistent with the ratio of the DM core radius to that of the gas.On the other hand, at large radii, near  shock , the gas density is much smaller.Cooling, therefore, is much less efficient, so we may assume that the density profile of the gas at these large distances continues to trace the DM, and thus asymptotically approaches the slope  ∝  −4 of an isothermal cylinder irrespective of the central temperature.This will produce a density profile consistent with Fig. 3 despite the gas not being isothermal.
We can use these insights to derive a prediction for the size of the isothermal filament core (the S zone) as a function of redshift and the mass of the halo that the filament is feeding.From Mandelker et al. (2018), the line-mass and virial radius of filaments are given by (their equations 10 and 14) where the inflow velocity of the filament towards the halo in units of the halo virial velocity, and  s,3 =  s /(1/3) ∼ 1 is the fraction of the total accretion onto the halo flowing along a given filament normalized by a fiducial value of 1/3.These predicted values are only a few tens of percent larger than those found in our simulation, in Section 3.4 and Section 5 respectively.For a filament in virial equilibrium perunit-length, eq. ( 9) relates the line-mass to the kinetic energy perunit-length.The latter is related to the virial temperature through eq. ( 5), K ∼ 1.5 B  v Λ gas /( p ), with Λ gas ∼  b Λ fil the line-mass of gas within the filament.Taken together with eq. ( 63), we obtain an expression for the filament virial temperature, Compared to the virial temperature of dark matter haloes (e.g., Dekel & Birnboim 2006),  v ∼ 2.5×10 6 K  2/3 12 (1+) 5 , the virial temperature of filaments increases more rapidly with redshift for a given halo mass.This is consistent with the results of Birnboim et al. (2016), who predicted more stable shocks in filaments feeding a given halo mass at higher redshifts.
We now assume that the initial gas core is comparable to the DM core, which is comparable to  v,fil (Section 3.3), and that this core cools isobarically from  v,fil to  s ∼ 2 × 10 4 K.The core density thus increases by a factor of  v,fil / s , and thus mass conservation dictates that the core radius is roughly with no dependence on halo mass,  s or M.This is consistent with the fact that the core radius in our simulation of  s ∼ 10 kpc is similar to the core radius in Ramsøy et al. (2021) despite the factor ∼ 30 difference in the halo mass.Associating the filament core (the S zone) with the cold stream that penetrates galaxy haloes, this model predicts a stream radius that scales differently from the predictions of both Mandelker et al. (2018), who assumed contraction from  v,fil until full angular momentum support with a constant spin parameter and obtained  s ∝  0.38 v (1 + ) −0.5 (their equation 23), and Mandelker et al. (2020b), who assumed thermal pressure equilibrium between the cold stream and the hot CGM at the halo virial temperature and obtained  s ∝ (1 + ) −1 independent of the halo mass (their equations 21 and 27).While all three of these models roughly agree on the radius of streams feeding haloes of  v ∼ 10 12 M ⊙ at  ∼ 4, future studies of the properties of cold streams as a function of halo mass and redshift will be able to distinguish these models and shed further light on what determines the size of cold streams.

Filament Dynamics
One of our main findings was the presence of a strong quadrupolar vorticity structure in the filaments, which dominates the filament dynamics and inertial forces in the S and V zones, at  < ∼ 0.6 shock (Fig. 12 and Section 6.3).Similar quadrupolar vorticity structures in intergalactic filaments have been found in several previous works using a variety of simulation and analysis methods (e.g., Pichon & Bernardeau 1999;Pichon et al. 2011;Codis et al. 2012Codis et al. , 2015b;  67) and ( 68), respectively, which are both reasonable matches to the simulation data.Laigle et al. 2015;Ramsøy et al. 2021;Xia et al. 2021).This vortical structure is thought to result from anisotropic accretion onto filaments from surrounding sheets, either due to a non-zero impact parameter or to shearing flow through a curved shock.While this seems consistent with our simulations, we defer a more detailed study of the origin and evolution of vorticity to future work.
Perhaps more intriguing, several recent works have claimed that filaments have a net coherent rotation or spin, giving rise to a helical flow along the filament, which can dominate over the quadrupolar vorticity structure discussed above.This has been claimed both using simulations (Xia et al. 2021) and stacked observations at low- (Wang et al. 2021).The filaments in our simulations show evidence for mean rotation which actually seems to be the primary support against gravity in the inner regions (Fig. 12).However, the strength of the mean rotation in our simulations is much weaker than the "residual" rotation resulting from the quadrupolar vorticity field.Future work with larger samples of filaments of different masses and across different environments and redshifts will be required to study the evolution of filament spin.
Finally, it is worth commenting on the similarities and differences between our force analysis presented in Section 6 and the recent force analysis of CGM gas in < ∼ 10 12 M ⊙ haloes at  ∼ 0 from the FOG-GIE simulations (Lochhaas et al. 2023).While our mathematical formalism, based on previous works focusing on the ICM (e.g., Lau et al. 2013), derives directly from the Euler equation and is hence exact, the interpretation of some of the inertial terms can be rather abstract (see the discussion in Section 6.1 and Section 6.2).On the other hand, Lochhaas et al. (2023) focused on five effective forces operating on the gas, which are all intuitive: gradients of thermal, turbulent, and ram pressure, centrifugal forces, and gravity.This results in a much more intuitive picture, as all the forces are well-defined and have clear physical meanings.However, while their thermal pressure gradients and gravitational forces have a one-to-one correspondence with ours, the other three forces are explicitly tied to a smoothing scale that separates mean from residual motions and serves only as approximations of the full shear tensor, which is accurately captured by the inertial forces in our method.Therefore, the exact force balance is not guaranteed in their method, even in the absence of temporal accelerations.Despite these differences and the very different systems studied (high- intergalactic filaments versus low- CGM), many of our results are similar.In particular, we both find that the outskirts of the system are primarily supported against gravity by thermal pressure gradients, while kinematic (inertial) forces dominate the inner regions and are not in a steady-state equilibrium.Moreover, we both find that locally, turbulent and ram pressure forces can act either outwards (opposing gravity) or inwards (aiding gravity).A more detailed comparison between these two mathematical methods of force decomposition is left for future work.This highlights the similarity between gaseous atmospheres in different cosmic-web elements and makes it tempting to describe the post-shock region in our filaments as a 'circumfilamentary medium' (CFM) surrounding the central cold stream within the potential well of the dark matter filament, much like the CGM surrounding the central galaxy within the potential well of the dark matter halo.

Implications for Cold Streams
The S zone of filaments is a dense, cold, isothermal core, where the density and temperature are roughly  H ∼ 10 −2 cm −3 and  ∼ 2 × 10 4 K (Section 3.1, Figs. 3 -4).This core represents the 'cold streams', predicted to be the main mode of gas accretion onto massive high- galaxies and to penetrate the virial accretion shocks around their dark matter haloes (Dekel & Birnboim 2006;Dekel et al. 2009a).As described above, the radius, density, and temperature of this region are consistent with predictions from various analytic models for the properties of cold streams as a function of the halo mass and redshift.
Several recent studies have attempted to model the evolution of cold streams as they travel towards massive central galaxies while interacting with the hot CGM, subject to various (magneto)hydrodynamic, thermal, and gravitational instabilities using analytic models and idealised numerical simulations (Mandelker et al. 2016(Mandelker et al. , 2018(Mandelker et al. , 2019a(Mandelker et al. , 2020a,b;,b;Padnos et al. 2018;Aung et al. 2019;Berlok & Pfrommer 2019).While these models are adding more and more physics and thus becoming more realistic, all of them have modelled the streams as initially laminar, hydrostatic systems supported predominantly by thermal pressure.While Mandelker et al. (2018) did account for mild turbulence with Mach number M turb < ∼ 1 driven by accretion onto the filament from the surrounding sheet, this is still a fundamentally different picture from the highly turbulent environment of cold streams in our simulations, which is dominated by vorticity and accretion and is surrounded by a cylindrical accretion shock.These must be accounted for in future work studying the stability and evolution of intergalactic filaments and cold streams.In particular, future models should account for internal filament dynamics such as rotation and turbulence, and for the confining ram pressure from the accreting gas, along with turbulence in the confining medium.
Finally, while the cold streams are dominated by non-thermal, turbulent, and vortical motions, they are in pressure equilibrium with the post-shock volume-filling CFM gas.This implies that as the streams penetrate the hot CGM of massive dark matter haloes, where both the densities and temperatures are larger, their confining pressure will increase by a factor of several.Indeed, we find in our simulations that the confining pressure of filaments jumps by a factor of ∼ (10 − 20) as they penetrate the virial shocks around > ∼ 10 12 M ⊙ haloes at  ∼ 4 (Lu et al., in prep.).First, it is unclear what happens to the hot CFM gas as the cold stream penetrates the halo.Does it stay behind and accumulate at the virial radius, or does it partially penetrate as well, dragged along by the cold stream?Second, this can have important implications for the morphology of cold gas in the outer CGM of such haloes, since such a sudden increase in confining pressure may cause the streams to 'shatter' into tiny fragments of order of the cooling length,  cool ∼  s  cool , with  s and  cool the sound speed and cooling time of gas with  > ∼ 10 4 K (e.g., McCourt et al. 2018;Gronke & Oh 2020, 2022).Similar phenomena have been predicted by other models of streams that penetrate virial shocks, although for different reasons (Cornuault et al. 2018).If the streams do shatter as they penetrate the virial shock, this may explain puzzling observations of large area covering fractions and small volume filling fractions of cold gas in the CGM of massive high- galaxies, which are consistent with a shattered mist-like collection of small cloudlets (e.g., Borisova et al. 2016;Cantalupo et al. 2019;Pezzulli & Cantalupo 2019).The process of how cold streams first penetrate the hot CGM around massive galaxies remains poorly understood and will be the subject of future work using idealised simulations of stream-shock interactions and shattering of filamentary systems (Yao, Mandelker, et al., in preparation).

SUMMARY AND CONCLUSIONS
Using a novel high-resolution cosmological simulation, IPMSim, that zooms in on a large patch of the cosmic-web in a large-scale protogroup environment, we study the structure and internal dynamics of sections of ∼ Mpc-scale intergalactic filaments feeding three ∼ 10 12 M ⊙ haloes at  ∼ 4. We select ten slices 30 kpc in length from three independent filaments that feed these three haloes and are embedded in the same cosmic sheet, so that no slice intersects any halo with  v > 10 10 M ⊙ .We then study these slices independently and in stacked form to reach the following conclusions: (i) Radial zones: The filaments can be broadly characterised by three radial zones (Fig. 2).We summarise the key properties of these zones here, while further below we offer a more detailed summary of each property over the entire filament.
• In the outer "thermal" (T) zone, a cylindrical accretion shock heats the gas to the virial temperature associated with the potential well of the dark-matter filament, and strong thermal pressure forces balance both inwardsgravitational and ram pressure forces resulting from the inflowing gas.Near  shock , excess thermal pressure forces cause the shock to expand, though the gas continues to flow with a roughly constant mass accretion rate throughout the T zone, both on and off the sheet, due to cooling of the gas behind the shock.The gas outside the sheet decelerates throughout this zone, while the sheet gas maintains a roughly constant inflow velocity.
• In the intermediate "vortex" (V) zone, the filament dynamics is dominated by a quadrupolar vortex structure in the velocity field, induced by the flow of gas through the sheet towards the central filament.The outwards centrifugal force due to the global net rotation is already comparable to the gravitational force inward, but it is small compared to the centrifugal forces due to the residual rotations associated with the quadrupolar vortices.The latter are largely balanced by shear and turbulent forces.These vortices lead to net outflows outside the sheet, and the induced mixing causes the gas inflowing along the sheet to decelerate.
• In the inner "stream" (S) zone, a dense isothermal core forms from an isobaric cooling-flow resulting from post-shock cooling in a freefall time, and is associated with a decrease in outwards forces.The size of this region expands with time, though the gas within this zone can move both inwards and outwards due to the impact parameter of the sheet and the resulting vorticity.The S zone represents the cold streams that feed massive galaxies from the cosmic web.
While this basic structure appears generic, the locations of the boundaries between these zones depend on filament line-mass (Section 7.1).In our case the T zone is at  > ∼ 0.65 shock and the S zone is at  < ∼ 0.25 shock .(ii) Thermal structure: The gas density profiles are well fit by the density profile of a self-gravitating isothermal filament with a temperature of ∼ 3×10 4 K.However, the filaments are not isothermal, and the temperature increases by a factor of ∼ (20 − 30) from the core in the S zone to the cylindrical accretion shock outside the T zone, at  shock ∼ 50 kpc (Figs.3-4).The isothermal core in the S zone is characterised by densities  H ∼ 0.01 cm −3 , temperatures  ∼ 3 × 10 4 K, and sizes ∼ 0.25  shock ∼ 10 kpc, consistent with predictions for cold streams that feed massive high- galaxies from cosmic-web filaments (Mandelker et al. 2018(Mandelker et al. , 2020b)).The thermal pressure in the filament is roughly constant within the S and V zones, suggesting that the confining pressure of cold streams in the IGM is ∼ (10 − 20) times smaller than their confining pressure once they penetrate the CGM of massive haloes.Outside  shock , the density becomes dominated by the underlying sheet, which penetrates and feeds the filament much the same way that filaments penetrate and feed haloes.
(iii) Dark matter structure: The gaseous filaments are embedded in dark matter filaments that set the gravitational potential wells.These are also well-fitted by models of self-gravitating isothermal cylinders, and dark matter is much more isothermal than gas with its effective temperature (i.e., kinetic energy or velocity dispersion squared) varying by a factor of < ∼ 2 over the range (0.1 − 1.0) shock (Fig. 6).The effective temperature of dark matter is very similar to the post-shock temperature of the gas,  ∼ 8 × 10 5 K, yielding a core radius that is ∼ 4 times larger than that of the gas (the S zone).However, this result is a function of the filament line-mass (Section 7.1).The baryon fraction near  shock approaches the universal value, while in the S zone  b > ∼ 0.5 (Fig. 5).(iv) Line-mass: The line-mass (mass per unit-length) of the dark matter filament is < 10% of the maximal line-mass for hydrostatic equilibrium in a self-gravitating isothermal cylinder at the temperature of the dark matter (Fig. 7).The line-mass for the gas is > 30% of this threshold, but still below it, which implies that intergalactic filaments are typically stable against gravitational fragmentation 17 .
(v) Thermal equilibrium: The cooling time near  shock is only slightly longer than the free fall time,  cool ∼ 3 ff , implying that the post-shock state is unstable to thermal condensation and the formation of multiphase gas.This prevents the gas in the T zone from maintaining hydrostatic equilibrium and results in a cooling flow towards the cold-stream in the S zone, where  cool <  ff (Fig. 8).
(vi) Virial equilibrium: The filament gas is in virial equilibrium per-unit-length, once we account for surface-pressure terms and remove the gas that is rapidly inflowing toward the filament along the sheet, with |  | >  vir = √ Λ (Figs.9-10).The virial radius of the filament is  vir ∼  shock .The associated virial temperature is approximately equal to the post-shock temperature, which is also roughly the dark matter effective temperature, suggesting that the dark matter is in virial equilibrium per unit-length as well.Throughout the T zone, the gas maintains the local virial temperature, with non-thermal motions offset by surface pressure.However, in the V and S zones these motions increase faster than the surface pressure with respect to the gravitational potential energy, resulting in net radial outflows induced by the vortices and the motion of the sheet.On the other hand, the filament is unbound along its axis, at all radii.
(vii) Filament dynamics: The component of vorticity along the filament axis has a characteristic quadrupolar structure, with two sets of counter-rotating vortices centred near the boundary between the T and V zones .However, despite a slight preference for the filament vorticity to align with its axis, there is substantial vorticity along the other axes as well (Fig. 14, right).In the T zone, the filament dynamics is dominated by accretion onto the filament (Fig. 13), with most of the mass flux flowing along the sheet (Fig. 9), whose impact parameter with respect to the filament centre seems related to the formation of the aforementioned vortices.
(viii) Dynamical equilibrium: We quantify the various force terms acting on the filament gas by decomposing the Euler equation into gravitational forces, thermal pressure forces and inertial forces.The latter is composed of ram, turbulent pressure, and shear forces.We find that over a large radial range, these forces balance each other, and the velocity field within the filament is in a steadystate (Fig. 11).However, near the location of the vortices in the outer V zone, and near the location of the shock in the outer T zone, the radial velocity field accelerates outward.Similarly, in the S zone, the velocity field accelerates inwards, damping transient outflows caused by the impact parameter of the sheet.In the T zone, the outwards thermal pressure forces roughly balance the inwardsgravity and the ram pressure forces (Fig. 12).In the V zone, where the dynamics is dominated by quadrupolar vortices, thermal pressure forces are comparable to centrifugal forces due to net average rotation, and these combine to provide support against gravity.At the same time, centrifugal forces due to residual rotation associated with the vortices are cancelledcanceled by shear forces and ram/turbulent pressure (Figs. 12,16).In the S zone, a drop in both centrifugal and thermal pressure forces, due to the development of a strong cooling flow, results in a net inwards acceleration of the velocity field.
While our analysis has been comprehensive and these results are enlightening, we stress that the properties of the filaments studied in 17 Though they may still become unstable in the inner CGM (Mandelker et al. 2018;Aung et al. 2019).
this work certainly depend on the mass of the filament, the redshift, and the environment (see, e.g., Ramsøy et al. 2021).Future studies using similarly high-resolution simulations and employing similar analysis methods to study a much larger sample of filaments across different redshifts and environments will be necessary to elucidate the properties of these extremely important elements of the cosmicweb and their implications for galaxy formation.

Figure B1.
Radial profile of the plasma  parameter, the ratio of thermal to magnetic pressure, stacked among our filament slices.We compute the value of  for each individual gas cell, and then average them over each radial bin weighted by volume (blue line, cyan shaded region), mass (orange line, yellow shaded region), and density (green line, pink shaded region).The three solid lines show the results after stacking the ten slices in log space, while the shaded regions show the 1 −  scatter. evolves from ∼ (10 5 − 10 6 ) at  ∼ 0.1  shock to ∼ (10 10 − 10 11 ) at  ∼  shock , implying that magnetic fields are not dynamically important in the filaments.This is reminiscent of speculations concerning how cold streams penetrate the virial shocks around the CGM in massive halos as seen in simulations (Dekel & Birnboim 2006;Dekel et al. 2009a;Bennett & Sĳacki 2020).However, a comprehensive discussion of this topic is beyond the scope of this paper.

APPENDIX B: MAGNETIC FIELDS AND HUBBLE DRAG
In Section 6.1, we analysed the gas dynamics assuming that this obeyed the Euler equation (eq.24).However, this was just an approximation, as it neglects the effect of magnetic fields and the expansion of the Universe.In this section, we justify both of these approximations.
The gas dynamics in the simulation follow the ideal magnetohydrodynamics (MHD) equations.In particular, the Eulerian momentum equation reads where  is the fluid velocity,  is its density,  its thermal pressure, Φ the gravitational potential, and  the magnetic field.Here,  represents the physical velocity, including the Hubble flow in an expanding Universe.In practice, the simulation stores the peculiar velocity,  pec =  −  ()  with  () the Hubble parameter at redshift , for which an additional term of (− ()) representing the so-called "Hubble drag" must be added to the right-hand side of eq.(B1).However, over the scales of interest,  < ∼  shock ∼ 50 kpc, the Hubble flow  () < ∼ 10 km s −1 is negligible compared to the peculiar velocity.Likewise, the Hubble drag is negligible compared to the other forces discussed in Section 6, at the level of a few percent at most, and has no effect on our results.We therefore have neglected this subtlety, and treat peculiar and physical velocities interchangeably.
The dynamical importance of the magnetic field is often quantified by the plasma  parameter, which is the ratio of thermal to magnetic pressure, where  B =  2 /(8) is the magnetic pressure.In Fig. B1, we show profiles of  stacked among our filament slices.We compute the value of  for each cell in the simulation and compute the mean value at each radius  weighted by volume (blue line, cyan shading), mass (orange line, yellow shading), and density (green line, pink shading).In all cases,  increases monotonically from the S zone to the T zone, consistent with magnetic field amplification due to gas condensation as it cools together with a turbulent dynamo effect.The volume-weighting, which highlights the hot and diffuse volumefilling gas, presents the highest values of , i.e. the lowest values of the magnetic fields, with  ∼ 10 7 in the S zone and  ∼ 10 11 in the T zone.The density-weighting, which highlights cold and dense gas, presents the strongest magnetic fields, but even here  ∼ 10 5 in the S zone and  ∼ 10 10 in the T zone.The magnetic pressure is thus negligible compared to the thermal pressure.
Furthermore, magnetic pressure gradients are negligible compared to thermal pressure gradients.Even in the extreme case where the thermal pressure varies over a characteristic length-scale of order  shock < ∼ 50 kpc, and the magnetic pressure varies over a characteristic length-scale of order the smallest cell size, Δ > ∼ 100 pc, magnetic pressure gradients would still be at least 100 times smaller than thermal pressure gradients.We conclude that magnetic fields are not dynamically important in our filament sample, and do not contribute to filament support against self-gravity.We, therefore, neglect them in our further analysis, simplifying eq.(B1) to the Euler equation (eq.24).
Observations of low- filaments in the vicinity of galaxy clusters do reveal the presence of magnetic fields (e.g., O'Sullivan et al. 2019;Vernstrom et al. 2021), though to our knowledge no such observational constraints exist for high- filaments.However, we note that values of  ∼ 10 5 , as found for cold and dense gas in the S zone, have been assumed for cold streams feeding massive high-z galaxies in idealized studies of stream evolution in the CGM (Ledos et al. 2023).A study of the growth of magnetic fields in intergalactic filaments over cosmic time is beyond the scope of this paper and is left for future work.

APPENDIX C: NUMERICAL DETAILS OF THE FORCE DECOMPOSITION
In this section, we provide further details regarding our numerical method for analysing the forces (Section 6.1-Section 6.2) and present several validation tests.Using our simulations, we numerically evaluate eqs.( 34)-( 49) as a function of radius, , for each filament slice.We begin by generating a uniform cylindrical grid in each slice, with   radial bins in the range  = (0.05 − 1.1) shock ,   azimuthal bins in the range  = (0 − 2), and   axial bins from  = −15 kpc to  = 15 kpc.We assign to each bin the density, velocity, and pressure of the gas cell whose centre is closest to the bin centre.While this does not conserve the total mass, momentum, or thermal energy in the slice, it accurately represents the fluid properties on the grid points, provided the bins are not much larger than the gas cells in the simulation.On the other hand, if the bins are much smaller than the gas cells, this can result in artificially null-gradients when using the grid to compute the terms in the integrals of eqs.( 34)-( 49).We, therefore, wish to use bin sizes as close as possible to the typical sizes of gas cells in the simulation.
In Fig. C1, the solid lines show the radial profiles of the average cell size, defined as the cubic root of the cell volume,  = Vol 1/3 , as a function of / shock stacked among our 10 filament slices.We show profiles for the volume-weighted average cell size (orange), the mass-weighted average cell size (blue), and the density-weighted average cell size (green).Since cells have the same mass to within a factor of ∼ 2, these averages roughly represent the largest, median, and smallest cells at a given radius, respectively.As in previous figures, the cyan-shaded regions represent the 1 −  standard deviations among the slices.The median cell size increases from ∼ 0.5 kpc in  ∼ 0.1 shock , to ∼ 1.5 kpc at  > ∼  shock , while the smallest densest cells grow from ∼ (0.3 − 0.7) kpc in the same range.We also show as dot-dashed lines in Fig. C1 the bin sizes, defined as the cubic root of the bin volume, for cylindrical grids with different numbers of radial, azimuthal, and axial bins as indicated in the legend.Note that in each case, the bin size scales as  1/3 , since the bin volume is ∼  ΔΔΔ.For our fiducial grid we use (  ,   ,   ) = (195, 180, 101), the highest-resolution grid shown in Fig. C1, corresponding to Δ ∼ 0.005  shock ∼ 0.23 kpc, Δ ∼ /90, and Δ ∼ 0.3 kpc.These dimensions are very similar to the density-weighted average cell size at all radii.However, all of our results are very similar when varying the number of bins as in Fig. C1

(see Figs. C2-C3 below)
To verify the validity of our method of depositing gas properties in the grid, which is not strictly conservative, we test the degree to which various gas quantities are conserved.In Fig. C2, we show the ratio of the enclosed gas mass profile inferred from our grid to that of the actual simulation cells.We see that for all but our lowest resolution grid, with (  ,   ,   ) = (87, 70, 43), the error on the enclosed mass is < ∼ 5% at all radii and approaches 0 at  shock .For our highest resolution grid, (  ,   ,   ) = (190,180,101), the error on both the total mass interior to  and the mass contained between  and  + Δ (not shown) is < ∼ 2% at all radii.The ratio of enclosed momenta and thermal energies is similar.We conclude that the fluid properties are sufficiently conserved when generating our highest-resolution grid.
We also experimented with other methods for depositing the fluid properties into the cylindrical grids.Specifically, we deposited the mass, momentum, and thermal energy of each gas cell onto the grid using either a top-hat or a cubic-spline smoothing kernel, with the kernel size set to the cell size (not shown).While this conserves the total mass, momentum, and thermal energy of gas within each slice, the unstructured nature of the mesh and the non-uniform sizes of gas cells result in many "empty" bins, into which no gas cell is deposited, unless we use a kernel size much larger than the cell size.For this reason, we prefer the deposition method described above, which guarantees that each bin has a well-defined non-zero value.However, we note that all three methods for depositing the gas data onto the cylindrical grid yield qualitatively similar results.
We compute the gradients, /, /, /, on the uni- We also experimented with second-order centred finite-difference derivatives: and obtained very similar results with no qualitative differences (see Figure C3 below).Note that in order to compute these derivatives at the grid boundaries self-consistently, we padded our grid with extra bins in both the radial and axial directions, with fluid properties deposited in the same way as described above.This is especially important for the  direction, since the axial terms (eqs. 35, 37, 39, and 49) are evaluated only at the boundaries of  = ±15 kpc.However, all of our results are presented only for the main grid, without the extra padding cells.
Prior to computing each spatial derivative, we smooth the grid with a one-dimensional Gaussian kernel along the direction of the derivative.For example, prior to computing /, we smooth the pressure values on the grid using a smoothing kernel ∝ exp − 2 /(2 2  ) , while prior to computing /, we smooth the pressure values on the grid using a smoothing kernel ∝ exp − 2 /(2 2  ) .This is similar to the smoothing method employed by Lau et al. (2013), though they used one-dimensional Savitzky-Golay filters rather than Gaussians.Our fiducial results use   =   =   = 0.5 kpc, ∼ (1 − 2) times the bin size in our highest resolution grid.We also present results for .Mass conservation test of our method for depositing gas cells into uniform cylindrical grids.We show the ratio of the enclosed gas mass inferred from our grid to the actual enclosed gas mass from the simulation.At all but our lowest resolution (blue line), the cumulative error in the gas mass in the grid is < 5% everywhere and approaches 0 at  shock .At our highest resolution (red line), both the cumulative and the local error in gas mass (not shown) are < 2% at all radii.
no-smoothing and for smoothing with   =   =   = 1.0 kpc in Fig. C3 below, and find overall good agreement between these different methods.In the same figure, we also present results for   = Δ, with  = , ,  and  = 1, 2, which yield very similar results as well.
Using the above procedures, we obtain all the terms in the integrands of eqs.( 36)-(49) (we address the temporal derivatives of eqs.( 34)-(35) below).We then perform the integrations using the "trapezoid-method", the standard second-order centred approximation for numerical integration of evenly-spaced data, We also experimented with a first-order, "rectangular", approximation, which yielded a slightly poorer fit to Gauss's theorem.We did not experiment with higher-order approximations involving higher-order interpolations between the grid points, as we found the trapezoidal method to be satisfactory and numerical integration is inherently less sensitive to noise than numerical differentiation.
We evaluated the robustness and validity of our methodology by directly computing ∇Φ on our grid and testing Gauss's law.The gravitational potential at the location of each gas cell in the simulation is provided in the simulation output and deposited into our uniform cylindrical grids in the same manner as the gas density, pressure, and velocity.We then compute  grav () following eq.( 25), and compare this to  tot (< ), the total mass of gas, stars, and dark matter enclosed within the cylinder of radius .The results of this test are shown in Fig. C3.The -axes show the ratio of  grav ()/ tot (< ), which should be unity everywhere according to Gauss's law.The -axes show the radius normalized to  shock .We computed the radial profile of this ratio for each slice and then stacked the slices by computing the mean and standard deviation of the ratio at each radius .In the left panel, we show results for different grid resolutions, focusing on those where the error on the enclosed mass is < ∼ 5% (Fig. C2).In each case, we used fourth-order derivatives and applied no smoothing prior to computing the derivatives.The three solid curves represent the mean ratio among our ten slices, while the cyanshaded region represents the 1 −  standard deviation of the highest resolution grid.The three different grid resolutions yield extremely similar results, with no systematic differences.Overall, our method for evaluating the enclosed mass using Gauss's law seems to systematically underpredict the true mass, with an average error of ∼ 5% at  > 0.3 shock , and < ∼ 10% at smaller radii.The scatter, however, is more symmetric about an error of 0 and is still limited to < ∼ 10% in most of the radial range.The scatter in the lower resolution versions (not shown) is somewhat larger, but qualitatively very similar with errors in the range of ∼ (10 − 20)%.In the right panel, we focus on the highest resolution grid and examine the effect of different smoothing widths (no smoothing, constant  = 0.5 or 1.0 kpc, and coordinate-dependent   = Δ, where Δ  is the bin width along the  direction with  = , ,  and  = 1, 2) and differentiation methods (second-order and fourth-order), as noted in the legend.At this resolution, the order of the numerical derivatives has hardly any impact on the results for each smoothing method.The version smoothed with  = 0.5 kpc seems to yield the best results, with effectively no error on the enclosed mass at 0.3 < ∼ / shock < ∼ 0.8, and errors of < 5% elsewhere.The version smoothed with  = 1.0 kpc yields the largest errors, of (5 − 10)%.Also, interestingly, unlike the no-smoothing and  = 0.5 kpc versions, which both systematically underestimate the enclosed mass, the  = 1.0 kpc version systematically overestimates the enclosed mass.As expected, given the bin sizes at this grid resolution,   = Δ  is very similar to no smoothing while   = 2Δ  is very similar to 0.5 kpc smoothing.We note that similar errors of ∼ 5% on estimates of the enclosed mass using Gauss's law were also reported in Lau et al. (2013) (the bottom panel of Fig. 1), which formed the inspiration for our methodology, but whose simulations were run on a much simpler Cartesian mesh.
We adopt a constant  = 0.5 kpc as our fiducial smoothing method, along with fourth-order centred finite differences and our highest resolution grid.The results presented in Figs.C2 and C3 suggest that this method produces reliable and robust results, both in terms of the conservation of fluid properties and in terms of accurately evaluating the necessary surface integrals in eqs.( 34)-( 49).
Finally, we address the computation of the two acceleration terms in eqs.( 34)-( 35), which represent a temporal change in velocity within a fixed volume.We attempted to compute these directly using multiple simulation snapshots, experimenting with first, second, and fourth-order finite-difference approximations for the time derivatives.When doing so, we attempted in each snapshot to centre ourselves on the same ten filament slices as in our fiducial snapshot.However, since the filaments themselves are all moving with respect to the simulation box, since the first stage of our filament selection process (Section 2.2) is done by-eye.Since there is some inherent uncertainty associated with defining the filament centre and rest-frame (as discussed in Section 2.2), this process is prone to errors and uncertainties.Furthermore, the time in between consecutive snapshots of the simulation, Δ ∼ 55 Myr near  ∼ 4, is within a factor of < ∼ 2 of both the cooling time and the eddy-crossing time for gas within the V and S zones,  < ∼ 0.6 shock .All of this makes the direct computation of the acceleration terms highly uncertain.As a result, in none of our experiments did the acceleration terms balance the other Testing our numerical methods for evaluating gradients and surface integrals on our cylindrical grids, using Gauss's law.We show the ratio of the gravitational mass at radius ,  grav from eq. ( 25), and the total mass interior to radius .Gauss's law states that this ratio must be unity everywhere.On the left, we show the results using fourth-order differentiation (eq.C1) with no smoothing applied.Different colour solid lines represent the mean ratio among the ten filament slices for different grid resolutions as shown in the legend.The cyan-shaded region represents the 1 −  standard deviation among the filament slices for the highest resolution.All resolutions yield similar results, with  grav systematically underpredicting the true enclosed mass by ∼ 5%.On the right, we show results for our fiducial (highest) resolution grid, but using different methods for smoothing the data and computing the gradients, as indicated in the legend.Using second or fourth-order derivatives has no effect on the results.The best results are obtained when smoothing the data with a Gaussian of  = 0.5 kpc prior to computing the derivatives.In this case, there is effectively no error on the enclosed mass at 0.3 < ∼ / shock < ∼ 0.8, and < ∼ 5% error elsewhere.The worst results are obtained with 1.0 kpc smoothing, but even here the errors are only ∼ (5 − 10)%.This confirms the validity and robustness of our methodology.terms in eq. ( 50).Since we have shown that our method reliably reproduces  grav () (Fig. C3), and since  therm () and  inertial () are computed similarly, we can assume that these are also computed reliably.We, therefore, use eq.( 50) and the corresponding axial equation to infer the acceleration terms from the difference between the gravitational, thermal, and inertial terms.
After computing the radial profiles of each mass-term from eqs. ( 36)-( 49) for each filament slice, we normalise these by  tot (< ) and stack the ten filament slices by taking the average of each normalised mass term in each radial bin.While it is the ratio of each mass term to  grav () that tells us the contribution of the respective force term to the support of the filament against gravity at that radius, we nonetheless opt to normalise the mass terms by  tot (< ).As shown in Fig. C3,  tot (< ) ≃  grav () to within ∼ 5%.However, while  grav () can fluctuate, especially for individual slices,  tot (< ) is a smooth, monotonically increasing function.We have verified that normalising by  grav () instead does not have a noticeable impact on our results.In Section 6.2, we use the notations  ... / tot ,  ... / grav , and  ... / grav interchangeably, and think of these as ratios of forces/accelerations rather than ratios of masses.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Hydrogen column density,  H , of the large-scale cosmic sheet in our simulation at  ∼ 4. Coordinates are given in proper physical distances.The left panel shows the "face-on" view of the sheet, while the right panel shows an "edge-on" view.The integration depth for both panels is ±60 kpc.Haloes with  v ≥ 10 10 M ⊙ within ±60 kpc of the sheet midplane are marked in both panels by black circles, whose radii correspond to the halo virial radii.The three largest haloes in the bottom-left, bottom-right, and top-right of the left-hand panel all have virial masses  v ∼ 10 12 M ⊙ .White boxes mark the ten filament slices selected for analysis, numbered from 1 to 10, selected not to contain any haloes with  v ≥ 10 10 M ⊙ (see Section 2.2).

Figure 2 .
Figure2.A schematic diagram showing the structure of the filament cross section, marking the key physical characteristics of the three zones described in Section 2.3 and detailed throughout the paper.The stream is surrounded by a cylindrical accretion shock, at a normalised radius of R = 1, which expands outwards as the filament grows.The shock-heated gas in the outer Thermal (T) zone, at R > ∼ 0.65, is at the virial temperature of the underlying dark matter filament, with low density and inwardsradial velocity driven by an isobaric cooling flow.In this region, the thermal pressure gradients balance both gravity and the ram pressure from the inflowing gas.In the intermediate Vortex (V) zone, at 0.65 > ∼ R > ∼ 0.25, the dynamics is dominated by a quadrupolar vortex structure which is induced by shearing motions between the sheet that feeds the filament and the post-shock filament gas.These vortices induce a complex combination of inflows, outflows, and rotation in this zone.Gravity here is balanced by a combination of thermal pressure gradients, centrifugal forces, and non-thermal random motions.The inner stream (S) zone, at R < ∼ 0.25, is a dense, cold, isothermal core which represents the cold streams that feed massive galaxies.While the stream boundary expands as the filament grows, the gas in this region is also characterised by both ithats and outflows due to non-radial inflow along the sheet.

Figure 4 .
Figure 4. Distribution of central densities,  0 on the -axis, central temperatures,  0 on the -axis, and shock radii,  shock shown in colour, for our ten filament slices.The central density is noramlised by the average baryonic density at  = 3.93,  b = Ω b  crit (1 + ) 3 = 1.00 × 10 −28 g cm −3 (labelled at the top).The corresponding slice numbers are denoted.Note that slices 03 and 04 in the lower right corner overlap with each other.

Figure 5 .
Figure 5. Baryon fraction as a function of radius.As in Fig. 3, red dots indicate the mean and the cyan-shaded region indicates the standard deviation among our ten slices, while the vertical lines/shaded regions mark the boundaries between the four filament zones.In the PS Zone,  >∼  shock ,  b converges to the universal baryon fraction of ∼ 0.16, marked by the horizontal dashed green line.The gas is more centrally concentrated than dark matter due to efficient cooling, so  b rises through the T and V zones, saturating at  b ∼ 0.5 in the S zone, where the gas and dark matter masses are comparable.

Figure 6 .
Figure6.Left panel: Projected map of the dark matter density averaged along the filament axis and stacked among our ten slices.The stacking procedure and slice orientation are the same as in the projected maps of gas properties shown in Fig.3.The shock radius,  shock , identified based on the gas temperature profiles, is marked by a white circle.As in Fig.3, we also mark the boundaries between the four radial zones with concentric circles.Middle and right panels: Structure of the dark matter filaments, stacked among our ten slices.We show the radial profiles of dark matter density (middle) and effective temperature (right, computed using the dark matter radial velocity dispersion following eq.5).These were normalised, stacked, and averaged in the same way as the gas profiles shown in Fig.3, with red points (cyan shaded region) representing the mean (standard deviation) in log space.The -axes in these two panels span the same range as the upper and middle right-hand panels of Fig.3.The black dashed line in the density profile shows a fit to the same parametric isothermal model used for the gas (eq.1), which matches the data well.The yellow and green stars mark the scale radii for the gas and dark matter density profiles, respectively.The latter is slightly larger than the shock radius,  shock , marked with a vertical blue line.The vertical shaded regions mark the boundaries between the S and V zones and the V and T zones.The effective dark matter temperature increases by a factor of ∼ 2 from the filament core to a peak near  0,DM , and is thus more isothermal than the gas, whose temperature increases by a factor of ∼ 20 over the same radial range.

Figure 8 .
Figure 8. Stacked profile of ratio between cooling and free fall times for gas with  >  0 = 3.6 × 10 4 K, the temperature of the cold stream.Near  shock ,  cool <∼ 3 ff implying that the shocked medium is likely unstable to condensation and precipitation.In the S Zone,  cool < ∼  ff consistent with the formation of a cold isothermal core at these radii.

Figure 9 .
Figure9.Profiles of gas radial velocities (left) and gas mass flow rate (right).Negative values refer to inflow towards the filament axis while positive values refer to outflowing gas.We distinguish between gas outside the sheet (blue lines) and within the sheet (red lines), where the sheet region is crudely defined as |  | < 10 kpc ∼ 0.22  shock (see Fig.3).The dashed-orange line in each panel shows the results for all gases regardless of their position relative to the sheet.The value of   in each bin is the mass-weighted average radial velocity, and the value of  in each bin is   /Δ, where Δ is the width of the bin.Prior to stacking, the   profiles for each slice are normalised by the corresponding virial velocity at  shock ,  vir ( shock ) = (Λ shock ) 1/2 (eq.(23)), while the profiles of  were normalised by Λ (  shock )  vir ( shock ).In the left-hand panel, the green lines show the profiles of ± vir ( )/ vir ( shock ).The sheet is inflowing at all radii, with |  | >  vir ( ) outside the S zone and with |  | <  vir ( ) inside the S zone.On the other hand, the off-sheet gas is inflowing in the T zone and outflowing in the V zone, but always with |  | <  vir ( ).

Figure 10 .
Figure10.Radial profile of the virial parameter,  vir (black line), along with the ratios of the various energy terms in the numerator of eq.(15) to the gravitational energy, | W ( ) |. Specifically, these are the volumetric kinetic and thermal terms, T dyn (green solid line) and T therm (solid magenta line), and the surface kinetic and thermal terms, T dyn, (green dashed line) and T therm, (solid magenta line).We see that  vir , T therm /| W |, and T therm, /| W | are all of order unity in the T and outer V zones, where the gas temperature is hot and comparable to the effective temperature of the dark matter.Therefore, the filaments are in approximate virial equilibrium per unit-length with  vir ∼  shock .

Figure
Figure11.Force equilibrium in the filaments.We compare the sum of thermal and inertial forces (solid lines) with the gravitational force (dashed lines), in both the radial and axial directions, as a function of / shock .The magenta dashed line and cyan shaded region show  grav ( )/ tot (<  ), which is ≃ 1 everywhere with a < ∼ 5% error at / shock < ∼ 0.3 and > ∼ 0.8 due to our numerical method.The navy and pink dashed lines show the ratios  grav,r ( )/ tot (<  ) and  grav,a ( )/ tot (<  ), respectively.The former is positive, indicating a radial gravitational field directed inwards, while the latter is negative, indicating an axial gravitational field directed outwards since the filament is unbound along its axis.Solid navy and pink lines show (  therm,r + inertial,r )/ tot (<  ) and (  therm,a + inertial,a )/ tot (<  ), respectively (eqs.36-39).Radii, where the navy or the pink solid lines are equal to the corresponding dashed lines, indicate regions in a dynamical steady state without temporal acceleration in the corresponding direction (eqs.(34) and/or (35) are equal to 0).If, however, the solid lines are greater (smaller) than the corresponding dashed lines, this indicates a negative (positive) temporal acceleration at this radius.Over a fairly large radial range in the V and T zones, / shock ∼ (0.25 − 0.40) and (0.60 − 0.80), there is a dynamical steady-state in the radial direction.  / > 0 at / shock ∼ (0.40 − 0.60) and (0.80 − 1.0), while   / < 0 in the S zone at / shock < ∼ 0.25.In the axial direction, there is no dynamical steady state at / shock > ∼ 0.2.

Figure 12 .
Figure 12.Detailed force balance in filaments.Left panel:We compare the radial thermal and inertial forces (solid red and green lines) with the radial gravitational force (navy dashed line, as in Fig.11), and the axial thermal and inertial forces (solid orange and yellow lines) with the axial gravitational force (pink dashed line).Centre panel: We show the four constituent forces comprising the radial inertial force (eq.38), namely the centrifugal force (purple lines), the ram/turbulent pressure forces (red lines) and the shear forces (orange and blue lines for azimuthal and axial shear, respectively).For each force, we further distinguish between their mean (dot-dashed lines) and residual (solid lines) components.The legend for this panel is in the box located below all of the panels.Right panel: We show the three constituent forces that comprise the axial inertial force (eq.39).These are the ram/turbulent pressure forces (blue lines) and the two shear forces (red and yellow lines for radial and azimuthal shear, respectively).Based on all three panels, we conclude: In the T Zone, the radial thermal force is directed outwards and is stronger than the gravitational force inwards, while the inertial forces are directed inwards and are dominated by the ram pressure, which mainly offsets the excess thermal pressure.In the V Zone, the radial thermal force is directed outwards and decreases from ∼ 1.5| grav,r | to ∼ 0.6| grav,r | as  decreases.Radial inertial forces are directed outwards with a roughly constant amplitude of ∼ 0.6| grav,r |, and are dominated by centrifugal forces.The mean rotation has a roughly constant amplitude of ∼ | grav,r |.The residual rotation is much larger, but is largely offset by ram/turbulent pressure and azimuthal shear.In the S Zone, the radial thermal and inertial forces become negative at  > ∼ 0.15  shock and  < ∼ 0.15  shock , respectively, resulting in a net inwards acceleration in this region.The axial gravitational force is directed outwards, representing the tidal stretching of the filament along its axis.It is roughly balanced by axial ram pressure forces at all radii, while the axial thermal pressure force is relatively small, directed outwards in the S zone and inwards the T and V zones.

Figure 13 .
Figure 13.The projected velocity field of slice 2, a representative example among our 10 slices.The colour map encodes the average gas density projected along the line-of-sight, and arrows represent the projected 2D velocity field perpendicular to the line-of-sight.The longest arrows correspond to 290 km s −1 , with lengths scaling as  = 30 sinh −1 /30 km s −1 , such that a 30 km s −1 arrow is ∼ 3.36 times shorter than the longest one.Four major vortices forming a quadrupolar pattern are present near the outer boundary of the V zone ( ∼ 0.6  shock , orange circle).These are highlighted with schematic eddies to guide the eye.
vir is the circular velocity of the cylin-drical filament.Despite slice-to-slice variation in the detailed velocity field and in the location and strength of the four vortices, the stacked projection nicely recovers the quadrupolar structure, with the vorticity being primarily positive (negative) in the upper-right and lower-left (upper-left and lower-right) quadrants 15 .

Figure 14 .
Figure 14.Vorticity along the filament axis.Left: Projection map of the vorticity along the filament axis, mass-weighted along the line of sight, normalized by  circ / shock , and stacked among the 10 filament slices.The same quadrupolar structure visible in the single slice shown in Fig. 13 is evident in the stacked map.The upper-right and lower-left (first and third) quadrants have positive vorticity, while the upper-left and lower-right (second and fourth) quadrants have negative vorticity.Right: The stacked profile of the contribution to the filament vorticity oriented along its axis.We calculate the mass-weighted average profiles of  2

Figure 15 .
Figure15.Our ideal vortex model is used to interpret the main features of the radial inertial forces from Fig.12.The red dashed circle represents the "shock" radius, , and we have added circles at 0.6 and 0.25, schematically representing the borders between the three zones in our simulated filaments.Each quadrant has a rotational vortex centred at ( , ) = (±, ±) ≃ (±0.42 , ±0.42 ).These form a quadrupolar structure, with positive (negative) vorticity in the first and third (second and fourth) quadrants.Within each quadrant, the rotation velocity scales as  ′ 0.5 , where  ′ is the distance to the vortex centre.

Figure 16 .
Figure16.Inertial (fictitious) force terms resulting from our toy model of an ideal quadrupolar vortex structure (Fig.15).Since our model is completely symmetric, the mean parts of   and   are, by definition, 0, and therefore only the residual parts contribute.Each force term of our model has been scaled by an additional   to account for the radial dependence of  g , the denominator in the ratios presented in Fig.12.Our simple model reproduces many of the trends seen in the simulation in the V zone, at 0.25 < ∼ / shock < ∼ 0.6 (Fig.12, middle panel, solid lines), where the ram pressure force has decreased and vortices dominate the dynamics of the filament (see the text for details).

Figure 17 .
Figure 17.Correlations between the total line-mass within  shock , Λ shock , and the shock radius,  shock (top) as well as the post-shock temperature,  shock (bottom).Points mark the different slices, which are numbered in each figure (see also Fig. 4).Dashed lines in each panel show our predictions from eqs. (67) and (68), respectively, which are both reasonable matches to the simulation data.

Figure C1 .
Figure C1.Cell-sizes in the simulation compared to bin-sizes of the uniform cylindrical grid we use to evaluate eqs.(34)-(49).Solid lines show the average cell size as a function of radius, weighted by volume (orange), mass (blue), or density (green), and stacked among our 10 filament slices.Dot-dashed lines show the bin sizes of our cylindrical grids, with different numbers of bins as indicated in the legend.At all radii, the bin size in our fiducial grid is comparable to the density-weighted average cell-size, representing the smallest and densest cells in the simulation.

Figure C2
Figure C2.Mass conservation test of our method for depositing gas cells into uniform cylindrical grids.We show the ratio of the enclosed gas mass inferred from our grid to the actual enclosed gas mass from the simulation.At all but our lowest resolution (blue line), the cumulative error in the gas mass in the grid is < 5% everywhere and approaches 0 at  shock .At our highest resolution (red line), both the cumulative and the local error in gas mass (not shown) are < 2% at all radii.

Figure C3.
Figure C3.Testing our numerical methods for evaluating gradients and surface integrals on our cylindrical grids, using Gauss's law.We show the ratio of the gravitational mass at radius ,  grav from eq. (25), and the total mass interior to radius .Gauss's law states that this ratio must be unity everywhere.On the left, we show the results using fourth-order differentiation (eq.C1) with no smoothing applied.Different colour solid lines represent the mean ratio among the ten filament slices for different grid resolutions as shown in the legend.The cyan-shaded region represents the 1 −  standard deviation among the filament slices for the highest resolution.All resolutions yield similar results, with  grav systematically underpredicting the true enclosed mass by ∼ 5%.On the right, we show results for our fiducial (highest) resolution grid, but using different methods for smoothing the data and computing the gradients, as indicated in the legend.Using second or fourth-order derivatives has no effect on the results.The best results are obtained when smoothing the data with a Gaussian of  = 0.5 kpc prior to computing the derivatives.In this case, there is effectively no error on the enclosed mass at 0.3 < ∼ / shock < ∼ 0.8, and < ∼ 5% error elsewhere.The worst results are obtained with 1.0 kpc smoothing, but even here the errors are only ∼ (5 − 10)%.This confirms the validity and robustness of our methodology.