Star formation scaling relations at ~100 pc from PHANGS: Impact of completeness and spatial scale

Aims: The complexity of star formation at the physical scale of molecular clouds is not yet fully understood. We investigate the mechanisms regulating the formation of stars in different environments within nearby star-forming galaxies from the PHANGS sample. Methods: Integral field spectroscopic data and radio-interferometric observations of 18 galaxies were combined to explore the existence of the resolved star formation main sequence (rSFMS), resolved Kennicutt-Schmidt relation (rKS), and resolved molecular gas main sequence (rMGMS), and we derived their slope and scatter at spatial resolutions from 100 pc to 1 kpc (under various assumptions). Results: All three relations were recovered at the highest spatial resolution (100 pc). Furthermore, significant variations in these scaling relations were observed across different galactic environments. The exclusion of non-detections has a systematic impact on the inferred slope as a function of the spatial scale. Finally, the scatter of the $\Sigma_\mathrm{mol. gas + stellar}$ versus $\Sigma_\mathrm{SFR}$ correlation is smaller than that of the rSFMS, but higher than that found for the rKS. Conclusions: The rMGMS has the tightest relation at a spatial scale of 100 pc (scatter of 0.34 dex), followed by the rKS (0.41 dex) and then the rSFMS (0.51 dex). This is consistent with expectations from the timescales involved in the evolutionary cycle of molecular clouds. Surprisingly, the rKS shows the least variation across galaxies and environments, suggesting a tight link between molecular gas and subsequent star formation. The scatter of the three relations decreases at lower spatial resolutions, with the rKS being the tightest (0.27 dex) at a spatial scale of 1 kpc. Variation in the slope of the rSFMS among galaxies is partially due to different detection fractions of $\Sigma_\mathrm{SFR}$ with respect to $\Sigma_\mathrm{stellar}$.


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. 2020b). 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;Tacconi et al. 2010;Bigiel et al. 2011;Schruba et al. 2011;Genzel et al. 2012).
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. 2020a;Morselli et al. 2020) and resolved Kennicutt-Schmidt relation (rKS;Bigiel et al. 2008;Leroy et al. 2008;Onodera et al. 2010;Schruba et al. 2011;Ford et al. 2013;Leroy 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 Article number, page 1 of 27 arXiv:2104.09536v3 [astro-ph.GA] 23 Jun 2021 A&A proofs: manuscript no. Paper_rSFMS 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. (2019aBacchini et al. ( , 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. (2020a), 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;Dib et al. 2017;Shi et al. 2018;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 galaxyto-galaxy variations. Additionally we study how the spatial scale of the data, fitting approach, and systematic assumptions, such as the CO-to-H 2 conversion factor, impact our results. This paper is structured as follows: In Sec. 2 we present our data set, and in Sec. 3 we describe the methods used in our analysis. In Sec. 4 we present our main results. In Sec. 5 we discuss our findings and their implications. Finally our conclusions are summarized in Sec 6.

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 1 ) survey (Leroy et al. 2021a). 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 1 http://phangs.org/ 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. (2021a) based on the distance compilation of Anand et al. (2021) as well as the inclinations determined by Lang et al. (2020).

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. 2021a). 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 to 15 individual MUSE pointings. Each pointing provides a 1 × 1 field of view sampled at 0 . 2 per pixel, with a typical spectral resolution of ∼2.5 Å (∼40 km s −1 ) covering the wavelength range of 4800−9300 Å. The total onsource 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 ∼0 . 5 to ∼1 . 0 for the targets with and without AO, respectively. Observations were reduced using a pipeline built on esorex and developed by the PHANGS team 2 (Emsellem et al. in prep.). The total area surveyed by each mosaic ranges from 23 to 692 kpc 2 . 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 Secs. 2.4 and 2.5.

ALMA
The 18 galaxies have CO(2-1) data from the PHANGS-ALMA Large Program (PI: E. Schinnerer; Leroy et al. 2021a). We used the ALMA 12m and 7m arrays combined with the total power antennas to map CO emission at a spatial resolution of 1 . 0−1 . 5. The CO data cubes have an rms noise of ∼0.1 K per 2.5 km s −1 channel (corresponding to 1 σ(Σ mol.gas ) ≈ 2 M pc −2 per 5 km s −1 interval). The inclusion of the Atacama Compact Array (ACA) 7m 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. (2021b). 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. (2020b), that is α CO = 4.35Z −1.6 M pc −2 (K km s −1 ) −1 , adopting a ratio CO(2-1)-to-CO(1-0) = 0.65 den Brok et al. 2021, T. 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 −2 (K km s −1 ) −1 , as the canonical value for our Galaxy (Bolatto et al. 2013) in Sec. 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.

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 4 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 4 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 4 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.

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 Sec. 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 3 (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 lev-3 https://irsa.ipac.caltech.edu/applications/DUST/ els 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.
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 degree 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.

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 4 K, and density n e = 100 cm 2 , following: Hα corr = Hα obs (Hα/Hβ) obs 2.86 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) (4), (5) and (6) correspond to those presented in Leroy et al. (2021a). Column (7) shows the vertical offset of the galaxy from the integrated main sequence of galaxies, as defined in Leroy et al. (2019). Distance measurements are presented in Anand et al. (2021) and inclinations in Lang et al. (2020). Uncertainties in columns (4), (5), (6) and (7) are on the order of 0.1 dex. Column (10) shows the area mapped by MUSE. 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 4 K), low density (10 −1 cm −3 ) 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, F. Belfiore et al. in prep). To this end, we followed the approach described in Blanc et al. (2009) where [S II] Hα DIG and [S II] Hα H II 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 0 parameters such that: 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 0 , 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 Sec. 3.2). In Sec. 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 Sec. 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, 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 Sec. 3.2), the higher uncertainty in the low Σ SFR regime has little effect in our analysis.

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 0 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 Sec. 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).

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 nonstar-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 starforming 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 dif-ferent 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.
When averaging in linear space either including nondetections 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 −2 ] ≈ 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 nondetections 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.

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 Section 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. Figure 4 shows the 2D histograms of the overall rSFMS (left), rKS (center), and rMGMS (right), and their corresponding bestfitting power laws (red dashed line and red dots, respectively) at 100 pc spatial scale.

