© I. Pessa et al. 2021

1. Introduction

In the current paradigm of evolution of galaxies, star formation occurs inside cold and dense molecular gas clouds. This process is regulated by complex small- and large-scale physics, such as feedback from recently born stars and supernovae resulting from the death of the most massive stars, magnetic fields or hydrostatic pressure exerted by the baryonic mass (Kennicutt & Evans 2012; Ostriker & Shetty 2011; Heyer & Dame 2015; Kruijssen et al. 2019; Krumholz et al. 2019; Chevance et al. 2020a). As a result of the interplay of such regulating mechanisms, different scaling relations at galactic scales arise between the total amount of star formation in a galaxy and the physical quantities that contribute to its regulation.

In this regard, the star formation main sequence (SFMS) is a tight (scatter of ∼0.3 dex) relation between the total star formation rate (SFR) of a galaxy and its total stellar mass. It consists of a power law, with a slope of ∼1, and it has been studied in the local Universe and at a higher redshift (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). Similarly, the Kennicutt–Schmidt relation has been extensively studied as it offers an alternative perspective on what drives the SFR in a galaxy. It correlates the total SFR with the total amount of gas and is consistent with a power law of order unity, even though the methodology does have an impact on the specific quantitative description (Schmidt 1959; Kennicutt 1998; Wyder et al. 2009; Genzel et al. 2010, 2012; Tacconi et al. 2010; Bigiel et al. 2011; Schruba et al. 2011).

Recent studies have demonstrated that these relations hold down to kiloparsec and sub-kiloparsec spatial scales, although their scatter is expected to raise below a critical spatial scale due to statistical undersampling of the star formation process (Schruba et al. 2010; Feldmann et al. 2011; Kruijssen & Longmore 2014; Kruijssen et al. 2018). The so-called resolved star formation main sequence (rSFMS; Cano-Díaz et al. 2016; Abdurro’uf & Akiyama 2017; Hsieh et al. 2017; Liu et al. 2018; Medling et al. 2018; Lin et al. 2019; Ellison et al. 2021; Morselli et al. 2020) and resolved Kennicutt–Schmidt relation (rKS; Bigiel et al. 2008; Leroy et al. 2008, 2013; Onodera et al. 2010; Schruba et al. 2011; Ford et al. 2013; Kreckel et al. 2018; Williams et al. 2018; Dey et al. 2019) correlate the local star formation rate surface density (Σ_SFR) with the local stellar mass surface density (Σ_stellar) and molecular gas surface density (Σ_mol. gas), respectively. However, there remains debate about the slope of these scaling relations, and it is known that their values depend on the specific approach used for their calculation (e.g., fitting technique; Calzetti et al. 2012; de los Reyes & Kennicutt 2019). Previous works have uncovered that their slope and scatter may link the physics from global to smaller scales. Moreover, the scatter of the rKS has been interpreted as individual star-forming regions undergoing independent evolutionary life cycles (Schruba et al. 2010; Feldmann et al. 2011), and its dependence on spatial scale provides insight into the timescales of the star formation cycle (Kruijssen & Longmore 2014; Kruijssen et al. 2018). Additionally, Bacchini et al. (2019a, 2020) found a tight correlation between the gas and the SFR volume densities in nearby disk galaxies, suggesting that the scatter of the rKS could be in part related to projection effects caused by the flaring scale height of gas disks. In addition to the rSFMS and the rKS, recently Lin et al. (2019), Ellison et al. (2021), and Morselli et al. (2020) reported the existence of a resolved molecular gas main sequence (rMGMS) as the correlation between the local Σ_mol. gas and local Σ_stellar, and several other studies have explored different forms of correlations between these locally measured quantities (Matteucci et al. 1989; Shi et al. 2011, 2018; Dib et al. 2017; Dey et al. 2019; Lin et al. 2019; Barrera-Ballesteros et al. 2021a)

The existence of these scaling relations suggests that their global counterparts are an outcome of the local mechanisms that drive star formation, and studying these relations and their respective scatter provides powerful insight into the local physical processes regulating the formation of stars. While the correlation between Σ_SFR and Σ_mol. gas is physically more intuitive, since molecular gas is the fuel for new stars, the origin of the relation between Σ_SFR and Σ_stellar is more uncertain. It can be understood as the interplay between local hydrostatic pressure of the disk and feedback mechanisms regulating the formation of new stars (Ostriker et al. 2010; Ostriker & Shetty 2011). On the other hand, using the ALMA-MaNGA QUEnch and STar formation (ALMaQUEST) survey (Lin et al. 2020), Lin et al. (2019) conclude that the rSFMS likely arises due to the combination of the rKS and the rMGMS, while the latter results from stars and gas following the same underlying gravitational potential.

In this paper we aim to probe these scaling relations to understand the mechanisms that locally regulate star formation, using ∼100 pc spatial resolution data from the PHANGS survey. We measured the slope and scatter of these resolved scaling relations in individual galaxies and in our combined sample. We study the universality of these relations and the processes driving galaxy-to-galaxy variations. Additionally we study how the spatial scale of the data, fitting approach, and systematic assumptions, such as the CO-to-H₂ conversion factor, impact our results.

This paper is structured as follows: In Sect. 2 we present our data set, and in Sect. 3 we describe the methods used in our analysis. In Sect. 4 we present our main results. In Sect. 5 we discuss our findings and their implications. Finally our conclusions are summarized in Sect. 6.

2. Data

We use a sample of 18 star-forming galaxies, all of them are close to the SFMS of galaxies. These galaxies represent a subsample of the Physics at High Angular resolution in Nearby GalaxieS (PHANGS¹) survey (Leroy et al., in prep.). The galaxies from the PHANGS survey have been selected to have a distance lower than 20 Mpc to resolve the typical scale of star-forming regions (50 − 100 pc) and to be relatively face on (i < 60°) to limit the effect of extinction and make the identification of clouds easier. These galaxies have been chosen to be a representative set of galaxies where most of the star formation is occurring in the local Universe. Our sample is summarized in Table 1 where we use the global parameters as reported by Leroy et al. (in prep.) based on the distance compilation of Anand et al. (2021) as well as the inclinations determined by Lang et al. (2020).

2.1. VLT/MUSE

We make use of the PHANGS-MUSE survey (PI: E. Schinnerer; Emsellem et al., in prep.). This survey employs the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2014) optical integral field unit (IFU) mounted on the VLT UT4 to mosaic the star-forming disk of 19 galaxies from the PHANGS sample. These galaxies correspond to a subset of the 74 galaxies from the PHANGS-ALMA survey (PI: E. Schinnerer; Leroy et al., in prep.). However, we have excluded one galaxy from the PHANGS-MUSE sample (NGC 0628) because its MUSE mosaic was obtained using a different observing strategy, leading to differences in data quality. The target selection for the PHANGS-MUSE sample focused on galaxies from the PHANGS parent sample that had already available ALMA data, as part of the ALMA pilot project, or from the ALMA archival.

The mosaics consist of 3–15 individual MUSE pointings. Each pointing provides a 1′×1′ field of view sampled at per pixel, with a typical spectral resolution of ∼2.5 Å (∼40 km s⁻¹) covering the wavelength range of 4800 − 9300 Å. The total on-source exposure time per pointing for galaxies in the PHANGS-MUSE Large Program is 43 min. Nine out of the 18 galaxies were observed using wide-field adaptive optics (AO). These galaxies are marked with a black dot in the first column of Table 1. The spatial resolution ranges from to for the targets with and without AO, respectively. Observations were reduced using a pipeline built on ESOREX and developed by the PHANGS team² (Emsellem et al., in prep.). The total area surveyed by each mosaic ranges from 23 to 692 kpc². Once the data have been reduced, we have used the PHANGS data analysis pipeline (DAP, developed by Francesco Belfiore and Ismael Pessa) to derive various physical quantities. The DAP will be described in detail in Emsellem et al. (in prep.). It consists of a series of modules that perform single stellar population (SSP) fitting and emission line measurements to the full MUSE mosaic. Some of these outputs are described in Sects. 2.4 and 2.5.

2.2. ALMA

The 18 galaxies have CO(2-1) data from the PHANGS-ALMA Large Program (PI: E. Schinnerer; Leroy et al., in prep.). We used the ALMA 12 m and 7 m arrays combined with the total power antennas to map CO emission at a spatial resolution of . The CO data cubes have an rms noise of ∼0.1 K per 2.5 km s⁻¹ channel (corresponding to 1σ(Σ_mol. gas)≈2 M_⊙ pc⁻² per 5 km s⁻¹ interval). The inclusion of the Atacama Compact Array (ACA) 7 m and total power data means that these maps are sensitive to emission at all spatial scales. For our analysis we use the integrated intensity maps from the broad masking scheme, which are optimized for completeness and contains the entirety of the galaxy emission. The strategy for observing, data reduction and product generation are described in Leroy et al. (2021). As our fiducial α_CO conversion factor we adopt the local gas-phase metallicity in solar units (Z′≡Z/Z_⊙) scaled prescription as described in Accurso et al. (2017) and Sun et al. (2020a), that is α_CO = 4.35Z′^−1.6 M_⊙ pc⁻² (K km s⁻¹)⁻¹, adopting a ratio CO(2-1)-to-CO(1-0) = 0.65 (Leroy et al. 2013; den Brok et al. 2021; Saito et al., in prep.). The radially-varying metallicity is estimated from the radial profile of gas-phase abundances in H II regions, as explained in Kreckel et al. (2020). Azimuthal variations in the metallicity of the interstellar medium have been previously reported (Ho et al. 2017; Kreckel et al. 2020), 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 = 4.35 M_⊙ pc⁻² (K km s⁻¹)⁻¹, as the canonical value for our Galaxy (Bolatto et al. 2013) in Sect. 5.4.2. IC 5332 has no significant detection of CO(2-1) emission and therefore has been excluded from the rKS and rMGMS analysis.

2.3. Environmental masks

We have used the environmental masks described in Querejeta et al. (2021) to morphologically classify the different environments of each galaxy and label them as disk, 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 we follow Herrera-Endoqui et al. (2015). Finally, spiral arms are only defined when they are clearly dominant features across the galaxy disk (excluding flocculent spirals). First, 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. For S⁴G, we rely on the analytic log-spiral segments from Herrera-Endoqui et al. (2015), and performed new fits for the remaining galaxies. These environmental masks allow us to examine the variations of the SFMS, rKS, and rMGMS relations across different galactic environments.

2.4. Stellar mass surface density maps

