Cold dust and stellar emissions in dust-rich galaxies observed with ALMA: a challenge for SED-fitting techniques

Over the past few years ALMA has detected dust-rich galaxies whose cold dust emission is spatially disconnected from the UV rest-frame emission. This represents a challenge for modeling their spectral energy distributions with codes based on an energy budget between the stellar and dust components. We want to verify the validity of energy balance modeling on a sample of galaxies observed from the UV to the sub-millimeter rest frame with ALMA and decipher what information can be reliably retrieved from the analysis of the full SED and from subsets of wavelengths. We select 17 sources at z~2 in the Hubble Ultra-Deep Field and in the GOODS- South field detected with ALMA and Herschel and for which UV to NIR. rest-frame ancillary data are available. We fit the data with CIGALE exploring different configurations for dust attenuation and star formation histories, considering either the full dataset or one that is reduced to the stellar and dust emission. We compare estimates of the dust luminosities, star formation rates, and stellar masses. The fit of the stellar continuum alone with the starburst attenuation law can only reproduce up to 50% of the total dust luminosity observed by Herschel and ALMA. This deficit is found to be consistent with similar quantities estimated in the COSMOS field and is found to increase with the specific star formation rate. The combined stellar and dust SEDs are well fitted when different attenuation laws are introduced. Shallow attenuation curves are needed for the galaxies whose cold dust distribution is very compact compared to starlight. The stellar mass estimates are affected by the choice of the attenuation law. The star formation rates are robustly estimated as long as dust luminosities are available. The large majority of the galaxies are above the average main sequence of star forming galaxies and one source is a strong starburst.


Introduction
Modeling the spectral energy distributions (SED) of galaxies is commonly used to measure some of the most important physical parameters that quantify galaxy evolution such as the star formation rate (SFR) and the stellar mass. The central idea behind this method is to reconstruct the stellar emission of a galaxy with population-synthesis models, assuming star formation histories (SFHs) of varying complexity (e.g., Walcher et al. 2011;Conroy 2013). Dust however plays a crucial role by strongly affecting and reshaping the SED: it absorbs and scatters stellar photons and thermally emits the absorbed energy in the infrared (IR; λ ∼ 1−1000 µm). To account for the impact of dust, modern physical models link the stellar and dust components: the energy radiated by dust corresponds to the absorbed stellar light.
This energy budget is implicit in radiation transfer models aimed at modeling nearby, spatially resolved galaxies with more or less complex geometries (e.g., De Looze et al. 2014). It is implemented explicitly in a handful of panchromatic SED fitting codes aimed at fitting the SED of galaxies such as MAGPHYS (da Cunha et al. 2008), CIGALE (Noll et al. 2009;Boquien et al. 2019), PROSPECTOR (Leja et al. 2017), and BAGPIPES (Carnall et al. 2018). These SED fitting codes are all based on simplified assumptions; for example, the emission is supposed to be isotropic at all wavelengths, and the effect of dust on the stellar continuum (and nebular emission for BAGPIPES and CIGALE) is accounted for by applying an effective attenuation law. This curve is either measured as for the commonly used starburst law (Calzetti et al. 2000) or built to reproduce average observed trends (Charlot & Fall 2000;Lo Faro et al. 2017;Cullen et al. 2018). These methods are very efficient at fitting very large samples of galaxies. Hayward & Smith (2015) showed that the SED modeling code MAGPHYS recovers physical parameters reasonably well for most of their simulated galaxies.
It has recently been shown that effective attenuation laws are far from universal but actually vary in local and more distant galaxies (e.g., Salmon et al. 2016;Salim et al. 2018;Buat et al. 2018). The attenuation law appears to flatten out when the amount of obscuration increases, a trend that is predicted by radiation transfer models in dusty media with idealized geometries (Chevallard et al. 2013, and references therein) or more realistic distributions from hydrodynamical simulations (Jonsson et al. 2010;Roebuck et al. 2019;Trayford et al. 2019). Lo Faro et al. (2017 found that a very flat attenuation law accurately describes a sample of ultra luminous infrared galaxies (ULIRGs) at z ∼ 2. Such a curve is able to explain the locus of these galaxies in the IRX-β diagram 1 , well above that of starburst galaxies obtained by Meurer et al. (1999). In this case high IRX values for a given β are reproduced by a global attenuation law that is grayer than the starburst curve rather than being due to a disconnection between fully obscured star-forming components and unobscured regions dominating the UV emission.
Energy balance methods assume a physical coupling between the dust and stellar components and should in principle only be valid in this case. Infrared-bright galaxies in which most of the dust emission is not coincident with the emergent stellar emission have long since been identified however (e.g., Charmandaris et al. 2004;Howell et al. 2010). More recently, observations of high-redshift massive dusty and submillimeter (submm) galaxies with the Atacama Large Millimeter Array (ALMA) report a very compact cold dust emission, with a rest-frame UV emission distributed in clumps that are spatially distant from the compact dust emission, and an extended H-band distribution (e.g., Hodge et al. 2016;Tadaki et al. 2017;Gómez-Guijarro et al. 2018;Puglisi et al. 2019). These dusty, compact objects could represent ∼50% of the massive (M > 10 11 M ) main sequence (MS) population (Puglisi et al. 2019). Conversely, by stacking ALMA images of less massive starforming galaxies at z ∼ 2, Lindroos et al. (2016) found typical sizes of ∼5 kpc for the ALMA continuum, similar to the median sizes measured in the near-infrared (NIR). Using observations of galaxies with vigorous star formation at z = 1−3 with both the Kark G. Jansky Very Large Array and ALMA, Rujopakarn et al. (2016) also report star formation distributed over an area similar to the distribution of stellar mass, although ALMA highresolution maps picture complex configurations for a few MS galaxies (Rujopakarn et al. 2019).
Very different dust and stellar distributions obviously challenge a local energy balance which is nevertheless expected to be valid at a global scale provided that we are dealing with a single physical object where dust emission comes from heating by stellar photons. However the way in which both emissions are linked, assuming a simple SFH and dust attenuation recipe, can have an impact on the physical parameters derived from the modeling of the global emission. In this paper, we analyze a sample of z ∼ 2 galaxies observed with ALMA and complemented by a rich multi-wavelength dataset from the rest-frame UV to the submm. We explore which SFHs and attenuation laws are needed to fit the SED of these galaxies and to what extent the SFR and stellar masses determined from the fit, with and without the IR emission emitted by dust, are reliable. With this aim in mind, we use data from the Hubble Ultra-Deep Field (HUDF; Dunlop et al. 2017) GOODS-South field  surveys. We use CIGALE  to fit the SED of these objects with three different attenuation laws that are representative of the very different shapes of effective attenuation laws used in the literature: the starburst law of Calzetti et al. (2000), and the recipes of Fall (2000) andLo Faro et al. (2017).
The paper is organized as follows. The galaxy sample and the multi-wavelength data are presented in Sect. 2. The different configurations of the SED fitting process are explained in Sect. 3 and the fits are performed and discussed in Sect. 4. In Sects. 5 and 6 we discuss the measures of the dust luminosity (L dust ), SFR, and stellar mass, and the locus of the galaxies on the MS. The dust and stellar relative distributions are related to the attenuation law in Sect. 7. We conclude in Sect. 8.

Target selection and multi-wavelength data
We consider the targets from the samples of Elbaz et al. (2018) and Dunlop et al. (2017), hereafter respectively labeled GS and UDF. We first select galaxies with ALMA observations corre-sponding to rest-frame wavelengths ≤500 µm. This corresponds to targets with a redshift z > 1.6 for the sources of Dunlop et al. (2017) observed at 1.3 mm and to z > 0.74 for the sample of Elbaz et al. (2018) observed at 870 µm. We are left with ten objects from Dunlop et al. (2017) and seven sources from Elbaz et al. (2018) (we do not include UDF12 listed with a redshift z = 5 and the optically dark galaxy GS8 at z = 3.2). All the sources but two (UDF10 and UDF15) are part of the extended sample analyzed by Elbaz et al. (2018) and ten sources are detected in X-ray. The redshifts range from 1.6 to 2.8 with a median value at 2.3. The identifiers, redshifts of the sources (with the references), and X-ray detections if any are listed in Table 2.
The objects have been observed from the rest-frame UV to the far-IR. From the shortest to the longest wavelengths: HST/ACS (F435W, F606W, F775W, F814W, F850LP) and HST/WFC3 (F105W, F125W, F160W), the VLT/ISAAC Ks band, Spitzer/IRAC (four bands), the Spitzer/MIPS band at 24 µm, and Herschel/PACS and SPIRE (five bands). The optical to IRAC data come from the multi-wavelength catalog of the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Guo et al. 2013). We cross-match the coordinates of the ALMA sources with those of the CANDELS catalog with a tolerance radius of 1 arcsec, confirming the associations performed by Elbaz et al. (2018), except for UDF 10 and 15 which are not included in the Elbaz et al. (2018) sample. The data from the Herschel GOODS-South field catalog were produced by the ASTRODEEP team 2 (Wang et al., in prep.). Postage-stamp HST images with ALMA contours for each source are shown in Dunlop et al. (2017), Rujopakarn et al. (2016), and Elbaz et al. (2018).

