VERTICO V: The environmentally driven evolution of the inner cold gas discs of Virgo cluster galaxies

Abstract

The quenching of cluster satellite galaxies is inextricably linked to the suppression of their cold interstellar medium (ISM) by environmental mechanisms. While the removal of neutral atomic hydrogen (H i) at large radii is well studied, how the environment impacts the remaining gas in the centres of galaxies, which are dominated by molecular gas, is less clear. Using new observations from the Virgo Environment traced in CO survey (VERTICO) and archival H i data, we study the H i and molecular gas within the optical discs of Virgo cluster galaxies on 1.2-kpc scales with spatially resolved scaling relations between stellar ($\Sigma_{\star}$), H i ($\Sigma_{\text{H}\,{\small\text{I}}}$), and molecular gas ($\Sigma_{\text{mol}}$) surface densities. Adopting H i deficiency as a measure of environmental impact, we find evidence that, in addition to removing the H i at large radii, the cluster processes also lower the average $\Sigma_{\text{H}\,{\small\text{I}}}$ of the remaining gas even in the central $1.2\,$kpc. The impact on molecular gas is comparatively weaker than on the H i, and we show that the lower $\Sigma_{\text{mol}}$ gas is removed first. In the most H i-deficient galaxies, however, we find evidence that environmental processes reduce the typical $\Sigma_{\text{mol}}$ of the remaining gas by nearly a factor of 3. We find no evidence for environment-driven elevation of $\Sigma_{\text{H}\,{\small\text{I}}}$ or $\Sigma_{\text{mol}}$ in H i-deficient galaxies. Using the ratio of $\Sigma_{\text{mol}}$-to-$\Sigma_{\text{H}\,{\small\text{I}}}$ in individual regions, we show that changes in the ISM physical conditions, estimated using the total gas surface density and midplane hydrostatic pressure, cannot explain the observed reduction in molecular gas content. Instead, we suggest that direct stripping of the molecular gas is required to explain our results.

1. Introduction

Galaxies experience different evolutionary processes depending on the environment in which they live, leading to variation in their observed properties as a function of environment. In particular, galaxies in denser environments contain a systematically lower mass of cold, gaseous interstellar medium (ISM) than isolated systems (e.g., Haynes & Giovanelli 1984; Giovanelli & Haynes 1985; Solanes et al. 2001; Fumagalli et al. 2009; Cortese et al. 2011; Catinella et al. 2013; Brown et al. 2017; Zabel et al. 2019; Stevens et al. 2019b, 2021). As this cold gas is the fuel for star formation (e.g., Schmidt 1959; Kennicutt 1989; Leroy et al. 2008; Bacchini et al. 2019; Saintonge & Catinella 2022), completing our understanding of the physics of galaxy evolution requires understanding how environmental mechanisms impact their cold gas content.

Galaxy clusters are one of the densest environments in the Universe, where galaxies are subject to numerous environmental mechanisms (e.g. Boselli & Gavazzi 2006). Their deep gravitational potentials cause tidal stripping of satellite galaxies near the cluster centre (Byrd & Valtonen 1990; Fujita 1998), and give galaxies large relative velocities ( $>$ $500\,\text{km} \, \text{s}^{-1}$ , Girardi et al. 1993; Struble & Rood 1999), causing them to undergo fast gravitational interactions such as galaxy harassment (Moore et al. 1996; Mihos 2004). Additionally, clusters are permeated by a hot ( $>$ $10^7\,$ K) intra-cluster medium (ICM, e.g., Sarazin 1986; Mohr, Mathiesen, & Evrard 1999; Boselli & Gavazzi 2006; Frank et al. 2013). Interaction with this ICM results in galaxies undergoing hydrodynamical processes such as thermal evaporation of the ISM (Cowie & Songaila 1977), turbulent stripping by viscous momentum transfer between the ICM and ISM (Nulsen 1982), and/or ram-pressure stripping (RPS, Gunn & Gott 1972). As the removal of gas results in the future cessation of star formation, the action of these numerous environmental processes is often invoked to explain the large fractions of quenched galaxies hosted by clusters (Dressler 1980; Kauffmann et al. 2004; Peng et al. 2010; Woo et al. 2013). While the timescale for quenching depends on the gas removal mechanism and the circularity of the in-fall orbit (Font et al. 2008; Jung et al. 2018), cluster satellites are consistent with their star formation being significantly reduced before or at first pericentre passage (e.g., Crowl & Kenney 2008; Ciesla et al. 2016; Boselli et al. 2016a; Vulcani 2020; Werle et al. 2022), suggesting a fast gas removal process (e.g., Boselli & Gavazzi 2006; Fossati et al. 2018; Boselli et al. 2020; Upadhyay, Oman, & Trager 2021; Wright et al. 2019, 2022).

Indeed, there is significant evidence that hydrodynamical stripping is a dominant mechanism for removing gas from cluster galaxies (see the reviews by Cortese, Catinella, & Smith 2021; Boselli, Fossati, & Sun 2022). In the canonical outside-in picture of stripping, the circumgalactic medium is first removed and cuts galaxies off from further accretion (Larson, Tinsley, & Caldwell 1980), before progressing to strip the ISM directly (e.g., Kenney, van Gorkom, & Vollmer 2004; Chung et al. 2007; Cortese et al. 2010). Further, cometary tails of stripped gas from the ISM have been observed with multi-wavelength tracers from the X-rays to the radio continuum (e.g. Sivanandam et al. 2014; Jáchym et al. 2014, 2017; Moretti et al. 2020b; Longobardi et al. 2020; Roberts et al. 2021; Ignesti et al. 2022), and reproduced in controlled and cosmological simulations of galaxy evolution (e.g. Tonnesen & Bryan 2010; Marasco et al. 2016; Yun et al. 2019).

Most of our understanding of how environmental processes impact the cold gas in galaxies comes from observations of the 21-cm emission line from neutral atomic hydrogen (H i), as its large radial extent makes it a sensitive tracer of external environmental mechanisms. In the case of hydrodynamical stripping, galaxies exhibit clear truncation of their H i distribution at a radius that depends on the strength of the ram pressure (Gunn & Gott 1972; Cayatte et al. 1990; Roediger & Hensler 2005; Tonnesen & Bryan 2009; Abramson et al. 2011; Yoon et al. 2017; Reynolds et al. 2021, 2022). The impact of these environmental processes on the molecular hydrogen (H $_{2}$ ) is less clear; however, truncation of the star formation in stripped galaxies (e.g., Koopmann & Kenney 2004; Boselli et al. 2006; Cortese et al. 2012; Fossati et al. 2013, 2018) indicates that the H $_{2}$ is also affected once stripping proceeds within the stellar disc. Further, the presence of star formation inside this truncation radius means that some gas can survive the stripping process. Whether the environment has also influenced this surviving gas, which would also affect the timescale for star formation to halt e.g. (e.g. Boselli et al. 2014; Cortese et al. 2021), remains unknown. As this surviving gas is often dominated by H $_{2}$ (Boselli et al. 2014; Cortese et al. 2016; Loni et al. 2021), understanding the impact of environmental processes on the molecular ISM in galaxies is essential to completing our picture of galaxy evolution in dense environments.

Under most ISM conditions, however, H $_{2}$ cannot be directly observed as it is a symmetric molecule without a permanent dipole. Instead, most of our understanding of the molecular ISM (from here, ‘molecular gas’ or subscript ‘mol’) comes from observing the rotational transitions of CO and its isotopologues (Bolatto et al. 2013; Saintonge et al. 2017). These observations have revealed that molecular gas is not affected by environmental processes as strongly as H i. This resistance is primarily because a large fraction of molecular gas exists in giant molecular clouds (GMCs), and it is typically more centrally concentrated where the gravitational restoring force is stronger, making it more resistant to direct stripping than the H i (e.g., Engargiola et al. 2003; Bolatto et al. 2008; Yamagami & Fujita 2011; Stevens et al. 2021). Instead, environmental processes typically lower the molecular gas content in the outskirts of galaxies first (Lee et al. 2020; Brown et al. 2021; Zabel et al. 2022), while significant reduction of the molecular gas mass of galaxies is only observed once galaxies have lost a large fraction of their total H i mass (Fumagalli et al. 2009; Boselli et al. 2014).

Although this picture is consistent with the outside-in stripping scenario, there is growing evidence that it is more complicated than just radial truncation. Simulations have shown that the multi-phase nature of the ISM allows the ICM to stream through low-density regions of the surviving gas inside the H i truncation radius, allowing environmental processes to affect the cold gas across the disc (Tonnesen & Bryan 2009; Choi, Kim, & Chung 2022). Disturbed molecular gas morphologies in cluster galaxies can resemble the H i (e.g., Lee et al. 2017), and high-resolution studies of individual objects have shown GMCs decoupled from large plumes of trailing gas (e.g., Cramer et al. 2020). Further, the mass of molecular gas detected in the tails of stripped gas trailing behind cluster galaxies requires some level of mass-loading from the ISM (e.g., Sivanandam, Rieke, & Rieke 2010; Moretti et al. 2018; Lee & Chung 2018; Jáchym et al. 2019), and there is evidence for the direct stripping of cold ISM from dust observations (e.g., Crowl et al. 2005; Longobardi et al. 2020). Alternately, stripping of the molecular gas may not be a dominant process, and the reduction of molecular gas mass could be due to the preferential removal of H i across the disc, leaving insufficient fuel to maintain efficient molecule formation (Fumagalli et al. 2009; Tonnesen & Bryan 2009; Zabel et al. 2022; Villanueva 2022). To further complicate the picture, there is evidence for increased molecular gas formation efficiency due to gas compression during the ICM–ISM interaction (e.g., Tonnesen & Bryan 2009; Moretti et al. 2020a,b). Clearly, our picture of how the environment impacts the molecular gas in galaxies is incomplete.

In this work, we present a new analysis of the cold gas properties within the stellar components of Virgo cluster galaxies using observations from the Virgo Environment Traced in CO(2-1) (VERTICO, Brown et al. 2021) ALMA ACA^a large program. The Virgo cluster is an excellent target for furthering our understanding of the environment’s impact on galaxies’ cold gas content. At a distance of 16.5 Mpc (Mei et al. 2007) it is close, allowing for detailed studies of individual galaxies, and it has been observed comprehensively across nearly the entire electromagnetic spectrum (York et al. 2000; Chung et al. 2009; Davies et al. 2010; Urban et al. 2011; Boselli et al. 2011, 2018a). The presence of infalling sub-structures also means that there are populations of galaxies at various stages of infall, and thus environmental impact (Gavazzi et al. 1999; Mei et al. 2007; Kim et al. 2016; Yoon et al. 2017). Additionally, the VERTICO sample was selected from the VLA Imaging of Virgo galaxies in Atomic gas survey (VIVA, Chung et al. 2009) of H i emission. Together, these data enable us to study how the environmental processes influence the cold gas properties of galaxies in the Virgo cluster on $\sim$ 1.2-kpc scales.

