A new census of dust and polycyclic aromatic hydrocarbons at $z=0.7-2$ with JWST MIRI

This paper is based on the MIRI imaging survey of Systematic Mid-infrared Instrument Legacy Extragalactic Survey (SMILES; Rieke et al. 2024; Alberts et al. 2024). SMILES is a MIRI US-GTO programme (PID 1207, PI: George Rieke) that surveys 34 arcmin² (using 15 MIRI pointings in a 3$\times$ 5 mosaic) centred on the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006). The survey is centred at 03h32m34.80s -27d48m32.60s and includes the eight broad bands of the MIRI Imager from 5.6 to 25.5 $\mu$m. A standard MIRI 4-point dither pattern was adopted to improve the point spread function (PSF) sampling and mitigate detector artifacts and cosmic-ray events. We used a FASTR1 readout pattern. Each pointing had a total science exposure time of 2.17 hours, divided over eight bands. The exposure times per band were designed to optimise the detection of typical active galactic nucleus (AGN) and star-forming galaxies down to an SFR of $\sim 10$ $M_{\odot}$yr^-1. The observations were carried out in December 2022 and January 2023, which achieved a 5$\sigma$ point source sensitivity of 0.11 $\mu$Jy (AB mag $=$ 26.3), 0.20 $\mu$Jy (25.6), 0.48 (24.7), 0.67 (24.3), 0.79 (24.2), 2.42 (22.9), 3.47 (22.6), and 16.91 (20.8) in F560W (754 s exposure time), F770W (866 s), F1000W (644 s), F1280W (755 s), F1500W (1121 s), F1800W (755 s), F2100W (2186 s), and F2550W (832 s)¹¹1Our derived sensitivities are $\sim 2-3\times$ deeper than the sensitivities predicted by ETC (version 3.0) for our observation setup.. The achieved sensitivities are an order of magnitude deeper than the confusion limit of Spitzer/MIPS 24  $\mu$m (56 $\mu$Jy; Dole et al. 2004) at similar wavelengths. The optimal combination of depth, spatial resolution, coverage area, and multi-wavelength capabilities makes SMILES well suited to investigating mid-IR dust emissions in galaxies during the cosmic noon epoch. In Figure 1, we show the advancement of this dataset compared to the existing IR data to study cosmic noon galaxies.

We started from the MIRI-selected sample (3$\sigma$ detected in F560W or F770W) and used the following criteria to define our final sample:

1. 1.

   $z=0.7-2.0$ to ensure adequate PAH coverage by MIRI data.

2. 2.

   $3\sigma$ detection in at least one MIRI band longwards of 15 $\mu$m.

3. 3.

   The removal of quiescent galaxies based on their rest-frame $U-V$ and $V-J$ colours (following Equation 3 in Schreiber et al. 2015). For a cleaner selection, we also removed galaxies with specific SFR (sSFR, SFR/$M_{*}$ from SED modelling, see Section 2.3) $<10^{-10}$ yr^-1.

4. 4.

   The removal of objects with unusual or incorrect photometry from visual inspections of the SED fits.

This selection results in 443 galaxies (‘full’ sample). The stellar mass and SFR distribution is shown in Figure 2. Specifically for the PAH analysis in this work, we further limited the sample to those with ‘robust’ PAH mass fraction (q_PAH, see Section 2.3.1). This category is defined as galaxies with either $>3\sigma$ q_PAH values or q_PAH error of less than 0.5%. The second condition allows the inclusion of near-zero PAH fractions that are still relatively robust. The rest of the sample is denoted as ‘low confidence’ q_PAH throughout the paper. The robust q_PAH sample has 216 objects. Figure 3 shows the distributions of redshift, F444W (tracing emission from old stellar populations), and F1500W (tracing PAH and dust emission) of the robust and full samples. The robust q_PAH sample is skewed towards brighter galaxies (and relatively more massive and dusty), by definition. Even though this sub-sample is biased, the separation is necessary to have a sample with reliable PAH fraction estimates. However, we show the results and trends for the low-confidence sample as well as the robust sample in the paper. The properties of the robust q_PAH sample are provided in an online repository²²2Online repository for supplementary data: https://zenodo.org/records/12671075 and a sample table is shown in Table 1.

