© I. Pessa et al. 2022

1. Introduction

The conversion of cold dense molecular clouds into stars ultimately occurs when the supporting pressures are insufficient to prevent gravitational collapse. While simple calculations only include gas pressure, both magnetic fields (e.g., Shu et al. 1987) and turbulence (Mac Low & Klessen 2004) have been proposed as mechanisms able to prevent gravitational collapse. In the case of magnetic fields, the collapse is prevented because the neutral hydrogen is coupled with ionized hydrogen, which is tied to the interstellar magnetic field and therefore resists collapse. For the latter proposed mechanism, supersonic turbulent motions act as an additional source of pressure. Once this pressure is no longer sufficient to support the self-gravity of the cloud, it collapses. The collapse of the molecular cloud, which might be triggered by an external source of pressure, leads to its fragmentation (Hoyle 1953), where individual fragments will form stars.

Despite the complexity of this process, several studies have reported correlations between the locally (approximately kpc and sub-kpc scales) measured star formation rate surface density (Σ_SFR), and spatially resolved quantities, such as the local stellar mass surface density (Σ_⋆), also known as the “resolved” star formation main sequence (rSFMS; Cano-Díaz et al. 2016; Abdurro’uf & Akiyama 2017; Hsieh et al. 2017; Lin et al. 2019; Morselli et al. 2020; Ellison et al. 2021), or the local molecular gas surface density (Σ_mol), also known as the resolved Kennicutt–Schmidt relation (rKS; Bigiel et al. 2008, 2011; Leroy et al. 2008, 2013; Blanc et al. 2009; Onodera et al. 2010; Schruba et al. 2011; Ford et al. 2013; Kreckel et al. 2018; Williams et al. 2018; Dey et al. 2019). The scatter of these correlations is expected to dramatically increase below a critical spatial scale due to statistical undersampling of the time evolution of the star formation process (Schruba et al. 2010; Feldmann et al. 2011; Kruijssen & Longmore 2014; Kruijssen et al. 2018). Pessa et al. (2021) recently confirmed such an increase in scatter with decreasing physical resolution down to ∼100 pc in a sample of 18 nearby galaxies. These local relationships also have well-studied global manifestations, the star-forming main sequence of galaxies (SFMS; Brinchmann et al. 2004; Daddi et al. 2007; Noeske et al. 2007; Salim et al. 2007; Lin et al. 2012; Whitaker et al. 2012; Speagle et al. 2014; Saintonge et al. 2016; Popesso et al. 2019), and the galaxy-integrated Kennicutt–Schmidt relation (Schmidt 1959; Kennicutt 1998; Wyder et al. 2009; Genzel et al. 2010; Tacconi et al. 2010) express relationships among the same quantities for entire galaxies. These galactic-scale relations are key to our current understanding of galaxy evolution and star formation across cosmic time. Studying their spatially resolved versions provides critical insights into their physical origin.

These resolved correlations arise as a result of some mechanism controlling the formation of stars. Many studies have investigated the dependence on Σ_⋆ (through the hydrostatic pressure exerted by the galactic potential) and by Σ_mol (being the fuel to form stars), as well as a combination of these parameters (e.g., Matteucci et al. 1989; Shi et al. 2011, 2018; Dey et al. 2019; Barrera-Ballesteros et al. 2021a). When approaching this topic empirically, some authors have studied a 2D plane in the 3D space spanned by Σ_SFR, Σ_⋆ and Σ_mol (Dib et al. 2017; Lin et al. 2019; Sánchez et al. 2021), which would imply a scenario where both Σ_⋆ and Σ_mol are responsible for modulating the formation of stars. While there is consensus that these are the main quantities that correlate with SFR to first order, there is debate about the universality of these correlations. Whereas Sánchez et al. (2021), using data from the CALIFA survey (Sánchez et al. 2012), concluded that the scatter in these relations is fully dominated by individual errors, Ellison et al. (2021) used data from ALMaQUEST (Lin et al. 2020) and found that the scatter in these correlations is dominated by galaxy-to-galaxy variations. Similarly, in Pessa et al. (2021), using high physical resolution data provided by MUSE, the authors found not only galaxy-to-galaxy variations, but also significant variations in these scaling relations across galactic environments. That is, they found different relations for spiral arms, disk, bars, centers and rings.

In this paper, we aim to test the universality of a two-dimensional (2D) planar star formation relation, defined as ${\mathrm{\Sigma }}_{\mathrm{SFR}}\propto {\mathrm{\Sigma }}_{\star }^a{\mathrm{\Sigma }}_{\mathrm{mol}}^b$. If the SFR in a given region is primarily driven by these two quantities, then, this ‘star-forming plane’ should not change significantly between individual galaxies or among different galactic environments. On the other hand, significant variations in the best-fit star-forming plane among different locations would indicate that such a relationship offers an incomplete description of the data and hint at additional quantities or functional forms that are key to regulating star formation in galaxies.

The Physics at High Angular resolution in Nearby GalaxieS (PHANGS¹) surveys (Leroy et al. 2021a; Emsellem et al. 2022; Lee et al. 2022) provide the opportunity to explore the relation between these three quantities at a high physical resolution (∼150 pc), making it possible to isolate galactic environments and study them separately. PHANGS targets probe a diversity of environments (centers and bulges, disks, bars, spiral arms) spanning a range of physical conditions (gas surface density, stellar surface density, dynamical pressure, orbital period, shear rate, streaming motions, gas phase metallicity; Meidt et al. 2018, 2020; Jeffreson et al. 2020; Kreckel et al. 2020; Emsellem et al. 2022; Querejeta et al. 2021). Recent works have shown that ∼100 − 150 pc scale surface density, velocity dispersion, and dynamical state of the molecular gas in these environments are sensitive to the local conditions (Colombo et al. 2014; Sun et al. 2020a; Rosolowsky et al. 2021) in a manner that appears to influence how molecular clouds form stars (Meidt 2016; Querejeta et al. 2021). Thus, investigating differences in the interplay of Σ_SFR, Σ_⋆ and Σ_mol across these galactic environments could be key to understand the galaxy-to-galaxy variations reported by Ellison et al. (2021) and Pessa et al. (2021), as relative contributions from different environments will vary from one galaxy to another.

The paper is structured as follows: in Sect. 2 we describe the data and data products used in this work. In Sect. 3 we describe in detail the methodology adopted to perform our analyses. In Sects. 4 and 5 we present our findings and discussions. Finally, our main conclusions are presented in Sect. 6.

2. Data

In this section, we cover only the main aspects of our data set. We refer the reader to Leroy et al. (2021b), Emsellem et al. (2022) and Pessa et al. (2021) for a more detailed description of our data. We use a sample of 18 star-forming galaxies from the overlap of the PHANGS–ALMA and PHANGS–MUSE samples. In order to resolve the typical separation between star-forming regions (∼100 pc), and limit the effect of extinction, the galaxies studied in this paper have been selected to have distances less than 20 Mpc and low inclinations (i < 60°). Our sample is summarized in Table 1 where we use the global parameters from Leroy et al. (2021a), which leveraged the distance compilation of Anand et al. (2021) and the galaxy orientations determined by Lang et al. (2020).

2.1. VLT MUSE

We make use of the PHANGS–MUSE survey (P.I.: E. Schinnerer; Emsellem et al. 2022). This survey employs the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2014) optical integral field unit (IFU) mounted on the Very Large Telescope (VLT) Unit Telescope 4 to mosaic the star-forming disk of 19 galaxies from the PHANGS sample. This sample corresponds to a subset of the 90 galaxies from the PHANGS–ALMA survey (P.I.: E. Schinnerer; Leroy et al. 2021a). For the sake of homogeneity of the data set, one galaxy from the PHANGS–MUSE sample (NGC 0628) has been excluded because its MUSE mosaic was obtained using a different observing strategy.

The mosaics consist of 3 to 15 individual MUSE pointings, each with a total on-source exposure time of 43 min. Nine out of the 18 galaxies were observed with adaptive optics (AO) assistance. These galaxies are marked with a black dot in the first column of Table 1. Each pointing provides a 1′×1′ field of view sampled at $0\stackrel{″}{.}2$ per pixel, with a typical spectral resolution of ∼2.5 Å (∼70 km s⁻¹) covering the wavelength range of 4800 − 9300 Å. Observations were reduced using recipes from the MUSE data reduction pipeline provided by the MUSE consortium (Weilbacher et al. 2020), executed with ESOREX using the python wrapper developed by the PHANGS team² (Emsellem et al. 2022). Once the data have been reduced, we have used the PHANGS data analysis pipeline (DAP) to derive various physical quantities, as described in detail in Emsellem et al. (2022). DAP is based on the GIST pipeline (Bittner et al. 2019), and consists of a series of modules that perform single stellar population (SSP) fitting and emission line measurements to the full MUSE mosaic. Some outputs from the pipeline relevant for this study are described in Sects. 2.4 and 2.5.

2.2. ALMA CO mapping

Our sample of 18 galaxies have CO(J = 2 − 1) [hereafter CO(2−1)] data from the PHANGS–ALMA survey (P.I.: E. Schinnerer; Leroy et al. 2021a). We used the ALMA 12 m and 7 m arrays combined with the total power antennas to map CO(2−1) emission at a spatial resolution of $1\stackrel{″}{.}0-1\stackrel{″}{.}5$ (version 4.0 of internal distribution, which is the first public data release). The molecular gas surface density maps have a typical uncertainty of ∼1.2 M_⊙ pc⁻² at a spatial resolution of 150 pc. We use integrated intensity maps and associated statistical uncertainty maps constructed using the PHANGS–ALMA “broad masking” scheme. These broad masks are designed to incorporate all emission detected at any scale in the PHANGS–ALMA data. As a result they have high completeness, that is to say, they include most of the flux in the galaxy, at the expense of also having more low-confidence pixels than a more stringent masking technique (see Leroy et al. 2021a,b, for details and completeness estimates at 150 pc). A signal-to-noise (S/N) cut of 1 is then applied to drop the most uncertain emission. The strategy for observation, data reduction and product generation are described in Leroy et al. (2021a,b).

