© The Authors 2025

1. Introduction

The hierarchical formation of large-scale structures in the Universe is driven by gravity, which competes with cosmic expansion (e.g. Peebles 2020). This leads to the anisotropic collapse of primordial density field fluctuations (Zel’Dovich 1970). Dark matter (DM) and baryons thus flow from underdense (voids) to overdense environments (galaxy clusters) that are connected by cosmic filaments. Together, they form the cosmic web (Bond et al. 1996). At the nodes of this network, galaxy clusters are the most massive gravitationally bound structures in the Universe; they are thus used to constrain, in particular, the cosmic matter density and distribution parameters, that is, Ω_m and σ₈, of the ΛCDM model. To constrain these parameters via the cluster number counts (e.g. Planck Collaboration XX 2014; Salvati et al. 2018; Aymerich et al. 2024) or the baryon fraction (e.g. Wicker et al. 2023), however, we need to estimate and calibrate their mass precisely and in an unbiased manner (see e.g. Nagai et al. 2007, Lau et al. 2009, 2013, Biffi et al. 2016, Salvati et al. 2019, Gianfagna et al. 2021, Lebeau et al. 2024a, Aymerich et al. 2024, Braspenning et al. 2025). In particular, when the mass of a cluster is dereived from observations of its intracluster medium (ICM) in the X-ray or sub-millimetre wavelengths, we usually assume hydrostatic equilibrium (see Ettori et al. 2013 or Kravtsov & Borgani 2012 for reviews). It was shown in several studies (observations and simulations; see e.g. Salvati et al. 2018, Gianfagna et al. 2021, Lebeau et al. 2024a and references therein), however, that this hypothesis tends to be too strong. Moreover, through the choices for the data analysis, for instance, accounting for substructures (e.g. Meneghetti et al. 2010) or instrumental effects (e.g. Schellenberger et al. 2015), the cluster mass estimation is biased. This is the so-called hydrostatic mass bias. On the side of physical sources of bias, it was shown that turbulence in the ICM might add to the non-thermal pressure support by a few percent to 30% (e.g., Rasia et al. 2004; Nelson et al. 2014; Vazza et al. 2018; Pearce et al. 2020; Angelinelli et al. 2020). When this is accounted for in the mass estimation, it might contribute to alleviating the hydrostatic mass bias.

In the context of precision cosmology, it has therefore become crucial to quantify turbulence in the ICM and understand how it develops. The increasing size and resolution of cosmological simulations in the past decade enabled an important effort that was dedicated to constraining turbulence in the ICM from simulations (e.g. Norman & Bryan 1999; Dolag et al. 2005; Nagai et al. 2007; Gaspari & Churazov 2013; Miniati 2014; Porter et al. 2015; Vazza et al. 2017; Vallés-Pérez et al. 2021a; Ayromlou et al. 2024; Groth et al. 2025). It was shown that vorticity, and then turbulence, is mainly generated by shocks as a result of mergers (e.g. ZuHone 2011; Vazza et al. 2012; Nagai et al. 2013), by accretion of cosmic gas (e.g. Iapichino et al. 2011), and by the feedback of active galactic nuclei (AGN) in cluster cores (e.g. Vazza et al. 2013; Bourne & Sijacki 2017), involving multiple injection scales. In addition, the ICM is multiphase (e.g. Wang et al. 2021; Mohapatra et al. 2022) and stratified (e.g. Shi et al. 2018; Mohapatra et al. 2020; Simonte et al. 2022; Wang et al. 2023), and other mechanisms such as galaxy motion (e.g. Ruszkowski & Oh 2011) and cosmic-ray pressure (e.g. Beresnyak et al. 2009) might contribute to turbulence. The combination of these processes leads to a more complex turbulence than in Kolmogorov’s idealised case (Kolmogorov 1941), which means that it is challenging to model.

At the same time, studies of the filament populations based on simulations (e.g. Cautun et al. 2014; Gheller & Vazza 2019; Pereyra et al. 2020; Kuchner et al. 2020; Galárraga-Espinosa et al. 2021; Rost et al. 2021; Gouin et al. 2022), including their detectability through HI 21 cm emission (Kooistra et al. 2017, 2019), and observations (e.g. Tanimura et al. 2022; Gouin et al. 2023; Gallo et al. 2024) have shown the impact of cosmic filaments on cluster shapes and dynamical states (e.g. Gouin et al. 2020, 2021; Santoni et al. 2024; Capalbo et al. 2025) and how cosmic gas is fed to clusters through filaments (Vallés-Pérez et al. 2020; Vurm et al. 2023; Kuchner et al. 2022). Some works studied shocks and heating of gas in filaments (e.g. Ryu et al. 2003, 2008; Pfrommer et al. 2006; Bennett & Sijacki 2020), the amplification of magnetic fields in filaments (Vazza et al. 2014), which is tightly related to turbulence and also velocity fields (e.g. Zhu et al. 2010, 2013) and the generation of turbulence at the interface of the intergalactic medium and galaxy haloes (e.g. Cornuault et al. 2018). To our knowledge, however, none studied deeply how vorticities develop along the way from the filaments that are connected to cluster to the ICM, however. The aim of this work is thus to conduct a qualitative and quantitative study of the gas dynamics in filaments that are attached to a specific cluster, Virgo.

To this end, we examined the structure of the velocity field and temperature along large filaments that are connected to a state-of-the-art high-resolution simulated replica of the Virgo cluster in its local environment (Sorce et al. 2021) at redshift z = 0. In the following, we introduce the simulation and the method we adopted to construct slices of temperature and velocity, and the derived quantities in Sect. 2. We qualitatively describe the gas flow and the turbulence sources using longitudinal and transverse cuts of the filaments connected to the Virgo cluster in Sects. 3 and 4. In Sect. 5 we decompose the velocity field into compressive and solenoidal components and compute their 2D power spectra and the spectrum of the total velocity field in the transverse slices. We discuss our results in Sect. 6 and conclude in Sect. 7.

Fig. 1. From left to right, we show longitudinal cuts of the temperature, electron density, Mach number, and velocity in the Virgo replica and the two filaments to which it is connected. The maps are about 22 Mpc wide, contain 157282 pixels, and are centred on the Virgo core. In the third panel, the Mach number and colour bands are used to display the wide range of values: The black-to-white colour band is on a linear scale in the [0.1,1] range, and the white-to-blue colour band is on a log scale in the [1100] range. In the right panel, the background map is the velocity along the x-axis of the simulation box, vx, and the foreground arrows represent the norm of the velocity field in the x − y plane, ||v||. The background (foreground) filament is in the left (right) part of each panel. The buoyant bubble is located roughly at a position (x = −5.0; y = 0.5) Mpc. Multiple galaxies are visible in the density map, for example, at positions ( − 2.5; −0.5) and ( − 6.0; 1.0) Mpc.

2. Method

2.1. The Virgo replica simulation

Using constrained initial conditions of the local Universe, reconstructed from the position and peculiar velocities of galaxies in the Cosmicflows-2 (Tully et al. 2013) dataset (see Sorce et al. 2016 for details), Sorce et al. (2019) simulated DM replicas of the Virgo cluster within its local cosmic environment. The most representative of the 200 DM simulations in terms of average properties and merging history compared to the full sample was used to run a high-resolution zoom-in hydrodynamical simulation of Virgo (Sorce et al. 2021). The DM and high-resolution hydrodynamical simulations agree well with observations. In particular, they accurately reproduce the observed filament in the foreground of Virgo along our observation sightline and also the group of galaxies that falls onto Virgo (more details about the simulation and its agreement with observations can be found in Sorce et al. 2021 and Lebeau et al. 2024a).