The PHANGS-MUSE DAP (Emsellem et al., in prep) includes a stellar population fitting module, a technique where a linear combination of SSP templates of known ages, metallicities, and mass-to-light ratios is used to reproduce the observed spectrum. 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 (which will be used in Sect. 5.4.3). 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 used a Voronoi tessellation (Cappellari & Copin 2003) to bin our MUSE data to a minimum signal-to-noise ratio (S/N) of ∼35, computed at the wavelength range of 5300 − 5500 Å. We chose this value in order to keep the relative uncertainty in our mass measurements below 15%, even for pixels dominated by a younger stellar population. To do this, we tried different S/N levels to bin a fixed region in our sample, and we bootstrapped our data to have an estimate of the uncertainties at each S/N level.

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, ranging from 30 Myr to 13.5 Gyr, logarithmically-spaced, and six metallicity bins [Z/H] = [ −1.49, −0.96, −0.35, +0.06, +0.26, +0.4]. 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 SSP fitting was done in two steps. First, we fitted our data assuming a Calzetti et al. (2000) extinction law to correct for internal extinction. We then corrected the observed spectrum using the measured extinction value before fitting it a second time, including a 12° multiplicative polynomial in this iteration in the fit. This two-step fitting process accounts for offsets between individual MUSE pointings. The different MUSE pointings are not necessarily observed under identical weather conditions, and therefore, even with careful treatment, we find some systematic differences between the individual MUSE pointings related to the different sky continuum levels. We studied regions in our mosaic where different pointings overlap, and found variations on the order of ∼3% (between an identical region in two different pointings). After inducing similar perturbations on a subset of spectra, we found that even these small differences could potentially cause systematic differences in measured stellar-population parameters. Therefore, in the first iteration of the SSP fitting, we measure a reddening value, and in the second iteration, we use a high-degree multiplicative polynomial to correct for those nonphysical features and homogenize the outcome of the different pointings. Additionally, we have identified foreground stars as velocity outliers in the SSP fitting and we have masked those pixels out of the analysis carried out in this paper.

2.5. Star formation rate measurements

As part of the PHANGS-MUSE DAP (Emsellem et al., in prep.), we fit 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 were able to construct emission lines flux maps for every galaxy. In order to calculate our final SFR rate measurement we use the Hα, Hβ, S II and O III emission line maps. We de-reddend the Hα fluxes, assuming that Hα_corr/Hβ_corr = 2.86, as appropriate for a case B recombination, temperature T = 10⁴ K, and density n_e = 100 cm², following:

(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 reddening in 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 determine whether the Hα emission comes from gas ionized by recently born stars or by a different source, such as active galactic nuclei (AGN), or low-ionization nuclear emission-line regions (LINER). We performed 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), to remove pixels that are dominated by AGN ionization from our sample. In the remaining pixels, we determined the fraction C_H II of the Hα emission actually tracing local star formation, and the fraction deemed to correspond to the diffuse ionized gas (DIG), 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., in prep.). To this end, we followed the approach described in Blanc et al. (2009) with the modifications introduced in Kaplan et al. (2016). Essentially, we first use the [S II]/Hα ratio to estimate C_H II in each pixel, following:

(2)

where and correspond to the typical [S II]/Hα ratio of DIG and H II regions as measured in the faint (10th percentile) and bright (90th percentile) end of the Hα distribution, respectively. Then, we perform a least squares fitting to find the best β and f₀ parameters such that:

(3)

Equation (3) represents the fraction of Hα emission tracing local star formation as a function of the Hα flux in each pixel. For Hα fluxes lower than f₀, the fraction is defined to be zero. This has been done for each galaxy separately. Bright pixels, dominated by star formation ionization have C_H II ∼ 1, while fainter and DIG-contaminated pixels have lower values. We include in Appendix A an example of the fitting of the parametrization defined in Eq. (3) to our data.

Finally, we compute the Hα emission tracing star formation (Hα_H II) as C_H II × Hα in each pixel, while the fraction 1 − C_H II is deemed to be DIG emission (Hα_DIG). We then calculate the total Hα_DIG to Hα_H II ratio (f_DIG) and rescale the Hα_H II flux of all pixels by (1 + f_DIG). This correction is performed because photons that ionize the DIG, originally escaped from H II regions. It represents, therefore, a spatial redistribution of the Hα flux. This approach permits us to estimate a star formation rate even in pixels contaminated by non-star-forming emission. A S/N cut of 4 for Hα and 2 for Hβ was then applied before computing the star formation rate surface density map using Eq. (4). However, most of the low S/N pixels have C_H II ≈ 0, and therefore, the 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 non-detections (see Sect. 3.2). In Sect. 3.2 we discuss the importance of non-detections (i.e., pixels with non-measured SFR) in our analysis. However, our main findings do not change qualitatively when we set to zero the SFR of these pixels (i.e., pixels that are not dominated by star-forming ionization according to the BPT criterion). Additionally, in Sect. 5.4.1 we discuss the impact of removing all Hα emission not associated with morphologically-defined H II regions.

To calculate the corresponding star formation rate from the Hα flux corrected for internal extinction and DIG contamination, we adopted the prescription described in Calzetti (2013):

(4)

This equation is scaled to a Kroupa universal IMF (Kroupa 2001), however, differences with the Chabrier IMF assumed for the SSP fitting are expected to be small (Kennicutt & Evans 2012). With these steps we obtain SFR surface density maps for each galaxy in our sample. We acknowledge that Eq. (4) assumes a fully-sampled IMF, and that the lowest SFR pixels (especially at high spatial resolution) may not form enough stars to fully sample the IMF. Hence, the measured SFR is more uncertain in this regime. However, due to our methodology to bin the data (see Sect. 3.2), the higher uncertainty in the low Σ_SFR regime has little effect in our analysis.

3. Methods

3.1. Sampling the data at larger spatial scales

In order to probe the relations under study at different spatial scales, we have resampled our native resolution MUSE and ALMA maps to pixel sizes of 100 pc, 500 pc and 1 kpc. The degrading of the data to larger spatial scales is done by only resampling into larger pixel, rather than performing a convolution before resampling. The 100 pc pixels are generally larger than the native spatial resolution of the maps. The resampling has been done for each one of the three relevant quantities: stellar mass surface density, star formation rate surface density, and molecular gas mass surface density. For stellar mass maps, the resampling was directly performed in the MUSE native resolution stellar mass surface density map, produced by the PHANGS-MUSE DAP. Calculating a resampled star formation rate surface density map required the resampling of each one of the line maps used for the BPT diagnostic and extinction correction. The parameters f₀ and β used for the DIG correction are calculated only at the native MUSE resolution, and do not change with spatial scale. This assumes that the typical DIG flux surface density does not vary with spatial scale, and the measurement at the native resolution is better constrained due to the higher number of pixels. Once the emission line maps were resampled, pixels dominated by AGN ionization were dropped from the analysis. The remaining Hα emission was then corrected by internal extinction and for DIG contamination as explained in Sect. 2.5. The DIG contamination correction is done independently at each spatial scale. For the molecular gas mass surface density, we have proceeded similarly to the stellar mass surface density. Additionally we have imposed a S/N cut of 1.5 for the molecular gas mass surface density map after the resampling, dropping the faintest and most uncertain pixels. For each one of the resampled quantities, we have also resampled the corresponding variance map to perform the S/N cut and report the corresponding uncertainty. Finally, we have corrected the maps by inclination, using a multiplicative factor of cos(i), where i corresponds to the inclination of each galaxy, listed in Table 1. The adopted inclinations correspond to those reported in Lang et al. (2020).

Figure 1 shows the star formation rate surface density map (top), molecular gas mas surface density map (center), and stellar mass surface density map (bottom) at each of the spatial scales probed in this work (100 pc, 500 pc, and 1 kpc from left to right) for one example galaxy (NGC 1512).

Fig. 1. Example of SFR, molecular gas, and stellar mass surface density (top, middle, and bottom row) at spatial scales of 100 pc, 500 pc and 1 kpc (left, middle, and right column) for one of the galaxies in our sample (NGC 1512). The black contour in each row encloses the pixels within the bar of the galaxy. Foreground stars have been masked in the stellar mass surface density maps (white pixels in bottom row).

3.2. Fitting technique

To fit each scaling relation, we binned the x-axis of our data in steps of 0.15 dex, and calculated the mean within each bin. Fitting bins rather than single data points avoids giving statistically larger weight to the outer part of a galaxy, where more pixels are available. A minimum of 5 data points with nonzero signal per bin was imposed at the 100 pc spatial scale for individual galaxies, in order to avoid sparsely sampled bins in the high or low end of the x-axis. For the full-sample relations, we imposed a minimum of 10 nonzero data points per bin at all spatial scales. We have also tested radial binning in the x-axis (i.e., calculating the means in bins defined as data points located at similar galactocentric radius), but this raised the scatter within each bin so we kept the x-axis binning, namely by stellar (gas) mass surface density.

The error in each bin has been calculated by bootstrapping: For a bin with N data points, we repeatedly chose 100 subsamples (allowing for repetitions), perturbed the data according to their uncertainties, and calculated their mean values. We adopt the standard deviation of the set of 100 mean values as the uncertainty of the mean in that bin. We opted for this approach instead of standard error propagation because the uncertainties of our SFR measurements are too small, and this method offers a more conservative quantification of the scatter within each bin. Similarly, for each binning resultant from the 100 iterations, we fit a power law and calculate a slope. To do this, we use a weighted least square fitting routine (WLS), where each bin is weighted by the inverse of its variance. This is justified since, by construction, the x-axis uncertainty of each bin is negligible, compared to its y-axis uncertainty. The final slope and its error correspond to the mean and standard deviation of the slopes distribution, respectively. This quantification of the error accounts for sample variance and statistical uncertainty of each data point, but it does not reflect the uncertainties induced by systematic effects. Finally, the scatter reported throughout the paper for the fitted power laws corresponds to the median absolute deviation of each data point, considering detections only, with respect to the best-fitting power law.

However, at physical resolutions of ∼100 pc, we deal with the issue that within each bin, we observe (in the case of the rSFMS) a bimodal distribution of SFR. At a given stellar mass surface density, a fraction of the pixels probed have a nonzero value of SFR that correlates with its corresponding stellar mass surface density value (with a certain level of scatter), while the remaining pixels do not show any SFR within our detection limits. This is either because these pixels are intrinsically non-star-forming, or because their SFRs are lower than our detection threshold. The fraction of these “non-detections” (N/D) is higher in the lower stellar mass surface density regime and close to zero at the high-mass end. This bimodality reflects the fact that star formation is not uniformly distributed across galactic disks – due to temporal stochasticity, or spatial organization of the star formation process due to galactic structure.

Thus, unlike studies done at ∼kpc resolutions, where star-forming regions smaller than the spatial resolution element will be averaged in a larger area, we need to properly account for the N/D fraction when investigating the scaling relations. This requires designing a fitting method such that our measurements at high resolution will be consistent to those obtained at larger spatial scales.

For the analysis presented in this paper, we are interested in (1) measuring how many stars are being formed per unit time on average at a given stellar mass surface density. This is different from asking (2) what is the typical SFR surface density at a given stellar mass surface density. The former requires us to include the non-detections as we are interested in averaging the SFR across the entire galactic disk, while the latter only tells us what are the most common SFR values to expect, and will depend on the probed spatial scale.

Measurements at lower spatial resolution are closer to addressing (1), as the physical interpretation of this is that in a given region of constant mass surface density, the mean SFR will be dominated by a few bright star-forming regions rather than by the much more numerous faint star-forming regions. Figures 2 and 3 exemplify this difference for the rSFMS derived using all the available pixels in our data at a spatial resolution of ∼100 pc. Figure 2 shows three different binning schemes for the rSFMS in the bottom panel, and the detection fraction (defined as 1 − f_N/D) of each bin in the top panel, where f_N/D quantifies the fraction of N/D in each bin.

Fig. 2. Bottom panel: 2D distribution of pixels in the overall rSFMS (i.e., including all pixels in our sample) and the three different binning schemes described in Sect. 3.2. As explained in the main text, different binning schemes address different questions. The red line shows the fiducial binning scheme adopted in this paper. Top panel: detection fraction of our SFR surface density tracer within each stellar mass surface density bin ranging from ∼20% in the low surface density regime to 100% at the high surface density end.