Scaling relations using the full sample
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 Sec. 3.2) that enhance the scatter of the rSFMS and the rKS, especially at high resolution (see Sec. 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 Sec. 5.2). Similarly, variations in the normalization of these relations are expected due to differences in the assumed CO-to-H 2 conversion factor (see Sec. 5.4.2), as well as in the binning methodology (see Sec. 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-inlog, 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.

Scaling relations in individual galaxies
Here we explore galaxy-to-galaxy variations for the three scaling relations. Figures 5, 6, and 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 Sec. 2.3). The black dashed lines show the correspond- 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).
ing 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 Sec. 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 Sec. 3.2). As explained in Sec. 4.1, the rMGMS is the relation that has the lowest total detection fraction in comparison.
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 Reference α rSFMS α rKS α rMGMS spatial scale Sánchez et al. (2021) 1.01 ± 0.015 0.95 ± 0.21 0.93 ± 0.18 ∼kpc Ellison et al. (2020a) 0.68 ± 0.01 0.86 ± 0.01 0.73 ± 0.01 ∼kpc Barrera-Ballesteros et al. (2021a) 0.92 0.54 -∼2 kpc Morselli et al. (2020) 0.74 ± 0.25 0.83 ± 0.12 0.94 ± 0.29 ∼500 pc Lin et al. (2019) 1.19 ± 0.01 1.05 ± 0.01 Table 2. Summary of some previously reported values for the slopes of the rSFMS, rKS law and rMGMS. The spatial scale at which each study was carried is indicated in the last column.

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 Sec. 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.
(σ ≈ 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. 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 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. 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 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.
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).
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. (2020b), where the authors found a systematic increase of Σ mol.gas in spiral arms with respect to interarm 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 Sec. 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 Figure 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 Sec. 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.

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 Sec. 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. 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.

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.  Table 3. Slope (α) and scatter (σ) using all the available pixels in our sample, for each one of the scaling relations probed, at spatial scales of 100 pc, 500 pc and 1 kpc.  Table 4. Slope (α) and scatter (σ) for each galaxy in our sample, for each one of the three relations probed, at 100 pc spatial scale. The overall measurement considering all the valid pixels is included in the last row.

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. (2020a).
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. 2020a). 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.

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 COto-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. 2008Leroy et al. , 2013 with the apparently shorter depletion times measured on ∼ 100 pc scales (∼40−500 Myr; Evans et al. 2009;Heiderman et al. 2010;Gutermuth et al. 2011;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 bestfitting power law for each galaxy. Evans et al. 2014;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 (largescale 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 starforming 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 Chevance et al. 2020a), 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 2 conversion factor. Similarly for the rSFMS, Ellison et al. (2020b) 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, 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.
particularly at 1 kpc spatial scale, where we are no longer in a stochastic regime, for the reasons described earlier.

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

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 2 -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 Fig. 9. Variations across different galactic environments for the rSFMS (left), rKS (center) and rMGMS (right) relations. The top row shows the 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. The bottom row shows, on the x-axis the slope measured for each individual environment. 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 Sec. 4.2 for discussion). 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.

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, H.-A. 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).
Additionally, we observe a flattening at the low end of the stellar mass surface density (log Σ stellar [M kpc −2 ] 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 nondetections 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.

Treatment of non-detections and its effect on the slope of the scaling relations
As explained in Sec. 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 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.
the SFR inside 1 kpc 2 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. 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 σ DF Hα ), 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.
We parametrize the difference in slope when non-detections are excluded from the fit as compared to the fiducial value as: 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. 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.
Relation  Table 5. Slope (α) and scatter (σ) using all the available pixels in our sample, for each one of the scaling relations probed, at spatial scales of 100 pc, 500 pc and 1 kpc. Non-detections are excluded when measuring the slope.
Target  Table 6. Slope (α) and scatter (σ) for the rSFMS measured in each galaxy in our sample at spatial scales of 100 pc, 500 pc and 1 kpc. The non-detections have been excluded from the slope measurement. Figure 15 shows the change of slope when non-detections are excluded in the calculation as a function of µ DF Hα for the galaxies in our sample. The figure shows a negative correlation (PCC ≈ −0.68) between µ DF Hα 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.
The scatter in the correlation shown in Figure 15 is likely because this diagnostic does not include any aspect of the un-derlying Σ 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. submitted), the second varies much more smoothly across galactic disks. The detection fraction distribution by its own cannot completely describe the flattening of the slope. Fig. 13. Sketch to exemplify the methodology to compute the distribution of the Hα detection fraction. In each 1 kpc 2 box we calculated the filling factor (DF Hα ) defined as the fraction of pixels with nonzero SFR.

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 σ DF Hα . The figure shows a clear positive correlation (PCC ≈ 0.72), meaning that galaxies with larger spread (i.e., larger σ DF Hα ) 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 σ DF Hα (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 σ DF Hα values show either a flatter DF Hα distribution or a concentration of their values around zero. In these scenarios, nondetections do not contribute to steepen the relation in the same way, since different mass ranges are similarly affected. 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.

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 ) and HERACLES (Leroy et al. 2009) surveys. Along the same line, Schruba et al. (2019) -studying 8 local galaxies-and Sun et al. (2020a) -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. (2020a) 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;Dib et al. 2017;Shi et al. 2018;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 Sec.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).
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. 2019aBacchini et al. ,b, 2020Sun et al. 2020b). 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).

