Understanding the Universal Dust Attenuation Scaling Relation of Star-Forming Galaxies

Star-forming galaxies (SFGs) adhere to a surprisingly tight scaling relation of dust attenuation parameterized by the infrared excess (IRX=$L_{\rm IR}/L_{\rm UV}$), being jointly determined by the star formation rate (SFR), galaxy size ($R_{\rm e}$), metallicity ($Z$/Z$_\odot$) and axial ratio ($b/a$). We examine how these galaxy parameters determine the effective dust attenuation and give rise to the universal IRX relation, utilizing a simple two-component star-dust geometry model in which dust in the dense and diffuse interstellar medium (ISM) follows exponential mass density profiles, connected with but not necessarily identical to the stellar mass profiles. Meanwhile, empirical relations are adopted to link galaxy properties, including the gas--star formation relation, the dust-to-stellar size relation, as well as the dust-to-gas ratio versus metallicity relation. By fitting a large sample of local SFGs with the model, we obtain the best-fitting model parameters as a function of metallicity, showing that the two-component geometry model is able to successfully reproduce the dependence of IRX on SFR, $R_{\rm e}$, $b/a$ at given $Z$/Z$_\odot$, as well as the dependence of power-law indices on metallicity. Moreover, we also retrieve constraints on the model geometry parameters, including the optical depth of birth clouds (BCs), BC-to-total dust mass fraction, BC covering factor of UV-emitting stars, and star-to-total dust disc radius ratio, which all evolve with galaxy metallicity. Finally, a consistent picture of how the star-dust geometry in SFGs evolves with galaxy metallicity is discussed.


INTRODUCTION
Dust, as a crucial component of the interstellar medium (ISM), plays a significant role in driving the cycles of baryons and energy.Firstly, it acts as a channel for converting radiation pressure into mechanical energy and outflows (Thompson et al. 2015;Naddaf & Czerny 2022).Secondly, dust serves as an important coolant for star and planet formation, and an interface for the formation of molecules (Galliano, Galametz, & Jones 2018).Additionally, dust has a profound influence ⋆ xzzheng@pmo.ac.cn on the spectral energy distributions (SED) of galaxies by absorbing the stellar ultraviolet (UV) to optical radiation and thermally re-emitting at the infrared (IR) wavelengths.This alteration of the SED significantly affects galaxy observables (Calzetti et al. 2000;Conroy 2013).Dust attenuation, defined as the effective sight-line absorption of light by dust in a galaxy (also referred to as 'obscuration'), 1 depends on both the dust content and the geometry between dust and stars.It is closely intertwined with star formation, chemical enrichment and structural growth of the galaxy.Understanding dust attenuation and the regulation by dust of physical processes across cosmic time is crucial not only for improving measurements of intrinsic galaxy properties, but also for unravelling the connections between dust, gas, metals, and stars in accordance with the structural evolution of galaxies (Salim & Narayanan 2020).
For star-forming galaxies (SFGs), various fundamental scaling relations involving gas, metals, and stars have been established up to a redshift of z ∼ 10, revealing their rapid evolution across this redshift range.One of these scaling relations is the star formation main sequence, which shows a nearly linear increase in star formation rate (SFR) with stellar mass (M * ) at a given redshift, while globally decreasing in SFR as redshift decreases (Noeske et al. 2007;Wuyts et al. 2011;Whitaker et al. 2012;Guo, Zheng, & Fu 2013;Speagle et al. 2014;Schreiber et al. 2015;Katsianis et al. 2020;Thorne et al. 2021;Leja et al. 2022;Kouroumpatzakis et al. 2023;Popesso et al. 2023).The density of star formation in SFGs is regulated by the density of cold gas, following the Kennicutt-Schmidt (K-S) Law (Σgas-ΣSFR) (Kennicutt 1998).Another important scaling relation is the stellar mass-metallicity (M * -Z) relation, which reveals an increase in gas-phase metallicity with stellar mass among SFGs.However, this relation also shows a rapid decline in metallicity at fixed M * as redshift increases (Tremonti et al. 2004;Erb et al. 2006;Mannucci et al. 2010;Zahid et al. 2012Zahid et al. , 2014;;Maiolino & Mannucci 2019;Sanders et al. 2021;Nakajima et al. 2023;Curti et al. 2023).In addition, SFGs are predominantly characterized by disc morphologies, and their sizes, quantified by the half-light radius (Re), typically increase with stellar mass.The stellar mass-size (M * -Re) relation of SFGs demonstrates a rapid size increase over time for SFGs with a given M * (van der Wel et al. 2014a;Suess et al. 2019;Mowla et al. 2019;Nedkova et al. 2021;Ono et al. 2023;van der Wel et al. 2024).These scaling relations provide a framework for understanding the evolution of and interplay between gas, metals, stars, and structural properties of SFGs from high redshifts to the present day.
On the other hand, dust attenuation serves as an indicator of effective dust columns in SFGs, which is closely related to gas density and metallicity.Observationally, dust attenuation in SFGs is often quantified using the IR to UV luminosity ratio, known as the IR excess (IRX≡LIR/LUV).Generally, SFGs with higher star formation rates (SFR) tend to exhibit higher infrared luminosities and dust temperatures (Dale & Helou 2002;Rieke et al. 2009;Schreiber et al. 2018), resulting in increased dust attenuation (Whitaker et al. 2014).The dust attenuation in SFGs has been found to positively correlate with stellar mass, SFR, surface density of SFR, UV continuum slope, and gas-phase metallicity (Kong et al. 2004;Martin et al. 2005;Johnson et al. 2007;Garn & Best 2010;Zahid et al. 2013;Battisti, Calzetti, & Chary 2016).However, investigations on local SFGs have revealed that the dependence of dust attenuation on SFR (as well as galaxy inclination angle) weakens with decreasing metallicity (Wild et al. 2011;Xiao et al. 2012).Similarly, Zahid et al. (2017) found distinct scaling relations of dust attenuation with SFR at different stellar masses.In the high-mass regime (M * > 10 10.2 M⊙), dust attenuation increases with SFR, while in the low-mass regime (M * < 10 10.2 M⊙), it decreases with SFR.Whitaker et al. (2014) found that the relation between galaxy stellar mass and attenuation evolves moderately at M * > 10 10.5 M⊙ over 0 < z < 2.5, but displays no or little evolution in the lowmass regime over the examined redshift range.Several studies using different indicators of attenuation have reported this puzzling 'no evolution' phenomenon in the stellar mass-dust attenuation relation for SFGs (Price et al. 2014;Whitaker et al. 2017;McLure et al. 2018;Shapley et al. 2022Shapley et al. , 2023;;Zhang et al. 2023).These correlations between dust attenuation and galaxy properties in SFGs demonstrate that dust attenuation is a complex process with multiple factors at play.
In the study by Qin et al. (2019a, hereafter referred to as Q19), we utilized a large sample of local SFGs to investigate the correlation between IRX and other galaxy parameters.By separating these parameters and determining independent correlations, we obtained several key findings: (1) Our analysis contradicted the previous understanding that galaxy stellar mass is the primary factor influencing IRX.We discovered that dust obscuration is not correlated with galaxy stellar mass, once controlling for other, more fundamental drivers.(2) We found that dust attenuation, as indicated by IRX, is determined by a combination of various parameters including SFR (or the indicator of infrared luminosity LIR), half-light radius (Re), gas-phase metallicity (Z), and galaxy inclination (axial ratio b/a).The relationship between these parameters and IRX can be described by a power-law function: IRX = 10 α L β IR R −γ e (b/a) −δ with the power-law slopes β, δ, and γ decreasing (i.e.getting closer to zero) as the metallicity decreases.
Furthermore, we verified that the empirical relation of dust attenuation obtained from local SFGs also holds out to z ∼ 2, demonstrating its universality.This universal IRX relation provides insights into the evolution of galaxies in terms of star formation, chemical enrichment, and morphological structure.However, the mechanism underlying how these three physical processes interact to make galaxies follow this relation remains to be explored.The determination of IRX is influenced by the dust content, star formation activity, and the geometry of dust and stars within a galaxy.It is now evident that the geometry of dust and stars is crucial for understanding the differences between various attenuation measurements and for quantifying the relative contributions that determine the overall dust attenuation of a galaxy.
In the past three decades, much effort has been devoted to studying galaxy dust attenuation using star-dust geometry models (Calzetti, Kinney, & Storchi-Bergmann 1994).Two commonly considered geometries are a uniform foreground dust screen and a homogeneous mixture of dust and emitting sources.However, in most cases, these simple geometries fail to adequately describe the observed attenuation patterns in real galaxies.Li et al. (2019) investigated the spatially resolved Balmer decrement distribution and its relation to other local and global properties, and concluded that neither the foreground screen nor the homogeneous mixture geometry can accurately represent the observed attenuation patterns (see also Kreckel et al. 2013).In addition, Q19 pointed out that while the simple homogeneous geometry can reproduce the IRX relations at Solar metallicity, it fails to reproduce the flat slopes observed at lower metallicities.This suggests that there is no universal star-dust geometry that applies to all types of galaxies.Instead, a more complex and flexible star-dust geometry model is required to better match the observations.
One approach to achieving this is by allowing the relative sizes of the dust and stellar discs to freely vary in the geometry model (Xiao et al. 2012).If the stars are more centrally concentrated than the dust disc, they will suffer an enhanced attenuation since the optical depth of the dust disc decreases with galactocentric radius.Another way to address the complexity of star-dust geometry is by considering a 'twocomponent' dust geometry, where the ISM of galaxies consists of dense birth clouds (BCs; associated with short-lived stars) and the surrounding diffuse ISM (Charlot & Fall 2000).This two-component geometry successfully explains the difference in attenuation between nebulae (birth clouds) and old stars in galaxies (Charlot & Fall 2000;Wild et al. 2011;Price et al. 2014;Koyama et al. 2019;Shivaei et al. 2020).The BCs in a galaxy follow certain distributions of ages and initial masses.If these distributions are similar for each galaxy, it would be expected that the average attenuation of BCs does not depend on either SFR or inclination.On the other hand, IRX measures the attenuation of the total UV light from young and intermediate-age populations.The former reside in the BCs while the latter are no longer surrounded by birth clouds due to stellar feedback.Then the steepness of IRX correlations might be influenced by the relative importance of dust attenuation caused by BC and diffuse dust.
In this work, we aim to build a new geometry model based on the two-component prescription with flexible model parameters, to reproduce the observational trends of IRX correlations highlighted in Q19.We focus on the scaling relations of dust attenuation, in particular, on the systematic trends of the steepness of IRX scaling relations with galaxy metallicity.In Section 2, we outline the details of our model as well as the impact of model parameters on dust attenuation.In Section 3, we describe how to fit the observed data with our geometry model.We show best-fitting results in Section 4 and discuss the implications in Section 5. Finally, we summarize our main results in Section 6.