We adopt the local gas-phase metallicity (in solar units) prescription for our fiducial α_CO conversion factor following Sun et al. (2020b) and partially based on Accurso et al. (2017), following α_CO = 4.35(Z/Z_⊙)^−1.6 M_⊙ pc⁻² (K km s⁻¹)⁻¹, adopting a ratio CO(2−1)-to-CO(1−0) = 0.65 (Leroy et al. 2013, 2021b; den Brok et al. 2021; Saito et al., in prep.). Metallicity gradients are measured from the gas-phase abundances in H II regions, as explained in Kreckel et al. (2020) and Santoro et al. (2022). Azimuthal variations in the gas-phase metallicity have been previously reported (Ho et al. 2017; Kreckel et al. 2020; Williams et al. 2022), however, these variations are small (0.04 − 0.05 dex), implying variations of ∼0.06 dex in α_CO and therefore do not impact our results. We test the robustness of our results against a constant α_{CO(1 − 0)} = 4.35 M_⊙ pc⁻² (K km s⁻¹)⁻¹, as the canonical value for our Galaxy (Bolatto et al. 2013) in Sect. 5.3.

2.3. Environmental masks

We used the environmental masks described in Querejeta et al. (2021) to morphologically classify the different environments of each galaxy and label them as disks, spiral arms, rings, bars and centers. This classification was done using photometric data mostly from the Spitzer Survey of Stellar Structure in Galaxies (S⁴G; Sheth et al. 2010). In brief, disks and centers are identified via 2D photometric decompositions of 3.6 μm images (see, e.g., Salo et al. 2015). A central excess of light is labeled as center, independently of its surface brightness profile. The size and orientation of bars and rings are defined visually on the NIR images; for S⁴G galaxies, the classification follows Herrera-Endoqui et al. (2015). Finally, spiral arms are only defined when they are clearly dominant features across the galaxy disk (excluding flocculent spirals). A log-spiral function is fitted to bright regions along arms on the NIR images, and assigned a width determined empirically based on CO emission.

2.4. Stellar mass surface density maps

The PHANGS–MUSE DAP (Emsellem et al. 2022) includes a stellar population fitting module, where a linear combination of single stellar population templates of specific ages and metallicities are used to reproduce the observed spectrum. We assume a Calzetti et al. (2000) extinction law to account for extragalactic attenuation in the fitting. This permits us to infer stellar population properties from an integrated spectrum, such as mass- or light-weighted ages, metallicities and total stellar masses, together with the underlying star formation history. Before doing the SSP fitting, we correct the full mosaic for Milky Way extinction assuming a Cardelli et al. (1989) extinction law and the E(B − V) values obtained from the NASA/IPAC Infrared Science Archive³ (Schlafly & Finkbeiner 2011).

In detail, our spectral fitting pipeline performs the following steps: First, we use a Voronoi tessellation (Cappellari & Copin 2003) to bin our MUSE data to a minimum S/N of ∼35, computed at the wavelength range of 5300 − 5500 Å. We use then the Penalized Pixel-Fitting (pPXF) code (Cappellari & Emsellem 2004; Cappellari 2017) to fit the spectrum of each Voronoi bin. To fit our data, we used a grid of templates consisting of 13 ages, logarithmically spaced and six metallicity bins. We fit the wavelength range 4850 − 7000 Å, in order to avoid spectral regions strongly affected by sky residuals. We used templates from the eMILES (Vazdekis et al. 2010, 2012) database, assuming a Chabrier (2003) IMF and BaSTI isochrone (Pietrinferni et al. 2004) with a Galactic abundance pattern.

The stellar mass map is then reconstructed using the weights assigned by pPXF to each SSP spectrum, given that the current stellar mass of each template is known. Finally, we have identified foreground stars as velocity outliers in the SSP fitting and we have masked those pixels for the analysis carried out in this paper.

2.5. Star formation rate measurements

As part of the PHANGS–MUSE DAP (Emsellem et al. 2022), we fit single Gaussian profiles to a number of emission lines for each pixel of the final combined MUSE mosaic of each galaxy in our sample. By integrating the flux of the fitted profile in each pixel, we construct emission line flux maps for every galaxy. We calculate SFR from extinction-corrected Hα. In detail, we dereddened the Hα fluxes, assuming that Hα_corr/Hβ_corr = 2.86, as appropriate for a case B recombination (Osterbrock 1989), temperature T = 10⁴ K and density n_e = 100 cm⁻³, following:

$$\begin{array}{c}\hfill \mathrm{H}{\alpha }_{\mathrm{corr}}=\mathrm{H}{\alpha }_{\mathrm{obs}}(\frac{{(\mathrm{H}\alpha /\mathrm{H}\beta )}_{\mathrm{obs}}}{2.86}{)}^{\frac{k_{\alpha }}{k_{\beta }-k_{\alpha }}},\end{array}$$(1)

where Hα_corr and Hα_obs correspond to the extinction-corrected and observed Hα fluxes, respectively, and k_α and k_β are the values of a given extinction curve at the wavelengths of Hα and Hβ. Opting for an O’Donnell (1994) extinction law, we use k_α = 2.52, k_β = 3.66 and R_V = 3.1.

Next, we remove pixels that are dominated by active galactic nuclei (AGN) or low-ionization nuclear emission-line regions (LINER) ionization from our sample, performing a cut in the Baldwin–Phillips–Terlevich (BPT; Baldwin et al. 1981) diagram using the [O III]/Hβ and [S II]/Hα line ratios, as described in Kewley et al. (2006).

For the remaining pixels, we determined the fraction of the Hα emission tracing local star formation (C_H II) and the fraction deemed to correspond to the diffuse ionized gas (DIG). DIG is a warm (10⁴ K), low density (10⁻¹ cm⁻³) phase of the interstellar medium (Haffner et al. 2009; Belfiore et al. 2015) produced primarily by photoionization of gas across the galactic disk by photons that escaped from H II regions (Flores-Fajardo et al. 2011; Zhang et al. 2017; Belfiore et al. 2022).

To this end, we use the [S II]/Hα ratio to estimate C_H II in each pixel, following Blanc et al. (2009) and Kaplan et al. (2016). After performing this correction and removing the DIG contribution of the Hα flux, we rescale the DIG-corrected Hα map by 1 plus the fraction of flux removed, to keep the total Hα flux constant. We perform this rescaling because photons that ionize the DIG are believe to originally have leaked out from H II regions⁴. This correction represents then a spatial redistribution of the Hα flux. We refer the reader to Pessa et al. (2021) for a detailed description of this procedure. This approach permits us to estimate a star formation rate in pixels contaminated by non-star-forming emission. A S/N cut of 3 for Hα and 1 for Hβ was then applied before computing the star formation rate surface density map using Eq. (2), effectively removing ∼13% of pixels in our sample. While 13% is not insignificant, we note that the majority have C_H II ≈ 0, and therefore our S/N cut does not largely impact our results. Pixels below this S/N cut, pixels with C_H II ≤ 0 or pixels where Hα_obs/Hβ_obs < 2.86, are considered nondetections (see Sect. 3.2). To calculate the corresponding star formation rate from the Hα flux, after correcting it for internal extinction and DIG contamination, we adopted the prescription of Calzetti (2013):

$$\begin{array}{c}\hfill \frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}=5.5\times {10}^{-42}\frac{\mathrm{H}{\alpha }_{\mathrm{corr}}}{{\mathrm{erg}\mathrm{s}}^{-1}}·\end{array}$$(2)

This equation is scaled to a Kroupa universal IMF (Kroupa 2001). Differences with the Chabrier IMF assumed for the SSP fitting are expected to be small (∼1.05; Kennicutt & Evans 2012). With these steps we obtain SFR surface density maps for each galaxy in our sample. We acknowledge that Eq. (2) assumes a fully sampled IMF, and that the lowest SFR pixels (especially at the ∼150 pc resolution of our data) may not form enough stars to fully sample the IMF. Hence, the measured SFR is more uncertain in this regime. We decided to exclude from our analysis the regions with low coverage of SFR or molecular gas surface density (see Sect. 3.2), hence minimizing the impact of the low SFR regime on our results.

2.6. Probing larger spatial scales

In Sect. 4.2, we investigate the effect of spatial resolution on our measurements. To do this, we first convolve our MUSE maps to a common fixed reolution of 150 pc, and then we rebin our 150 pc resolution MUSE and ALMA maps to have pixel sizes of 150 pc, 200 pc, 300 pc, 500 pc and 1 kpc. Then we replicate our measurements using these rebinned maps. We conduct an identical rebinning process for all three relevant quantities: stellar mass surface density, star formation rate surface density and molecular gas mass surface density (Pessa et al. 2021, for a detailed description). After beginning with matched resolution 150 pc MUSE and ALMA data, we favor this rebinning approach rather than, for example, convolution to a coarser Gaussian beam, because it minimized pixel-to-pixel covariance and yields approximately statistically independent measurements. Our core results use the 150 pc pixels, which are larger than the native spatial resolution of the maps for all galaxies, except for NGC 1672, which has an ALMA native resolution of ∼180 pc.

The rebinning step is followed by an inclination correction, simply using a cos(i) multiplicative term, where i is the inclination of each galaxy as listed in Table 1 (adopted from Lang et al. 2020). All following results and conclusions in the next Sections pertain to a fixed spatial resolution of 150 pc, unless specifically stated (see Sect. 4.2).

3. Methods

In this section, we present our methodology to fit a 2D plane (i.e., power law) predicting log₁₀Σ_SFR as function of log₁₀Σ_⋆ and log₁₀Σ_mol. We fit planes separately to the data for each distinct galactic environment, as well as to the full sample. Because the fitting can be potentially biased by the influence of nondetections in either Σ_SFR or Σ_mol, we explain our approach to nondetections first and then describe the full fitting method.

3.1. Nondetections in our data

In Pessa et al. (2021), we found that the adopted treatment of nondetections can have a considerable impact on the derived slope of fitted scaling relations. Here, nondetections (N/Ds) are pixels with a value of Σ_SFR or Σ_mol lower than our defined detection threshold. As stellar mass is detected essentially everywhere in our maps at high significance, nondetections in Σ_⋆ are essentially a nonissue.

This impact of N/Ds becomes particularly strong when performing measurements at high spatial resolution, where a larger fraction of pixels are deemed to be N/Ds. For more discussion of how such sparse maps emerge from timescale effects or other stochasticity, see Pessa et al. (2021) and references therein.

In Pessa et al. (2021), we overcame this issue by binning the data before fitting a power law. Here we are trying to fit a 2D plane, and discretizing the data in higher dimensions becomes problematic, so we do not proceed in the same way. Instead, we opt here to focus our analysis only on those Σ_⋆ ranges that have high detection fractions in both Σ_SFR and Σ_mol. Here the detection fraction is defined as the fraction of pixels with a measurement of Σ_SFR (or Σ_mol) above our detection threshold in a given Σ_⋆ bin. As a reminder, we required S/N > 1 for the ALMA CO intensity, S/N > 3 for the Hα and S/N > 1 for the Hβ at the native resolution, and N/Ds are a nonissue for stellar mass.