Systematic effects
In this section, we test how some of our assumptions impact our measurements of the scaling relations probed. In Sec. 5.4.1, we derive the rSFMS and rKS relations after removing diffuse ionized gas (DIG) emission. In Sec. 5.4.2, we use a constant α CO to measure the rKS and the rMGMS.

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 Sec. 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 F. Santoro et al. (in prep.). In short, H II regions are identified 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 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. 2020a). 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. 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 clas-I. Pessa et al.: Resolved scaling relations with PHANGS b 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 Sec. 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 Sec. 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). sified 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. 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. 5.4.2. Impact of the chosen CO-to-H 2 conversion factor As explained in Sec. 2.2, our fiducial CO-to-H 2 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 −2 (K km s −1 ) −1 (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 2 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 2 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 2 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).

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 Sec. 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 7 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.
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 Sec. 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). 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 Sec. 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.

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 starforming cycle. The same was recently reported in Ellison et al. (2020a) 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 I. Pessa et al.: Resolved scaling relations with PHANGS Fig. 21. rKS (top) and rMGMS (bottom) measured using all pixels from our sample and assuming a constant CO-to-H 2 conversion factor. This conversion factor leads to a flatter rKS and a steeper rMGMS (see Sec. 5.4.2 for a discussion). The best fit power law obtained with the fiducial approach is indicated with the gray line. 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. 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.

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 2 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 Sec. 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 Secs. 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 Sec. 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 nondetections than the inner region with a higher detection fraction, resulting in a steeper rSFMS (see Secs. 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 Sec. 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. 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 (µ DF CO ). 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 (σ DF CO ). 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 Section 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 µ DF CO for all galaxies. On the other hand, the impact in σ DF CO 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 µ DF CO instead. Figure D.3 shows the difference in the rMGMS slope with the overall slope as a function of µ DF CO , 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.  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. 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.
Article number, page 25 of 27