The slippery slope of dust attenuation curves: Correlation of dust attenuation laws with star-to-dust compactness up to z = 4

Aims. We investigate dust attenuation of 122 heavily dust-obscured galaxies detected with the Atacama Large Millimeter Array (ALMA) and Herschel in the COSMOS field. We search for correlations between dust attenuation recipes and the variation of physical parameters, mainly the effective radii of galaxies, their star formation rates (SFR), and stellar masses, and aim to understand which of the commonly used laws best describes dust attenuation in dusty star-forming galaxies at high redshift. Methods. We make use of the extensive photometric coverage of the COSMOS data combined with highly-resolved dust continuum maps from ALMA. We use CIGALE to estimate various physical properties of these dusty objects, mainly their SFR, their stellar masses and their attenuation. We infer galaxy effective radii (Re) using GALFIT in the Y band of HSC and ALMA continuum maps. We use these radii to investigate the relative compactness of the dust continuum and the extension of the rest-frame UV/optical Re(y)/Re(ALMA). Results. We find that the physical parameters calculated from our models strongly depend on the assumption of dust attenuation curve. As expected, the most impacted parameter is the stellar mass, which leads to a change in the"starburstiness"of the objects. We find that taking into account the relative compactness of star-to-dust emission prior to SED fitting is crucial, especially when studying dust attenuation of dusty star-forming galaxies. Shallower attenuation curves did not show a clear preference of compactness with attenuation, while the Calzetti attenuation curve preferred comparable spatial extent of unattenuated stellar light and dust emission. The evolution of the Re(UV)/Re(ALMA) ratio with redshift peeks around the cosmic noon in our sample of DSFGs, showing that this compactness is correlated with the cosmic SFR density of these dusty sources.