We adopted the PROSPECTOR SED fitting code (Johnson et al., 2021) to fit the UV to mid-IR photometry. PROSPECTOR is a highly flexible SED fitting code that is based on Bayesian forward modelling and Monte Carlo sampling of the parameter space using FSPS stellar populations (Conroy, 2013). We have adopted the ‘current’ or ‘surviving’ stellar mass throughout this paper, which takes into account losses throughout the evolution of stars in the galaxy via AGB winds, supernovae, etc. This is different from the mass of stars that are formed in the lifetime of the galaxy. In our sample, the ratio of the surviving to the formed stellar mass is on average 0.69 with a scatter of 0.04 (for a Chabrier IMF). Below, we first discuss the assumptions on modelling the dust emission (Section 2.3.1), then discuss the validity of the derived PAH fractions (Section 2.3.2). Finally, in Section 2.3.3, we discuss the remaining assumptions in the SED fitting. We note that the main results of this paper (i.e. PAH and dust emission) are mainly constrained by the assumed dust models (Draine & Li, 2007), and therefore using a different SED fitting code is not expected to significantly change the results.

Dust sizes and surface areas were estimated from MIRI data as a proxy for where the bulk of stars are forming in galaxies. These measurements were later used to calculate SFR surface densities. We calculated dust sizes from MIRI F1500W images with 0.5” resolution, tracing rest-frame $5-9$  $\mu$m. The choice of the F1500W filter to measure dust sizes was made as this filter provides a good compromise between wavelength and spatial resolution.

We used WebbPSF (Perrin et al., 2014) to model the F1500W PSF. In order to estimate sizes, we used GALFIT (Peng et al., 2002, 2010) to fit the spatial profile of each source. Each one was modelled twice: once with a point-source model and once with a single-Sersic model. For a complete morphological analysis, more complicated models are required but for our current purposes these simple models provide reasonable sizes for the majority of the sample. Galaxies that fit better (smaller reduced $\chi^{2}$) with a Sérsic model than the point-source model are considered spatially resolved. Sources that are best fit with the point source model or with a Sérsic profile but whose half-light radius is smaller than the 50% encircled energy radius for the appropriate filter are considered unresolved. For the unresolved sources, we used the 50% encircled energy size of the PSF as an upper limit on their effective radius. We visually inspected all fits and removed the bad fits (e.g. if the source had signs of merger or nearby bright companions). This results in 105 F1500W good fits for extended objects in our sample.

In Figure 7, we show the relationship between stellar mass and q_PAH (left panel) and the 7.7  $\mu$m PAH luminosity fraction relative to total IR luminosity emission ($F_{7.7}/F_{\rm{TIR}}$; right panel). $F_{7.7}$ is the integrated PAH flux ($\int(F_{\lambda}-F^{\rm{continuum}}_{\lambda})\leavevmode\nobreak\ {\rm d}\lambda$) between 6.9 and 9.7 $\mu$m and is corrected for the continuum emission, following the recipe of Draine et al. (2021), and the total IR (TIR) power is the integrated emission ($\int F_{\lambda}\leavevmode\nobreak\ d\lambda$) from 3 to 1100 $\mu$m of the best-fit SED model.