The top panels of Fig. 1 show the detection fraction of Σ_SFR and Σ_mol in each bin of Σ_⋆, and the bottom panels show rSFMS and rMGMS binned for each individual environment (colors) and for the full sample (black). Two types of lines are shown, solid lines represent the binned trend accounting for N/Ds, and dashed lines represent the binned trend neglecting N/Ds. It is easy to see in the bottom panels that for all environments the two lines deviate, that is, neglecting N/Ds from the analysis systematically leads to flatter slopes. This bias in the slope is due to varying fractions of N/Ds across the Σ_⋆ range, usually being higher at the low Σ_⋆ end. To minimize the impact of N/Ds on our measurements, we confine our analysis only to those Σ_⋆ ranges, where the detection fraction (of both Σ_SFR and Σ_mol) is higher than 60%. In Appendix B, we discuss how using a different detection fraction threshold impacts our results, and we show that a threshold of 60% provides a good compromise between minimizing the impact of N/Ds in our analysis, while maintaining the statistical significance of our data.

Fig. 1. rSFMS (bottom left) and rMGMS (bottom right) measured at 150 pc resolution for all galaxies in our sample. For reference and to guide the eye, the various solid colored lines represent the binned trends (bins of 0.15 dex) obtained for each environment, accounting for nondetections as detailed in Pessa et al. (2021). The black solid line shows the same measurement for all environments simultaneously. The dashed lines show the binned trends obtained when nondetections are neglected from the analysis. The panels on the top row show the detection fraction of ΣSFR (left) and Σmol (right), defined as the fraction of pixels with a measurement above our detection threshold in each bin of Σ⋆. The horizontal dotted line marks the detection fraction level of 60%.

We impose the detection fraction threshold independently for each environment, that is, the range of Σ_⋆ values used for ‘spiral arms’ is different to that used for ‘disks’. This approach allows us to maximize the number of data points used in the analysis.

We do caution that our approach to the different galactic environments might impose some bias on the results. The spiral arms defined by Querejeta et al. (2021) tend to have more gas and more star formation at fixed Σ_⋆. Because they were identified partially based on multiwavelength data that trace gas and SFR, this is somewhat by construction. It is not a priori certain that spiral arms must behave distinctly from other environments in the Σ_⋆ − Σ_SFR − Σ_mol space, but some bias might be expected from how we set up the analysis. This also manifests in the detection threshold cuts: because spiral arms have higher Σ_mol and Σ_SFR at fixed Σ_⋆, detections extend to lower Σ_⋆ for these environments and we might also expect a different best-fitting plane. This is already somewhat evident in Fig. 1, where the “Sp. arm” environment shows higher detection fractions at fixed Σ_⋆ than other cases.

It is worth noting that we are not dropping a large fraction of our data by imposing these thresholds. Indeed, only ∼20% of our 150 pc pixels are removed from the entire sample across all environments. Given that different galactic environments sample different parts of the galactic disk and therefore different (typical) surface densities, the adopted N/D thresholds are most relevant for the disk environment, which extends to the largest galactic radii and therefore (also) samples low surface densities. As a result, ∼38% of its original pixels do not satisfy the threshold. For all other environments, the fraction of pixels dropped is only < 5%. After performing this cut, the disk still remains the galactic environment with the largest number of pixels. Thus, we are only mildly increasing the statistical uncertainty of our results with this approach, and it permits us to reduce any bias in the reported quantities by considering only those Σ_⋆ regimes not dominated by nondetections.

3.2. Fitting technique

In this work, we aim at finding the best set of parameters C_⋆, C_mol and C_norm that describe the star-forming plane conformed by the quantities logΣ_SFR, logΣ_⋆ and logΣ_mol, such that:

$$\begin{array}{c}\hfill log{\mathrm{\Sigma }}_{\mathrm{SFR},i}=C_{\star ,i}log{\mathrm{\Sigma }}_{\star }+C_{\mathrm{mol},i}log{\mathrm{\Sigma }}_{\mathrm{mol}}+C_{\mathrm{norm},i},\end{array}$$(3)

where the subindex i stands for each galactic environment. If Σ_⋆ and Σ_mol are the only quantities that determine Σ_SFR in a given region of a galaxy, then this plane should not change between different galactic environments. However, given that galactic environments represent a diversity of other physical conditions including gas phase metallicity, stellar age, stellar geometry (flattened vs. triaxial), shear rate, radial flows and the gas structure and organization, then we might expect there to be additional factors that influence gas stability and star formation at fixed Σ_mol and Σ_⋆ (i.e., Hunter et al. 1998; Martig et al. 2009; Krumholz et al. 2018; Meidt et al. 2018, 2020; Gensior et al. 2020; Gensior & Kruijssen 2021). Thus, the question we aim to address in this paper is whether there is a single relationship between Σ_SFR, Σ_⋆ and Σ_mol which is valid in all environments, or if this relation differs between (some) environments. In the latter case, can we also aim at identifying the parameter(s) setting the level of Σ_SFR, besides Σ_⋆ and Σ_mol.

To minimize the covariance between model parameters, we normalize the distributions of the three involved variables, logΣ_SFR, logΣ_⋆ and logΣ_mol, by their (full sample) mean before fitting the data (i.e., centering). After the data are normalized, we find the best set of parameters C_⋆, C_mol and C_norm such that:

$$\begin{array}{cc}\hfill log{\mathrm{\Sigma }}_{\mathrm{SFR}}-⟨log{\mathrm{\Sigma }}_{\mathrm{SFR}}⟩=& C_{\star }(log{\mathrm{\Sigma }}_{\star }-⟨log{\mathrm{\Sigma }}_{\star }⟩)\hfill \\ \hfill & +C_{\mathrm{mol}}(log{\mathrm{\Sigma }}_{\mathrm{mol}}-⟨log{\mathrm{\Sigma }}_{\mathrm{mol}}⟩)+C_{\mathrm{norm}},\hfill \end{array}$$(4)

where the quantities inside the ‘⟨ ⟩’ brackets represent the mean of each distribution, across the full sample. Collecting these terms, we can define a new recentered normalization ${\widehat{C}}_{\mathrm{norm}}$ so that:

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathrm{norm}}=C_{\mathrm{norm}}+⟨log{\mathrm{\Sigma }}_{\mathrm{SFR}}⟩-C_{\star }⟨log{\mathrm{\Sigma }}_{\star }⟩-C_{\mathrm{mol}}⟨log{\mathrm{\Sigma }}_{\mathrm{mol}}⟩.\end{array}$$(5)

For consistency we apply a fixed normalization to the data (⟨logΣ_⋆⟩ = 8.11, ⟨logΣ_mol⟩ = 7.25 and ⟨logΣ_SFR⟩= − 2.48) throughout this work, that is, we use the same average values at different spatial resolutions (Sect. 4.2), and when exploring environments in individual galaxies (Sect. 5.2).

We choose to fit a single C_norm value for all different galactic environments to avoid degeneracies between the normalization (C_norm) and the slopes of the planes (C_⋆ and C_mol). This choice focuses our analysis specifically on variations of the dependency on Σ_⋆ and Σ_mol. However, we stress that the main conclusions of this paper are not affected by this approach.

We designed our methodology (i.e., unique normalization and detection threshold) to test for environment-driven differences in the plane relating logΣ_⋆, logΣ_SFR, and logΣ_mol. The existence of such differences is of considerable interest, because it can reveal the degree to which physics covariant with the environment definitions impact the star formation process. However, we caution that our methodology is not necessarily the optimal approach to determine the ‘real’ slopes of the star-forming plane. Environment may also affect the normalization (C_norm). Because we keep C_norm fixed, our fitting method may allow some of the signal associated with normalization variation to the derived slopes (C_⋆ and/or C_mol).

In order to find the best-fitting parameters that describe the star-forming plane in each environment, as defined in Eq. (3), we perform a Bayesian hierarchical linear regression, based on routines from the PYMC3 python package (Salvatier et al. 2016). The hierarchical linear regression represents the middle ground between assuming that different environments are completely independent populations, and assuming that they all are identical and described by the same model. Instead, it assumes that the parameters of the models that describe the data sets of different galactic environments have some underlying similarity. That is, the coefficients that describe the star-forming plane in different galactic environments are assumed to follow the same underlying ‘hyperprior’ distribution, defined by a common set of ‘hyperparameters’. We choose weakly informative hyperpriors in order to ensure that our results are not dominated by the priors. The hyperprior distributions considered are the following:

$$\begin{array}{cc}& C_{\star ,\mu }\sim \mathcal{N}(1,{10}^2),\hfill \end{array}$$(6)

$$\begin{array}{cc}& C_{\star ,\sigma }\sim \mathcal{H}(5^2),\hfill \end{array}$$(7)

$$\begin{array}{cc}& C_{\mathrm{mol},\mu }\sim \mathcal{N}(1,{10}^2),\hfill \end{array}$$(8)

$$\begin{array}{cc}& C_{\mathrm{mol},\sigma }\sim \mathcal{H}(5^2),\hfill \end{array}$$(9)

where 𝒩(μ, σ²) and ℋ(σ²) stand for Normal and Half-Normal distributions, respectively, with mean μ and variance σ². The prior distribution of each coefficient is then defined as

$$\begin{array}{c}\hfill C_j\sim \mathcal{N}(C_{j,\mu },C_{j,\sigma }),\end{array}$$(10)

for j in { ⋆ ,mol}. In addition to these hyperpriors, we adopt the following priors for the normalization coefficient.

$$\begin{array}{cc}& C_{\mathrm{norm},\mu }\sim \mathcal{N}(0,{10}^2),\hfill \end{array}$$(11)

$$\begin{array}{cc}& C_{\mathrm{norm},\sigma }\sim \mathcal{H}(5^2).\hfill \end{array}$$(12)

Finally, we include an additional term to account for intrinsic dispersion, which is common for all environments, and its prior is defined as

$$\begin{array}{c}\hfill ϵ\sim \mathcal{N}(0,{\sigma }_{\mathrm{intr}}),\end{array}$$(13)

where σ_intr corresponds to the intrinsic dispersion of the data with respect to the model, and for its prior distribution we use:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{intr}}\sim \mathcal{HC}(5^2),\end{array}$$(14)

where ℋ𝒞(σ²) stands for Half-Cauchy distribution, a common choice for a prior distribution of a scaling parameter like the intrinsic scatter. We find that σ_intr describes well the standard deviation in the fit residuals of our data.

