We analyze the internal cross-section structure and dynamics of inter-halo cosmic-web filaments that connect to massive high-$z$ galaxies. Our analysis is based on a high-resolution  cosmological simulation zooming in on a volume encompassing three filaments of $\sim 1 \, \Mpc$ length feeding three massive haloes of $\sim 10^{12}\,\msun$ at $z \sim 4$, embedded in a large-scale sheet. Each filament is surrounded by a strong cylindrical accretion shock of radius $r_{\rm shock} \sim 50 \kpc$, within which the gas is in virial equilibrium in the potential well set by an isothermal dark-matter filament, with a virial temperature of $\sim 7 \times 10^5\,{\rm K}$. The total filament line-mass is $\sim 9\times 10^8\,\msun\,\kpc^{-1}$ and the gas to total mass ratio within $r_{\rm shock}$ is the universal baryon fraction. These values all match the expectations from analytical models for filament properties as a function of halo-mass and redshift.
The cross-section of the filament can be divided into three radial zones. In the outer “thermal” (T) zone, $r \geq 0.65 \, r_{\rm shock}$, the inward gravity and ram-pressure forces from the inflowing gas are over-balanced by outward thermal pressure forces, causing the inflowing gas to decelerate and the shock to expand outwards. In the intermediate “vortex” (V) zone, $0.25 \leq r/  r_{\rm shock} \leq 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 outward force is provided mainly by centrifugal forces associated with these vortices and, to a lesser extent, by global rotation and thermal pressure, while shear and turbulent forces that are also associated with the vortex structure add to the force inward. The inner “stream” (S) zone, $r < 0.25 \, r_{\rm shock}$, is a dense isothermal core defining the cold streams that feed the galaxies, with $T\sim 3 \times 10^4 \, \Kdegree$ and $n_{\rm H}\sim 0.01 \cmc$. The core is formed by an isobaric cooling flow that results from post-shock cooling on a free-fall time, and is associated with a decrease in the outward forces, though the gas within the core exhibits both inflows and outflows.

intergalactic medium – large-scale structure of Universe – hydrodynamics – methods: numerical – methods: analytical

§ INTRODUCTION

In the last two decades, our understanding of galaxy formation in a cosmological context has undergone a paradigm shift. Galaxies, 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 <cit.>, and can be seen in the distribution of galaxies in observational surveys such as 2dFGRS <cit.>, SDSS <cit.>, and the 2MASS redshift survey <cit.>. Perhaps the most prominent feature of the cosmic-web is a network of intergalactic filaments, which contain roughly half of the mass in the Universe <cit.>.

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 slices below). At high redshift, $z\sim (2-6)$ near the peak of cosmic star- and galaxy-formation, these filaments manifest as streams of cold, dense gas ($T\sim 10^4\,\Kdegree$, $n\sim 10^{-2}\,{\rm cm}^{-3}$), predicted to be the main mode of accretion onto high-$z$ galaxies, feeding them directly from the cosmic-web. In massive haloes with virial mass $\Mv \gsim 10^{12}\,\msun$ above the critical mass for the formation of a stable virial shock <cit.>, the streams are predicted to penetrate the hot circumgalactic medium (CGM), maintaining temperatures of $T_{\rm s}\gsim 10^4\Kdegree$ and reaching the central galaxies in roughly a virial crossing time <cit.>. This gas supply sustains the large observed star formation rates (SFRs) of $\sim 100\,\msun\yr^{-1}$ without
galaxy mergers <cit.>.

Though their narrow size and low density make the streams difficult to detect observationally directly, numerous observational studies of the CGM and the intergalactic medium (IGM) around massive high-$z$ galaxies, in both absorption <cit.> and
emission <cit.>, reveal large quantities of cold gas with spatial and kinematic properties consistent with the predictions for cosmic filaments and cold streams. Cold streams are ubiquitous in cosmological simulations <cit.>,
where they are found to be the primary source of gas accretion onto both the halo and the central galaxy <cit.>.

Besides supplying galaxies with cold gas to fuel ongoing star formation, additional effects of filaments and streams on galaxies and DM haloes have been extensively explored theoretically. These include halo spin, galaxy angular momentum growth and disk formation <cit.>, the shape and alignment of galaxies and haloes <cit.>, driving and sustaining turbulence and disk instabilities <cit.>, and the emission and absorption signatures of the CGM <cit.>. Gravitational instabilities in and fragmentation of filaments and streams have also been speculated to contribute the formation of stars and even globular clusters outside of galaxies <cit.>.