THE STAR-DUST GEOMETRY MODEL
2.1 The distribution of dust and stars, and dust opacity While star formation does occur in various types of galaxies, including irregular galaxies and elliptical galaxies, the highest rates of ongoing star formation are typically observed in the most abundant disc galaxies which provide better conditions for the birth of stars.For example, ellipticals were found to contribute only 10-13% of the total ongoing star-formation budget (Zheng et al. 2009;Wuyts et al. 2011;Nelson et al. 2013;van der Wel et al. 2014b;Kaviraj 2014;Lofthouse et al. 2017).We thus focus on disc SFGs and address how dust attenuation is affected by different quantities/processes.We aim to examine the scaling relations of dust attenuation among galaxy populations.It is reasonable to assume that stars and dust are well mixed in the ISM, and the size distribution of dust grains and their chemical compositions are identical everywhere across a galaxy.
The ISM in a galaxy can be briefly described with two components: the dense birth clouds (BCs, or star-forming regions) with short-lived young stars, as well as more extended large-scale distributions of diffuse ISM, dust and old stars (Charlot & Fall 2000;Wild et al. 2011).The large-scale distribution of diffuse dust plays a major role in mediating the propagation of photons in galaxy discs and, at least in nearby galaxies, dominates the bolometric output of dust emission.To simplify, the total dust optical depth can be considered as the sum of the optical depth caused by the diffuse dust (diff ) plus the dense BC dust (dense): (1) This calculation relies on the assumption that the BCs and young stars share the same geometry distribution and no overlapping effects between BCs along the line of sight are taken into account in our model.A double exponential form is commonly used in astrophysics to capture the observed distribution of matter in galactic discs (Misiriotis et al. 2000;Smith et al. 2015).For the diffuse components, the distribution of diffuse dust ρ d and stars ρ⋆, both residing in a disc configuration, can be described by a double exponential, respectively where r and h are the radius and height in cylindrical coordinates, while R and H are the scale-length and scale-height of the dust or stellar disc, respectively.We use subscript ⋆ to represent the stellar disc and use subscript d to represent the diffuse dust disc hereafter.We note that the subscript d always refers to the diffuse dust disc.To not be confused, we add superscript tot if denoting the total dust disc (including also BCs).Then the integrated diffuse dust mass is given as Similarly, the integrated intrinsic luminosity of stellar emission is The observed luminosity of stellar emission due to the dust attenuation is where τ rh is the diffuse optical depth along the line of sight towards a given (r, z).For a galaxy observed under a face-on orientation, where Σ d is the diffuse dust surface density defined as M d /(πR 2 d ), and κ is the dust mass extinction coefficient that converts the dust surface/column density into dust optical depth.κ is wavelength dependent because different wavelengths of light have different absorption and scattering cross sections for the same dust grain, e.g.UV light has higher κ than optical.Here we omit the subscript λ to avoid the complexity of the notation hereafter.By definition, the effective optical depth caused by the diffuse dust is where We then do the following variable substitutions: Equation 8 can be rewritten as where The diffuse dust attenuation τ diff in our geometry model is only determined by the diffuse dust surface density Σ d , star-to-diffuse dust scale-length ratio R, star-to-diffuse dust scale-height ratio Ĥ, and the dust mass extinction coefficient κ.We notice that the thickness of both the dust and stellar discs (i.e., H d /R d and H⋆/R⋆) have vanished.
The above equations allow translating a given diffuse dust mass and geometry to the effective optical depth seen by a face-on observer.Under an inclined angle, sightlines to individual disc regions will cross dust located at different radii, but as we will demonstrate in Appendix A, the associated net attenuation of the galaxy light can still be approximated adequately using Equations 8-10, simply by boosting Σ d by a factor 1/(b/a).
For the dense (or BC) component, we first assume that the BCs are spherical and have similar sizes.We return to this assumption in Section 5.1.Additionally, for two similar size BCs, more metal-rich BCs with higher dust-to-gas ratio are expected to have higher BC optical depth (τ bc ).The galaxy dust attenuation probed by IRX measures the attenuation of the total UV light from young and intermediate-age populations (up to a few ×10 8 yr).The young stars are surrounded by BCs while the intermediate-age stars no longer live in BCs (due to feedback) but still contribute the UV emission (Kennicutt & Evans 2012).Given the evolutionary picture of BC dispersal, it is natural to introduce a BC covering factor C bc , defined as the total fraction of UV emission that is subject to BC attenuation.We note that our covering factor C bc is not the same as the clumpiness factor in the literature (e.g.Tuffs et al. 2004;van der Giessen et al. 2022), they assume the UV light is totally obscured in BCs (i.e., τ bc ≫ 1) whereas we do not.The optical depth caused by the dense BCs follows where I⋆,int denotes the UV light intrinsically emitted by the stellar population.Another important parameter in our two-component model is the BC dust mass fraction F bc , defined as the BC-to-total dust mass fraction.It controls the proportion of each component in the two-component model, and hence the evolution of dust geometry.For example, at F bc = 0 and F bc = 1, it reverts back to the pure BC-dominated and diffuse dustdominated geometry, respectively.Naturally, F bc also impacts the relative contribution of dust attenuation caused by BC and diffuse dust.Using the BC mass fraction F bc , we are able to link the diffuse dust surface density Σ d in Equation 10with the observationally more relevant total dust surface density Σ tot d .The latter is defined as , where M d is the total dust mass and R tot d is the scale-length of the total dust disc.Then we have where Rtot = R tot d /R d .In Appendix B, we explain how the size ratio of the total to diffuse dust disc can be written as a function of R and F bc .Then the total dust optical depth can be obtained with In summary, the effective dust optical depth in our twocomponent model depends on the dust extinction coefficient κ λ , total dust surface density Σ tot d , BC-to-total dust mass fraction F bc , star-to-diffuse dust disc scale length ratio R and scale height ratio Ĥ, BC optical depth τ bc , and BC covering fraction C bc .Other parameters have vanished during the derivation of the above formula.
Although the derivation of the above equations is based on young/intermediate-age stellar discs, it also holds for stellar discs at different ages (or a nebular disc) with appropriate changes to the model parameters.The old stellar disc has longer wavelength emission than young stars and therefore has smaller dust opacity (smaller κ and τ bc ).On the other hand, the old stars are not expected to be surrounded by BCs, and their BC coverage fraction C bc should be set to 0. Conversely, the nebular emission traces < 10 7 yr stellar emission, thus younger than those traced by the UV (or IRX), and is expected to have a higher C bc .
While in this paper, we focus on the IRX diagnostic of attenuation, the framework outlined in this section will thus also be applicable to the interpretation of, for example, Balmer decrement measurements, and how they contrast to IRX (see e.g.Qin et al. 2019b).The dashed lines represent the dust attenuation caused by diffuse dust, while the solid lines consider both diffuse and BC dust attenuation (with a fixed τ bc,V = 0.5).All model parameters (F bc , C bc , R, Ĥ) are fixed to a typical value of 0.5 except for the one indicated in the diagram.