Fig. 3. Distribution of the SFR surface density values within the mass bin of log Σstellar [M⊙ pc−2]≈8.5. Around 40% of the pixels are not detected in ΣSFR and we observe a bimodality in the distribution. These non-detections have been highlighted in brown and moved to an artificial value (N/D). The vertical lines correspond to the derived average value for the three binning schemes in Fig. 2. The mean in log-space (green dashed line) matches the peak of the distribution of detections, while the mean in linear space (red and blue solid lines) correspond to the average SFR at this stellar mass surface density, including and excluding N/D, respectively.

When averaging in linear space either including non-detections as zeros (red line) or excluding them from the average calculation (blue line) we address question (1), whereas averaging in log-space (green line) provides information on (2) and hence goes through the bulk of the 2D distribution in the log Σ_SFR versus log Σ_stellar plane. Figure 3 shows the SFR distribution of pixels at an average log Σ_stellar [M_⊙ pc⁻²]≈8.5, and the average log Σ_SFR for this stellar mass surface density bin computed with each one of the three binning schemes is shown by the vertical lines (using the same color scheme as for Fig. 2). Non-detections in this bin are highlighted in brown. Here it becomes clear that while the mean in log-space (green line) matches the peak of the nonzero distribution, the means in linear space including or excluding non-detections (red and blue lines) are shifted toward higher values. Additionally the mean in linear space excluding non-detections (blue line) will always be greater or equal to the mean in linear space when non-detections are accounted for, as zeros in the calculation of the average (red line).

For our analysis we are interested in understanding if and how the star formation scaling relations (i.e., rSFMS, rKS and rMGMS) vary with measurement scale, and therefore we adopt the mean measured in linear space as our fiducial approach (i.e., red line in Fig. 2). To probe the effect of excluding non-detections on the slope determination, we use the blue binning scheme (i.e., mean in linear space excluding N/D) in order to perform a fair comparison to the fiducial case. However, our results (in terms of impact of spatial scale in the measured slope) are qualitatively unchanged if we use the mean in log space (i.e., green line) in this case instead.

Finally, we highlight here that for our fiducial approach, we are interested in probing the scaling relations across the entire galactic disk. Hence, we consider all the pixels from a given galaxy in the computation of the scaling relations, including faint and more uncertain SFR measurements, as well as N/D. We note that a pixel is defined as N/D (in a given scaling relation) if it has a nonzero measurement in the quantity shown on the x-axis, and a value below our detection limits in the quantity shown on the y-axis, i.e., a non-detection in molecular gas will lead to the omission of this pixel in the rKS, while it will be treated as an N/D in the rMGMS.

4. Results

In this section we present our measurements of the three relations under study. We recover all three scaling relations at a spatial scale of 100 pc and we explore how these relations change when the data are degraded to lower spatial resolutions. We present each relation when considering all available pixels from the full galaxy sample, and derived for each individual galaxy to study galaxy-to-galaxy variations. As mentioned in Sect. 2.2, IC 5332 has not been detected in CO(2-1) emission and therefore, the rKS and rMGMS have been measured for the remaining 17 galaxies only.

4.1. Scaling relations using the full sample

Figure 4 shows the 2D histograms of the overall rSFMS (left), rKS (center), and rMGMS (right), and their corresponding best-fitting power laws (red dashed line and red dots, respectively) at 100 pc spatial scale.

Fig. 4. 2D distribution of the overall resolved star formation main sequence (left), resolved Kennicutt–Schmidt relation (center), and molecular gas main sequence (right) at 100 pc spatial scale. The red points show the mean binned data and the red dashed line is the best-fitting power law. The colored region shows the 98% confidence interval of the linear fit.

The number of data points used is 313 227, 110 084 and 309 569 for the rSFMS, rKS, and rMGMS, respectively. This number is smaller for the rMGMS than for the rSFMS as IC 5332 is missing in the former. In the same order, the total detection fraction, defined as the fraction of pixels with a nonzero detection in the y-axis (1 − f_N/D), is 0.50, 0.84, and 0.35. This means that about 16 − 65% of the data points with a valid measurement of the property on the x-axis are not detected in the property on the y-axis. This reflects the high level of stochasticity dominating these relations at 100 pc spatial scale.

Out of these three relations, we find that the rMGMS has the lowest scatter (σ ≈ 0.34 dex), followed by the rKS and the rSFMS with scatter σ ≈ 0.41 dex and σ ≈ 0.51 dex, respectively (see bottom row in Table 4). However, one must be careful when interpreting these numbers, since the scatter is computed using the nonzero pixels only. Due to our methodology, we include faint nonzero SFR pixels (see Sect. 3.2) that enhance the scatter of the rSFMS and the rKS, especially at high resolution (see Sect. 5.1).

We have measured a slope of 1.05 ± 0.01 for the rSFMS, 1.06 ± 0.01 for the rKS, and 1.18 ± 0.01 for the rMGMS. In the same order, the intercepts that set the normalization of these relations are −10.63, −9.96, −2.23. These numbers are within the range of previously reported values, using lower resolution data. Differences in the slope possibly arise due to the use of different fitting methods (Hsieh et al. 2017) or differences in the treatment of non-detection (see discussion in Sect. 5.2). Similarly, variations in the normalization of these relations are expected due to differences in the assumed CO-to-H₂ conversion factor (see Sect. 5.4.2), as well as in the binning methodology (see Sect. 3.2). Figure 2 illustrates that the latter can lead to differences of up to ∼0.5 dex, with respect to the standard approach of average-in-log, depending on the fraction of N/D in a given bin. We summarize values from recent studies in Table 2. Some studies reported more than one value, using different fitting methods, often orthogonal distance regression (ODR) and ordinary least squares (OLS). In that case we show here the OLS value, as we obtained similar numbers for both fitting techniques in our data, when non-detections are excluded from the analysis. The best-fitting slopes and scatter for the overall sample are provided in the bottom row of Table 4.

4.2. Scaling relations in individual galaxies

Here we explore galaxy-to-galaxy variations for the three scaling relations. Figures 5–7 show the rSFMS, rKS, and rMGMS, respectively, for each individual galaxy in our sample. Different colors represent different galactic environments, namely disk, spiral arms, bar, rings (inner and outer), and centers (see Sect. 2.3). The black dashed lines show the corresponding overall best-fitting power law from Fig. 4 as reference and the magenta dashed line show the best-fitting power law for each individual galaxy. In Sect. 4.1 we find that the rMGMS is the relation with the lowest scatter in our overall sample. Figure 7 shows that the distribution of data points in this relation is much more compact (on a logarithmic scale) than in the rSFMS or the rKS. As mentioned earlier, this is due to the inclusion of very faint SFR pixels (see Sect. 3.2). As explained in Sect. 4.1, the rMGMS is the relation that has the lowest total detection fraction in comparison.

Fig. 5. rSFMS for each galaxy at 100 pc resolution. Different environment in the galaxies are colored according to the color scale next to the bottom right panel. The dashed black line represents the best-fitting power law to the overall measurement and the magenta line shows the best-fitting power law for each galaxy.

Fig. 6. rKS for each galaxy at 100 pc resolution. Different environments in the galaxies are colored according to the color scale next to the bottom right panel. The dashed black line represents the best-fitting power law to the overall measurement and the magenta line shows the best-fitting power law for each galaxy.

Fig. 7. rMGMS for each galaxy at 100 pc resolution. Different environments in the galaxies are colored according to the color scale next to the bottom right panel. The dashed black line represents the best-fitting power law to the overall measurement and the magenta line shows the best-fitting power law for each galaxy.

On the other hand, prominent differences between galaxies are seen for the rSFMS (Fig. 5). For instance, in some galaxies such as NGC 1365 or NGC 1566, their disk and spiral arm pixels agree well with the overall rSFMS, while in others, like NGC 3351 or NGC 4321, disk, spiral arms, and outer ring pixels tend to exhibit a constant SFR surface density independent of their local stellar mass surface density (albeit probing only a limited range here). Figure 8 highlights the diversity in the slopes measured for individual galaxies for the three relations probed. The rSFMS relation shows the largest dispersion in slopes among galaxies (σ ≈ 0.34), followed by the rMGMS (σ ≈ 0.32), and finally the rKS (σ ≈ 0.17), where galaxies show slopes generally closer to the overall value (black diamond). The two galaxies with larger uncertainties in the rKS best-fitting power law slopes correspond to NGC 2835 and NGC 5068, two low-mass and low-SFR galaxies. Figure 6 shows that these galaxies have only a small number of data points, hence their uncertain measurements.

Fig. 8. Slopes measured for the rSFMS, rKS and rMGMS for the individual galaxies in our sample. Each circle represents a galaxy. The size and color of each point scales with total SFR and stellar mass of the galaxy, respectively. The dispersion of the slope values for each relation is indicated on the top of the panel. The horizontal shift is arbitrary, with galaxies ordered from left to right by NGC number. The two galaxies with larger error bars in the KS law correspond to NGC 2835 and NGC 5068. This is due to the low number of data points available, as can be seen in Fig. 6. The black diamond with the magenta error bar indicates the global measurement for each relation. The horizontal gray dashed line show the median slope for each relation, in a galaxy basis.

Finally, we note that the shape of the rKS relation across different galaxies does not seem to change. While it has a larger scatter than the rMGMS at high resolution (0.41 dex versus 0.34 dex, respectively), the relation is indeed more uniform across different environments. On the contrary, certain galactic environments do impact the shape of the rSFMS and the rMGMS. When present, inner structures (bar, inner ring, and center) often “bend” these scaling relations. In the case of the rSFMS, this “bend” might be related to an increase of SFR in an inner ring or a suppression of SFR in the bulge or bar, which is not reflected by a change in the stellar mass surface density but it is reflected in the molecular gas surface density.

Figure 9 highlights variations across different morphological environments. The top row shows binned scaling relations for each environment separately, considering all the pixels in our sample, compared to the overall relation we derived for all environments together. We use the offset (Δ_env) defined as the vertical offset between the binned scaling relation from each environment and the overall best-fitting power law, and the slope measured for each environment to quantify deviations from the overall relation. It can be seen that the slope and normalization of the rKS are similar across different galactic environments while the variations for the rSFMS and rMGMS relations are substantially larger (up to ∼0.4 dex).

Fig. 9. Variations across different galactic environments for the rSFMS (left), rKS (center) and rMGMS (right) relations. Top row: binned data for each galactic environment separately, following the same color scheme from Fig. 7, while the black dashed line shows the overall best-fitting power law to all data points together. Bottom row: slope measured for each individual environment (on the x-axis). The vertical offset between each binned environment and the overall best-fitting power law is reported on the y-axis. It has been calculated by fitting a power law with a fixed slope (that of the overall relation) to each binned environment, and computing the intercept difference respect to the overall best fitting power law. Its error-bars represent the standard deviation of the difference between the binned environment and the overall best fit. The black diamond show the slope of the overall best-fitting power law (see Sect. 4.2 for discussion).

The disk environment (being the largest in area) is dominating the overall slope for rSFMS and rMGMS. The spiral arms share a similar slope, but are systematically offset above the overall relations by up to ∼0.4 dex. This is consistent with the findings reported in Sun et al. (2020a), where the authors found a systematic increase of Σ_mol. gas in spiral arms with respect to inter-arm regions analyzing 28 galaxies from the PHANGS-ALMA survey with identified spiral arms. Querejeta et al. (2021) also report higher average Σ_mol. gas and Σ_SFR in spiral arms compared to inter-arm regions, using data from the PHANGS-ALMA survey and narrow-band Hα imaging as SFR tracer. Outer rings (with low Σ_{stellar, mol. gas}) appear in the rSFMS and rMGMS as flatter features laying below the overall relation, while inner rings (with high Σ_{stellar, mol. gas}) show up above it. Finally, bars lie systematically below the overall relations, especially in the case of the rSFMS (see discussion in Sect. 5.1). To our knowledge this is the first time these three scaling relations are probed separately across different morphological environments of galaxies, at a spatial scale comparable to individual star-forming regions.