The PAH ‘strength’, defined as either the PAH mass or emission fraction, is positively correlated with stellar mass with a correlation coefficient of $\rho_{\rm Pearson}=0.5$ for the robust q_PAH sample and 0.6 for the full sample (with p-values of $\ll 1$). We also estimate the gas-phase metallicity ($12+\log{\rm{(O/H)}}$) based on stellar masses and the $z\sim 1.5$ O3N2-based mass-metallicity relation (MZR) of Topping et al. (2021) and the Pettini & Pagel (2004) calibration, shown on the top axis of Figure 7. The median trend of q_PAH with the MZR-estimated metallicities closely resembles the local galaxies trend from Draine et al. (2007) (the horizontal lines in Figure 7-left). The two q_PAH median values of the local galaxies were reported at metallicities below and above $12+\log{\rm{(O/H)}}$$=8.1$. The rapid change in q_PAH happens at $\sim 0.1-0.2$ dex higher oxygen abundances in our $z\sim 1.3$ sample compared to that at $z\sim 0$. However, given the uncertainties in metallicity calibrations (see also the next section), we conclude that the behaviour of the $z\sim 1.3$ sample is consistent with that at $z\sim 0$.

There is a large scatter in q_PAH at intermediate masses, such that in the mass range of $\log(M_{*}/M_{\odot})\sim 9.5-10.5$, there are galaxies with q_PAH of $<0.5\%$ to $>4\%$. The majority of these galaxies seem to have young ages and high specific SFRs (sSFR). To better understand the physical process driving the scatter, we also calculate partial correlation coefficients between q_PAH, $M_{*}$, and a third galaxy parameter, including sSFR, SFR, age, optical depth, obscured SFR fraction, and SFR surface density. However, we do not find a significant secondary dependence on any of these parameters in the q_PAH-$M_{*}$ correlation of the robust sample. If we limit the sample to only q_PAH values with $>3\sigma$ significance, the specific SFR (sSFR) and age reduce the q_PAH-$M_{*}$ correlation coefficient of 0.4 to a partial coefficient of 0.25, indicating a potential secondary influence on the q_PAH-$M_{*}$ correlation, such that at a given stellar mass, galaxies with higher sSFRs and younger ages have lower PAH strengths.

The relation between the PAH fraction and metallicity has been extensively studied in the local Universe (Marble et al., 2010; Engelbracht et al., 2005; Madden et al., 2006; Smith et al., 2007; Galliano et al., 2008; Draine et al., 2007; Rémy-Ruyer et al., 2015; Aniano et al., 2020; Li, 2020; Chastenet et al., 2023) and also at cosmic noon (Shivaei et al., 2017, using Spitzer data). To investigate this relation with the deeper MIRI data at $z>0$, we estimate metallicities from the fundamental metallicity relation (FMR) between stellar mass, SFR, and metallicity. The FMR is shown to be established and invariant up to $z\sim 3$ (Cresci et al., 2019; Sanders et al., 2021). Sanders et al. (2021) estimated the scatter in $12+\log{\rm{(O/H)}}$ from FMR to be 0.06 dex at $z\sim 2.3$, and 0.04 dex at $z\sim 0$. This is typically smaller (or at most comparable) to the systematic uncertainty in strong-line metallicity calibrations (e.g. $\sim 0.2$ dex average systematic uncertainty between the Pettini & Pagel (2004) and the Pilyugin & Thuan (2005) calibrations in Aniano et al. (2020), see their Figure 1). We adopt the FMR presented in Sanders et al. (2021) to estimate metallicities from our mass and SFR values. The adopted FMR is in agreement with those derived from direct-method metallicities (Andrews & Martini, 2013) and also with the best-fit for local galaxies in Curti et al. (2020). We adopt a Monte Carlo simulation by randomly perturbing stellar mass and SFR within their uncertainties for 500 times, and derive the metallicities according to the FMR relation. The final metallicity and uncertainty of each object is the median and standard deviation of the 500 realisations. The average measurement uncertainty of the inferred $12+\log{\rm{(O/H)}}$ is 0.04 dex, comparable to the FMR scatter in $12+\log{\rm{(O/H)}}$.