Dependence of IRX on model parameters
The equations derived in Section 2.1 give the optical depth τ tot of a galaxy, while the attenuation indicator of infrared excess (IRX≡ LIR/LUV), defined as the ratio between infrared and UV luminosity, is used in our observational analysis.The intrinsic UV emission is dominated by the young (< 10 Myr) and intermediate-age (10-500 Myr) stellar populations.For a constant star formation history (SFH), the young stellar populations contribute 56% of the total UV luminosity (1216 Å-3000 Å), 19% from these at ages of 10 − 30 Myr, and the remaining 25% stems mainly from stellar populations of ages 30-500 Myr.Note that dust heating by old stars can contribute an increasing fraction to the far-IR luminosity (i.e., cold dust emission) at decreasing specific SFR but the fraction is marginal for normal SFGs (Nersesian et al. 2019).Moreover, the total IR luminosity of local SFGs used in our analysis is derived from the mid-IR luminosity (see Section 3.1) in which the contribution from the dust heated by old stars is negligible.We thus assign the total IR luminosity to the dust-reprocessed radiation of young and intermediateage stellar populations, and take IRX to reflect their dust attenuation.
A conversion between IRX and FUV optical depth should be made.Following Hao et al. (2011), the IRX can be obtained by IRX = (e τ FUV − 1) /0.46,where τFUV is the dust optical depth at the FUV according to Equation 13.This conversion is based on the energy balance principle in the sense that the absorbed UV-to-optical radiation is totally transformed into the IR through dust thermal emission.Throughout our analysis, we adopt a fixed κV of 0.62×10 −5 M −1 ⊙ kpc 2 following Kreckel et al. (2013) (see Section 3.4).We further adopt the Calzetti et al. (2000) attenuation curve, which has τFUV/τV = 2.5.Here the dust optical depths are associated with the projected dust surface density and thus dependent on viewing angle.Since our IRX estimates are based on the computed τFUV, they too will tightly correlate with the projected dust column density Σ dust,proj (see also Appendix A).We caution that in detail, the application of the energy balance principle is only warranted when comparing the UVto-optical emission versus IR emission as integrated over 4π steradian.For individual viewing angles, deviations from energy balance may occur as the IR radiation tends to be emitted isotropically, whereas the emerging UV light is not.The latter effect is not captured in our modelling.For disc SFGs, the average attenuation over 4π steradian equals that of an inclined SFG with project axial ratio b/a = 0.6 (see Q19).Accounting for the fact that the IR emission is anticipated to emerge isotropically would modestly increase (decrease) IRX for galaxies with higher (lower) b/a than 0.6 compared to the estimates obtained with the approach adopted in this work.Variations in the effective attenuation law, that in practice may arise from star-dust geometry and viewing angle effects, are also not captured.For instance, for the fixed κV adopted in this paper, a greyer (i.e., flatter) attenuation law would give rise to a lower τFUV and thus reduced IRX.We return to this point in Section 5.4.Fig. 1 shows the model-predicted IRX as a function of Σ dust for a set of model parameters.If τ bc = 0, i.e., in the case of diffuse dust only, IRX increases monotonically with dust surface density.This is consistent with the prediction of a homogeneous mixture star-dust geometry given by Q19, although the quantitative form of the relation shows modest differences due to the non-uniform density distribution of dust and stars (see also figure 1 in Zhang et al. 2023).Once BC attenuation is considered (where for illustrative purposes we here adopt a V -band BC optical depth of τ bc,V = 0.5), the Σ dust -IRX relation flattens at the low-Σ dust end.This is because the BC attenuation is assumed to be constant and independent of the opacity of diffuse dust.From Panel (b) we can see that the attenuation caused by cloud dust depends not only on the adopted BC opacity but also on the BC covering fraction C bc .The latter controls the fraction of UV emission that arises from within BCs.The more UV stars are located in BCs, the more dust attenuation.
The attenuation caused by diffuse dust is determined by F bc , R, and Ĥ but independent of C bc at a given total dust surface density Σ dust .F bc controls the remaining fraction of diffuse dust.The higher F bc , the less diffuse dust and hence the lower the overall dust attenuation.R and Ĥ control the relative spatial distribution between UV-emitting stars and dust in the galaxy.With decreasing R and Ĥ, the UV emission becomes more centrally concentrated within the dust disc, both radially and vertically.For an exponential dust disc, the central region always has a higher dust column density.We note that R and Ĥ have a very similar impact on the resulting Σ dust -IRX relation, and consequently might not be very distinguishable in subsequent model fitting.We will impose R = Ĥ in our model as discussed in Section 3.4.
In conclusion, in the BC dominated regime, IRX increases with τ bc and C bc , and has a flat Σ dust -IRX relation; in the diffuse dust dominated regime, IRX decreases with F bc , R and Ĥ, and has a steep Σ dust -IRX relation.The slope of the Σ dust -IRX relation is thus not constant over the full range of dust surface densities, but is controlled by the relative importance between BC and diffuse dust attenuation.We notice that the (projected) dust surface density is expected to correlate with the metallicity, SFR, Re and b/a (Wuyts et al. 2011;Li et al. 2019;Shapley et al. 2022).This gives a hint that the systematic trends between the power-law slope of IRX as a function of LIR or SFR, Re and b/a with gas-phase metallicity as shown in Q19 may arise from systematic changes in the relative importance of BC attenuation.