Despite their prominence in the large-scale structure of the cosmic-web, 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 works have suggested that this could be due to thermal pressure gradients <cit.>, coherent filament rotation <cit.>, helical vortices <cit.>, turbulence <cit.> 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 which have measured filament properties in large-volume cosmological simulations <cit.> commonly involve either DM-only N-body simulations or low-resolution hydrodynamic simulations which 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 <cit.>, 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 ionization 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-$z$ using cosmological zoom-in simulations with relatively high resolution was made by <cit.>. These authors used a cosmological zoom-in simulation of a Milky Way (MW)-mass halo, $\Mv\sim 5\times 10^{11}\msun$ at $z=0$, taken from the  suite <cit.> and run with the adaptive mesh refinement (AMR) code  <cit.>. They studied the properties of the main filament feeding the central halo in the redshift range $z\sim 3.5-8$, focusing primarily on $z\sim 4$ when the filament extended $\lsim 200\kpc$ on either side of the halo and did not seem to connect to any other massive object. The spatial resolution in the analyzed filament was $\sim 1.2\kpc$, with small patches around haloes embedded in the filament reaching a resolution of $\sim 0.6\kpc$. These authors found the filament to be well described by a model of an isothermal, self-gravitating, infinite cylinder <cit.>, embedded inside an isothermal self-gravitating sheet that dominated the mass distribution at $r\gsim (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 post-shock gas was found to cool rapidly, and the shock width was limited to a few kpc. The azimuthally-averaged radial profile of gas temperature increased by $\lsim 60\%$ between the central value and the peak of the shock. While significant vorticity was identified within the filament, consistent with previous results <cit.>, the filament was found to be primarily supported by thermal pressure.

The analysis and results of <cit.> are both novel and extremely detailed. However, since their work focused on a single filament feeding an isolated high-$z$ progenitor of a MW-mass halo, it is unclear whether their results describe filaments feeding massive, $\Mv\gsim 10^{12}\msun$, haloes as cold streams at high-$z$, 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 in order to draw more general conclusions about the filament population.
Finally, while the resolution of $\sim 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 $\sim 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., ).

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 <cit.> and <cit.>, which zoom in on a large patch of the IGM in between two massive haloes, which at $z\sim 2$ are $\Mv\sim 5\times 10^{12}$ M$_{\odot}$ each and connected by a $\sim$ Mpc scale cosmic filament. We hereafter refer to this suite of simulations as , where “IPM” refers to the “intra-pancake medium” of multiphase gas within cosmological sheets (or “Zel'dovich pancakes,” see ; ). The large filament present at $z\sim 2$ was formed by the merger of several smaller filaments at $z \sim 2.5$. Our analysis focuses on $z\sim 4$, when the system contained several well-defined filaments connecting three haloes with $\Mv\sim 10^{12}\,\msun$ each (see slices 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 Sims), the typical cell size in the IGM at $z\sim 4$ is $\sim 1\kpc$ near the mid-plane of the sheet and $\sim 300\pc$ near the cores of the filaments.
A detailed convergence study presented in <cit.> suggests that this appears sufficient to resolve the formation of multiphase gas in the turbulent IPM, though the morphology and distribution of the cold component are likely only marginally resolved.
As shown in thermal below, the characteristic scale of the accretion shock surrounding our filaments is $r_{\rm shock}\sim 50\kpc$, while the cold stream radius is $r_{\rm stream}\sim 0.25r_{\rm shock}\sim 12\kpc$. We thus have $\sim 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 organized as follows: In Methods, 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 Radial_prop, we analyze the radial structure of gas and dark matter in filaments, including their thermal structure. In tcool, we discuss the thermal stability of gaseous filaments by comparing their cooling and free-fall time scales. In virial, we address whether the gaseous filaments are in virial equilibrium within the potential well set by the dark matter filament. In fil_dyn, we decompose the forces that are supporting 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 disc. Finally, in conclusions, we conclude.

§ METHODS

§.§ Simulation Method

We use the highest resolution version of   introduced in <cit.> and <cit.>. 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  <cit.>. 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 $z\sim 2.3$ snapshot of the Illustris TNG100[http://www.tng-project.org] magnetohydrodynamic cosmological simulation <cit.>. These have virial masses $\Mv$ (defined using the  spherical overdensity) in the range $\sim (1-40)\times 10^{12}\msun$. We then selected from these all halo pairs with a comoving separation in the range $(2.5-4.0)h^{-1}\Mpc \sim (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 $\Mv\sim 5\times 10^{12}\msun$ each, separated by a proper distance of $D\sim 1.2\Mpc$ at $z\sim 2.3$. By $z=0$, the two haloes evolve into mid-size groups with $\Mv\sim (1.6-1.9)\times 10^{13}\msun$ separated by $\sim 2.7\Mpc$, such that their comoving distance has decreased by $\lsim 30\%$.

The zoom-in region is the union of a cylinder with length $D$ and radius $R_{\rm ref}=1.5\times R_{\rm v,max}\sim 240\kpc$ extending between the two halo centers, and two spheres of radius $R_{\rm ref}$ centred on either halo, where $R_{\rm v,max}$ is the larger of the two virial radii at $z=2.3$. We trace all dark matter particles within this volume back to the initial conditions of the simulation at $z=127$, refine the corresponding Lagrangian region to higher resolution, and rerun the simulation to a final redshift, $z_{\rm fin}$. In this paper, we use the highest resolution version of this simulation, dubbed ZF4.0 in <cit.>, which has a dark matter particle mass of $m_{\rm dm}=8.2\times 10^4\,\msun$ and a Plummer-equivalent gravitational softening of $\epsilon_{\rm dm}=250\pc$ comoving. Gas cells are refined such that their mass is within a factor of 2 of the target mass, $m_{\rm gas}=1.5\times 10^4\msun$. Gravitational softening for gas cells is twice the cell size, down to a minimal gravitational softening $\epsilon_{\rm gas}=0.5\epsilon_{\rm dm}=125\pc$. The simulation was run until $z_{\rm fin}\sim 2.9$, though our analysis here focuses on $z\sim 4$.

The simulations were performed with the same physics model as used in the TNG simulations, described in detail in <cit.> and <cit.>. This includes radiative cooling down to $T=10^4{\rm K}$ from Hydrogen and Helium <cit.>, metal-line cooling <cit.>, radiative heating from a spatially constant and redshift-dependent ionizing ultraviolet background (UVB) <cit.>, partial self-shielding of dense gas from the UVB <cit.>, and additional heating from the radiation field of active galactic nuclei (AGN) within $3\Rv$ of haloes containing actively accreting supermassive black holes <cit.>. Star-formation occurs stochastically in gas with densities greater than $n_{\rm thresh}=0.13\cmc$, which is placed on an artificial equation of state meant to mimic the unresolved multiphase ISM <cit.>. 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 $b = 0.2$ to the dark matter particles. We then assign gas and stars to FoF groups based on their nearest-neighbour dark matter particle. Finally, we apply  <cit.> to the total mass distribution in each FoF group. The most massive  object in each FoF group is identified as the central halo, whose virial radius $\Rv$ is defined using the <cit.> spherical overdensity criterion, and virial mass $\Mv$ is the total mass of dark matter, gas, and stars within $\Rv$.

§.§ Filament Selection

Hydrogen column density, $N_{\rm H}$,  of the large-scale cosmic sheet in our simulation at $z\sim 4$. All the coordinates are denoted by proper physical distances. The left panel shows the “face-on” view of the sheet, while the right panel shows an “edge-on” view. Haloes with $\Mv\ge 10^{10}\,\msun$ are marked by black circles whose radius corresponds to the halo virial radius. The three largest haloes in the bottom-left, bottom-right, and top-right of the left-hand panel all have virial masses $\Mv\sim 10^{12}\,\msun$. White boxes mark the ten filament slices selected for analysis, numbered from $1$ to $10$, selected not to contain any haloes with $\Mv\ge 10^{10}\,\msun$ (see Fil_selec).

We focus our analysis on $z\sim 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 $\Mv\sim 10^{12}\,\msun$, and along which nearly all haloes with $\Mv>10^9\,\msun$ are located. The sheet is formed by the collision of two smaller sheets at $z\gsim 5$, while most of the filaments merge by $z\sim 2.5$ <cit.>. slices shows the Hydrogen column density, $N_{\rm H}$, of the sheet at $z\sim 4$ in face-on (left) and edge-on (right) projections. In order to minimize contamination of filament properties by haloes, we avoid analyzing filament regions that contain haloes with $\Mv\ge 10^{10}\,\msun$, marked by black circles in slices.
We selected ten such slices from three separate filaments, marked with white rectangles and numbered from $1-10$ in slices. Each slice is $30\kpc$ thick along the filament axis, slightly larger than the typical characteristic scale radius of the filaments, $r_0$, as described in thermal below. However, we note that using a slice thickness of $(20-40)\kpc$ does not change any of our key results.

The orientation of the filament axis in each slice is determined by eye by approximating a straight line to the ridge-line of maximal $N_{\rm H}$ in each slice, as seen in the face-on projection through the sheet (left-hand panel of slices). 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 $z$-axis to refer to the local filament axis in each slice, while the $(x,y)$ plane represents the filament cross-section with the $y$-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 $R_0=50\kpc$ and length $L=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\sim 1.6\kpc$, while keeping $L$ 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 $\sim (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 Radial_prop-fil_dyn 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 cartoon.

* Outer thermal (T) zone, $0.65 r_{\rm shock}\lsim r \lsim r_{\rm shock}$: The post-shock region, where the gas is hot (at the virial temperature, see virial) and the dynamics are dominated by thermal and ram pressure.

* Intermediate vortex (V) zone, $0.25 r_{\rm shock}\lsim r\lsim 0.65 r_{\rm shock}$: The gas is cooling and the dynamics are dominated by a quadrupolar vortex structure and increasing centrifugal forces (see fil_dyn).

* Inner stream (S) zone, $0\lsim r \lsim 0.25 r_{\rm 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 $r\gsim r_{\rm shock}$. This is the region outside the accretion shock surrounding the filament, where cold and diffuse gas is free-streaming towards the filament.
While $r_{\rm shock}$ is well-defined (thermal), the boundaries between the T and V zones and between the V and S zones are not as sharp. We assign a width of $\sim 0.05r_{\rm shock}$ to each of them, $r/r_{\rm shock}\sim (0.6-0.65)$ and $(0.2-0.25)$ respectively, marked in the figures below.

A schematic diagram showing the structure of the filament cross section, marking the key physical characteristics of the three zones described in zones and detailed throughout the paper. The stream is surrounded by a cylindrical accretion shock, at a normalized radius of ${\rm R}=1$, which expands outwards as the filament grows. The shock-heated gas in the outer Thermal (T) zone, at ${\rm R}\gsim 0.65$, is at the virial temperature of the underlying dark matter filament, with low density and inwards radial velocity driven by an isobaric cooling flow. In this region, thermal pressure gradients balance both gravity and ram pressure from inflowing gas. In the intermediate Vortex (V) zone, at $0.65\gsim {\rm R}\gsim 0.25$, the dynamics are 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 ${\rm R}\lsim 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 characterized by both inflows and outflows due to non-radial inflow along the sheet.

§ FILAMENT RADIAL STRUCTURE

In this and the following sections, we analyze the radial structure of various filament properties. These include the gas density, temperature, and thermal pressure (thermal), the baryon fraction (fbar), the dark matter density and radial velocity dispersion (DM_prop), the mass-per-unit-length of gas and total mass (Line_Mass), the ratio of cooling time to free-fall time (tcool), the virial parameter (virial), and the internal kinematics and dynamics (fil_dyn). We compute all profiles as a function of the cylindrical radius, i.e., the distance from the local filament axis as defined in Fil_selec, with each radial bin representing a cylindrical shell extending the thickness of the slice, $L=30\kpc$.

§.§ Thermal Properties and Characteristic Radii

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). See the text for details regarding the averaging and stacking procedures. The temperature profiles reach a well-defined maximum that defines a shock radius, $r_{\rm shock}$, marked by white circles 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, $r\lsim 0.25 \, r_{\rm shock}$, where the gas is dense and cold. In the V zone, $r/r_{\rm shock}\sim (0.25-0.65)$, the pressure is roughly constant while the temperature increases by a factor of $\sim 3$. In the T zone, $r/r_{\rm shock}\sim (0.65-1.0)$, the pressure drops by a factor of $\sim 2$, the temperature increases by a factor of $\sim 5$, and the density develops a slope of $\rho \propto r^{-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 \, r_{\rm shock}$ in the radial profiles. At $r<r_{\rm shock}$, the density profile is well-fit by that of an infinite self-gravitating isothermal cylinder in hydrostatic equilibrium, shown by the black dashed line, with a scale radius of $r_0\sim 0.5\,r_{\rm 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.

In thermal profiles, 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 $z$ direction) across the slice thickness of $L=30\kpc$. The average density and pressure are weighted by volume, while the temperature is weighted by mass. Prior to stacking, we align the $y$-axes such that they lie in the sheet containing the filaments (see slices), while the $x$-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 $y$-direction breaking the circular symmetry represents the sheet. There are well-defined accretion shocks around the filament and the sheet, visible in both the temperature and pressure maps. While the accretion-shock around the sheet is discussed in <cit.>, those works did not explicitly identify filamentary accretion shocks as seen here. Such filamentary accretion shocks have been both predicted theoretically <cit.> and seen in other cosmological simulations <cit.>. While the sheet clearly breaks the circular symmetry of the filament cross-section, for most of our analysis (except radial inflows 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 $R=100\kpc$ for each individual slice, namely a bin size of $\Delta r = 100\kpc/150\approx0.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 bin sizes and cell sizes below). However, due to the unstructured nature of the grid, all of our bins are sufficiently populated for good statistics.
The density and pressure (temperature) profiles represent the volume-weighted (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, $r_{\rm shock}$. Once $r_{\rm shock}$ is identified, we recompute each profile using 150 linearly spaced bins from $r/r_{\rm shock}=0.1-2.0$ and normalize 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 thermal profiles. Hereafter, all stacked profiles and projection maps presented used bins proportional to $r_{\rm shock}$ of each slice. The location of the shock seen in the temperature and pressure maps is indeed at $r\simeq r_{\rm shock}$, marked by a dashed-white circle in each projection. We show the values of central density, $\rho_0$, central gas temperature, $T_0$, and shock radius, $r_{\rm shock}$, for our ten slices in rho_0 and T_0 scatterings.

At small radii, the filaments are characterized by a core of constant density and temperature extending to roughly $\sim (0.2-0.25)\,r_{\rm shock}$, namely the S zone defined in zones. 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 $\sim (10-20)$ times larger than the core temperature. Within the V zone, $r\sim (0.25-0.65)r_{\rm shock}$, the pressure remains roughly constant, decreasing by only $\sim 10\%$ from its central value. The temperature increases by a factor of $\sim 3$ from its central value, often considered the threshold for `cool' gas in studies of multiphase gas mixing <cit.>, 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, $r\sim (0.65-1.0)r_{\rm shock}$, the pressure declines by a further factor of $\lsim 2$, the temperature rapidly increases by a factor of $\sim 5$, and the density reaches a slope of roughly $\rho\propto r^{-3}$ at $r\lsim r_{\rm 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 <cit.>. 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.

Distribution of central densities, $\rho_0$ on the $x$-axis, central temperatures, $T_0$ on the $y$-axis, and shock radii, $r_{\rm shock}$ shown in colour, for our ten filament slices. The central densities are also normalized by the average baryonic density at $z=3.93$, $\overline{\rho}_{\rm b}=5.00\times 10^{-29}{\rm 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.

Most of the slices have central densities $\rho_0 \sim (10^{-26.5}-10^{-25.5})\cmc$, corresponding to a Hydrogen number density of $n_{\rm H}\sim 0.01\cmc$ and to a baryonic overdensity of $\sim (100-300)$. Their central temperatures are $T_0\sim 2\times 10^4\K$, placing the filament core in approximate thermal equilibrium with the UVB <cit.>. The shock radii are typically $r_{\rm shock}\sim (35-50)\kpc$, with an average value of $\sim 45\kpc$ among the ten slices. The typical core size is thus $r_{\rm core}\sim 0.25\,r_{\rm shock}\sim 12 \kpc$. These values of $\rho_0$, $T_0$, and $r_{\rm core}$ are in good agreement with those predicted for cold streams feeding massive high-$z$ galaxies from the cosmic-web <cit.>. 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 $T_0\gsim 10^5\K$. Both of these slices lie close to the intersection between two filaments (see slices 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 $\lsim 3\times 10^4\K$. Likewise, slices 01 and 02 have unusually extended shocks, with $r_{\rm shock}\sim (70-80)\kpc$. Slice 01 is adjacent to a massive, $\sim 10^{12}\msun$ halo, which likely influences the extent of the shock. 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.

To further characterize the filament radial structure, we fit the stacked density profile shown in the upper-right hand panel of thermal profiles to the density profile of an infinitely long, self-gravitating, isothermal cylinder in hydrostatic equilibrium <cit.>,

\begin{equation}
\label{eq:isothermal_density}
\rho_{\rm isothermal}(r)=\rho_0\left[1+\left(\dfrac{r}{r_0}\right)^2\right]^{-2},
\end{equation}

\begin{equation}
\label{eq:r0}
r_0=\sqrt{\frac{2k_{\rm B}T_0}{\mu m_{\rm p}\pi G \rho_0}},
\end{equation}

where $\rho_0$ is the central density, $G$ is the gravitational constant, $k_{\rm B}$ is the Boltzmann constant, $m_{\rm p}$ is the proton mass, $\mu$ is the mean molecular weight, and $T_0$ is the isothermal temperature. While the filament is clearly not isothermal within $r_{\rm shock}$, several previous works have attempted to model intergalactic filaments as such <cit.>. 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 thermal profiles. As we can see, it actually fits the density profile quite well, except at large radii, $r\gsim r_{\rm shock}$, where the density becomes dominated by the sheet rather than the filament. Similar results were found in <cit.>, who attempted to fit the density profile with a combination of an isothermal cylinder and an isothermal sheet, though we make no such an attempt here.

When fitting our stacked density profile to isothermal_density, the only free parameter is $r_0/r_{\rm shock}$, since the density has been normalized to its central value so effectively $\rho_0=1$[Fitting for $\rho_0$ as well produces a value extremely close to $1$.]. From the fit, we derive a characteristic scale radius of $r_0\sim 0.5\,r_{\rm 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

\begin{equation}
\label{eq:isothermal_temperature}
T_{\rm isothermal}=\frac{r_0^2\mu m_{\rm p}\pi \rho_0 G}{2k_{\rm B}},
\end{equation}

and the isothermal pressure profile is given by

\begin{equation}
\label{eq:isothermal_pressure}
P_{\rm isothermal}(r)=\frac{\pi G\rho_0^2r_0^2}{2}\left[1+\left(\frac{r}{r_0}\right)^2\right]^{-2}.
\end{equation}

Plugging the average value of $\rho_0\sim 10^{-26}\cmc$ and the derived value of $r_0\sim 0.5r_{\rm shock}\sim 22\kpc$ into isothermal_temperature returns a value of $T_0\sim 3.6\times 10^4 \, \K$, which is roughly the average value of $T_{\rm isothermal}$ measured among our ten slices (see rho_0 and T_0 scatterings). The density profile thus resembles that of an isothermal filament at the core temperature. We address this further in comparison. 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 sim_dynamics.

Baryon fraction as a function of radius.
As in thermal profiles, 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, $r\gsim r_{\rm shock}$, $f_{\rm b}$ converges to the universal baryon fraction of $\sim 0.16$, marked by the horizontal dashed green line. The gas is more centrally concentrated than the dark matter due to efficient cooling, so $f_{\rm b}$ rises through the T and V zones, saturating at $f_{\rm b}\sim 0.5$ in the S zone, where the gas and dark matter masses are comparable.

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 thermal profiles. The shock radius, $r_{\rm shock}$, identified based on the gas temperature profiles, is marked by a white circle. As in thermal profiles, 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 DM_temp). These were normalized, stacked, and averaged in the same way as the gas profiles shown in thermal profiles, with red points (cyan shaded region) representing the mean (standard deviation) in log space. The $y$-axes in these two panels span the same range as the upper and middle right-hand panels of thermal profiles. The black dashed line in the density profile shows a fit to the same parametric isothermal model used for the gas (isothermal_density), 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, $r_{\rm 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 $\sim 2$ from the filament core to a peak near $r_{\rm 0, DM}$, and is thus more isothermal than the gas, whose temperature increases by a factor of $\sim 20$ over the same radial range.

§.§ Baryon Fraction

In fbar we show the baryon fraction as a function of radius, stacked among our ten slices as described in thermal. The baryon fraction in each radial bin, $f_{\rm b}(r)$, is given by the ratio of baryonic mass (gas and stars) to total mass (gas, stars, and dark matter) interior to radius $r$. In practice, the stellar mass is negligible, and $f_{\rm b}$ can be thought of as the gas fraction. In the PS zone, $r\gsim r_{\rm shock}$, $f_{\rm b}(r)$ converges to the universal baryon fraction, $f_{\rm b}\sim 0.16$, with very little scatter among the different slices. Within the filament core, the S zone, $f_{\rm b}$ saturates at roughly $0.5$ on average, with a large scatter of $\sim (0.3-0.7)$. The gas is thus more centrally concentrated than the 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 $\sim 0.5$.

§.§ Dark Matter Properties

Following the discussion in thermal and fbar, we now examine the radial structure of the dark matter in our filament slices. In the left panel of DM prop, 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 thermal profiles. Similar to gas, the dark matter filament has a roughly circular cross-section if ignoring the contribution of the overdense pancake along the $y$ direction. However, unlike the gas density which is very centrally concentrated, the dark matter density appears approximately constant within $r_{\rm 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 $\Mv<10^{10}\msun$ (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 all disappears once haloes with $\Mv\gsim 10^8 \msun$ 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 DM prop, we present the radial profile of dark matter density, normalized by its central value and stacked among our ten filament slices, using the same procedure described in thermal for the gas profiles. We fit the dark matter density profile to the same parametric isothermal model used for the gas (isothermal_density). The fit is shown in the middle panel of DM prop with a black-dashed line and is a good representation of the data. From fbar we know that the central dark matter density is comparable to the central gas density, namely an overdensity of $\sim 30$ with respect to the mean matter density[Corresponding to an overdensity of $\sim 200$ with respect to the mean baryonic density, as shown in rho_0 and T_0 scatterings.] at $z\sim 4$. From the fit, we extract the scale radius of the DM density $r_{0,{\rm DM}}\simeq 1.17\, r_{\rm shock}$, which is marked in DM prop along with the scale radius of the gas density profile, $r_0\lsim 0.5\,r_{\rm shock}$. The radius of the DM density core is $r_{\rm core,\,DM}\sim 0.8r_{\rm shock}$, roughly 4 times larger than the gas density core. This is consistent with the baryon fraction profile (fbar) which showed that the gas was much more centrally concentrated than the dark matter. This is unlike the case studied in <cit.>, where the gas and dark matter were found to have similar concentrations and density profiles. We address this in comparison.

In the right-hand panel of DM prop, we present the radial profile of dark matter effective temperature, normalized by its central value and stacked among our ten filament slices as described above. This is derived from the radial velocity dispersion of the dark matter, which we compute in the standard way as $\sigma_r^2 = \langle v_r^2\rangle-\langle v_r\rangle^2$, where $v_r$ is the radial velocity with respect to the filament axis and $\langle \dots\rangle$ represents a mass-weighted average[Note that since in our simulation every dark matter particle has the same mass this is simply the arithmetic variance.] over all dark-matter particles in a given cylindrical shell at radius $r$. The effective dark matter temperature[By assuming a mean particle mass of $\mu m_{\rm p}$, DM_temp defines the temperature the gas would have if it were in equipartition with the dark matter, similar to the way the halo virial temperature is typically defined.] is then defined via

From the profile, we see that $T_{\rm eff}$ increases by a factor of $\sim 2$ from the filament core to a peak just outside $r_{\rm shock}$, near $r_{\rm 0,DM}$. This is much less than the factor $\sim 20$ increase in gas temperature, showing that the isothermal approximation is more valid for dark matter than for gas. The average central value of $T_{\rm eff}$ in the filament core is $T_{\rm eff, C}\sim 9.6\times 10^5 \Kdegree$. Using the average central dark matter density, our fit to $r_{0, {\rm DM}}$, and isothermal_temperature, we calculate the effective “isothermal” temperature of the dark matter, $T_{\rm iso,DM}=1.2\times 10^6\, \Kdegree$. 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, $\rho\propto r^{-2}$, rather than a central core as we find for filaments.

§.§ Mass per unit length (line-mass)

Profiles of the filament line-mass, $\Lambda$, 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 normalized by the corresponding, temperature-dependent, maximal line-mass for hydrostatic equilibrium, $\Lambda_{\rm max}$ (Lmax). While the total line mass is clearly larger than the line-mass of gas only, the larger temperature of the dark matter results in a larger value of $\Lambda_{\rm max}$, leading to a lower normalized line-mass. While the normalized line-mass for gas approaches unity near the filament edge, the total line-mass is well below the critical value.

In line mass we show radial profiles of the mass per-unit-length (hereafter line-mass), $\Lambda(<r)$, normalized and stacked among our ten filament slices. For a given slice, $\Lambda(<r)$ is simply defined as the total mass (of gas and/or dark matter, see below) interior to radius $r$, divided by the slice thickness, $L=30\kpc$.
Motivated by our fits of the gas and dark matter density to isothermal profiles (thermal, DM_prop), prior to stacking we normalize the line mass in each slice by the maximal possible value of a self-gravitating isothermal cylinder in hydrostatic equilibrium <cit.>,

where $c_s^2=k_{\rm B}T/(\mu m_{\rm p})$ is the isothermal sound speed squared, and $T$ is the effective isothermal temperature given by the model fit using the average central density, the characteristic scale radius, $r_0$, and isothermal_temperature.
This is similar to the measured central temperature in the filament core. In line mass, we show separately the normalized line-mass profiles of gas alone (red) and total mass (blue, primarily gas and dark matter). Note that while the total line-mass is obviously larger than the line-mass of gas alone, the gas has a much lower effective
temperature than the dark matter, $\sim 3.6\times 10^4\,{\rm K}$ and $\sim 1.2\times 10^6\,{\rm K}$ respectively, while the effective
temperature derived from the density profile of total mass
is $\sim 8.6\times 10^5\,{\rm K}$. This leads to a much smaller value of $\Lambda_{\rm max}$ for the gas alone compared to the total mass — $2.3\times 10^8\, {\rm M}_\odot \, \kpc^{-1}$
and $5.5\times 10^9\, {\rm M}_\odot \, \kpc^{-1}$ respectively.
As a result, the normalized line-mass of gas alone is larger than the normalized total line-mass.

For the gas, $\Lambda/\Lambda_{\rm max}\sim 0.2$ near the outer boundary of the S zone at $r\sim r_{\rm core}\sim 0.25 \, r_{\rm shock}$, $\sim 0.7$ near $r_{\rm shock}$, and approaches unity at $2\,r_{\rm shock}$. This is consistent with the results of <cit.> (see their equation 40), who predicted that cold streams feeding massive high-$z$ galaxies should have $\Lambda\sim \Lambda_{\rm 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 at small radii (see sim_dynamics below), and the non-isothermality of the gas which is important at large radii (see thermal profiles). 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 $\Lambda/\Lambda_{\rm max}\sim 0.02$ and $0.15$ within $r_{\rm core}$ and $r_{\rm 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, $t_{\rm cool}$, to the free-fall time, $t_{\rm ff}$.

The gas cooling time, $t_{\rm cool}$, is given by

\begin{equation}
\label{eq:tcool}
t_{\rm cool}=\frac{u\rho}{\mathcal{L} n_{\rm H}^2},
\end{equation}

where $u$ is the gas internal energy per unit mass, $\rho$ is the gas density, $n_{\rm H}$ is the hydrogen number density, and $\mathcal{L}=\mathcal{C}-\mathcal{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 $t_{\rm cool}$ for each cell. We then compute the average cooling time as a function of radius, $t_{\rm cool}(r)$, by taking the mass-weighted average of tcool among all cells in each radial bin with temperatures $T>T_0=3.6\times 10^4\K$, the temperature of the cold, isothermal core. Gas with $T\lsim T_0$ is near thermal equilibrium with the UVB and is often undergoing net heating, while virtually all gas with $T>T_0$ is undergoing net cooling.

For an infinite cylinder, the free-fall time is given by

\begin{equation}
\label{eq:tff}
t_{\rm ff}(r)=\sqrt{\frac{1}{4G\rho_{\rm mean}(<r)}},
\end{equation}

where $\rho_{\rm mean}(<r)=\Lambda(<r)/(\pi r^2)$ is the mean density interior to $r$.

For isobaric cooling, $u\rho\propto P \sim {\rm const}$, so $t_{\rm cool}\propto \rho^{-2}$, while $t_{\rm ff}\propto \rho^{-1/2}$. The ratio $t_{\rm cool}/t_{\rm ff}$ is thus smaller at smaller radii closer to the filament axis. When $t_{\rm cool}/t_{\rm ff}\lsim 1$,
the system cannot maintain the shock-heated component and the bulk of the gas cools and falls to the centre <cit.>. The outer radius where $t_{\rm cool}/t_{\rm ff}= 1$ is often referred to as the “cooling radius”. Even if $t_{\rm cool}/t_{\rm 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 galaxy groups or the ICM in galaxy clusters, this process is often referred to as precipitation and occurs when $t_{\rm cool}/t_{\rm ff}\lsim 10$ <cit.>. The analogous threshold for cylindrical geometry may be different due to the different radial dependence of stabilizing buoyancy forces. A detailed study of this is beyond the scope of the current paper, and we instead adopt the threshold of $t_{\rm cool}/t_{\rm ff}\sim 10$ as representative of the threshold for condensation to occur.

Stacked profile of ratio between cooling and free fall times for gas with $T>T_0=3.6\times 10^4 \K$, the temperature of the cold stream. Near $r_{\rm shock}$, $t_{\rm cool}\lsim 3t_{\rm ff}$ implying that the shocked medium is likely unstable to condensation and precipitation. In the S Zone, $t_{\rm cool}\lsim t_{\rm ff}$ consistent with the formation of a cold, isothermal core at these radii.

tcool_tff shows the ratio $t_{\rm cool}(r)/t_{\rm ff}(r)$ as a function of $r/r_{\rm shock}$, stacked among the ten slices by taking the average and standard deviation.
At $r\lsim r_{\rm shock}$, $t_{\rm cool}/t_{\rm ff}\sim 3$, implying that while the shock may be stable to monolithic cooling, it is likely unstable to precipitation and condensation. The cooling radius, where $t_{\rm cool}/t_{\rm ff}\sim 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 thermal profiles).

§ 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 <cit.>. In this section, we wish to test this hypothesis in our simulations.

§.§ The Virial Theorem per-unit-length

An infinite, self-gravitating, collisionless and isolated cylinder in virial equilibrium (per-unit-length) obeys <cit.>

where $\mathcal{K}$ is the kinetic energy per-unit-length and $\Lambda$ is the line-mass. 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 will show in fil_dyn below (plasma_beta), magnetic support is negligible in the filament gas. We, therefore, use the virial theorem for an unmagnetized gaseous system. For such a system in time-independent equilibrium, the virial theorem states <cit.>

where $\mathcal{T}$ is the total thermal plus kinetic energy in the system, $\mathcal{T}_s$ is the confining pressure on the surface, including both thermal and ram pressure, and $\mathcal{W}$ is the total gravitational energy of the system. In integral form, these terms can be expressed as

𝒯 = ∫_V (1/2ρv^2+3/2P)dV,

𝒯_s = ∮_S (r⃗ ·π) ·dS⃗,

𝒲 = -∫_V (ρr⃗·∇⃗Φ) dV.
In vir_tau-vir_grav, $\rho$ is the gas density, ${\vec{v}}$ is its velocity, $P$ is its thermal pressure, $\Phi$ is the gravitational potential[Note that vir_grav depends only on the gravitational acceleration, $\vec{\nabla}\Phi$, and therefore does not depend explicitly on the zero point of the potential, which cannot be defined at $r=\infty$ for an infinite cylinder as we are assuming.], and $\overleftrightarrow{\pi}$ is the fluid pressure tensor given by

π_ij=ρv_iv_j + δ_ijP,

where $\delta_{ij}$ is the Kronecker delta.
We define the virial parameter

\begin{equation}
\label{eq:vir para}
\alpha_{\rm vir}=\frac{2(\mathcal{T}-\mathcal{T}_s)}{|\mathcal{W}|}.
\end{equation}

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 systems[See <cit.> for the importance of including an analogous surface pressure term for collisionless particles when estimating the virial equilibrium of dark matter haloes.]. When $\alpha_{\rm vir}>1$, the system is expanding due to kinetic plus thermal energy, while for $\alpha_{\rm vir}<1$ the system is collapsing due to the combined effects of surface pressure and gravity.

When evaluating vir_tau-vir_grav 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 the 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.
We thus approximate the gravitational acceleration as

∇⃗ Φ≃2GΛ_tot(<r)/rr̂,
with $\Lambda_{\rm tot}(<r)$ being the total line-mass interior to radius $r$. Likewise, the gas density and line-mass are related by

2πr ρ≃d Λ_gas/d r.
We, therefore, approximate the gravitational potential energy as

𝒲(r)≃-2GL∫_0^rΛ_totd Λ_gas/dr' dr',

where $L=30\kpc$ is the filament slice thickness. We note that this factor $L$ is only included for consistency with the formalism presented in <cit.>, and in practice, we normalize all the energy terms in vir_tau-vir_grav by $L$ to obtain energy per-unit-length before computing $\alpha_{\rm vir}$. We numerically evaluate the integral in W(r) by performing the following double-sum:

where the sum $\sum_{<r'}$ in the square brackets is taken over all particles (dark matter and stellar) and gas cells interior to radius $r'$, and the sum $\sum_{r'}$ is taken over all gas cells within the cylindrical shell of thickness $\Delta r'$ at $r'$.

We numerically evaluate the volumetric kinetic and thermal terms (vir_tau) as

𝒯(r)=∑_<r1/2m_gasv^2+∑_<r3/2k_B T≡𝒯_dyn(r)+𝒯_therm(r),
where $v=\sqrt{v_x^2+v_y^2+v_z^2}$ is the total velocity. Note that we include $v_z^2$ in $\mathcal{T}$ despite our assumption of filaments as infinite cylinders when computing $\mathcal{W}$. This is because even when evaluating the energy per-unit-length, such that there is no net velocity or acceleration along the $z$ direction, we still expect the velocity dispersion to be three-dimensional locally, such that $v_z^2$ remains an important contribution to the kinetic energy density of the system. We comment below on the effect of removing $v_z^2$ from the kinetic term.

The surface term, vir_taus, is composed of three surfaces for each radius $r$, one on the cylindrical shell and two others on either end of the cylindrical slice. However, the latter two surfaces are irrelevant for evaluating the virial equilibrium per-unit-length and we, therefore, do not consider them here. We comment below on the effect of including them in the calculation.
The surface term is thus given by

𝒯_s(r) = ∫_-L/2^L/2 ∫_0^2π (ρr v_r^2 + ρz v_rv_z + Pr ) r dϕdz.

The second term, namely the integral of $\rho z v_rv_z$, is not compatible with our assumption of an infinite cylinder, nor with the evaluation of $\mathcal{T}_{s}(r)$ per-unit-length. For an infinite cylinder in equilibrium per-unit-length, $v_{\rm z}$ would be constant (or zero). Likewise, grad_phi-W(r) are only valid for $\rho$ which does not depend on $z$, and in such cases $v_{\rm r}$ should also be independent of $z$. This term thus reduces to the integral of $z{\rm d}z$ 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 $\rho r v_r^2$ at all radii $r\gsim 0.15\,r_{\rm shock}$, because $\rho v_r v_z$ is indeed roughly independent of $z$, consistent with our assumed symmetry.

In practice, when numerically evaluating $\mathcal{T}_{s}$ we perform the following sum over gas cells,

𝒯_s=∑_rρr v_r^2V/Δr + ∑_r PrV/Δr≡𝒯_dyn,s + 𝒯_therm,s,
where $V$ is the volume of each gas cell and $\Delta r$ is the bin width. We have thus approximated ${\rm d}S \simeq V/\Delta r$. We have verified that this approximation is roughly independent of our choice of $\Delta r$, so long as this is within a factor of $\sim (2-3)$ of the typical cell size in the simulation.

§.§ Rapidly Flowing Gas on and off the Sheet

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 defined as $|x|<10\kpc\sim 0.22\,r_{\rm shock}$ (see thermal profiles). The dashed-orange line in each panel shows the results for all gas regardless of its position relative to the sheet. The value of $v_r$ in each bin is the mass-weighted average radial velocity, and the value for $\dot{M}$ in each bin is $\sum mv_r/\Delta r$, where $\Delta r$ is the width of the bin.
Prior to stacking, the $v_r$ profiles for each slice are normalized by the corresponding virial velocity at $r_{\rm shock}$, $v_{\rm vir}(r_{\rm shock})=(G\Lambda_{\rm shock})^{1/2}$ (v_vir), while the profiles of $\dot{M}$ were normalized by $\Lambda_(r_{\rm shock})v_{\rm vir}(r_{\rm shock})$.
In the left-hand panel, the green lines show profiles of $\pm v_{\rm vir}(r)/v_{\rm vir}(r_{\rm shock})$.
The sheet is inflowing with $|v_r|>v_{\rm vir}(r)$ outside the S zone, and is outflowing with small velocity within the S zone. The off-sheet material is inflowing in the T zone and outflowing in V zone, but always with $|v_r|<v_{\rm vir}(r)$.

Finally, before we calculate the corresponding profiles, we wish to remove gas that is rapidly accreting onto or outflowing from the filaments. Such gas is not expected to be in equilibrium within the filament, but rather to either be unbound or on its first infall, and will therefore bias our results <cit.>. We do this by applying a cut in the radial velocity of gas cells included in the energy terms described above. At each radius $r$, we calculate the corresponding “virial” velocity (from virial_ideal)

and only consider gas cells where $|v_r|<v_{\rm vir}(r)$. We apply this threshold to all terms appearing above, except for $\Lambda_{\rm 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 radial inflows 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, $\dot{M}$. 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 slices as everything within $|x|<10\kpc\sim 0.22\,r_{\rm shock}$ in the frame of the projection maps in thermal profiles and DM prop. In the left-hand panel, we also show the profile of $\pm v_{\rm 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 $|v_r|<v_{\rm vir}(r)$ for the virial parameter predominantly removes rapidly inflowing gas (with $v_r<0$) along the sheet in the T and V zones. The mean inflow velocity along the sheet exceeds $v_{\rm vir}$ at $r\gsim 0.3\,r_{\rm shock}$, while outside the sheet it only exceeds $v_{\rm vir}$ at $r>1.2r_{\rm 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 $r\lsim 2r_{\rm 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 sim_dynamics). At $r>2r_{\rm shock}$ the off-sheet gas is either accelerating or maintaining constant velocity. However, note that the gas does not completely stall at $r_{\rm 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 $r\lsim r_{\rm shock}$ (tcool_tff). 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 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-$z$ through their shock-heated haloes are seen in simulations to maintain roughly constant inflow velocities throughout the hot CGM <cit.>.

Within the V zone, at $0.2\lsim r/r_{\rm shock}\lsim 0.6$, the off-sheet gas is outflowing with an average outflow velocity smaller than $v_{\rm vir}$. This seems to be caused by a strong quadrupolar vortex structure that develops in this region, as discussed in toy_model. These vortices cause gas outside the sheet to swirl around and can lead to outward 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 <cit.>, 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, the mass accretion onto the S zone is dominated by flow outside the sheet in the PS zone, has equal contribution from sheet and off-sheet material in the T zone, and is dominated by the sheet in the V zone.

Within the S zone, the off-sheet component is not well defined due to our crude definition of the sheet as everything within $|x|<10\kpc\sim 0.22\,r_{\rm shock}$. The sheet gas in this region is outflowing, albeit at very small velocities, such that the outflow is better seen in the mass flux panel. This seems to be primarily due to the impact parameter of the sheet which can generate positive radial velocities very close to the filament centre, along with vortical motions that develop in the central regions (see toy_model). However, we note that this outflow seems transient and is not seen in preliminary examinations of additional snapshots roughly $55\Myr$ before and after our fiducial snapshot (see additional discussion of this point in sim_dynamics).

Radial profile of the virial parameter, $\alpha_{\rm vir}$ (black line), along with the ratios of the various energy terms in the numerator of vir para to the gravitational energy, $|\mathcal{W}(r)|$. Specifically, these are the volumetric kinetic and thermal terms, $\mathcal{T}_{\rm dyn}$ (green solid line) and $\mathcal{T}_{\rm therm}$ (magenta solid line), and the surface kinetic and thermal terms, $\mathcal{T}_{{\rm dyn},s}$ (green dashed line) and $\mathcal{T}_{{\rm therm},s}$ (magenta dashed line). We see that $\alpha_{\rm vir}$, $\mathcal{T}_{\rm therm}/|\mathcal{W}|$, and $\mathcal{T}_{{\rm therm},s}/|\mathcal{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. The filaments are thus in approximate virial equilibrium per unit-length with $r_{\rm vir}\sim r_{\rm shock}$.

§.§ Virial Equilibrium and the Virial Radius

Using vir para, W(r)-tau_sum, and taus_sum, we evaluate $\alpha_{\rm vir}$ for each slice as a function of radius for gas with $|v_r|<v_{\rm vir}(r)$. We discuss the sensitivity of our results to this threshold below. We then stack all ten slices by weighting $\alpha_{\rm vir}(r)$ by $|\mathcal{W}(r)|$. This is equivalent to treating all ten slices as one long cylinder when computing the virial parameter per-unit-length, in accordance with our modelling of the filaments as infinite cylinders.

We present the stacked profile of $\alpha_{\rm vir}(r)$ in virial parameter in black, along with the ratios of the four energies in the numerator of $\alpha_{\rm 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 $\alpha_{\rm vir}(r)$ is useful to locate the radius where $\alpha_{\rm vir}\sim 1$, which may be considered the “virial radius” of the system. The profile monotonically decreases from $\alpha_{\rm vir}\sim 10$ in the S zone where there is indeed a net outflow (radial inflows), to $\alpha_{\rm vir}\sim 1$ in the T and outer V zones at $0.5\lsim r/r_{\rm shock}\lsim 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 $r_{\rm vir}\sim r_{\rm shock}$ for our filament sample.

We further note that the volumetric and surface thermal terms are both $\sim 1$ in the T zone, where the gas temperature is roughly equal to the effective temperature of the dark matter as defined by its velocity dispersion (DM_prop). This implies that the dark matter is also expected to be in approximate virial equilibrium (per-unit-length) within $r_{\rm shock}$. This is consistent with previous estimates of the virial radii of filaments feeding massive haloes at high-$z$ <cit.>. For a filament feeding a $10^{12}\,\msun$ halo at $z\sim 4$, these authors predict $r_{\rm vir}\sim 55\kpc$, slightly larger than our average $r_{\rm shock}\sim 45\kpc$ and within the range of $r_{\rm shock}$ values among our ten slices (rho_0 and T_0 scatterings). Note that this is different from the case of galaxy clusters, where the shock radius is found to be larger than the virial radius <cit.>.

Removing $v_z^2$ from $\mathcal{T}$ and $\mathcal{T}_s$ reduces $\alpha_{\rm vir}$ by $\sim 30\%$, such that $\alpha_{\rm vir}\sim 0.6$ near $r_{\rm shock}$. This is consistent with approximate equipartition between the three dimensions of kinetic energy. Our inclusion of $v_z^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 $\mathcal{T}_s$ results in $\alpha_{\rm vir}<0$ at all radii. The exact numerical value of $\alpha_{\rm 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 gas with $|v_r|>2v_{\rm vir}$ has no affect on our results showing that the exclusion is dominated by very rapidly flowing gas. Changing our threshold in $v_r$ to only exclude inflowing gas with $v_r<-v_{\rm vir}$ while including rapidly outflowing gas retains $\alpha_{\rm vir}\sim 1$ in the
T and outer V zones,

showing that our estimation of a virial radius at $r_{\rm vir}\sim r_{\rm shock}$ is robust. At smaller radii, where the system is far from equilibrium anyway based on virial parameter, $\alpha_{\rm vir}$ becomes negative due to a much larger surface kinetic term. Including rapidly inflowing gas as well results in negative $\alpha_{\rm vir}$ everywhere, consistent with our interpretation that gas currently inflowing towards the filament along the sheet is not yet in equilibrium.

In summary, the gas within $r_{\rm 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 $r_{\rm vir}\sim r_{\rm shock}$. Within the T zone, the gas is roughly 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 radial inflows.

§ 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 summation by describing the mathematical formalism of our force decomposition, and present results from our simulation in sim_dynamics. In toy_model, we interpret the simulation results using an analytic toy model for the internal filament dynamics. In Appendix sim_numerics, we discuss the robustness of the results presented in sim_dynamics to our numerical method.

§.§ Force Decomposition

§.§.§ Basic Framework

Radial profile of the plasma $\beta$ parameter, the ratio of thermal to magnetic pressure, stacked among our filament slices. We see that $\beta$ evolves from $\gsim 10^4$ at $r\sim 0.1\,r_{\rm shock}$ to $\gsim 10^7$ at $r\sim r_{\rm shock}$, implying that magnetic fields are not dynamically important in the filaments.

The gas dynamics in the simulation evolves according to the ideal magnetohydrodynamics (MHD) equations. In particular, the Eulerian momentum equation reads

v⃗t+(v⃗·∇⃗ ) v⃗ = -1/ρ ∇⃗ (P+B^2/8π) + 1/4πρ(B⃗·∇⃗)B⃗ - ∇⃗ Φ,
where $\vec{v}$ is the fluid velocity[$\vec{v}$ in MHD represents the physical velocity, including the Hubble flow in an expanding Universe. In practice, the simulation stores the peculiar velocity, $\vec{v}_{\rm pec}=\vec{v}-H(z)\vec{r}$ with $H(z)$ the Hubble parameter at redshift $z$, for which an additional term of $(-H(z)\vec{v})$ representing the so-called `Hubble drag' must be added to the right-hand side of MHD. However, over the scales of interest, $r\lsim r_{\rm shock}\sim 50\kpc$, the Hubble flow $H(z)r\lsim 10\kms$ is negligible compared to the peculiar velocity. Likewise, the Hubble drag is negligible compared to the other forces, at the level of a few percent at most, and has no affect on our results. We hereafter neglect this subtlety and treat peculiar and physical velocities interchangeably.], $\rho$ is its density, $P$ its thermal pressure, $\Phi$ the gravitational potential, and $\vec{B}$ the magnetic field. The dynamical importance of the magnetic field is often quantified by the plasma $\beta$ parameter, which is the ratio of thermal to magnetic pressure,

β≡P_th/P_B = 8πP/B^2,
where $P_{\rm B}=B^2/(8\pi)$ is the magnetic pressure. We show the profile of $\beta$, stacked among our ten filament slices, in plasma_beta. Near the filament centre, $\beta\gsim 10^4$, while near the shock radius $\beta\gsim 10^7$. Magnetic pressure is thus negligible compared to 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 $r_{\rm shock}\lsim 50\kpc$, and the magnetic pressure varies over a characteristic length-scale of order the smallest cell size, $\Delta\gsim 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 MHD to the Euler equation,

v⃗t+(v⃗·∇⃗ ) v⃗ = -1/ρ ∇⃗ P - ∇⃗ Φ.
We note that observations of low-$z$ filaments in the vicinity of galaxy clusters do reveal the presence of magnetic fields <cit.>. 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.

We follow <cit.> in using euler 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 here generalize their method to study the dynamics of gas in intergalactic filaments assuming cylindrical symmetry. We describe the method below, both general considerations as well as formulae 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,

M_grav≡1/4πG∮_S∇⃗Φ·dS⃗ = M_tot,
where $M_{\rm tot}$ is the total mass enclosed by an imaginary closed surface $S$ over which the surface integral is taken. Combining euler and gauss, we have

\begin{equation}
\label{eq:summation}
M_{\rm grav}=-\frac{1}{4\pi G}\oint_{S}\left(\pdv{\vec{v}}{t}+\left(\vec{v}\cdot\vec{\nabla}\right)\vec{v}+\frac{1}{\rho}\vec{\nabla} P\right)\cdot{\rm d}\vec{S}.
\end{equation}

Each of the terms in the integral on the right-hand-side of summation represents an acceleration. When combined, these terms must balance the gravitational acceleration. From the left-hand side of summation we learn that these terms can be thought of as “mass terms”, contributing to the total mass interior to the surface $S$. We thus define

\begin{equation}
\label{eq:Mtot}
M_{\rm grav}=M_{\rm accel}+M_{\rm therm}+M_{\rm inertial},
\end{equation}

\begin{equation}
\label{eq:acceleration term}
M_{\rm accel} \equiv -\frac{1}{4\pi G}\oint \pdv{\vec{v}}{t}\cdot {\rm d}\vec{S},
\end{equation}

is the acceleration term representing the temporal change of the velocity field normal to the bounding surface, [Note that in our Eulerian representation, this represents a temporal change of the velocity field inside the filament, rather than the local acceleration experienced by a test particle as it moves through the filament.]

\begin{equation}
\label{eq:thermal term}
M_{\rm therm}\equiv -\frac{1}{4\pi G}\oint \frac{1}{\rho}\vec{\nabla} P\cdot {\rm d}\vec{S},
\end{equation}

is the thermal term representing the support of the filament against gravity by thermal pressure gradients, and

\begin{equation}\label{eq:inertial term}
M_{\rm inertial}\equiv -\frac{1}{4\pi G}\oint \left(\vec{v}\cdot\vec{\nabla} \right)\vec{v}\cdot{\rm d}\vec{S},
\end{equation}

is the inertial term. As we demonstrate below and expand upon in toy_model, 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

In order to estimate the overall dynamical equilibrium of the filament and the relative contribution of each term to the filament support against gravity, we analyze acceleration term-inertial term as a function of radius for each filament slice and compare their sum (Mtot) to the actual mass interior to radius $r$ in the simulation. For each radius $r$, we assume a Gaussian surface which is a cylinder of radius $r$ and length $L=30\kpc$, the thickness of each filament slice. We comment on the impact of $L$ on our results in sim_dynamics. Unlike the discussion of virial equilibrium presented in virial, 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 $S_0$ with ${\rm d}\vec{S}_0={\rm d}S_0{\hat{r}}$, and two others on either end of the cylinder, hereafter $S_{\pm}$, with ${\rm d}\vec{S}_{\pm}={\rm d}S_{\pm}{\pm\hat{z}}$. We then have

\begin{equation}
\label{eq:S1}
\int_{S_0}\text{d}S=r\int_{0}^{2\pi}\text{d}\phi\int_{-L/2}^{L/2}\text{d}z,
\end{equation}

\begin{equation}
\label{eq:S2}
\int_{S_{\pm}}\text{d}S=\left(\int_0^r r'\text{d}r'\int_{0}^{2\pi}\text{d}\phi\right)\bigg|_{z=\pm L/2}.
\end{equation}

Now let

∫_ℛdS≡∫_S_0dS,   ∫_𝒜dS≡∫_S_+dS-∫_S_-dS,
where ${\mathcal R}$ stands for radial and ${\mathcal A}$ stands for axial. Each of the integrals in acceleration term-inertial term is thus decomposed into two parts, $\oint_{\rm cylinder}\text{d}S=\int_{\mathcal R}\text{d}S+\int_{\mathcal A}\text{d}S.$

Using this notation, we can explicitly write out all the terms in acceleration term-inertial term as follows:

M_accel,r(r)=-1/4πG∫_ℛv_rt dS,

M_accel,a(r)=-1/4πG∫_𝒜v_zt dS,

M_therm,r(r)=-1/4πG∫_ℛ1/ρPr dS,

M_therm,a(r)=-1/4πG∫_𝒜1/ρPz dS,

M_inertial, r(r)=-1/4πG∫_ℛ(v_rv_rr+v_ϕ/rv_rϕ+v_zv_rz-v_ϕ^2/r) dS,

M_inertial, a(r)=-1/4πG∫_𝒜(v_rv_zr'+v_ϕ/r'v_zϕ+v_zv_zz) dS.
Note that all the integrals in accel-axial inertial term are functions of radius, $r$, which we have highlighted on the left-hand-side of each of these equations. accel, Mtherm_r, and radial inertial term represent radial accelerations towards or away from the filament axis, while accel_2, Mtherm_a and axial inertial term represent axial accelerations along the filament axis.

§.§.§ Further decomposing the Inertial Terms

The last term in the integrand of radial inertial term represents the centrifugal acceleration, and is hereafter referred to as the rotation term,

M_rot(r)=1/4πG∫_ℛv_ϕ^2/r dS.
When evaluating rotation, we further decompose the azimuthal velocity at radius $r$ into mean and residual components. The mean rotation is given by

v_ϕ,mean(r)=1/2πr L∫_-L/2^L/2∫_0^2π v_ϕrdϕdz.
Note that mean_rot represents the volume-weighted average azimuthal velocity at radius $r$, rather than the mass- or density-weighted average. This is necessary for consistency with rotation. The residual azimuthal velocity in each cell is then given by

\begin{equation}
\label{eq:rand_rot}
v_{\phi, \rm res}(r,\phi,z)=v_\phi(r,\phi,z)-v_{\phi,\rm mean}(r).
\end{equation}

With these definitions, it is straightforward to show that

M_rot(r)=14πG_ℛ(v_ϕ,mean^2r + v_ϕ,res^2r)dS

≡M_rot, mean(r)+M_rot, res(r),

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 radial inertial term can be interpreted in two ways. The first, which we will demonstrate in detail in toy_model, is as additional fictitious forces resulting from the fluid motion. The second is as contributing 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 radial inertial term represent 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 mean_rot-rand_rot, such that the mean is only a function of $r$, then it is straightforward to show that these terms can be written as follows:

M_ram, r(r)=-1/4πG∫_ℛv_r,meanv_r,meanr dS,

M_turb, rr(r)=-1/4πG∫_ℛv_r,resv_r,resr dS,

M_shear, rϕ(r)=-1/4πG∫_ℛ (v_ϕ, mean/r+v_ϕ, res/r)v_r,resϕ dS ,

M_shear, rz(r)=-1/4πG∫_ℛ (v_z, mean+v_z, res)v_r,resz dS.
We expect the first terms in Mturb_p and Mturb_z to be negligible compared to the second, because $v_{\phi, {\rm mean}}(r)$ and $v_{z, {\rm mean}}(r)$ can be taken outside the integral and the average of the derivatives of $v_{r,{\rm res}}$ should be small. Analogous divisions[We note that this distinction between ram and turbulent pressure is not unique. If one defines the mean motion with respect to some local smoothing kernel around each fluid element (e.g., as was done in the study of CGM dynamics by ), then the mean motion is no longer only a function of $r$ and all three terms contribute to both ram and turbulent pressure.] into ram, turbulent, and shear forces in the axial direction can be applied to the three terms in axial inertial term.

Finally, we also divide the gravitational term (Mtot) into its radial and axial components:

M_grav, r(r) = 1/4πG∫_ℛ ∂Φ/∂r dS,

M_grav, a(r) = 1/4πG∫_𝒜 ∂Φ/∂z dS.
These represent the radial and axial gravitational fields in the filament, respectively. Note that since the gravitational field is $\vec{g}=-\vec{\nabla}\Phi$, $M_{\rm grav, r}$ and $M_{\rm grav, a}$ are defined to be positive when the radial and axial gravitational fields are directed inwards towards the filament centre. However, the thermal and inertial terms in Mtherm_r-axial inertial term and decomposition of the rotation term-Mturb_z are all 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

M_grav, r = M_accel, r + M_inertial, r + M_therm, r,
with an analogous equation for the axial direction.

§.§.§ Numerical Calculation

Using our simulations, we numerically evaluate accel-Mgrav_a as a function of radius, $r$, for each filament slice. Full details of how this is done, including convergence and robustness tests, are provided in Appendix sim_numerics. 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)r_{\rm shock}$ with a typical cell-size of $\Delta\sim 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 $\lsim 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 $\sigma=0.5\kpc$. Evaluating both sides of gauss, we find that our numerical estimate of $M_{\rm grav}$
deviates from $M_{\rm tot}$ by at most $\lsim 5\%$ at each radius.

After computing the radial profiles of each mass-term for each filament slice, we normalize these by $M_{\rm tot}(<r)$ and stack the ten filament slices by taking the average of each normalized mass term in each radial bin. While it is the ratio of each mass term to $M_{\rm grav}(r)$ that tells us the contribution of the respective force to the filament support at that radius, we nonetheless normalize the mass terms by $M_{\rm tot}(<r)$, which is a smooth, monotonically increasing function that is within $\lsim 5\%$ of $M_{\rm grav}(r)$, which exhibits small fluctuations. We hereafter use the notations $M_{...}/M_{\rm tot}$, $M_{...}/M_{\rm grav}$, and $F_{...}/F_{\rm 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 accel-accel_2, which represent a temporal change to the velocity field in the filament. Directly evaluating these terms requires using multiple snapshots. As detailed in Appendix sim_numerics, 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 summation_radial 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 dynamical stabilities by analyzing where the velocity field in the filaments is in steady state, with $M_{\rm accel}\propto \partial \vec{v}/\partial t\simeq 0$, and where it is changing in time. We then analyze in detail the various components of the thermal and inertial terms in detailed force decomposition, to better characterize the force balance and dynamical state in the filaments. Throughout, we emphasize key trends and physical processes in the three radial zones.

§.§.§ Overall Force Balance

In summation method, we assess to what extent thermal and inertial forces balance gravity within the filament, as a function of $r/r_{\rm shock}$. The dashed dark-blue (hereafter navy) and pink lines show the stacked profiles of $M_{\rm grav, r}/M_{\rm tot}$ and $M_{\rm grav, a}/M_{\rm 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 Mgrav_r-Mgrav_a, 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 virial), 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-$\sigma$ standard deviation among our ten slices of the radial profile of $M_{\rm grav}/M_{\rm tot} = (M_{\rm grav, r}+M_{\rm grav, a})/M_{\rm tot}$. This ratio is $\simeq 1$ everywhere, as Gauss's theorem tells us it must be (gauss). While our numerical method introduces a $\sim 5\%$ error at $r/r_{\rm shock}\lsim 0.3$ and $\gsim 0.8$ (see also Grad Phi together), this is negligible compared to the variations in $M_{\rm grav, r}$ and $M_{\rm grav, a}$.

Force equilibrium in the filaments. We compare the sum of thermal and inertial forces (solid lines) to the gravitational force (dashed lines), in both the radial and axial directions, as a function of $r/r_{\rm shock}$. The magenta dashed line and cyan shaded region show $M_{\rm grav}(r)/M_{\rm tot}(<r)$,
which is $\simeq 1$ everywhere with a $\lsim 5\%$ error at $r/r_{\rm shock}\lsim 0.3$ and $\gsim 0.8$ due to our numerical method. Navy and pink dashed lines show the ratios $M_{\rm grav,r}(r)/M_{\rm tot}(<r)$ and $M_{\rm grav,a}(r)/M_{\rm tot}(<r)$, 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. The solid navy and pink lines show $(M_{\rm therm,r}+M_{\rm inertial,r})/M_{\rm tot}(<r)$ and $(M_{\rm therm,a}+M_{\rm inertial,a})/M_{\rm tot}(<r)$, respectively (Mtherm_r-axial inertial term). 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 (accel and/or accel_2 are equal to 0). If, however, the solid lines are greater (smaller) than the corresponding dashed lines, this then indicates a negative (positive) temporal acceleration at this radius. Over a fairly large radial range in the V and T zones, $r/r_{\rm shock}\sim (0.25-0.40)$ and $(0.60-0.80)$, there is a dynamical steady state in the radial direction. $\partial v_r/\partial t > 0$ at $r/r_{\rm shock}\sim (0.40-0.60)$ and $(0.80-1.0)$, while $\partial v_r/\partial t < 0$ in the S zone at $r/r_{\rm shock}\lsim 0.25$. In the axial direction, there is no dynamical steady state at $r/r_{\rm shock}\gsim 0.2$.

The solid navy and pink lines show the stacked profiles of $(M_{\rm therm, r}+M_{\rm inertial,r})/M_{\rm tot}$ and $(M_{\rm therm, a}+M_{\rm inertial,a})/M_{\rm tot}$, respectively (Mtherm_r-axial inertial term). In a steady state where thermal and inertial forces balance gravity and the filament velocity field is constant in time, these must equal $M_{\rm grav, r}/M_{\rm tot}$ (summation_radial) and $M_{\rm grav, a}/M_{\rm tot}$, respectively.
If they are larger (smaller) than the corresponding $M_{\rm grav,r}$ and $M_{\rm grav,a}$ terms at some radius $r$, then the velocity field at that radius accelerates outwards (inwards) in the corresponding direction.

In the radial direction, the filament exhibits a dynamical steady-state with no temporal acceleration throughout the inner halves of both the V zone and the T zone, $0.25\lsim r/r_{\rm shock}\lsim 0.4$ and $0.6\lsim r/r_{\rm shock}\lsim 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 $r\gsim 0.8r_{\rm shock}$, there is an outward 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 thermal to identify the shock radius in our ten filament slices at $z\sim 4.05$ and $z\sim 3.82$, roughly $55\Myr$ before and after our fiducial snapshot at $z\sim 3.93$. The typical $r_{\rm shock}$ increases from $\sim 36\kpc$ to $\sim 51\kpc$ during this interval,
corresponding to a shock velocity of $v_{\rm shock}\sim 130\kms$, while the sound speed in the post-shock gas is $c_{\rm s}\sim 85\kms$. The outer T zone near $r_{\rm shock}$ has thus not reached equilibrium after being hit by the shock, explaining the positive acceleration in this region.

In the V zone, there is an outward radial acceleration at $0.4\lsim r/r_{\rm shock}\lsim 0.6$. This is the region where the radial velocity of sheet material begins decreasing and the off-sheet material begins moving outwards (radial inflows). As we will show in toy_model 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.

Detailed force balance in the filaments. Left panel: We compare the radial thermal and inertial forces (red and green solid lines) to the radial gravitational force (navy dashed line, as in summation method), and the axial thermal and inertial forces (orange and yellow solid lines) to the axial gravitational force (pink dashed line). Centre panel: We show the four constituent forces comprising the radial inertial force (radial inertial term), 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 the panels.
Right panel: We show the three constituent forces comprising the axial inertial force (axial inertial term). 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 inwards gravitational force, while the inertial forces are directed inwards and dominated by ram pressure, which mostly offsets the excess thermal pressure. In the V Zone, the radial thermal force is directed outwards and decreases from $\sim 1.5|F_{\rm grav,r}|$ to $\sim 0.6|F_{\rm grav,r}|$ as $r$ decreases. The radial inertial forces are directed outwards with a roughly constant amplitude of $\sim 0.6|F_{\rm grav,r}|$, and are dominated by centrifugal forces. The mean rotation has a roughly constant amplitude of $\sim |F_{\rm 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 $r\gsim 0.15\, r_{\rm shock}$ and $r\lsim 0.15\,r_{\rm shock}$, respectively, resulting in a net inwards acceleration in this region which damp the outflows seen in radial inflows. 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.

In the S zone, at $r\lsim 0.25r_{\rm shock}$, there is an inwards acceleration. This corresponds to the outer boundary of the isothermal core (thermal profiles), where the thermal pressure decreases locally due to strong cooling (tcool_tff). However, we do not interpret this as evidence that the cold-stream is collapsing with no support against gravity. Rather, recall that both the mean radial velocity and the mass flux are directed outwards in the S zone (radial inflows), which appears to be a manifestation of the impact parameter of the flow within the sheet with respect to the filament centre. This outwards velocity is also consistent with the very large virial parameter in this region (virial parameter). A negative radial acceleration term in this region would then imply that this outwards motion is decelerating as the cold stream relaxes. Indeed, at $z\sim 3.82$ roughly $55\Myr$ after our fiducial snapshot, the velocity in this region is no longer outwards and is close to 0. The detailed dynamics and kinematics of the central streams may thus be highly time-dependent.

Finally, we note that in the axial direction, there is no steady-state at $r\gsim 0.2\,r_{\rm shock}$, as the filament gets continually stretched by tidal forces and material accelerates towards the massive haloes on either end. Note that the acceleration term is positive in the S zone and negative outside of it. This suggests that the velocity of the cold streams towards the haloes increases with time, while that of the hot filament gas decreases with time. We discuss this further below.

§.§.§ Detailed Force Decomposition

In summation method another we show the individual profiles of all the forces that went into the curves shown in summation method, to asses their individual contributions to the dynamical state of the filament. On the left, the solid lines show $F_{\rm therm, r}/F_{\rm grav}$ (red), $F_{\rm inertial, r}/F_{\rm grav}$ (green), $F_{\rm therm, a}/F_{\rm grav}$ (orange) and $F_{\rm inertial, a}/F_{\rm grav}$ (yellow). The dashed lines are as in summation method, showing $F_{\rm grav,r}/F_{\rm grav}$ (navy), $F_{\rm grav,a}/F_{\rm grav}$ (pink), and $M_{\rm grav}/M_{\rm tot}$ (magenta). First, we note that the four non-gravitational forces obtain both positive and negative values. As discussed following Mgrav_r-Mgrav_a, a positive value implies an acceleration outward counteracting collapse, while a negative value implies an inwards acceleration. In the centre and right panels, we show the profiles of all the individual components that comprise the radial and axial inertial terms, respectively (radial inertial term-axial inertial term). 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 $r\gsim 0.6r_{\rm shock}$, $F_{\rm therm,r}$ is directed outwards and is stronger than the inwards gravitational force, monotonically increasing from $\sim |F_{\rm grav,r}|$ at $r\sim 0.6\,r_{\rm shock}$ to $\sim 4|F_{\rm grav,r}|$ at $r\sim r_{\rm shock}$. This excess in the outwards thermal pressure force is partially balanced by a strong inwards inertial force, $F_{\rm inertial,r}$, which grows in amplitude from $\sim (0.3-2)|F_{\rm grav,r}|$ in the same radial range. The inwards inertial force in this region is dominated by the $-v_r \partial v_r/\partial r$ term (centre panel, red lines), which grows in amplitude from $\sim (1.5-3)|F_{\rm grav,r}|$. As discussed in summation, this term contributes to the radial components of both the ram pressure force (Mram, dot-dashed red line) and the turbulent pressure force (Mturb_r, solid red line). The ram pressure force can be visualized by examining the radial velocity profiles in radial inflows. While the velocity in the sheet is roughly constant in the T zone, the inwards velocity outside the sheet decreases from $\sim v_{\rm vir}$ at $r_{\rm shock}$ to $\sim 0$ at $0.6\, r_{\rm shock}$, resulting in $-v_r \partial v_r/\partial r<0$.
In the outer T zone, $0.8\lsim r/r_{\rm shock}$ where there is a positive radial acceleration term (summation method), the ram pressure force dominates over the turbulent pressure force which actually becomes negligibly small at $r\gsim r_{\rm shock}$.
In the inner T zone, $0.6\lsim r/r_{\rm shock}\lsim 0.8$ where the radial velocity field is in steady-state (summation method), the ram pressure and turbulent pressure forces are comparable. The two shear forces, $r\phi$ and $rz$ (Mturb_p-Mturb_z, orange and blue lines in the centre panel, respectively), are dominated by residual motions (solid lines), while 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. To summarize, strong thermal pressure gradients in the T zone support the filament against both gravity and ram pressure from inflowing gas. At $r\lsim 0.8\, r_{\rm shock}$, the gas is roughly in a steady state, while at $r\gsim 0.8\, r_{\rm shock}$ the velocity field is accelerating outwards in response to the shock.

The V zone, $0.25 \lsim r/r_{\rm shock} \lsim 0.6$, is the only zone where both $F_{\rm therm,r}$ and $F_{\rm inertial,r}$ are directed outwards. In the outer V zone, $0.4 \lsim r/r_{\rm shock} \lsim 0.6$ where there is a positive radial acceleration (summation method), the thermal pressure is dominant and still stronger than the radial gravitational force and comparable in amplitude to its value in the inner T zone. In the inner V zone, $0.25 \lsim r/r_{\rm shock} \lsim 0.4$ where the radial velocity field is in steady-state (summation method), the thermal pressure term decreases and becomes comparable to the inertial term, which maintains a value of $\lsim 0.6|F_{\rm 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 primarily due to a decreasing ram pressure force inwards and an increasing centrifugal force outwards. The mean rotation term increases from $\sim 0$ at $r\sim r_{\rm shock}$ to a constant $\sim |F_{\rm grav,r}|$ throughout the S and inner V zones, at $r\lsim 0.4 r_{\rm shock}$. The residual rotation term is $\sim 1.5|F_{\rm grav,r}|$ at $r\sim r_{\rm shock}$, and increases monotonically towards smaller $r$. Throughout the V and S zones, at $r < 0.6 r_{\rm shock}$, this is mostly offset by the inwards $r\phi$ shear force, $-(v_{\phi}/r)\partial v_r/\partial \phi$, with a smaller contribution from ram/turbulent pressure. In toy_model, 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 summarize the situation in the outer V zone, at $0.4\,r_{\rm shock}\lsim r\lsim 0.6\,r_{\rm 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 toy_model below).

In the inner V zone, at $0.25 r_{\rm shock}\lsim r \lsim 0.4 r_{\rm shock}$ where the filament velocity field is in a steady state with no acceleration (summation method), the thermal force continues to be directed outwards but is now weaker than the inwards gravitational force, $\sim (0.3-0.6)|F_{\rm grav,r}|$. The inertial term continues to be directed outwards as well, maintaining a constant amplitude of $\sim 0.6|F_{\rm 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 $\sim |F_{\rm grav,r}|$. The residual rotation term continues rising towards smaller $r$, with amplitudes in the range $\sim (4-7)|F_{\rm grav,r}|$, though this is again compensated for by shear and turbulent pressure, (see toy_model below). This trio yields a slight net inwards acceleration, that lowers the total inertial term compared to the mean rotation. To summarize the situation in the inner V zone, at $0.25\,r_{\rm shock}\lsim r\lsim 0.4\,r_{\rm shock}$, the filament is kept in force equilibrium with roughly equal contributions from thermal pressure forces and centrifugal forces in the support against gravity.

In the S zone, at $r\lsim 0.25r_{\rm shock}$, there is a net inwards acceleration term (summation method). This is primarily driven by the thermal pressure force becoming negative at $r\sim (0.15-0.25)r_{\rm shock}$, near the outer boundary of the isothermal core where cooling is maximal, and the inertial term becoming negative at $r\lsim 0.15r_{\rm shock}$ due to a slight drop in the centrifugal forces. As stated in dynamical stabilities, this inwards acceleration coincides with the region where the mean radial velocity and the total mass flux are directed outwards (radial inflows), and damps the outflow. This suggests that the dynamics of the isothermal cores of high-$z$ filaments, representing the cold streams feeding massive high-$z$ galaxies (thermal), 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 summarize the radial force balance throughout the filament, emphasizing 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, thermal pressure forces become comparable to 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 which damps outflows in this region.

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. $F_{\rm inertial,a}$ is dominated by axial ram pressure forces encapsulated in the term $-v_z\partial v_z/\partial z$, 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 cold_streams.

§.§ A Toy Model for the Filament Velocity Field

In sim_dynamics we saw that, in addition to thermal pressure gradients, the filament's 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, $(v_\phi/r)\partial v_r/\partial \phi$ and $v_r\partial v_r/\partial r$. 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 sim_dynamics. We begin in kepler by discussing the very simple case of an elliptical orbit in a Keplerian potential as an instructive example. In filament's vorticity we discuss the vorticity structure of the filaments in our simulation, highlighting a characteristic quadrupolar structure. Finally, in vortex, 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 $M$. The equation of motion in the radial direction, with respect to the central mass, is

where $r$ is the distance to the central mass and $\phi$ is the azimuthal angle. The azimuthal and radial velocities are $v_{\phi}=r\dot{\phi}$, and $v_r=\dot{r}$. The first term on the left-hand side of radial motion is thus $v_\phi^2/r$ while the second is $d v_r/dt$. The latter term can be expressed in two possible ways, depending on whether the orbit is parametrized with respect to the radius, $r$, or azimuthal angle[An ellipse is a one-dimensional curve with one independent parameter.], $\phi$.

If we parametrize the orbit in terms of $r$ we have

\begin{equation}
\label{eq:r derivative}
\dv{}{t}=\dot{r}\dv{}{r}=v_r\dv{}{r}.
\end{equation}

Inserting this into radial motion yields

\begin{equation}
\label{eq:r parametrization}
\frac{v_\phi^2}{r}-v_r\dv{v_r}{r}=\frac{GM}{r^2},
\end{equation}

showing that in the frame of the particle, gravity is indeed balanced by the two relevant components of the inertial term in radial inertial term, which act as fictitious forces along the orbit. We can explicitly evaluate r parametrization for the Keplerian orbit by recalling that both the specific angular momentum, $l=v_{\phi}r$, and the specific energy, $\mathcal{E}=-GM/(2a)$ with $a$ the semi-major axis, are conserved. Therefore,

\begin{equation}
\label{eq:vphi in r}
\end{equation}

\begin{equation}
\label{eq:vr in r}
\end{equation}

One can easily verify that vphi in r-vr in r obey r parametrization.

If we parametrize the orbit in terms of $\phi$, r derivative-vr in r become

\begin{equation}
\dv[]{}{t}=\dot{\phi}\dv{}{\phi}=\frac{v_\phi}{r}\dv{}{\phi},
\end{equation}

\begin{equation}
\label{eq:phi parametrization}
\frac{v_\phi^2}{r}-\frac{v_{\phi}}{r}\dv{v_r}{\phi}=\frac{GM}{r^2},
\end{equation}

\begin{equation}
\label{eq:vphi in phi}
\end{equation}

\begin{equation}
\label{eq:vr in phi}
\end{equation}

where $e$ is the eccentricity of the orbit, which is a conserved quantity. Again, one can easily verify that vphi in phi-vr in phi obey phi parametrization and that gravity is balanced by the two relevant components of inertia term in radial inertial term, 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 radial inertial term act as fictitious forces that balance gravity.

§.§.§ The Velocity Structure of Filaments

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 \kms$, with lengths scaling as $l=30 \sinh^{-1}\left(v/30\kms\right)$,
such that a $30 \kms$ arrow is $\sim 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 ($r\sim 0.6\,r_{\rm shock}$, orange circle). These are highlighted with schematic eddies to guide the eye.

Vorticity along the filament axis. Left: Projection map of the vorticity along the filament axis, mass-weighted along the line of sight, normalized by $v_{\rm circ}/r_{\rm shock}$, and stacked among the 10 filament slices. The same quadrupolar structure visible in the single slice shown in streamplot 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 $\omega_z^2$ and $\omega_{\rm tot}^2$ 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 show a preference for vorticity to be aligned with their axis, with a typical ratio of $\gsim 0.6$, compared to $0.5$ which is characteristic of a random distribution.

In order to model the inertia forces seen in summation method another using a similar analysis to that presented in kepler, we must first characterize the gas velocity field within filaments. In streamplot, we show the velocity field in slice 2 (see slices), 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 thermal profiles. The arrows represent the 2D gas velocity field perpendicular to the filament's axis, $v_x{\hat{x}}+v_y{\hat{y}}$, mass-weighted-averaged 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 $y$-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. Sheet gas, on the other hand, maintains its approximate trajectory due to its higher density which destabilizes the shock there <cit.>. Inside $r_{\rm shock}$, the velocity is broadly characterized by four vortices centred near the transition between the T and V zones at $r\sim 0.6\,r_{\rm 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 <cit.>, 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 the dark matter, an additional source of vorticity for the gas is the 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 vorticity (left), we show a stacked projection map of $\omega_z$, the component of the vorticity parallel to the filament axis. In each slice, we compute the mass-weighted average value of $\omega_z$ along the line-of-sight, and weight this by $v_{\rm circ}(r_{\rm shock})/r_{\rm shock}$ prior to stacking, where $v_{\rm circ}(r)=\sqrt{2G\Lambda_{\rm tot}(r)}=\sqrt{2}v_{\rm vir}$ is the circular velocity of the cylindrical 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[Note that while all of our previous analysis was symmetric under the transformation $z\rightarrow -z$, not so $\omega_z$. While the vorticity in each slice located in the same filament has the same orientation, we take care to orient the three different filaments (see slices) such that $\omega_z$ is mostly positive (negative) in the first and third (second and fourth) quadrants prior to stacking.].

In the right-hand panel of vorticity, we address the fraction of the vorticity aligned with the filament axis. For each slice, we compute the mass-weighted average of both $\omega_{z}^2$ and $\omega_{\rm tot}^2=(\omega_x^2+\omega_y^2+\omega_z^2)$ at each radius $r$, 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 was randomly oriented, this ratio would be $\sim 0.5$, which is the average of ${\rm cos}(\theta)$ uniformly distributed between $0$ and $1$. Instead, the typical ratio is $\gsim 0.6$ within the shock and increases slightly towards smaller radii, displaying a preference for the vorticity to be aligned with the filament axis. This is broadly consistent with the results of <cit.> who found a preference for the filament vorticity to be aligned with its axis, with an excess probability of $\sim 20\%$ to have an angle less than $60^{\circ}$, corresponding to a cosine greater than $0.5$.

§.§.§ Ideal Vortex Model

Our ideal vortex model, used to interpret the main features of the radial inertial forces from summation method another. The orange dashed circle represents the “shock” radius, $R$, and we have added circles at $0.6R$ and $0.25R$, schematically representing the borders between the three zones in our simulated filaments. Each quadrant has a rotational vortex centred at $(x,y)=(\pm a,\pm a)\simeq(\pm 0.42\, R, \pm 0.42\, R)$. 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 $r'\,^{0.5}$, where $r'$ is the distance to the vortex centre.

While the gas dynamics in filaments is inherently three-dimensional and complex, the quadrupolar vortex structure seen in streamplot-vorticity appears to be the most robust feature common to all filament slices. Therefore, we here examine to what extent such a velocity field alone can reproduce the robust features of the radial inertial forces seen in summation method another (centre panel). To this end, we construct a simple two-dimensional toy model for the filament velocity field, characterized 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 streamplot-vorticity. This model is shown in vortex velocity. 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. $\alpha$ sets the slope of the rotation curve with respect to the vortex centres while $\beta$ is a normalization constant. $r'$ is the distance to the centre of each vortex, which is assumed to be at $(x,y)=(\pm a,\pm a)$.

The vorticity corresponding to such a velocity field, in the vortex frame, is:
ω⃗=∇×v⃗=1/r'(r'v_ϕ')r'ẑ'=±β(α+1) r'^α-1ẑ'.
For any $\alpha\not=-1$, the vorticity is non-zero[$\alpha=-1$ corresponds to a constant specific angular momentum profile.] and finite at all $r'>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 $r'$ and $\phi'$ to the unprimed coordinates $x$ and $y$ depends on the quadrant. For example, in the first quadrant $r'=[(x-a)^2 + (y-a)^2]^{1/2}$ and $\phi'=\arctan[(y-a)/(x-a)]$. We can then write out the full velocity field in terms of the polar coordinates in the filament frame, $r$ and $\phi$ using $x=r\cos(\phi)$ and $y=r\sin(\phi)$. We use this model within the boundaries $|x|,|y|\le R$, where $R$ can be thought of as analogous to the shock radius of the filament, $r_{\rm shock}$. We hereafter place the vortices at a radius of $r=0.6\,R$, motivated by streamplot, so $a/R=0.6/ \sqrt{2}\simeq 0.42$. The trends at $r<a$ and $r>a$ are independent of this choice.

A visual inspection of streamplot indicates that, the rotation velocity increases away from the vortex centres. We, therefore, take $\alpha>0$. We further require that the ram/turbulent pressure term, $-v_r\partial v_r/\partial r$ be negative, based on the middle panel of summation method another, which turns out to require $\alpha<1$. We use $\alpha=0.5$ as our fiducial value, but note that we obtain similar trends with any value of $\alpha\sim(0.1-0.9)$. The value of the normalization constant $\beta$ is arbitrary and we hereafter take $\beta=1$.

Given parameters, $a$, $\beta$, and $\alpha$, we compute the velocity field and the resulting fictitious forces using the same method we used in summation. 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 two-dimensional, all axial derivatives vanish, i.e., $\partial/\partial z=0$. Furthermore, in our highly idealized model, the mean azimuthal and radial velocities are both zero as well, $v_{\phi,{\rm mean}}=v_{r,{\rm mean}}=0$.
This leaves only the three components in the radial inertial term which depend only on the residual velocities, namely the turbulent pressure (Mturb_r), residual azimuthal shear (Mturb_p), and residual centrifugal force (decomposition of the rotation term) (red, orange, and purple solid lines in the centre panel of summation method another). We evaluate these following our methodology described in Appendix sim_numerics.
Finally, we normalize all forces by $(r/R)^{\gamma}$, with $\gamma=1.48$. This accounts for the radial dependence of the gravitational force, $F_{\rm g}$, the denominator in all the profiles shown in summation method another which is given by the total mass interior to radius $r$ (blue line in line mass).

Inertial (fictitious) force terms resulting from our toy model of an ideal quadrupolar vortex structure (vortex velocity). Since our model is completely symmetric, the mean parts of $v_r$ and $v_\phi$ are by definition $0$, and therefore only the residual parts contribute. Each force term from our model has been scaled by an additional $r^{\gamma}$ to account for the radial dependence of $F_{\rm g}$, the denominator in the ratios presented in summation method another. Our simple model reproduces many of the trends seen in the simulation in the V zone, at $0.25\lsim r/r_{\rm shock}\lsim 0.6$ (summation method another, middle panel, solid lines), where the ram pressure force has decreased and vortices dominate the filament dynamics (see text for details).

The resulting force terms are plotted in ideal vortex results. These reproduce many of the trends seen in the corresponding curves in the centre panel of summation method another. At small radii ($r\lsim 0.5 R$), the residual rotation is mostly balanced by the azimuthal shear, both of which appear to diverge as $r\rightarrow 0$ as $r^{-\gamma}$. The turbulent pressure contributes to the balance of the rotation term as well, but saturates at a finite value as $r\rightarrow 0$. Even the sign-flip in this term at $r\lsim 0.6\,R$, 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 summation method another near the location of the vortices at $r\sim 0.6\,r_{\rm shock}$. The corresponding sign flip of the shear force at $r=a$ is not seen in the data. However, in the simulations, the inertial forces in the T zone at $r>0.6\,r_{\rm shock}$ are dominated by ram pressure from inflowing gas, not included in our ideal vortex model, which is therefore only expected to represent the data at $r<a=0.6\,R$. Furthermore, unlike our simple two-dimensional toy model, the actual filament dynamics are complex and three-dimensional, with vorticity that is not even truly aligned with the filament axis (vorticity). Despite these caveats, our model confirms that the quadrupolar vorticity structure in intergalactic filaments does yield fictitious force terms that resemble the inertial terms dominating the filament dynamics at $r\lsim 0.6\,r_{\rm shock}$. More detailed theoretical studies of cosmic velocity fields in filaments can be found in, e.g., <cit.>.

§ DISCUSSION

§.§ Comparison to Previous work

§.§.§ Filament Thermal Structure and Size

We compare our results to those of <cit.>, who used an AMR simulation to study the properties of an intergalactic filament feeding the progenitor of a MW-mass halo ($\Mv\sim 10^{11.5}\msun$ at $z=0$). While they examined the filament in the redshift interval $z\sim (3.5-8)$, their main analysis focused on $z\sim 4$ as ours does, 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, self-gravitating, isothermal cylinder, out to the radius where the sheet begins dominating the mass distribution. This was true for both the gas and the 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, $r_0$, for the gas and dark matter were within a few tens of percent of each other, while these differed by a factor of $\gsim 3$ in our simulations. Likewise, the dark matter effective temperature in <cit.> was only a factor of $\lsim 2$ larger than the central gas temperature (see their figure 7), while it was $\sim 30$ times larger in our simulations (DM_prop).
Furthermore, the shock surrounding the filament in their simulation was both weaker and narrower than in our own. They measured an azimuthally-averaged gas-temperature increase of $\lsim 60\%$ between the filament core and the shock,
while the temperature in our simulation increases by a factor of $\sim (20-30)$. Likewise, they found the shock to be nearly isothermal, with the gas temperature reaching its core value $\sim 3\kpc$ from the shock front, while in our case the gas remains hot outside the cooling radius of $\sim 0.5\,r_{\rm shock}\sim 25\kpc$ (tcool) and only reached the core temperature $\sim (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 only be true 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 <cit.>.