We use four Markov chain Monte Carlo (MCMC) chains to sample the posterior distributions, each one of them having 2000 iterations plus 2000 additional burn-in iterations. The posterior is sampled using the NUTS algorithm (Homan & Gelman 2014). The convergence of the MCMC chains is ensured by using the $\widehat{R}\approx 1$ criterion (Gelman & Rubin 1992). Essentially, $\widehat{R}$ corresponds to the ratio of the between-chain variance and the within-chain variance. In Appendix A, we use a toy model to test the accuracy of the hierarchical modeling, and find that this approach correctly recovers the test values in mock data sets.

4. Results

4.1. Star-forming plane across galactic environments at 150 pc resolution

Figure 2 shows the posterior distributions of each one of the parameters that define the star-forming plane in each galactic environment, following the same color code as in Fig. 1. The mean and standard deviation of each one of these distributions are reported in Table 2, as well as the number of pixels used in the fit for each environment. We firstly note that the marginalized posteriors of C_⋆ and C_mol of different environments are significantly different (> 1σ in most of cases). The posterior distributions of C_⋆ and C_mol are significantly broader for galaxy centers than for the rest of environments. This could either be statistical (i.e., due to the smaller number of pixels probing this environments) or intrinsic (i.e., different galaxy centers have different scaling relations). The posterior of C_norm is only shown for the full sample, since we model a single normalization for all environments as explained in Sect. 3.2. As a single value of C_norm is fitted for the full sample, the posterior distribution of C_norm is considerably narrower than that of C_⋆ and C_mol. Thus, the x-axis in the right-most panel has smaller bins, and the distribution of C_norm has been renormalized for an easier visualization. After correcting for the centering of the data (Eq. (5)), this value of C_norm implies an absolute normalization ${\widehat{C}}_{\mathrm{norm}}$ of −8.83. This value represents a characteristic depletion time of ∼1.5 Gyr, well within the range of values reported in Querejeta et al. (2021).

Fig. 2. Posterior distributions for the coefficients C⋆, Cmol and Cnorm that define the star-forming plane in each galactic environment, measured at a spatial resolution of 150 pc. As the posterior distribution of Cnorm is considerably narrower than that of C⋆ and Cmol, the x-axis has been binned in smaller bins and the distribution renormalized for an easier visualization.

We note that C_⋆ parametrizes the rate of change of the rKS normalization, associated with Σ_⋆, and the index of the power law spans a wide range of values, changing sign significantly across the sample, from about −0.41 to 0.45. Interestingly, bars have a very negative value of C_⋆ (∼ − 0.4), meaning that they show longer depletion times at higher Σ_⋆ values. While this low value of C_⋆ is certainly an indication of lower star formation efficiencies, compared to other environments, we remind the reader that due to our adopted approach of using a single normalization across different galactic environments, a negative C_⋆ could be interpreted also as the slope trying to capture a vertical offset between environments.

In Sect. 5.2, we discuss the possibility of this low value of C_⋆ being linked to lower star formation efficiencies (SFE) driven by radial and/or turbulent motion of gas in bars. In contrast, spiral arms show the highest values of C_⋆ (∼0.45), leading to shorter depletion times at higher Σ_⋆ values. As Σ_⋆ varies strongly as a function of radius, C_⋆ expresses a radial trend in the normalization of the molecular gas scaling relation. Analogous variations were reported in Muraoka et al. (2019), where the authors report radial variations of a factor 2 − 3 in the SFE of individual nearby galaxies from the CO Multi-line Imaging of Nearby Galaxies (COMING) project (Sorai et al. 2019).

On the other hand, C_mol shows a more homogeneous behavior with values in the range of 0.83 to 1.01, which implies that higher Σ_mol leads to higher Σ_SFR everywhere in the galaxy. Subtle variations are found for bars having the lowest values (C_mol ∼ 0.83) or rings having the largest value (C_mol ∼ 1.01). It is interesting to note that while rings show nearly a linear slope in molecular gas, other galactic environments exhibit a sublinear behavior. While a linear relation implies a constant depletion time (defined as τ = Σ_mol/Σ_SFR), a sublinear slope leads to a depletion time that increases with gas surface density (neglecting changes in C_⋆). Shetty et al. (2014) examined different possibilities to physically interpret this sublinearity. On one hand, differences in the molecular gas properties such as star formation efficiency or volume density (see Bacchini et al. 2019, 2020) could lead to variations of the depletion time. Alternatively, variations in the depletion time could be produced by a diffuse component of the molecular gas (i.e., molecular gas not actively forming stars). Schinnerer et al. (2019) and Pan et al. (2022) found that a significant fraction of molecular gas is indeed decoupled from the Hα emission in PHANGS galaxies. Although for this analysis, we choose only pixels that show both molecular and ionized gas emission, quiescent gas could still be present in our data due to projection effects. In this scenario, the sublinearity is an indication that the fraction of diffuse molecular gas grows with Σ_mol. We highlight that this sublinearity persists for bars, disks and spiral arms, when we let C_norm free for each individual environment.

A third possibility is that variations in the depletion times are driven by variations in the real underlying α_CO conversion factor or CO(2−1)-to-CO(1−0) line ratios. Leroy et al. (2022) studied variations in the R₂₁ ≡ CO(2−1)-to-CO(1−0) line ratio within PHANGS galaxies, and reported a central enhancement of R₂₁ of a factor of ∼0.18 dex. Such anticorrelation of R₂₁ with galactocentric radius (and thus, correlation with Σ_SFR) implies that our fiducial choice of R₂₁ = 0.65 could potentially lead to an overestimation of the total molecular gas content in the inner and denser regions of the galaxy, and to an underestimation of the gas content in the outer regions. This can effectively depress the measured C_mol slope by up to ∼0.15. This effect likely plays a more relevant role in the global measurement, considering all environments, due to the larger dynamic range of Σ_⋆. However, it can not be ruled out that it is also, partially, driving the sublinearity found in the single-environment measurements.

Figure 3 shows the covariance between the posterior distributions of the parameters C_⋆ and C_mol measured for each galactic environment. As C_⋆ and C_mol are the coefficients of the linear combination to predict logΣ_SFR of a given pixel from its logΣ_⋆ and logΣ_mol measurements, they show a negative covariance, that is, a higher C_⋆ implies a lower value of C_mol. The covariance is much larger for centers likely due to the lower number of pixels in these environments (see Table 2), but could also be connected to a physically more heterogeneous behavior of centers. On the other hand, the covariance in disks is significantly smaller. In Appendix A, we show that this covariance arises, in part, as an artifact of the fitting procedure. However, the covariance could be also partially due to the fact that Σ_⋆ and Σ_mol are physically highly correlated quantities, through the rMGMS, thus the covariance reflects the slope between them.

Fig. 3. Posteriors distributions of the coefficients C⋆ and Cmol for each separate environment, measured at an spatial resolution of 150 pc, to show the covariance between them. The posteriors have been smoothed using a Gaussian kernel. The color scale for each environment indicates the 1-, 2- and 3-sigma confidence intervals.

Figure 4 provides an alternative visualization of the 2D planes obtained with the fitting procedure. It shows the partial residual for each independent variable after removing the dependency on the second independent variable. Specifically, we define δ_⋆ and δ_mol as

$$\begin{array}{cc}& {\delta }_{\star }=log{\mathrm{\Sigma }}_{\mathrm{SFR}}-C_{\mathrm{norm}}+C_{\mathrm{mol}}log{\mathrm{\Sigma }}_{\mathrm{mol}},\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\delta }_{\mathrm{mol}}=log{\mathrm{\Sigma }}_{\mathrm{SFR}}-C_{\mathrm{norm}}+C_{\star }log{\mathrm{\Sigma }}_{\star },\hfill \end{array}$$(16)

Fig. 4. Partial residual for each independent variable, dropping the dependency from the second independent variable, found for each individual environment and for all environments together as a function of Σ⋆ (left) and Σmol (right), measured at an spatial resolution of 150 pc. The binned data are shown for each corresponding environment. The trend for the full sample is overplotted as a black dotted line.

to visualize the partial residuals as a function of Σ_⋆ and Σ_mol, respectively. The figure shows in each row the residuals for a given environment as a 2D histogram. The bottom row shows the residuals for the full sample. The slope in each panel corresponds to the best-fitting value of C_⋆ (left column) and C_mol (right column) for each environment.

The sharp cut at low Σ_⋆ in the x-axis of the left panels is a consequence of the adopted detection threshold for each galactic environment as described in Sect. 3.2. The figure shows the range of slopes derived for the different environments probed. Spiral arms are the environment that show the most positive trend for Σ_⋆, while centers and disks show more subtle positive trends. In contrast, bars show the most negative C_⋆ values of all environments. Similarly, bars also exhibit the lowest values of C_mol, which can be understood as a lower efficiency at converting molecular gas into stars toward higher values of Σ_mol. On the other hand, star-forming rings show the largest C_mol value (≈1.01), which indicates a nearly constant depletion time across the Σ_mol range. We inspect these trends per environment and find that they are not driven by individual galaxies and that, indeed, the trends of the multiple galaxies are consistent with each other. These trends reveal that the galaxy-to-galaxy variations in the scaling relations reported in Ellison et al. (2021) and Pessa et al. (2021) could plausibly be explained by a different relative contribution of the different galactic environments across the sample galaxies.

4.2. Effect of spatial resolution

In this section, we explore how the spatial scale of the data impacts our measurements. We degraded our data, as explained in Sect. 2.6, and repeated the measurement of C_⋆, C_mol and C_norm, following the same procedure as detailed in Sect. 3.2. Figure 5 displays the mean of the distribution of each parameter, for each galactic environment, as a function of the spatial scale of the data. The error-bars show the standard deviation of the corresponding posterior distribution. The bottom panel shows the ‘real’ normalization value ${\widehat{C}}_{\mathrm{norm}}$, as defined in Eq. (5). It is clear that different environments have different coefficients in their scaling relations at ∼150 pc resolution, but the differences are reduced when looking at larger spatial scales. This is likely due to a combination of two effects: (i) At larger spatial scales, the light from different environments is blended, and (ii) the number of available pixels decreases drastically at larger spatial scales, from a total sample size of ∼50 000 to ∼1200 making the measurement statistically less certain. Whereas (ii) would have an effect across the whole range of spatial scales probed, (i) will be relevant only at spatial scales larger than that of the typical structural size of each environment.

Fig. 5. Change of the coefficients that define the star-forming planes in each galactic environment with spatial scale. The color code for each environment is the same as used in Fig. 1. The cyan diamond and a magenta hexagon indicate the results reported by Sánchez et al. (2021) (EDGE–CALIFA sample) and Lin et al. (2019), respectively.