The code CIGALE
The SED fitting is performed with CIGALE 3 . We refer to Boquien et al. (2019) for a detailed description of the code. CIGALE combines a stellar SED with dust attenuation and emission components. A critical aspect is that it conserves the energy between dust absorption in the UV-to-NIR domain and emission in the mid-IR and far-IR. The quality of the fit is assessed by the best χ 2 (and a reduced best χ 2 defined as χ 2 r = χ 2 / (N − 1), with N being the number of data points). The value of the physical properties and their corresponding uncertainties are estimated as the likelihood-weighted means and standard deviations. In this work, we introduce different attenuation models. To compare them we use the Bayesian inference criterion defined as BIC = χ 2 + k × ln(N), with k being the number of free parameters and N the number of data points used for the fit (Ciesla et al. 2018). The BIC is an approximation of the Bayes factor (Kass & Raftery 1995). Below, we describe the assumptions and choices adopted to run CIGALE for our current study. The values of the input parameters are listed in Table 1.

Star formation histories, dust, and nebular components
We assume a delayed SFH with the functional form SFR ∝ t × exp(−t/τ). A current burst of star formation is introduced if it provides better fits. This burst is modeled as a constant star formation that has been ongoing for at most 100 Myr and that is superimposed to the delayed SFH. Its amplitude is measured by the stellar mass fraction produced during the burst episode. The age of the main stellar population and of the burst are also free parameters as well as the e-folding timescale τ; the input values are listed in Table 1. We adopt a Chabrier (2003) initial mass function (IMF) and the stellar models of Bruzual & Charlot (2003 For the sources detected in X-ray an AGN component is added to the fit following the method developed in Buat et al. (2015): the Fritz et al. (2006) library of CIGALE is reduced to two AGN models with an angle between the AGN axis and the line of sight equal to 80 • and either a small or high optical depth at 9.7 µm. The nebular emission is added from the Lyman continuum photons produced by the stellar component. The nebular continuum and the emission lines are calculated from a grid of nebular templates generated with CLOUDY 08.0. In this work the ionization parameter is chosen to be log(U) = −2.

Attenuation laws
We expect the diverse and complex structures of our targets to translate to different effective attenuations. Three different dust attenuation recipes already used in the literature are considered. We do not introduce full flexibility: the large majority of studies using SED fitting methods assume a fixed law. Our aim in this work is to discuss the reliability of physical parameter estimations under simple classical approaches. First, we consider the attenuation law measured for the stellar continuum of nearby starbursts by Calzetti et al. (2000), for which the amount of attenuation is quantified with the color excess A(λ) = E(B − V) star k (λ). The nebular emission is attenuated with a simple screen model and a Milky Way extinction curve with a color excess E(B − V) line . The ratio between the color excess for the stellar continuum and the nebular lines was originally measured to be equal to 0.44 (Calzetti 2001) in local starbursts. At z ≥ 1, Kashino et al. (2013), Puglisi et al. (2016) obtained more similar attenuations for stellar and nebular emissions and Buat et al. (2018) concluded that there is great variability among galaxies with an average value of 0.55. Our targets are highly obscured, with a compact and dusty star formation. Young stars emitting in the UV are expected to experience very high attenuation and we do not have the measurement of emission lines to constrain the differential attenuation. We first tested three different values of the ratio, namely 0.44, 0.6, and 1, with no significant change in the results. We decided to use the same color excess for stellar and nebular emission. The recipe is referred to as C00 hereafter.
Here, we introduce the attenuation recipe of Charlot & Fall (2000), which differs from C00 in its philosophy. A differential attenuation between young (age <10 7 years) and old (age >10 7 years) stars is assumed. Light from both young and old stellar populations is attenuated in the interstellar medium (ISM) but light from young stars is also affected by an extra attenuation in the birth clouds (BCs). Both attenuation laws are modeled by a power law and the amount of attenuation is quantified by the attenuation in the V band, Charlot & Fall (2000) fixed both exponents of the power laws n BC and n ISM to −0.7, although a value of −1.3 is initially introduced in their model and further adopted by da Cunha et al. (2008). Here, we use the original scenario which is more consistent with the studies of HII regions (Lo Faro et al. 2017). The µ parameter is defined as the ratio of the attenuation in the V band experienced by old and young stars Charlot & Fall (2000) obtained µ = 0.3 from their study of nearby starburst galaxies. From a study of Shapes of the attenuation laws introduced for the modeling. The ISM power-laws for CF00 (black) and LF17 (blue) are plotted with dotted lines. The solid lines represent the C00 (red) law and the average attenuation laws obtained for our sample with CF00 and LF17; see text for details. several thousand galaxies detected by Herschel as part of the HELP project, Małek et al. (2018) found that a µ value slightly higher than 0.3 led to better fits. The recipe defined with the exponents of the power laws n BC and n ISM equal to −0.7 and µ equal to 0.5 is referred to as CF00 hereafter.
We further introduce the modification to CF00 as proposed by Lo Faro et al. (2017) to fit ULIRGs at z 2 in agreement with radiation transfer modeling. This attenuation law is flatter with n ISM = −0.48. We refer to this third law as LF17 hereafter. We use the same value of µ = 0.5 in this case, in agreement with the results of Lo Faro et al. (2017).
We compare these three laws in Fig. 1. While the shape of C00 is fixed and applies uniformly to the stellar continuum, the effective shape of CF00 and LF17 instead depends on the SFH since the young and old stars experience a different attenuation. In Fig. 1 we report the average law for our galaxy sample from the fits performed in Sect. 4. We see that CF00 combined with the SFH obtained for our sample leads to an attenuation similar to C00 at short wavelengths. This is expected since Charlot & Fall (2000) built their model to reproduce the properties of the starburst galaxies analyzed by Calzetti et al. (1994). However a significant difference is visible for λ > 0.5 µm: CF00 is much flatter than C00 (Chevallard et al. 2013;Lo Faro et al. 2017;Buat et al. 2018). The average LF17 law is flatter than the two previous laws over the full wavelength range.

Fitting the SEDs
The usual way to model galaxies is to use all the photometric data available over the broadest feasible wavelength range in order to measure the physical parameters as reliably as possible. In the particular case of the galaxies studied in this work, the different locations of dust and stellar emissions lead us to explore different configurations for the fits. We perform three fits here: stellar continuum only, dust emission only, and finally the full SED.

Fitting the stellar continuum
We model the stellar continuum using all the HST, VLT, and IRAC data available. We run CIGALE with the three attenuation laws and with and without a burst added to the delayed SFH. First we test the requirement to add a burst to the delayed SFH. The comparison of the BIC is used for this purpose. Two additional free parameters are introduced when a burst is added to the delayed SFH. The difference between the two BICs corresponding to the models with and without a burst is calculated as ∆BIC = χ 2 (burst) + 2 × ln(N) − χ 2 (no burst), with N being the number of fitted bands. We adopt the limit of |∆BIC| > 6 to put a preference on a model: this represents strong evidence against the model with the higher BIC (Kass & Raftery 1995;Salmon et al. 2016). In all cases, ∆BIC is found between +5 and −3. There is no strong evidence for a recent burst and a simple delayed SFH is adopted for all the fits of the stellar continuum.
We then compare the three attenuation scenarios. They have the same number of free parameters and ∆BIC is now defined as the difference between the best χ 2 of each law: ∆BIC = χ 2 (C00) − χ 2 (CF00 or LF17). The χ 2 are plotted in Fig. 2 with limits corresponding to |∆BIC| > 6. We clearly see that C00 gives the best fits. In all cases, C00 is acceptable with ∆BIC < 6 (∆BIC is even lower or equal to zero except in one case). We adopt the C00 dust attenuation scenario for all our galaxies. The best fits of the stellar continuum for the whole galaxy sample are shown in Appendix A.

Fitting the dust IR SED
We model the dust emission using all the IR data available (Spitzer/MIPS, Herschel/PACS and SPIRE, and ALMA) to measure L dust independently of the stellar continuum. This quantity is commonly used to measure obscured star formation using simple assumptions about the SFH and dust absorption (e.g., Kennicutt & Evans 2012). In this case no energy budget is performed and the IR data are fitted to the Draine & Li (2007) models. No AGN is introduced. Without IRAC data which are not considered for the fits, the constraint on the AGN component would be very weak. We plot the resulting values in Fig. 3 and compare them to the measurements of Elbaz et al. (2018) for the 15 galaxies in common. We find very good agreement between both values. The L dust values are listed in Table 3.

Fitting the full SED
We now fit the full datasets with the same models and parameters defined in Sect. 3. First, as for the stellar continuum, we test the A79, page 4 of 13 requirement to introduce a recent burst by calculating the BIC difference between the scenarios with and without a burst. We apply the same criterion: there is strong evidence for a burst when ∆BIC < −6. Using C00 and CF00 we find strong evidence for a burst for respectively 13 and 11 of the 17 galaxies. We find ∆BIC 6 (i.e., close to strong evidence for no burst) for UDF13 with both C00 and CF00 and also for UDF7 and GS4 with CF00. With LF17 we find strong evidence for a burst for only three galaxies (UDF3, GS3, and GS6). For all the other cases, ∆BIC 6. The requirement for a burst for C00 and CF00 is essentially due to the fact that the delayed SFH model alone is most of the time unable to produce enough IR emission. With LF17, the attenuation law is much grayer and the amount of attenuation is high enough to produce the dust emission observed without a burst, or with a burst of very low amplitude. There is no strong evidence in favor of a simple delayed SFH without a burst. As expected when ∆BIC 6, the burst fractions returned by the code are always very small and never higher than 1%. Therefore for the sake of simplicity we adopt a delayed SFH on which a recent burst is superimposed for all the targets and dust attenuation laws. The average (median) values of the burst fraction found for C00, CF00, and LF17 are respectively 0.18 (0.20), 0.10 (0.10), and 0.03 (10 −6 ).
We plot in Fig. 4 the best χ 2 for the three attenuation laws. This time, no model appears to be clearly better than the other ones for the entire sample. In Table 2 we report the best attenuation law for each target based on the best χ 2 .
The other acceptable laws on the basis of a BIC difference lower than six are also listed. The C00, CF00, and BLF7 models correspond to the best fit for seven, six, and four objects, and are acceptable for two, one, and eight other cases respectively. In Fig. 5 we show two examples (UDF4 and GS6) for which different laws are preferred 4 .
The AGN fractions (defined as the fraction of the total IR luminosity coming from the AGN) obtained for X-ray detected objects are small ( f AGN = 0.05 ± 0.07) and do not have a significant  impact on the results. In the following we use the best attenuation law for each target according to Table 2. The best fits for the full SED of the whole galaxy sample are shown in Appendix A.

Dust emission
L dust is an output of the three fits performed in Sect. 4. Here we compare the different estimates. CIGALE outputs L dust due to absorption of starlight, as the result of the energy budget. We correct the measure of dust emission with IR data (Sect. 4.2) for the AGN contribution estimated with the full SED in order  Examples of best fits of the full SED with the three attenuation laws for two galaxies: UDF4 (lower row) and GS6 (upper row). x-axis: observed wavelength in µm, y-axis: observed flux in mJy. From left to right: C00, CF00, and LF17 laws. The black line represents the best fitted spectrum composed of attenuated stellar emission, nebular lines, and dust emission. The data points are plotted with empty blue squares and their 3σ error with blue vertical lines. The intrinsic, unabsorbed stellar continuum is plotted with a dashed blue line, the red shaded area indicates the amount of absorbed light, reprocessed in IR. The best fit is found with C00 for UDF4 (χ 2 min = 42.7, ∆BIC = 23.3 and 29.81 for LF17 and CF00 respectively): the fit with C00 is better in the NIR which is underestimated with CF00 and LF17. Here, GS6 is best fitted with LF17 (χ 2 min = 54.06, ∆BIC = 17.8 and 117 for CF00 and C00 respectively) and C00 and, to a lesser extent CF00, produce an overly high attenuation in the UV, a grayer attenuation law as LF17 better reproduces the UV with a substantial attenuation in the NIR.
to compare the three measurements of L dust . While we expect a large difference between the results from the fits of the stellar continuum and the other measurements, the L dust values measured from IR data and from the full SED should be similar. We show in Fig. 6 that the agreement between both measurements is indeed very good.

Dust emission from the fit of the stellar continuum
We find a substantial attenuation when fitting the stellar continuum alone with an average value of A V = 1.6 ± 0.5 mag. In , and with the fit of the stellar continuum only (filled squares for GS galaxies and empty squares for UDF galaxies, y-axis). The 1σ errors given by CIGALE for the fit of the stellar continuum are plotted as vertical and horizontal lines. The blue crosses correspond to L dust values obtained with the slope of the UV continuum. The solid line represents equal quantities on both axes, and dashed and dotted lines represent a L dust respectively three and ten times lower on the y-axis.
the luminosities measured with IR data. The median value of the ratio between both measurements is 0.41 for the whole sample, and 0.56 and 0.28 for the UDF and GS galaxies respectively. In one case (UDF13) the IR emission is over-predicted from the fit of the stellar continuum alone, with an excess of 30%. In Sect. 4 we show that both the stellar and full SEDs of UDF13 can be fitted without adding a burst, for any attenuation law. No extra attenuation is needed to explain L dust measured from Herschel and ALMA data. Conversely, for two galaxies (UDF3 and GS6), less than 10% of the dust emission can be explained from the rest-frame UV to NIR fit. We found that a burst must be added to fit the full SED of these galaxies with any of our attenuation laws. We see in Sect. 6 that both are the most active galaxies of our sample in star formation. Both Dunlop et al. (2017) and Elbaz et al. (2018) also report a substantial attenuation from the analysis of the stellar continuum. This result is puzzling: in these objects the ALMA emission does not match the UV rest-frame emission and in some cases there is even a full disconnection between the two. Our fits are averaged over the whole of the galaxy in each case and are not necessarily valid for the specific regions where the UV emission is detected. The light distribution in the HST/ACS F160W filter (corresponding to visible rest frame) extends much farther than what is seen in the HST/ACS F606W and F814W filters, which correspond to rest-frame UV (Rujopakarn et al. 2016;Elbaz et al. 2018). These different spatial distributions may be of consequence for energy balance models even when only the stellar continuum is considered. One way to check the energetic link between the UV and dust emissions is to estimate L dust only from the rest-frame UV emission as we do in the following section. respectively, allowing us to estimate the slope of the UV continuum, which is expected to follow a power law ( f λ ∝ λ β ). We exclude the most distant galaxies of the sample from this analysis, namely UDF1 and UDF2. At z 2.7, the F435W filter samples the emission below 120 nm where the attenuation curve is not well constrained and is likely to depart from the Calzetti law (Leitherer et al. 2002;Buat et al. 2002;Reddy et al. 2016). We calculate β 5 using the fits of the stellar continuum alone after ensuring that the Bayesian estimates of the fluxes in the F435W and F606W filters are in very good agreement with the observed values (Fig. 7). With β at hand, we estimate the IRX ratio using the IRX-β relation valid for local starburst galaxies. We use the inner relation of Overzier et al. (2011), close to the original fit of Meurer et al. (1999). L dust is then deduced by multiplying IRX with the UV rest-frame luminosity estimated by CIGALE and corresponding to the GALEX/FUV filter.

Dust emission from the slope of the UV continuum
We plot the L dust obtained from IRX-β in Fig. 6 for comparison with the results obtained in the previous subsection. We find similar trends considering either the fit of the stellar continuum or using only β. Given the large uncertainties on the measure of β with only 2 bands and on the validity of the IRX-β relation for these objects, we do not consider the values based on the IRXβ relation as sufficiently reliable. We only use them confirm the rather large attenuation found for the UV emission, which corresponds to a substantial fraction of the total dust emission of these galaxies, at least in terms of energy. Gómez-Guijarro et al. (2018) also found that the dust emission is always associated with the reddest stellar component.
In the case of UDF13, L dust calculated from β is much higher than the total observed dust luminosity from the dust only model. Similarly, L dust measured from its stellar continuum is also larger than when measured from the dust emission only, although with a much lower excess. This galaxy is not starbursting as we showed before, and we cannot interpret its relatively red UV slope (β = 1.2) as being exclusively due to dust attenuation.

Comparison of the HUDF, GOODS-S, and
COSMOS/HELP surveys As part of the HELP project 6 , Małek et al. (2018) analyzed the SED of more than 40 000 galaxies of the pilot field ELAIS-N1 5 β = −6.9 × log( f F435W / f F606W ) − 2 with fluxes expressed in mJy. 6 The Herschel Extragalactic Legacy Project (https://herschel. sussex.ac.uk/). HELP data products are distributed on HeDaM (http://hedam.lam.fr/HELP/).   Fig. 8. Comparison between L dust estimated from the fit of the stellar continuum (x-axis) and with IR data (y-axis). Filled squares are used to represent GS galaxies and UDF galaxies are plotted with empty squares. The COSMOS sample is plotted with dots color coded with redshift. The black solid and dashed lines represent the linear regression and 2σ dispersion for the full COSMOS sample. The red dot-dashed and dotted line the 2 and 3σ dispersion for the COSMOS galaxies with 1.5 < z < 2.5 and L dust > 10 11 L . and performed SED fitting with CIGALE. They estimated L dust from the stellar continuum for 686 galaxies in the ELAIS-N1 field detected with Herschel with a signal-to-noise ratio (S/N) larger than two. They found a very good agreement and a tight relation between both estimations (their Fig. 12.c). The ELAIS-N1/HELP field is not very deep with a mean redshift of 0.93 and not suitable for a comparison with our current study. In order to compare samples in the same range of redshift and L dust , we consider the COSMOS/HELP field for which similar data products are available. As for the ELAIS-N1 study, we select galaxies with Herschel detections (at least two SPIRE fluxes with an S /N > 3 and at least one PACS flux with S /N > 2), at least six photometric measurements of the stellar continuum, and good-quality fits of the stellar continuum and IR SED (for more details of the selection see Małek et al. 2018). We select 3365 galaxies with z < 3 and 694 with 1.5 < z < 2.5. In Fig. 8 we plot the two estimations of L dust (from the stellar continuum alone and with IR data) for the COSMOS sample, the correlation is good with a correlation coefficient of 0.90. We can see on the plot that the dispersion around the linear fit increases with redshift and luminosity. The dispersion for the full sample is σ = 0.26 and increases to σ = 0.36 for the subsample of galaxies with 1.5 < z < 2.5 and L dust > 10 11 L . The values for the GS and UDF galaxies are over-plotted. All of them lie on or above the linear regression line found for the COSMOS sample. All but one of the UDF galaxies (UDF3) lie within the 2σ limit for 1.5 < z < 2.5 and L dust > 10 11 L . The GS galaxies are close to or above this 2σ limit but they remain below 3σ. This difference of behavior between UDF and GS objects is expected since HUDF is a blind survey whereas the sources of Elbaz et al. (2018) were selected to be very massive and IR bright. We return to this comparison in Sect. 6.3 after SFR and stellar masses estimations and in relation to the position of the galaxies in the MS. Nevertheless, a full comparison between ALMA-and Herschel-selected samples is beyond the scope of this paper, and would include a careful analysis of the selection bias for each sample. Multiwavelength analyses, including SED fitting, of larger ALMA blind surveys like the GOODS-ALMA and ASAGAO surveys (Franco et al. 2018;Hatsukade et al. 2018) Fig. 9. Comparison between stellar masses (M ) measured with the fit of the stellar continuum and C00 (x-axis) and the fit of the full SED with the best attenuation law (y-axis). The different attenuation laws adopted for the fit of the full SED are indicated with red (C00), blue (LF17), and black (CF00) filled circles. The solid line represents the 1:1 relation, the dashed (resp. dotted) line a stellar mass two (resp. three) times lower on the y-axis.
are also needed to perform a robust statistical comparison with the large Herschel surveys.

Star formation rate and stellar mass determinations
The main parameters extracted from the SED fitting are usually stellar masses and SFRs. We list in Table 3 the stellar masses obtained from the fit of the stellar continuum together with the SFR and stellar masses obtained using the full SED. We discuss their reliability here.

Stellar mass
In Fig. 9 we compare the stellar masses measured from the stellar continuum alone (with C00 and a simple delayed SFH) and the full SED (with the best attenuation law and a delayed SFH with a recent burst). As expected the choice of the attenuation law has a strong impact on the stellar mass.
The flatter attenuation laws are always favored for the most massive systems (>∼10 11 M ). There is a good agreement between the masses when C00 is used for both fits of the full SED and stellar continuum. In this case, we obtain only slightly lower values with the full SED (factor of 0.8 on average). We attribute this difference to the presence of a burst of young stars in the case of the full SED fitting, which over-shines older stellar populations (Pforr et al. 2012). We find higher stellar masses with CF00 and LF17 with an average systematic difference of close to a factor of two (reaching 2.8 for LF17 against 1.5 for CF00) in agreement with the findings of Małek et al. (2018). This is likely due to the substantial attenuation affecting restframe stellar emission and even the NIR with CF00 and LF17.
The stellar masses obtained by Dunlop et al. (2017) from their fit of the UV-to-NIR SED of UDF galaxies agree within 0.15 dex with our estimations from the stellar continuum. This agreement is expected since both fits are performed assuming C00.

Star formation rate
In Fig. 10, we compare the SFR measured from the full SED with the simple approach of adding the SFR obtained from the total dust and UV emissions: SFR(IR)+SFR(UV).  Fig. 10. Comparison between the SFRs measured from SED fitting (y-axis) and the sum of SFR(IR) and SFR(UV), (x-axis). The filled circles correspond to the fit of the full SED as in Fig. 9. The empty circles correspond to the SFRs obtained with the fit of the stellar continuum alone. The solid line represents equal quantities on both axes; the dashed and dotted lines represent a quantity respectively 3 and 10 times lower on the y-axis.
We adopt here the calibrations of Buat et al. (2008) between luminosities and SFR. We use L dust values obtained with the fit of the IR data and corrected for the AGN contribution as explained in Sect. 5. The SFRs(IR) agree within a factor two with the same quantity measured by Dunlop et al. (2017) for the UDF galaxies (SFR FIR1 in their Table 4) except for UDF3 for which our measurement is 2.5 higher than the value they reported. The restframe UV luminosity (without dust correction) is an output of CIGALE already used in Sect. 5.2: it comes from the fit of the stellar continuum alone.
We find a satisfying agreement between the two estimates of the SFRs. The SFR(IR)+SFR(UV) calculation is based on very simple assumptions on the SFH and dust heating. We therefore do not expect a perfect agreement with our measurements coming from SED fitting. In Fig. 10 we also plot the SFRs measured from the stellar continuum alone. We find they are much lower than the ones based on the full SEDs by a factor similar to the one found for the estimation of L dust .
Once again, UDF13 departs from the general trend. Its SFR measured from the full SED is lower than both the addition of the IR and UV SFRs and the SFR measured from the stellar continuum.The difference with the former comes from the SFH used for the fit (a delayed SFH with a peak 1 Gyr before and no recent burst with a burst fraction equal to 10 −7 ) versus the simple assumption of a constant SFH over 100 Myr to calculate SFR(IR)+SFR(UV) (e.g., Kennicutt & Evans 2012;Buat et al. 2008). The discrepancy with the latter is explained by higher L dust inferred from the stellar continuum alone.

Position on the main sequence of star forming galaxies
We define the MS with the relations of Schreiber et al. (2015) as a function of stellar mass and redshift. For each source we consider the two stellar masses previously computed either from the stellar continuum or from the full SED and we calculate the two SFRs corresponding to these stellar masses for a galaxy on the MS. In Fig. 11, for each target, we compare these two SFRs with the SFR measured previously from the fit of the full SED and plot the ratio of these values as a function of the ratio of the two stellar masses. We assume a dispersion of 0.3 dex for the MS (Schreiber et al. 2015). With the stellar masses measured with the stellar continuum, six galaxies are found above the MS and one below. If we define a starburst as a galaxy with an SFR exceeding the MS by a factor at least four, only one galaxy is found to be starbursting (UDF3). Another one (GS6) has an SFR 3.5 times higher than the MS. If we now use the stellar masses measured from the full SED, five galaxies are found above the MS but close to the boundaries of the relation, and no galaxy is starbursting. UDF13 is still well below the MS and UDF3 very close to starbursting. Dunlop et al. (2017) concluded that their data were completely consistent with the published MS relations at z 2. Elbaz et al. (2018) also introduced the variation of the MS with redshift by using the relations of Schreiber et al. (2015) as we did above. They also found that the UDF and GS galaxies are located in the upper part of the MS and only UDF3, GS6, and GS7 are starbursting, with a SFR more than four times above the MS. In our analysis, GS6 and GS7 appear less extreme but still well above the MS for a stellar mass estimated with the stellar continuum.
The difference found between L dust estimations with the fit of either IR data or stellar continuum discussed in Sect. 5.3 is found to increase with star forming activity. In Fig. 12 the same quantities as in Fig. 8 are plotted but this time color coded with the specific SFR defined as the SFR from the fit of the full SED divided by the stellar mass from the fit of the stellar continuum. For COSMOS, UDF, and GS galaxies the excess of dust emission measured with IR data increases with sSFR. The fact that all our UDF and GS galaxies are found in the upper part of the MS explains why they depart from the mean trend found by Małek et al. (2018) between both estimations of L dust for the ELAIS-N1 field and confirmed here in the COSMOS field, but all remaining with the 3σ limit. Comparison between L dust estimated from the fit of the stellar continuum (x-axis) and with IR data (y-axis) as in Fig. 8. Diamonds are used to represent GS galaxies and UDF galaxies are plotted with squares. The COSMOS sample is plotted with dots. The symbols are color coded with sSFR.

Dust attenuation laws and compactness
7.1. ALMA, H-band radii, and attenuation laws The compactness of the most intense sources detected by ALMA leads us to compare the effective radii measured in the H and ALMA bands respectively. We use the circularized half-light radii of Elbaz et al. (2018) and we define their ratio: R H ALMA = R(H)/R(ALMA). Rujopakarn et al. (2016) also measured radii corresponding to the star formation and stellar mass distributions for the UDF galaxies. Although these radii correspond to the distribution of physical quantities that we are discussing in this paper, we prefer to adopt the Elbaz et al. (2018) radii values to have homogeneous measurements for both GS and UDF galaxies. The radii measured by Rujopakarn et al. (2016) for the UDF galaxies correspond to slightly smaller radius ratios, that is, more similar star formation and stellar mass distributions; using them would not change our conclusions. We calculate R H ALMA for 15 of our 17 galaxies (UDF10 and UDF15 are not included in the Elbaz et al. 2018 sample and are not measured by Rujopakarn et al. 2016). In Fig. 13 we split the sample in two parts: the galaxies which can be fitted with C00 (from Table 2 with C00 in Col. 4 or 5) and the galaxies which are only well fitted with a power-law curve (CF00 or LF17). We plot the histogram of R H ALMA for the two subsamples. All the galaxies which can be fitted with C00 correspond to R H ALMA < 2. This is the case of all the UDF galaxies except for UDF3, for which R H ALMA = 1 and the best law is LF17. Six out of the seven galaxies whose SEDs are only well fitted with a power-law curve (CF00 or LF17) exhibit R H ALMA ≥ 1.6. All the GS galaxies are in this configuration except GS1, with R H ALMA = 0.9 and a best fit with C00. Our current study is limited by the small number of sources. The need for flatter attenuation laws for sources with a very compact dust distribution has to be confirmed with a larger sample of sources with measured morphological parameters from optical and ALMA images (e.g., Fujimoto et al. 2018, and references therein). This is planned for a future study.

Validity of attenuation laws
The galaxies of our sample with the largest difference between dust and stellar distributions correspond to a power-law attenuation curve, which leads to a substantial attenuation up to  Fig. 13. Histogram of the ratio of the H-band radius to the ALMA radius, R H ALMA . The filled histogram represents galaxies for which C00 gives a satisfying fit, and the empty histogram represents galaxies for which a shallower power-law attenuation (CF00 or LF17) is needed. the rest-frame NIR. This attenuation is particularly large when LF17 is used (UDF3, GS3, GS5, and GS6, as seen in Fig. A.1). The better fits obtained with a shallower attenuation law reflect the need to produce a large dust emission with a substantial to strong attenuation of the whole stellar continuum. Variable attenuation laws from the UV to the NIR are found from dust radiative transfer calculations, either for simplified models of galactic disks or coupled with hydrodynamical simulations and a substantial to high optical depth experienced by all the stellar populations leads to a flattening of the effective attenuation law. (Chevallard et al. 2013;Trayford et al. 2019). Roebuck et al. (2019) performed radiation transfer on a library of hydrodynamic simulations for isolated disk and mergers. They found flat attenuation curves for the dusty phases of each simulation. These curves are well reproduced by the LF17 law and are found consistent with a heavily obscured young stellar population with a contribution of an unattenuated older stellar population, with a similar stellar mass in both components.
Variations of the effective UV attenuation curve are reported by several studies based on observations and the IRX-β scatter plot is recognized as a diagnostic for these variations (e.g., Salmon et al. 2016;Lo Faro et al. 2017;Popping et al. 2017;Narayanan et al. 2018;Salim & Boquien 2019). A flattening at longer wavelengths than the UV is much more difficult to identify on the observed SED because of the well-known degeneracy between age and attenuation. The analysis of the UV to submm SED puts some constraints on this flattening from an energy balance approach, as we did in this study and in Lo Faro et al. (2017). As mentioned above, radiation transfer models link the flattening of the effective attenuation law to an increase of the effective optical depth at the corresponding wavelength. It is not clear if this is really valid for our targets with R H ALMA > 2 as their rest-frame visible stellar light extends well beyond the ALMA emission. In this case energy balance SED modeling may over-estimate the stellar mass by adding extra attenuation up to the NIR rest frame if a substantial part of this emission comes from unattenuated old stars.
When the stellar and dust size distributions are more similar, the C00 law gives satisfying fits. The stellar masses from the full SED and the stellar continuum only are similar. We note that we cannot exclude a moderate flattening in the NIR: the steep decrease of C00 for λ > 0.5 µm (Fig. 1) is not easy to reproduce, even when a single stellar population is assumed (Seon & Draine 2016;Buat et al. 2019). In these galaxies fitted with C00, the ALMA and HST distributions can also be very different at small scales with a smooth cold dust distribution and a patchy/disturbed stellar emission as reported by Rujopakarn et al. (2019) with high ALMA spatial resolution imaging of three UDF galaxies. Nevertheless their dust emission is fitted with the "classical" C00.
The case of UDF3 is not consistent with the general trend. Whereas the ALMA and H-band radii are similar, LF17 fits this galaxy better. In this case, a gray attenuation curve is likely to be representative of the stars and dust distributions inside the galaxy and the stellar mass estimated from the fit of the stellar SED with C00 underestimates the stellar mass by a factor 2.6. We note however that this is our worst fit with any attenuation law considered in this work (Table 2 and Fig. A.1) Our analysis shows that some indirect information can be retrieved from the UV to submm (i.e stellar and dust) energy balance modeling. According to Fig. 13, when we can fit the full SED with C00, a similar radius is found for the stellar population in the NIR and the cold dust detected by ALMA. Conversely, for galaxies with a very compact dust distribution, some additional obscuration is needed to produce enough dust emission. In this work we model it with a grayer attenuation law. This modification of the attenuation law seems to reflect heavily obscured objects with different dust and observed stellar distributions, as also identified in the nearby universe (Charmandaris et al. 2004;Howell et al. 2010) and in various studies at higher redshift based on ALMA data (e.g., Hodge et al. 2016;Gómez-Guijarro et al. 2018). The ULIRGs studied by Lo Faro et al. (2017) and the galaxies of the ELAIS-N1 field studied by Małek et al. (2018) and fitted with flat attenuation laws could also be in this configuration.

Conclusions
Here, we analyze the SEDs of 17 dust-rich galaxies at z ∼ 2 detected with ALMA from the blind survey conducted by Dunlop et al. (2017) in the HUDF and the massive Herschel galaxies in the GOODS-South field by Elbaz et al. (2018). We performed several fits with CIGALE with different attenuation laws, considering only the stellar continuum, the dust emission, and the full dataset.
The fit of the stellar continuum with the starburst law of Calzetti et al. (2000) reveals an obscuration which corresponds on average to one-third (UDF galaxies) to half (GS galaxies) of the total dust emission detected with Herschel and ALMA. Qualitatively, the analysis of the slope of the UV continuum confirms this obscuration of the stellar continuum. We compared our estimations of dust emission to the same quantities measured for a sample of galaxies of similar L dust and redshift in the COSMOS field. The deficit of L dust found from the fit of the stellar continuum is on average absent for the COSMOS galaxies, but is consistent within 2σ and 3σ respectively for the UDF and GS galaxies. In all samples, UDF, GS, and COSMOS, the deficit is found to strongly correlate with sSFR.
The fit of the full SED (stellar and dust components) with CIGALE reveals that different attenuation laws have to be considered. The Calzetti et al. (2000) law is valid for most of the UDF galaxies. Shallower power-law attenuations (Charlot & Fall 2000;Lo Faro et al. 2017) are preferred for most of the GS galaxies. A requirement for flatter attenuation laws is found for galaxies with a very compact dust distribution detected by ALMA, much less extended than the stellar emission observed by HST. The stellar mass estimates are dependent on the choice of the attenuation law, with differences reaching a factor of approximately three when either the stellar continuum or the full SED is used for the fit, assuming different attenuation laws. The SFRs are robustly estimated as long as the dust emission is considered. We find that the SFR from the fit of the full SED are consistent with the sum of the SFR from the IR and UV.  Table 2, Col. 4. The dashed blue line represents the spectrum corresponding to the best fit of the stellar continuum (Sect. 4.1) extending to IR and submm to show the dust emission predicted. The AGN component, if any, is plotted with a thin green line. The intrinsic, unabsorbed, stellar continuum of each best fit are plotted with a thin blue line, the shaded areas indicate the corresponding amount of attenuation (in red for the fit of the full SED, in blue for the fit of the stellar continuum). The data points are plotted with empty blue squares and their 3σ error with blue vertical lines.