This simulation used the Planck Collaboration XVI (2014) cosmological parameters; the Hubble constant, H₀ = 67.77 km s⁻¹ Mpc⁻¹, dark energy density, Ω_Λ = 0.693, total matter density, Ω_m = 0.307, baryonic density, Ω_b = 0.048, amplitude of the matter power spectrum at 8 Mpc h⁻¹, σ₈ = 0.829, and spectral index, n_s = 0.961. It was produced using the adaptive mesh refinement (AMR) hydrodynamical code RAMSES (Teyssier 2002). In this run, Euler equations are solved following the MUSCL-Hancock method using the Riemann solver of Toro et al. (1994), which is second-order accurate in space, with a variable time step based on the Courant Friedrich Levy condition (see Teyssier 2002, for more details). Within the 500 Mpc h⁻¹ local Universe box, the zoom region is a 30 Mpc diameter sphere with an effective resolution, i.e. only in the zoom region, of 8192³ DM particles of mass m_DM = 3 × 10⁷ M_⊙. The finest cell size of the AMR grid is 0.35 kpc following a pseudo-Lagrangian refinement criterion, i.e. based on DM particles density (see Sorce et al. 2021, for more details). In the high-resolution run, the maximum resolution of 0.35 kpc is reached only in the core of the cluster; in the filaments, the resolution of the cells is mostly 22.5 and 45 kpc. In the following, we used a best resolution of 22.5 kpc to compute the Hodge-Helmholtz decomposition and the power spectra. In addition, it is worth mentioning that the development of turbulent motion at small scales in this simulation is limited by numerical diffusion due to the accuracy order of the Riemann solver and the maximum resolution of the AMR grid. Therefore, the results presented in this study are to be taken as a lower limit of the actual intensity of turbulent motion in the ICM that could be better resolved with a higher-order Riemann solver and a finer grid resolution.

In addition, the simulation implements sub-grid models for radiative gas cooling and heating, star formation, and kinetic feedback from the AGN and type II supernova (SN) similar to the Horizon-AGN implementation of Dubois et al. (2014, 2016) with an AGN feedback model improved by orienting the jet according to the black hole (BH) spin (see Dubois et al. 2021 for details). The simulation was run from redshift z = 120 to z = 0. In this work, we only used the snapshot at z = 0.

2.2. Longitudinal and transverse cuts along the Virgo filaments

We used the rdramses code¹, to extract the DM particles and gas cell properties from the simulation. To identify the Virgo DM halo and its galaxies, we used the Tweed et al. (2009) halo finder, respectively, on the DM particles and the star particles. In this work, the gas is defined as the baryonic component in the simulation cells, in which hydrodynamics equations are solved following an Eulerian approach; it thus does not include the stars. We extracted very thin slices from the simulation to study variations of temperature and velocity at small scales in the filaments connected to Virgo. We adapted the code used in Lebeau et al. (2024a) to do so. We set the integration depth along a given sightline to 22.5 kpc, which is 64 times the finest cell size of 0.35 kpc. In underdense regions, cells can be larger than 22.5 kpc but still partly belong to the selected slice; in that case, only their fraction inside the slice was accounted for.

The Virgo replica is connected to two cosmic filaments aligned along the same direction, almost perfectly in the xy simulation box plane. We thus extracted a slice in this plane centred on the Virgo cluster core; we call it the longitudinal cut throughout the paper; it is presented in Fig. 1. Moreover, we show a zoom into the filament with Δx_cen = x − x_{virgo centre} < 0 (Δx_cen > 0) in Fig. 2 (Fig. 3), x_{virgo centre} being the centre of the Virgo DM halo. Given that this simulation replicates the Virgo cluster in its local environment, including the Milky Way, we can put ourselves in the sightline of the observer, as in Lebeau et al. (2024a). Thus, we can call them background filament (left side of panels of Fig. 1) and foreground filament (right side of panels of Fig. 1).

Fig. 2. Zoom into the background filament in the longitudinal cut. The maps are 11.061 × 12 Mpc large with the same resolution as those presented in Fig. 1. The left panel shows the velocity field, similar to the right panel of Fig. 1, the central panel shows the divergence, and the right panel shows the z component of the vorticity. The Virgo R500 and virial radius, Rvir, are shown as black circle arcs on each map. The merger-accelerated accretion shock is located roughly in the range x = [ − 5, −2] Mpc and y = [ − 5, −1] Mpc.

Fig. 3. Zoom into the foreground filament in the longitudinal cut. The maps are 11.061 × 12 Mpc large with the same resolution as those presented in Fig. 1. The left panel shows the velocity field, similar to the right panel of Fig. 1, the central panel shows the divergence, and the right panel shows the z component of the vorticity. The Virgo R500 and virial radius, Rvir, are shown as black circle arcs on each map.

Then, we extracted perpendicular slices to the filaments to study the temperature and velocity fields in their core. We extracted slices in planes parallel to the xz simulation box plane instead of computing slices with a normal vector aligned with the major axis of the filaments (the angle between the two is roughly 25 degrees) to avoid non-physical small-scale blurring due to the misalignment with the AMR grid. This choice preserves the resolution at very small scales, allowing a qualitative analysis of the structure of both temperature and the velocity fields in the filaments. These slices are presented in Figs. 4 and 5 for the background and foreground filaments, respectively. To visualise the position of the slices, we refer the reader to Fig. A.1 in Appendix A.

Fig. 4. Transverse cuts along the background filament. From left to right, we present cuts at distances of 8, 6, 4 and 2 Mpc from the Virgo centre. From top to bottom, we present the temperature, the velocity field with the background map being the velocity along the z simulation box axis, vz, and the foreground arrows representing the norm of the velocity field in the zy plane, ||v||, the divergence of the velocity field and the x component of its vorticity. Similarly to longitudinal cuts, the maps are 22.122 Mpc wide, contain 15 7282 pixels, and are centred on the Virgo centre.

Fig. 5. Transverse cuts along the foreground filament. From left to right, we present cuts at distances of 2, 4 and 6 Mpc from the Virgo centre. From top to bottom, we present the temperature, the velocity field with the background map being the velocity along the z simulation box axis, vz, and the foreground arrows representing the norm of the velocity field in the zy plane, ||v||, the divergence of the velocity field and the x component of its vorticity. Similarly to longitudinal cuts, the maps are 22.122 Mpc wide, contain 157282 pixels, and are centred on the Virgo centre.

The slices presented in Figs. 1, 4 and 5 are (22.122 Mpc)² and contain 15728² pixels; each pixel is thus 1.4 kpc wide. The pixel size is the same for the zooms presented in Figs. 2 and 3. For each of the longitudinal and transverse cuts, we present the mass-weighted mean temperature and electron density (the transverse cuts can be found in Appendix B), and the mean velocity of the cells integrated along the sightline in the Virgo restframe. For the velocity, the blue-to-red field in the background is the velocity along the abscissa of the figure, that is, v_x (v_z) for the longitudinal (transverse) cut in Figs. 1, 2 and 3 (Figs. 4 and 5), and the black-to-white arrows in the foreground represent the velocity field in the plane of the slice. In addition, we present the 2D divergence and the z (x) component of the curl, that is, the vorticity of the velocity field in the longitudinal (transverse) cuts in the central and right panels of Figs. 2 and 3 (third and bottom panels of Figs. 4 and 5). For the longitudinal cut, the divergence in the xy plane and the z component of the vorticity are defined as

$$\begin{array}{c}\hfill \mathbf{\nabla }·\mathit{v}=\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}\text{and}{\omega }_z=\frac{\partial v_y}{\partial x}-\frac{\partial v_x}{\partial y}·\end{array}$$(1)

For the transverse cut, the divergence in the yz plane and the x component of the vorticity are defined as

$$\begin{array}{c}\hfill \mathbf{\nabla }·\mathit{v}=\frac{\partial v_y}{\partial y}+\frac{\partial v_z}{\partial z}\text{and}{\omega }_x=\frac{\partial v_z}{\partial y}-\frac{\partial v_y}{\partial z}·\end{array}$$(2)

In Sect. 5, we decompose the velocity field in compressive and solenoidal components following the Hodge-Helmholtz procedure as in Vazza et al. (2017). This can only be done on a regular 3D grid, however, whereas our simulation of Virgo, based on the AMR code RAMSES, has an irregular grid. We thus adapted the simulation output data specifically for this section. Given the range of the turbulent scales considered in this study, roughly from a few Mpc to a few tens of kpc, we did not need to use the simulation at its finest resolution to study statistical properties of the velocity field through the power spectrum; we thus set the maximum resolution level to cells of 22.5 kpc while reducing the output data with rdramses. It means that we did not use smaller cells than this resolution limit in the AMR dependence tree; the outputted physical properties are thus the mean of all the smaller so-called children cells, as defined in Teyssier (2002), associated with a given larger so-called parent cell. It is worth noting that this cell size equals the thickness of the slices considered in this work. Cells with a lower resolution, i.e. larger cells in underdense regions, were divided into smaller ones up to the resolution limit with the same physical properties except for the mass equally distributed among all the new smaller cells to keep the same density. Choosing this resolution limit also avoided generating a large amount of data, given that the number of cells is multiplied by eight each time we increase the resolution level by one. It also significantly reduced the calculation time.