Figure 6 explicitly shows the posterior distributions of C_⋆ and C_mol, for each environment and spatial scale probed. It shows that at spatial scales larger than 300 pc, the differences between the posteriors of different environments become smaller than 1σ (with respect to the full-sample measurement), especially for C_mol. This spatial scale roughly matches the width of the spiral arms in our Galaxy (Reid et al. 2014). In our sample, rings, spiral arms and centers are found to have sizes on the order of a few hundreds parsecs (Querejeta et al. 2021), and therefore, their emission is expected to blend at spatial scales larger than this. On the other hand, bars have larger typical sizes, on the order of several kpc, which is consistent with finding that differences in the posterior distributions of their scaling relation parameters persist up to scales of ∼500 pc in both, C_⋆ and C_mol. This homogeneization toward lower spatial resolutions agrees well with the absence of systematic galaxy-to-galaxy variations reported in Sánchez et al. (2021), using data from the CALIFA (Sánchez et al. 2012) and EDGE (Bolatto et al. 2017) surveys.

Fig. 6. Posterior distribution of the coefficients that define the star-forming planes in each galactic environment, at each one of the spatial scales probed, from 150 pc to 1 kpc. The spatial scale and the number of pixels used for each measurement are indicated in the top left corner for each row.

It is also worth noting the general trends of C_⋆ and C_mol toward larger spatial scales. While the posteriors of C_⋆ are shifted toward more negative values, the posteriors of C_mol move toward steeper slopes. These changes are likely a combination of the effect of N/Ds in our data (even though we remove the Σ_⋆ ranges more critically affected by N/Ds to minimize their impact), and the C_⋆ − C_mol covariance. The former could induce a steepening in these slopes, while the latter would drive their (a)symmetry.

Finally, we find that galactic environments are not only different in terms of the coefficients describing their star-forming plane, but also in terms of the scatter around it. The top panel of Fig. 7 shows the scatter (defined as the standard deviation of the residuals) of each environment with respect to its modeled star-forming plane, at the different spatial scales probed (solid lines).

Fig. 7. Comparison between plane and single power law models. Top: scatter with respect to the star-forming 2D plane (solid) and rKS (dashed) determined for each galactic environment, at the different spatial scales probed. Bottom: difference of the LOO score measured for the plane and the rKS models, where a > 0 indicates a preference for the plane model. The error bars are the 3σ standard errors. The black horizontal line at 0 marks where there would be no preference between the models. The plane model is preferred over the rKS at all spatial scales, increasingly so at finer spatial scales.

The first noticeable feature is how the scatter decreases toward larger spatial scales. The same trend is reported in Pessa et al. (2021) for the 1D scaling relations. This effect has been associated with the decoupling of different stages of the star-forming cycle at high spatial resolution, in other words, a given aperture can be either dominated by a peak in the CO emission (i.e., early in the star-forming cycle) or a peak in the Hα emission (i.e., late in the star-forming cycle) (e.g., Schruba et al. 2010; Feldmann et al. 2011; Leroy et al. 2013; Chevance et al. 2020a,b). At larger spatial scales, these peaks are averaged, which diminishes the sampling effect. (Schruba et al. 2010; Kruijssen & Longmore 2014; Semenov et al. 2017; Kruijssen et al. 2018). The fraction of pixels containing only molecular gas or only ionized gas has been empirically characterized in Schinnerer et al. (2019) and Pan et al. (2022). Both works found that ∼500 pc is a critical spatial scale, below which distinguishing particular stages of the star formation process is possible, in nearby star-forming galaxies. In this line, Kruijssen et al. (2019) and Chevance et al. (2020a) found that the typical distance between independent star-forming regions is 100 − 300 pc. Kim et al. (in prep.) and Machado et al. (in prep.) report an average distance of 250 − 300 pc in the full PHANGS sample.

The figure also shows that at ∼150 pc, the scatter in certain environments is significantly lower than in others. Centers and spiral arms show a particularly low scatter, whereas disk is the only environment that shows a scatter larger than the full sample. It is also worth noticing that the rate at which the scatter decreases with increasing spatial scale also varies across environments, with bars showing the flattest trend between 150 and 300 pc, while rings and centers showing the steepest trend. Together with the scatter, this indicates variations in the properties of the molecular clouds, such as differences in their lifetimes, or a different typical separation between clouds (as suggested by Chevance et al. 2020a). In this line, Henshaw et al. (2020) found evidence that fragmentation indeed occurs on smaller size scales in centers and potentially also (but to a lesser extent) in spiral arms. This idea would be consistent with the higher scatter, and its flatter decrease toward larger spatial scale that we see in bars and disks, resulting from a lower spatial density of star-forming regions, as compared to rings, spiral arms or centers. Additional variations in the intrinsic properties of molecular clouds have been reported in previous works (Sun et al. 2018, 2020b,a; Rosolowsky et al. 2021). Sun et al. (2020b), using data from PHANGS–ALMA, found that cloud-scale molecular gas surface density, velocity dispersion and turbulent pressure of molecular clouds depend on local environmental conditions. Similarly, Colombo et al. (2014) used data from the PdBI Arcsecond Whirlpool Survey (PAWS; Schinnerer et al. 2013) to characterize variations of molecular gas properties across galactic environments, and concluded that this environmental variations are a consequence of the combined action of large-scale dynamical processes and feedback from high-mass star formation. A similar conclusion was obtained by Renaud et al. (2015), where the authors used hydrodynamical simulations to study the regulation of the star formation in bars.

4.3. Plane versus single power law

Here we explore whether a 2D plane offers a more accurate prediction of Σ_SFR than the rKS. For this comparison, we find the best-fitting power law in each environment following the same prescription described in Sect. 3.2 (i.e., centering the data and using a single C_norm value), dropping the dependency on Σ_⋆ (i.e., removing C_⋆ from the model). The top panel of Fig. 7 shows as dashed lines the scatter with respect to the rKS measured for each galactic environment, at each one of the spatial scales probed. The figure shows that when we separate the data by environment, the additional dependency of Σ_⋆ does not reduce the scatter significantly. Moreover, rKS shows a slightly smaller scatter than the 2D plane in some environments. This is a consequence of the low predictive power found for Σ_⋆, with respect to Σ_mol (i.e., C_mol > C_⋆), to predict Σ_SFR. However, it should be kept in mind that the lower predictive power found for Σ_⋆ could be, at least partially, related to its limited dynamic range, as shown by Fig. 4. Nevertheless, this finding is consistent with Lin et al. (2019), where the authors find that the ‘extended’ version of the rKS (i.e., with an additional dependency of Σ_⋆) has its scatter only slightly reduced with respect to the conventional rKS.

Despite offering only a modest improvement in the scatter, we find that the plane model is formally favored over the rKS at all spatial scales. We compare the models using leave-one-out (LOO; Sammut & Webb 2010) cross-validation, a standard model comparison method (Vehtari et al. 2017). The LOO cross validation is computationally more expensive than other forms of cross validations (e.g., K-fold, random subsampling), but it offers a more robust estimate of the predictive accuracy of a model, and it is suitable for relatively small datasets.

The bottom panel of Fig. 7 shows the difference of the LOO statistics for both models (the environmental separation is included in each model), where a higher LOO statistic indicates the preferred model. The figure shows that the LOO statistic for the plane model is higher than for the rKS model at all spatial scales at the > 3σ level, and that the preference for the plane model becomes stronger at high spatial resolution. For this comparison, we checked that the measured σ_intr at each spatial scales captures the different levels of scatter of the data across the range of spatial scales probed, making the decrease in the LOO statistic not primarily driven by the intrinsic scatter.

This result is consistent with Fig. 6, which shows that the width of the posterior distributions of C_⋆ increases toward larger spatial scales in all environments, and thus nonzero C_⋆ are more constrained on finer spatial scales. This means that including the C_⋆ in the model is more important at smaller scales, but less critical at larger spatial scales. An explanation for this is that in the high spatial resolution measurements, the high Σ_SFR (and Σ_mol) regions are often located in the inner regions of the galaxy (and vice-versa). However, at lower spatial resolutions, Σ_SFR (and Σ_mol) is diluted in a larger regions, due to its intrinsic patchy configuration. As this does not occur with Σ_⋆, this results in Σ_⋆ providing less information to predict Σ_SFR, compared to Σ_mol, toward lower resolutions.

5. Discussion

5.1. Comparison to previous findings

Here we compare our results with those in the literature. Our results agree relatively well with those reported in Lin et al. (2019), where the authors used ALMaQUEST (Lin et al. 2020) data to explore the extended version of the rKS relation (Shi et al. 2011, 2018). In terms of the coefficients we use in this paper, they found C_⋆ = −0.29, C_mol = 0.97 and C_norm = −8.27. These numbers agree with those we report here (see Table 2), not only in their absolute values but also in terms of the relative predicting power of logΣ_⋆ and logΣ_mol to predict logΣ_SFR. Shi et al. (2018) performed similar measurements using data from The H I Nearby Galaxy Sample (THINGS; Walter et al. 2008) and from the GALEX data archive (using total gas surface density rather than molecular gas). Interestingly, in terms of the same coefficients, they found C_⋆ = 0.55, C_mol = 1.09 and C_norm = −10.47. This is roughly consistent with the idea that mid-plane pressure helps to regulate current star-formation, as pressure scales as ${\mathrm{\Sigma }}_{\star }^{0.5}{\mathrm{\Sigma }}_{\mathrm{gas}}$ (Shi et al. 2018). Furthermore, they found that outer disks of dwarfs galaxies and local luminous infrared galaxies show the largest offset in this relation, whereas local spirals show almost no offset. This is an intriguing result, as we also find a relatively similar behavior for spiral arms. However, we stress that this comparison has the caveats of (i) local spiral galaxies host different environments (not only spiral arms) and (ii) the study carried out in Shi et al. (2018) uses data with different spatial resolutions. Nevertheless, this is an interesting comparison that could indicate that mid-plane hydrostatic pressure plays a more relevant role in spiral arms than in other environments.

On the other hand, in Sánchez et al. (2021), the authors performed a similar measurement using data from the CALIFA and EDGE surveys. For the EDGE–CALIFA sample, which represents a better reference for comparison in terms of molecular gas tracer and spatial scale, they found C_⋆ = 0.44 ± 0.01, C_mol = 0.57 ± 0.01 and C_norm = −9.5 ± 0.01. These coefficients partially agree with our results, in terms of the lower relative importance of logΣ_⋆ to predict the logΣ_SFR value with respect to logΣ_mol. However, significant quantitative differences between the measured coefficient exist. The differences persist even if we exclude the detection-fraction-threshold step from the fitting (this does not change our results significantly at a spatial resolution of ∼1 kpc).