We limit the $z\sim 0$ comparison samples to two studies that adopted, similar to this study, the Draine & Li (2007) models to derive q_PAH and have homogeneous metallicity calculations: Aniano et al. (2020) (hereafter, A20) and Rémy-Ruyer et al. (2015) (hereafter, R15). A20 used two strong-line metallicity calibrations for the KINGFISH sample: 1) the ‘PT’ or Pilyugin & Thuan (2005) method, adopted from Moustakas et al. (2010), and 2) the ‘PP04’ or Pettini & Pagel (2004) based on [Nii]/H$\alpha$. They showed that the two metallicities can differ by as much as 0.4 dex (on average $\sim 0.2$ dex), reflecting the large uncertainty in the calibrations. R15 incorporated the Dwarf Galaxy Survey (DGS; Madden et al. 2013) into the KINGFISH which extended the sample to lower metallicities. They adopted the PT metallicity calibration.

Figure 8 shows the q_PAH trend with metallicity. There is a positive correlation (Pearson $\rho=0.55$ with p-value $\ll 1$) with a large scatter at high metallicities. The trend is very similar to the $z\sim 0$ trend from A20 (assuming their preferred PP04 metallicity³³3A20 and Hunt et al. (2016), who used the same metallicities in their analysis, preferred the PP04 calibration as it shows tighter scaling relations with other calibrations and with q_PAH.). A20 approximated the q_PAH-metallicity behaviour by both a linear function and a step function. The median values below and above their threshold metallicity, $12+\log{\rm{(O/H)}}$$\sim 8.25$, and the slopes of the linear fits ($3.6\pm 0.7$ and $4.0\pm 0.5$ at $z\sim 0$ and 1.3, respectively) for the two samples suggest a non-evolving PAH fraction-metallicity relation from $z\sim 0-2$.

In Figure 9 (top panel), we compare our results with the $z\sim 0$ data over a wider range of metallicities by combining both the A20 and R15 data (using the PT metallicity). The q_PAH values of R15 extend to larger values compared to A20 and ours, because they set a higher upper limit for the q_PAH prior in their fittings. We also include A20 data for the PP04 metallicity, for completeness. To interpret and compare the datasets, we adopt a moving (running) median to smooth the data without assuming a parametric functional form a priori. All three trends show a constant low ($\lesssim 1\%$) q_PAH at low metallicities ($12+\log{\rm{(O/H)}}$$\lesssim 8.1-8.2$ or $\sim 0.3\,Z_{\odot}$ assuming $12+\log{\rm{(O/H)}}$$=8.69$ for Solar metallicity; Asplund et al. 2009), then a linear increase in q_PAH with metallicity that ends with a plateau (q_PAH $\sim 3.4\%$) at $12+\log{\rm{(O/H)}}$$\gtrsim 8.4$ ($\gtrsim 0.5\,Z_{\odot}$). We show the fits to the moving medians in the lower panel of Figure 9. The constant q_PAH at high metallicities confirms the successful usage of PAH strengths as IR luminosity and obscured SFR indicators for moderately massive and metal-rich galaxies based on local calibrations (e.g. Calzetti et al., 2007; Shipley et al., 2016). Our findings are also in general agreement with Marble et al. (2010) based on the $z\sim 0$ SINGS data, who found a linear relation between PAH luminosity fraction and metallicity. These observations show that, at least up to $z\sim 2$, metallicity is a good indicator of the ISM properties that affect the balance between the formation and destruction of PAHs in star-forming galaxies, and the responsible processes maintain a constant PAH abundances of $\sim 3.4$% above a metallicity of $\sim 0.5\,Z_{\odot}$.

Below $\sim 0.5\,Z_{\odot}$, we see that there is a linear decrease in the PAH fraction, which has a similar slope between the $z\sim 0$ and 1.3 samples. This universal PAH fraction-metallicity relation can potentially be calibrated to estimate the metallicities of metal-poor galaxies at high redshifts using IR imaging surveys alone. The source of the decreasing PAH fraction with decreasing metallicity can be explained by the destruction of PAHs at low metallicities, either through radiative destruction in less shielded low-metallicity environments (Madden et al., 2006; Hunt et al., 2010; Narayanan et al., 2023), or supernova shocks (Seok et al., 2014). Some studies also explain it by insufficient production of PAHs due to a delayed injection into the ISM by AGB stars (Galliano et al., 2008) or more efficient production at high metallicities through shattering of large grains (Seok et al., 2014). Future cosmological simulations with on-the-fly dust evolution simulations can shed light on the complex interplay of the formation and destruction processes of these grains in the ISM.