We present our study using the framework of spatially resolved scaling relations (e.g., Bigiel et al. 2008; Morselli et al. 2020; Ellison et al. 2021a,b; Pessa et al. 2022; Abdurro’uf et al. 2022a,b; Baker et al. 2022), in particular those between stellar, H i, and molecular gas mass surface densities. While azimuthally averaged radial profiles are useful for increasing sensitivity to faint emission, in this work, we focus on the bright emission within the optical components of galaxies. Thus, comparing the physical quantities of kpc-scale regions within galaxies increases our statistics, and reduces the effects of averaging over larger areas that include many different physical regions. For analyses of how the environment impacts the radial profiles of cold gas content in VERTICO galaxies, we refer the reader to the works of Brown et al. (2021), Zabel et al. (2022), and Villanueva 2022).

This paper is structured as follows. In Section 2 we describe our VERTICO sample, control sample, the homogenisation of the data, and ancillary data products. In Section 3 we investigate the resolved scaling relations with all galaxies combined and discuss some of the caveats of this approach. We present our main results in Section 4, where we study how the environment impacts the spatially resolved H i and molecular gas content of galaxies across their stellar discs and within the truncation of their gas discs. In Section 5 we investigate the physical mechanisms responsible for our observations, and in Section 6 we summarise our results.

As described in Brown et al. (2021), all physical quantities used in this work were computed using a Kroupa (2001) stellar initial mass function, a $\Lambda$ CDM cosmology with $h=0.7$ , $\Omega_\mathrm{M}=0.3$ and $\Omega_\Lambda=0.7$ , and we adopt a constant distance for Virgo galaxies of 16.5 Mpc (Mei et al. 2007).

2. Datasets

2.1 VERTICO sample

We selected our main galaxy sample from the VERTICO survey (Brown et al. 2021), which obtained CO(2-1) observations ^b for 51 Virgo cluster galaxies extracted from the VIVA survey. We refer the reader to Brown et al. (2021) for a detailed description of the parent sample selection and data reduction. In this work, we focus on massive ( $\log(M_{\star}/\mathrm{M}_{\odot})\geq10$ ) galaxies to minimise metallicity-dependent uncertainties in tracing the molecular gas content of galaxies with CO(2-1) emission, and we enforce an optical inclination cut of $i\leq70^\circ$ to avoid large de-projection uncertainties. We also removed two galaxies that met these criteria, NGC 4772 and NGC 4698, as they are passive, early-type discs that have recently accreted their H i gas (Haynes et al. 2000; Cortese & Hughes 2009; Yoon et al. 2017), and thus their H i and molecular gas properties do not reflect the action of the Virgo environment on their cold gas content. Our final VERTICO sample consists of 18 galaxies, and the properties of the galaxies used in this work are summarised in Table 1, along with their sources and derivations that are also described below.

Table 1. Properties of the galaxies used in this work. VERTICO galaxies are contained in the top section, while the HERACLES control sample is in the bottom segment.

$^{\mathrm{a}}$ VERTICO data are 7-m array only, the rest have 7-m and Total Power array data.

Columns: (1) Galaxy name; (2) Luminosity distance (Mei et al. 2007; Walter et al. 2008); (3) Position angle (Brown et al. 2021); (4) inclination (Brown et al. 2021); (5) B-band optical radius (Chung et al. 2009; Walter et al. 2008); (6) $\log_{10}$ integrated stellar mass (Leroy et al. 2019); (7) $\log_{10}$ integrated star-formation rate Leroy et al. 2019; (8) H i deficiency parameter (Chung et al. 2009), this work; (9) average $\log\Sigma_{\star}$ at the H i-radius $\Sigma_{\text{H}\,{\small\text{I}}}=1\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ (this work).

The native spatial resolution of VERTICO data is $\sim$ 8-arcsec, and data products with 9-arcsec and 15-arcsec resolution were also created using the imsmooth routine of the Common Astronomy Software Applications package (CASA, McMullin et al. 2007). We used the coarser 15-arcsec data products sampled with 8-arcsec spaxels, which are closer to the spatial resolution of the H i data. The VIVA H i data used in this work were re-imaged by author A. Chung to circularise the synthesised beam and match the VERTICO data products as closely as possible, resulting in datacubes with 8 arcsec spaxels and synthesised beams varying from 15.2 to 20 arcsec, with an average of 17 arcsec. At the assumed 16.5 Mpc distance to Virgo, the 15/8-arcsec beam/spaxels of the VERTICO data correspond to a 1 200/650 pc physical resolution.

Signal masks were created for the CO(2-1) and H i datacubes as described in Brown et al. (2021, but see also Sun et al. 2018 and Leroy et al. 2021), and moment 0 maps created for each galaxy. The H i mass in each spaxel was calculated assuming optically thin emission (e.g., Meyer et al. 2017),

(1) \begin{equation}M_{\text{H}\,{\small\text{I}}}\, [\mathrm{M}_{\odot}] = 2.35\times10^5 \ D_\mathrm{L}^2(1+z)^{-2}\ S_{\text{H}\,{\small\text{I}}},\end{equation}

where $D_\mathrm{L}$ is the luminosity distance to the object (assumed to be $16.5\,$ Mpc for all VERTICO galaxies), z is the redshift, and $S_{\text{H}\,{\small\text{I}}}$ is the velocity-integrated spectral flux density (in Jy $\text{km} \, \text{s}^{-1}$ ). Unless otherwise specified, we do not correct our H i masses for the Helium abundance. The velocity-integrated CO(2-1) flux density ( $S_\mathrm{CO}$ ) was used to calculate the total CO(2-1) luminosity using the method of Solomon & Vanden Bout (2005),

(2) \begin{equation}L_\mathrm{CO} [\mathrm{K}\, \text{km} \, \text{s}^{-1}\, \mathrm{pc}^{-2}] = 3.25\times10^7\, \nu^{-2}\, D_\mathrm{L}^2\, (1+z)^{-3}\, S_\mathrm{CO},\end{equation}

where $\nu$ is the observed frequency. The molecular gas mass was then computed using

(3) \begin{equation}M_\mathrm{mol}\, [\mathrm{M}_{\odot}] = \frac{\alpha_\mathrm{CO}}{R_{21}} L_\mathrm{CO},\end{equation}

where $R_{21}$ is the integrated flux ratio of the CO(2-1) and CO(1-0) lines and $\alpha_\mathrm{CO}$ is the molecular gas mass to CO(1-0) luminosity correction factor. We used a constant $R_{21}=0.8$ , as derived by Brown et al. (2021), and a constant $\alpha_\mathrm{CO}=4.35\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ based on observations of the Milky Way disc by Bolatto et al. (2013), which includes a 36% contribution from Helium. The H i and molecular gas mass of spaxels were converted to surface densities $\Sigma_{\text{H}\,{\small\text{I}}}$ and $\Sigma_{\text{mol}}$ in units of $\mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ using the 650-pc physical pixel size.

Global stellar masses were drawn from the $z=0$ Multi-wavelength Galaxy Synthesis (z0MGS, Leroy et al. 2019) sample, a collection of (Galaxy Evolution Explorer, Martin et al. 2005), and (Wide-field Infrared Survey Explorer, Wright et al. 2010) data for $\sim$ $11\,000$ galaxies with $D<50\,$ Mpc and brighter than $M_\mathrm{B}=-18\,$ mag. Stellar-mass estimates were computed from WISE band-1 observations using the WISE band-1 and -3 colour to determine the mass-to-light ratio, as described in Leroy et al. (2019). The same prescriptions were applied on a pixel-by-pixel basis to create resolved maps of stellar surface density ( $\Sigma_{\star}$ ) at the same resolution as the H i and molecular gas data. Photometric position angles (PA) and inclinations (i) for each galaxy were calculated by fitting Kron apertures to SDSS ^c r-band images (Brown et al. 2021) and assuming an intrinsic axial ratio of $q=0.2$ , and B-band 25th $\mathrm{mag}\, \mathrm{arcsec}^{-2}$ optical radii ( $R_{25}$ ) were sourced from Chung et al. (2009). All surface densities used in this work were de-projected to face-on values using a factor of $\cos(i)$ .

To quantify the impact of the environment on VERTICO galaxies, we used the type-independent H i deficiency parameter (def-H i) for each galaxy from Chung et al. (2009), which compares the observed H i content to the average content in isolated galaxies with the same optical size. It is defined as

(4) \begin{equation} \text{def-H}\,{\small\text{I}} = \overline{\log\Sigma_{\text{H}\,{\small\text{I}}}} - \log\Sigma_{\text{H}\,{\small\text{I}}, obs}\end{equation}

where $\Sigma_{\text{H}\,{\small\text{I}}, obs} = \text{S}_{\text{H}\,{\small\text{I}}}/ (2\, R_{25})^2$ is the mean surface density within the optical disc. Following Chung et al. (2009), this definition uses $\text{S}_{\text{H}\,{\small\text{I}}}$ in units of Jy $\text{km} \, \text{s}^{-1}$ and $R_{25}$ in arcminutes, and $\overline{\log\Sigma_{\text{H}\,{\small\text{I}}}} = 0.37$ is the average type-independent value of $\log(M_{\text{H}\,{\small\text{I}}}/D_{25}^2)$ from Haynes & Giovanelli (1984, Table 4, 6.81) after a unit conversion from $\mathrm{M}_{\odot}\,\mathrm{kpc}^{-2}$ . A larger def-H i value corresponds to a smaller H i mass compared to isolated galaxies, and the typical scatter in the def-H i distribution for isolated galaxies is $\sim$ $0.2-0.3$ dex (Haynes & Giovanelli 1984; Jones et al. 2018, thus automatically encapsulating the variations in the spatial distribution of H i surface density observed in isolated galaxies). In this work, we define ‘H i-deficient’ galaxies as those with def-H i $\geq0.3 \, \mathrm{dex}$ , namely, a factor of $\gtrsim$ 2 more H i-poor than expected.

We are primarily interested in the properties of the gas within the optical disc of galaxies, so we masked each galaxy to only include pixels/spaxels within the B-band $R_{25}$ , where the radius of each spaxel was calculated from the photometric PA and inclination. While this allows us to study the impact of the environment across the discs of galaxies, we are also interested in how the environment impacts the surviving cold gas disc. The extent of the H i gas is a good tracer of the gas truncation radius, outside of which stripping is usually assumed to be effective. As our analysis uses the relationships between $\Sigma_{\star}$ , $\Sigma_{\text{H}\,{\small\text{I}}}$ , and $\Sigma_{\text{mol}}$ , it is useful to quantify the size of the H i reservoir in terms of $\Sigma_{\star}$ .

We computed radial profiles of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ and $\log\Sigma_{\star}$ using 4 arcsec annuli, which we found gave the best consistency with the Chung et al. (2009) VIVA measurements. Defining the H i truncation radius ( $R_{\text{H}\,{\small\text{I}}}$ ) to be where the radial average $\Sigma_{\text{H}\,{\small\text{I}}} = 1\,\mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ (e.g.,Wang et al. 2014, 2016; Stevens et al. 2019a), we used the value of $\log\Sigma_{\star}$ at this radius to determine the average $\log\Sigma_{\star}$ at $R_{\text{H}\,{\small\text{I}}}$ , $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ . In the most H i-rich galaxies, $R_{\text{H}\,{\small\text{I}}}$ can fall outside of the B-band $R_{25}$ , but this does not affect any of our results as these galaxies are not the focus of our analysis. Further, in all galaxies, less than 5% of spaxels with $\Sigma_{\star}\geq\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ have radius $>R_{\text{H}\,{\small\text{I}}}$ , meaning trends with higher (lower) $\Sigma_{\star}$ trace trends with smaller (larger) galactocentric radius. We also note here that while the H i observations cover the optical disc of each galaxy, the size of the VERTICO CO(2-1) observation footprints were based on the $R_{\text{H}\,{\small\text{I}}}$ of each target. Subsequently, there is some variation in the physical coverage of the molecular gas in galaxies ( $0.44 - 1\times R_{25}$ ), and we accounted for this in our analysis of the molecular gas content.

