Investigating the Intracluster Medium Viscosity Using the Tails of GASP Jellyfish Galaxies

To investigate the properties of the stripped ISM, we select the most extreme ram pressure stripped galaxies from the GASP sample, which is composed of galaxies at z = 0.04–0.07 located in different environments, from the field to clusters. Specifically, we select the 10 galaxies in clusters with the most extended Hα tails in the plane of the sky in the MUSE images. Their properties, and those of their hosting clusters, which are presented in B. M. Poggianti et al. (2017a) and M. Gullieuszik et al. (2020), are summarized in Table 1. The resulting sample spans more than two orders of magnitude in stellar mass, from log (M_star/M_⊙) = 9.05 (JW56) to log (M_star/M_⊙) = 11.47 (JW100) and covers different regions in the phase-space diagram (Figure 1), which indicates that galaxies are in different ICM and RPS conditions.

Zoom In Zoom Out Reset image size

These galaxies, and their stripped tails, have been the subject of previous studies. They all show Hα tails with projected lengths of a few tens of kpc that host extraplanar star formation (B. M. Poggianti et al. 2019b; N. Tomičić et al. 2021). In the MUSE spectra, JW56 and JW100 show extended regions in which the optical line ratios are consistent with LINERS emission (B. M. Poggianti et al. 2019b), indicating that the main ionization mechanism could be the photoionization due to the ICM thermal emission (M. G. Campitiello et al. 2021). For JW100, Chandra observations detected extended X-ray emission in the tail, which correlates spatially with the Hα emission (B. M. Poggianti et al. 2019a) and supports the hypothesis of ongoing ISM–ICM mixing (M. Sun et al. 2021). Finally, in JO206, JW39 and JW100 show nonthermal radio tails (A. Müller et al. 2021; A. Ignesti et al. 2022a), which is evidence of large-scale magnetic fields extending in the tails.

This work is based on the kinematics of the ionized gas derived from the Hα emission. We use the MUSE data from the GASP survey; the observations, data reduction, and analysis are presented in B. M. Poggianti et al. (2017a). The spatial resolution of these seeing-limited observations is ∼1'', which at the redshift of the targets corresponds to ∼1 kpc. Emission line fluxes and gas kinematics were measured using the IDL software KUBEVIZ (M. Fossati et al. 2016) on the data cubes obtained by subtracting the stellar component from the observed data and corrected for dust extinction (for details see B. M. Poggianti et al. 2017a). KUBEVIZ also provides the 1σ error for each fitted parameter, including the line-of-sight velocity. We report the Hα line velocity maps in the galaxy rest-frame by normalizing them for the average velocity measured at the center of the stellar disk.

The stripped tails are complex, filamentary structures in which the properties of the plasma are affected by the presence of star-forming complexes (B. M. Poggianti et al. 2019b; E. Giunchi et al. 2023), and by other large-scale dynamical processes which imprint coherent velocity differences in the VSF. Thus, it is necessary to further process the data to recover the turbulence signature in the Hα VSF.

2.3.1. Step I: Locating the Diffuse ISM

2.3.2. Step II: Filtering Out Rotation and Advection

The stripped ISM dynamic outside the stellar disk is primarily set by the combination of (1) advection induced by the ram pressure, (2) rotation inherited from the disk, and (3) ISM small-scale motions (e.g., M. Gullieuszik et al. 2017; S. Tonnesen & G. L. Bryan 2021; N. Akerman et al. 2023; R. Luo et al. 2023; M. Sparre et al. 2024). To study the ISM small-scale motions, and hence the ICM turbulence signature in the VSF, it is crucial to remove the other two factors. However, their removal is complicated by the intrinsic filamentary structure of the stripped ISM, which makes it difficult to define an unambiguous rotation axis to conduct an analysis of the 3D rotation (e.g., C. Bacchini et al. 2023, for the analysis of the rotating molecular gas in the disks of the same GASP galaxies), and by projection effects, which can lead gas originating from different sides of the disk to overlap along the line of sight. Therefore, we propose a novel technique to remove these components from the velocity field, which we show here for JO206, and for the other galaxies in Appendices A and B.