Part of the differences we see here are in the statistical questions we are exploring with respect to the variation of these scaling relationships with environment. Sánchez et al. (2021) demonstrate the importance of the statistical framework for interpreting the results of these regressions. Our hierarchical Bayesian approach for this analysis provides a framework for analyzing the data within well-resolved environments that could not be addressed at coarser resolution by the complementary EDGE–CALIFA sample, despite its much larger size, compared to our PHANGS–MUSE sample. A direct comparison is not possible since PHANGS–MUSE galaxies all have distances lower than 20 Mpc, and an overlap between both samples does not exist. Nevertheless, it is worth mentioning that Sánchez et al. (2021) also report measurements using their CALIFA–only sample, where the molecular gas estimates are based on the ISM dust attenuation prescription reported in Barrera-Ballesteros et al. (2021a). For this sample, they obtained C_⋆ = 0.66 ± 0.02, C_mol = 0.38 ± 0.01 and C_norm = −9.9 ± 0.05. These numbers are in stronger disagreement with our results, not only due to the differences in the values of the measured slopes, but also in the relative importance of logΣ_⋆ and logΣ_mol to predict the logΣ_SFR Hence, any comparison of results should be made carefully, as differences in the sample, such as the molecular gas tracer, could lead to different conclusions. The results reported in Sánchez et al. (2021) (EDGE–CALIFA sample) and in Lin et al. (2019) are indicated in Fig. 5 as a cyan diamond and a magenta hexagon, respectively.

Other studies have also explored systematic difference across galactic environments. In Querejeta et al. (2021), the authors used 74 galaxies of the PHANGS sample to measure the distributions of Σ_mol, Σ_SFR and depletion times across galactic environments at a fixed spatial scale of ∼1.5 kpc. They found a strong correlation between molecular gas and SFR surface densities, with a global slope of N = 0.97, but little variation across galactic environments. However, they reported a slight offset toward shorter depletion times for centers (τ_dep = 1.2 Gyr), and longer depletion times for bars (τ_dep = 2.1 Gyr). This result is consistent with our measurement of the lowest C_⋆ in this environment at a spatial resolution of ∼1 kpc. On the other hand, we do not find that C_⋆ in centers is higher than in the other environments. However, this specific measurement is highly uncertain (due to the low number of ‘center’ pixels at this resolution in our study) and thus, a proper comparison is not possible. Overall, they did not find evidence of strong variations in how efficiently different environment form stars. Not finding strong differences across environments is also consistent with the homogenization of galactic structure at large spatial scales which we report in this work (see Sect. 4.2).

5.2. What drives the environmental variations we see?

In previous sections, we reported the coefficients that describe the star-forming plane for each individual galactic environment. Furthermore, we find significant differences in the star-forming planes associated with these different environments, particularly for bars, spiral arms and rings. In this section, we explore which parameter(s) could be driving the differences we report in Sect. 4.

To perform this exploration, we measure C_⋆, C_mol and C_norm in each environment of each individual galaxy in our sample (following the same procedure detailed in Sect. 3.2), and search for correlations between the coefficients obtained in each environment, and a set of additional parameters:

• Molecular gas fraction; logf_mol, calculated as log(Σ_mol/Σ_⋆)

• Star formation efficiency; logSFE, calculated as log(Σ_SFR/Σ_mol)

• Free fall time (t_ff), calculated following Utomo et al. (2018), i.e., $t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G}(\frac{H}{{\mathrm{\Sigma }}_{\mathrm{mol}}})}$, where G is the gravitational constant and H is the vertical scale height of the molecular gas layer which can be estimated as $H=\sqrt{\frac{{\sigma }_{\mathrm{mol}}^2{h}_{\star }}{G{\mathrm{\Sigma }}_{\star }}}$, with σ_mol corresponding to the velocity dispersion of the molecular gas component, and h_⋆ is the typical stellar scale height, for which we adopt a value of 300 pc (Utomo et al. 2018).

• Depletion time in units of free fall time; τ/t_ff, where τ = Σ_mol/Σ_SFR and t_ff is calculated as described above.

• Mid-plane hydrostatic pressure; logP_h, calculated following Elmegreen (1989) as $P_{\mathrm{h}}=\frac{\pi }{2}G{\mathrm{\Sigma }}_{\mathrm{mol}}({\mathrm{\Sigma }}_{\mathrm{mol}}+\frac{{\sigma }_{\mathrm{mol}}}{{\sigma }_{\star }}{\mathrm{\Sigma }}_{\star })$. This approach neglects the contribution from atomic gas, however, the relative contribution of atomic gas with respect to molecular gas in these galaxies and this regime of Σ_mol values is expected to be subdominant (Bigiel et al. 2008; Schruba et al. 2019; Leroy et al., in prep.). Barrera-Ballesteros et al. (2021b) explore the role of P_h in regulating Σ_SFR, and they found a tight (scatter ∼0.2 dex) correlation between these two quantities. Sun et al. (2020a) also found that Σ_SFR correlates with the dynamical equilibrium pressure of the ISM. This is encouraging to explore if differences between environments could be driven by differences in the mid-plane hydrostatic pressure, which has the additional dependencies on σ_mol and σ_⋆.

• Gas-phase metallicity; [Z/H]_gas, modeled from emission-line measurements in the MUSE data, allowing for azimuthal variations, as detailed in Williams et al. (2022).

• Stellar velocity dispersion; σ_⋆, resulting of SSP fitting of the MUSE data (see Sect. 2.4).

• Hα velocity dispersion; σ_Hα, measured from the MUSE data (see Sect. 2.5).

• Molecular gas velocity dispersion; σ_mol, measured from the ALMA data (Leroy et al. 2021b).

We have additionally considered stellar population parameters measured by SSP fitting of the MUSE data (Emsellem et al. 2022; Pessa et al., in prep.).

• light-weighted stellar age; logAGE_LW

• mass-weighted stellar age; logAGE_MW

• light-weighted stellar metallicity; [Z/H]_LW

• mass-weighted stellar metallicity; [Z/H]_MW.

A similar exploration of parameters was carried out by Dey et al. (2019), using data from the EDGE–CALIFA sample, performing a data-driven approach to investigate what shapes the SFR, finding that Σ_SFR scales primarily with Σ_⋆ and Σ_mol. Conversely, we perform this exploration in order to understand what is driving the environmental differences we see in the plane spanned by these three quantities.

We compute the mean of the distribution of these quantities in each environment (considering only those pixels that were used for the fitting of the plane), and measure the weighted Pearson’s correlation coefficient (ρ) between these environmental-averaged quantities and the corresponding coefficients. For this latter step, we consider only the subset of galaxies that satisfy the following criteria:

1. Probe at least three different environments.

2. Host galaxy has a bar and spiral arms.

3. Fraction of pixels removed due to the imposed detection threshold is < 50%.

These selection criteria ensure that we select galaxies in which we can simultaneously probe several environments, including bars and spiral arms, which are the environments that exhibit the largest differences, especially in terms of C_⋆ (see Fig. 2), and that most of the pixels of these galaxies are actually used in the fitting. Out of our sample, 8 galaxies satisfy the listed conditions (NGC 1300, NGC 1365, NGC 1512, NGC 1566, NGC 1672, NGC 3627, NGC 4303, NGC 4321). As a result of these conditions, we drop from our sample galaxies with total stellar mass of logM_⋆ [M_⊙]< 10.6. We perform this exploration on an individual-galaxy basis because we acknowledge that, while different environments usually exhibit different properties in a given galaxy, the same environments are not necessarily identical across different galaxies.

For this subset of galaxies, we quantify the level of correlation between each of the parameters listed, and the coefficients C_⋆ and C_mol obtained for each environment, using the ‘overall’ Pearson’s correlation coefficient, defined as:

$$\begin{array}{c}\hfill ⟨\rho ⟩=\frac{|{\sum }_{i=1}^{n_{\mathrm{gal}}}{\rho }_i|}{n_{\mathrm{gal}}},\end{array}$$(17)

which corresponds to the average ρ across the subset of galaxies. For each galaxy, the uncertainty in its correlation coefficient ρ_i (i.e., between the plane coefficients and each one of the parameters listed) is estimated by performing 50 Monte-Carlo simulations, perturbing the measured C_⋆ and C_mol for each environment in each iteration. The uncertainty in ⟨ρ⟩ is then calculated using standard error propagation.

Figure 8 shows the overall Pearson’s correlation coefficient, that quantifies the average level of correlation between the median of each one of the quantities listed on the x-axis for a given environment and its corresponding C_⋆ (blue) and C_mol (orange) value. It is clear that some parameters show a high level of correlation with C_mol, the highest being the stellar mass-weighted age of the underlying stellar populations. However, we acknowledge that a possible correlation with the average age of the underlying stellar population would likely be a consequence of differences in the star-forming process, rather than driving it. For this reason, we focus on the following high-ρ parameter, ⟨τ/t_ff⟩ (⟨ρ⟩∼0.77).

Fig. 8. Overall Pearson coefficient to quantify the level of correlation between the C⋆ and Cmol parameters that define the star-forming plane in each environment, and the set of different parameters explored.

As ⟨τ/t_ff⟩ corresponds to the depletion time normalized by the free fall time, that is, the characteristic time that it would take a gas cloud to collapse under its own gravitational attraction, it is a metric that describes how efficiently the gas is collapsing and forming stars. A high value can be interpreted as a relative excess of molecular gas with respect to SFR, in other words, a larger fraction of the gas is not forming stars (quiescent), compared to an environment with a shorter average depletion time. Figure 9 shows explicitly the correlation between ⟨τ/t_ff⟩ measured in each environment, and the corresponding C_mol coefficient for each individual galaxy in our subset. For most of the subset of galaxies, we see a negative trend (except for NGC 1365). Bars tend to have longer depletion times and lower values of C_mol, as opposed to spiral arms, that show shorter depletion times and higher values of C_mol. However, this correlation is not surprising, as C_mol quantifies how efficiently molecular gas forms stars, and ⟨τ⟩ corresponds to the inverse of ⟨logSFE⟩ in a given environment. Thus, environments with on average longer depletion times (or lower SFE) will be naturally described by a lower C_mol. Nevertheless, this correlation shows that quantifiable differences in depletion time across environments exists, and that these differences lead to variations in the star-forming plane.

Fig. 9. Cmol − ⟨τ/tff⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.