The local SFGs and universal IRX relation
We carry out our investigation of dust attenuation using the sample and data from Q19, which consists of 32,354 local SFGs.More details about the sample selection and data extraction can be found in Q19.We note that in Q19 the IR luminosities (8-1000 µm) were estimated from the single WISE W4 band.However, we find that the estimated total LIR from the combination of W3 and W4 bands have a smaller dispersion than using a single W4 band (see Appendix C).In this work, we therefore update the total IR luminosity by combining the WISE W3 and W4 bands.The IR-to-UV luminosity ratio is referred to as IRX and SFR is also estimated from the combination of IR and UV luminosities following Bell et al. (2005).
The original universal IRX relation states that IRX is determined by a combination of the infrared luminosity LIR, half-light radius (Re), gas-phase metallicity (Z) and galaxy inclination (axial ratio b/a).In this work, we replace the observational parameter of IR luminosity by the physically more relevant parameter of star formation rate (SFR).The two parameters are tightly correlated with each other (see figure 2   Here, log(Z/Z⊙) = 12 + log (O/H) − 8.69 represents the oxygen abundance and Z⊙ refers to Solar metallicity.The powerlaw slopes decrease with decreasing metallicity.Moreover, galaxies at z ∼ 2 also follow this empirical relation for dust obscuration, supporting its universality.
3.2 Determining projected Σ dust from metallicity, SFR, Re and b/a We have shown that IRX in the model is controlled by the dust surface density and other model parameters.However, our empirical universal IRX relation shows that IRX depends on metallicity, LIR or SFR, size and axial ratio.We infer dust surface density from Z, SFR, Re, and b/a in the following steps.At first, dust surface density (Σ dust = M dust /(πR2 e )) can be inferred from gas surface density via the dust-to-gas ratio (DGR) as The gas surface density is found to be tightly correlated with SFR surface density, known as the Kennicutt-Schmidt Law where n ∼1.4 as suggested by Kennicutt (1998).Here we do not use a fixed value for K-S slope but let it vary freely in our model.We notice that the 'surface density' presented here is defined using a half-light radius, 2 while the Equations in Section 2.1 use disc scale-length instead.The two radii relate as Re = 1.678R * for an exponential disc, and when applying Equation 10 we account for this difference in definition of surface density.Motivated by the power-law relation between stellar and dust radii discussed by ?, we parametrize On the other hand, the dust-to-gas ratio is found to be correlated with gas-phase metallicity in the form of a powerlaw (Issa, where µ ≡ ϵ/[(πϕψ 2 ) 1/n ] can be seen as a new normalization Power-law slope of DGR ∝ Z q (also τ bc ∝ Z q ) [0, 3] 0.98 parameter since these sub-parameters are degenerate and indistinguishable.Finally, we interpret the axial ratio as a projection effect in determining projected dust surface density, i.e., Here we do not distinguish between the axial ratios of stellar and dust discs, as the difference between the two is not significant.We notice that there is a possible risk to interpreting the axial ratio as a simple projection effect if the galaxy is not an ideal thin disc.We will show that it will not significantly affect our results in Appendix A. Finally, for brevity, we omit the 'proj' superscript hereafter.
In summary, the (projected) Σ dust of galaxies in our observed sample are calculated from four observed parameters: Z, SFR, Re and b/a, that used to define the universal dust attenuation relation in Q19.This calculation involves four free parameters: µ, p, n and q.
3.3 Metallicity dependence of τ bc , F bc , and C bc Metallicity is the most important parameter in our attenuation model, as it is believed to affect not only the dust surface density but also the star-dust geometry.At first, the BCs in this model are assumed to be identical, i.e., have the same radii and gas masses.Their dust properties are determined by the dust-to-gas ratio of the BC, which is found to be tightly correlated with gas-phase metallicity.Following Equation 18, the optical depth τ bc can be written as On the other hand, how to set the forms of F bc and C bc as a function of metallicity is a non-trivial task because there are no direct observational constraints to guide us.Some hints can be obtained from the scaling relations of other observed quantities.The specific SFR (sSFR= SF R/M * ) is found to be anti-correlated with gas-phase metallicity in the form of a power law (e.g.Lara-Lopez et al. 2013; Calabrò et al. 2017).
Generally speaking, SFR traces recent star formation in dense ISM (associated with BC) while the M * traces old stars in the diffuse ISM (associated with diffuse dust).Therefore, we can assume the BC-to-diffuse dust ratio follows a power-law and by definition, the BC-to-total mass ratio then follows to be At Solar metallicity, F bc,Z ⊙ = 1/(1 + 1/ζ).We rewrite Equation 23 by replacing the ζ with F Z ⊙ bc as On the other hand, C bc describes the fraction of UV light emitted from within BCs (relative to the total UV emission).
Similarly, we have C bc,Z ⊙ = 1/(1 + 1/ξ) at Solar metallicity and thus have obtaining dust surface density from the observed metallicity, SFR, Re, and b/a, we have the model parameters µ, n, p, and q; for obtaining IRX from Σ dust , we have parameters κ, τ bc,Z ⊙ , q, C bc,Z ⊙ , η, F bc,Z ⊙ , ν, R, and Ĥ.The definition, units and ranges of these model parameters are summarized in Table 1.
Of them, κV is fixed to 0.62 × 10 −5 M −1 ⊙ kpc 2 according to Kreckel et al. (2013).Kreckel et al. (2013) convert dust mass surface density to AV by assuming the observed Milky Way ratio of visual extinction to hydrogen column, and a fixed dust-to-gas mass ratio.They obtain Calzetti et al. (2000) attenuation curve is adopted with τFUV/τV = 2.5.
The normalization parameter µ in determining Σ dust from metallicity, SFR, Re and b/a can be constrained in observations.The DustPedia project uses the Herschel Space Telescope to make far-infrared high-resolution observations of 810 nearby galaxies, which is an ideal sample for our Σ dust calibrations from observations.By definition in Equation 19, the normalization parameter µ is defined as the Σ dust at SF R = 1 M⊙ yr −1 , Re,star = 3 kpc, and Z = Z⊙.We select a sub-sample with a limited range of Re around 3 kpc (from 2 to 4 kpc) and 12+log(O/H) around 8.69 (from 8.6 to 8.85) and fit the sample to obtain the best-fitting power-law relation between Σ dust and SFR.At SF R = 1 M⊙ yr −1 , we finally obtain the normalization parameter µ ≈ 10 5.16 M⊙ kpc −2 .This value is fixed in our model fitting.
We have two structural parameters in the model: the starto-diffuse dust scale-length ratio R and scale height ratio Ĥ.By analysing the dust profile of edge-on galaxy discs, ?found that both the 500-to-100 µm scale-height ratio (H500/H100) and the scale-length ratio (R500/R100) are on average larger than unity.Given that the 100 µm and 500 µm emission roughly trace SFR and dust mass, respectively, it suggests that the galaxy star-to-diffuse dust scale-height ratio Ĥ < 1 and the star-to-diffuse dust scale-length ratio R < 1 For simplicity, we let Ĥ ≡ R in our model fitting.
Then, our model has the remaining nine free parameters to be constrained: n, p, q, τ bc,Z ⊙ , C bc,Z ⊙ , η, F bc,Z ⊙ , ν, and R. We make use of the Monte Carlo Markov chain (MCMC) sampling method to explore the 9-dimensional parameter space and efficiently characterize the uncertainties associated with the derived parameters.In detail, we consider a Gaussian likelihood function and assume flat priors for our model parameters in the following ranges: With the large sample of our local SFGs, we can place considerable constraints on the geometry model with these free parameters.The Python package emcee4 (Foreman-Mackey et al. 2013) is used to perform the MCMC sampling using an affine invariant algorithm proposed by Goodman & Weare (2010).The algorithm behind emcee has several advantages over traditional MCMC sampling methods and it has excellent performance as measured by the autocorrelation time.The sampler is initialised around the position in parameter space that gives the best fit to the data with a maximum likelihood method.
As shown in Fig. 2, we obtain an excellent match to the observed IRX.The 1 σ dispersion of the IRX residual is ∼ 0.19 dex, consistent with the dispersion of empirical powerlaw fitting applied by Q19.Moreover, this relation derived from local SFGs also holds for galaxies up to z = 3. Fig. 2 also shows the corner plot of the MCMC fit in order to examine the robustness of the best-fitting parameters.We can see that, although there are some degeneracies (e.g. a lower C bc tends to couple with a higher τ bc in the fitting and vice versa), the model parameters are well constrained.Even in the worst case (e.g. an anti-correlation between inferred C bc,Z ⊙ and τ bc,Z ⊙ ), the relative error does not exceed 30%.We point out that our sample includes 37 local SFGs with Z <∼ 1/3 Z⊙ to provide reasonable constraints on these BC model parameters, although the vast majority of our sample galaxies cover the high-metallicity regime where these parameters are hardly constrained.The estimated values of each parameter are listed in Table 1.

RESULTS
We show in this section the results of fitting the IRX of local SFGs with our geometry model.We first consider the best-fitting power-law indices of the ISM scaling relations in Section 4.1 and the best-fitting geometry parameters as a function of metallicity in Section 4.2.Then, in Section 4.3, we show how our geometry model reproduces the universal IRX relation presented in Q19.
4.1 The best-fitting scaling relations in determining Σ dust We have three free parameters n, p, and q in determining Σ dust based on the galaxy scaling relations.The best-fitting KS-Law power-slope n is 1.84, which is slightly steeper than the value of 1.4 by Kennicutt (1998) or the revised version of 1.5 by Kennicutt & De Los Reyes (2021).However, the precise KS-Law slope depends on the sample selection and the method used to measure the physical quantities.For example, Kennicutt & De Los Reyes (2021) found that the slope increases to 1.92 if a metallicity dependent XCO (following the prescription of Bolatto, Wolfire, & Leroy 2013) is adopted, consistent with our best-fitting result.Using a somewhat different but physically motivated prescription of XCO, Narayanan et al. ( 2012) also found a similarly steep slope of 1.95.In a modelling of observed Balmer decrements using simplified dust geometries, sharing several ingredients in common with our approach, Li et al. ( 2019) likewise concluded that a significantly superlinear KS-slope (of ∼ 1.5) was required to reproduce the observed Hα/Hβ ratios.
The best-fitting relation between dust and stellar (r-band) radius follows R e,dust ∝ R p e,star with the power slope p = 0.69.Using the carefully selected sample from the DustPedia project, the best-fitting relation between dust and stellar radius follows R e,dust ∝ R 0.8 e,star , very close to our model prediction.This indicates that the growth of the dust disc is slower than the stellar disc.In an inside-out disc growth scenario, the materials that build a stellar disc are primarily accreted from the surrounding environment.In this sense, the properties of a galaxy would be to a large extent shaped by the properties of the accreted material (such as specific angular momentum, accretion rate, etc.).Pan et al. (2021) found that more gas-rich galaxies have a higher HI-to-stellar-size ratio.Our low-mass, gas-rich galaxies have a higher dust-to-stellar size ratio.
A near linear Z-DGR relation (i.e., q ∼ 1) is suggested by our model fitting.This is consistent with the observational results that the DGR is well represented by a power law with a slope of ∼1 (at least at 12+log(O/H)>8.2, e.g.James et al. 2002;Draine & Li 2007;Leroy et al. 2011;Sandstrom et al. 2013;Giannetti et al. 2017;Galliano, Galametz, & Jones 2018), and with theoretical expectations due to the dust grain growth (e.g.Asano et al. 2013;Zhukovska 2014;Feldmann 2015;Aoyama et al. 2017;McKinnon et al. 2018) Different colours denote the optical depth of a single BC in the V -band τ bc,V (blue), the BC-to-total mass fraction F bc (black), the fraction of UV light emitted from within BCs C bc (red), the star-to-total dust scale-length ratio Rtot (green).Here Rtot differs from the star-to-diffuse dust disc scale-length ratio R which is metallicity-independent in our model.The vertical dashed line marks Solar metallicity.
ing on the adopted metallicity calibration) is sufficient to describe the observed Z-DGR relation at all metallicities.This is based on the chi-square (χ 2 ) of a linear regression to the sample and does not necessarily imply that the true relationship follows this functional form.If we only focus on the high metallicity regime, as shown in figure 10 of their paper, the dust-to-metal ratio more or less flattens with metallicity, indicating that there is a linear Z-DGR relation in the high metallicity regime.

The best-fitting geometry parameters
For the geometry parameters, including F bc , τ bc,Z ⊙ , C bc , and Rtot , we show them as a function of metallicity in Fig. 3.We find a birth-cloud optical depth of τ bc,V = 0.31 at Solar metallicity, which is comparable to the value of 0.33 at same metallicity as given in Lu et al. (2023).τ bc,V increases with metallicity following τ bc ∝ Z, i.e., a linear relation.This comes from the best-fitting linear Z-DGR relation, and the τ bc in our model is only determined by the dust-to-gas ratio.
The BC mass fraction F bc increases with decreasing metallicity.It equals 0.14 at Solar metallicity and increases to 0.53 at 1/4 Solar metallicity.This is consistent with the lowermass/metal-poor galaxies that have younger stellar populations and more clumpy geometry (Baes et al. 2020).By modelling the inclination dependence of stellar and nebular attenuation together with a dust geometry model, Lu et al. (2022) suggested that the typical BC mass fraction is 0.3 for an MW-like galaxy.Determining F bc in observations is a non-trivial task, as it is not a directly observed quantity.The observed dust SED can be divided into two parts: the hot dust around young stars and the diffusely distributed cold dust (Draine & Li 2007), very similar to our definition of BC and diffuse dust, respectively.By decomposing the infrared SED of nearby galaxies, Draine et al. (2007) found that the mass fraction of hot dust is in the range of 0-6%.On the other hand, the molecular gas is often associated with BCs and the atomic gas is associated with the diffuse ISM.The molecular gas fraction (∼ F bc ) is ∼20% for our local galaxies (Catinella et al. 2018).We do not intend to compare our F bc to either hot dust fraction or molecular gas fraction, but note that they span a similar range as our best-fitting F bc .
The BC covering fraction also decreases with metallicity, with values close to unity at lower metallicity but dropping to 0.84 at Solar metallicity.Massive and metal-rich galaxies are relatively older and have a higher fraction of intermediate-age stars that still contribute to the UV emission but had their surrounding BCs dispersed due to stellar feedback.Similar to F bc , there are no direct observational constraints on C bc .The dense clouds at low metallicity not only have a higher fraction of their dust comprised in BCs, but also a higher fraction of UV stars embedded in BCs.In contrast, if there is no BC dust (i.e., F bc = 0), there are no UV stars surrounded by BCs (i.e., C bc = 0).
The radiation and stellar winds from the newly-born massive stars may clean up the leftover gas in BCs over a timescale of only a few Myr, before the occurrence of supernova (SN) explosions from these stars (e.g.Kruijssen et al. 2019).It is reasonable to expect that intermediate-age stellar populations could exist in metal-poor galaxies although the best-fitting C bc ∼ 1 at low metallicity suggest the contribution from intermediate-age populations may be negligible.Rather, star formation across a galaxy can reasonably be seen as a constant process over the most recent 500 Myr, possibly independent of the metallicity of the ISM in the galaxy and thus giving metallicity dependence of ν ∼ 0, instead of the obtained best-fitting ν = −3.9, which is quite steep.We argue that the best-fitting C bc and ν at low metallicity should not be thought of a measure of the UV radiation fraction from within BCs; instead it may encode physics that impacts what fraction of UV light will be subject to BC attenuation.Alternatively, the expelling of gas and dust driven by stellar feedback might become inefficient in metal-poor BCs and the expelled clouds could remain clumpy and act similarly to BCs in attenuating the UV light.We notice that at low metallicity these model parameters might be biased by our assumption of R = Ĥ (see Section 5.4 for more details).
The best-fit star-to-diffuse dust radius ratio R is ∼0.53 in our model, which suggests that the stellar disc has a smaller radius than diffuse dust.The star formation often takes place in more central regions that feature a higher dust/gas density, while the diffuse dust extends further into the galaxy outskirts (Kennicutt & Evans 2012;Smith et al. 2016;Casasola et al. 2017).We notice that the observed dust emission consists of both dense/BC and diffuse dust.To compare with observations, we calculate the star-to-total dust radius ratio Rtot based on R and F bc (see Appendix B).We can see that the star-to-total dust radius Rtot decreases with decreasing metallicity.In the case of extremely metal-poor environments (below the metallicity coverage of our sample), F bc is expected to be close to unity and the 'total dust' is dominated by BC dust, which traces the stars.In other words, BC dust and total dust discs are the same one, i.e, Rtot = 1.Conversely, in the high-metallicity limit, F bc approaches zero, diffuse dust plays a dominant role and therefore Rtot = R = 0.52.By examining the radial distribution of dust, stars, gas, and SFR in a sub-sample of 18 face-on spiral galaxies extracted from the DustPedia sample, Casasola et al. (2017) found the average Rtot = RSFR/R tot dust is ∼0.57(see their table 7).We found the average metallicity of their sample is 8.6 and a consistent Rtot of 0.53 is predicted in our model at that metallicity. 5We provide a PHTHON package IRX TAU TOT to calculate the total dust optical depth of a galaxy with given metallicity and best-fitting geometry parameters.6

Comparing with the Universal IRX Relation presented in Q19
The main goal of this paper is to understand the physical origins of the universal IRX relation presented in Q19, particularly the systematic flattening of power slopes towards lower metallicity.Given a set of best-fitting model parameters, we have IRX (Z, SFR, Re, b/a).Then the IRX-Z relation can be predicted with our model at given SFR, Re and b/a.As shown in Q19, besides the metallicity term, the power-law indices of SFR, Re and b/a in the power-law relations also contain a dependence on metallicity.One way to eliminate the metallicity dependence of other parameters and obtain a 'pure' IRX-Z relation is by normalising SFR, Re and b/a to unity.For consistency, we apply the same normalisation of SFR,Re and b/a to unity in our model.We compare our thus obtained model-predicted IRX-Z relation against the relation given by Q19 (who similarly normalised SFR, Re and b/a to unity) in the first panel of Fig. 4. We can see that our model perfectly reproduces the IRX-Z relation given by Q19.A power-law relationship adequately describes the observed IRX-Z relation.
Similarly, the SFR-IRX relation can be obtained from the best-fitting model by fixing metallicity, Re and b/a.We notice that the slope of the empirical SFR-IRX relation is not constant but decreases with decreasing SFR.7 One way to reduce this to a single measure of slope is adopting the slope of the line tangent to the curve at some representative SFR.We divide our sample galaxies into continuous metallicity bins with a constant bin width of 0.1 dex.In each metallicity bin, we calculate the average metallicity, SFR, Re and b/a.With these averaged values, we obtain a model-predicted SFR-IRX relation and then measure the slope of the tangent line to the curve at the average SFR (i.e., β).Repeating the calculation to all metallicity bins, the β as a function of metallicity is obtained (top-right panel of Fig. 4).With a similar approach, the power slope of the Re-IRX relation (γ) and b/a-IRX relation (δ) as a function of metallicity can be also determined, respectively.As shown in the remaining panels of Fig. 4, all the derived power-law slopes β, γ, and δ from our best-fitting two-component geometry model agree perfectly with that in Q19 as a function of metallicity.We notice that the dispersion of IRX residuals decreases slightly from 0.188 dex in the empirical power-law fitting to 0.185 dex in the geometry model fitting.That is, the IRX relations are described slightly better by our geometry model than by the power-law forms.The latter should be regarded as an approximate description over a certain dynamic range of observed parameters.

Constant τ bc at a given metallicity
While the IRX caused by diffuse dust monotonously increases with Σ dust , adding a BC component with constant optical depth can be very effective in flattening the slope of IRX relations as shown in Fig. 1.In other words, the key point to capture the systematic change in IRX relation slopes is that we assume the optical depth of BCs is uniquely set by metallicity but is independent of other galaxy properties, such as SFR, Re and b/a.It is natural to question whether the average τ bc may depend on other galaxy properties than metallicity instead.
A galaxy consists of numerous (of the order of hundreds to millions) BCs, that follow certain distributions of ages and initial masses (determined by the star-forming process).If these distributions are universal, then the average optical depth of the BCs should be similar in each galaxy.In contrast, for example, if the BC mass function and stellar IMF depend on the global properties of galaxies, such as their SFR or gas surface density, then one might expect that the average dust opacity is correlated with such global properties.Lu et al. (2023) argued that it is reasonable to assume the optical depth of BCs (they call it 'clumps') is only decided by the metallicity or dust-to-gas ratio, since the average electron density and size distribution of H II regions (associated with BCs) are similar for all types of local SFGs (Oey et al. 2003;Liu et al. 2013;Santoro et al. 2022).However, Santoro et al. (2022) found that the power-law slope of the Hα luminosity function of extragalactic H II regions decreases with the galaxy star formation rate surface density ΣSFR.This would correspond to an enhanced clustering of young stars at high gas surface densities.If the star formation efficiency is similar among the BCs, the higher SFR surface density galaxies would then have more massive BCs on average.We note that more massive BCs do not necessarily have a higher column density or optical depth.The masses and sizes of giant molecular clouds (GMC) in our Milky Way are found to follow a relation of MGMC ∝ R 2 , i.e., a constant surface/column density (Larson 1981;Lombardi, Alves, & Lada 2010;Ballesteros-Paredes, D'Alessio, & Hartmann 2012;Chen et al. 2020;Lada & Dame 2020).Mapping CO emission at the cloud-scale resolution for ten nearby galaxies, Rosolowsky et al. (2021) likewise found the extragalactic GMCs to follow a mass-size relation of MGMC ∝ R 2 .Using the data of average ΣGMC of nearby galaxies given by Rosolowsky et al. (2021), we find that, within uncertainties, the average ΣGMC changes little with either global SFR or gas surface density.In brief, there is no evidence that the average τ bc of nearby galaxies should be correlated with the galaxy properties, particularly the SFR or gas surface density.It is a good approximation that the optical depth of BCs is only determined by metallicity without additional dependencies on other galaxy properties (see also Lu et al. 2023).Exceptions may occur in the more extreme environments of local ultra-luminous infrared galaxies (ULIRGs, Krumholz & McKee 2005) or higher-redshift galaxies which are substantially more gas-rich and surface dense in a galaxy-averaged sense, and also feature enhanced local densities (Davies et al. 2021).

Physical origins of the universal IRX relation
One main goal in this work is to explore the physical origins of the universal IRX relation, particularly the systematic changes in the slopes of IRX as functions of SFR, Re, and b/a with metallicity.To do so, we built up a new 'twocomponent' dust geometry model with flexible parameters and performed a model fitting to the observed data.Our model fitting shows consistent results with the power-law fitting presented in Q19, including the systematic changes in the power-law slope of IRX with metallicity, SFR, Re and b/a (see Fig. 4).We present a schematic diagram in Fig. 5 to illustrate our result of how the dust geometry evolves with metallicity.
At first, the BC attenuation is assumed to be constant (only determined by metallicity) while the diffuse dust attenuation increases with Σ dust (see bottom panels).The latter is determined by the combination of Z, SFR, Re, and b/a in the forms of power-laws.The changes in the slopes of the total Σ dust -IRX relation are controlled by the competition of the BC and diffuse dust attenuation, which is regulated by the changes in metallicity.At low metallicity, there is less diffuse dust and the BC attenuation dominates.In this case, one can expect that the total dust attenuation is only controlled by the gas-phase metallicity (changing dust-to-gas ratio) and independent of either SFR, disc size or inclination.This is why we see a flat Σ dust (SFR, Re, b/a)-IRX relation at lower metallicity in the bottom-left panel of Fig. 5.
At high metallicity, there is a large amount of diffuse dust, both due to an increasing total amount of dust and a higher diffuse dust fraction.A fraction of UV stars lose their surrounding BCs due to the stellar feedback (the empty circle), and the attenuation caused by the BCs will be somehow reduced.Although the BC dust opacity moderately increases due to the increased dust-to-gas ratio at higher metallicity, the diffuse dust attenuation still plays a dominant role.The total dust attenuation now becomes sensitive to the changes in global dust surface density that are decided by SFR, metallicity, size and inclination.This is why we see a steep Σ dust (SFR, Re, b/a)-IRX relation at high metallicity in the bottom-right panel of Fig. 5.
Unlike SFR, Re, and b/a, the variations in galaxy metallicity do not only change the galaxy-averaged (projected) dust surface density through the associated dust-to-gas ratio, but also alter the dust geometry in galaxies.Lowmetallicity galaxies have a lower dust surface density and a BC-dominated or clumpier dust geometry.The diffuse starlight embedded in a clumpy geometry will more easily escape and BC clouds themselves feature a lower opacity, altogether leading to a smaller net attenuation.In other words, the positive Z-IRX relation in observations is shaped by the combination of a positive Z-Σ dust relation as well as an anticorrelation between metallicity and clumpiness of the geometry.
While galaxy stellar mass is a crucial physical parameter governing the galaxy formation process, gas-phase metallicity directly measures the gaseous metal abundance and indirectly probes the dust abundance of the ISM, where star formation and dust obscuration are intimately connected.Therefore, a measure of dust attenuation such as IRX is more closely linked to the specific property of gas-phase metallicity than the absolute amount of stellar mass.As mentioned in Sec-tion 3.1, the universal IRX relation involves SFR, Re, Z/Z⊙, and b/a, but not M * .The dependence on b/a is solely attributed to the inclination effect, while SFR, Re, and Z/Z⊙ are all related to the ISM of SFGs.Here, Re is measured from the stellar light profile but is used as a proxy for the spatial distribution of the ISM and young stellar populations (see Q19 for more discussion).It is important to emphasize that the dust attenuation parameter IRX is derived from the fraction between the absorbed and unabsorbed bolometric luminosity of the young and intermediate-age stellar populations that are tightly associated with the ISM, rather than the older stars that dominate the stellar mass in SFGs.The geometry parameters in our model, such as the star-to-total dust scale-length ratio R tot , are set to be related to Re of a galaxy and thus are connected with stellar mass.Our model fitting results show that R tot moderately decreases with metallicity.Given that metallicity is correlated with stellar mass, the best-fitting R tot as a function of metallicity might also reflect its dependence on stellar mass.For general populations of disc-dominated SFGs, the ISM is globally associated with older stellar populations in terms of its geometry.In conclusion, our study shows that our model parameters as a function of metallicity are able to unveil the processes that regulate the universal IRX relation.

Explaining the no evolution of mass-attenuation relation
Recent studies found that the relation between dust attenuation and stellar mass evolves mildly with cosmic time up to z ∼ 3 (e.g.Price et al. 2014;Whitaker et al. 2017;McLure et al. 2018;Shapley et al. 2022Shapley et al. , 2023;;Zhang et al. 2023).However, it introduces a 'non-evolution puzzle' since either the SFR, ISM content, size, and metallicity evolve significantly with redshift at fixed stellar mass.Quantitatively, Shapley et al. (2022) found the dust surface density (estimated either via direct or indirect methods) to increase by a factor of ∼4 from z = 0 to z = 2.3 at given M * while the dust attenuation remains unchanged.This means that the dust attenuation per unit of Σ dust (hereafter 'dust attenuation efficiency') should decrease towards higher redshifts.They propose several potential origins at z ∼ 2.3: 1) a smaller dust-to-metal ratio (steeper Z-DGR relation); 2) a dust distribution that is relatively more extended compared to the stars; 3) clumpier dust distributions; and 4) a lower dust mass extinction coefficient κ λ .
Of them, 2) is firstly excluded since high-z galaxies usually have smaller far-IR sizes than their local counterparts (e.g.Fujimoto et al. 2017;Tadaki et al. 2020;Gómez-Guijarro et al. 2022).8In our geometry model, both the Z-DGR relation and κ are fixed, i.e., 1) and 4) are not included in our model.In other words, through the modelling study in this work, only the clumpier geometry at high-z is favoured.Highz galaxies usually have lower metallicity than local galaxies, and they are expected to have a clumpier geometry as shown in Fig. 5.In this case, more dust is in the form of BCs rather than smoothly distributed in the ISM.The decrease of diffuse dust causes a decrease in dust attenuation, while the increasing BC dust only increases the number of BCs but does not affect the global dust attenuation.On the other hand, if the BC covering fraction is not 100%, the more clumpy geometry means that a fraction of starlight will escape directly from the galaxy without encountering clumps and suffers little dust attenuation.Both effects lead to a smaller dust attenuation at a given dust surface density.The interpretation in terms of a more clumpy (i.e.BC-dominated) dusty geometry is also in line with the reduced (or even absent) inclination dependence of attenuation in high-z SFGs as observed out to cosmic noon by Lorenz et al. (2023) and Zhang et al. (2023), and beyond cosmic noon by Gómez-Guijarro et al. (2023).
Motivated by the universality of our IRX relation, we predict the mass-Σ dust relation and mass-attenuation relation at different redshifts based on the best-fitting model as well as the evolution of galaxy scaling relations, including the mass-SFR relation (Popesso et al. 2023), the mass-metallicity relation (Sanders et al. 2023), and the mass-size relation (van der Wel et al. 2014a).The results are shown in Fig. 6.We remind the reader that such predictions are based on the assumption that all of the K-S Law, Z-DGR relation, R d -R⋆ relation and metallicity-dependence of dust geometry changes only mildly across cosmic time.We can see that, although a significant evolution of Σ dust at given M * is present, the dust attenuation evolves mildly at low and intermediate masses.This indicates that a clumpy geometry at high-z is sufficient to explain the puzzle of significant evolution in Σ dust but no evolution in dust attenuation (see also Zhang et al. 2023).Besides, our model predicts a moderate evolution of dust attenuation at M * > 10 10.3 M⊙.This is consistent with Whitaker et al. (2014) who found the IRX evolves mildly at M * < 10 10.5 M⊙ but evolves moderately at higher masses.However, Shapley et al. (2022) did not find a moderate evolution of dust attenuation traced by the Balmer Decrement at the high-mass end.We notice that the sample used for Shapley et al. (2022) is not 'very' massive, and the majority of galaxies have M * less than ∼ 10 10.5 M⊙ (see also Shapley et al. 2023).Another possible reason might be attributed to the difference between the Balmer decrement and IRX indicators.The former traces nebular emission attenuation, while the latter traces the stellar emission attenuation (Qin et al. 2019b).More efforts should be made to give a self-consistent explanation.Such an analysis, however, goes beyond the scope of this work.
Given that both distant and nearby SFGs conform to the universal IRX relation, it is intriguing to explore the evolutionary trajectory of a galaxy from high redshift to the present day.Extracting such an evolutionary track from the IRX-mass-redshift relation depicted in the middle panel of Fig. 6 requires connecting data points representing low-mass SFGs at z = 3 with more massive SFGs at lower redshifts.Following the methodology outlined by Förster Schreiber & Wuyts (2020) to retrieve the evolutionary paths of a Milky-Way-like galaxy from z = 3 to the present day, we show in Fig. 7 their derived histories of star formation rate (SFR), stellar mass (M * ), effective radius (Re), and metallicity (Z/Z⊙).Additionally, for comparison, we present the corresponding best-fitting parameters (τ bc,V , F bc and C bc ) of our two-component star-dust geometry model, as well as the model-predicted IRX.It can be observed that, in general, the IRX curve initially follows the SFR curve, with deviations likely resulting from metallicity-and structure-related processes.A striking transition occurs roughly around z ∼ 1.5 from the metal-poor and compact progenitor to the metalrich and extended disky galaxy.