To visualize the structure of the velocity maps, we compose the phase-space diagram for each galaxy, where we compare the velocity measured in each pixel with its minimum distance from the stellar disk edge, as defined in M. Gullieuszik et al. (2020). The diagrams are then converted into density plots to highlight the substructures in the velocity maps (Figure 2, top panel). First, we correct the velocity field for the advection component along the line of sight. The ram pressure can induce a velocity increase with the distance from the stellar disk as the gas clouds gradually approach the ICM wind speed (R. Luo et al. 2023; P. Serra et al. 2023). So, we compute the median velocities at increasing distances from the disk and use them to fit a linear model of the velocity gradient along the tail (Figure 2, top panel, orange points). Then, we subtract this systematic component from each pixel depending on their distance from the disk. A similar procedure is presented also by R. Luo et al. (2023) to correct the systematic velocity in ESO137-001. This correction can also remove the systematic velocity shift derived from the galaxy velocity in the cluster. The result is a velocity distribution that is typically centered around zero without clear gradients with the distance (Figure 2, central panel). However, we recognize the presence of clear filamentary patterns in the residual phase-space intensity distribution. We assume that these structures are produced by the rotation of the stripped gas inherited from the disk kinematics.

Zoom In Zoom Out Reset image size

To remove the rotation, we first identify the different filaments in the phase-space density plots as regions where the points are clustered around certain velocity values. To do so, we bin the phase-space plots in several velocity bins and then for each bin we compute the median velocity values at increasing distances from the disk. The numbers of velocity and distance bins are different for each galaxy and they are set by visual inspection. For each filament (colored points in Figure 2, central panel), we compute the average velocity, weighted by the number of points in each bin (colored lines in Figure 2, central panel). Then, we correct each filament by subtracting the corresponding average velocity from each pixel in their bins. At this step, we mask out those pixels located outside the filaments. This operation independently corrects each filament for its velocity, resulting in a velocity distribution centered around zero and without any clear symmetric structure (Figure 2, bottom panel). Finally, we further mask out the strong outliers, which typically are residuals left by rotation correction or contaminations from the [N ii] line. In Figure 3, we show the original velocity map of the diffuse, stripped ISM (top panel) and the model obtained by combining the advection and rotation components (bottom panel).

Zoom In Zoom Out Reset image size

2.3.3. Step III: Computing the VSF

The VSFs are computed from the pixels of the residual velocity map as follows: for each pair of pixels, we measure their projected physical distance l, which we derive by multiplying their angular separation by the kpc-to-arcsec scale at the hosting cluster redshift, and the velocity difference δ V. The VSF is derived by computing the average absolute value of the velocity differences 〈∣δ V∣〉 within bins of l from 0.2 to 100 kpc with a step of 0.2 kpc, which is approximately the physical scale corresponding to the pixel size of the MUSE images. Correspondingly, by propagating the uncertainties on the velocity measurement of each spaxel, the uncertainties on the fit of the systematic velocity, and the rotating filaments, we compute the error on each bin of the VSF that results in the order of a few km s⁻¹. We show in Figure 4 the residual velocity map (left-hand panel) and the corresponding VSF (blue line). For comparison, we use the original velocity map not corrected for the rotation and advection to compute the VSF in the tail (VSF_NC, orange-dashed line), and inside the stellar disk (VSF_Disk, green-dashed–dotted line). The VSFs show different slopes, with the first one resembling the slope l^1/3 (gray-dashed line), and the other ones being generally steeper and more similar to the linear scaling (black line). We note that as a consequence of atmospheric seeing the effective angular resolution is limited to ∼1'' (≃1 kpc marked in gray in Figure 4), and thus below that scale the VSF signal is less reliable because the pixels are correlated within the resolution elements. Therefore, in the following analysis we conservatively consider l = 1 kpc as the lower limits of the VSFs. As a caveat, it is important to remember that this analysis is subject to projection effects because filaments that are physically distant can appear adjacent along the line of sight. We refer to Section 4 for a detailed discussion of the role of projection effects.