Therefore, it is interesting that other parameters that show a high level of correlation with C_mol are σ_⋆ and σ_Hα. Velocity dispersion encodes information about noncircular motion of gas and stars in the galaxy, such as radial motions along the bar, or turbulence induced by star formation feedback or by AGN activity, for galaxies hosting an AGN (NGC 1365, NGC 1566, NGC 1672, NGC 3627, NGC 7496). Thus, the correlations between C_mol and σ_⋆, Hα could be an imprint of how the noncircular motion of the gas prevents collapse and efficient star formation (leading to longer depletion times). The correlation with σ_Hα is presented in Appendix D for completeness.

Of the remaining parameters, only two show a potential correlation with C_mol, ⟨logSFE⟩ and ⟨f_mol⟩. The first one is not surprising, as it corresponds to the inverse of ⟨τ⟩. For ⟨f_mol⟩, we find hints of a positive correlation, that is, a higher gas fraction leads to a higher C_mol. This implies that stars are formed more efficiently in galactic environments which are more gas-rich.

For t_ff, logP_h, [Z/H]_gas, σ_mol, log[Z/H]_LW and log[Z/H]_MW, we do not find a correlation with C_mol.

On the other hand, C_⋆ shows in general lower levels of correlation for most of the parameters probed. This is not unexpected, due to the lower predictive power of Σ_⋆ to predict Σ_SFR as compared to Σ_mol. Nevertheless, we still see a correlation of C_⋆ with logSFE (⟨ρ⟩∼0.68). Figure 10 shows how environments with higher average star formation efficiency exhibit higher measured C_⋆ values. A similar correlation is found for ⟨f_mol⟩, which is shown in Appendix D for completeness.

Fig. 10. C⋆ − ⟨log SFE⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.

⟨τ/t_ff⟩ and logAGE_LW show hints of a mild correlation with C_⋆ (⟨ρ⟩∼0.5), but as discussed earlier, the first is expected as depletion time correspond to the inverse of SFE, and the second, if exists, is likely a consequence (rather than a cause) of differences in the star-forming process. None of the remaining parameters show hints of correlation with C_⋆.

The poor level of correlation of C_⋆ and C_mol with ⟨logP_h⟩ can be explained by the fact that P_h scales primarily with Σ_mol and Σ_⋆. These dependencies are already captured by C_⋆ and C_mol, thus, P_h contributes little additional information.

Finding these correlations with physical conditions in a given galactic environment suggests that an additional parameter is indeed causing the variations we measure across galactic environments (and also can explain galaxy-to-galaxy variations), and thus, playing a relevant role in regulating the SFR. However, determining what is (or are) the key parameter(s) is beyond the scope of this paper, as we would need a larger sample of galaxies that show a variety of galactic environments, and have, at the same time, high coverage of Hα and CO emission.

Finally, we stress here that due to the relatively low number of galaxies where we can reliably probe several galactic environments (8), and since each one of the correlations has only as many points as different environments the galaxy has, we can only speculate about the influence of an additional parameter in the coefficients of the star-forming plane, rather than robustly establishing a causal relation. Thus, we present this as a possible line of exploration, with intriguing dependencies that need confirmation with larger samples.

5.3. Choice of α_CO conversion factor

We have tested our main conclusions against a constant α_CO conversion factor of 4.35 M_⊙ pc⁻² (K km s⁻¹)⁻¹ (Bolatto et al. 2013), instead of the metallicity-dependent prescription (described in Sect. 2). We found that qualitatively, none of our conclusions are affected by our choice of α_CO. Figure 11 shows the posterior distributions of the coefficients determined at a spatial scale of 150 pc, under the assumption of a constant α_CO conversion factor. It is easy to see that even though some specific coefficients may change, the relative difference between environments, and the relative predictive power of C_⋆ and C_mol remains qualitatively the same.

Fig. 11. Posterior distributions for the coefficients C⋆, Cmol and Cnorm that define the star-forming plane in each separate environment, using a constant αCO. The vertical dashed lines show the centers of the posteriors distribution obtained under our fiducial choice of αCO. As the posterior distribution of Cnorm is considerably narrower than that of C⋆ and Cmol, the x-axis has been binned in smaller bins and renormalized for an easier visualization.

6. Summary

We have investigated the star-forming plane, conformed by Σ_SFR, Σ_⋆, and Σ_mol in 18 galaxies from the PHANGS sample at a physical resolution of 150 pc, and explored potential variations driven by galactic environment. Our main conclusions are as follows:

1. We found significant differences (> 1σ) in the coefficients that describe the star-forming plane across galactic environments. These differences are particularly significant for bars, spiral arms and rings. We interpret these variations as evidence for an additional regulator mechanism that is not captured by neither Σ_⋆ nor Σ_mol.

2. These variations between environments homogenize toward lower spatial resolutions. Differences are no longer significant at spatial scales larger than ∼500 pc. This is because the combined effect of blending of environments at lower spatial resolutions, and the lower number of pixels per environment at large spatial scales.

3. We find a good agreement with similar measurements done using data from the ALMaQUEST survey. On the other hand, we find a moderate agreement with measurements reported using data from the EDGE–CALIFA sample. Differences with coarser resolution studies are not surprising, as our statistical framework is specifically designed to quantify differences that would be measurable at high spatial resolution only.

4. We used a subset of galaxies from our sample, where we can probe (≥3) several galactic environments, and search for correlations between variations in the coefficients that define the star-forming plane in each environment, and the median of the distribution of a number of variables measured in those environments. We find a strong correlation of C_mol with depletion time and Hα velocity dispersion (⟨ρ⟩∼0.77 and 0.72, respectively), which could be an imprint of longer depletion times driven by increased turbulent or radial motion in some galactic environments, leading to the environmental differences we observe for the star-forming plane.

Our results are consistent with the existence of additional physics being at play in the regulation of the star formation. However, a larger sample of galaxies where we can simulataneously probe different galactic environments is required to confirm the correlations explored in this work and, thus, to provide a definitive answer as to what is (or are) the additonal parameter(s) modulating the formation of stars.

Acknowledgments

This work was carried out as part of the PHANGS collaboration. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes IDs 094.C-0623, 098.C-0484 and 1100.B-0651. This paper also makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.01161.S, ADS/JAO.ALMA#2015.1.00925.S, ADS/JAO.ALMA#2015.1.00956.S and ADS/JAO.ALMA#2017.1.00886.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. M.C. and J.M.D.K. gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through an Emmy Noether Research Group (grant number KR4801/1-1), as well as from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). J.G. gratefully acknowledges financial support from the Swiss National Science Foundation (grant no CRSII5_193826). K.K. gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, P.I.: Kreckel). E.W.K. acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. E.R. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2017-03987. M.Q. acknowledges support from the Spanish grant PID2019-106027GA-C44, funded by MCIN/AEI/10.13039/501100011033. G.A.B. gratefully acknowledges support by the ANID BASAL project FB210003. J.Pe. acknowledges support by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP, co-funded by CEA and CNES. J.Pe acknowledges support from the ANR grant ANR-21-CE31-0010. R.S.K. and S.C.O.G. acknowledge support from DFG via the Collaborative Research Center (SFB 881, Project-ID 138713538) “The Milky Way System” (subprojects A1, B1, B2 and B8) and from the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex dat”, funded by the German Excellence Strategy. They also thank for funding form the European Research Council in the ERC Synergy Grant “ECOGAL – Understanding our Galactic ecosystem: From the disk of the Milky Way to the formation sites of stars and planets” (project ID 855130). F.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 726384/Empire) P.S.B. acknowledges support from grant PID2019-107427GB-C31 from the Spanish Ministry of Science and Innovation A.U. acknowledges support from the Spanish grants PGC2018-094671-B-I00, funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, and PID2019-108765GB-I00, funded by MCIN/AEI/10.13039/501100011033. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: Toy model to test the hierarchical fitting

In this section, we use a toy model to test the performance of the hierarchical fitting algorithm and also to explore the origin of the covariance between the parameters C_⋆ and C_mol. To this end, we define 20 data subsets, each one having 3000 randomly selected data points from our full data set (Σ_⋆ and Σ_mol). For each subset, a pair of C_⋆ and C_mol is randomly generated from normal distributions centered at 0.5 and 1, respectively, with standard deviation of 1. C_norm is fixed to the value measured in our data at 150 pc for all subsets, to replicate the procedure we used for our data. The three coefficients are then used to generate simulated Σ_SFR values for each subset separately. A scatter of ∼0.6 dex (comparable to that found in our data) is added to the simulated Σ_SFR values.

Once the data subsets are created, we use the hierarchical fitting routine to find the best set of C_⋆ and C_mol to describe each subset, whereas C_norm is fitted for all subsets simultaneously. In this test, the hyperprior distributions used are C_⋆, μ ∼ 𝒩(0.5, 1²), C_⋆, σ ∼ ℋ(1²), C_mol, μ ∼ 𝒩(1, 1²), C_mol, σ ∼ ℋ(1²), C_norm, μ ∼ 𝒩(0, 2²) and C_norm, σ ∼ ℋ(2²).

Figure A.1 shows the plane conformed by the posterior distributions obtained for C_⋆ and C_mol for each data subset. The color code represents the posterior of C_norm. The black diamond marks the median of the posterior, whereas the blue square shows the real pair of C_⋆ − C_mol used to generate each simulated data set. This test primarily demonstrates the accuracy of the hierarchical fitting, as in all cases the real values of the coefficients fall somewhere within the posterior distribution, typically within 1σ of the recovered one. Secondly, this test shows that the covariance between C_⋆ and C_mol discussed in Sect. 4 is partially an artifact resulting from the fitting, as we recover it even though the real coefficients were chosen independently. However, we can not rule out that the covariance arises as a result of an intrinsic correlation between Σ_⋆ and Σ_mol.

Fig. A.1. Results of the toy model fitting to test accuracy of the hierarchical fitting as well as origin of the C⋆ − Cmol covariance. Each panel shows the posterior distribution obtained for a simulated data set. The black diamond and its error bar marks the median and 1-sigma of the posterior distribution. The blue square marks the real pair of C⋆ − Cmol used to generate each simulated data set. The color code scales with the posterior value of Cnorm.

Nevertheless, this brings up the question of how many data points are required to break this covariance? And thus, how many galaxies would we require to break (or at least reduce) the covariance in the probed galactic environments?