While the VERTICO sample contains some H i-normal systems, environmental processes can begin to impact galaxies before they become H i-deficient. Thus, we chose to draw our control sample from a different survey to ensure that we had an H i-normal reference that was as free from environmental effects as possible. We include the H i-normal VERTICO systems in some of our figures (Figures 1, 3, and 6) and discuss them briefly in Section 4, but the main focus of our analysis is the H i-deficient galaxies.

2.2 HERACLES control sample

We extracted a control sample of galaxies from the HERA ^d CO(2-1) Legacy Survey (HERACLES, Leroy et al. 2009). HERACLES targeted the CO(2-1) emission line in 33 galaxies from The H i Nearby Galaxy Survey (THINGS, Walter et al. 2008), a sample of nearby ( $D<15\,$ Mpc), late-type galaxies spanning from dwarfs to massive spirals. After excluding non-detections and five galaxies that overlap with VERTICO, there are 25 HERACLES galaxies with public data (for more details see Brown et al. 2021). Physical quantities for HERACLES galaxies were calculated using the same methods and sources as described for VERTICO, and B-band $R_{25}$ were drawn from Walter et al. (2008). We did not compute $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ for HERACLES galaxies as the $R_{\text{H}\,{\small\text{I}}}$ was nearly always outside $R_{25}$ .

The HERACLES native resolution datacubes have a synthesised beam-width of 13 arcsec, while the THINGS data are more varied depending on the VLA configuration. Before calculating physical quantities, the HERACLES and THINGS data were smoothed using the imsmooth routine of CASA to match the physical spaxel size and resolution of the VERTICO data, and circularise the synthesised beam. As the distances to HERACLES sample galaxies vary, this left 11 galaxies that could be matched to a 1 200-pc physical resolution in both CO(2-1) and H i.

To ensure that our HERACLES control sample is a reasonable representation of star-forming galaxies, we restricted this sample to galaxies more star-forming than $1\sigma$ below the star-forming main-sequence using the parametrisation from the z0MGS sample ( $\sigma=0.36$ , with star-formation rates (SFR) from the same catalogue, Leroy et al. 2019). After the stellar mass and inclination selection cuts, our final HERACLES control sample consists of eight galaxies, listed in Table 1. We compared the integrated properties of the final control sample galaxies to a subset of star-forming galaxies from the xCOLD GASS ^e (Saintonge et al. 2017) sample, selected using the same SFR threshold. While the HERACLES galaxies lie on or slightly above ( $<1\sigma$ ) the locus of the star-forming main-sequence, they span the same range in $\log\,M_{\text{H}\,{\small\text{I}}}/M_{\star}$ and $\log\,M_{\mathrm{MG}}/M_{\star}$ as the xCOLD GASS sample, indicating that there is no bias towards H i- or molecule-rich systems in our control sample.

Figure 1. Spatially resolved $\Sigma_{\text{H}\,{\small\text{I}}}-$ and $\Sigma_{\text{mol}}-\Sigma_{\star}$ scaling relations for the combined samples. The left panel shows $\log\Sigma_{\text{H}\,{\small\text{I}}}$ vs $\log\Sigma_{\star}$ , the resolved H i–stellar mass scaling relation, and the right panel shows $\log\Sigma_{\text{H}\,{\small\text{I}}}$ vs $\log\Sigma_{\star}$ , the resolved molecular gas–stellar mass scaling relation. Contours for both samples enclose 5, 16, 50, 84, and 95% of the H i or molecular gas detected spaxels, while non-detected spaxels are shown with downward triangles at $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=\log(\Sigma_{\text{mol}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=-0.5$ . Density-normalised histograms of $\log\Sigma_{\star}$ , $\log\Sigma_{\text{H}\,{\small\text{I}}}$ , and $\log\Sigma_{\text{mol}}$ are shown opposite their respective axes, with filled histograms representing H i- or molecular gas-detected spaxels while open histograms represent non-detections. The HERACLES $\Sigma_{\star}$ histogram for H i non-detections uses 0.25 dex bins compared to the 0.1 dex used for VERTICO as it has fewer points (331 compared to 5 717). Compared to HERACLES, the combined VERTICO sample has fewer $\Sigma_{\text{H}\,{\small\text{I}}}$ detected spaxels, particularly at lower $\log\Sigma_{\star}$ , and a lower average $\log\Sigma_{\text{H}\,{\small\text{I}}}$ . Similarly, but less clearly, the combined VERTICO sample has fewer $\Sigma_{\text{mol}}$ detections at lower $\Sigma_{\star}$ , while the histogram of detected $\Sigma_{\text{mol}}$ has a larger tail towards lower values.

Last, we briefly mention the detection limits of each dataset after the homogenisation of the spatial resolutions and the masking procedures. The HERACLES:THINGS data are shallower than the VERTICO:VIVA data in both H i ( $\Sigma_{\text{H}\,{\small\text{I}}}\geq1-3\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ within B-band $R_{25}$ compared to $\Sigma_{\text{H}\,{\small\text{I}}}\geq1\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ ) and molecular gas ( $\Sigma_{\text{mol}}\geq 1-3\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ compared to $\Sigma_{\text{mol}}\geq0.4-1\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ ). Given that we are interested in environment-driven reductions in both H i and molecular gas surface densities, the fact that the VERTICO:VIVA data are deeper ensures that our findings are not due to differences in sensitivity between the two datasets. The shallower nature of HERACLES means that it can have an increased fraction of non-detected spaxels near the edge of the disc, and we account for this in our analysis.

3. Combined $\Sigma_{\text{H}\,{\small\text{I}}}$ – $\Sigma_{\star}$ and $\Sigma_{\text{mol}}$ – $\Sigma_{\star}$ relations

In the left panel of Figure 1 we plot $\Sigma_{\text{H}\,{\small\text{I}}}$ against $\Sigma_{\star}$ , in other words, the resolved H i–stellar mass (rH i SM) relation. Detected spaxels for HERACLES (blue) and VERTICO (red) are shown with contours enclosing 5, 16, 50, 84, and 95% of the population, and non-detections are shown at $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=-0.5$ with downward triangles. Density-normalised histograms of $\log\Sigma_{\star}$ and $\log\Sigma_{\text{H}\,{\small\text{I}}}$ are shown opposite to their respective axes, computed with 0.25 dex bins. Filled histograms show detected spaxels, while the open histograms on the $\Sigma_{\star}$ axis show the histogram of $\log\Sigma_{\star}$ for H i non-detected spaxels.

Focusing on HERACLES galaxies, the rH i SM relation of H i-normal systems occupies a $\sim$ 0.75 dex range in $\Sigma_{\text{H}\,{\small\text{I}}}$ across the $\sim$ 2.5 dex in $\Sigma_{\star}$ . H i is detected across the whole $\Sigma_{\star}$ range, with non-detected spaxels having preferentially lower $\Sigma_{\star}$ compared to detections, also shown in the $\Sigma_{\star}$ histograms. This distribution of HERACLES points reflects the known result that H i surface densities typically saturate around $\sim10\, \mathrm{M}_{\odot}\, \mathrm{pc}^{-2}$ within the stellar components of galaxies (Martin & Kennicutt 2001; Krumholz, McKee, & Tumlinson 2009; Bigiel et al. 2008; Bigiel & Blitz 2012). Focusing now on VERTICO, it is clear that the VERTICO combined rH i SM relation inhabits a different region of the parameter space compared to HERACLES. First, there is a decreased number of H i-detected spaxels across the range of $\Sigma_{\star}$ , but preferentially at lower $\Sigma_{\star}$ . This change is also visible in the $\Sigma_{\star}$ histograms where both the average VERTICO detected and non-detected spaxels are shifted to higher $\Sigma_{\star}$ compared to HERACLES. Second, the typical $\Sigma_{\text{H}\,{\small\text{I}}}$ of VERTICO galaxies is lower. It is missing the largest $\Sigma_{\text{H}\,{\small\text{I}}}$ values, there is a cloud of red points at $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})>1.5$ with $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2}) \sim 0-0.5$ that does not exist in HERACLES, and the $\log\Sigma_{\text{H}\,{\small\text{I}}}$ histogram of VERTICO galaxies is broader, extending to lower $\Sigma_{\text{H}\,{\small\text{I}}}$ . The median detected spaxel has $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.58$ for VERTICO galaxies, while HERACLES has $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.93$ , and a two-sample Kolmogorov–Smirnov (K-S) test rejects the null hypothesis that the two $\Sigma_{\text{H}\,{\small\text{I}}}$ distributions are drawn from the same parent distribution with high significance ( $p<10^{-4})$ .

Given that the VERTICO sample spans a range of def-H i values ( $-0.3-2.25\, \mathrm{dex}$ ) and $\Sigma_{\star}$ increases towards the centres of galaxies, these changes are broadly consistent with galaxies undergoing ram-pressure stripping from the ICM. Namely, the ISM at large radii is affected first, while lowering the average gas surface density within the central regions where the gas is more gravitationally bound requires stronger stripping that has proceeded within the stellar disc (Roediger & Brüggen 2007; Fumagalli et al. 2009; Steinhauser, Schindler, & Springel 2016; Yoon et al. 2017), and these two effects clearly emerge in this figure.

How are these differences reflected in the molecular gas surface density of galaxies? In the right panel of Figure 1, we show the resolved molecular gas–stellar mass (rMGSM) relation ^f for the combined HERACLES and VERTICO samples using the same colouring and non-detection labelling as the rH i SM relation. Unlike $\Sigma_{\text{H}\,{\small\text{I}}}$ , $\Sigma_{\text{mol}}$ shows a positive correlation with $\Sigma_{\star}$ , and differences between the combined rMGSM relation of HERACLES and VERTICO are less distinct. Non-detected spaxels in both samples exist preferentially at lower $\Sigma_{\star}$ , although the typical $\Sigma_{\star}$ of VERTICO non-detected spaxels is larger than HERACLES.