Zoom In Zoom Out Reset image size

To further corroborate our analysis, we explore the case of the VSF in the absence of ICM turbulence. To do so, we make use of the simulation setup presented in N. Akerman et al. (2023) (see Appendix D). The simulation models stripping of a massive Milky Way-like galaxy with a radius of ∼30 kpc and which rotates at ∼200 km s⁻¹. In this set-up, the gas dynamics are set solely by rotation and ram pressure, and the simulated ICM wind is initially laminar at the galaxy scale. Furthermore, there are no magnetic fields. We produced maps of the gas velocity along the line of sight for the gas phase in the temperature interval 3.5 < logT <6.5, each one at four different snapshots in time (100, 200, 300, and 400 Myr after the beginning of the stripping). The maps capture the gas in the stripped tails, defined as 2 kpc above the galaxy plane in the direction of stripping. Specifically, we analyze the case of face-on stripping (W0) and inclined stripping (W45, wind hits the galaxy at 45°). In W0, the wind blows along the z-axis, while in W45 the wind has components along both the y- and the z-axes. The latter ensures that the line of sight is aligned with the stripping direction and allows us to further investigate the role of projection effects. For the simulated VSFs, we take the gas projections along the x-axis for both W0 and W45 galaxies, and an additional projection along the y-axis for W45 ("W45-y" in Figure 5). We set the pixel size of the velocity maps to 0.2 kpc to match the physical scale of the MUSE images. The VSFs are computed with the same method adopted for the real MUSE observations, i.e., by considering only the gas outside of the stellar disk. We note that at this stage we do not attempt to remove the rotation and advection. The resulting VSFs are shown in Figure 5, divided by simulation run (rows). The velocity maps themselves can be found in Appendix D.

Zoom In Zoom Out Reset image size

In the early stages (100–200 Myr), when the tail is still rich in rotating gas (see Figure 14), the VSFs show a recurrent shape characterized by an almost constant slope from 0.2 to 2–3 kpc, a strong inflection between 3 and 7 kpc, and linear slope expected by the rotation up to the disk radius scale. As a cautionary tale we note that for the late-stage stripping, when the tail has already lost most of the gas and what is left does not keep track of the disk rotation anymore, the VSFs are flatter than the initial-stage stripping case. This result suggests that an advection-dominated tail, in which the gas particles move with a uniform velocity set by the ram pressure wind, results in a sublinear VSF, which is similar to a turbulent-dominated one.

Finally, we use the simulated tails to demonstrate the effectiveness of the correction procedure described in Section 2.3.2. Specifically, we perform the procedure on the W45-y run, at 300 Myr after the beginning of the stripping because of the complex morphology composed of filamentary structures visually resembling the observed tails. We show each step and the final results in Figure 6. After the removal of rotation and advection, the residual VSF (blue line) results are significantly flattened, which denotes the fact that the residual velocity map is devoid of large-scale coherent velocity differences, thus confirming the efficacy of the procedure. We note that the residual VSF is also flatter than the turbulent slope, which indicates that the simulated tails lack small-scale turbulent motions. We conclude that this is likely due to the limits of the simulation setup (see Appendix D), which cannot fully resolve the turbulence development. Although done intentionally because such a configuration is designed to save on computational time and resources, and furthermore the simulation was designed for a different scientific analysis (see N. Akerman et al. 2023), this means that future studies aiming to focus on the development of stripped tails will require tailored numerical simulations.

Zoom In Zoom Out Reset image size

In Figure 7, we show the resulting, corrected VSF for each galaxy compared with the Kolmogorov spectrum 〈∣δ V∣〉 ∝ l^1/3 (A. Kolmogorov 1941). The individual VSFs are presented in Appendix C. Generally, the galaxies approximately follow the ∝l^1/3 scaling below 5 kpc, whereas they are flat on larger scales due to the removal of rotation and advection, which dominate the large-scale dynamic. This pattern is consistent with the one shown by Y. Li et al. (2023).