We investigated if any global galaxy parameter could be related to the diversity in the measured slopes. In particular, we explored any dependence of the slope on total stellar mass (M_⋆), total SFR, specific SFR (sSFR = SFR/M_⋆), and vertical offset of the galaxy from the global star formation main sequence of galaxies (ΔMS) as computed in Leroy et al. (2019). Only correlations for ΔMS are shown here, as the other parameters show even weaker or no correlations. Figure 10 shows the difference between the slope derived from the full sample and that of each individual galaxy (Δα_overall) as a function of ΔMS of each galaxy. Specifically, if α_gal corresponds to the slope measured for one individual galaxy, and α_overall is the slope measured for the overall sample (see Table 4), we define

(5)

Fig. 10. Differences in the slope measured for each individual galaxy with respect to the global measurements in Fig. 4 for the rSFMS. Each dot represent a galaxy in our sample. The gray dots are NGC 2835 and NGC 5068, two low-mass galaxies. The PCC of the correlation is indicated in the top-left corner.

Fig. 10 shows the slope variations with respect to ΔMS for the rSFMS (An analog figure for the rKS and the rMGMS is included in the Appendix B for completeness). NGC 2835 and NGC 5068 are marked as gray points. A weak correlation can be identified between the slope difference in the rSFMS and ΔMS (Pearson correlation coefficient (PCC) of ∼0.49). A potential origin of this correlation is discussed in Sect. 5.1. No correlation is found for the slope difference in the rKS or rMGMS with any of the tested global parameters. Table 4 summarizes the slopes measured for each galaxy and each scaling relation probed.

4.3. Variation of slope and scatter as a function of spatial scale

In this section, we explore how varying the spatial scale of the data affects the three relations. We resample our data to spatial scales of 500 pc, and 1 kpc as explained in Sect. 3.1. Figure 11 shows the three relations probed at spatial scales of 100 pc, 500 pc, and 1 kpc. The measured slopes and scatter are reported in Table 3. At all spatial scales, the slopes show no evidence of systematic dependence with spatial resolution. Uncertainties in the slope measurements are larger at 1 kpc due to the smaller number of data points, while the scatter is systematically lower at larger spatial scales, consistent with the findings reported in Bigiel et al. (2008), Schruba et al. (2010), Leroy et al. (2013), Kreckel et al. (2018). This results from averaging small scales variations of regions at different stages of their evolutionary cycle (e.g., Kruijssen & Longmore 2014). In fact, at 1 kpc resolution the rKS shows the lowest level of scatter.

Fig. 11. 2D histograms of the overall rSFMS (top row), rKS (center row), and rMGMS (bottom row) using all available pixels from our sample and probed at spatial scales of 100 pc (left column), 500 pc (middle column), and 1 kpc (right column). The x-axis binning and the best-fitting power law are indicated with red dots and a red dashed lines, respectively. The measured slope and its error are stated in each panel.

Finally, the number of pixels used at 1 kpc resolution drops to 2860, 1510, and 2820, and the fraction of pixels with detected signal increases to 0.62, 0.90, and 0.52 for the rSFMS, rKS, and rMGMS, respectively. This represents an increase of 24%, 7%, and 49% with respect to the 100 pc spatial scale.

5. Discussion

5.1. Main findings of our analysis

In the previous section, we presented our results on the scaling relations for both, our overall sample and individual galaxies. Our main findings are summarized in the following.

5.1.1. Variations with spatial scale between and within galaxies

We find in our overall sample that the rMGMS is the relation with the smallest scatter (at a spatial scale of 100 pc), followed by the rKS, and the rSFMS (see Table 3). However, the rKS shows the most consistency across different galaxies and environments (i.e., similar values of slope and normalization). Its larger scatter is due to the inclusion of low SFR pixels that deviate more from the overall relation. The consistency of the rKS across different galaxies and environments reflects the direct connection between molecular gas and SFR, since stars form out of molecular gas. Hence, SFR “follows” the molecular gas distribution within short timescales (Hα emission traces star formation in the last ∼10 Myr; Calzetti 2013; Leroy et al. 2012; Catalán-Torrecilla et al. 2015; Haydon et al. 2020) and this produces some correlated features in the shape of the rSFMS and the rMGMS for individual galaxies. NGC 1300, NGC 1512, and NGC 3351 are clear examples of these features. These similar features suggest that differences in the rSFMS across different environments are partially driven by the availability of molecular gas to fuel the formation of stars. However, bars, for example, generally appear below the overall best-fitting power law in the rKS as well, which suggests that variations in star formation efficiency may also play a role in driving the scatter. This consistency (across different galaxies and environments) suggest that physically, it is the more fundamental relation, compared to the rSFMS and the rMGMS, in agreement with what has been found by Lin et al. (2019) and Ellison et al. (2021).

In contrast, from the perspective of individual galaxies, we see that the galactic environment has a stronger impact on the rSFMS and the rMGMS than on the rKS. In particular, in the rSFMS, bars display a small range of SFR values, all systematically below the overall sequence. This apparent suppression of SFR in bars is consistent with the “star formation desert” (James & Percival 2018) in barred galaxies. Furthermore, variations in the rSFMS between different galaxies as those presented here have been previously reported (Hall et al. 2018; Liu et al. 2018; Vulcani et al. 2019; Ellison et al. 2021). Thus, a single power law across the full stellar mass surface density range does not provide a good representation of the Σ_stellar versus Σ_SFR plane; instead, different galactic environments may define different relations. However, here we aim at measuring these scaling relations to first order, i.e., a “simple” power law, and finding a more realistic description of the data goes beyond the scope of this paper.

We also investigated if any global galaxy parameter relates to the galaxy-to-galaxy variations in the slope of the scaling relations studied, and we found a potential correlation between the slope of the rSFMS of a galaxy and its ΔMS parameter. A positive correlation with a globally enhanced SFR could be explained as the increase of SFR happening preferably in the inner and more dense regions of the galaxy. This would lead to a steeper rSFMS in SFR-enhanced galaxies and it is consistent with the findings reported by Ellison et al. (2018), where a global enhancement or suppression of SFR was found to impact the inner region of a galaxy stronger than its outskirts, studying a sample of galaxies from the MaNGA survey (Bundy 2015). However, we stress here that even though the rSFMS shows hints of a correlation, the scatter and the limited number of data points make it hard to draw a robust conclusion.

Regarding spatial scale, we do not find evidence for a systematic dependence of the slope of these relations on the spatial resolution of the data, and the measurements are consistent within 2σ of their respective uncertainties. On the other hand, the scatter decreases at coarser spatial resolution as small scale variations are averaged out. The rKS is the relation that shows the least scatter at 1 kpc spatial scale.

5.1.2. Origin of observed scatter

Numerous authors have investigated the physical origins of the scatter in these resolved scaling relations and determined how this scatter evolves with spatial scale. According to Schruba et al. (2010) and Kruijssen & Longmore (2014), the scatter in the CO-to-Hα ratio at small spatial scales is dominated by the fact that 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). At larger spatial scales, these peaks are averaged and this sampling effect diminishes.

Along the same line, Semenov et al. (2017) presented a simple model to conciliate the long molecular gas depletion times measured at galactic scales (∼1 − 3 Gyr; Kennicutt 1998; Bigiel et al. 2008; Leroy et al. 2008, 2013) with the apparently shorter depletion times measured on ∼100 pc scales (∼40 − 500 Myr; Evans et al. 2009, 2014; Heiderman et al. 2010; Gutermuth et al. 2011; Schruba et al. 2017). In the proposed scenario, the difference in these time scales originates from the fact that not all the gas is going through the star formation process at the same time. It is suggested that the fraction of molecular gas, that is actively forming stars, is regulated by local (stellar feedback, turbulence, gravitational instabilities) and global (large-scale turbulence, differential rotation) mechanisms that continuously turn on and off the formation of stars. Consequently, several of these star formation cycles must occur in order to process all the molecular gas in a given volume. Thus, at high spatial resolution, a larger scatter in the distribution of depletion times (and thus in the rKS relation) results from the decoupling of the actively star-forming clouds from the quiescent ones.

The local mechanism is discussed more generally in Kruijssen & Longmore (2014) and Kruijssen et al. (2018), where the authors define a critical spatial scale on which the rKS breaks down from its integrated version, due to an incomplete sampling of the star-forming cycle. This spatial scale is defined as a function of the typical separation between independent star-forming regions, and the duration of the shortest phase of the star-forming process. This critical spatial scale equals the smallest region in which the star formation process is statistically well sampled. In this regard, having found the scatters at 100 pc σ_rSFMS > σ_rKS > σ_rMGMS is consistent with the expectation from this evolutionary scenario, provided that τ_Hα < τ_CO < τ_stars, where τ corresponds to the duration each tracer is visible across the star formation cycle. Due to the fact that the period of time in which young stars ionize their surrounding interstellar medium is shorter than the lifetime of molecular clouds, i.e., τ_Hα < τ_CO (Kruijssen et al. 2019; Chevance et al. 2020b), the impact of time evolution on the rMGMS will be smaller than on the rSFMS or the rKS. Consequently, the minimum spatial scale needed to properly sample the rMGMS is smaller than that needed for the rSFMS or the rKS. At spatial scales of ≳1 kpc, we are no longer in this “undersampling” regime, and the rKS shows the least scatter.

However, statistically incomplete sampling is not the only source of scatter for these resolved scaling relations. The slope of the rKS is also sensitive to the conditions present in the interstellar medium, such as metallicity and ultraviolet radiation field (Feldmann et al. 2011). Further Leroy et al. (2013) reported that variations in the measured depletion time (τ_dep = Σ_mol. gas/Σ_SFR) across different environments at ∼kpc spatial scales are consistent with variations of the CO-to-H₂ conversion factor. Similarly for the rSFMS, Ellison et al. (2020) found that at ∼kpc spatial scales, the scatter in this relation is likely driven by local variations in star formation efficiency. Here, we find large variations in these resolved scaling relation for different galaxies and galactic environments, in terms of both slope and normalization, which are particularly strong in the case of the rSFMS and rMGMS. These differences are key in setting their scatter, particularly at 1 kpc spatial scale, where we are no longer in a stochastic regime, for the reasons described earlier.

5.2. Role of non-detections in our analysis

In this section, we discuss the impact of excluding non-detections from the analysis. At high spatial resolution the non-detection fraction is quite large and their treatment can drastically alter the results.

5.2.1. Insights from simulations

In Calzetti et al. (2012), the authors studied variations in the slope of the observed rKS with different spatial resolution using simulated galaxies. They found that variations with spatial scale depend mainly on the slope of the true underlying rKS. In particular, for a slope of 1, they found a nearly constant slope in the observed relation (variation of ∼0.05) from 200 pc to 1 kpc resolution. However, due to the assumed underlying close H₂–SFR relation, both quantities share the same spatial distribution and hence have basically identical detection fractions which is not necessarily the case for our data.