The clearest difference, however, is that the VERTICO contours extend to lower $\Sigma_{\text{mol}}$ at all $\Sigma_{\star}$ , and this is true even at higher $\Sigma_{\star}$ where all $\Sigma_{\text{mol}}$ are above the detection limit. This difference is also visible in the VERTICO $\Sigma_{\text{mol}}$ histogram, which extends to lower $\Sigma_{\text{mol}}$ . Quantitatively, the median detected VERTICO spaxel has $\log(\Sigma_{\text{mol}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.73$ , lower than the HERACLES value of $\log(\Sigma_{\text{mol}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.85$ , and we can exclude that the distributions are drawn from the same parent sample ( $p<10^{-4}$ ). Some of these differences will be due to the shallower HERACLES detection limit (causing more non-detections at lower $\Sigma_{\star}$ ) and the smaller spatial coverage of VERTICO galaxies with smaller H i discs (causing an under-sampling of lower $\Sigma_{\star}$ ). However, they also suggest that the molecular gas is being lost preferentially at lower $\Sigma_{\star}$ (at least to our detection threshold), while the remaining, detectable molecular gas has lower surface density. However, the magnitude of these effects appears to be much smaller, and the affected region(s) are less well defined compared to the rH iSM relation.

Figure 2. Resolved $\Sigma_{\text{H}\,{\small\text{I}}}-$ and $\Sigma_{\text{mol}}-\Sigma_{\star}$ relations for HERACLES control galaxies. Left two panels: The combined rH i SM relation (top) and rMGSM (bottom) are shown with grey points in the background, the median relation in bins of $\log\Sigma_{\star}$ for each galaxy is shown with thick, coloured lines, and the 16th– $84\mathrm{th}$ percentiles are shaded in the same colour, which matches the coloured name of each galaxy in the top left corner. Galaxies are coloured from red through to black increasing integrated molecular gas mass and alternate between panels for better visual separation. Third panel: The combined rH i SM relation (top) and rMGSM (bottom) are shown with grey points in the background, overlaid are the medians (solid lines) and 16th– $84\mathrm{th}$ percentiles (dashed lines) treating all points equally (unweighted, blue) and each galaxy equally (weighted, red). The rightmost panels show the density-normalised histograms for the unweighted and weighted treatments. The key message of this Figure is that galaxies lie along well-defined, individual sequences within the rH i SM and rMGSM relations. As the location of an individual sequence is best described by global galaxy properties, we choose to weight each galaxy equally. HERACLES galaxies do not show a large difference between the weighted and unweighted statistics, as they are restricted to the star-forming regime.

How, then, do we extract the comparatively weaker signal (if any) of the impact of environmental processes on the rMGSM relation? Using a sample of galaxies from the ALMaQUEST survey (ALMA-MaNGA QUEnching and STar Formation, Lin et al. 2020) with similar $M_{\star}$ and SFR ranges to our VERTICO sample, Ellison et al. (2021a,b) showed that, to first order, the scatter within the rMGSM relation is not determined by the properties of individual spaxels. Instead, galaxies form well-defined individual sequences in the parameter space, and their scatter from the median relation is best described with global galaxy properties (e.g., the location of an individual rMGSM is best described by specific SFR, Ellison et al. 2021a). We show these individual sequences for both the rH i SM and rMGSM relations of our HERACLES control sample in the left panel of Figure 2, ordered from red to black by increasing total molecular gas mass. The distribution of HERACLES points is shown in the background, with the median relation (thick coloured lines) and $16\mathrm{th}$ – $84\mathrm{th}$ percentiles (coloured, shaded regions) for each galaxy in bins of $\log\Sigma_{\star}$ overlaid. Clearly, the galaxy–galaxy variation is larger than the scatter within individual objects, and this is particularly evident in the higher $\Sigma_{\star}$ regions of the rH iSM relation, where galaxies might have a central H i depression.

As the VERTICO sample contains both H i-normal and H i-deficient galaxies, combining the entire sample into one relation acts to wash out the signature of environmental effects, particularly in the molecular gas. Thus, to understand how environmental processes impact the resolved properties of galaxies, we need to control for their global properties. Further, we must also ensure that any histograms and statistics we compute are not biased towards larger galaxies and galaxies that are less affected by the environment, as these systems contain more spaxels. To do this, whenever we compute a histogram such as $\log\Sigma_{\text{H}\,{\small\text{I}}}$ or $\log\Sigma_{\text{mol}}$ we weight each spaxel by $1/(N_\mathrm{sp}N_\mathrm{G})$ , where $N_\mathrm{sp}$ is the number of spaxels contributed by its host galaxy and $N_\mathrm{G}$ the number of contributing galaxies. Under this definition, the sum of the weights is always 1, and it ensures that each galaxy contributes to the histogram (and any subsequently computed statistics) equally. All statistics presented in this paper from now on are weighted values.

In the middle panels of Figure 2 we compare the weighted (red) and unweighted (blue) median (solid line), and 16th and 84th percentiles (dashed lines),^g of the rH i SM relation and rMGSM relation for the HERACLES control sample. The rightmost panels show the weighted and unweighted histograms of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ and $\log\Sigma_{\text{mol}}$ across all $\Sigma_{\star}$ . The advantage of this weighting is clear in $\Sigma_{\text{mol}}$ at $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})\lesssim 1.5$ , where the HERACLES data begin to reach their detection threshold. While the unweighted statistics become non-detections, the weighted statistics better trace the relation as the typical galaxy is still detected. However, to avoid biasing our statistics, we will restrict some of our quantitative $\Sigma_{\text{mol}}$ comparisons to $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})\geq1.25$ , as this is the typical $\Sigma_{\star}$ where the HERACLES galaxies transition to being dominated by non-detections. At these higher $\Sigma_{\star}$ , and at all $\Sigma_{\star}$ for $\Sigma_{\text{H}\,{\small\text{I}}}$ , there are negligible differences between the two weighting methods. This is expected as our control sample is defined to be star-forming systems, which exhibit less galaxy-to-galaxy variation in their spatially resolved cold gas scaling relations (Ellison et al. 2021a). As we see below, this variation is much larger for VERTICO galaxies.

Figure 3. Individual rH i SM relations for VERTICO galaxies, with the galaxy name given in the top-right of each panel. The panels are ordered in rows by def-H i (noted in the top-left corner), which also sets the colour of each panel, and a colour bar is displayed on the right-hand side. Solid points indicate detected spaxels, while downward triangles at $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=-0.5$ indicate non-detections. The median rH i SM for each galaxy computed in 0.25 dex $\Sigma_{\star}$ bins is shown with large white markers, where downward triangles indicate the median is a non-detection, and error bars show the 16th and 84th percentiles of the $\log\Sigma_{\text{H}\,{\small\text{I}}}$ distribution in each bin. The HERACLES control sample is shown with thick (median) and thin (16th and 84th percentiles), and the vertical dashed, coloured lines show $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ , the average $\log\Sigma_{\star}$ at the H i truncation radius ( $R_{\text{H}\,{\small\text{I}}}$ ). The two galaxies marked with stars are interacting, and their H i observations are partially confused. The impact of the Virgo cluster on the rH i SM relation is two-fold. As def-H i increases there are more non-detected spaxels and $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ moves to larger values; and there is an increasing reduction of the typical $\Sigma_{\text{H}\,{\small\text{I}}}$ within $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ .

4. Impact of environment on the spatially resolved cold gas properties of galaxies

Environmental effects on the H i content of galaxies are a well-studied topic, though many previous works have focused on the gas at large radii due to the extended nature of H i. In Section 4.1, we present an analysis of the H i content in the central regions of galaxies using the framework of the spatially resolved scaling relations. Not only does this present a new view of how environmental processes impact the resolved H i content of galaxies, but it is also useful to outline the visualisations and techniques that we apply to the molecular gas in Section 4.2.

4.1 Atomic hydrogen

In Figure 3, we present the rH i SM relations of VERTICO galaxies in individual panels. H i-detected spaxels are shown with coloured points and non-detections with downward triangles. The median $\log\Sigma_{\text{H}\,{\small\text{I}}}$ in $0.25$ dex bins of $\Sigma_{\star}$ (which are combined with adjacent bins if they contain less than 10 points) is shown with white, bordered symbols, where a circle denotes a median detection and a downward triangle denotes that the median is a non-detection. The median (thick, grey line) and 16th and 84th percentiles (thin, grey lines) of the HERACLES control sample are shown in the background, and the average $\log\Sigma_{\star}$ measured at $R_{\text{H}\,{\small\text{I}}}$ ( $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ , as calculated in Section 2.1) is shown with a vertical dashed line. Panels are ordered in rows from left to right by increasing def-H i (the numerical value is displayed in the top left of each panel, in-line with the galaxy name), which also determines the colour of the individual points (with a colour bar on the right). Figure 3 is rich in physical information, and we dissect this in the following analysis.

We first focus on the top row, which contains VERTICO galaxies considered to be H i-normal based on their def-H i ( $-0.3<\text{def-H}\,{\small\text{I}}\ [dex] <0.3$ ). Three galaxies are consistent with the HERACLES control sample, while NGC 4567 appears to be a slight outlier, particularly at lower $\Sigma_{\star}$ . NGC 4567 is interacting with NGC 4568 (row 2, column 2) and the galaxies cannot be separated in velocity space (Chung et al. 2009), causing an overlap of the rH i SM relation of each galaxy and meaning that it has experienced some environmental mechanisms.

Table 2. Weighted median $\log\Sigma_{\text{H}\,{\small\text{I}}}$ (Figure 4) and $\log\Sigma_{\text{mol}}$ (Figure 7) for the control sample, and the moderate and large def-H i samples, in bins of $\log\Sigma_{\star}$ . The $\Sigma_{\star}$ range used to select spaxels are given in the first column, the second to fourth columns use all spaxels, while the remaining columns use only detected spaxels for each gas phase. Uncertainties are calculated with bootstrap re-sampling using $10^4$ iterations where entire galaxies are selected in each subsample to maintain equal weighting.

^aMedian is a non-detection.

^bMedian consistent with control sample within uncertainty.

The rH iSM relations of individual galaxies show two clear trends as a function of increasing def-H i that occur distinctly at lower and higher $\Sigma_{\star}$ . First, the value of $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ (measuring $R_{\text{H}\,{\small\text{I}}}$ , vertical dashed lines), on average, increases with increasing def-H i. In other words, the H i is stripped (at least to our sensitivity limit) to smaller radii within the disc as def-H i increases. Second, inside $R_{\text{H}\,{\small\text{I}}}$ , the individual rH i SM relations increasingly fall below the median and then the 16th (lower) percentile of the HERACLES control sample at fixed $\Sigma_{\star}$ . That is, there is an increasing reduction of the typical $\Sigma_{\text{H}\,{\small\text{I}}}$ within the $R_{\text{H}\,{\small\text{I}}}$ of galaxies as def-H i increases. This effect is largest in the highest def-H i systems, but several moderate def-H i galaxies also begin to show evidence for this reduction (NGC 4321, NGC 4501, NGC 4535). These trends show that the impact of environmental processes on the H i is two-fold: not only is the H i stripped from the outskirts, but the typical surface density of the remaining gas is also reduced. Additionally, we observe no evidence for an environment-driven enhancement of $\Sigma_{\text{H}\,{\small\text{I}}}$ , at least at the 1.2-kpc spatial scale that we trace.

To quantify these effects, and because this is a complex parameter space, we use two different strategies for computing and comparing $\log\Sigma_{\text{H}\,{\small\text{I}}}$ histograms of VERTICO galaxies to the control sample. In all comparisons, we split the VERTICO galaxies into two bins of global def-H i,^h def-H i $=[0.3,0.7)\, \mathrm{dex}$ and def-H i $\geq0.7 \, \mathrm{dex}$ , which we refer to these as the ‘moderate’ and ‘large’ def-H i samples, respectively. These thresholds were defined as approximately 1 and 2 times the scatter in the def-H i of isolated galaxies ( $\sigma$ =0.3 dex), with a slight adjustment to ensure an equal number of galaxies in each sample. Histograms are computed by weighting each galaxy equally (see Section 3), and we also computed weighted medians for each def-H i sample and the control sample, which are presented in Table 2.