Zoom In Zoom Out Reset image size

The residual VFSs are different from the simulated ones (Figure 5), lacking both the steep, high-velocity part at large scales, induced by the combination of rotation and advection, and the flattening at small scales. On the contrary, these features are present in the not-corrected VSFs (Appendix C, orange lines). Therefore, we argue that we have removed most of the rotation and advection components, and thus the residual VSFs trace the small-scale, turbulent motions of the stripped ISM.

The turbulent cascade we observe could either result from the "frozen-in" ISM turbulence inherited from the disk or the ICM turbulence. To test the first scenario, in Figures 4 and 13 we compare the VSF computed inside (green lines) and outside (orange lines) of the stellar disk. We observe two cases, galaxies in which the disk VSF is lower than the tail one (JO113, JO162, JO175, JO206, and JW39), and galaxies in which we observe the opposite (JO147, JO171, JO204, JW56, and JW100). The first case naturally challenges the frozen-in ISM turbulence scenario because it indicates that the ISM turbulence can be weaker and the tail turbulence is likely dominated by the growth of instabilities induced by the interaction with the ICM. Regarding the second case, we observe that the ISM turbulence can be significantly strong (e.g., the case of JW56 and JW100), but once the gas is stripped it decays quickly. For JW56 and JW100, the flattening of the VSF takes place at the 2–3 kpc scale with δ v ∼ 70 km s⁻¹. The driving scale is typically a factor 2–3 larger than the flattening scale (C. Zhang et al. 2024), thus we estimate an eddy turnover time (t ≃ l/v) of ∼50–100 Myr. Assuming an average stripped plasma velocity of 300 km s⁻¹, which is in line with the constraints provided by the decay of the radio continuum emission in ram pressure stripped tails (A. Ignesti et al. 2023; I. D. Roberts et al. 2024), then within one eddy turnover time the stripped ISM has moved about ∼15–30 kpc. This length scale is consistent with the observed tail length. So, we argue that what we observe in the tail is never pure, frozen-in ISM turbulence but we are tracing fully mixed ISM or, at least, stripped ISM in a transition stage.

To further understand how turbulence develops in the tails of GASP galaxies, we study the variation of the residual VSF at increasing distances from the stellar disk. We perform this analysis on JO147, JO204, JO206, and JW100, the galaxies with the longest tails in the sample. Similar to the approach presented by Y. Li et al. (2023), we define an adaptive moving frame, with a width corresponding to one-third of the tail length, in which we compute the VSF in eight bins at increasing distance from the stellar disk. We also compute the average value of each VSF between 2 and 3 kpc to quantify the trend with the distance. The results are shown in Figure 8.

Zoom In Zoom Out Reset image size

The resulting VSFs broadly follow the l^1/3 slope below 5 kpc, as Figure 7 shows. We observe a gradient with the distance, with the farthest bins showing a higher level of turbulence with respect to those close to the disk. This result is in line with a scenario where the mixing with the ICM is increasing the stripped ISM turbulence velocity with respect to the initial stellar disk conditions (F. Fraternali et al. 2002; G. Iorio et al. 2017; C. Bacchini et al. 2020). Therefore, the observed VSFs support our hypothesis that environmental turbulence ultimately drives the stripped gas small-scale motions.

Our working hypothesis, based on the observational evidence of ISM–ICM mixing (see Section 1), is that the stripped ISM clouds can trace the local ICM turbulence. The increase in turbulence power with the distance from the disk and the fast dissipation of the small-scale turbulence with respect to the dynamical time, both discussed in Section 3.2, indicate that the turbulence observed the residual VSFs is generated outside of the disk, which is consistent with our assumption. Therefore, similarly to the case discussed by Y. Li et al. (2023), the observed, residual VSFs can provide us with new insights into the ICM properties. The ICM turbulent cascade is supposed to extend down to a dissipation scale beyond which the viscous forces dominate over the inertial ones, and thus the kinetic energy gets dissipated into heating. Therefore, in the classical treatment of the ICM as a fluid driven by Coulomb collision, the smallest VSF scale would be the so-called Kolmogorov microscale η:

$$\begin{eqnarray}&&\eta ={\left(\displaystyle \frac{{\nu }^{3}}{\epsilon }\right)}^{\tfrac{1}{4}}\ \mathrm{cm},\end{eqnarray} \tag{ 1 }$$

where is the specific energy flux and ν is the kinematic viscosity. The latter is related to the dynamic viscosity, μ, which can be computed as:

$$\begin{eqnarray}&&\mu =\rho \nu =5500{\left(\displaystyle \frac{{T}_{e}}{{10}^{8}\ {\rm{K}}}\right)}^{\tfrac{5}{2}}{\left(\displaystyle \frac{\mathrm{log}{\rm{\Lambda }}}{40}\right)}^{-1}\ {\rm{g}}\ {\mathrm{cm}}^{-1}{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 2 }$$

where logΛ is the Coulomb logarithm (C. L. Sarazin 1988, Equation (5.33)), and ρ and T_e are, respectively, the ICM mass density and electron temperature (e.g., L. Spitzer 1962). The specific energy flux can be computed as:

$$\begin{eqnarray}&&\epsilon \simeq \displaystyle \frac{{v}^{3}}{l}\ {\mathrm{cm}}^{2}{{\rm{s}}}^{-3}.\end{eqnarray} \tag{ 3 }$$

To compute η it is thus necessary to know the local ICM electron density and temperature near the galaxies, which can be provided by the X-ray spectral analysis of the ICM thermal emission. We collect these values for eight out of ten galaxies, for which an X-ray spectral analysis is reported in the literature (see Table 2). For three of them, namely JO113, JO147, and JO162, only the ICM temperature is available, thus we compute a range for ν for the typical range of electron density n_e = (0.1 – 1) × 10⁻³ cm⁻³. Finally, for each galaxy, we estimate by using the VSF value at the reference scale of 3 kpc, i.e., the scale at which every galaxy shows a Kolmogorov-like slope. We note that for v ∝ l^1/3, is scale-invariant. We estimate η to be <100 kpc, which is in line with the previous estimates (I. Zhuravleva et al. 2019; Y. Li et al. 2023). For comparison, we also compute the electron mean free path λ_e following L. J. Spitzer (1962):

$$\begin{eqnarray}&&{\lambda }_{e}=\displaystyle \frac{{3}^{3/2}{\left({{kT}}_{e}\right)}^{2}}{4{\pi }^{1/2}{n}_{e}{e}^{4}\mathrm{log}\,{\rm{\Lambda }}}\ \mathrm{cm},\end{eqnarray} \tag{ 4 }$$

where k is the Boltzmann constant, T_e and n_e are the ICM electron temperature and particle density, and e is the electron charge. For JO113, JO147, and JO162, we compute the λ_e range corresponding to n_e = (0.1– 1) × 10⁻³ cm⁻³. We show the values of temperature, kT_e , and electron density, n_e , at the cluster-centric distance of each galaxy, VSF at 3 kpc, and the resulting η and λ_e in Table 2.

To test if the VSFs stop at the expected dissipation scale, we rescale each VSF for the corresponding η, and, for JO113, JO147, and JO162, we use the lower limit of the interval of η (Table 2). The resulting VSFs are presented in Figure 9, in which we also show their components below l = 1 kpc (dashed lines). It emerges that the VSFs easily extend below the expected dissipation scale, which corresponds to l/η = 1 (red-vertical line), with JO113, JO171, and JO175 completely developing in the regime l/η < 1. This result is consistent with the previous study by Y. Li et al. (2023), and it confirms that the ICM dissipation scale is smaller than predicted by Equation (1).