3. Matter flows across the filaments

In this section, we focus on the longitudinal cut to study the gas across the filament. First, in the left panel of Fig. 1, the temperature reveals the components of the Virgo local environment. In the left part of the slice, a filament of almost 5 Mpc diameter extends over more than 10 Mpc from the limit of the simulation box to very deep in Virgo. The temperature in the core of the filament is in the range [10^6.5, 10⁷] K, with an electron density below 10⁻⁴ cm⁻³, it is thus a Warm-Hot Intergalactic Medium (WHIM) given the classification in the literature (e.g. Martizzi et al. 2019; Gouin et al. 2022). At the borders of the filament, the temperature can reach up to 10^7.5 K due to the filament accretion shocks, as discussed below. In the right part of the slice, we observe a smaller, hotter and less collimated filament connected to Virgo. Then, underdense regions surround Virgo and the filaments, with temperatures below 10⁶ K. Finally, the temperature in the Virgo ICM ranges from 10⁷ K to more than 10⁸ K, as expected for galaxy clusters.

This thin slice of 22.5 kpc width reveals small-scale temperature fluctuations in the Virg filaments and core. In the background filament, numerous wavy structures are visible; these might be the signature of instabilities leading to rotation and potentially turbulence. They might be generated by Kelvin-Helmholtz instabilities due to streams with different relative velocities and by galaxies acting as obstacles generating trailing turbulence. In the Virgo core, temperature fluctuations are even at smaller scales and show a very complex gas stirring, which traces complex turbulence.

In the second panel, we present the electron density. The densest region is the Virgo ICM, as expected due to the gravitational potential well. Then, in the filaments, by zooming in on the density map, we observe elongated density peaks with a size of about a hundred kpc, coinciding with cold regions in the temperature maps. These cold density peaks are, in fact, galaxies, given that we could match their positions with the positions of galaxies detected by the halo finder. For example, two galaxies are positioned roughly at the coordinates (x = −2.5; y = −0.5) and ( − 6.0; 1.0) Mpc. On the other hand, we observe an underdense and hot region of similar scale in the background roughly at coordinates ( − 5.0; 0.5) Mpc; this buoyant bubble might be due to SN or AGN feedback.

Then, in the central panel in Fig. 1, we present the dimensionless flow Mach number, M = v/c_s, with v the magnitude of the velocity field and c_s the sound speed defined as

$$\begin{array}{c}\hfill c_{\text{s}}=\sqrt{\frac{\gamma k_{\text{B}}T}{\mathrm{\mu }{\mathrm{m}}_{\text{p}}}},\end{array}$$(3)

with k_B the Boltzmann constant, μ the mean molecular weight, m_p the proton mass and γ the adiabatic index set to 5/3 in the simulation, given that we assume for simplicity a monoatomic perfect gas. We use two colour bands to display the large range of values: the black-to-white is on a linear scale in the [0.1,1] range, and the white-to-blue is on a log scale in the [1100] range.

Similarly to previous works (e.g. Gaspari & Churazov 2013; Vazza et al. 2017), we can see that the gas flow is subsonic to transonic in the ICM with a flow Mach number below M = 2. Then, in the filaments, the gas flow is supersonic, between M = 2 and M = 20, which is due both to the temperature and density drop and to the higher overall velocity in the filament, as we can see in the right panel of Fig. 1. In the upper underdense region, where the temperature is lower, but the density is not much lower than in the filaments, we found flow Mach numbers similar to the latter. Finally, in the lower underdense region, where the density and temperature are much lower than in the upper one, the flow Mach number is of the order of M = 100, which is consistent with Ryu et al. (2003). This is because the low temperature in this medium induces a small sound speed, thus increasing the flow Mach number for a velocity similar to that in the upper underdense region; it is thus a highly supersonic flow.

Finally, we focus on the velocity field in the slice presented in the right panel of Fig. 1. We observe that most of the cosmic gas is first accreted from underdense regions to the cosmic matter sheet, particularly around the background filament (see Sect. 3), then to filaments and finally to Virgo. It is accelerated up to 1000 km s⁻¹ in the underdense regions and matter sheet before being slowed down and heated when entering the filament and then re-accelerated up to 1200 km s⁻¹ before being strongly shocked and heated again, at the termination shock as called by, e.g., Vurm et al. (2023) when entering the cluster. We see how filaments funnel the gas to Virgo and act as cosmic highways, as shown in, e.g., Gouin et al. (2022) and Vurm et al. (2023). Virgo is thus the accretion endpoint where the two major flows from the background and foreground filaments meet and stir. Indeed, the x component of the velocity field, v_x, is entirely positive (red) in the background filament. In contrast, it is entirely negative (blue) in the foreground filament, and there is a mix of pale red and blue in Virgo. Finally, a strong outflow is located below Virgo on this slice, located roughly in the range x = [ − 5, −2] Mpc and y = [ − 5, −1] Mpc, which seems to generate vorticity at its boundaries; see left panel of Fig. 2 for more details. This outflow has already been identified in Lebeau et al. (2024b) and might result from the recent major merger identified in the simulation by Sorce et al. (2021). It ressembles a merger-accelerated accretion shock as studied in Zhang et al. (2020), which could generate vorticity via the baroclinic mecanism (see e.g. Vazza et al. 2017; Wittor et al. 2017).

To highlight the velocity field structure in the filaments, particularly shocks and rotation, we show zooms of the background and foreground filaments of the velocity field, their divergence and the vorticity, respectively, in Figs. 2 and 3. The Mach number and density-weighted divergence and vorticity in these zoom regions can also be found in Appendix C. On each map, the Virgo R₅₀₀ = 1.09 Mpc and virial radius, R_vir = 2.15 Mpc, are shown as black circle arcs. For the background filament, in Fig. 2, we observe that the divergence (central panel) highlights the filament boundary. Before entering the filament, the gas flow is essentially laminar with a velocity between 800 and 1000 km s⁻¹ in a single coherent flow, and the divergence is null. When entering the filament, it experiences an oblique shock, or filament accretion shock, as shown by the divergence, and is thus slowed down and heated, as shown in Fig. 1. Moreover, the gas starts swirling, as shown by the vorticity (right panel). Then, in the filament, the gas flow is accelerated by its collapse on Virgo, from less than 200 km s⁻¹ at Δx_cen = −10 Mpc to more than 1200 km s⁻¹ at Δx_cen = −2 Mpc. It is a more complex flow with multiple parallel streams, as shown by the divergence that transports a mildly swirling gas with relatively large eddies from Δx_cen = −10 Mpc to Δx_cen = −2 Mpc. Finally, the gas enters the cluster; there is an orthogonal shock, named termination shock in, e.g., Vurm et al. (2023), at Δx_cen = −2 Mpc, which is of about the virial radius, as shown by the high divergence at this distance. The gas is once again slowed down and heated. We also observe that in the Virgo core, the gas highly rotates at much smaller scales, partly due to this termination shock.

The picture is even more complex for the foreground filament, presented in Fig. 3. We observe a gas flow, but the filament from the bottom right corner of the map to the Virgo core is less distinctly defined. It is visible on the divergence map as this tracer shows multistreams inside the filament, but its boundaries are not as clear as for the background filament. It is confirmed by the vorticity, where we observe slow and large-scale rotation in the same region, and by the temperature in the left panel of Fig. 1. Nonetheless, we observe the same features as for the background filament: the gas is shocked twice, first when entering the filament and then in the cluster close to the virial radius. This generates eddies, first quite slow and at large scales in the filament, then much quicker and at small scales in the cluster.

From this longitudinal cut analysis, the picture that emerges is that Virgo is fed by two aligned filaments with different structurations, the collimated background filament embedded in a cosmic sheet and the funnel-shaped foreground one, but comparable gas flows, both contributing to the complex gas stirring in the Virgo core.