Figure 3 provides some clues to these questions. It shows a larger covariance in the environments with fewer data points (centers) than in those environments probed with more pixels (disks and spiral arms), while the lowest covariance is seen in the full sample. Then, if we would aim to reduce the covariance of the measurement in rings or centers (∼500 pixels) to the levels of spiral arms (∼15 000 pixels), we would need a sample on the order of 30 times larger. However, if we restrict the observations exclusively to galaxies that show all different environments (such as NGC 1672 or NGC 4321), the number of additional galaxies would need to be about ∼12 times larger. However, this is only an estimation based on the area identified as ring or center in these two galaxies. The exact number will depend on the area occupied by these features in each galaxy.

Appendix B: Choice of detection fraction threshold

Through this paper we adopt a detection fraction threshold of 60% to our data, which means that we confine our analysis to those Σ_⋆ ranges in which we have a detection fraction of Σ_SFR and Σ_mol higher than this threshold (see Sect. 3.1 for a detailed description). Here we show how our results are influenced by the choice of detection fraction threshold.

Figure B.1 shows the posterior distributions obtained for the coefficients C_⋆ and C_mol under different assumptions of detection fraction threshold. The top panel shows the distributions obtained without applying any detection fraction threshold (i.e., using the full sample). The bottom panels show our fiducial adopted value (60%), 50%, and 70%.

Fig. B.1. Posterior distributions obtained for the coefficients C⋆ and Cmol, following the same methodology described in Sect. 3.1, but applying different levels of detection fraction threshold. The top two panels show the base case scenario, using the full sample without applying any detection fraction threshold. The following panels show our fiducial adopted value (60%) and thresholds at 50% and 70%. The threshold level and the number of pixels used for each measurement are indicated in the top left corner for each row.

Firstly, we note that the posteriors of nondisk environments only change slightly from using the full sample to applying a threshold of 60% in the detection fraction (first three rows). This is because at the 60% level, the fraction of pixels dropped from these environments is very small (< 5%), as explained in Sect. 3.1. On the other hand, disks show a clear steepening of the C_mol slope in the same range of detection fraction. This behavior is consistent with the expected role of N/Ds in the measurement on the slope (Fig. 1).

We also note that only 5000 data points are lost when going from 50% to 60%, while going from 60% to 70% leads to drop ∼22 000 additional data points. Furthermore, at the 70% threshold level, the fraction of data points removed from spiral arms and rings rises to ∼11% and ∼35%, respectively. As a result of this higher number of dropped data points, and the significant truncation of the Σ_⋆ range in which these environments are probed, their posteriors change drastically at this detection fraction level, becoming unstable with respect to less strict threshold levels. Thus, we conclude here that a 60% threshold (fiducial value) provides a good balance between minimizing the impact of N/Ds on our analysis, and keeping the integrity of our data.

Appendix C: Global versus per-environment detection fraction threshold

In Sect. 3.1 we present our methodology to minimize the impact of N/Ds in the slope measurement, by restricting our analysis to the Σ_⋆ in which each environment has a detection fraction higher than 60%. Here, we present our results under a slightly different approach. Instead of defining the used Σ_⋆ range for each environment separately, we use the same Σ_⋆ range for all environments, defined on the detection fraction of the full sample (black line in Fig. 1).

As our full sample is statistically dominated by the disk environment, this implies that the Σ_⋆ range at which the rest of the environments are probed is truncated to roughly match the Σ_⋆ range in which disks satisfy the detection fraction threshold. In practice, this has an impact mainly in the sampling of spiral arms and rings, which due to their inherently higher detection fractions (Fig. 1), can be probed in a larger range of Σ_⋆ values following our fiducial approach.

Figure C.1 shows the obtained posterior distributions of the C_⋆ and C_mol coefficients in each environment, for four different levels of detection fraction threshold (similar to Fig. B.1). Due to the truncation of the Σ_⋆ range, the posterior of the spiral arms becomes unstable, and highly dependant on the threshold level applied, reaching almost a value of C_mol ≈ 1 for a detection fraction threshold of 70%. Thus, we opt for a “dynamic” Σ_⋆ range, chosen for each environment separately, as our fiducial approach.

Fig. C.1. Posterior distributions obtained for the coefficients C⋆ and Cmol, restricting the Σ⋆ range of different galactic environments according to the detection fraction of the full sample, for different levels of detection fraction threshold. The top two panels shows the base case scenario, using the full sample without applying any detection fraction threshold. The following panels show our fiducial adopted value (60%) and thresholds at 50% and 70%. The threshold level and the number of pixels used for each measurement are indicated in the top left corner for each row.

We note that the total number of data points in the sample at this threshold is drastically increased with respect to the fiducial scenario. This is because additional data points from the disk environment are included in the sample (as the detection fraction of the full sample is slightly higher than that for the disk only), at the cost of excluding data points from other environments.

Nevertheless, we find that environmental differences, in terms of C_⋆ and C_mol, persist under this different methodology, and thus, qualitatively our conclusions are robust against this different treatment of our data.

Appendix D: Correlation of ⟨σ_Hα⟩ and ⟨f_mol⟩ with coefficients calculated in individual environments of single galaxies

Here we show the correlation between the C_⋆ and C_mol coefficients derived for individual environments in single galaxies (see Sect. 5.2) and the second-highest Pearson coefficient parameters in Fig. 8. Figure D.1 shows the C_mol calculated for each environment, as a function of the average σ_Hα within that environment, for each galaxy in the subsample defined in Sect. 5.2. This correlation could be an imprint of how the dynamical state of the gas might play role in regulating the level of local SFR.

Fig. D.1. Cmol − ⟨σHα⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.

For comparison, Fig. D.2 shows the correlation between C_⋆ and ⟨f_mol⟩ for each environment, for the same galaxies. A positive correlation implies that environments with higher gas fractions also show typically lower depletion times.

Fig. D.2. C⋆ − ⟨fmol⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.

All Tables

All Figures

Fig. 1. rSFMS (bottom left) and rMGMS (bottom right) measured at 150 pc resolution for all galaxies in our sample. For reference and to guide the eye, the various solid colored lines represent the binned trends (bins of 0.15 dex) obtained for each environment, accounting for nondetections as detailed in Pessa et al. (2021). The black solid line shows the same measurement for all environments simultaneously. The dashed lines show the binned trends obtained when nondetections are neglected from the analysis. The panels on the top row show the detection fraction of ΣSFR (left) and Σmol (right), defined as the fraction of pixels with a measurement above our detection threshold in each bin of Σ⋆. The horizontal dotted line marks the detection fraction level of 60%.

In the text

	Fig. 2. Posterior distributions for the coefficients C⋆, Cmol and Cnorm that define the star-forming plane in each galactic environment, measured at a spatial resolution of 150 pc. As the posterior distribution of Cnorm is considerably narrower than that of C⋆ and Cmol, the x-axis has been binned in smaller bins and the distribution renormalized for an easier visualization.
In the text

	Fig. 3. Posteriors distributions of the coefficients C⋆ and Cmol for each separate environment, measured at an spatial resolution of 150 pc, to show the covariance between them. The posteriors have been smoothed using a Gaussian kernel. The color scale for each environment indicates the 1-, 2- and 3-sigma confidence intervals.
In the text

Fig. 4. Partial residual for each independent variable, dropping the dependency from the second independent variable, found for each individual environment and for all environments together as a function of Σ⋆ (left) and Σmol (right), measured at an spatial resolution of 150 pc. The binned data are shown for each corresponding environment. The trend for the full sample is overplotted as a black dotted line.

In the text

	Fig. 5. Change of the coefficients that define the star-forming planes in each galactic environment with spatial scale. The color code for each environment is the same as used in Fig. 1. The cyan diamond and a magenta hexagon indicate the results reported by Sánchez et al. (2021) (EDGE–CALIFA sample) and Lin et al. (2019), respectively.
In the text

	Fig. 6. Posterior distribution of the coefficients that define the star-forming planes in each galactic environment, at each one of the spatial scales probed, from 150 pc to 1 kpc. The spatial scale and the number of pixels used for each measurement are indicated in the top left corner for each row.
In the text

Fig. 7. Comparison between plane and single power law models. Top: scatter with respect to the star-forming 2D plane (solid) and rKS (dashed) determined for each galactic environment, at the different spatial scales probed. Bottom: difference of the LOO score measured for the plane and the rKS models, where a > 0 indicates a preference for the plane model. The error bars are the 3σ standard errors. The black horizontal line at 0 marks where there would be no preference between the models. The plane model is preferred over the rKS at all spatial scales, increasingly so at finer spatial scales.

In the text

	Fig. 8. Overall Pearson coefficient to quantify the level of correlation between the C⋆ and Cmol parameters that define the star-forming plane in each environment, and the set of different parameters explored.
In the text

	Fig. 9. Cmol − ⟨τ/tff⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.
In the text

	Fig. 10. C⋆ − ⟨log SFE⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.
In the text

Fig. 11. Posterior distributions for the coefficients C⋆, Cmol and Cnorm that define the star-forming plane in each separate environment, using a constant αCO. The vertical dashed lines show the centers of the posteriors distribution obtained under our fiducial choice of αCO. As the posterior distribution of Cnorm is considerably narrower than that of C⋆ and Cmol, the x-axis has been binned in smaller bins and renormalized for an easier visualization.

In the text

Fig. A.1. Results of the toy model fitting to test accuracy of the hierarchical fitting as well as origin of the C⋆ − Cmol covariance. Each panel shows the posterior distribution obtained for a simulated data set. The black diamond and its error bar marks the median and 1-sigma of the posterior distribution. The blue square marks the real pair of C⋆ − Cmol used to generate each simulated data set. The color code scales with the posterior value of Cnorm.

In the text

Fig. B.1. Posterior distributions obtained for the coefficients C⋆ and Cmol, following the same methodology described in Sect. 3.1, but applying different levels of detection fraction threshold. The top two panels show the base case scenario, using the full sample without applying any detection fraction threshold. The following panels show our fiducial adopted value (60%) and thresholds at 50% and 70%. The threshold level and the number of pixels used for each measurement are indicated in the top left corner for each row.

In the text

Fig. C.1. Posterior distributions obtained for the coefficients C⋆ and Cmol, restricting the Σ⋆ range of different galactic environments according to the detection fraction of the full sample, for different levels of detection fraction threshold. The top two panels shows the base case scenario, using the full sample without applying any detection fraction threshold. The following panels show our fiducial adopted value (60%) and thresholds at 50% and 70%. The threshold level and the number of pixels used for each measurement are indicated in the top left corner for each row.

In the text

	Fig. D.1. Cmol − ⟨σHα⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.
In the text

	Fig. D.2. C⋆ − ⟨fmol⟩ correlation across different galactic environments, for the galaxies that satisfy our single-galaxy selection criteria. The color code for each environment is the same as used in Fig. 1.
In the text