In Figure 10 we show q_PAH as a function of other galaxy parameters: age, SFR, stellar mass, obscured luminosity fraction ($f_{\rm obscured}$), and SFR surface density ($\Sigma_{\rm SFR}$). Obscured luminosity fraction is calculated as the ratio of the dust-obscured UV (1550Å) luminosity to the intrinsic UV luminosity (see Section 4). Sizes are calculated in the F1500W filter, tracing rest-frame 5 to 9  $\mu$m dust emission, for galaxies with spatially resolved profiles (see Section 2.4). In the literature, $\Sigma_{\rm SFR}$ is often calculated based on the area derived from a shorter-wavelength size measurement that is not a good proxy for where the bulk of the stars are forming. Here, surface area is based on dust sizes (from F1500W) that closely reflects the star-forming region, and hence, the resulting SFR surface densities are more reliable.

We show the moving medians for the robust (black curve) and the full sample (grey curve), which includes the robust sample, in Figure 10. In the robust sample, q_PAH correlates with age, stellar mass, and obscured fraction. The q_PAH-mass correlation is the strongest among these, and is likely a byproduct of the mass-metallicity relation (as was discussed in previous sections). The correlation with age can be expected as a stellar evolutionary effect if the AGB stars are the primary production channel of PAHs (Galliano et al., 2008; Li, 2020) and, in fact, it is also seen at $z\sim 2$ from MIPS data (Shivaei et al., 2017). However, this correlation disappears when the full sample is considered. The correlation with obscured fraction is interesting as it may suggest PAH destruction due to reduced shielding by dust grains in systems with low obscured fractions. We also note the lack of a significant correlation with SFR surface densities for the spatially resolved galaxies, which may suggest that the fraction of PAH-to-total dust mass is insensitive to the star formation intensity. Even including lower limits on SFR surface density for the unresolved sources (triangles in Figure 10) does not suggest a correlation with SFR surface density. Further analysis of this trend using size measurements from other filters that have higher resolution but still trace the dust emission will be done in a future work.

As most of the metals in the ISM are depleted onto dust grains that attenuate starlight, it is expected that obscured luminosity fraction (ratio of attenuated to intrinsic luminosity) correlates with metallicity, and hence stellar mass (e.g. Reddy et al., 2010; Whitaker et al., 2017), which is also manifested in the relationship between dust attenuation and mass (e.g. Garn et al., 2010; Shivaei et al., 2020b; Bogdanoska & Burgarella, 2020).

In Figure 11, the obscured fraction is shown as a function of stellar mass. We investigate the average shape and scatter of the relationship at $M_{*}>4\times 10^{8}$ $M_{\odot}$, where the sample is complete. The mass completeness is estimated based on the detection limits of 5.6  $\mu$m and 15  $\mu$m data and assuming a mass-to-light ratio of 0.6 at 2.2 $\mu$m (McGaugh & Schombert, 2014)⁴⁴4The exact completeness varies with redshift from 0.7 to 2.. To study the average trend, we fit the data in three ways: a moving median function⁵⁵5The moving weighted-average fit shows an identical trend to the moving median curve., the functional form of $f_{\rm obsc}$-$M_{*}$ suggested by Whitaker et al. (2017),

where $f_{\rm obscured}$ is the obscured fraction, and $a$ and $b$ are the fit parameters, and the following function that best matches the moving median trend:

where erf is the error function, and $f_{\rm max}$, $f_{\rm min}$, $\mu$, and $\sigma$ are the fitting parameters. Using an orthogonal distance regression (ODR) fitting procedure (python scipy.odr), considering the uncertainties on both the obscured fraction and mass, the fit using Equation 2 matches the Whitaker et al. (2017) curve closely. However, our fit is highly biased by the few highly obscured galaxies with very small uncertainties. The moving (running) weighted-average curve, where the weight is the inverse square of uncertainty in $f_{\rm obsc}$, also shows how this small population of highly obscured galaxies bias the fit (dotted curve in Figure 11). The running median (dashed curve) and the Equation 3 fit (black line) better represent the behaviour of $f_{\rm obsc}$ as a function of mass. The best fit parameters are listed in Table 2.