4. Matter flows from the cosmic matter sheet to the filaments

We now study the gas flowing from underdense regions to the filaments, including through a cosmic matter sheet in between for the background filament. We use transverse cuts across at different distances from the cluster centre along the filaments.

4.1. Background filament

In Fig. 4, we show four transverse cuts across the background filament at distances Δx_cen = { − 8, −6, −4, −2} Mpc from the cluster centre. The top row presents the temperature, the second the velocity field, the third the divergence, and the bottom row the x component of the vorticity, ω_x. The joint study of these four tracers helps us understand the structuration of this filament and how gas is funnelled towards Virgo. When focusing on the Δx_cen = −8 Mpc slice, we can see that the filament is embedded in a cosmic matter sheet, as we have already shown in Lebeau et al. (2024b). The gas goes from underdense regions to the matter sheet. It is then accreted by the filament, its typical baryonic radius, which we could define by the maximum of divergence of the velocity field, being only about 2 Mpc at this distance. We can already observe multistreams, heating and stirring in the matter sheet, that divergence separates from underdense regions, and one large eddy around the spine of the filament.

Then, in the cut across at Δx_cen = { − 8, −6, −4} Mpc, we can see that the two opposite gas flows fueling the filament are not perfectly aligned, which generate gas coiling around the spine of the filament and thus a large eddy, it is particularly visible on the Δx_cen = −4 Mpc slice (third column of Fig. 4). We also observe that when approaching Virgo, the gas in the cosmic sheet and the filament heats up and accelerates. Moreover, if defined with the velocity maximum of divergence, its envelope goes from about 3–4 Mpc, while vortices gain intensity. Finally, when the gas enters the cluster (right column of

Fig. 4), as we have already shown with the longitudinal cut, it is heated and stirred, leading to a large amount of turbulence. We also observe that the gas reaches its highest velocity at Δx_cen = −2 Mpc and that the divergence defines the cluster baryonic boundary fairly well, which we could relate to the accretion radius. This radius has about the same value as the R_vir but is much smaller than the splashback radius estimated from baryons in Lebeau et al. (2024b). If observable, the divergence of the velocity field might thus be a boundary definition of clusters that could be comparable with R_vir.

4.2. Foreground filament

In Fig. 5, we show three transverse cuts across the foreground filament at distances Δx_cen = {2, 4, 6} Mpc from the cluster centre. At larger distances, the foreground filament is much more diluted than the background one. The figure is organised similarly to Fig. 4. Coherently with what we observed on the longitudinal cut (see Fig. 3), we can see that this filament is much more diffuse than the background one. Indeed, at Δx_cen = 6 Mpc, the filament can be identified as the converging point of the velocity field (second row) and the hottest area, but it does not stand out much from its environment. Moreover, there are multiple flows towards it, whereas the background filament is the confluence of the two gas flows from each part of the matter sheet. The gas around this filament directly flows from underdense regions, and we cannot distinguish a cosmic matter sheet in this area. At Δx_cen = 4 Mpc, multiple gas streams are also visible. Still, the filament is more defined as it is hotter, better defined by the divergence, and the gas inside is rotating, as shown by the vorticity. It has no cylindrical shape, however, which is certainly due to the multiple streams that might be considered as smaller-scale filaments (e.g. Aragón-Calvo et al. 2010; Feldbrugge & van de Weygaert 2025). Then, at Δx_cen = 2 Mpc, we are roughly at the Virgo virial radius, and the picture is indeed roughly similar to the cut at Δx_cen = −2 Mpc in the background filament.

In conclusion, we can see that the gas swirling starts in this foreground filament, similarly to the background filament, even though it is not as large and collimated. The eddies seem to be transported along the gas flow towards Virgo while becoming more intense and numerous due to shocks and gas accretion at the boundaries of the filament.

5. Velocity fields in filaments

In the previous sections, we presented a qualitative study of the gas flows in the filaments connected to Virgo, showing how eddies are generated and grow towards Virgo. We now decompose the velocity field into its compressive and solenoidal components to quantify their contribution to the total velocity field and characterise turbulence through the 2D total power spectrum.

5.1. Hodge-Helmholtz velocity field decomposition

The velocity field is a mixture of compressive and rotating flows. The Hodge-Helmholtz theorem stipulates that any vector can be separated into a curl-free, a divergence-free, and a harmonic component (e.g. Arfken et al. 2011). In the case of fluid mechanics, the first is the compressive component, the second is the solenoidal component, and we consider the third one to be null in this work, following Vallés-Pérez et al. (2021b). Following Vazza et al. (2017), we apply this decomposition to the Fourier transform of the velocity field written as v_k = v_k, c + v_k, s with

$$\begin{array}{c}\hfill {\mathit{v}}_{k,c}=\mathit{k}(\mathit{k}·{\mathit{v}}_k)/{|k|}^2\text{and}{\mathit{v}}_{k,s}=-\mathit{k}\times (\mathit{k}\times {\mathit{v}}_k)/{|k|}^2.\end{array}$$(4)

We can thus efficiently calculate v_k, c by projecting v_k on k using Fast Fourier Transform (FFT)², then return to real space and deduce the solenoidal component as v_s = v − v_c.

The decomposition was applied on a rectangular 12 × 12 × 20 Mpc box encompassing Virgo and the two filaments. We present slices of this box in Fig. 6. The slices are similar to those presented in Figs. 4 and 5, that is Δx_cen = { − 8, −6, −4, −2, 0, 2, 4, 6} Mpc. The top panels show the total velocity field norm, $||\mathit{v}||=\sqrt{v_x^2+v_y^2+v_z^2}$, the middle panels show the compressive component norm, ||v_c||, and the bottom panels show the solenoidal component norm, ||v_s||.

In every slice and particularly for central slices at Δx_cen = { − 2, 0, 2} Mpc, we can see that the compressive component (middle row) highlights the boundaries of the baryonic matter sheet, the filaments and the cluster and more generally regions where accretion shocks occur. These features traced by the compressive component are coherent with the ones traced with the divergence of the velocity field presented in the central panels of Figs. 2 and 3, and third rows of Figs. 4 and 5, which is expected given that the compressive component is, by definition, curl-free. On the other hand, we note that the solenoidal component traces rotation occurring mostly inside large-scale structures. This feature agrees with that traced by the vorticity presented in the right panels of Figs. 2 and 3 and bottom rows of Figs. 4 and 5, which is also expected given that this component is, by definition, divergence-free. To quantify turbulence in these slices, we analyse their power spectrum in the following subsection.

Fig. 6. Slices of the norm of the 3D velocity field. The top panel presents the total velocity field, and the second (bottom) panel presents the compressive (solenoidal) component obtained via the Hodge-Helmholtz decomposition. The slices are extracted from the 12 × 12 × 20 Mpc box used for the decomposition; the maps are thus 12 × 12 Mpc large, and each pixel is 22.5 kpc large. The slices presented are similar to those in Figs. 4 and 5, that is at Δxcen = { − 8, −6, −4, −2, 0, 2, 4, 6} Mpc.

Fig. 7. Comparison of the total velocity field power spectra. The top and bottom rows present spectra computed from 5 and 1 Mpc square maps, respectively. The left column presents the total velocity field power spectra in slices in the background filament, that is, from red to black at Δxcen = { − 8, −6, −4, −2} Mpc, and in the Virgo core for comparison. The right column presents the total velocity field power spectra in slices in the foreground filament, that is, from black to blue at Δxcen = {2, 4, 6} Mpc, and also in the Virgo core for comparison. For comparison, we show the slopes predicted by Kolmogorov’s (P(k)∝k−11/3) subsonic turbulence theory and Burgers’ (P(k)∝k−4) supersonic turbulence theory.

5.2. 2D power spectra

We now quantitatively characterise the velocity field in the slices and compare the contribution of each component to the total velocity field throughout the filaments. To do so, we compute the power spectrum, P(k), that we define for the total velocity field as