In our first comparison, shown in Figure 4, we divide the rH i SM relation into 0.5 dex bins of $\log\Sigma_{\star}$ and compute the histogram of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ in each bin. The control sample is shown as a grey, solid histogram, the moderate def-H i population as a blue, open histogram, and the large def-H i as a red, open histogram. All H i non-detected spaxels are included in the bin located at $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=-0.5$ , and if the value of this bin exceeds the figure limits, its value is written in the top left.

Figure 4. $\log\Sigma_{\text{H}\,{\small\text{I}}}$ histograms in bins of $\log\Sigma_{\star}$ . Each panel shows density-normalised histograms of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ for spaxels in the $\Sigma_{\star}$ range quoted in the top right corner. HERACLES control sample galaxies are shown with a grey, filled histogram, VERTICO moderate def-H i galaxies with a blue, open histogram, and large def-H i galaxies with a red, open histogram. Non-detected spaxels are included in the bin centred at $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=-0.5$ , and if the value of this bin exceeds the y-axis limit then its value is quoted in the top-left of the panel in the same colour as its histogram. At fixed $\Sigma_{\star}$ , galaxies with larger def-H i have a higher fraction of non-detected spaxels and their detected spaxels have lower $\Sigma_{\text{H}\,{\small\text{I}}}$ . The difference is larger at lower $\Sigma_{\star}$ .

Figure 4 again highlights the trends seen in Figure 3. First, both def-H i samples show a larger fraction of non-detected H i spaxels than the control sample, and the fraction of non-detected spaxels is higher in the larger def-H i galaxies. Namely, the radius at which H i is below the detection limit gets closer to the galaxy centre as def-H i increases. This increased fraction of non-detections is also evident in Table 2, where the weighted median for detections only is often different to the value for all points in both def-H i samples. Second, at all $\Sigma_{\star}$ , the typical detected $\Sigma_{\text{H}\,{\small\text{I}}}$ is also lower as def-H i increases. Table 2 shows that the moderate def-H i sample has a weighted median for detected $\Sigma_{\text{H}\,{\small\text{I}}}$ spaxels that is 0.34–0.83 dex lower than the control sample, while it is 0.53–0.85 dex lower for the large def-H i sample, and these differences are greater than the estimated uncertainties. Last, neither def-H i sample shows higher $\Sigma_{\text{H}\,{\small\text{I}}}$ than the control sample, highlighting that we observe no environment-driven enhancement of $\Sigma_{\text{H}\,{\small\text{I}}}$ .

We note that there is some overlap between the 1 $\sigma$ scatter in the $\Sigma_{\text{H}\,{\small\text{I}}}$ distributions. For example, in the highest $\Sigma_{\star}$ bin the 16th percentile of the control sample (0.42) is smaller than the 84th percentiles of both the moderate (0.73) and large (0.49) def-H i samples. As we do not expect environmental mechanisms to operate across the whole gas disc with the same magnitude (e.g., the centre or the trailing edge will be less affected), some overlap between the distributions is expected. Nevertheless, in all $\Sigma_{\star}$ bins we find that the two distributions cannot be drawn from the same sample ⁱ ( $p<10^{-4}$ ), confirming that we are witnessing the effect of environmental processes already in the moderate def-H i regime.

Considering the K-S test results, Figure 4 also contains interesting insights into how environmental processes change the distribution of $\Sigma_{\text{H}\,{\small\text{I}}}$ in the remaining gas. While the bins with $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})<1.25$ are dominated by non-detections, in the higher $\log\Sigma_{\star}$ bins, the $\log\Sigma_{\text{H}\,{\small\text{I}}}$ histogram of the moderate def-H i galaxies (blue) is similar to the control sample. Namely, there is a peak at higher $\Sigma_{\text{H}\,{\small\text{I}}}$ and a tail towards lower surface densities. The main difference is a decrease in the average $\Sigma_{\text{H}\,{\small\text{I}}}$ , while the increase in non-detections comes from the low surface density tail falling below the detection limit. The large def-H i galaxies show a larger decrease in $\Sigma_{\text{H}\,{\small\text{I}}}$ values, but also, the shape of the $\log\Sigma_{\text{H}\,{\small\text{I}}}$ histogram becomes more gaussian rather than skewed. This change in shape suggests that in addition to the decrease in the typical surface density, the densest H i regions are no longer present.

In our second comparison, we consider only H i within the $R_{\text{H}\,{\small\text{I}}}$ of each galaxy. For both def-H i samples, we used all spaxels with $\log\Sigma_{\star} \geq \langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ to compute equally weighted density-normalised $\log\Sigma_{\text{H}\,{\small\text{I}}}$ histograms. Achieving a fair comparison with the control sample required an extra step, however, as the $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ of each galaxy is different. For each VERTICO galaxy, we computed a unique, equally weighted control sample histogram using all spaxels from each control galaxy with $\log\Sigma_{\star} \geq \langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ . The control histograms for each galaxy in a given def-H i sample were then stacked and re-normalised using the same procedure.^j This dependence on $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ results in slightly different control samples for each def-H i sample^k and is visible when comparing the grey histograms in the two panels of Figure 5. Weighted medians for these histograms are listed in Table 3.

Figure 5. Histograms of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ within the $R_{\text{H}\,{\small\text{I}}}$ of galaxies. We compare density-normalised histograms of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ for spaxels with $\log\Sigma_{\star}\geq \langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ for galaxies with moderate def-H i (top panel; blue, open histogram) and large def-H i (bottom panel; red, open histogram) to their matched HERACLES control samples (both panels, grey filled histogram). As def-H i increases, the fraction of non-detected spaxels within $R_{\text{H}\,{\small\text{I}}}$ increases, and the typical detected $\Sigma_{\text{H}\,{\small\text{I}}}$ is lower. In galaxies with the largest def-H i the histogram of detected $\log\Sigma_{\text{H}\,{\small\text{I}}}$ is gaussian, rather than skewed, indicating that the densest H i is no longer present. Neither def-H i sample shows evidence for elevated $\Sigma_{\text{H}\,{\small\text{I}}}$ compared to the control sample.

In Figure 5, we compare histograms of $\log\Sigma_{\text{H}\,{\small\text{I}}}$ for spaxels within $R_{\text{H}\,{\small\text{I}}}$ between the control sample (grey, solid), the moderate def-H i sample (top panel, blue) and large def-H i sample (bottom panel, red). Here, we can see the impact of environmental processes on the remaining H i gas that has not been removed. The moderate def-H i galaxies have median detected $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.52\pm0.07$ , while their matched control sample has $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.95\pm0.02$ , indicating that the typical $\Sigma_{\text{H}\,{\small\text{I}}}$ is 0.43 dex lower (Table 3). The difference is larger in the large def-H i galaxies, which have a median $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.28\pm0.06$ that is 0.66 dex lower than the median of their matched control sample, $\log(\Sigma_{\text{H}\,{\small\text{I}}}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=0.94\pm0.02$ . We also see the same changes in the shape of the $\Sigma_{\text{H}\,{\small\text{I}}}$ histogram described above in Figure 4, and we can rule out the control and def-H i samples being drawn from the same parent distribution ( $p<10^{-4}$ ). Thus, when considering the remaining gas inside $R_{\text{H}\,{\small\text{I}}}$ , environmental processes typically lower the average $\Sigma_{\text{H}\,{\small\text{I}}}$ of moderate def-H i galaxies, while large def-H i galaxies can no longer support the densest H i gas.

Our results so far have provided a unique view of how environmental mechanisms impact the resolved H i content of galaxies, and agree with results from previous studies (e.g. Vollmer et al. 2001; Tonnesen & Bryan 2009; Chung et al. 2009; Yoon et al. 2017; Cortese et al. 2021; Boselli et al. 2022; Lee et al. 2022). However, how these mechanisms impact the molecular gas in galaxies on these spatial scales has so far remained unquantified, and is discussed next.

4.2 Molecular hydrogen

In Figure 6, we present the individual rMGSM relations of VERTICO galaxies in the same way as the rH i SM relations. Focusing first on the top row of H i-normal galaxies, aside from NGC 4567 (interacting, as discussed in Section 4.1), NGC 4651 seems to be molecule-poor for an H i-normal galaxy. The kinematic position angles of the H i inside and outside the optical disc of NGC 4651 are different (Chung et al. 2009), and it has stellar and H i tidal features (Chung et al. 2009; Martínez-Delgado et al. 2010; Morales et al. 2018), suggesting that much of its H i disc could have been recently accreted. Thus, the molecular gas content likely does not reflect that of an H i-normal galaxy. NGC 4254 appears slightly molecule-rich at fixed $\Sigma_{\star}$ . There is evidence that it recently underwent a fly-by interaction (Vollmer et al. 2005; Haynes, Giovanelli, & Kent 2007; Boselli et al. 2018b), and it has a molecular bar and asymmetric spiral arms containing dense molecular gas (Sofue et al. 2003). These differences between the control sample and NGC 4651 and NGC 4254 highlight how environmental mechanisms can influence galaxies before they become significantly H i-poor, and the importance of selecting a control sample from field galaxies.

Figure 6. Individual rMGSM relations for VERTICO galaxies, presented in the same way as Figure 3. As def-H i increases, the clearest deviation of $\Sigma_{\text{mol}}$ from the HERACLES control sample occurs at lower $\Sigma_{\star}$ , typically close to or outside of $R_{\text{H}\,{\small\text{I}}}$ ( $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ ). Only some of the galaxies with the largest def-H i have reduced $\Sigma_{\text{mol}}$ at higher $\Sigma_{\star}$ and within their $R_{\text{H}\,{\small\text{I}}}$ (e.g., NGC 4579 NGC 4380, NGC 4450).

Before dissecting Figure 6, it is useful to point out some interesting features. First, the median rMGSM relation of four galaxies (NGC 4298, NGC 4689, NGC 4457, NGC 4450) transition from detections to non-detections coincidently with $R_{\text{H}\,{\small\text{I}}}$ , $\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ , indicating that these galaxies show co-spatial H i and molecular gas truncation. Second, several galaxies show changes in the slope of their rMGSM relation. Some of these changes are close to $R_{\text{H}\,{\small\text{I}}}$ and are likely the signature of the environment (e.g. NGC 4298, NGC 4689). Other galaxies show changes in their more central regions, such as the superposition of two slopes in NGC 4535, a flattening and then upturn in NGC 4548, or the zig-zag structure of NGC 4579. These are not solely due to changes in the $\Sigma_{\text{mol}}$ content but also reflect changes in the central stellar structure of galaxies. NGC 4535 and NGC 4548 are barred with central bulges, causing a steeper $\Sigma_{\star} - \Sigma_{\text{mol}}$ relation in their centres and bar regions and a flatter relation in the disc. NGC 4579 has a bar and a ring, causing an alternation between elevated $\Sigma_{\text{mol}}$ and $\Sigma_{\star}$ in the centre, lowered $\Sigma_{\text{mol}}$ in a higher $\Sigma_{\star}$ region between the centre and the ring, and higher $\Sigma_{\text{mol}}$ in a lower $\Sigma_{\star}$ region in the ring and disc. Beyond having interesting effects on the shape of the rMGSM relation, this suggests that we must use caution when comparing $\Sigma_{\text{mol}}$ properties at fixed $\Sigma_{\star}$ , particularly in the ranges where stellar structures could be important.

Table 3. Weighted median $\log\Sigma_{\text{H}\,{\small\text{I}}}$ (Figure 5) and $\log\Sigma_{\text{mol}}$ (Figure 8) within $R_{\text{H}\,{\small\text{I}}}$ for the moderate and large def-H i samples, and their matched control samples. As the control sample histogram is constructed by selecting spaxels with $\log\Sigma_{\star}\geq\langle\log\Sigma_{\star}\rangle_{\text{H}\,{\small\text{I}}}$ in each galaxy, the two def-H i samples have unique control samples, and thus unique weighted control sample medians. The first column gives the def-H i sample, the second and third columns the control sample and def-H i sample medians using all spaxels, and the remaining columns use detected spaxels only. Uncertainties are calculated using jackknife re-sampling and rejecting an entire galaxy in each iteration.

Comparing the rMGSM relations of VERTICO galaxies to the control sample, there is less difference as a function of def-H i than was visible for the rH i SM relation. In the second and third rows, namely at moderate def-H i, the rMGSM relations of individual galaxies typically agree well with the median of the control sample, though there are some small signatures of deviation near $R_{\text{H}\,{\small\text{I}}}$ in some galaxies. Instead, it is primarily the largest def-H i systems (the bottom two rows) where the individual rMGSM relations begin to deviate below the control sample median and occasionally below the 16th percentile, at fixed $\Sigma_{\star}$ . These differences are largest at lower $\Sigma_{\text{mol}}$ and closer to $R_{\text{H}\,{\small\text{I}}}$ , though there is galaxy-galaxy variation and at high $\Sigma_{\star}$ some galaxies are consistent with the control sample (e.g. NGC 4457, NGC 4569), while others show reduced $\Sigma_{\text{mol}}$ (e.g. NGC 4380, NGC 4450).

These results imply that the $\Sigma_{\text{mol}}$ at lower $\Sigma_{\star}$ (outer parts of the disc) is more strongly impacted than at higher $\Sigma_{\star}$ (the inner regions). Further, as $\Sigma_{\text{mol}}$ is typically lower at lower $\Sigma_{\star}$ , the environmental processes preferentially impact the lower surface density molecular gas within galaxies. Interestingly, this can begin even at moderate def-H i (e.g. NGC 4298 and its $\Sigma_{\text{mol}}$ slope-change and $\Sigma_{\text{mol}}$ truncation) and effectively changes the slope of the individual rMGSM relations, making them steeper. These results agree with previous studies, which found steeper radial profiles of $\Sigma_{\text{mol}}$ in galaxies with more disturbed H i morphology and larger def-H i (Fumagalli et al. 2009; Zabel et al. 2022; Villanueva 2022). Additionally, like $\Sigma_{\text{H}\,{\small\text{I}}}$ we observe no environment-driven elevation of $\Sigma_{\text{mol}}$ at the spatial scale that we are sensitive to.

In Figure 7 we compare the $\log\Sigma_{\text{mol}}$ histograms of VERTICO galaxies to the control sample in 0.5 dex $\log\Sigma_{\star}$ bins. We do not expect stellar structures, as mentioned above, to cause any biases in these histograms as they are primarily contained in the $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})\geq2.25$ bin. These histograms reflect the smaller offsets between the individual VERTICO rMGSM relations and the control sample.