Similarly, in Hani et al. (2020) the authors used a sample of Milky Way-like galaxies from FIRE-2 simulations to study the rSFMS measured at spatial scales of 100 pc, 500 pc, and 1 kpc. They reported a systematic steepening of this relation at larger spatial scales. It is important to note that only pixels with nonzero SFR values were considered in their analysis. They interpreted that this effect was primarily caused by differences in the detection fraction of the SFR across the galactic disk. A lower detection fraction in outer regions, where the stellar mass surface density is lower, would cause a stronger dilution of the SFR in the low stellar mass surface density regime than in the inner and denser regions when spatially averaging the data. As a consequence, the rSFMS has been measured to be flatter at 100 pc and to steepen at larger spatial scales, where the sparser intrinsic SFR distribution in the low mass surface density regime is leading to lower SFR surface density values.

This steepening of the rSFMS is mainly a consequence of the methodology used (i.e., excluding N/D in the fitting process). In the following, we investigate this effect in our data, and measure how much our slope measurements change under this alternative approach.

5.2.2. Slope changes in our scaling relations

Figures 12, C.1 and C.2, show the rSFMS, rKS, and rMGMS, respectively, when non-detections are excluded from the fit. We define α_D/O and σ_D/O as the measured slope and scatter considering pixels with nonzero detections only. Each relation studied here is presented at three spatial scales (100 pc, 500 pc, and 1 kpc from left to right). The slopes and scatter measured in the overall sample at each spatial resolution are summarized in Table 5, while Tables 6, C.1, and C.2 list the slopes and scatter measured for the rSFMS, rKS, and rMGMS, respectively, in each galaxy at each spatial scale. When we exclude non-detected pixels from our analysis, we see a systematic steepening of the slope at larger spatial scales in the rSFMS and the rMGMS. The steepening is particularly strong in the rMGMS due to lower (CO) detection fraction in this relation (∼35% at 100 pc). The rKS is less affected as molecular gas and SFR share a more similar spatial distribution (Schinnerer et al. 2019; Pan et al., in prep.) and non-detections in SFR frequently coincide with a non-detection in molecular gas, resulting in a high detection fraction of SFR for a given molecular gas value in the rKS relation (∼84% at 100 pc).

Fig. 12. 2D distribution of the resolved star formation main sequence considering all galaxies in our sample at three spatial resolutions (100 pc, 500 pc and 1 kpc from left to right). Non-detection pixels have been excluded from the fit. This produces a systematic flattening at higher spatial resolution. The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

Additionally, we observe a flattening at the low end of the stellar mass surface density (log Σ_stellar [M_⊙ kpc⁻²] ≲ 7.5) in the rSFMS and the rMGMS at 100 pc scale. This flattening is qualitatively similar to that reported by Barrera-Ballesteros et al. (2021a) or Cano-Díaz et al. (2019), attributed to Hα pollution from non-SF emission at low SFR values and due to the Hα detection threshold (see also Salim et al. 2007). As we correct for non-star-forming contributions in our Hα derived SFR maps, we interpret this flattening is caused by the larger fraction of non-detections in these low stellar mass surface density bins. Therefore, excluding non-detections leads to an overestimation of the mean measured in these bins, which results in the observed flattening at low stellar mass surface densities.

Even though we see flatter relations at 100 pc than at 1 kpc for the overall sample when non-detections are excluded, the behavior of individual galaxies is more diverse. Differences between the lowest and highest resolution measurements in each galaxy mainly depend on the spatial distributions of stellar mass and ionized gas (in the case of the rSFMS) across the galactic disk within each galaxy. However, for all galaxies the slopes obtained at 100 pc are flatter when non-detections are excluded, compared to when they are included.

5.2.3. Treatment of non-detections and its effect on the slope of the scaling relations

As explained in Sect. 5.2, the flattening of the relations toward small spatial scales arises when non-detections are excluded from the analysis. This is because “empty” regions at smaller spatial scales are averaged at larger spatial scales with regions including detections. This effectively dilutes the signal when the data are degraded, and the magnitude of this dilution will depend on the local detection fraction. The spatial distribution of Hα with respect to the stellar mass surface density is complex, and varies not only from galaxy to galaxy, but also between different environments within galaxies. Therefore, the local detection fraction follows an equally complex distribution.

Here we explore if we can recover a relation between the detection fraction and the amount of flattening, despite of the complex distribution of these quantities, using the slope measured at 100 pc spatial scale in a given galaxy, when N/D are excluded from the measurement. We measured the detection fraction of the SFR inside 1 kpc² boxes in the 100 pc spatial scale maps. For each box we calculated the corresponding detection fraction (DF_Hα), defined as the fraction of pixels with nonzero emission, and we computed the distribution of the DF_Hα within these boxes for each galaxy. Figure 13 shows a schematic representation of the approach for one example galaxy.

Fig. 13. Sketch to exemplify the methodology to compute the distribution of the Hα detection fraction. In each 1 kpc2 box we calculated the filling factor (DFHα) defined as the fraction of pixels with nonzero SFR.

Figure 14 shows the normalized distribution of DF_Hα for each galaxy. The corresponding mean (μ_DF) and standard deviation (σ_DF) are indicated in each panel. The latter contains information about the spatial configuration of the SFR. A galaxy with a compact (spatially concentrated) configuration will have a rather bimodal distribution of DF_Hα (and consequently higher σ_DFHα), where boxes will have mostly ∼0 or ∼1 values, whereas a more clumpy (assembled as individual separated clumps) configuration will lead to a flatter and uniform distribution.

Fig. 14. Normalized distribution of the Hα detection fraction (DFHα) inside 1 kpc2 size boxes in the galactic disk. Each panel corresponds to a galaxy from our sample. The mean (μ) and standard deviation (σ) are indicated in each panel.

We parametrize the difference in slope when non-detections are excluded from the fit as compared to the fiducial value as:

(6)

where α is the slope measured for a galaxy using the fiducial approach and α_D/O is the slope measured for the same galaxy when non-detections are excluded in the fit.

Figure 15 shows the change of slope when non-detections are excluded in the calculation as a function of μ_DFHα for the galaxies in our sample. The figure shows a negative correlation (PCC ≈ −0.68) between μ_DFHα and the change in the measured slope when non-detections are excluded. This can be interpreted as galaxies with higher average detection fraction are less affected by dilution of their signal, which implies that excluding non-detections has a smaller impact on the measured slope. So, despite of the complexity of the 2D distribution of ionized gas across the galactic disk, we identify a trend between the flattening and the distribution of the detection fraction.

Fig. 15. Difference between the slope of the rSFMS of each galaxy measured with the fiducial approach and the slope measured when non-detections are excluded, as a function of the mean detection fraction of the SFR tracer in each galaxy (μDFHα). All values are positive in the y-axis because excluding non-detections always leads to a flatter relation. The PCC of the correlation is indicated in the top-right corner of the panel.

The scatter in the correlation shown in Fig. 15 is likely because this diagnostic does not include any aspect of the underlying Σ_stellar (or Σ_mol. gas) distribution. Ultimately, the change of the slope will not only depend on the detection fraction in a given region, but also on its underlying Σ_stellar (for the rSFMS). Whereas the first is highly driven by local environment (Querejeta et al. 2021; Meidt et al. 2021), the second varies much more smoothly across galactic disks. The detection fraction distribution by its own cannot completely describe the flattening of the slope.

5.2.4. Role of detection fraction in the dispersion of slopes at 100 pc

In this section, we test if the dispersion that we see in the slopes measured for the rSFMS could be linked to differences in the SFR detection fraction distributions. Figure 16 shows Δα_overall (as defined in Eq. (5), the difference in slope of a given galaxy and the full sample) as a function of σ_DFHα. The figure shows a clear positive correlation (PCC ≈ 0.72), meaning that galaxies with larger spread (i.e., larger σ_DFHα) tend to have a steeper rSFMS (i.e., more positive Δα_overall) and galaxies with smaller spread tend to have a flatter rSFMS (i.e., more negative Δα_overall). This can be explained as galaxies with larger σ_DFHα (i.e., with a more bimodal distribution of DF_Hα) are associated with a spatial configuration in which the impact of non-detections is stronger in some regions (usually in the low stellar mass surface density regime) and less significant in others (central regions with higher stellar mass surface densities). This leads to a systematically steeper rSFMS by pulling down the low mass surface density end of the rSFMS, where non-detections dominate, with respect to the high mass surface density end. Galaxies with smaller σ_DFHα values show either a flatter DF_Hα distribution or a concentration of their values around zero. In these scenarios, non-detections do not contribute to steepen the relation in the same way, since different mass ranges are similarly affected.

Fig. 16. Difference between the slope of the rSFMS of each galaxy measured with the fiducial approach and the overall measurement as a function of the standard deviation of the detection fraction distribution of the SFR tracer in each galaxy (σDFHα). The PCC of the correlation is indicated in the top-right corner of the panel.

We find similar correlations between the detection fraction distributions of the molecular gas tracer and changes in the slope of the rMGMS (see Appendix D). Hence, the detection fraction of Hα (or molecular gas) is a relevant aspect to set the rSFMS (rMGMS) slope when measured at high resolution. This is also consistent with finding a similar dispersion in the slopes of the rSFMS and the rMGMS, and a significantly smaller dispersion in the slopes of the rKS. Due to the more similar spatial distribution between Hα and the molecular gas tracer, differences in the detection fraction are also less extreme.

5.3. Role of baryonic mass surface density in regulating SFR

Previous studies have demonstrated that the mid-plane pressure of the interstellar medium impacts its physical properties and consequently the local SFR surface density. Leroy et al. (2008) reported a good correlation between the mid-plane hydrostatic gas pressure and the local SFR surface density and gas depletion time in 23 nearby galaxies in the THINGS (Walter et al. 2008) and HERACLES (Leroy et al. 2009) surveys. Along the same line, Schruba et al. (2019) – studying 8 local galaxies – and Sun et al. (2020b) – studying 28 star-forming galaxies from the PHANGS-ALMA sample – found that the average internal pressure of molecular clouds tends to balance the sum of their own weights and the external interstellar medium pressure. Sun et al. (2020b) also reported a tight relationship between the mid-plane hydrostatic pressure and SFR surface density across their sample. These observations are in line with feedback-driven star formation models (Ostriker et al. 2010; Ostriker & Shetty 2011), in which the dynamical equilibrium of the interstellar medium in the galaxy gravitational potential regulates the local SFR surface density and vice versa.

Since both stars and molecular gas define the gravitational potential of a galaxy (neglecting the contribution from atomic gas), previous studies have also explored the existence of a relation between Σ_SFR and a combination of Σ_stellar and Σ_mol. gas (Matteucci et al. 1989; Shi et al. 2011, 2018; Dib et al. 2017; Dey et al. 2019; Lin et al. 2019; Barrera-Ballesteros et al. 2021a; Sánchez et al. 2021).

In Barrera-Ballesteros et al. (2021a), the authors computed the “baryonic” mass surface density as Σ_b = Σ_stellar + Σ_mol. gas at ∼kpc scales and found that Σ_b tightly correlates with SFR, with a slope of 0.97 and a residual scatter lower than that measured for the rSFMS or the rKS. We derived the Σ_b versus Σ_SFR relation at 100 pc spatial scale (see Fig. 17) and measured a slope of 1.21 ± 0.01 and a scatter of 0.49 dex. At 500 pc we obtain a slope of 1.22 ± 0.02 and a scatter of 0.46 dex, and at 1 kpc a slope of 1.14 ± 0.03 and a scatter of 0.42 dex. These values of scatter are higher than computed for the rKS and lower than those of the rSFMS at all spatial scales. Thus, we find that Σ_mol. gas is a better predictor of Σ_SFR than Σ_b. Differences with Barrera-Ballesteros et al. (2021a) are probably related to the treatment of N/D in the fit (which may directly impact the measured slope, as shown in Sect. 5.2), and the use of Balmer decrement as molecular gas surface density tracer (which is expected to introduce some level of scatter to the correlation).