$$\begin{array}{c}\hfill P(\mathit{k})=\frac{1}{2}{|\stackrel{\mathbf{\sim }}{\mathit{v}}(\mathit{k})|}^2=\frac{1}{2}(|{\stackrel{\sim }{v}}_x{(\mathit{k})|}^2+|{\stackrel{\sim }{v}}_y{(\mathit{k})|}^2+|{\stackrel{\sim }{v}}_z(\mathit{k}){|}^2),\end{array}$$(5)

with $\stackrel{\mathbf{\sim }}{\mathit{v}}(\mathit{k})$ the Fourier transform of the total velocity field. The formula is similar for the compressive and solenoidal components of the velocity field. It is worth noting that we compute the power spectrum in the rest frame of the Virgo centre; thus, we do not subtract the coherent flow, in other words, the average gas velocity in each slice, in the filaments from the total velocity, as done in, e.g. Vazza et al. (2017). Tests have shown that, thanks to the properties of the Fourier Transform, it yields exactly the same results. We thus compare the contribution of each component to the total velocity field and do not use the terms solenoidal turbulence and compressive turbulence as in Vazza et al. (2017). We only talk about turbulence when discussing the total power spectrum.

To investigate the velocity field and its components in the core of the filaments and the cluster, but also in their direct surrounding, we compute the 2D power spectrum on 1 and 5 Mpc square maps centred on the core of the filament. The selected regions are displayed as white squares on the top row of Fig. 6 and their centre can be found in Table 1. Since our study focuses on the gas flows, we defined manually the centre of the selected regions as the point of convergence of the velocity flows, as shown in the second row of Figs. 4 and 5. We thus purposefully did not use any filament finders, which usually rely on the positions of DM particles or galaxies. Using a filament finder to define the centre of the maps as the spine of the filament could give a slightly different position compared to the hand-selected centres of the slices, which could marginally modify the velocity power spectrum of the 1 Mpc large maps, but should not modify the velocity power spectra computed from larger maps. The 2D power spectrum is calculated similarly to the 3D case³ (see Eq. (5)), but only on one plane extracted from the 3D regular grid after the Hodge-Helmoltz decomposition, as if computed on a map but without projecting the components of the velocity field along a sightline as in Lebeau et al. (2025).

For the sake of clarity, we only analyse here the power spectra and their ratios computed from 1 and 5 Mpc square maps; all the 2D spectra, including those computed from 2 Mpc square maps, can be found in Fig. G.1 in Appendix G. We chose a 2D power spectrum on slices over 3D to avoid overlap, and thus possible correlations, between studied regions given that the slices are only 2 Mpc apart; we discuss this choice in Sect. 6. To simplify the analysis and enable comparisons of the spectra, we only show the 2D power spectra of the total velocity field in Fig. 7, the solenoidal-to-total power spectra ratios in Fig. 8 and the compressive-to-solenoidal power spectra ratios in Fig. 9. The same figures for the spectra computed from 2 Mpc square maps can be found in Fig. D.1 in Appendix D. In these three figures, the top and bottom panels present the spectra, or ratios, computed from 5 and 1 Mpc square maps, respectively. The left column presents the total power spectrum, or ratios, in the background filament at distances Δx_cen = { − 8, −6, −4, −2, 0} Mpc from red to black and the right column present those in the foreground filament at distances Δx_cen = {0, 2, 4, 6} Mpc from black to blue.

Fig. 8. Ratio of the power spectrum of the solenoidal component over the total velocity field. The top and bottom rows present ratios of spectra computed from, respectively, 5 and 1 Mpc square maps. The left column presents the ratios in slices in the background filament, that is, from red to black at Δxcen = { − 8, −6, −4, −2} Mpc, and in the Virgo core for comparison. The right column presents the ratios in slices in the foreground filament, from black to blue at Δxcen = {2, 4, 6} Mpc, and the Virgo core for comparison.

Fig. 9. Ratio of the power spectrum of the compressive over the solenoidal component of the velocity field. The top and bottom rows present ratios of spectra computed from, respectively, 5 and 1 Mpc square maps. The left column presents the ratios in slices in the background filament, that is, from red to black at Δxcen = { − 8, −6, −4, −2} Mpc, and in the Virgo core for comparison. The right column presents the ratios in slices in the foreground filament, from black to blue at Δxcen = {2, 4, 6} Mpc, and the Virgo core for comparison.

First, we can visually notice, in Fig. 7, that all the total power spectra have a negative slope, typical of a turbulent cascade (Richardson 1922); this is also the case for each component, with relatively similar slopes compared to their associated total power spectrum (see Fig. G.1 in Appendix G). We can also see that, for both filaments and independently of the map size, although less visible for the 5 Mpc square maps, there is an amplitude increase when approaching Virgo, meaning more intense eddies due to more energy injection by shocks at multiple scales, as shown by the divergence of the velocity field in the third panel of Figs. 4 and 5. In the Virgo core, however, the amplitude is lower than at Δx_cen = { − 2, 2} Mpc, because kinetic energy is dissipated and converted into heat, and the eddies are therefore slowed down when they enter the cluster because of the termination shock at R_vir≈ 2 Mpc.

We now compare the power spectra slopes. For simplicity, we only discuss the total velocity and do not compare the slopes of each component, given that they have similar behaviours. To do so, we fitted the power spectra to a power law P(k) = Ak^α, the best-fit values for the slope α can be found in Table 2. The slopes were fitted for k < 0.095 (L > 66 kpc); which is chosen to be slightly smaller in k than the Nyquist frequency, f_N, which is twice the best cell resolution of 22.5 kpc, to avoid contamination in the fit from numerical flattening close to f_N (see Fig. E.1). Each row represents a given square map size from which the power spectra were computed. We observe that in the range Δx_cen = [ − 8, −4] Mpc for the background filament, whatever the map size, the slope does not vary. It gets steeper at Δx_cen = −2 Mpc, likely due to the termination shock front and the associated small-scale shocks within R_vir transferring more kinetic energy into heat. Then, in the foreground filament, the slope gets smoothly steeper when approaching Virgo and is the steepest at Δx_cen = 2 Mpc, most likely for the same reason as in the background filament. In the Virgo core, the slope is less steep than at Δx_cen = −2 and 2 Mpc, particularly on the spectra computed from the 1 Mpc square map. Visually, however, the slope is steeper at the smallest scales in the core for the power spectrum computed from the 5 Mpc square map (black line in the top left panel of Fig. 7), and a small flattening at the largest scales, which could be the injection scale, but since we do not compute the power spectra on larger scales we cannot confirm this speculation. When restricting the fitting interval to k= [0.006, 0.095] kpc⁻¹, we still observe an increase in slope along the filament, but the steepest slope is in the Virgo core. It is certainly due to the termination shock at R_vir that is covered in the 5 Mpc square map.

We then discuss the contribution of each component to the total velocity field by computing the ratios of power spectra in the slices, from the Virgo core to far into the filaments. In the Virgo core, Δx_cen = 0 Mpc is displayed in black in Figs. 8 and 9, whatever the map size, the solenoidal to total ratio is almost flat and equal to one, meaning that the solenoidal component largely dominates the velocity field and that there is indeed almost no contribution to the overall velocity from the compressive component. Moreover, the compressive to the solenoidal ratio in Fig. 9 is generally the lowest among the slices, whatever the map size. It is in agreement with Vazza et al. (2017) who also found a solenoidal-dominated turbulence in the core of their simulated cluster.

At 2 Mpc from the Virgo core, Δx_cen = −2 and 2 Mpc, that is, at the Virgo connection with both the background and foreground filaments and at about the virial radius, the solenoidal component contains most of the energy in the flow of gas again. Still, the energy of the compressive component is slightly higher than in the Virgo core, as we see in Fig. 9, whatever the map size. In particular, for the 1 Mpc square map in the background filament (bottom left panels of Figs. 8 and 9), the solenoidal to total ratio is the lowest, and the solenoidal component spectrum is almost equal to that of the compressive component at small scales. It might be because the compressive component traces the gravitational accretion around the spine of the filament.