We do not see any significant redshift evolution in this relation, as is shown in the right panel of Figure 11, in agreement with the findings of Whitaker et al. (2017) that at a fixed mass the obscured fraction is constant from $z\sim 0.5$ to 2.5. Based on our best-fit curve (Eq. 3) that matches the moving median, on average, the UV luminosity is dominantly obscured ($f_{\rm obsc}$$>50\%$) above the stellar mass of $\sim 5\times 10^{9}$  $M_{\odot}$. Interestingly, owing to the large scatter in $f_{\rm obsc}$ at any mass, there are galaxies with dominantly obscured emission even down to $M_{*}\sim 5\times 10^{8}$  $M_{\odot}$. We investigate the scatter in $f_{\rm obsc}$ at fixed mass in the next section. Above the stellar mass of $2\times 10^{10}$  $M_{\odot}$, more than 80% of the emission is obscured. However, the SMILES survey is likely not representative of the population of heavily obscured (and rare) galaxies owing to its relatively small area.

Whitaker et al. (2017) studied the Spitzer/MIPS samples and found that the $f_{\rm obscured}$-$M_{*}$ relation does not change from $z=0$ to $2.5$. They fit the UV-near-IR photometry with Bruzual & Charlot (2003) models and used the Dale & Helou (2002) template to convert 24 $\mu$m flux densities to total IR luminosities. Then assuming the Kennicutt (1998) and Bell et al. (2005) calibrations, they converted UV and IR luminosities to SFRs to estimate $f_{\rm obsc}=\frac{\rm SFR_{IR}}{\rm SFR_{IR}+SFR_{UV}}$. While the two relationships are in relatively good agreement within the scatter of the relation in Figure 11, there are a few differences between the Whitaker et al. (2017) study and this work that contribute to some disagreement between the two relationships. Whitaker et al. (2017) used a single log-averaged IR template from Dale & Helou (2002) (i.e. a fixed IR SED shape) to convert PAH luminosity to total IR luminosity for their sample at $z\sim 0.5-2.5$. While at high metallicities and masses, galaxies may have a more uniform IR SED shape, and hence, the assumption of a single IR SED would be valid, that is not always the case at low metallicities. As we showed in the previous section, below $\sim 0.5Z_{\odot}$ the PAH fraction decreases uniformly from $z\sim 0-2$. Therefore, using the high-metallicity end calibrations would overestimate the IR luminosities at lower metallicities (and masses; see also e.g. Shivaei et al. 2017, 2022). Therefore, the single IR SED assumption in Whitaker et al. (2017) might have overestimated their IR luminosity values at lower masses. Additionally, the different area and sensitivity of SMILES compared to the MIPS/Spitzer surveys can play a role in the observed discrepancy between this study and the Whitaker et al. (2017) study. SMILES is more sensitive to less obscured galaxies at a fixed mass relative to the MIPS samples. In fact, the fit to our data using the Whitaker et al. (2017) formalism and an ODR fitting procedure (i.e. taking into account the uncertainties on $f_{\rm obsc}$) closely follows the Whitaker et al. (2017) curve. This is because a few highly obscured galaxies with very low $f_{\rm obsc}$ uncertainties dominate the fit – a similar bias that a MIPS-based sample would have. Furthermore, SMILES has a smaller area compared to the Spitzer/MIPS surveys, and hence, is likely missing a (rarer) population of massive and highly attenuated galaxies.