Caveats
For simplicity, we have utilized the Calzetti Law as the fixed galaxy dust attenuation curve in our modelling.Additionally, we have assumed that the star-to-diffuse dust disc scalelength and scale-height ratios follow each other as R = Ĥ.The Calzetti Law has been successfully applied in our modelling of the universal IRX relation for the overall population of SFGs.It is important to note that the dust attenuation curve may differ for various galaxy populations.For example, the Milky Way, Large Magellanic Cloud, and Small Magellanic Cloud exhibit different curves.Previous studies have reported variations in the slope of the dust attenuation curves for different galaxy samples, which can be either steeper than the Calzetti Law (e.g.Shivaei et al. 2020) or shallower (e.g.Salim & Narayanan 2020).The slope of the dust attenuation curves and dust attenuation were reported to be anti-correlated (Salim & Narayanan 2020), likely influenced by degeneracies in galaxy spectral energy distribution (SED) fitting (Qin et al. 2022).Furthermore, the inclination of galaxies has been found to affect the shape of the observed attenuation curve, with edge-on galaxies exhibiting a grayer attenuation curve compared to face-on galaxies (see e.g., Appendix of Zhang et al. 2023).These findings suggest that incorporating flexibility to the dust attenuation curve could further improve the realism with which observations are reproduced by our modelling, and could provide insights into the factors influencing the slope of the effective dust attenuation curve.
A greyer curve, as reported for highly-inclined (i.e., edgeon) SFGs (e.g.Lu et al. 2022), may result in a lower τFUV at a given IRX, as a relatively larger portion of the dust-absorbed energy that is reprocessed into the IR would originate from stellar emission at wavelengths longward of the FUV regime.In such cases, the use of the Calzetti Law in our modelling could lead to an overestimation of τFUV.However, we anticipate the model parameters best-fitting the observed IRX, metallicity, SFR, and Re to remain largely unchanged.We note that the highly-inclined galaxies represent only a small fraction (2.5% for SFGs with b/a < 0.2) of our sample of SFGs, leaving the results of our model fitting unaffected.We acknowledge that our model does not have the capability to determine the slope of the dust attenuation curve as enabled for example by the radiative transfer approaches employed by Zhang et al. (2023), but note that this also falls beyond the scope of this work.
Our assumption that R = Ĥ, meaning that the stellar disc and the dust disc mirror each other in shape with a scaling factor, ensures that the axial ratios of the stellar and dust components are identical.This assumption simplifies the calculation of inclination-related effects and is a crucial parameter for our two-component star-dust geometry model.Our results demonstrate that this assumption holds well on a global level.However, in real galaxies, the values of R and Ĥ can vary within a range and may not necessarily mirror each other on an individual basis.For example, massive disc galaxies like NGC 891 often exhibit a thick stellar disc and a thin dust disc, resulting in an apparent dust lane when viewed edge-on.Such variations in R and Ĥ can introduce scatter in the dust surface density (Σ dust ) relative to star formation and, consequently, in dust attenuation.For metalpoor galaxies, their morphologies tend to be spherical or rounder compared to regular discs, resulting in a higher Ĥ relative to R. Such star-dust geometries significantly reduce the inclination-dependent effect on dust attenuation (Q19).Metal-poor galaxies may thus require distinct model parameters R and Ĥ compared to those determined predominantly by the metal-rich galaxies in our sample.The fixed R = Ĥ in our model fitting may bias C bc towards unity in order to dramatically reduce the inclination-dependent effect at low metallicity.On the other hand, the dust attenuation indicator, IRX, is primarily influenced by the radiation from young and intermediate-age stellar populations, which are closely associated with the dust disc.IRX is less sensitive to the contribution from old stellar populations.It is important to note that a comprehensive analysis of dust attenuation should consider the geometric differences between stars and dust, accounting for possible variations in both R and Ĥ.
Dust attenuation in our model accounts for absorption of light and does not take into consideration the effect of scattering.Light scattering would introduce additional light for galaxies seen face-on and less for galaxies seen edge-on.We believe this effect plays a secondary role in affecting the relationship between IRX and effective dust surface density.In contrast, the orientation effect of increasing projected dust columns with increasing inclination, which can be accurately quantified, has a more significant influence on this relationship.