When focusing on the background filament, Δx_cen = { − 8, −6, −4, −2} Mpc, the amplitude ratio between the spectra of the two components depends on the map size, as we see in the left column of Fig. 9. Looking at the spectrum computed from the 1 Mpc square maps (bottom row), we observe that the solenoidal component dominates the compressive one, like in the Virgo core. It is due to the large eddy in the core of the filament (see second row of Fig. 4, particularly the second panel at Δx_cen = −6 Mpc). The spectrum computed from the maps with a size of 5 Mpc shows a different trend, however. The compressive component dominates the solenoidal component at Δx_cen = −8 Mpc, and we observe a smooth progressive inversion as the ratio goes from between 1 and 5 at Δx_cen = −8 Mpc to around 0.2 at Δx_cen = −2 Mpc. It indicates that, in this large and collimated filament, the gas accretion, highlighted by the compressive component of the velocity field, contributes to the generation of rotation far away and is then superseded by the contribution of large-scale eddies in the filament as the latter grows in radius. In addition, the slope of the solenoidal component is steeper than that of the compressive component on a 5 Mpc square size map at Δx_cen = −8 Mpc. It is because at this distance there is not much rotation at small scales, only one large eddy, see the vorticity on transverse slices in Fig. 4. Most of the power comes from the laminar flows accreting matter that the 5 Mpc square size map encompasses. Consequently, the power spectrum of the solenoidal component shows power at large scales and then drops rapidly. We observe the same feature in the foreground filament on the 5 Mpc square size map at Δx_cen = 6 Mpc.

Now focusing on the foreground filament (Δx_cen = {2, 4, 6} Mpc, right column of Fig. 9), we observe that, whatever the map size, the solenoidal component dominates the compressive one at Δx_cen = 2 Mpc. Then, the two components have a more or less equivalent amplitude at Δx_cen = 4 Mpc, and the compressive component dominates the solenoidal one at Δx_cen = 6 Mpc. In particular, at this distance, the ratio of compressive to solenoidal is of the order of 5–10 from the 5 Mpc square maps and 1–2 from 1 Mpc square maps. It shows that there is indeed a higher contribution from the compressive component further away from the core of the filament due to gas accretion, as also visible in the right column of Fig. 5. The right column of Fig. 9 shows that regardless of the map size, however, the compressive-to-solenoidal ratio decreases continuously while approaching Virgo. It shows that, in the case of a less collimated filament with multi-stream accretion, there is no obvious eddy around the spine of the filament. The cosmic gas looks to be funnelled towards Virgo, thus generating turbulence in a less structured manner than in the well-collimated background filament, where the gas is accreted from two streams coming from both sides of the matter sheet.

Overall, the conjugate increase in amplitude and steepening of the slope of the velocity power spectrum up to the Virgo outskirts shows that more and more energy is injected into the turbulent cascade because of shocks, which induce more conversion of kinetic energy into heat. In addition, the velocity field transits from a compressive- to solenoidal-dominated regime. All these findings tend to show that turbulence in the filaments develops towards the cluster, finally transiting at the termination shock from a compressible supersonic state, as modelled by Burgers (1948), in filaments to a more developed subsonic and solenoidal-dominated turbulence in Virgo, close to the ideal fully developed, subsonic and incompressible case of Kolmogorov (1941). We compare our findings to the aforementioned models in the following Sect. 6.

6. Discussion

In the following discussion, we compare this work to pioneering ones on gas temperature and velocity in large-scale structures, discuss the impact of computing 2D power spectra instead of 3D and propose interpretations of the results. Finally, we propose future extensions to this work given studies on turbulence in galaxy clusters.

The search for the origin of the acceleration of cosmic rays in the cosmic web and the thermalisation of cosmic baryons led to the study of shocks during large-scale structure formation. Ryu et al. (2003) and Pfrommer et al. (2006) produced cosmological simulations in which they both found signatures of accretion shocks at the boundaries of cosmic sheets, filaments and clusters. These so-called external shocks (Ryu et al. 2003), with a Mach number at shocks estimated from the gas temperature jump following the Rankine-Hugoniot conditions, were up to M = 100 or higher. They also found internal shocks in filaments and haloes with a Mach number around M = 2. They did not extend their studies to the eddies generated by these shocks, however.

Ryu et al. (2008), on the contrary, studied the cascade of vorticity to estimate the generated magnetic field strength. In an MHD simulation reproducing the conditions in the ICM, Porter et al. (2015) tested the generation of turbulence from solenoidal-only and compressive-only forcing; they have shown that both these forcing processes led to a mixture of solenoidal and compressive turbulent motions; in particular, compressive flows generate shocks which induce vorticity. Our results generally agree with theirs as we observe the contribution of both compressive and solenoidal components to the velocity field. Our Virgo replica is a hydrodynamics-only simulation, however, and therefore does not allow us to compare this in detail because the magnetic fields and turbulence in their work are correlated.

In a series of works, turbulence was simulated in the inter-galactic medium (IGM) (Iapichino et al. 2011) and clusters’ outskirts (Iapichino et al. 2017), both using the AMR code ENZO (Bryan et al. 2014) and including a subgrid model for unresolved turbulence, computing the energy content on unresolved length-scales. In Iapichino et al. (2011), they have shown that turbulence is mostly generated by shocks in the WHIM and merger-induced shear flows in the ICM, which agrees with our work. In Iapichino et al. (2017), they particularly studied the shocked ICM in cluster outskirts, focusing on the energy budget, which, in their case, is dominated by the major merger. This both heats the gas and induces turbulence, which we also observed at Δx_cen = { − 2, 2} Mpc, although we are observing the shocked gas accreted from filaments and not the consequences of a merger.

Recently, Lu et al. (2024) studied filaments connecting massive high-redshift galaxies following a similar tomographic approach. In a zoom-in simulation encompassing three ∼1 Mpc filaments in-between ∼10¹² M_⊙ haloes at z ∼ 4, they have shown that their filament could be separated in an outer thermal, intermediate vortex and inner stream zones, the stream zone being the cold gas stream feeding galaxies. Although the methodology is similar to ours, we find very different results; in particular, we observe a large eddy around the spine of the background filament. It comes from the fact that we study much hotter, longer and thicker filaments connected to a ∼4 × 10¹⁴ M_⊙ galaxy cluster at z = 0.

We now compare our results with theoretical predictions. Some studies (e.g. Vazza et al. 2017) focused on the 3D power spectrum, allowing for a direct comparison. We computed 2D power spectra, however, and we therefore need to account for the slope difference between 2D and 3D power spectra. Other works, for instance Zhuravleva et al. (2012) in mock X-ray observations, predicted a proportionality relation, thus a similar slope, between 2D and 3D power spectra. Because we computed the velocity power spectrum on a slice and not in a projection, however, we expect that P_3D(k)k³ ∼ P_2D(k)k², which rewrites as P_3D(k)/P_2D(k)∼k⁻¹. To verify this, we computed the 3D power spectrum in a 1 Mpc side cube in the cluster core region and compared it with its 2D counterpart. Results are presented in Fig. 10. We also compared the spectra in the other slices, which can be found in Appendix F. We did not compute the 3D power spectra on larger cubes to avoid overlap and possible correlations, as we discussed above.

Fig. 10. Comparison between the 2D (solid black line) and 3D velocity power spectra (dashed black line) computed from a 1 Mpc square map and a 1 Mpc side cube, both centred on the core of Virgo. The best-fit value of the slope is displayed as dotted coloured lines.

First, by fitting the slope of the power spectra as P(k) = Ak^α, we have P_3D(k)/P_2D(k)∝k^{α_3D − α_2D}. Computing the difference between the 3D slope and the 2D slope in each case (see Table F.1), we can observe that the difference is always close to −1, which agrees well with the theoretical prediction. Then, we can observe that the 3D power spectrum in the core of Virgo is in excellent agreement with subsonic Kolmogorov (1941) turbulence theory, predicting a −11/3 ≈ −3.67 slope. Moreover, the 3D power spectrum computed around the spine of the filaments is in quite good agreement with the supersonic Burgers’ (Burgers 1948) turbulence theory predicting a −4 slope (see Fig. F.1). It seems to show that, at least in 1 Mpc cube boxes in the core of Virgo and around the spine of the filaments, turbulence has reached the stationary regime. This does not imply that turbulence is in a stationary regime within R_vir ∼ 2 Mpc, however, because of the transition from the supersonic to the subsonic state at this distance. In the filaments, stationary supersonic turbulence might be maintained by the steady accreting flows of matter from the outskirts of the filaments, themselves fed by flows arriving from a cosmic sheet of matter and local voids.