Fig. 17. ΣSFR as a function of Σb, defined as Σmol. gas + Σstellar. We explored the existence of a tighter correlation with Σb as it is a better tracer for the total hydrostatic pressure exerted by the baryonic mass. We measured a slope of 1.18 ± 0.005 with a scatter, and a scatter lower than the rSFMS and higher than the rKS.

However, we stress here that the mid-plane hydrostatic pressure is not necessarily directly proportional to the baryonic mass surface density. Although Σ_b is sometimes used as a proxy for the hydrostatic pressure, it does not accurately reflect the gravitational potential felt by the gas disk in galaxies, because the stellar disk and gas disk usually have quite different scale heights (e.g., Ostriker et al. 2010; Bacchini et al. 2019a,b, 2020; Sun et al. 2020a). More careful investigations of the quantitative relation between Σ_mol, Σ_SFR, and the mid-plane hydrostatic pressure in the future will shed more light on the role of hydrostatic pressure in regulating star formation (e.g., Barrera-Ballesteros et al. 2021a,b).

5.4. Systematic effects

In this section, we test how some of our assumptions impact our measurements of the scaling relations probed. In Sect. 5.4.1, we derive the rSFMS and rKS relations after removing diffuse ionized gas (DIG) emission. In Sect. 5.4.2, we use a constant α_CO to measure the rKS and the rMGMS. Finally, in Sect. 5.4.3 we show the resultant rSFMS and rKS when a longer timescale SFR tracer is used instead of Hα.

5.4.1. Role of diffuse gas emission

In our fiducial approach, we aim at measuring scaling relations consistently at different spatial scales. Thus, we want to include any detectable SFR emission at high angular resolution (as well as N/D), that would then be averaged in the lower resolution measurement. However, even though we followed the approach described in Sect. 2.5 in order to correct the Hα emission by DIG contamination, we could still be left with some amounts of contaminant emission in our analysis. Here, we explore an alternative way to identify the star formation-associated Hα emission at the native MUSE resolution that is suitable to identify individual H II regions. For this, we use the H II region catalog from Santoro et al. (in prep.). In short, H II regions are identified in the Hα map using HIIphot (Thilker et al. 2000) for resolved sources and DAOStarFinder (Bradley et al. 2020) for point-like sources. A BPT diagnosis using the [O III]/Hβ and [N II]/Hα line ratios of the integrated spectra of each region was then used to select the Hα-emitting regions dominated by star formation, as described in Kewley et al. (2006). Once H II regions have been identified at the native MUSE resolution, we mask all remaining emission before computing our maps at 100 pc, 500 pc, and 1 kpc scales. The fraction of emission removed is on the order of ∼30% (∼60% of the pixels with nonzero emission in the fiducial approach are masked), consistent with the diffuse fractions obtained by decomposing the narrowband Hα images of these galaxies into diffuse and compact components in Fourier space (Hygate et al. 2019; Chevance et al. 2020b). This demonstrates that our procedure reliably removes all potentially contaminant emission before resampling the SFR map.

Figure 18 shows the resulting SFR map from this different approach at a spatial scale of 100 pc for NGC 1512. Individual H II regions are much more easily distinguished in the resultant map. Quantitatively, this creates large differences, mostly at the edges of H II regions or molecular clouds, where the signal is now averaged with neighboring zero-emission pixels during the resampling and pulled to lower levels of SFR. At the same time, it removes all faint, DIG-contaminated pixels from the relation. Figure 19 highlights this effect for the rSFMS of NGC 4254. The top panel shows the rSFMS from the fiducial approach, and the bottom panel shows the same relation using this alternative approach. The color scheme is the same as used in Fig. 5, indicating the morphological environment of each pixel. The black dashed line corresponds to the overall fit from the fiducial approach and has been overplotted as reference. Faint DIG-contaminated pixels are removed, which reduces the number of data points in the relation. However, the remaining pixels in the low SFR regime are pulled to even lower values resulting in a qualitatively similar rSFMS.

Fig. 18. Example of SFR surface density at spatial scales of 100 pc, for one of the galaxies in our sample (NGC 1512). All the Hα emission not associated with morphologically-defined H II regions has been excluded before re-sampling the maps. The methodology used to exclude this emission is described on the main text in Sect. 5.4.1.

Fig. 19. rSFMS of one galaxy in our sample (NGC 4254) to illustrate the effect of removing Hα emission not associated with morphologically defined H II regions before resampling the SFR surface density map. Top: fiducial methodology. Bottom: alternative methodology as described in Sect. 5.4.1. The color scheme is the same as in Fig. 6. The black line correspond to the global rSFMS measured using all pixels with the fiducial approach (see Fig. 4).

Figure 20 shows the 2D distributions for the rSFMS and the rKS probed at 100 pc, 500 pc, and 1 kpc scale with their corresponding binned trends and best-fitting power laws. The inferred slopes agree relatively well with those obtained using the fiducial method (see Fig. 11). Differences are expected since we are now fitting a power law to a subset of pixels (while the rest is classified as N/D). However, the scatter is systematically enhanced in both relations by a factor of ∼1.1 − 1.4, depending on the spatial scale considered (see Table 7). This increase is mainly driven by the effect previously described: the dilution of emission at the edges of H II regions or molecular clouds. Thus, we conclude that we could be underestimating the SF-associated emission by removing too much signal and, thus, introducing additional scatter to the relation. The true scatter of these scaling relations is probably somewhere in between our fiducial approach and applying the strict H II region mask. Hence, in our fiducial approach the scatter of these relations is possibly underestimated due to contaminant Hα emission boosting the SFR of the faintest pixels.

Fig. 20. Effect of excluding diffuse gas emission before resampling: 2D distributions of the overall rSFMS (top row) and rKS (bottom row) using all the available pixels in our sample, probed at spatial scales of 100 pc (left column), 500 pc (middle column) and 1 kpc (right column). The x-axis binning and the best-fitting power law are indicated with red dots and a red dashed line respectively. The measured slope and its error are indicated in each panel (see Sect. 5.4.1 for details). The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

5.4.2. Impact of the chosen CO-to-H₂ conversion factor

As explained in Sect. 2.2, our fiducial CO-to-H₂ conversion factor scales with local gas-phase metallicity. In this section, we explore the impact of assuming a constant conversion factor of 4.35 M_⊙ pc⁻² (K km s⁻¹)⁻¹ (Bolatto et al. 2013) on the measured rKS and rMGMS.

All of the galaxies in our sample exhibit either a negative or flat metallicity profile (Ho et al. 2015). This leads to a higher CO-to-H₂ conversion factor in the outer part of the galactic disk than in the central part. When the radial gradient is removed, differences between central (usually denser) clouds and outer (usually fainter) clouds are enhanced. This widens the range of molecular gas surface densities probed. Figure 21 shows the rKS (top) and rMGMS (bottom) at 100 pc spatial scale using a constant α_CO. The best-fitting power law obtained with the fiducial approach is indicated by the gray line. As a consequence of this widened range of molecular gas surface densities from applying a constant CO-to-H₂ conversion factor, the slope of the rKS decreases to 1.01 ± 0.01 and that of the rMGMS increases to 1.28 ± 0.01, representing a change of ∼5% and ∼8%, respectively, as compared to the fiducial scenario. Finally, the normalization of these relations is also affected by the CO-to-H₂ prescription. Under the assumption of a constant α_CO, we find an intercept of −9.29 for the rKS and −3.36 for the rMGMS (for reference, the intercepts computed under our fiducial α_CO prescriptions are −9.96 and −2.23, respectively).

Fig. 21. rKS (top) and rMGMS (bottom) measured using all pixels from our sample and assuming a constant CO-to-H2 conversion factor. This conversion factor leads to a flatter rKS and a steeper rMGMS (see Sect. 5.4.2 for a discussion). The best fit power law obtained with the fiducial approach is indicated with the gray line.

5.4.3. Probing SFR on a 150 Myr timescale with spectral fitting

So far, we have used Hα emission as our SFR tracer. Hα is known to trace the formation of stars during the last ∼10 Myr (Calzetti 2013; Leroy et al. 2012; Catalán-Torrecilla et al. 2015; Haydon et al. 2020). Here we explore how the slope measured for the rSFMS and rKS varies when we use a SFR tracer sensitive to longer star formation timescales. We use the resolved star formation histories (see Sect. 2.4) to map the fraction of total stellar mass assembled during the last 150 Myr (i.e., the youngest 4 age bins in our age–metallicity grid). Multiplying this fraction with the stellar mass surface density and dividing it by 1.5 × 10⁷ yr results in a SFR surface density map that probes longer timescales.

A longer timescale SFR tracer smooths the SFR values, i.e., high SFR values from the short timescale are nearly unaffected and low SFR values are pushed to higher values. Figure 22 shows how this change impacts the rSFMS (top) and the rKS (bottom). Consequently, both relations are significantly flattened. In particular, the slope of the rSFMS decreases to 0.61 ± 0.01 and the slope of the rKS to 0.60 ± 0.01.

Fig. 22. rSFMS (top) and rKS (bottom) measured at 100 pc resolution, using all available pixels from our sample and adopting a SFR tracer with a longer timescale. We have used the derived SFH to compute a SFR tracer sensitive to SF episodes in the last 150 Myr. This smooths the SFR measurements, flattens both relations, and reduces their scatter. The red dashed line show best-fitting power law to the binned data (red points). The best fit obtained with the fiducial approach is indicated with the gray line.

This also drastically reduces the scatter in both relations to 0.27 dex and 0.24 dex, respectively. A decrease in the scatter is consistent with the uncertainty principle reported in Kruijssen & Longmore (2014) and Kruijssen et al. (2018; also see discussion in Sect. 5.1). Averaging the SFR over a longer time scale increases the period of time in which young stars can be detected (i.e., τ_young-stars > τ_Hα), and thus reduces the critical spatial scale at which the relation breaks due to statistical undersampling of the star formation process.

We summarize that the timescales involved in the SFR determination strongly influence the resulting slope and scatter of the rSFMS and rKS relations. This agrees with a similar finding presented in Barrera-Ballesteros et al. (2021a).

5.5. Implication of our results

While the relation between Σ_SFR and Σ_mol. gas is direct, since star formation occurs in molecular clouds, the relation between Σ_SFR and Σ_stellar is understood as the interplay between the hydrostatic pressure exerted by the stellar (and cold gas) disk, together with feedback processes (such as stellar winds and supernovae) regulating the star formation. However, given that we found the rKS was tighter than the rSFMS and the Σ_SFR versus Σ_b correlation, we conclude that it is mainly the amount of available molecular gas that regulates the star formation rate rather than the amount of stellar or baryonic mass. This agrees well with what was reported by Lin et al. (2019), where the authors also conclude that the rSFMS could originate from the existence of the rKS and the rMGMS, where the latter can be explained either as the molecular gas following the gravitational potential of the stellar disk or both, stars and gas following the same underlying potential defined by the total mass. The similarities between the rSFMS and the rMGMS in our data across different galactic environments (see Figs. 5 and 7) and the high scatter of the rSFMS as compared to the other two scaling relations suggest that this might be the case, and that the substantial variations in the rSFMS are driven by a combination of abundance/lack of molecular gas to fuel star formation as well as variations in star formation efficiency.