SUMMARY
In this work, we have developed a new two-component dust model to describe the universal dust attenuation relation in local galaxies.We fit our model to the same SFG data exploited in Q19 using a Bayesian MCMC sampling.Our main findings are summarized as follows: i) Our model produces a good fit to the observational data, and successfully reproduces the systematic changes in IRX scaling relations over the entire metallicity range as presented in Q19.ii) We give constraints on three galaxy scaling relations through model fitting, including the star-formation law of ΣSFR ∝ Σ 1.84 gas , the dust-stellar radius relation of R e,dust ∝ R 0.69 e,star , and the metallicity-dust/gas relation of DGR ∝ Z.All of these relations are consistent with the observational studies.
iii) The evolution of dust geometry as a function of gasphase metallicity is also constrained quantitatively.As metallicity increases from 1/3 Solar to Solar, the starto-total dust disc scale-length ratio decreases from 0.69 to 0.53, the optical depth of birth clouds increases from 0.11 to 0.34, the birth cloud mass fraction decreases from 0.42 to 0.14, and the BC covering fraction of UV light decreases from ∼1 to 0.75.Low-metallicity galaxies not only have a smaller dust surface density but also a clumpier geometry.iv) We find the variations in the slopes of IRX with SFR, Re and b/a stem from the competition between BC and diffuse dust attenuation in galaxies, which is controlled by the galaxy metallicity.When a galaxy is metal-rich, there is a large amount of diffuse dust in its ISM, and IRX increases with the global dust surface density; At low-metallicity, the birth cloud attenuation becomes important and even dominant, and then the IRX becomes insensitive to the changes of either SFR, Re or b/a.