Finally, we could better characterise turbulence in filaments by calculating the Froude number and other stratified turbulence parameters such as the Brunt–Väisälä (BV) frequency (e.g. Mohapatra et al. 2020), conduct a more thorough study of the velocity field by analysing the enstrophy (e.g. Vazza et al. 2017; Vallés-Pérez et al. 2021a) or use continuous wavelets transform (e.g. Wang & He 2024). We plan to conduct such a study in future work where we will also extend our analysis to a larger number of filaments, the Virgo replica and its filaments not being representative of the cluster population as shown in Sorce et al. (2021), and study the time evolution of turbulence in filaments through multiple simulation snapshots (e.g. Wang et al. 2023) when we only studied the Virgo replica at redshift z = 0. This will be done with the upcoming LOCALIZATION⁴ cosmological simulation that reproduces the Local Universe.

7. Conclusion

We studied velocity fields in filaments that are connected to the replicated Virgo cluster by investigating longitudinal and transverse slices, first by visual inspection, and then by decomposing the velocity field and computing the power spectra.

In addition, we compared our work with previous studies on the generation of turbulence in the large-scale structures of the Universe. They agree rather well. Because of the differences in size, resolution, and methods of these works, we were unable to conduct detailed comparisons, however. We also discussed a comparison of the slopes of 2D power spectra to 3D power spectra and theoretical predictions. We showed that the slope differs between our 3D and 2D velocity power spectra, which seems to agree well with theoretical expectations. We found that the slopes of the 3D power spectra agree excellently with Kolmogorov’s subsonic turbulence theory in the core of Virgo and agree well with Burgers’ supersonic turbulence theory around the spine of the filaments. This most likely shows that turbulence reached the stationary regime in these regions. We concluded by proposing improvements and extensions to this work, which we might conduct in the near future.

On the one hand, the study of the longitudinal slice showed how cosmic gas streams in filaments towards Virgo. On the other hand, the investigation of transverse slices revealed the inner structure of the filaments and how gas flows from underdense regions to matter sheets and then filaments. In particular, the background filament is long, well-structured, and connected to a large matter sheet, whereas the foreground filament is much smaller, more diffuse, and funnels the gas towards Virgo. Despite these differences, gas swirls inside both, as shown through the vorticity. Accretion shocks at the filament boundary mostly generate this swirling, which we traced through the divergence of the velocity field. Moreover, the gas accelerates when approaching Virgo and is shocked and heated once again when entering it, roughly at the virial radius. This accelerates the gas swirling.

By computing the 2D total power spectrum, we found an increase in amplitude and steepness towards Virgo in both filaments, with the maxima at Δx_cen = −2 and 2 Mpc, that is, at roughly the virial radius. This means that increasingly more kinetic energy is injected into turbulence by shocks, which is then partly transferred into smaller eddies in the cascade and is partly dissipated into heat. In the Virgo core, the amplitude of the power spectrum is lower than at the virial radius, and the slope is fainter. This tends to show that turbulence reaches a more developed state as it approaches the cluster, and transits from a supersonic and compressible state in the filaments to a more developed subsonic and almost incompressible state in the ICM. This is close to the ideal model of Kolmogorov (1941).

Then, we decomposed the velocity field into compressive and solenoidal components and studied the compressive-to-solenoidal and solenoidal-to-total power spectrum ratios. We found a transition from a compression-dominated to a solenoidal-dominated velocity field when approaching Virgo for both filaments. This reinforced the scenario of transition from Burgers-like to Kolmogorov-like turbulence. Moreover, the power spectra we extracted from the square maps highlighted the different turbulent scales and eddies in these two filaments. In the background filament, the spectra we computed from the 1 Mpc square maps showed that turbulence is dominated by one large eddy in the core of the filament, and the spectra we computed from the 5 Mpc square maps showed that the swirling is generated by gas accretion from the cosmic matter sheet. In the foreground filament, in contrast, the power spectra show the same smooth compressive-to-solenoidal-dominated transition regardless of the map size.

To conclude, the case study of turbulence in the Virgo filaments, the most advanced to date to our knowledge, showed that accretion shocks and instabilities generate gas rotation in cosmic filaments and clusters and that the velocity field develops from a compressive- to a solenoidal-dominated regime. This was possible because of the superior resolution of the simulation. It demonstrates that the velocity field in filaments is a precursor of the velocity field in clusters. It paves the way for future observations of gas dynamics in the core of the clusters and potentially also in outskirts or bridges, in the X-ray with XRISM or Athena, or even in filaments via 21 cm HI with LOFAR or SKA.

Acknowledgments

The authors thank the referee for their helpful comments to improve this article. The authors acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gauss-centre.eu) for providing computing time on the GCS Supercomputers SuperMUC at LRZ Munich. This work has been supported by the grant agreements ANR-21-CE31-0019 / 490702358 from the French Agence Nationale de la Recherche / DFG for the LOCALIZATION project. This work has been supported as part of France 2030 program ANR-11-IDEX-0003. The authors thank Florent Renaud for sharing the rdramses RAMSES data reduction code. The authors thank Eugene Churazov for his precious comments on the velocity power spectrum discussion. SZ thanks Philipp Busch for the early discussions on this topic and acknowledges the hospitality of the Institute d’Astrophysique Spatiale, at the initial phases of this project. SZ also acknowledges support from the Israel Science Foundation grant No. 1388/24.

References

Appendix A: 3D visualisation of the slices

Fig. A.1. 3D visualisation of the temperature in Virgo and its filaments. The green plane represents the longitudinal cut presented in Figs. 1,2 and 3, similarly to the slice shown in the left panel. The red to black to blue planes represent the transverse cuts presented in Figs. 4 and 5. In the right panel, the sheet of matter embedding the filaments is in the bottom-left to top-right diagonal. This clipping and slicing visualisation was made using the PyVista Python library (Sullivan & Kaszynski 2019). A video showing the simulation at different angles is available at this link: https://youtu.be/VZwacpKZNH4.

Appendix B: Electron density in transverse slices

Fig. B.1. Transverse slices of the electron density along the background filament (top) and foreground (bottom) filaments. On the top row, from left to right, we present cuts at Δxcen = { − 8, −6, −4, −2} Mpc. On the bottom row, from left to right, we present cuts at Δxcen = {2, 4, 6} Mpc. Similarly to longitudinal cuts, the maps are 22.122 Mpc wide, contain 157282 pixels, and are centred on the Virgo centre.

Appendix C: Additional quantities in the zoomed regions

Fig. C.1. Zoom into the background (top) and foreground (bottom) filaments in the longitudinal cut. The maps are 11.061 × 12 Mpc large with the same resolution as those presented in Fig. 1. The left panel shows the Mach number, the central panel shows the electron density-weighted divergence, and the right panel shows the z component of the electron density-weighted vorticity. The Virgo R500 and virial radius, Rvir, are shown as black circle arcs on each map.

Appendix D: Power spectra and their ratios computed from 2 Mpc square maps

Fig. D.1. Total power spectra (top) and the ratios of the solenoidal component over the total (middle) and solenoidal over compressive component power spectra (bottom) computed from 2 Mpc square maps, similarly to Figs. 7, 8 and 9.

Appendix E: Range of validity of the power spectra

Fig. E.1. Full scale range of the 3D and 2D power spectra computed in the 5 Mpc cubed size box (squared size map) centred on the core of Virgo, similarly to Fig. 10. The chosen range of validity and the Nyquist frequency are shown by vertical dashed lines.

Appendix F: Comparison of the 3D and 2D power spectra slopes

For each best-fit α parameter, the error, i.e. standard deviation, on the fit is of the order of 1e-3. It is worth noting that the difference between the 3D and 2D slopes is roughly -1, as shown in Tab. F.1 for all the slices-cube pairs.

Fig. F.1. Comparison of the 2D and 3D power spectra slope computed from a 1 Mpc side map or cube with the same centre.

Appendix G: All 2D power spectra