Moderate def-H i galaxies (blue) show little differences from the control sample, except in the lowest $\Sigma_{\star}$ bin that has a significant contribution from non-detections. In this bin, the relatively larger number of non-detected spaxels in the control sample is due to the shallower detection limit of HERACLES. Overall, the main difference between the control sample and moderate def-H i sample is the presence of a slightly more extended tail towards lower $\Sigma_{\text{mol}}$ . This tail imprints on the median detected $\Sigma_{\text{mol}}$ being up to $\sim$ 0.18 dex lower than the control sample, but the distributions peak in similar locations and in the two higher $\Sigma_{\star}$ bins the medians are consistent within their estimated uncertainties (Table 2). The presence of this tail, however, is sufficient to rule out the control sample and moderate def-H i sample being drawn from the same distribution in all $\Sigma_{\star}$ bins ( $p<4\times10^{-4}$ ).

In the large def-H i sample (red), the impact of environmental processes is stronger. Non-detections dominate the smallest $\Sigma_{\star}$ bin, and despite the VERTICO data being deeper, the relative number of non-detected spaxels is larger than the control sample in this bin. Considering the detected spaxels, the weighted median $\log\Sigma_{\text{mol}}$ is $0.34-0.59$ dex lower than the control sample across all $\log\Sigma_{\star}$ bins, all of which are significant considering their uncertainties (Table 2), and we can exclude that the distributions are drawn from the same parent sample ( $p<10^{-4}$ ).

These results support the interpretation that it is primarily the lower surface density molecular gas at lower $\Sigma_{\star}$ that is affected by the environment, although as galaxies become significantly H i-deficient ( $\gtrsim$ $0.7$ ), the impact is larger and can also affect the gas at higher $\Sigma_{\star}$ . Further, in the histograms of Figure 7 we see no evidence for any environment-driven elevation of $\Sigma_{\text{mol}}$ , although in the largest $\Sigma_{\star}$ bin the moderate def-H i sample appears to show slightly higher $\Sigma_{\text{mol}}$ . However, the medians are consistent within the uncertainties, and this enhancement disappears if we instead consider spaxels with $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})\geq2$ or $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})=[2,2.5)$ , meaning that it is due to the slightly different $\Sigma_{\star}$ coverage of galaxies in this bin.

In Figure 8, we compare the histograms of $\log\Sigma_{\text{mol}}$ within $R_{\text{H}\,{\small\text{I}}}$ for moderate def-H i galaxies (top, blue) and large def-H i galaxies (bottom, red) to their respective control sample $\log\Sigma_{\text{mol}}$ histograms^l (grey, solid). The $\Sigma_{\text{mol}}$ histogram of moderate def-H i galaxies is similar to the control sample, and the median detected $\Sigma_{\text{mol}}$ is $0.15$ dex lower (a factor of $\sim$ 1.4, Table 3), which is significant considering the uncertainty. Much like Figure 7, the main difference is the tail towards lower $\Sigma_{\text{mol}}$ values, and there is also a reduction of spaxels with surface densities near the peak of the histogram. The control sample has slightly more non-detected spaxels due to the shallower detection limit, and we verified that the result remains if consider only spaxels with $\log(\Sigma_{\star}/\mathrm{M}_{\odot}\, \mathrm{pc}^{-2})\geq1.5$ . We also note that the tail is primarily driven by NGC 4394, which sits below the control sample rMGSM relation at all $\Sigma_{\star}$ . We re-computed the $\log\Sigma_{\text{mol}}$ histogram with NGC 4394 excluded and found that a smaller tail towards lower surface density $\Sigma_{\text{mol}}$ still remains.

Figure 7. Histograms of $\log\Sigma_{\text{mol}}$ in bins of $\Sigma_{\star}$ , presented in the same way as Figure 4. At fixed $\Sigma_{\star}$ , moderate def-H i galaxies have a similar range of detected $\Sigma_{\text{mol}}$ to HERACLES galaxies, with a slight extension to lower $\Sigma_{\text{mol}}$ , while large def-H i galaxies show a lower typical $\Sigma_{\text{mol}}$ .

Figure 8. Histograms of $\log\Sigma_{\text{mol}}$ within the $R_{\text{H}\,{\small\text{I}}}$ of galaxies, presented in the same way as Figure 5. Inside the $R_{\text{H}\,{\small\text{I}}}$ of moderate def-H i galaxies, environmental processes affect primarily the molecular gas with low surface density, while in large def-H i galaxies even the typical $\Sigma_{\text{mol}}$ is reduced.

The $\log\Sigma_{\text{mol}}$ histogram of the large def-H i galaxies shows a different behaviour. Unlike the moderate def-H i galaxies, there is no extended tail of lower $\Sigma_{\text{mol}}$ , but instead, it has a similar shape to the control sample, only shifted to a lower average $\Sigma_{\text{mol}}$ . The weighted median detected spaxel is $0.45$ dex lower (a factor of 2.8 or 64%) than the control sample and significant considering the uncertainties (Table 3), reflecting a larger change, which would be expected in larger def-H i galaxies. Thus, in the largest def-H i galaxies, environmental processes act to lower the typical $\Sigma_{\text{mol}}$ within $R_{\text{H}\,{\small\text{I}}}$ . These results support the interpretation that at moderate def-H i, the lower surface density molecular gas within $R_{\text{H}\,{\small\text{I}}}$ is impacted, while at large def-H i the typical $\Sigma_{\text{mol}}$ is lowered.

To summarise, we have presented clear evidence that environmental processes impact the molecular gas content of galaxies on 1.2-kpc physical scales, even within the H i truncation radius. The effect is primarily on the lower surface density molecular gas and can begin even at moderate def-H i, while in the largest def-H i systems, the denser molecular gas shows a reduction in the average molecular gas surface density.

5. Physical drivers of $\Sigma_{\text{mol}}$ reduction

The observed reduction of $\Sigma_{\text{mol}}$ at fixed $\Sigma_{\star}$ suggests that cluster environmental processes are impacting the molecular gas in VERTICO galaxies. This reduction could be due to direct stripping of the molecular gas, similar to the H i (e.g., Moretti et al. 2018; Cramer et al. 2020). Alternately, it could result from inefficient molecule formation caused by a change in the ISM physical conditions due to the stripping of the H i (e.g., Blitz & Rosolowsky 2006; Krumholz et al. 2009; Fumagalli et al. 2009). In this case, the suppressed molecular gas content is due to consumption by star formation and the lack of cold ISM replenishment (akin to strangulation/starvation). However, for this last scenario to work requires the assumption that the ISM in a galaxy is able to respond to any changes on a timescale equal to, or shorter than, the stripping timescale of the H i ( $\lesssim$ $200\,$ Myr; Vollmer et al. 2001; Roediger & Brüggen 2007; Tonnesen & Bryan 2009; Choi et al. 2022). Indeed, this might be expected due to the short lifetimes of GMCs ( $\lesssim$ 30 Myr, Chevance et al. 2020b). In this case, VERTICO galaxies should follow the scaling relations used to interpret the resolved molecular-to-atomic gas mass ratio ( $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ ) in normal galaxies, which is sensitive to the physical conditions within the ISM. Instead, we show here that these scaling relations break down when applied to VERTICO galaxies, and they cannot describe their $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ , suggesting that changes to the ISM physical conditions alone cannot explain the observed reduction of $\Sigma_{\text{mol}}$ .