The rMGMS is the relation with the least scatter among the three relations at a spatial scale of 100 pc, consistent with the expectation from the perpective of the time-scales of the star-forming cycle. The same was recently reported in Ellison et al. (2021) in an analysis carried out at kpc spatial scales. Interestingly, the rKS exhibits the most homogeneous behavior across different environments and between galaxies when measured at 100 pc scales. Variations in the slope of the rSFMS and the rMGMS across different galaxies seem to be related to differences in the detection fraction of either the SFR tracer or the molecular gas tracer. This effect is less important in the rKS, since molecular and ionized gas share more similar spatial distributions.

6. Summary

We have used VLT/MUSE and ALMA data from the PHANGS survey to derive SFR, molecular gas mass, and stellar mass surface densities across the galactic disks and to study the rSFMS (Σ_SFR versus Σ_stellar), rKS (Σ_SFR versus Σ_mol. gas), and rMGMS (Σ_mol. gas versus Σ_stellar) relations in a sample of 18 star-forming galaxies at spatial scales of 100 pc, 500 pc, and 1 kpc. We tested for systematic differences induced by spatial scales considered, fitting approaches used, and assumptions made. Additionally, we have explored the Σ_{stellar + mol. gas} − Σ_SFR correlation in our data to probe the effect of the mid-plane hydrostatic pressure of the disk as a regulator of the local SFR surface density. We applied a different approach to remove non-star-forming (diffuse) emission contaminating our SFR tracer before measuring the scaling relations, a different CO-to-H₂ conversion factor prescription, and a SFR tracer probing a longer timescale. Our main findings can be summarized as follows:

1. We have recovered all three scaling relations at a spatial scale of 100 pc. Of the three relations, rMGMS shows the least scatter (0.34 dex) for our global data set, whereas the rKS is the relation that shows the highest level of consistency between different galaxies and across environments. Its higher scatter in the high resolution (100 pc) measurement is related to the inclusion of very low surface density SFR data points, following our methodology to recover all potential SFR emission. When probed at 1 kpc scales, these data points are averaged over a larger region and the rKS shows the least scatter among the three relations (see Sect. 4.1).

2. At 100 pc, we found that the scatter of the scaling relations follows σ_rSFMS > σ_rKS > σ_rMGMS. This is consistent with the expectation from the evolutionary scenario perspective, given τ_Hα < τ_CO < τ_stars, where τ corresponds to the duration each tracer is visible across the star formation cycle (see Sects. 4.1 and 5.1.2).

3. We found significant variations in the studied scaling relations across different galactic environments. These variations are particularly strong in the case of the rSFMS and the rMGMS. The disk is the dominant feature in setting the slope, being the largest in area, while spiral arms share a similar slope, but are offset above the overall relations by up to 0.4 dex. Bars lie systematically below the overall relations (see Sect. 4.2).

4. We searched for global parameters that could be driving the dispersion in the measured slopes between the galaxies in our sample. We found a correlation with ΔMS, which implies that a global enhancement of SFR could change the slope of the scaling relations measured in a galaxy. However, we found a tighter correlation between the standard deviation of the distribution of the SFR tracer detection fraction with the deviation from the overall rSFMS slope. This can be explained by the fact that in galaxies with a larger spread in the detection fraction distribution, outer regions have typically a lower detection fraction values. Thus, the low stellar mass surface density regime will be more affected by non-detections than the inner region with a higher detection fraction, resulting in a steeper rSFMS (see Sects. 4.2 and 5.2.4).

5. As long as non-detections are included in the measurement of the slope, the spatial scale of the data do not greatly or systematically impact the measured slope. The scatter on the other hand decreases at larger spatial scales (see Sect. 4.3).

6. Excluding non-detections from the analysis artificially flattens the relation at smaller spatial scales, resulting in a steepening when the analysis is carried out at larger spatial scales. This is because pixels with nonzero signal are averaged with the non-detection pixels at larger spatial scales. Furthermore, this effectively causes an artificial flattening of the rSFMS and the rMGMS at the low mass surface density end in the 100 pc scale measurement (see Sect. 5.2.2).

7. How much the slope of a galaxy is flattened when non-detections are excluded at 100 pc spatial scale depends on the 2D distribution of ionized gas with respect to the stellar mass surface density (in the case of the rSFMS). We computed the distribution of the Hα detection fraction (DF_Hα) and we found a correlation between the mean of the detection fraction distribution and the change in slope when non-detections are not included as compared to the fiducial approach (see Sect. 5.2.3).

8. At all spatial scales, the scatter in the Σ_b versus Σ_SFR relation is higher than that seen for the rKS. We interpret this behavior such that the amount of available gas plays a primary role in locally regulating the SFR (see Sect. 5.3).

9. We removed all non-star-forming (diffuse) emission before resampling the SFR surface density maps, in order to reduce the level of contamination in our data. This strongly affects the level of emission at the edges of H II regions and molecular clouds, and increases the scatter of the rSFMS and the rKS. This suggests that with our fiducial approach, we may underestimate the real scatter of these relations due to diffuse non-SF-related flux boosting the emission of the faintest SF pixels (see Sect. 5.4.1).

10. We have recomputed the rKS and the rMGMS under the assumption of a constant CO-to-H₂ conversion factor. Removing the metallicity-dependent radial gradient from the conversion factor leads to a slightly flatter rKS (∼5% flatter) and a slightly steeper rMGMS (∼8% steeper; see Sect. 5.4.2).

11. We have recomputed the rSFMS and the rKS at 100 pc scale using a longer timescale SFR tracer derived from the MUSE star formation histories. This longer timescale tracer smoothens the SFR values with respect to the shorter timescale tracer. Both relations are significantly flattened and their scatter is drastically reduced (see Sect. 5.4.3).

Studying star formation in nearby galaxies at high physical resolution is a powerful tool to connect extragalactic observations with measurements in our own Galaxy. Assessing star formation scaling relations across different galaxy populations and quantifying systematic variations in different galactic environments will provide valuable insights in how galaxies grow and evolve in the local universe.

Acknowledgments

We thank the anonymous referee for their insightful comments which helped to improve the paper. This work was carried out as part of the PHANGS collaboration. DL, ES, HAP, TS, and TGW acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). FB acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 726384/Empire). The work of AKL and JS is partially supported by the National Science Foundation (NSF) under Grants No. 1615105, 1615109, and 1653300. E.C. acknowledges support from ANID project Basal AFB-170002. JMDK and MC gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through an Emmy Noether Research Group (grant number KR4801/1-1) and the DFG Sachbeihilfe (grant number KR4801/2-1). JMDK gratefully acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program via the ERC Starting Grant MUSTANG (grant agreement number 714907). KK gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI Kreckel). MQ acknowledges funding from the research project PID2019-106027GA-C44 from the Spanish Ministerio de Ciencia e Innovación. AU acknowledges support from the Spanish funding grants AYA2016-79006-P (MINECO/FEDER), PGC2018-094671-B-I00 (MCIU/AEI/FEDER), and PID2019-108765GB-I00 (MICINN). SCOG and RK acknowledge funding from the European Research Council via 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), from the 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 data”, funded by the German Excellence Strategy. 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.

References

Appendix A: Fitting of the C_H II parametrization

Figure A.1 shows the C_H II values (i.e., the fraction of Hα flux tracing star formation as defined in Eq. (2)) of the pixels in the NGC 1512 map, as a function of their Hα flux surface density. The solid black line shows the paramatrization defined in Eq. (3), and the best-fit parameters f₀ and β are listed in the top of the panel.

Fig. A.1. 2D distribution of the CH II fraction values of the pixels at the MUSE native resolution in NGC 1512 as defined in Eq. (2), as a function of the log Hα flux surface density of each pixel, in units of 10−20 erg s−1 cm−2 pc2. The black solid line shows the best-fitting parametrization, as described in Eq. (3). The obtained f0 and β parameters are shown in the top of the panel.

Appendix B: Δα_overall as a function of ΔMS for the rKS and the rMGMS

Figure B.1 shows the dependence of Δα_overall (as defined in Eq. (5)) from ΔMS, for the rKS and the rMGMS relations. We do not report any level of correlation between these parameters, and these figures are included here for completeness only.

Fig. B.1. Differences in the slope measured for each individual galaxy with respect to the global measurements in Fig. 4 for the rKS (left) and the rMGMS (right). Each dot represent a galaxy in our sample. The gray dots are NGC 2835 and NGC 5068, two low-mass galaxies. The PCC of the correlation is indicated in the top-left corner of each panel.

Appendix C: Measurement of the rKS and rMGMS at different spatial scales excluding non-detections

We include here the figures and tables with the measurement of the rKS and the rMGMS when non-detections are excluded from the analysis. Figures C.1 and C.2 show the obtained overall rKS and rMGMS respectively, at spatial scales from 100 pc to 1 kpc. Tables C.1 and C.2 show the corresponding slope and scatter obtained for the same relations in each individual galaxy at the different spatial scales.

Fig. C.1. 2D distribution of the resolved Kennicutt–Schmidt relation considering all galaxies in our sample at three spatial resolutions (100 pc, 500 pc and 1 kpc from left to right). Non-detection pixels have been excluded from the fit. The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

Fig. C.2. 2D distribution of the molecular gas main sequence considering all galaxies in our sample at three spatial resolutions (100 pc, 500 pc and 1 kpc from left to right). Non-detection pixels have been excluded from the fit. This produces a systematic flattening at higher spatial resolution. The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

Appendix D: Role of the detection fraction of the molecular gas tracer in the slope of the rMGMS

Figure D.1 shows how much the slope of the rMGMS in each galaxy varies when non-detections are excluded from its calculation, as a function of the mean in the distribution of the molecular gas tracer (μ_DFCO). We find a negative correlation (PCC ≈ −0.77), similar to that of the rSFMS. The outlier is the galaxy NGC 2835, which has a large uncertainty on its measured slope. Figure D.2 shows the difference between the slope of the rMGMS for each galaxy and the global measurement as a function of the standard deviation of the distribution of the molecular gas tracer (σ_DFCO). The level of correlation we find here is modest (PCC ≈ 0.50), in contrast to that of the rSFMS slope differences and the SFR tracer distribution. This could be due to the completeness limit of our molecular gas tracer. As stated in Sect. 4.1, the rMGMS is the relation with the lowest total detection fraction (∼35%). The detection insensitivity to a fainter component should not impact our slope measurements, as our methodology is robust against the detection threshold. However, the detection fraction distribution will be affected. This leads to systematically lower values of μ_DFCO for all galaxies. On the other hand, the impact in σ_DFCO is less systematic, and will depend on the mean of the distribution in each galaxy. We found a better correlation (PCC ≈ −0.70) of Δα_overall with μ_DFCO instead. Figure D.3 shows the difference in the rMGMS slope with the overall slope as a function of μ_DFCO, excluding the two galaxies with the more uncertain measurements (NGC 2835 and NGC 5068). The negative correlation implies that galaxies with a lower mean detection fraction of the molecular gas tracer have a steeper rMGMS.