Introduction
In its earlier years, the Universe did not only witness more star formation rates (SFRs), but a significant fraction of this SFR was obscured by dust (e.g., Blain et al. 2002;Takeuchi et al. 2005;Chapman et al. 2005;Bouwens et al. 2012;Madau & Dickinson 2014;Magnelli et al. 2014;Bourne et al. 2017;Whitaker et al. 2017;Gruppioni et al. 2020;Khusanova et al. 2020). The heavily attenuated in the ultraviolet (UV) dusty star-forming galaxies (DSFGs, e.g., Viero et al. 2013;Weiß et al. 2013;Casey et al. 2014;Béthermin et al. 2015;Strandet et al. 2016;Casey et al. 2017;Reuter et al. 2020) have contributed notably to the cosmic SFR, making them crucial to the general comprehension of galaxy evolution. These dust-rich galaxies have managed to stack up their stellar masses in relatively short timescales, while efficiently depleting their gas reservoirs. Therefore, they might be the direct progenitors to ultramassive passive red galaxies that are often encountered at high redshift (e.g., Daddi et al. 2005;Whitaker et al. 2013;Nayyeri et al. 2014;Toft et al. 2014;Carnall et al. 2020).
Our understanding of the nature of DSFGs have benefited from the plethora of multiwavelength photometry over the last two decades. With a decade worth of observation programs of the Atacama Large Millimeter/submillimeter Array (ALMA), it has become possible to enrich our datasets of infrared (IR) galax-ies at higher redshifts, by their multiplicity from their blended, lower resolution detections with Spitzer and Herschel. This has made modeling dust attenuation at high redshifts a hot topic (e.g., Popping et al. 2017;Wang et al. 2017;McLure et al. 2018;Buat et al. 2019;Salim & Boquien 2019;Fudamoto et al. 2020).
Interstellar dust is highly efficient in absorbing the short wavelength photons, predominantly originating from young UVbright massive stars, rendering the star-forming cold regions of the DSFGs virtually inaccessible. This attenuated light can be successfully reproduced by assuming a dust attenuation law (e.g., Burgarella et al. 2005;Buat et al. 2012Buat et al. , 2014Lo Faro et al. 2017;Salim et al. 2018;Salim & Narayanan 2020). In reality however, a single attenuation law cannot mimic dust extinction in the interstellar medium (ISM) of a large and diverse sample of galaxies (e.g., Wild et al. 2011;Kriek & Conroy 2013;Buat et al. 2018;Małek et al. 2018;Salim & Narayanan 2020). Different approaches appear to work when modeling dust attenuation in reproducing spectral energy distributions (SEDs) of galaxies. Calzetti et al. (2000) derived a universal effective attenuation law which is used as a screen model, by measuring the extinction in local starburst galaxies. The attenuation curve of Calzetti et al. (2000) succeeds in modeling dust reddening even at high-redshift metal-poor galaxies with bright cold dust component. The double component dust attenuation law of Charlot & Fall (2000) assumes a more complex, physical mixing of dust Article number, page 1 of 18 arXiv:2304.13713v1 [astro-ph.GA] 26 Apr 2023 A&A proofs: manuscript no. aanda and stars. With this approach, newly formed stars are placed in the cold molecular clouds, and experience double attenuation by dust of the molecular clouds and the ISM. Older stars are attenuated by the dust grains of the ISM alone. These attenuation laws, together with their different recipes, are often used in the literature when modeling the SED of galaxies. Such recipes include a steeper curve than that of Calzetti et al. (2000) with a UV bump around 0.217µm (Buat et al. 2011(Buat et al. , 2012, and a shallower ISM attenuation than that of Charlot & Fall (2000) (e.g., Lo Faro et al. 2017. However, when it comes to DSFGs, a proper dust attenuation curve that is able to mimic accurately the absorbed photons is crucial for any SED analysis, especially that in these objects, as dust plays a massive role in their evolution, and is an important agent in their SEDs. On the other side of the SED, the long-wavelength cutoff of the Rayleigh-Jeans approximated tail is extensively covered by ALMA. Therefore, its cold dust continuum provides a vital element in the far-infrared (FIR) SEDs of DSFGs. However, high-resolution ALMA-detected cold dust emission maps were often found to disagree spatially with detections at shorter wavelengths, such as the UV-emitting star-forming regions or the stellar populations of DSFGs (e.g., Dunlop et al. 2017;Elbaz et al. 2018;Buat et al. 2019). This disagreement can be either in a form of different dust compactness relative to the higher energy detections, or in some cases a complete physical dissociation of these two components. Such offsets are often found to be more significant than the systematic offsets that arise from instrumental uncertainties or large beam sizes . In most cases, radio detections with the Karl G. Jansky Very Large Array (VLA) confirmed the ALMA-detected physical dissociation (e.g., Rujopakarn et al. 2016;Dunlop et al. 2017;Elbaz et al. 2018;Hamed et al. 2021). This separation challenges a local energy balance which could be an issue for SED modeling that takes into account a global conservation of energy. This problem was highlighted lately in Buat et al. (2019), since such energetic balance is the core of widely used SED fitting tools (e.g., da Cunha et al. 2008;Noll et al. 2009;Boquien et al. 2019).
Reverse engineering the spectral distribution of DSFGs does not come without obstacles, despite the advent in current understanding of the physical and chemical processes that such galaxies undergo. The most commonly confronted obstacle in SED fitting is the degeneracy problem. Such degeneracy arises from overlapping different physical contribution in one specific wavelength domain, such as the dust and stellar age degeneracy (e.g., Hirashita et al. 2017). Although this can be overcome by assuming an energetic balance and using the FIR emission as an additional constrain for dust attenuation, and by choosing an appropriate star formation history (SFH) (see Ciesla et al. 2016), a well-constrained attenuation curve will help limit such degeneracies, leading to better estimated physical parameters.
Dust attenuation curves significantly alter the stellar mass determination of galaxies in general and of DSFGs in particular. As flatter and geometrically complex attenuation curves can dim the light coming from the older stellar populations in the ISM more efficiently than steeper ones, they naturally result in a significant hidden older stellar population in the optical to nearinfrared (NIR) range. This, along with other assumptions such as the initial mass function (IMF) and the SFH, lead to large differences in the resulting stellar masses (e.g., Zeimann et al. 2015;Małek et al. 2018;Buat et al. 2019).
The uncertainty in stellar mass determination of DSFGs (see the stellar mass controversy in Hainline et al. 2011;Michałowski et al. 2012;Casey et al. 2014), remarkably influences the posi-tion of these objects along the commonly named main-sequence (MS) of star-forming galaxies (e.g., Brinchmann et al. 2004;Noeske et al. 2007;Daddi et al. 2010;Rodighiero et al. 2011;Lo Faro et al. 2015;Schreiber et al. 2015;Hamed et al. 2021). Most galaxies seem to follow the tight scatter of the MS independently of the redshift, whose outliers are typically referred to as starbursts. However, the strong dependence of the stellar mass estimation from SED fitting techniques on the assumed attenuation curve, varies massively this scatter and affects the "starburstiness" of already active DSFGs. It is therefore crucial to choose suitable attenuation laws in order to limit biases on the stellar mass determination.
Recently, Donevski et al. (2020) investigated the dust and gas contents of a large sample of DSFGs, linking dust abundance to other physical characteristics such as their SFRs. Despite the growing understanding of these objects, we still lack a complete picture of how dust attenuates their stellar radiation. Properly quantifying dust attenuation of DSFGs is crucial for quantifying the cosmic SFR and get a better grip on galaxy evolution. To achieve this goal, we study closely the dust attenuation curves in DSFGs and investigate the possible physical properties that they might depend on. With the ever-growing understanding of the nature of these dusty galaxies, in recent years many works studied the relation between dust attenuation and other physical properties (e.g., Fudamoto et al. 2020;Lin et al. 2021;Lower et al. 2022;Boquien et al. 2022). In these studies, links between dust attenuation properties and other physical observables were studied, such as the dust grain sizes, star-dust geometry, and the star formation activity in these galaxies. The most widely used attenuation laws are that of Calzetti et al. (2000) and Charlot & Fall (2000). Both of these laws are used interchangebly when modeling galaxies' photometry at high redshifts. However, we still lack a complete knowledge of the preference of attenuation laws in DSFGs at different redshift range. In this paper, we aim at answering the question of which of the commonly used laws best describes dust attenuation in DSFGs at high redshift. Measuring the slopes of attenuation is beyond the scope of this paper. We make use of a large statistical sample with available dust continuum maps and their UV/optical counterparts. We study the effect of UV/optical to dust continuum compactness on the preferred attenuation law for our sample. This paper is structured as follows: in Section 2 we describe the data analyzed in this work, both the photometry and the images. In Section 3.1 we detail the method we use to estimate the circularized effective radii of our sources, and Section 3.2 provides a description of the SED fitting procedures we used to achieve the physical properties of our sources. The results and their respective discussions are presented in Section 4, and the summary is concluding this paper in Section 5. Throughout this paper, we adopt the stellar IMF of Chabrier (2003) and a ΛCDM cosmology parameters (WMAP7, Komatsu et al. 2011): H 0 = 70.4 km s −1 Mpc −1 , Ω M = 0.272, and Ω Λ = 0.728.

Sample selection
The large two square degree COSMOS field, centered on R.A., Dec. = (10h00m27.9171s, +02d12m35.0315s) (Scoville et al. 2007;Ilbert et al. 2013;Laigle et al. 2016), has been observed with an unmatched commitment from different instruments, covering a wide range of wavelength observation of galaxies up to redshift of six. This unique survey design offers rich data sets  Table 1: Summary of available photometric data in each band with its centered wavelength, the mean of S:N, and the number of detections in our sample. The detections of different ALMA bands (6 and 7) concern different galaxies.
of millions of identified galaxies, allowing deep investigation of galaxy evolution at various redshifts.
The choice of COSMOS field galaxies in this work is motivated by the abundance of multiwavelength data spanning across a wide range of redshifts, and the significant number of ALMA detections that this field enjoys. This set of optical, infrared and submillimeter detections allows us to build a statistical sample of DSFGs. This sample is ideal to investigate the evolution of DSFGs properties and the crucial role that dust attenuation plays in their evolution.

ALMA data
Since the main science objective in this study is to quantify the effect of the distribution of dust emission relative to the stellar continuum, at different redshift ranges, the core of our sample was built around ALMA detections. For that we used ALMA fluxes and continuum maps from the A 3 COSMOS automated ALMA data mining in the COSMOS field (Liu et al. 2019). This data set assembles hundreds of identified galaxies from the ALMA archive into a single catalogue. In our work we use the primary-beam-corrected ALMA maps.
The main advantage in our work is having access to dust continuum morphology relative to the spatial distribution of the star formation region and the stellar population in our sample. For this reason, we carefully select a sample of galaxies characterized by good quality detections.
For submillimeter images, we use the A 3 COSMOS generated maps (Liu et al. 2019). These images were deconvolved using a robust cleaning with a Brigg's parameter of 2, i.e. a natural weighting of visibilities. This results in significantly better signal-to-noise (S:N) of the innermost regions of a source in the uv plane, thus defining its outermost "borders" in the im-age plane (see section 2.1 in Liu et al. 2019 for a more detailed description of produced ALMA continuum images). To study dust attenuation through cosmic time, we select ALMA-detected galaxies with S:N higher than 5. Our preliminary sample is composed of 1,335 ALMA detected individual galaxies.

Ancillary maps
For the shorter wavelength (rest-frame UV and optical) images, we used the third data release of the deep field continuum maps detected with the y band of the Hyper Suprime-Cam (HSC) of Subaru (Miyazaki et al. 2018;Aihara et al. 2022). These images have high angular resolution (0.64 ) which allows a physical comparison with their ALMA counterparts. Comparison between the photometric redshifts of 43 galaxies from the final sample, for which spectroscopic redshifts are available from the literature (Liu et al. 2019). Spearman's rank correlation coefficient is shown as ρ, with the photometric redshift accuracy given by the σ MAD (Ilbert et al. 2009). Red dotted lines correspond to z p = z s ± 0.15(1 + z s ) (Ilbert et al. 2009).

Auxiliary photometric data
We used the photometric data from the Herschel Extragalactic Legacy Project (HELP) panchromatic catalogue (Shirley et al. 2019). This catalogue was achieved taking into account visible to mid-infrared (MIR) range surveys as a prior, homogenizing them and extracting fluxes from Herschel maps whose angular resolution is significantly lower than their short wavelength counterparts.
To build the HELP catalogue, Herschel fluxes were extracted using the probabilistic deblender XID+ (Hurley et al. 2017), which was used on SPIRE maps, taking into into account the positions of sources observed with the high-resolution detections from Spitzer at 24µm. This technique was shown to increase the accuracy of photometric redshift estimations (Duncan et al. 2018).
We perform a positional cross-match between ALMA catalogue and the HELP catalogue with a rather conservative 1 search radius. Although galaxies whose dust continuum or molecular gas emission is significantly dissociated from the shorter wavelength continua are not uncommon (e.g., Elbaz et al. 2018;Hamed et al. 2021), due to various factors such as astrometry problems, positional errors as a result of the beam size, or otherwise physical factors, this timid search radius avoids false matches (e.g., Buat et al. 2019;Liu et al. 2019). This crossmatching procedure resulted in 383 individual galaxies ranging from z=0.3 to z=5.5.
To better constrain the physical properties of our statistical sample through SED modeling, we discarded galaxies that have less than three detections in the UV-NIR wave range. This requirement rejected 30% of the sources. Additionally, we required minimum six detections in the MIR-FIR bands (8-1000 µm) out of which at least five detections having a S:N>3. As a result of these selection criteria, the finally-selected sources have at least 10 detections in the UV to NIR (0.3-8 µm) range with S:N>5. For objects that had detections in similar bandpass filters, we used the detections coming from the deeper survey. This has a high importance in SED fitting procedures especially because measurements at similar wavelengths may be order of magnitudes different from each other, which will negatively affect the quality of their fitted spectrum. Moreover, dense coverage of a very short part of the SED could add too much weight during the SED fitting process and bias the final fit ). Our UV-NIR photometric data as well as the FIR counterparts have overall high S:N (mean of 60 for the former and 13.23 for the latter), while in the IR side of the spectrum, MIPS measurements at 24 µm, as well as SPIRE detections at 250 and 350 µm, and evidently all ALMA detections have high S:N (mean of 12.68 for all bands). All Herschel's SPIRE fluxes (at 250, 350 and 500 µm) are essential since they cover the thermal part of the total SED up to z ∼ 4, which contains the Rayleigh-Jeans dust emission tail. To even better constrain the IR-submilimiter part of the SED fits, we appended our sources with VLA detections at 3 GHz from Smolčić et al. (2017).

Final sample
The above described selection yields in the final sample of 122 galaxies with panchromatically high S:N, covering a redshift range of 1 < z < 4. Forty-three galaxies of our sample have spectroscopic redshifts from Liu et al. (2019), and for the rest of the sample we use the reliable photometric redshifts provided by the HELP catalogue. Figure 1 shows a comparison between the spectroscopic and photometric redshifts of the galaxies in our sample that possess both measurements. We calculate the photometric redshift accuracy (Ilbert et al. 2009) as σ MAD = 1.48 × median(|z p − z s |/(1 + z s )) where MAD is the mean absolute deviation. This redshift accuracy resulted in reliable photometric redshifts of our sample with σ MAD = 0.014. Table 1 shows the photometric bands and the associated S:N for the final sample of 122 DSFGs. Almost half of the final sample, 60 galaxies, have a S:N of Y band detection higher than 5, and all of the sources had a VLA detection. We want to stress that only 15 DSFGs from our sample (12%) do not have u band detection, and we miss five detections in g and V bands. With the exception of those galaxies, the rest of the sample have full set of 23 photometric bands, assuring excellent spectral coverage, essential for a detailed SED fitting.

Size measurements
To study the spatial extent of dust emission and that of the stellar populations and the star-forming regions in our sample, we derived homogeneous effective radii (R e ) of the dust continuum maps and their short (UV-optical) wavelength counterparts. To achieve this, we used GALFIT (Peng et al. 2002), parametrically fitting two-dimensional Sérsic profiles in the primary-beamcorrected images of our sample. With the Sérsic index (n Sérsic ), obtained from the fitting procedure, it is possible to quantify the concentration of light in a galaxy, which can provide important information about its morphology. Moreover, GALFIT provides the values of the ratio between the minor and major axis, and this allows for calculating the effective circularized radii (hereafter simply effective radii or R e ). With spectroscopic and reliable photometric redshift information, we can convert the real sizes of our objects. In this work, we analyze the evolution of IR radii between redshift four and one, but also simultaneously changes in the rest-frame UV-optical.
For our analysis we used the Y band of HSC, and ALMA (bands 6 and 7) images to compute physical sizes of our final sample of DSFGs. The mean angular resolution of the ALMA detections in our sample is 0.75 , and equivalently 0.64 for the HSC Y band detection. This renders valid a comparison between sizes estimated from each of these bands directly, without the need to degrade them on resolution.
In computing the Sérsic profiles, we adopted a similar approach than Elbaz et al. (2018), by leaving (n Sérsic ) free. For the sake of comparison, we also fixed n Sérsic = 1 (e.g., Hodge et al. 2016;Elbaz et al. 2018). While 5% of galaxies in our sample did not reach convergence when fixing n Sérsic , R e was found to be in good agreement in both cases (on average 7% larger when fixing n Sérsic for ALMA detections), and we found an agreement within 23% in resulted n Sérsic with a fixed and a free n Sérsic .
In our profile fittings, we initially used an automated computation with GALFIT, and carefully checked the resulted models and their residuals. This was performed to get a general understanding of our sample's range of effective radii in each of the two used wavelength domains. This allowed us to test different input values of n Sérsic , including a Gaussian model. This is an important step when using this tool, since it permits validating the initial parameters needed for GALFIT. We found that varying n Sérsic between an exponential disk profile (n Sérsic = 1) and a Gaussian profile (n Sérsic = 0.5) leads to a slight change in the models and their consequent effective radii, with R e n=0.5 ALMA being 26% smaller than R e n=1 ALMA .
After the aforementioned tests, we individually fitted Sérsic profiles to each of our sources in the two bands (UV/optical and IR), while monitoring the resulted models, the residuals and the profile parameters. This individual fitting was especially required for galaxies that did not fit in the automated computation (∼ 5% of the total sample), in which case a simple and slight parameter adaptation managed to fit these sources. The distribution of the computed n Sérsic was rather narrow, ranging Our technique in fitting Sérsic profiles in two bands of our sample was performed homogeneously in the same methodical approach. This is important to accurately quantify ratios of different R e at varying wavelengths. Our primary effort was to calculate our effective circularized radii in a homogenized method -with the same tool and approach, which reduces possible biases for the final physical interpretation. To test the reliability of our size measurements, we also computed the minimum possible size that can be accurately measured using the formula by Martí-Vidal et al. (2012) and Gómez-Guijarro et al. (2022): where the minimum size for each source (θ min ) is in units of the synthesized beam FWHM (θ beam ), depending on the S:N of the source. All size measurements of our sample were above that limit. Figure 2 shows an example of a galaxy HELP-J095953.305+014250.922 at redshift z=2.15 seen in the HSC Y and ALMA band 6. Figure 2 shows the original image of that galaxy at two different wavelengths, with their light profiles fitted with GALFIT, and their subsequent residuals.  Table 2: Summary of the derived effective radii of our sample from the available detections at two different wavelengths.
The median radii and their errors (the median absolute deviation) at different redshift ranges are presented in Table 2, and the redshift evolution of these derived sizes along with the evolution of their star-to-dust compactness will be shown in Figure 11 (left panel). This Figure shows the change in the radii sizes (also related to the sample selection of DSFGs) of ALMA and HSC y as well as a comparison with similar work done by Buat et al. 2019 based on galaxies at z∼2. The evolution of these derived radii with the dust luminosity and stellar mass of our sample are shown in Appendix B.

SED fitting method
To derive the physical properties of our well-constrained multiwavelength sample, we use the Code Investigating GALaxy Emission CIGALE 1 , an energy balance technique of the SED fitting (Noll et al. 2009;Boquien et al. 2019). This technique of SED fitting takes into account the balance between the energy absorbed in the rest-frame UV-NIR part of the total galaxy emission, and its rest-frame IR emission. The mediator agent in this energy balance is the dust, since it will absorb a significant part of the short wavelength photons emitted by the stars, and emanate in the form of thermal emission in the FIR. Reverse-engineering the total spectrum of a galaxy is not an easy quest. Some physical processes are completely unrelated, such as the synchrotron emission of accelerated electrons, which dominates the radio part of the SED, and the UV photons whose origin is traced directly to the young stars. However, some physical processes release photons of the same frequency range, such as the MIR range which can have different contributors like Active Galactic Nuclei (AGN) and polycyclic aromatic hydrocarbons (PAHs), resulting in degeneracies. Therefore, carefully choosing physically-motivated templates and parameters is crucial in order to deduce key physical properties that galaxies are experiencing, since these parameters depend on the assumptions made (e.g., Ciesla et al. 2015;Leja et al. 2018;Carnall et al. 2018).
In the following subsections, we describe the different aspects of our SED fitting strategy, and motivate our choice of certain laws and parameters. The SED modules description used in our work are presented below with the stellar part and the dust part.

Stellar SED
To forge the SED of a galaxy, we first assume a stellar population that is behind its direct and indirect emission. This means that we should take into account a stellar population library 1 https://cigale.lam.fr/ and its spectral evolution, stars with different ages and a certain metallicity. In this work, we use the stellar population library of Bruzual & Charlot (2003), a solar metallicity and an initial mass function (IMF) of Chabrier (2003), which takes into account a single star IMF as well as a binary star systems.
The stellar population models are then convoluted with an assumed star formation history (SFH). SFHs are sensitive to many complex factors including galaxy interaction, merging, gas accumulation and its depletion (e.g. Elbaz et al. 2011;Ciesla et al. 2018;Schreiber et al. 2018;Pearson et al. 2019). The SFHs have a significant effect on fitting the UV part of the SED, and consequently affecting the derived physical parameters such as the stellar masses and the SFRs. Ciesla et al. (2016Ciesla et al. ( , 2017 showed that simple SFH models (such as a delayed model) are not enough to reproduce a precise fit of the UV data, especially for galaxies that are undergoing a starburst or quenching activity.
To model the SEDs of our IR-bright sample, we use a delayed SFH with a recent exponential burst (e.g. Małek et al. 2018;Buat et al. 2018;Donevski et al. 2020). This recent burst is motivated by the ALMA detection which makes very likely a numerous populations of young stars manifesting their presence through dust. And in such scenario, a galaxy will build the majority of its stellar population in its earlier evolutionary phase, then the star formation activity slowly decreases over time. This is followed by a recent burst of SFR. The SFR evolution over time is hereby modeled with: where the first term translates into a delayed SFH slowed by the factor of τ 2 , which is the e-folding time of the main stellar population, and extended over the large part of the age of the galaxy. The second term is the exponential decrease of recent SFR, where t and τ are the age of the burst and the e-folding  (Calzetti et al. 2000) Colour    time of the burst episode respectively. This SFH provided better fits compared to the simpler delayed SFH, especially for nonquiescent galaxies. We vary τ as shown in Table 3, to give a comprehensive flexibility of the delayed formation of the main stellar population. We discuss the choice of this SFH over a truncated version in Appendix A.

Dust SED
The dust content of our sample of DSFGs is presumed to be the driving component of the shape of their SEDs. This makes the modeling of dust attenuation important to extract accurate physical properties.
To model the effect of dust, we use two different attenuation laws for our SED fitting: the approach of Calzetti et al. (2000, henceforth C00) and that of Charlot & Fall (2000, henceforth Article number, page 7 of 18 A&A proofs: manuscript no. aanda CF00). While these two attenuation laws are relatively simple, they differ on how to attenuate a given stellar population. The attenuation curve of C00 was tuned to fit a sample of starbursts in the local Universe. This curve attenuates a stellar population assuming a screen model: where k(λ) is the attenuation curve at a given wavelength λ, A(λ) is the extinction curve, and E(B-V) is the color excess, which is the difference between the observed B-V color index and the intrinsic value for a given population of stars. Despite its simplicity, this attenuation curve, with its modifications, is widely used in the literature (e.g. Burgarella et al. 2005;Buat et al. 2012;Małek et al. 2014Małek et al. , 2017Pearson et al. 2017;Elbaz et al. 2018;Buat et al. 2018;Ciesla et al. 2020). However, it does not always succeed in reproducing the UV extinction of galaxies at higher redshifts (Noll et al. 2009;Lo Faro et al. 2017;Buat et al. 2019).
Another approach is to also consider dust present in birth clouds. This is the core of the attenuation curve of CF00. In this approach, dust is considered to attenuate the dense and cooler molecular clouds (hereafter MCs) differently than ambient diffuse interstellar media (ISM). This configuration is expressed by the following analytical expression: where δ IS M and δ MC are the slopes of attenuation in the ISM and the MCs respectively. Young stars that are in the MCs will therefore be attenuated twice: by the surrounding dust and additionally by the dust in the diffuse ISM. A ratio of Av IS M is also considered to account for the attenuation of young stars residing in the birth clouds, and the older stars residing in the ISM. CF00 found that δ IS M = δ MC = -0.7 satisfied dust attenuation in nearby galaxies, however, this curve is frequently used at higher redshifts (e.g. Małek et al. 2018;Buat et al. 2018;Pearson et al. 2018;Salim & Narayanan 2020). By attenuating at higher wavelengths (until the NIR) more efficiently than C00, this approach considers a more attenuated older stellar population.
Lo Faro et al. (2017, henceforth LF17) have found that a shallower attenuation curve reproduces the attenuation in ultraluminous and luminous IR galaxies (ULIRGs and LIRGs) at z∼2. For their sample, LF17 have found δ IS M = -0.48. This curve was used in Hamed et al. (2021) for a heavily dust-obscured ALMA-detected galaxy at z∼2, and provided an overall better fit that other steeper attenuation laws.
To model dust attenuation of the galaxies of our sample, we use the aforementioned laws, with the parameters presented in Table 3. The mean normalized attenuation of our sample as a result of the three attenuation curves are shown in Fig. 3

left panel.
To obtain the left panel of Fig. 3, we computed the attenuation values in the UV-NIR bands. Then we averaged the attenuation in each band for the whole sample. The curve of C00 is steeper than the double-component power-laws of CF00 and LF17, especially in the NIR domain.

Hot and cold dust components
Dust grains heated by AGN, along with the vibrational modes of polycyclic aromatic hydrocarbons (PAHs), dominate the MIR part of the SED of a galaxy. Thus, it is important to include AGN modeling in our SED fitting procedures, as well as taking into account PAH contribution to the overall dust emission. Our initial analysis of the IRAC photometry of Spitzer did not suggest Article number, page 8 of 18 M. Hamed et al.: The slippery slope of dust attenuation curves AGN candidates. We also included AGN-heated dust templates of Fritz et al. (2006) in our SED fitting procedure, but found no AGN contribution in our sample.
To model the IR emission in our SED models, we use the templates of Draine et al. 2014. These templates take into consideration different sizes of grains of carbon and silicate, hence, allowing different temperatures of dust grains. These templates rely on observations and are widely used in the literature to fit FIR SEDs.

SED quality, model assessment
In assessing which SED provides the best fit for modeling the galaxies of our sample, we adopt a similar methodology as in Buat et al. (2019). We compare the reduced χ 2 of the resulted fits with the three attenuation recipes used, and in the case of different attenuation laws providing a good reduced χ 2 , we checked the Bayesian inference criterion (BIC) defined as BIC = χ 2 + k × ln(N), where k is the number of free parameters and N is the number of data points. The mean reduced χ 2 was found to be 2.4 for our sample (2 left panel, where we show the reduced χ 2 for the best fits of our sample).
To test the robustness of the method we use to select the best attenuation law for each galaxy, we perform the following test: -We fit all galaxies of our sample with the attenuation law of CF00. We then take the best fit values for each filter of each galaxy. -We perturb these "best" fluxes (obtained from the fit) using the initial photometry errors to obtain a mock catalog. -This mock catalog was then fitted with: Calzetti law, CF00 law and LF17 law, in the same way as the initial real photometry was treated.
The reduced χ 2 obtained from those fits are then compared for each source. We show this comparison in Fig. 6. The obtained reduced χ 2 for CF00 fits of the mock sample are consistently smaller than that obtained using the two other attenuation recipes. Precisely, 93% of mock galaxies preferred CF00 attenuation law using this test. The other 7% had a reduced χ 2 of ∼0.1 lower using the other attenuation laws.
The choice of the best attenuation law that better describes the observed fluxes is crucial for our study. Therefore we used an additional method in reliably attributing the best attenuation law for each galaxy. We introduced perturbations to the fluxes by applying a Gaussian distribution with a standard deviation that corresponds to the uncertainties of each band. To make the calculation time faster, and since the main task was to check the reliability of the best attenuation law for each galaxy, we generated mock fluxes until the mid-IR bands of IRAC. We created ten mock catalogues with this method, and applied the same approach of SED fitting for our original sample to the mock samples. Fig. 7 shows the ratio of mock fluxes to real fluxes of our sample.
The results of these tests are presented in Tab. 4, where the vast majority of our mock galaxies (> 94%) preferred the same attenuation law of the initial real galaxies. This shows that the usage of the reduced χ 2 in assessing the best attenuation law in our case is valid. This is directly linked with the good S:N ratios of our photometric data, that were shown in Tab 1. Fig. 6. Comparison between reduced χ 2 of the mock sample obtained with CF00, C00, and LF17 attenuation laws.

Real
Mock C00 CF00 LF17 C00 94% 5% 1% CF00 0% 97% 3% LF17 0% 3% 97% Table 4: Summary of the results of the best fits using the same SED method on the fluctuated photometry preferring the same attenuation law (average of ten realizations) and the best attenuation laws for the real galaxies.
A&A proofs: manuscript no. aanda An example of a computed SED of our sample is given in Fig. 5, where we attenuated the same galaxy in three different ways. In this example, the shallower attenuation of Lo Faro et al. 2017 was preferred since it provided a significantly better fit. To test the reliability of our SED models, with CIGALE we generated a mock galaxy sample and fitted SEDs with the same methods applied to our sample. We show the comparison between the real physical properties that we derived for our sample and its mock equivalent in Fig. D.1.

Galaxy properties and dust attenuation
We applied a bouquet of attenuation slopes in fitting our sample of DSFGs. The attenuation curve of C00 results in lower stellar masses compared to the ones obtained with the shallower double component attenuation laws of CF00 and LF17 (with a mean stellar mass of 10 10.87 M for C00, 10 11.22 M and 10 11.52 M for CF00 and LF17 respectively). The distribution of the obtained stellar masses was shown in Fig. 3 right panel with the mean M * for the whole sample portrayed for every attenuation law used. Stellar masses computed using shallower attenuation slopes are higher than the one produced with steeper curves.
Star formation rates computed from the panchromatic SED fitting using the three aforementioned attenuation laws do not change, similar to the results found in Małek et al. (2018). The mean values of the log 10 (SFR) of the sample fitted with C00, CF00, and LF17 are 2.75, 2.63, and 2.60 M yr −1 respectively, which is of a similar range of galaxies studied in Buat et al. 2019. The dust masses computed with the three attenuation laws are invariant (with a mean of 1.80×10 9 M for the whole sample). This is mainly due to the strong constrain of the FIR part of the SED provided by the ALMA detections as well as the good fitting of the spectrum. The significant difference in produced stellar masses using the different attenuation law slopes result in a clear distinction in the "starburstiness" of galaxies, and also affects the quiescent systems (Lo Faro et al. 2015). In our sample, the number of starburst galaxies decrease with a shallower attenuation curve (60% for C00, 25% for CF00, and 14% for LF17). Despite its simplicity, the attenuation law given by Calzetti et al. (2000) provided good fits in building the SEDs and was favored over the shallower curves of Fall (2000) andLo Faro et al. (2017) in 49% of the whole sample, that is 61 sources (by comparing the resulted reduced χ 2 ). This was mainly noticed below redshift of z = 2. The attenuation curve adapted by Lo Faro et al. (2017) provided better fits for 38 galaxies in total, but 79% of these galaxies fell in the redshift range of 1.5 < z < 2.5, which supports the initial tuning of the ISM attenuation at -0.48 (Lo Faro et al. 2017).
We show the preference of the attenuation laws for our sample based on the V band attenuation and the SFR in Figs. 8 and 9 respectively. We find no clear correlation between the attenuation in the V band and the preference of the attenuation laws in our sample. We also checked the correlation with the stellar masses. But since these masses are directly a byproduct of the attenuation law used, as it was shown in Figs. 3, we cannot tell if this correlation is physical. Galaxies that are preferring CF00 double-component attenuation law and its shallower version LF17, by construction, result in a significantly higher older stellar population, therefore increasing the stellar mass (e.g., Małek et al. 2018;Buat et al. 2019;Hamed et al. 2021;Figueira et al. 2022). We also checked the preference of attenuation laws used with SFR for our galaxies. This is shown in Fig. 9. We found that towards higher SFRs, there are no preferences in attenuation laws for our sample. However, in the lower limit of SFR, a double-component attenuation was slightly preferred, but still within the error bars. In the sample we had only 18 galaxies with log(SFR)< 2.4 and 38% of them were fitted with C00 while the rest preferred CF00/LF17. The small statistics at this lower end of SFR for our sample does not allow us to make a strong statement about the correlation of the attenuation laws for low-SFR galaxies.

1.0
Fraction C00 CF00/LF17 Fig. 8. Preference of attenuation laws of our sample according to the attenuation in the V band. To facilitate the reading of this plot, we shifted the bins of CF00/LF17 slightly to the right (+0.1).

Stellar vs. dust components
To analyze the energy balance of our sample of DSFGs, especially with the available dust and stars emission and images, we follow the method introduced in Małek et al. (2018) and Buat et al. (2019) by dissecting the stellar continuum of our sample of galaxies without taking into account the FIR detections. Equivalently, we also fit the FIR continua of our SFGs. To model the stellar continuum, we use the photometric bands of CFHT, Subaru, VISTA and available IRAC bands. This is done in order to ensure the attenuation curve requirements without adding the energy balance constrain to the global SED fitting method. As shown in e.g., Buat et al. 2019 andHamed et al. 2021, the dichotomy between the stellar SED and its dust counterpart is important in testing the validity of the energy balance concept that is the basis of most of panchromatic SED fitting tools. Moreover, this method is critical in cases where the dust continuum maps are not centered on their short wavelength counterparts. When fitting the short wavelength part of the SEDs, we model the stellar light taking into account the delayed SFH boosted by a recent burst, a stellar population library of Bruzual & Charlot (2003), and dust attenuation laws that are discussed in 3.2. We also modeled the IR emission using Draine et al. (2014) dust emission templates, but without taking into account the IR photometry, allowing the energy balance to dictate dust luminosities and masses based on the amount of the stellar light that is attenuated.
A total of 61% of our sample (75 sources) provided better fits with the simple power law of C00. While a shallower attenuation was needed to reproduce the spectra with the least χ 2 for the rest (14% with CF00 and 25% with an even grayer LF17 slope).
The steep law of C00 provided better fits for more galaxies when taking into account the stellar emission only. This result was also found in Buat et al. (2019) for a smaller sample. This tendency will be confirmed in future studies based on the new generation of IR datasets from JWST and well constrained short wavelength counterparts from LSST.
Equivalently, we estimated dust luminosities based on the rest frame UV to NIR photometry of our sample, assuming an energy balance between dust absorption and dust emission. We compare these IR luminosities with the ones calculated from the IR photometric points with Draine et al. (2014) templates. The results are shown in Fig. 10. We confirm the scatter initially found in Buat et al. (2019) and Hamed & Małek (2022) with dust-rich galaxies which significantly differs from the normal SFGs found in Małek et al. (2018). We find that the dust emission evoked from a pure energy balance based on the short wavelength does not always explain the one calculated from the FIR photometry. Galaxies that are fitted with C00 attenuation law are found to inhibit more attenuation and they produce higher L IR from their dust content rather than the stellar parts. The star-to-dust compactness did not seem to play a role in this trend. This shows that the dust luminosity concluded from the direct and attenuated UV photons based on the energy balance is not enough to reproduce the dust luminosity observed from the actual IR photometry.

Dust attenuation and sizes
We study the effects of star-to-dust compactness of the unobscured star-forming regions and the stellar population emission to the extent of dust emission detected by ALMA. We define the ratio of the short wavelength radii to their FIR counterparts as R e (UV)/R e (ALMA) where R e (UV) is the circularized effective radii measured from the HSC Y bands of our sources. The UV radii of our sample decrease with redshift, while the ALMA counterparts are rather constant accross the studied redshift range in this work. This is shown in Fig. 11 and is explained by the bright IR ALMA-selected sources, and their active star-forming/stellar population regions.
Despite the fact that the attenuation curve of Calzetti et al. (2000) managed to fit the largest sub-sample of our DSFGs around the cosmic noon and at lower redshift ranges, the distribution of these galaxies for which a steeper attenuation curve provided the best fit is relatively compact. A relative star-to-dust compactness was found to have an average of 1.7, while for the shallower attenuation curve recipes were found to not follow a specific preference, and they are rather scattered accross all the studied redshift bins. We show the redshift distribution of the preferred attenuation laws in Appendix B. The range of the resulted (n Sérsic ) was too small (0.4 to 1.6), and in that narrow range we did not find any correlation with other observables. We find that the ratio of R e (UV) to R e (ALMA) changes accross the redshift in our sample of ALMA-detected DSFGs. This distribution is found to peak at z = 2 around the cosmic noon, as shown in the right panel of Fig. 11. This change is more prominent at higher redshift with the decrease of the star-to-dust compactness is directly connected to the rapidly decreasing restframe UV sizes of these DSFGs at higher redshift. This might be explained by the more intense star formation around that cosmic epoch, especially that DSFGs contributed significantly to the total star formation activity in the Universe. Moreover, this peak is found to be stronger for galaxies that require a shallower attenuation curve. This might correlate with the higher need of cold star forming regions in these galaxies to explain the higher SFRs. This result shows that dust attenuation and the star-to-dust compactness of dust in DSFGs might be correlated with the cosmic SFR density.
These findings partly agree with the smaller sample size studied in Buat et al. (2019). However, our statistically larger sample size allows for an extrapolation of this correlation not only at different redshift ranges, but also showed that shallower attenuation curves do not favor higher star-to-dust compactness, but rather a scattered trend was found, unlike the steeper curves which clearly preferred relatively smaller R e (UV)/R e (ALMA) (∼ 2) at different redshift ranges. One possible explanation for this can be the fact that for galaxies with relatively smaller sizes of their unobscured star-forming regions and their dust content, might be translated by a screen model of attenuation, due to a non-effective mixing of stars and dust. On the other hand, a very compact dust emission requires a more complex mixing of dust and stars which is translated in a shallower attenuation curve. A correlation is visible between the fraction of our sample that is fitted with a specific attenuation law and the relative compactness. We show this in Fig. 12.
We find that the galaxies with smaller opt/UV sizes relative to dust radii largely prefer the C00 attenuation law in our sample. For larger R e (Y)/R e (ALMA) ratio (>3) only CF00 and its shallower modification LF17 fitted the galaxies better. This shows that taking relative compactness into account is highly important when performing SED fitting. One can infer the most likely attenuation law for a given galaxy by taking into account its sizes in the unattenuated SF region and the dust continuum. This in turn enables a reliable procedure for fitting its spectral energy distribution (SED).

Summary
In this paper, we studied a statistical sample of 122 DSFGs not hosting AGN accross a wide range of redshift (1 < z < 4). We derived their circularized effective radii in two different bands, HSC's Y band when available, and equivalently their dust components' radii with ALMA detections. On the other hand, we carefully analyzed their SEDs, modeling them with varied dust attenuation laws, particularly that of Calzetti et al. (2000) and the shallower curves of Charlot & Fall (2000) and their recipes.
We also dissected the stellar SEDs alone and their IR counterparts as done in (e.g., Buat et al. 2019;Hamed et al. 2021) to investigate the validity of the energy balance when taking into account the ALMA detections of our sources. We found that even if most of our sources seem to produce the same dust emission when relying on an energetic balance from the short wavelengths, and when fitting the IR photometry apart, some galaxies expressed dimmer star formation when the Calzetti et al. (2000) attenuation law was favored. This was translated into an under-supply of dust emission from the stellar population alone.
We found that knowing the information of the sizes of DSFGs especially their stellar and dust contents in the analysis, constrains the used attenuation curves to fit the photometry of these galaxies. We found that a starburst curve of Calzetti et al. (2000) was favored in the reproduction of the SEDs of DSFGs with comparable star-to-dust radii ratios, precisely at 1.2 < R e (UV)/R e (ALMA) < 2.5. However, despite the seemingly irrelevant star-to-dust ratios with the attenuation curves of Charlot & Fall (2000) and the shallower version (Lo Faro et al. 2017), we found that compact dust emission and extended stellar radii needed shallow curves and double exponential attenuation laws to account for the missing photons absorbed by dust. This shows that when fitting SEDs using broad band photometry, a careful analysis of the radii of different components should be taken into account, before using a unique attenuation law which will result in wrongly-estimated stellar masses.
We stress that recent ALMA studies of smaller samples of z∼2 galaxies suggested that compact IR sizes could be connected to rapid growth of supermassive black holes during the SMG star-formation phase (Ikarashi et al. 2017). Semianalytical models (i.e. Lapi et al. 2018) predict that such sources would experience forthcoming/ongoing AGN feedback, which is thought to trigger the morphological transition from star-forming discs to early-type galaxies. However, our galaxies do not show AGN activity, and are suggestive of SF galaxies caught in the compaction phase characterised by clump/gas migration toward the galaxy center, where the highly intense dust production takes place and most of the stellar mass is accumulated (Pantoni et al. 2021a andPantoni et al. 2021b). Interestingly, dust compaction phase is suggested to play role in metallicity enrichment efficiency, which further affects dust growth in ISM (Pantoni et al. 2019;Donevski et al. 2020). We expect that our dusty galaxies have relatively wide range of gas metallicities and dust growth efficiencies, which would reflect on different attenuation slopes, rather than favoring the single one. This expectation is also in line with results from recent cosmological simulations that found dependence of attenuation on dust compactness and/or geometry (e.g., Schulz et al. 2020) and the ratio between the small and large dust grains (e.g., Hou & Gao 2019).
In our sample, we have observed that the C00 attenuation law is mostly favored by galaxies with opt/UV sizes two times larger than the dust radii. However, for galaxies with R e (Y)/R e (ALMA) >3 the CF00 and its shallower modification LF17 were found to fit better. These findings suggest significant importance of considering the relative compactness when conducting SED fitting. By taking into account the sizes of a galaxy in the unattenuated SF region and the dust continuum, one can deduce the most probable attenuation law for that galaxy, thus providing a dependable approach for fitting its SED.
We conducted a test to ensure that the trend observed between the preferred attenuation law and relative compactness is not influenced by other physical properties. Specifically, we investigated the potential relationship between the relative compactness and two other properties -attenuation in the V band and SFR ( Fig. 8 and Fig. 9). Fig. 13 presents the results of this analysis. We find no significant correlation between these properties and relative compactness. This finding supports the conclusion that the observed trend between relative compactness and preferred attenuation law can be robust and not influenced by these physical properties.
We find that the star-to-dust compactness of the unobscured star-forming regions/stellar population regions to dust emission of these DSFGs peaks around the cosmic noon (z∼2). This is notable at higher redshift. A possible correlation might be with the cosmic SFR density for which the DSFGs were a major contributor in the early Universe.
These results are promising in the era of highly resolved deep-field detections with the LSST and JWST, where dust attenuation and size measurements are getting more precise. Combining these detections with the FIR information especially with ALMA, is unparalleled to deal with the dust attenuation curve problem at different redshift ranges. Histogram of the ratio of HSC's Y band radius to the ALMA radius. The red-filled histograms at different redshift ranges represent galaxies for which a shallow attenuation curve of CF00 or LF17 was critical in order to reproduce the observed stellar light and resulted in better fits. The dashed histogram shows galaxies for which the C00 attenuation curve gave a satisfying fit.