Fig. G.1. 2D power spectrum of the total velocity field (black) and its compressive (blue) and solenoidal (red) components computed from 1 (bottom), 2 (middle) and 5 Mpc (top) square maps extracted from the slices. The squares are displayed in the top row of Fig. 6, and their centre’s position can be found in Tab. 1. From left to right, each column presents the power spectra extracted from a given slice at a distance of between -8 and 6 Mpc from the Virgo centre; they are each 2 Mpc apart.

All Tables

All Figures

Fig. 1. From left to right, we show longitudinal cuts of the temperature, electron density, Mach number, and velocity in the Virgo replica and the two filaments to which it is connected. The maps are about 22 Mpc wide, contain 157282 pixels, and are centred on the Virgo core. In the third panel, the Mach number and colour bands are used to display the wide range of values: The black-to-white colour band is on a linear scale in the [0.1,1] range, and the white-to-blue colour band is on a log scale in the [1100] range. In the right panel, the background map is the velocity along the x-axis of the simulation box, vx, and the foreground arrows represent the norm of the velocity field in the x − y plane, ||v||. The background (foreground) filament is in the left (right) part of each panel. The buoyant bubble is located roughly at a position (x = −5.0; y = 0.5) Mpc. Multiple galaxies are visible in the density map, for example, at positions ( − 2.5; −0.5) and ( − 6.0; 1.0) Mpc.

In the text

Fig. 2. Zoom into the background filament in the longitudinal cut. The maps are 11.061 × 12 Mpc large with the same resolution as those presented in Fig. 1. The left panel shows the velocity field, similar to the right panel of Fig. 1, the central panel shows the divergence, and the right panel shows the z component of the vorticity. The Virgo R500 and virial radius, Rvir, are shown as black circle arcs on each map. The merger-accelerated accretion shock is located roughly in the range x = [ − 5, −2] Mpc and y = [ − 5, −1] Mpc.

In the text

Fig. 3. Zoom into the foreground filament in the longitudinal cut. The maps are 11.061 × 12 Mpc large with the same resolution as those presented in Fig. 1. The left panel shows the velocity field, similar to the right panel of Fig. 1, the central panel shows the divergence, and the right panel shows the z component of the vorticity. The Virgo R500 and virial radius, Rvir, are shown as black circle arcs on each map.

In the text

Fig. 4. Transverse cuts along the background filament. From left to right, we present cuts at distances of 8, 6, 4 and 2 Mpc from the Virgo centre. From top to bottom, we present the temperature, the velocity field with the background map being the velocity along the z simulation box axis, vz, and the foreground arrows representing the norm of the velocity field in the zy plane, ||v||, the divergence of the velocity field and the x component of its vorticity. Similarly to longitudinal cuts, the maps are 22.122 Mpc wide, contain 15 7282 pixels, and are centred on the Virgo centre.

In the text

Fig. 5. Transverse cuts along the foreground filament. From left to right, we present cuts at distances of 2, 4 and 6 Mpc from the Virgo centre. From top to bottom, we present the temperature, the velocity field with the background map being the velocity along the z simulation box axis, vz, and the foreground arrows representing the norm of the velocity field in the zy plane, ||v||, the divergence of the velocity field and the x component of its vorticity. Similarly to longitudinal cuts, the maps are 22.122 Mpc wide, contain 157282 pixels, and are centred on the Virgo centre.

In the text

Fig. 6. Slices of the norm of the 3D velocity field. The top panel presents the total velocity field, and the second (bottom) panel presents the compressive (solenoidal) component obtained via the Hodge-Helmholtz decomposition. The slices are extracted from the 12 × 12 × 20 Mpc box used for the decomposition; the maps are thus 12 × 12 Mpc large, and each pixel is 22.5 kpc large. The slices presented are similar to those in Figs. 4 and 5, that is at Δxcen = { − 8, −6, −4, −2, 0, 2, 4, 6} Mpc.

In the text

Fig. 7. Comparison of the total velocity field power spectra. The top and bottom rows present spectra computed from 5 and 1 Mpc square maps, respectively. The left column presents the total velocity field power spectra in slices in the background filament, that is, from red to black at Δxcen = { − 8, −6, −4, −2} Mpc, and in the Virgo core for comparison. The right column presents the total velocity field power spectra in slices in the foreground filament, that is, from black to blue at Δxcen = {2, 4, 6} Mpc, and also in the Virgo core for comparison. For comparison, we show the slopes predicted by Kolmogorov’s (P(k)∝k−11/3) subsonic turbulence theory and Burgers’ (P(k)∝k−4) supersonic turbulence theory.

In the text

Fig. 8. Ratio of the power spectrum of the solenoidal component over the total velocity field. The top and bottom rows present ratios of spectra computed from, respectively, 5 and 1 Mpc square maps. The left column presents the ratios in slices in the background filament, that is, from red to black at Δxcen = { − 8, −6, −4, −2} Mpc, and in the Virgo core for comparison. The right column presents the ratios in slices in the foreground filament, from black to blue at Δxcen = {2, 4, 6} Mpc, and the Virgo core for comparison.

In the text

Fig. 9. Ratio of the power spectrum of the compressive over the solenoidal component of the velocity field. The top and bottom rows present ratios of spectra computed from, respectively, 5 and 1 Mpc square maps. The left column presents the ratios in slices in the background filament, that is, from red to black at Δxcen = { − 8, −6, −4, −2} Mpc, and in the Virgo core for comparison. The right column presents the ratios in slices in the foreground filament, from black to blue at Δxcen = {2, 4, 6} Mpc, and the Virgo core for comparison.

In the text

	Fig. 10. Comparison between the 2D (solid black line) and 3D velocity power spectra (dashed black line) computed from a 1 Mpc square map and a 1 Mpc side cube, both centred on the core of Virgo. The best-fit value of the slope is displayed as dotted coloured lines.
In the text

Fig. A.1. 3D visualisation of the temperature in Virgo and its filaments. The green plane represents the longitudinal cut presented in Figs. 1,2 and 3, similarly to the slice shown in the left panel. The red to black to blue planes represent the transverse cuts presented in Figs. 4 and 5. In the right panel, the sheet of matter embedding the filaments is in the bottom-left to top-right diagonal. This clipping and slicing visualisation was made using the PyVista Python library (Sullivan & Kaszynski 2019). A video showing the simulation at different angles is available at this link: https://youtu.be/VZwacpKZNH4.

In the text

Fig. B.1. Transverse slices of the electron density along the background filament (top) and foreground (bottom) filaments. On the top row, from left to right, we present cuts at Δxcen = { − 8, −6, −4, −2} Mpc. On the bottom row, from left to right, we present cuts at Δxcen = {2, 4, 6} Mpc. Similarly to longitudinal cuts, the maps are 22.122 Mpc wide, contain 157282 pixels, and are centred on the Virgo centre.

In the text

Fig. C.1. Zoom into the background (top) and foreground (bottom) filaments in the longitudinal cut. The maps are 11.061 × 12 Mpc large with the same resolution as those presented in Fig. 1. The left panel shows the Mach number, the central panel shows the electron density-weighted divergence, and the right panel shows the z component of the electron density-weighted vorticity. The Virgo R500 and virial radius, Rvir, are shown as black circle arcs on each map.

In the text

	Fig. D.1. Total power spectra (top) and the ratios of the solenoidal component over the total (middle) and solenoidal over compressive component power spectra (bottom) computed from 2 Mpc square maps, similarly to Figs. 7, 8 and 9.
In the text

	Fig. E.1. Full scale range of the 3D and 2D power spectra computed in the 5 Mpc cubed size box (squared size map) centred on the core of Virgo, similarly to Fig. 10. The chosen range of validity and the Nyquist frequency are shown by vertical dashed lines.
In the text

	Fig. F.1. Comparison of the 2D and 3D power spectra slope computed from a 1 Mpc side map or cube with the same centre.
In the text

Fig. G.1. 2D power spectrum of the total velocity field (black) and its compressive (blue) and solenoidal (red) components computed from 1 (bottom), 2 (middle) and 5 Mpc (top) square maps extracted from the slices. The squares are displayed in the top row of Fig. 6, and their centre’s position can be found in Tab. 1. From left to right, each column presents the power spectra extracted from a given slice at a distance of between -8 and 6 Mpc from the Virgo centre; they are each 2 Mpc apart.

In the text