More detailed modelling of the physical state of the cold gas reservoir, and its relationship to $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ , is outside the scope of this work and will be investigated in future VERTICO papers. For an analysis of $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ across the discs of VERTICO galaxies using radial profiles, we refer the reader to Villanueva 2022), who found elevated ratios in galaxies with perturbed and truncated H i morphologies. Here, we compare $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ in individual spaxels to two parameterisations of the ISM physical conditions to try and understand what physics might be responsible for the reduction in molecular gas content of the VERTICO galaxies. To be consistent with previous works, in this section, we correct our H i masses for a 36% contribution from Helium. Our analysis is restricted to detected spaxels only, and we also only consider spaxels with $\log\Sigma_{\star}\geq1.25$ so that we only study the ISM parameterisations where the control sample is dominated by detected spaxels. We also denote the total cold gas surface density as $\Sigma_\mathrm{gas}=\Sigma_{\text{H}\,{\small\text{I}}}+\Sigma_{\text{mol}}$ .

It is also useful to quantify how the removal of H i alone manifests in the scaling relations that we will investigate. To do this, we computed the offset between the median $\log\Sigma_{\text{H}\,{\small\text{I}}}$ of each VERTICO galaxy and the median of the control sample at fixed $\Sigma_{\star}$ . We then subtracted this offset (as a function of $\Sigma_{\star}$ ) from the spaxels of each HERACLES galaxy to create mock ‘H i-stripped’ galaxies. Using these mocks, we re-computed weighted histograms and statistics to estimate how the HERACLES control sample would shift if it had experienced H i reduction equal to what is observed for each VERTICO galaxy.

We use these mock galaxies in the following analysis to show where deviations from the ISM scaling relations can be explained by H i-only removal, and where we also require extra reduction of the molecular gas content. Indeed, it is our ability to perform this kind of spatially resolved analysis that allows us to identify these small differences, which are on the order of 0.1–0.2 dex, that would likely be washed out in studies considering only integrated quantities such as global H i and molecular gas deficiencies. For conciseness, we only show individual scaling relations for a subset of VERTICO galaxies in this section. Individual relations for all VERTICO galaxies, and the combined VERTICO and HERACLES relations are shown in Appendix 1.

First, we use a proxy for the hydrostatic pressure at the midplane of the disc, which was predicted to determine the molecular gas fraction by Elmegreen (1993, but see also Blitz & Rosolowsky 2006 and Leroy et al. 2008). We adopt the prescription from Blitz & Rosolowsky (2006), which we estimate as $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ . $\Sigma_{\star}$ and $\Sigma_\mathrm{gas}$ are the dominant terms in the Blitz & Rosolowsky (2006) prescription, with the stellar scale height and gas velocity dispersion having minor contributions. If this prescription is an accurate predictor of physics determining $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ , then we might expect that stripping of the ISM (predominantly H i) would reduce the pressure and cause $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ to decrease correspondingly.

In Figure 9, we show $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ as a function of $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ for two moderate def-H i galaxies (left two panels) and two large def-H i galaxies (right two panels), presented similarly to Figures 3 and 6. We also show, with red lines, how the HERACLES control sample median would shift if it had experienced H i reduction equal to what is observed for the displayed VERTICO galaxy. All four VERTICO galaxies have larger $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ than expected at fixed $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ showing that, despite a reduction in pressure due to the stripping of the H i, the $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ is larger. Thus, internal hydrostatic pressure cannot explain the $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ ratio in these VERTICO galaxies, suggesting that stripping of the gas (at least of the H i) is more dominant than the hydrostatic pressure, at least over the timescale traced by VERTICO galaxies. Further, there are regions within Figure 9 where the mock H i-stripped control sample over-estimates the $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ of VERTICO galaxies. This over-estimation is typically larger at lower pressures (corresponding to larger radii) and in the larger def-H i galaxies, which are the regions and galaxies where we observed molecular gas reduction in Section 4.2. Thus, changes in the hydrostatic pressure plus H i removal alone cannot explain the observed molecular gas reduction in these galaxies.

Figure 9. Ratio of molecular-to-atomic gas mass of individual spaxels as a function of the midplane pressure proxy $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ for four VERTICO galaxies, coloured and ordered by def-H i. Galaxies are presented in the same format as in Figures 3 and 6, with the left two belonging to the moderate def-H i sample and the right two to the large def-H i sample. The red lines show how the median and 16th and 84th percentiles of the HERACLES control sample would shift if it had experienced the same magnitude of H i reduction at fixed $\Sigma_{\star}$ as each VERTICO galaxy. At fixed midplane pressure, VERTICO galaxies have a larger $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ than expected for the same pressure in a control sample galaxy. However, this difference is not as large as the H i-only removal scenario, particularly at lower pressures and larger def-H i.

In Figure 10, we show that the story told by these four galaxies holds for our whole VERTICO sample, using normalised histograms of $\log(\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}})$ in two bins of $\log \Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ . The blue and red empty histograms correspond to the moderate and large def-H i samples, and the solid grey histogram shows the $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ histogram for the control sample. The solid blue and red histograms show the mock H i-stripped control samples, matched to the moderate and large def-H i samples respectively. In both bins, the VERTICO galaxies have elevated $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ compared to the HERACLES control sample; but not so elevated that they match with the histograms of the mock H i-stripped control samples.

Figure 10. Histograms of $\log(\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}})$ in bins of midplane pressure, presented in the same way as Figures 4 and 7. Additionally, solid blue and red histograms show how the $\log(\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}})$ histogram of the control sample would move if it had experienced H i removal of equal magnitude as the moderate and large def-H i samples, respectively. Moderate and large def-H i galaxies have systematically larger $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ than control sample galaxies at fixed $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ . Further, the H i-only removal scenario over-estimates the observed $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ , except in moderate def-H i galaxies in the largest $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ bin.

Quantitatively, in the lowest $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ bin, an H i-only removal scenario would result in $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ distributions with medians of $0.24\pm0.03$ dex and $0.65\pm0.03$ dex for the moderate and large def-H i populations, which over-estimate their observed values of $0.09\pm0.07$ dex and $0.31\pm0.04$ dex by 0.15 dex and 0.34 dex, respectively. In the highest $\Sigma_{\star}^{0.5}\Sigma_\mathrm{gas}$ bin, the large def-H i galaxies similarly show 0.2 dex smaller $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ elevation (median $0.61\pm0.04$ dex) than expected from H i removal alone ( $0.81\pm0.04$ dex). A K-S test between the def-H i samples and their respective mock control samples rules out that they are drawn from the same distribution ( $p<10^{-4}$ ). Conversely, the moderate def-H i sample in this bin matches its mock control sample, which is reflected in the expected and observed $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ medians of $0.47\pm0.03$ dex and $0.48\pm0.04$ dex, respectively. This match between the mock control sample and the moderate def-H i sample agrees with our results in Section 4.2, where no reduction in $\Sigma_{\text{mol}}$ was observed in these galaxies. Further, we can rule out that the elevated $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ is due to ram pressure providing additional pressure support and increasing the conversion rate of H i to molecular gas, as we found no evidence for elevated molecular gas content in any of our VERTICO galaxies (at least at our 1.2-kpc spatial scale). All these points indicate that some reduction in $\Sigma_{\text{mol}}$ is required to explain the $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ of VERTICO galaxies, but this parameterisation with hydrostatic pressure alone does not provide a physical explanation.

An alternate scenario based on ISM physics is that, due to the removal of H i, there is insufficient fuel to form new molecular gas. We test this hypothesis using the total gas surface density ( $\Sigma_\mathrm{gas}$ ), suggested by Krumholz et al. (2009, see also Fumagalli et al. 2009) to be the primary regulator of $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ with metallicity and the interstellar radiation field playing secondary roles. If the stripping of the ISM lowers the total surface density such that the formation of molecular gas is inefficient, we would expect VERTICO galaxies to lie along the same relation between $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ and $\Sigma_\mathrm{gas}$ as the control sample, as claimed in Fumagalli et al. (2009). In Figure 11, we show $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ as a function of $\Sigma_\mathrm{gas}$ in the same style as shown in Figure 9 for the pressure relation. The displayed moderate (left two panels) and large (right two panels) def-H i galaxies do not lie along the same relation as the control sample. Further, the removal of H i alone would over-estimate the observed $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ at fixed $\Sigma_\mathrm{gas}$ , particularly in lower $\Sigma_\mathrm{gas}$ regions.

Figure 11. Ratio of molecular-to-atomic gas mass of individual spaxels as a function of $\Sigma_\mathrm{gas}$ for the same galaxies as Figure 9. The $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ of VERTICO galaxies is systematically larger than in control sample galaxies at fixed $\Sigma_\mathrm{gas}$ , but not as large as the H i-only removal scenario shown with the red lines.

This trend remains for the whole sample, and in Figure 12, we show histograms of $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ in two bins of $\log\Sigma_\mathrm{gas}$ in the same style as Figure 10. In the lower $\Sigma_\mathrm{gas}$ bin, the expected median $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ from H i-only removal scenarios, $0.17\pm0.04$ dex and $0.63\pm0.04$ dex, over-estimate the observed values of $0.05\pm0.07$ dex $0.47\pm0.03$ dex for the moderate and large def-H i samples by 0.12 dex and 0.16 dex respectively. In the higher $\Sigma_\mathrm{gas}$ bin, the elevation in $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ of the large def-H i sample is smaller by 0.11 dex, with an expected median $0.74\pm0.03$ dex compared to $0.63\pm0.04$ dex observed, while the moderate def-H i galaxies agree more closely with their mock control sample, having $0.38\pm0.03$ dex expected and $0.41\pm0.05$ dex observed. We also can exclude that both def-H i samples are drawn from the same parent distribution as their respective mock control samples ( $p<10^{-4}$ ). Thus, in VERTICO galaxies it is not $\Sigma_\mathrm{gas}$ that sets $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ , and explaining the observed $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ requires a reduction in their molecular gas content that is not related to insufficient fuel or an inefficient conversion of H i into molecular gas.

Figure 12. Histograms of $\log(\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}})$ in bins of total cold gas surface density, presented in the same way as Figure 10. Moderate and large def-H i galaxies have systematically larger $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ than expected if molecular gas formation was set primarily by $\Sigma_\mathrm{gas}$ . The H i-only removal scenario over-estimates $\Sigma_{\text{mol}} / \Sigma_{\text{H}\,{\small\text{I}}}$ , except at large $\Sigma_\mathrm{gas}$ in moderate def-H i galaxies.

We must also consider if there are any systematic differences between how molecular gas is traced in VERTICO and HERACLES. As the data have been homogenised and both use the same constant $R_{21}$ and $\alpha_\mathrm{CO}$ to convert CO(2-1) luminosity to molecular gas mass, this scenario is unlikely. Although $R_{21}$ and $\alpha_\mathrm{CO}$ have their own uncertainties, this homogenisation means that our results are actually re-scaled differences between the observed CO(2-1) emission in our samples. Thus, to explain the difference between the typical $\Sigma_{\text{mol}}$ inside the $R_{\text{H}\,{\small\text{I}}}$ of large def-H i galaxies would require VERTICO galaxies to preferentially have CO(2-1)-dark gas such that the observed CO(2-1) luminosity under-estimates the actual amount of molecular gas present by a factor of 2.8 (0.45 dex).