As was mentioned in the previous section, there is a large scatter in $f_{\rm obsc}$ at a given mass, such that there are galaxies with significant obscured UV emission down to stellar masses of $\sim 5\times 10^{8}$  $M_{\odot}$. The scatter in the obscured fraction at a given mass increases from 0.1 dex at $M_{*}\sim 10^{11}$  $M_{\odot}$ to $\sim 0.5$ dex at $M_{*}\sim 10^{9}$  $M_{\odot}$. We calculate secondary correlations in the $f_{\rm obsc}-M_{*}$ relation with sSFR, age, and SFR surface density using partial correlation coefficients. The surface density is measured using dust sizes from F1500W data (Section 2.4). In those calculations the sample is limited to galaxies with spatially resolved discs in F1500W. The spatially resolved sample is limited to galaxies with $M_{*}>10^{9.5}$  $M_{\odot}$, and inevitably biased against the extremely compact star-forming galaxies. The stellar mass distribution of the spatially resolved sample above $M_{*}>10^{9.5}$  $M_{\odot}$ is very similar to that of the parent sample with 50% of the spatially resolved sample at $M_{*}<10^{10.2}$  $M_{\odot}$. Among the sSFR, age, and SFR surface density parameters, we find a relatively strong correlation between the obscured fraction and SFR surface density at a fixed stellar mass, but no significant correlations with age and SFR at constant mass.

Figure 12 shows the obscured fraction versus SFR surface density ($\Sigma_{\rm SFR}$, left) and F1500W effective radii (middle). $\Sigma_{\rm SFR}$ is calculated using dust surface area (at rest-frame $8-9$  $\mu$m) that closely traces the active star-forming region. We divide the sample in three stellar mass bins to investigate the trend at a fixed mass. The average values are shown with large symbols in the figure.

There is a clear correlation between $f_{\rm obsc}$ and $\Sigma_{\rm SFR}$ at all masses with Pearson correlation coefficient of $\sim 0.5$ (p-value $\ll 1$). We fit the average trends in each mass bin with a linear function. As was expected, the overall normalisation increases with increasing mass; that is, as was shown in the previous figures, $f_{\rm obsc}$ increases with increasing mass. However, not only does $f_{\rm obsc}$ clearly increase with $\Sigma_{\rm SFR}$ (and decreases with size) at a given mass, but the rate at which $f_{\rm obsc}$ increases with $\Sigma_{\rm SFR}$ (i.e. the slope) appears to be similar across the mass range of the sample. The fit parameters are presented in Table 3. We discuss the possible explanation of this behaviour in the next Section. While the overall trend with size (effective radius) is consistent, the slope of the $f_{\rm obsc}$-size relation varies among the three mass bins and the evolutionary behaviour is not as clear as with $\Sigma_{\rm SFR}$. We also note that the width of the $f_{\rm obsc}$ distribution decreases as the mass increases, as is shown with the error bars on the average values in Figure 12. The standard deviation of $f_{\rm obsc}$ in the highest $\Sigma_{\rm SFR}$ bin from low to high mass is 0.20, 0.15, 0.07, respectively.

The right panel of Figure 12 confirms that the $f_{\rm obsc}$ increasing trend with $\Sigma_{\rm SFR}$ is not a byproduct of increasing SFR. While there is a general correlation between SFR and size for the full sample, the average SFR does not vary with size in the fixed mass ranges. Therefore, the underlying cause of $f_{\rm obsc}$ dependence on $\Sigma_{\rm SFR}$ at a given mass is the compactness of the star-forming region, and not only an increase in the rate of star formation.

Using submm or mm data for samples with generally much higher stellar masses and SFRs than presented here, studies have shown that star formation occurs in more compact regions compared to the stellar mass sizes (e.g. Murphy et al., 2017; Jiménez-Andrade et al., 2021; Hodge et al., 2016; Elbaz et al., 2018; Fujimoto et al., 2017; Rujopakarn et al., 2019; Hamed et al., 2023a). Additionally, it has been argued that starburst galaxies on the main-sequence with particularly short depletion timescales exhibit higher SFR surface densities with more compact mm dust emission (Elbaz et al., 2018; Gómez-Guijarro et al., 2022) and radio emission (Jiménez-Andrade et al., 2019), and hotter dust (Gómez-Guijarro et al., 2022) compared to other main-sequence galaxies at the same mass and redshift. In the regime of lower SFR surface density ($\log(\Sigma_{\rm SFR}/M_{\odot}{\rm yr^{-1}}\leavevmode\nobreak\ {\rm kpc^{-2}})\sim-1$ to 1 ) and mass ($\log(M_{*}/M_{\odot})\sim 9.5$ to 11), we observe the same relation for purely star-forming galaxies (no AGN contamination) that PAH dust sizes are more compact with higher SFR surface densities at a fixed mass (right panel of Figure 12).