APPENDIX A: THE INFLUENCE OF THE DISC THICKNESS
As mentioned in Section 3.2, the derivation of our geometry model equations is based on the assumption that galaxies are oriented face-on.We include the axial ratio to our model by interpreting it as a simple projection effect, i.e., an edge-on galaxy (axial ratio of b/a) with an observed (i.e., projected) dust surface density of Σ dust is in our model treated as a face-on galaxy with a dust surface density of Σ dust ÷ (b/a).This simplified approach ignores the fact that for an inclined observer the sightline to stars at any galactocentric radius will pierce through dust located at a range of galactocentric radii, an effect that will be more pronounced for thicker discs and/or closer to edge-on viewing angles.In other words, disc thickness may bias the projection effect.Does it affect our results?We start by considering the τ rh in Equation 7. At first, we assume a galaxy has a certain inclination θ with respect to the observer, defined such that θ = 0 • corresponds to a face-on viewing angle.We further treat the galaxy as an axi-symmetric system.Rewriting τ rh in Equation 7 as τ xyhθ , and introducing l as the integration variable capturing distance along the line of sight, we obtain where h ′ = h H d and Σ diff d refers to the galaxy-averaged surface density of the diffuse dust disc.The expression for the effective optical depth due to diffuse dust (Equation 8) should likewise be rewritten as We can see τ diff is a function of κ, Σ dust , R⋆, H⋆, R d , H d and θ.We notice that we mainly focus on the influence of relative thickness but do not care about the absolute value of the galaxy radius.Here we let R d = 5 kpc.We have checked that using another value of radius does not alter our conclusion.According to the our best-fitting model, we let R = Ĥ ≈ 0.5, i.e., R⋆ = 2.5 kpc, H⋆ = H d /R d × 2.5 kpc.The κ is fixed as 0.68 × 10 −5 M −1 ⊙ kpc 2 in the V -band following Section 3.4.Now, the optical depth of diffuse dust τ diff depends on the face-on Σ dust , thickness H d /R d , and inclination θ.We integrate the equations numerically to obtain the dust opacity at different inclinations.Besides, we project the three-dimensional galaxy disc onto a plane with a viewing angle of θ for a set of different thicknesses.We then use the python package photutils9 to fit the projected ellipse to obtain the axial ratio.
The top-left panel of Fig. A1 shows IRX as a function Σ dust for a disc-like galaxy with a thickness of 0.1.We can see that at a given Σ dust , more inclined galaxies have higher IRX.If the thickness increases to 1 (H d = R d , the top-right panel), the shape of the galaxy is close to a sphere, and therefore the dust opacity becomes no longer sensitive to the changes of inclination.The bottom panels show the IRX as a function of projected dust surface density, determined by dividing Σ dust by b/a.Surprisingly, at a given Σ dust,proj , the dependence of IRX on inclination disappears for different thicknesses.Although our derivation of the model formulas is based on the face-on galaxy orientation, including the axial ratio via a simple projection approach enables our results not to be affected.We here derive the relation between the size ratio of SFR and dust disc when accounting for all dust ( Rtot ) or only diffuse dust ( R).As we will demonstrate, the relation between the two varies with changing birth cloud dominance (F bc ).Again, we use no superscript, superscript bc, and superscript tot to denote the diffuse dust, BC dust and total dust components, respectively.We assume both BC and dust disc follow exponential profiles as and the integrated mass of BC and diffuse dust are F bc and 1 − F bc , respectively (i.e., for simplicity we use a total dust mass of unity in our derivation).Then the total dust (BC+diff) profile is For a given F bc and R, we thus have a profile for the total dust disc and the Rtot can be obtained by fitting the profile (within r ′ < 3).In other words, Rtot is a function of R and F bc .We show the numerically derived Rtot as a function of R at different F bc in Fig. B1.We can see that Rtot increases with R and decreases with F bc .When F bc = 1, i.e., there is no diffuse dust, Rtot = 1; When F bc = 0, i.e., there is no BC dust, Rtot = R.It is worthwhile noting that the total dust disc can be well described by an exponential profile in most cases.In some extreme cases, e.g. at R < 0.3, it slightly deviates from an exponential profile.This will not affect our results since our model fitting suggests the R ≈ 0.5.