Metallicity variations cannot explain a reduction in CO(2-1) emission, as while at smaller metallicities there is more CO(2-1) dark gas, cluster galaxies show similar or elevated metallicity over isolated systems (e.g., Skillman et al. 1996; Gupta et al. 2018; Franchetto et al. 2020). Changes to the UV radiation field due to reduction in star formation also are unlikely to cause any differences in CO(2-1) emission. VERTICO galaxies show reduced SFR surface densities up to $\sim$ 1 dex compared to star-forming galaxies at fixed $\Sigma_{\star}$ (Brown et al., in preparation); however, in the same regions, we observe reduced H i and molecular gas densities with similar magnitude, which reduces the gas column’s ability to self-shield. Thus, any changes in the UV field are likely balanced by changes in the gas density (Krumholz et al. 2009). VERTICO galaxies also show similar $R_{21}$ values to normal galaxies from xCOLD GASS (Brown et al. 2021), suggesting that there is no significant difference in the CO(2-1) excitation that could under-estimate the amount of molecular gas. Additionally, Smith et al. (2014) and Glover & Smith (2016) found that CO(2-1)-dark gas is typically diffuse and low surface density, making it prone to turbulent compression and gravitational collapse and thus less likely to exist in the central regions of galaxies that we have studied.

Further evidence that VERTICO galaxies do not have substantial CO(2-1)-dark gas comes from studies of the Kennicutt–Schmidt (Schmidt 1959; Kennicutt 1998) star formation law. If the CO(2-1) emission line underestimated the total molecular gas mass, then galaxies would deviate from the Kennicutt–Schmidt law and appear to have elevated SFR at fixed molecular gas mass. However, VERTICO galaxies show similar galaxy-to-galaxy variation in their resolved Kennicutt–Schmidt relations to what is observed in star-forming galaxies (Vollmer et al. 2012; Jiménez-Donaire et al. 2022). This consistency with unperturbed systems implies that the physics of star formation, and by extension the ISM physical conditions giving rise to the CO(2-1) emission, do not show sufficient difference between VERTICO galaxies and the control sample to explain our results. Thus, we need to invoke a process that is unrelated to the ISM physical conditions and can reduce molecular gas content on timescales similar to the H i gas. Direct stripping of the molecular gas meets this requirement.

There is evidence that not all molecular gas in galaxies resides in bound GMCs, with 20–50% (Polk et al. 1988;Wilson & Walker 1994; Pety et al. 2013; Roman-Duval et al. 2016; Chevance et al. 2020a; Saintonge & Catinella 2022) existing in a more diffuse phase or a molecular thick disc that has velocity dispersion similar to the H i e.g. (e.g. Caldú-Primo et al. 2013, 2015; Mogotsi et al. 2016), making it much easier to strip than GMC gas. Our data cannot differentiate between diffuse and dense gas; however, the 64% reduction (Section 4.2) in the typical molecular gas content within $R_{\text{H}\,{\small\text{I}}}$ observed in our large def-H i galaxies is greater than the ‘typical’ diffuse gas fractions mentioned above. This is interesting, as there is also evidence that GMCs can be disrupted by the ICM. Cramer et al. (2020) found that only 12% of the molecular gas on the leading edge of NGC 4402 exists in GMCs that have the survived the ongoing stripping, while the remaining gas exists in a trailing plume with lower surface density. The authors suggest that this large fraction of lower surface density molecular gas implies that an indirect stripping process, by which GMCs are converted to lower surface density molecular gas or H i before being directly stripped, could be taking place.

This interpretation is also supported by simulations, which have shown that the multi-phase structure of the ISM allows stripping to affect gas within the truncation radius (e.g., Quilis, Moore, & Bower 2000; Tonnesen & Bryan 2009). Further, momentum-transfer from the ICM to the ISM through turbulent mixing can effectively ablate the dense, cold ISM into warmer phases (e.g., Gronke & Oh 2018; Sparre, Pfrommer, & Vogelsberger 2019) and accelerate it out of the galaxy (Tonnesen & Bryan 2021; Choi et al. 2022). The typical ICM pressures expected for Virgo galaxies are within the ranges used in these simulations (e.g., Vollmer et al. 2001; Roediger & Brüggen 2007; Chung et al. 2007; Abramson et al. 2011), making stripping of the molecular gas the most feasible explanation for our observed molecular gas reduction. However, future work with higher spatial resolution observations is needed to allow quantification of how the diffuse and dense phases of the ISM are impacted by the Virgo environment, and search for signatures of direct stripping of the molecular gas.

6. Summary

In this paper, we used spatially resolved scaling relations between $\Sigma_{\star}$ , $\Sigma_{\text{H}\,{\small\text{I}}}$ , and $\Sigma_{\text{mol}}$ to study how cluster environmental processes impact the cold gas content within the optical extent of massive ( $\log(M_{\star}/\mathrm{M}_{\odot})\geq10$ ) galaxies on 1.2-kpc scales. We presented a new view of how the H i gas distribution changes with increasing H i deficiency, and we have shown, statistically, how environmental processes impact the spatial distribution of the molecular gas, the distribution of $\Sigma_{\text{mol}}$ , and molecular-to-atomic gas mass ratio using observations from the VERTICO survey (Brown et al. 2021). Our main results are as follows:

• • While the stripping of gas from galaxies impacts the H i gas at large radii, it has the additional effect of lowering the average H i surface density of the remaining gas within the H i-stripping radius, $R_{\text{H}\,{\small\text{I}}}$ . This reduction in surface density increases as the gas removal proceeds inward until the densest H i regions can no longer form (Figures 3, 4, 5).

• • Some galaxies show co-spatial truncation of their H i and molecular gas (Figure 6), indicating that environmental processes can effectively remove both phases on similar timescales in the same spatial region.

• • Overall, the molecular gas is less affected than the H i, and we find evidence that environmental mechanisms first reduce the lower surface density molecular gas that exists preferentially at lower $\Sigma_{\star}$ (and thus typically larger radii) in galaxies (Figure 7). However, we find that once the removal of H i proceeds sufficiently far within the optical disc, the average surface density of the remaining molecular gas within $R_{\text{H}\,{\small\text{I}}}$ is reduced (Figure 8).

• • At the spatial scales we are sensitive to, we see no evidence for environment-driven elevation of the H i or molecular gas surface density outside the ranges observed in star-forming galaxies (Figures 4, 7).

• • VERTICO galaxies have systematically larger molecular-to-atomic gas mass ratios in individual regions compared to the control sample, but smaller values than expected from the removal of H i alone (Figures 9, 11). Using proxies for the physical conditions in the ISM, we showed that changes to the physical state of the gas could not explain the observed molecular-to-atomic gas mass ratios (Figures 10, 12). Thus, ablation and/or stripping of the molecular gas on a similar timescale as the H i in a given 1.2-kpc region is the best explanation for the reduction of molecular gas content in VERTICO galaxies.

Our results agree with, and extend upon, previous literature focused on gas stripping and environmental effects in cluster galaxies (Cortese et al. 2021; Boselli et al. 2022). Combined with the recent, complementary studies from the VERTICO team (Brown et al. 2021; Villanueva 2022; Zabel et al. 2022), a coherent evolutionary scenario for Virgo cluster galaxies begins to emerge. Compared to isolated systems, statistically, Virgo galaxies have normal molecular gas content (Boselli et al. 2014; Cortese et al. 2016) and elevated molecular-to-atomic gas mass ratios (Mok et al. 2016; Villanueva 2022) when considering their integrated properties only. Subsequently, there is poor correlation between H i and molecular gas deficiency measurements (Zabel et al. 2022), but this does not mean that the molecular gas is not affected by environmental processes. Instead, major changes in the integrated molecular gas mass occur only once a significant fraction of the cold gas within the stellar disc has been affected.

Our analysis reveals how the molecular gas content of galaxies is affected by cluster environmental processes in a spatially resolved way. We have shown that the reduction of molecular gas begins at the edge of the stellar disc even in moderately H i-deficient galaxies, while in the highest H i-deficiency systems the molecular gas becomes significantly reduced across the disc. While this result has been shown using integrated measurements and radial profiles (Villanueva 2022; Zabel et al. 2022), our approach allows us to quantify this effect as a function of $\Sigma_{\star}$ with exquisite detail, and has allowed us to disentangle between environmentally driven and secular/internal processes.

Interestingly, previous works have found evidence for the compression of molecular gas due to the ICM–ISM interaction (e.g., Lee et al. 2017; Moretti et al. 2020a,b), particularly in the early stages of infall. We find no evidence for elevated H i or molecular gas surface densities, suggesting that in our sample, which spans a larger (and later) range of infall stages and at the 1.2-kpc scales we are sensitive to, environment-driven elevation of cold gas content is not significant.

Our results also emphasise that the gas removal process is more complex than a simple outside-in stripping scenario, and it is clear that the environmental processes affect both the H i and molecular gas within the H i truncation radius (see also Lee et al. 2022), which is typically assumed to be the boundary of where stripping is efficient. Even at moderate H i deficiency, we find that the average surface density of H i remaining within the disc is reduced, indicating that the environment does not just affect the H i gas in the outskirts. This result can be understood as the multi-phase ISM enabling the ICM to stream through low-density regions and affect the gas across the disc (e.g., Tonnesen & Bryan 2009; Lee et al. 2020). However, we have also shown that the molecular gas is affected similarly, although it requires a stronger ICM–ISM interaction (as traced by more H i lost). Thus, the ICM can ablate and strip gas across the entire discs of galaxies and does not discriminate which phase is affected so long as the interaction is strong enough (e.g., Tonnesen & Bryan 2021; Choi et al. 2022). Indeed, this agrees well with observations of larger molecular gas masses in the cometary tails of stripped galaxies when the associated $\mathrm{H}\,{\alpha}$ emission is brighter or the stripping is more advanced (Moretti et al. 2018). In these cases, it is likely that the stronger stripping has led to a larger mass-loading of the outflowing gas. Future work with larger samples of field and cluster galaxies, but in particular, higher resolution observations that can separate dense and diffuse CO(2-1) emission will be key to understanding the details of gas suppression in cluster galaxies.

Last, these results have interesting implications for the quenching of star formation in cluster galaxies. In particular, the reduction of molecular gas within the inner regions of galaxies will cause star formation to shut down earlier, and with a different distribution, than just radial truncation. Comparison to deep $\mathrm{H}\,{\alpha}$ observations of VERTICO galaxies from the VESTIGE (A Virgo Environmental Survey Tracing Ionised Gas Emission, Boselli et al. 2018a) survey would enable the study of changes to recent, massive star formation in regions where molecular gas is reduced. Further, future observations with sensitive optical integral-field spectroscopic instruments such as MUSE (Multi Unit Spectroscopic Explorer, Bacon et al. 2017) will allow the study of the physics of galaxy quenching in unprecedented detail, including the additional impact of feedback from active galactic nuclei (e.g., Boselli et al. 2016b; George et al. 2019; Peluso et al. 2022). In particular, employing spectro-photometric fitting techniques (e.g., Fossati et al. 2018) will enable high spatial resolution studies of star formation and star-formation history reconstruction across the discs of Virgo cluster galaxies. These results can then be connected to cold gas reduction on sub-kpc scales where ISM physics operates.