We also see that the obscured emission is more dominant in compact galaxies with higher $\Sigma_{\rm SFR}$. This observed trend is illustrated by a simple cartoon in Figure 13. For an extended star-forming galaxy the dust-star geometry is more complex: some population of stars are only affected by the diffuse ISM dust, while others can be trapped in their birthcloud dust and highly attenuated. If the star-forming region is compressed into a small area, the galaxy will have a more homogeneous dust distribution with higher covering fraction that attenuates most of the light emitted from the stars. Hence, the obscured fraction of the emitted light is higher in galaxies with compact star-forming regions. This explanation works for closed-box systems. However, strong outflows and feedback may clear pathways through the ISM of the compact star-forming regions and expose the starlight without attenuation. This may explain the larger variation in $f_{\rm obsc}$ in the low mass bins at any $\Sigma_{\rm SFR}$: low-mass galaxies have highly stochastic star formation histories (Hopkins et al., 2014), and massive stars produced in the vigorous starbursts can destroy their birth clouds via feedback (e.g. winds, supernovae; Naidu et al. 2022) to make sightlines clear of dust. This will generate a diversity of $f_{\rm obsc}$ from galaxy to galaxy at low masses, while in massive galaxies dust stays in the compact core with high covering fraction. In our data, the scatter (standard deviation) of $f_{\rm obsc}$ in the lowest mass bin is $15-20$% at a fixed $\Sigma_{\rm SFR}$ and decreases to $8-12$% in the highest mass bin.

The relationship of obscured fraction with compactness of star-forming region can also explain the spread in the locus of galaxies in the IRX-$\beta$ diagram. The IRX-$\beta$ relation is an empirical calibration of UV stellar continuum slope ($\beta$) based on the ratio of IR to unobscured UV luminosity (IRX; Meurer et al. 1999). The location of galaxies in this diagram depends on the slope of the attenuation curve in the UV and the intrinsic UV slope (e.g. Reddy et al., 2018; Shivaei et al., 2020a; Salim & Narayanan, 2020; Hamed et al., 2023b). Figure 14 shows the distribution of our sample in the IRX-mass and IRX-$\beta$ diagrams. Galaxies with higher $\Sigma_{\rm SFR}$ at a given mass or $\beta$ have higher IRX, and preferentially shallower attenuation curves similar to that of the Calzetti et al. (2000) curve (middle panel; similar trends are also seen at $z\sim 0$, Johnson et al. 2007). The preference for shallower curves is potentially due to both a higher average optical depth and a lower fraction of unobscured evolved stars. Radiative transfer models have shown that a high dust optical depth flattens the dust attenuation curve (Chevallard et al., 2013) and the right panel of Figure 14 shows that $\Sigma_{\rm SFR}$ in our sample is correlated with the dust optical depth (Pearson $\rho=0.65$ with p-value $\ll 1$). Additionally, as is shown in the schematic Figure 13, galaxies with more concentrated star-forming regions (higher $\Sigma_{\rm SFR}$) have lower fraction of evolved stars unobscured by dust, as these stars would not have sufficient time to migrate away from their dusty birth clouds. This dust-star geometrical effect also tends to flatten the dust attenuation curve of galaxies in which the evolved stars dominate the optical luminosity (Narayanan et al., 2018). The inclination (axis ratio) of galaxies can also play a role in the observed trends, which will be investigated in future studies.