Fig. D.1. Difference between the slope of the rMGMS of each galaxy measured with the fiducial approach and the slope measured when non-detections are excluded, as a function of the mean detection fraction of the molecular gas tracer in each galaxy (μDFCO). All values are positive in the y-axis because excluding non-detections always leads to a flatter relation. The PCC of the correlation is indicated in the top-right corner of the panel.

Fig. D.2. Difference between the slope of the rMGMS of each galaxy measured with the fiducial approach and the global measurement as a function of the standard deviation of the detection fraction distribution of the molecular gas tracer in each galaxy (σDFCO). The PCC of the correlation is indicated in the top-right corner of the panel.

Fig. D.3. Difference between the slope of the rMGMS of each galaxy measured with the fiducial approach and the global measurement as a function of the mean of the detection fraction distribution of the molecular gas tracer in each galaxy (μDFCO). The PCC of the correlation is indicated in the top-right corner of the panel.

All Tables

All Figures

Fig. 1. Example of SFR, molecular gas, and stellar mass surface density (top, middle, and bottom row) at spatial scales of 100 pc, 500 pc and 1 kpc (left, middle, and right column) for one of the galaxies in our sample (NGC 1512). The black contour in each row encloses the pixels within the bar of the galaxy. Foreground stars have been masked in the stellar mass surface density maps (white pixels in bottom row).

In the text

Fig. 2. Bottom panel: 2D distribution of pixels in the overall rSFMS (i.e., including all pixels in our sample) and the three different binning schemes described in Sect. 3.2. As explained in the main text, different binning schemes address different questions. The red line shows the fiducial binning scheme adopted in this paper. Top panel: detection fraction of our SFR surface density tracer within each stellar mass surface density bin ranging from ∼20% in the low surface density regime to 100% at the high surface density end.

In the text

Fig. 3. Distribution of the SFR surface density values within the mass bin of log Σstellar [M⊙ pc−2]≈8.5. Around 40% of the pixels are not detected in ΣSFR and we observe a bimodality in the distribution. These non-detections have been highlighted in brown and moved to an artificial value (N/D). The vertical lines correspond to the derived average value for the three binning schemes in Fig. 2. The mean in log-space (green dashed line) matches the peak of the distribution of detections, while the mean in linear space (red and blue solid lines) correspond to the average SFR at this stellar mass surface density, including and excluding N/D, respectively.

In the text

	Fig. 4. 2D distribution of the overall resolved star formation main sequence (left), resolved Kennicutt–Schmidt relation (center), and molecular gas main sequence (right) at 100 pc spatial scale. The red points show the mean binned data and the red dashed line is the best-fitting power law. The colored region shows the 98% confidence interval of the linear fit.
In the text

	Fig. 5. rSFMS for each galaxy at 100 pc resolution. Different environment in the galaxies are colored according to the color scale next to the bottom right panel. The dashed black line represents the best-fitting power law to the overall measurement and the magenta line shows the best-fitting power law for each galaxy.
In the text

	Fig. 6. rKS for each galaxy at 100 pc resolution. Different environments in the galaxies are colored according to the color scale next to the bottom right panel. The dashed black line represents the best-fitting power law to the overall measurement and the magenta line shows the best-fitting power law for each galaxy.
In the text

	Fig. 7. rMGMS for each galaxy at 100 pc resolution. Different environments in the galaxies are colored according to the color scale next to the bottom right panel. The dashed black line represents the best-fitting power law to the overall measurement and the magenta line shows the best-fitting power law for each galaxy.
In the text

Fig. 8. Slopes measured for the rSFMS, rKS and rMGMS for the individual galaxies in our sample. Each circle represents a galaxy. The size and color of each point scales with total SFR and stellar mass of the galaxy, respectively. The dispersion of the slope values for each relation is indicated on the top of the panel. The horizontal shift is arbitrary, with galaxies ordered from left to right by NGC number. The two galaxies with larger error bars in the KS law correspond to NGC 2835 and NGC 5068. This is due to the low number of data points available, as can be seen in Fig. 6. The black diamond with the magenta error bar indicates the global measurement for each relation. The horizontal gray dashed line show the median slope for each relation, in a galaxy basis.

In the text

Fig. 9. Variations across different galactic environments for the rSFMS (left), rKS (center) and rMGMS (right) relations. Top row: binned data for each galactic environment separately, following the same color scheme from Fig. 7, while the black dashed line shows the overall best-fitting power law to all data points together. Bottom row: slope measured for each individual environment (on the x-axis). The vertical offset between each binned environment and the overall best-fitting power law is reported on the y-axis. It has been calculated by fitting a power law with a fixed slope (that of the overall relation) to each binned environment, and computing the intercept difference respect to the overall best fitting power law. Its error-bars represent the standard deviation of the difference between the binned environment and the overall best fit. The black diamond show the slope of the overall best-fitting power law (see Sect. 4.2 for discussion).

In the text

	Fig. 10. Differences in the slope measured for each individual galaxy with respect to the global measurements in Fig. 4 for the rSFMS. Each dot represent a galaxy in our sample. The gray dots are NGC 2835 and NGC 5068, two low-mass galaxies. The PCC of the correlation is indicated in the top-left corner.
In the text

Fig. 11. 2D histograms of the overall rSFMS (top row), rKS (center row), and rMGMS (bottom row) using all available pixels from our sample and probed at spatial scales of 100 pc (left column), 500 pc (middle column), and 1 kpc (right column). The x-axis binning and the best-fitting power law are indicated with red dots and a red dashed lines, respectively. The measured slope and its error are stated in each panel.

In the text

Fig. 12. 2D distribution of the resolved star formation main sequence considering all galaxies in our sample at three spatial resolutions (100 pc, 500 pc and 1 kpc from left to right). Non-detection pixels have been excluded from the fit. This produces a systematic flattening at higher spatial resolution. The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

In the text

	Fig. 13. Sketch to exemplify the methodology to compute the distribution of the Hα detection fraction. In each 1 kpc2 box we calculated the filling factor (DFHα) defined as the fraction of pixels with nonzero SFR.
In the text

	Fig. 14. Normalized distribution of the Hα detection fraction (DFHα) inside 1 kpc2 size boxes in the galactic disk. Each panel corresponds to a galaxy from our sample. The mean (μ) and standard deviation (σ) are indicated in each panel.
In the text

Fig. 15. Difference between the slope of the rSFMS of each galaxy measured with the fiducial approach and the slope measured when non-detections are excluded, as a function of the mean detection fraction of the SFR tracer in each galaxy (μDFHα). All values are positive in the y-axis because excluding non-detections always leads to a flatter relation. The PCC of the correlation is indicated in the top-right corner of the panel.

In the text

	Fig. 16. Difference between the slope of the rSFMS of each galaxy measured with the fiducial approach and the overall measurement as a function of the standard deviation of the detection fraction distribution of the SFR tracer in each galaxy (σDFHα). The PCC of the correlation is indicated in the top-right corner of the panel.
In the text

	Fig. 17. ΣSFR as a function of Σb, defined as Σmol. gas + Σstellar. We explored the existence of a tighter correlation with Σb as it is a better tracer for the total hydrostatic pressure exerted by the baryonic mass. We measured a slope of 1.18 ± 0.005 with a scatter, and a scatter lower than the rSFMS and higher than the rKS.
In the text

	Fig. 18. Example of SFR surface density at spatial scales of 100 pc, for one of the galaxies in our sample (NGC 1512). All the Hα emission not associated with morphologically-defined H II regions has been excluded before re-sampling the maps. The methodology used to exclude this emission is described on the main text in Sect. 5.4.1.
In the text

Fig. 19. rSFMS of one galaxy in our sample (NGC 4254) to illustrate the effect of removing Hα emission not associated with morphologically defined H II regions before resampling the SFR surface density map. Top: fiducial methodology. Bottom: alternative methodology as described in Sect. 5.4.1. The color scheme is the same as in Fig. 6. The black line correspond to the global rSFMS measured using all pixels with the fiducial approach (see Fig. 4).

In the text

Fig. 20. Effect of excluding diffuse gas emission before resampling: 2D distributions of the overall rSFMS (top row) and rKS (bottom row) using all the available pixels in our sample, probed at spatial scales of 100 pc (left column), 500 pc (middle column) and 1 kpc (right column). The x-axis binning and the best-fitting power law are indicated with red dots and a red dashed line respectively. The measured slope and its error are indicated in each panel (see Sect. 5.4.1 for details). The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

In the text

	Fig. 21. rKS (top) and rMGMS (bottom) measured using all pixels from our sample and assuming a constant CO-to-H2 conversion factor. This conversion factor leads to a flatter rKS and a steeper rMGMS (see Sect. 5.4.2 for a discussion). The best fit power law obtained with the fiducial approach is indicated with the gray line.
In the text

Fig. 22. rSFMS (top) and rKS (bottom) measured at 100 pc resolution, using all available pixels from our sample and adopting a SFR tracer with a longer timescale. We have used the derived SFH to compute a SFR tracer sensitive to SF episodes in the last 150 Myr. This smooths the SFR measurements, flattens both relations, and reduces their scatter. The red dashed line show best-fitting power law to the binned data (red points). The best fit obtained with the fiducial approach is indicated with the gray line.

In the text

	Fig. A.1. 2D distribution of the CH II fraction values of the pixels at the MUSE native resolution in NGC 1512 as defined in Eq. (2), as a function of the log Hα flux surface density of each pixel, in units of 10−20 erg s−1 cm−2 pc2. The black solid line shows the best-fitting parametrization, as described in Eq. (3). The obtained f0 and β parameters are shown in the top of the panel.
In the text

	Fig. B.1. Differences in the slope measured for each individual galaxy with respect to the global measurements in Fig. 4 for the rKS (left) and the rMGMS (right). Each dot represent a galaxy in our sample. The gray dots are NGC 2835 and NGC 5068, two low-mass galaxies. The PCC of the correlation is indicated in the top-left corner of each panel.
In the text

	Fig. C.1. 2D distribution of the resolved Kennicutt–Schmidt relation considering all galaxies in our sample at three spatial resolutions (100 pc, 500 pc and 1 kpc from left to right). Non-detection pixels have been excluded from the fit. The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.
In the text

Fig. C.2. 2D distribution of the molecular gas main sequence considering all galaxies in our sample at three spatial resolutions (100 pc, 500 pc and 1 kpc from left to right). Non-detection pixels have been excluded from the fit. This produces a systematic flattening at higher spatial resolution. The fiducial overall best-fitting power law at each spatial scale is marked with the gray dashed line for reference.

In the text

Fig. D.1. Difference between the slope of the rMGMS of each galaxy measured with the fiducial approach and the slope measured when non-detections are excluded, as a function of the mean detection fraction of the molecular gas tracer in each galaxy (μDFCO). All values are positive in the y-axis because excluding non-detections always leads to a flatter relation. The PCC of the correlation is indicated in the top-right corner of the panel.

In the text

	Fig. D.2. Difference between the slope of the rMGMS of each galaxy measured with the fiducial approach and the global measurement as a function of the standard deviation of the detection fraction distribution of the molecular gas tracer in each galaxy (σDFCO). The PCC of the correlation is indicated in the top-right corner of the panel.
In the text

	Fig. D.3. Difference between the slope of the rMGMS of each galaxy measured with the fiducial approach and the global measurement as a function of the mean of the detection fraction distribution of the molecular gas tracer in each galaxy (μDFCO). The PCC of the correlation is indicated in the top-right corner of the panel.
In the text