APPENDIX C: CALIBRATING THE TOTAL INFRARED LUMINOSITY FROM WISE MID-IR BANDS
We calibrate the total infrared luminosity by making use of the high-quality data from the DustPedia survey (Davies et al. 2017;Clark et al. 2018). 10The total IR luminosities are taken from Nersesian et al. (2019), in which CIGALE was adopted to perform the energy balance fitting and the THEMIS model was used to derive the dust properties.Galaxies that have large errors on either stellar mass, SFR or IR luminosity (i.e., > 0.3 dex) were excluded from our analysis.Our final sample contains 324 SFGs with secure detections in multiple bands.
To better sample the IR SED range and obtain reliable LIR measurements, we require detections with a signal-to-noise ratio S/N > 3 in at least three Herschel bands (70-500 µm).We give a new calibration of total infrared luminosity [LIR (8-1000 µm)] as a function of the combination of WISE 12 µm and 22 µm luminosities, log LIR = 1.19 + 0.97 × (log L22 + 0.7(log L12 − log L22)), (C1) where L12 and L22 are the monochromatic luminosities given by L12 = νLν (12 µm) and L22 = νLν (22 µm) in units of Solar luminosity.Fig. C1 demonstrates that the scatter significantly decreases if two mid-IR bands are used to constrain the total IR luminosity.

Figure 1 .
Figure 1.IRX as a function of Σ tot dust at different F bc (a), C bc (b), R (c), and Ĥ (d), indicated by blue and orange colours.The dashed lines represent the dust attenuation caused by diffuse dust, while the solid lines consider both diffuse and BC dust attenuation (with a fixed τ bc,V = 0.5).All model parameters (F bc , C bc , R, Ĥ) are fixed to a typical value of 0.5 except for the one indicated in the diagram.

Figure 2 .
Figure 2. Corner plot showing the best-fitting results of posterior PDFs of nine model parameters: n, p, τ bc,Z ⊙ (V -band), q, F bc,Z ⊙ , η,C bc,Z ⊙ , ν and R. The 50th percentile value (slightly different from the likelihood-weighted mean value in Table1) and 1 σ upper and lower uncertainties (obtained based on the 16-50-84th percentile values) are given at the top of each histogram.

Figure 4 .
Figure 4.The power-law slope of the universal IRX relation IRX = 10 α SF R β R −γ e (b/a) −δ as a function of metallicity (dashed lines).Results from our best-fitting two-component geometry model are overplotted as solid lines.

Figure 5 .
Figure 5. Top: Schematic diagrams of the two-component dust geometry model in galaxies at low (top-left) and high (top-right) metallicities.The smooth background with density gradient represents the distribution of diffuse dust, while the discrete filled clumps represent the dense clouds and open circles represent the dispersed clouds due to stellar feedback from newly-born stars.The blue stars represent the young and intermediate-age stellar populations with UV emission, and red stars represent old stellar populations.Bottom: IRX as a function of Σ dust at low (bottom-left) and high (bottom-right) metallicity.The blue and green lines represent the model-predicted dust attenuation caused by diffuse and dense (BC) dust, respectively, while the red curve represents the total dust attenuation.The grey points are our local SFGs and the magenta circles with error bars represent the high-z SFGs in an average sense.Here the (projected) dust surface density is calculated from metallicity, SFR, Re, and b/a according to Equation 20.

Figure 6 .
Figure6.Top: the predicted Σ dust as a function of stellar mass at different redshifts over 0 < z < 3. The prediction is made with our two-component star-dust geometry model with the best-fitting parameters to the local SFGs in combination with the empirical mass-SFR, mass-metallicity, mass-size scaling relations, the K-S Law and the Z-DGR relation from the literature.Middle: The predicted IRX as a function of stellar mass for redshifts 0 < z < 3, obtained by pairing our dust-star geometry model with the aforementioned empirical scaling relations.Bottom: Solid symbols show the predicted IRX from our star-dust geometry model, contrasted to the IRX predicted by the universal IRX relation (Equation14).The open and solid symbols of the same size and shape correspond to the same galaxy subpopulations.Clearly, the observed mass-IRX relations given inWhitaker et al. (2014) can be well reproduced in the sense that IRX mildly evolves with redshift at the high-mass end but shows no evolution at low-or intermediate stellar masses.When replacing M * with the predicted IRX by the universal IRX relation, the scattered data points (the open symbols) at the high-mass end in the middle panel redistribute to form a tight relation (the solid symbols) in the bottom panel.

Figure 7 .
Figure7.The growth history of a Milky-Way like galaxy in given parameters/quantities: star formation rate (SFR), stellar mass (M * ), half-light radius (Re), and gas-phase metallicity [12+log(O/H)] adopted from FörsterSchreiber & Wuyts (2020).The predicted IRX (black) and the best-fitting parameters τ bc,V (cyan), F bc (orange) and C bc (green) from our two-component star-dust geometry model are shown for comparison.

Figure A1 .
Figure A1.Top panels: IRX a function of Σ dust with disc thickness H d /R d of 0.1 (top-left), 0.5 (top-middle), and 1 (top-right).The solid lines with blue-to-red colour represent the increasing value of inclination θ.Bottom panels: similar to the top but replacing Σ dust with projected dust surface density Σ dust,proj = Σ dust ÷ (b/a).

Figure B1 .
Figure B1.Rtot as a function R at different F bc .Here, Rtot represents the size ratio of total dust and diffuse dust disc, R the size ratio of SFR (or BC) and diffuse dust disc, and F bc the fraction of dust comprised in BCs.
APPENDIX B: DETERMINING THE NUMERICAL SOLUTION Rtot ( R, F bc ) Figure C1.L IR as a function of L 12 , L 22 and their combination for the DustPedia sample.The red lines depict the corresponding bestfitting relations and the 1 σ scatter of the sample around this relation is also given.

Table 1 .
The model parameters of our geometry model