Zoom In Zoom Out Reset image size

The observed and predicted ICM viscosity discrepancy, ν/ν_S , where ν_S is the Spitzer viscosity (Equation (2)), can be inferred by comparing the observed dissipation scale, l_d , with the predicted one, η. Specifically, from Equation (1) it follows that:

$$\begin{eqnarray}&&\displaystyle \frac{\nu }{{\nu }_{S}}={\left(\displaystyle \frac{{l}_{d}}{\eta }\right)}^{\tfrac{4}{3}}\ .\end{eqnarray} \tag{ 5 }$$

We note that in our work we cannot reliably sample the VSF below ∼1 kpc due to the atmospheric seeing (see Section 2), thus we conservatively estimate the upper limit for ν/ν_S by setting l_d = 1 kpc. We report the resulting log ₁₀(ν/ν_S ) in Table 2. Our results indicate that the real ICM viscosity is smaller than the Spitzer estimate, with a difference that spans, at least, two orders of magnitude from log ₁₀(ν/ν_S ) = −0.6 (JW100) to log ₁₀(ν/ν_S ) = −2.7 (JO171), in line with the previous results (I. Zhuravleva et al. 2019; Y. Li et al. 2023).

This discrepancy could be evidence that the ICM particles' mean free path is smaller than the predicted λ_e (Table 2) due to plasma effects. Since we are measuring the VSFs at scales that are similar to, or smaller than, the Coulomb scale λ_e , for completeness we note that it is also possible that self-organization effects in collisionless turbulence (e.g., magneto-immutability, J. Squire et al. 2019) may allow a magnetohydrodynamic-like turbulent inertial range to also form below the putative viscous scale (S. Majeski & M. W. Kunz 2024). In this case, rather than a reduced effective λ_e , we are probing effects due to collisionless physics in the ICM turbulence. We also note that the viscosity can be further suppressed of a factor ∼3 by magnetic fields (e.g., F. Mogavero & A. A. Schekochihin 2014; S. Komarov et al. 2018). The presence of magnetic fields in the jellyfish galaxies' tails is demonstrated by the existence of nonthermal radio continuum emission, which constrains the typical intensity to be between 4 and 7 μG (e.g., H. Chen et al. 2020; I. D. Roberts et al. 2021a, 2021b; A. Ignesti et al. 2022a, 2022b, 2023).

The presence of turbulence in the stripped ISM and the detection of a minimum scale can have implications for local star formation (e.g., M.-M. Mac Low & R. S. Klessen 2004; D. M. Elmegreen et al. 2014; C. Federrath et al. 2021) because below the dissipation scale turbulence does not significantly contribute to the gas internal pressure that contrasts the gravitational collapse. For JO175, JO204, JO206, JW39, and JW100, E. Giunchi et al. (2023) report the size distribution of star-forming clumps, finding that those in the tails of these galaxies are indeed limited to sizes of 0.3 kpc, which would be in agreement with the dissipation scales derived in our work. Remarkably, the complexes in the tails are smaller than those in the stellar disk, hinting at different underlying processes that shape the origin of star-forming complexes in these different environments.

Additionally, in Figure 7 we observe that the VSF, measured at the reference value of 3 kpc, ranges from ∼30 to 90 km s⁻¹. This gradient could be due to the work impressed by the ram pressure during the galaxies' orbits, which can be tentatively outlined by the position of these galaxies in the phase-space diagram. In this diagram, which we show in Figure 1, galaxies in the top left-hand corner show the highest velocity and the smallest cluster-centric distance, thus they should be subject to the strongest ram pressure. We note, however, that the phase-space coordinates are only lower limits of the corresponding quantities. JW100 and JW56, which show the highest line-of-sight velocities and the smallest projected distances, also show the highest 〈∣δ V∣〉_{3 kpc} values, namely 79 and 63 km s⁻¹ (Table 2). Therefore, the data tentatively support the hypothesis of a connection between ram pressure and VSF amplitude, but a larger sample is required to confirm it.