# ALMA Lensing Cluster Survey: Dust mass measurements as a function of redshift, stellar mass, and star formation rate from $z = 1$ to $z = 5$

Jean-Baptiste Jolly, Kirsten Knudsen, Nicolas Laporte, Andrea Guerrero, Seiji Fujimoto, Kotaro Kohno, Vasily Kokorev, Claudia del P. Lagos,

Thiébaut-Antoine Schirmer, Franz Bauer, et al.

► To cite this version:

Jean-Baptiste Jolly, Kirsten Knudsen, Nicolas Laporte, Andrea Guerrero, Seiji Fujimoto, et al. ALMA Lensing Cluster Survey: Dust mass measurements as a function of redshift, stellar mass, and star formation rate from $z = 1$ to $z = 5$. Astronomy & Astrophysics - A&A, 2025, 693, pp.A190. 10.1051/0004-6361/202346239. [hal-05469649](https://hal.science/hal-05469649v1)

HAL Id: hal-05469649

https://hal.science/hal-05469649v1

Submitted on 21 Jan 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# ALMA Lensing Cluster Survey: Dust mass measurements as a function of redshift, stellar mass, and star formation rate from $z = 1$ to $z = 5$

Jean-Baptiste Jolly<sup>1,2,*</sup>, Kirsten Knudsen<sup>2</sup>, Nicolas Laporte<sup>3</sup>, Andrea Guerrero<sup>4</sup>, Seiji Fujimoto<sup>5,6</sup>, Kotaro Kohno<sup>7,8</sup>, Vasily Kokorev<sup>9</sup>, Claudia del P. Lagos<sup>10</sup>, Thiébaut-Antoine Schirmer<sup>2</sup>, Franz Bauer<sup>11,12,13</sup>, Miroslava Dessauge-Zavadsky<sup>14</sup>, Daniel Espada<sup>15,16</sup>, Bunyo Hatsukade<sup>7</sup>, Anton M. Koekemoer<sup>17</sup>, Johan Richard<sup>18</sup>, Fengwu Sun<sup>19</sup>, and John F. Wu<sup>20,21</sup>

Received 24 February 2023 / Accepted 15 November 2024

## ABSTRACT

**Context.** Understanding the dust content of galaxies, its evolution with redshift and its relation to stars and star formation is fundamental for our understanding of galaxy evolution. Dust acts as a catalyst of star formation and as a shield for star light. Advanced millimeter facilities like ALMA have made dust observation ever more accessible, even at high redshift. However, dust emission is typically very faint, making the use of stacking techniques is instrumental in the study of dust in statistically sound samples.

**Aims.** Using the ALMA Lensing Cluster Survey (ALCS) wide-area band-6 continuum dataset (~110 arcmin<fnref n="2" /> across 33 lensing clusters), we constrain the dust-mass evolution with redshift, stellar mass, and star formation rate (SFR).

**Methods.** After binning sources according to redshift, SFR, and stellar mass as extracted from an HST-IRAC catalog, we performed a set of continuum-stacking analyses in the image domain using LINESTACKER on sources between $z = 1$ and $z = 5$, which further improved the depth of our data. The large field of view provided by the ALCS allowed us to reach a final sample of ~4000 galaxies with known coordinates and SED-derived physical parameters. We stacked sources with an SFR between $10^{-3}$ and $10^3$ M☉ per year and a stellar mass between $10^8$ and $10^{12}$ M☉, and we split them into different stellar mass and SFR bins. Through stacking, we retrieved the continuum 1.2 mm flux, which is a known dust-mass tracer. This allowed us to derive the dust-mass evolution with redshift and its relation to the SFR and stellar mass.

**Results.** We clearly detect the continuum in most of the subsamples. From the nondetections, we derive $3\sigma$ upper limits. We observe a steady decline in the average dust mass with redshift. Moreover, sources with a higher stellar mass or SFR have a higher dust mass on average. This allows us to derive scaling relations. Our results mostly agree well with models at $z \sim 1-3$, but they indicate a typically lower dust mass than predicted at higher redshift.

**Key words.** dust, extinction – galaxies: evolution – galaxies: ISM – galaxies: statistics

<footnote n="1">¹ Max-Planck-Institut für extraterrestrische Physik, 85748 Garching, Germany</footnote>

<footnote n="2">² Department of Space, Earth and Environment, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden</footnote>

<footnote n="3">³ Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, 13388 Marseille, France</footnote>

<footnote n="4">⁴ Astronomy Department, Universidad de Concepción, Barrio Universitario S/N, Concepción, Chile</footnote>

<footnote n="5">⁵ Cosmic Dawn Center (DAWN), Copenhagen, Denmark</footnote>

<footnote n="6">⁶ Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark</footnote>

<footnote n="7">⁷ Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan</footnote>

<footnote n="8">⁸ Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="9">⁹ Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="10">¹⁰ International Centre for Radio Astronomy Research (ICRAR), M468, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia</footnote>

<footnote n="11">¹¹ Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Campus San Joaquín, Av. Vicuña Mackenna 4860, 7820436 Macul Santiago, Chile</footnote>

<footnote n="12">¹² Centro de Astroingeniería, Facultad de Física, Pontificia Universidad Católica de Chile, Campus San Joaquín, Av. Vicuña Mackenna 4860, 7820436 Macul Santiago, Chile</footnote>

<footnote n="13">¹³ Millennium Institute of Astrophysics, Nuncio Monseñor Sótero Sanz 100, Of 104, Providencia, Santiago, Chile</footnote>

<footnote n="14">¹⁴ Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland</footnote>

<footnote n="15">¹⁵ Departamento de Física Teórica y del Cosmos, Campus de Fuentenueva, Edificio Mecenas, Universidad de Granada, E-18071 Granada, Spain</footnote>

<footnote n="16">¹⁶ Instituto Carlos I de Física Teórica y Computacional, Facultad de Ciencias, E-18071 Granada, Spain</footnote>

<footnote n="17">¹⁷ Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="18">¹⁸ CRAL, Observatoire de Lyon, Université Lyon 1, 9 Avenue Ch. André, F-69561 Saint Genis Laval Cedex, France</footnote>

<footnote n="19">¹⁹ Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson 85721, USA</footnote>

<footnote n="20">²⁰ Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="21">²¹ Department of Physics & Astronomy, Johns Hopkins University, 3400 N Charles St, Baltimore, MD 21218, USA</footnote>

<footnote>* Corresponding author; jbjolly@mpe.mpg.de</footnote>
# 1. Introduction

Dust impacts the direct evolution of galaxies and their observations in multiple ways. It is thought to be the main catalyst of H₂ formation (Wakelam et al. 2017), which is one of the main components of molecular clouds in the interstellar medium (ISM), and hence, it is one of the principal drivers of star formation (e.g. Scoville 2013). In addition to the direct effect on the physical properties of galaxies, dust also plays a critical role in astrophysical observations. The dust grains attenuating the starlight are consequently heated, and they in turn reradiate the energy at longer wavelengths. This is the so-called dust-obscured star formation (see for e.g. Casey et al. 2014; Hodge & da Cunha 2020; Zavala et al. 2021). As a consequence, the light from star-forming regions is typically dominated by dust continuum emission. The relation between dust mass and other physical characteristics of galaxies have been intensely studied, for example, through the dust-to-gas (atomic and molecular) or dust-to-metal ratios (e.g. Hunt et al. 2005; Draine et al. 2007; Engelbracht et al. 2008; Galametz et al. 2011; Magrini et al. 2011; Saintonge et al. 2013; Rémy-Ruyer et al. 2014; Combes 2018; Li et al. 2019; Shapley et al. 2020; Tacconi et al. 2020; Birkin et al. 2021; Popping & Péroux 2022; Popping et al. 2023), or through the ratio of the dust mass versus star formation rate (SFR) (e.g. da Cunha et al. 2010; Casey 2012; Santini et al. 2014; Dudzevičiūtė et al. 2020, 2021).

Recent studies have shown that the proportion of dust-obscured star formation has evolved with redshift (Bouwens et al. 2020). It has been estimated to be a dominant fraction during the peak of the cosmic star formation history. However, it appears to have accounted for only ~20%–25% of the total star formation at redshift $z = 6–7$ (see Zavala et al. 2021). Similarly, the overall cosmic dust density has been shown to peak around $z \sim 1–2$ and to rapidly decline at higher redshift (see for example Driver et al. 2018; Magnelli et al. 2020; Pozzi et al. 2020; Dudzevičiūtė et al. 2021). Furthermore, massive dust reservoirs have been detected in massive high-redshift galaxies (e.g. Bertoldi et al. 2003; Valiante et al. 2009; Venemans et al. 2012; Watson et al. 2015; Laporte et al. 2017; Tamura et al. 2019, 2023). While this is probably not the case for the more typical less massive high-redshift galaxies, it highlights the importance of studying the evolution of dust mass with redshift.

Millimeter (mm) or submillimeter (submm) emission was proposed to be a tracer of dust mass (Scoville et al. 2014, 2016, 2017), with the assumption that the emission is optically thin and measured far from the peak of the spectral energy distribution of the dust. The Atacama Large Millimeter/Submillimeter Array (ALMA) has become instrumental in the quest to study the dust mass in lower-mass galaxies and its evolution in the mm and submm bands. With large surveys and high-redshift single-target observations becoming progressively more accessible (e.g. Knudsen et al. 2016; Walter et al. 2016; González-López et al. 2017, 2020; Laporte et al. 2017; Béthermin et al. 2020; Aravena et al. 2020), the evolution of dust is studied in increasingly statistically sound samples.

However, observations of high-redshift galaxies are typically biased by construction toward the observation of the brightest sources (Malmquist bias). To draw a complete picture of the dust evolution with redshift, it is necessary to also study galaxies with lower intrinsic luminosities. Gravitational lensing can be used as a tool to enhance the signal from faint galaxies

without the need for excessive integration time. Similarly, tools such as stacking can statistically drastically improve the signal-to-noise ratio (S/N) for large samples. With the help of gravitational lensing and stacking, it is hence possible to push the limit of observations toward sources with a lower intrinsic dust luminosity.

We present the stacking analysis of 10386 gravitationally lensed galaxies from the ALMA Lensing Cluster Survey (ALCS) at $z > 1$ (cluster and field sources at $z \le 1$ are studied in a separate paper; Guerrero et al. (2023). By binning galaxies by redshift and further splitting them according to their stellar mass or SFR, we study the evolution of the dust mass with redshift and its scaling relation with stellar mass and SFR. We also integrate the contribution from all galaxies in each redshift bin to assess the evolution of the cosmic dust density.

In Sect. 2 we describe the overall data set, the catalog, and the different subsamples. In Sect. 3 we describe the stacking method, which we performed using LINESTACKER (Jolly et al. 2020), and the processes involved in the dust-mass calculation. In Sect. 4 we present our results, and we discuss them in Sect. 5. Finally, in Sect. 6 we summarize this paper. In the appendix (available on Zenodo), we show alternative stacking procedures as well as some additional details for the main analysis.

Throughout the paper, we assume a ΛCDM cosmology with $\Omega_m = 0.3$, $\Omega_\lambda = 0.7$ and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-3}$. All magnitudes are quoted in the AB system, such that $M_{AB} = 23.9–2.5 \log_{10}(S_v [\mu Jy])$.

# 2. Data and sample

## 2.1. ALCS

The ALCS is a large ALMA program that was accepted in cycle 6 (Project ID: 2018.1.00035.L; PI: K. Kohno). It observed 33 lensing clusters in band 6 ($\lambda \sim 1.2$ mm), which spanned ~110 arcmin² (primary beam (PB) > 0.5). The clusters are distributed as follows: 16 clusters were taken from RELICS (the Reionization Lensing Cluster Survey; Coe et al. 2019), 12 clusters were taken from CLASH (the Cluster Lensing And Supernova survey with Hubble; Postman et al. 2012), and 5 clusters were taken from the Frontier Fields survey (Lotz et al. 2017). The observations were carried out between December 2018 and December 2019 (cycles 6 and 7) in compact-array configurations C43-1 and C43-2 in a double-frequency window setup. The observations were made from 250.0 to 257.5 GHz and from 265.0 to 272.5 GHz for a total bandwidth of 15 GHz. When available, the ALCS data were concatenated with existing ALMA data, notably the ALMA Frontier Fields Survey (Project ID: 2013.1.00999.S, PI: Bauer and Project ID; 2015.1.01425.S; PI: Bauer). The data were reduced and calibrated using the Common Astronomy Software Applications (CASA, McMullin et al. 2007) package version 5.4.0 for the 26 clusters observed in cycle 6 and with v5.6.1 for the remaining clusters that were observed in cycle 7. Throughout this paper, we use natural weighted primary-beam-corrected and uv-tapered continuum maps, with a tapering parameter of 2 arcsec (the full width at half maximum of the synthesized beam is 2 arcsec; with a corresponding pixel size of 0.16 arcsec). uv-tapered maps were chosen over natural resolution maps to ensure the beam-size homogeneity of the different images when stacking. The average root mean square (RMS) of the maps is ~63 µJy/beam (see Table 1 for the detailed RMS of each map). A full description of the survey can be found in Kohno et al. (2023).
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Distance from the main sequence $\Delta(\text{MS}) = \log(\text{SFR}_{\text{MS}})$ – $\log(\text{SFR})$ (from Speagle et al. 2014) as a function of redshift for the galaxies in the full sample. The points are color-coded according to stellar mass. The dashed line at $\Delta(\text{MS}) = -0.5$ highlights the region from which quiescent galaxies were excluded (see the appendix, which is available on Zenodo).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Comparison between the SMF of the COSMOS2020 sample (Weaver et al. 2023, 2022) and the SMF of the sources studied in this paper. The SMF was obtained by counting the galaxies in each redshift bin separately for each cluster and dividing by the corresponding volume, corrected for the mean magnification of the sources in each bin. The error bars correspond to the 16th to 84th percentile of the stellar mass of the sources in the sample, and the impact of the 20% error on the magnification on the total volume.</figcaption></figure>

### 2.3. Subsamples

To study the evolution of the dust mass across redshift, we separated the sources into four redshift bins: $1 \le z < 2$, $2 \le z < 3$, $3 \le z < 4$, and $4 \le z < 5$. We additionally split sources according to either their SFRs or stellar masses into the following bins: (i) SFR (upper and lower limits in $M_{\odot} \text{yr}^{-1}$): $0.001 \le \text{SFR} \le 1$, $1 < \text{SFR} \le 10$, $10 < \text{SFR} \le 50$, $50 < \text{SFR} \le 100$, and $100 < \text{SFR} \le 1000$, and (ii) stellar mass ($M_{*}$, upper and lower limits in $M_{\odot}$): $10^8 \le M_{*} \le 10^9$, $10^9 < M_{*} \le 10^{10}$, $10^{10} < M_{*} \le 10^{11}$, and $10^{11} < M_{*} \le 10^{12}$. The total number of sources in each subsample is listed in Table 2. Furthermore, the distribution of the SFR and stellar mass in each subsample can be found in the appendix (available on Zenodo).

<table><caption>**Table 2.** Number of stacked sources in each subsample.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Redshift range</th></tr><tr><th>1-2</th><th>2-3</th><th>3-4</th><th>4-5</th><th>Total</th></tr></thead><tbody><tr><td>8 &le; log(M<sub>*</sub>) &le; 9</td><td>1004</td><td>521</td><td>432</td><td>375</td><td>2332</td></tr><tr><td>9 &lt; log(M<sub>*</sub>) &le; 10</td><td>508</td><td>295</td><td>309</td><td>256</td><td>1368</td></tr><tr><td>10 &lt; log(M<sub>*</sub>) &le; 11</td><td>154</td><td>43</td><td>84</td><td>60</td><td>341</td></tr><tr><td>11 &lt; log(M<sub>*</sub>) &le; 12</td><td>27</td><td>4</td><td>11</td><td>20</td><td>62</td></tr><tr><td>0.001 &le; SFR &le; 1</td><td>1046</td><td>441</td><td>306</td><td>250</td><td>2043</td></tr><tr><td>1 &lt; SFR &le; 10</td><td>529</td><td>377</td><td>455</td><td>394</td><td>1755</td></tr><tr><td>10 &lt; SFR &le; 50</td><td>113</td><td>43</td><td>67</td><td>58</td><td>281</td></tr><tr><td>50 &lt; SFR &le; 100</td><td>5</td><td>1</td><td>1</td><td>5</td><td>12</td></tr><tr><td>100 &lt; SFR &le; 1000</td><td>-</td><td>1</td><td>7</td><td>4</td><td>12</td></tr></tbody></table>

Sources that were individually detected in the ALCS data were not excluded from the samples. We decided to distinguish solely based on the properties mentioned above to derive the average properties of the population, and the 1.2 mm continuum flux density of individual sources should therefore not be a criterion of exclusion from our samples. To assess the impact of the 1.2 mm flux distribution in our samples, we performed median stacking and bootstrapping analyses, alongside the main mean-stacking analyses (see Sects. 3.4 and 4.2).

## 3. Methods

### 3.1. Stacking

We performed the stacking with LINESTACKER (Jolly et al. 2020) in single-channel mode (for the continuum stacking) on uv-tapered images (see Sect. 2.1). Sources were stacked pixel to pixel using mean stacks without weights (median stacks were also performed; see Section 4.2). The stacking was operated with stamp sizes of $9.76 \times 9.76 \text{ arcsec}^2$ (61 $\times$ 61 pixels). The sources were spatially aligned using the position extracted from the Kokorev et al. (2022) catalog. The stacked flux in Jy was retrieved by integrating the continuum flux in a central circular region of the stack stamp, with a radius $R_{\text{integ}} = 2''$ (12.5 pixels). This corresponds to the synthesized beam size. We then divided it by the beam size in pixel units.

To compute the RMS associated with each stacked cube, a random set of source-free coordinates was drawn for each stacking position. An empty stack was then generated from the set of random coordinates (see Jolly et al. 2020) with the same characteristics as the normal stacks (i.e., stamp size and number of targets). This process was performed 1000 times for each subsample. The standard deviation across all stack stamps was then computed and used to derive the RMS associated with each stack cube. When the flux in the central region of the stack stamp was lower than three times its associated RMS, an upper limit of $3\sigma$ was used in place of the integrated flux to compute $M_{\text{dust}}$.

### 3.2. Magnification corrections

The averaged 1.2 mm fluxes extracted from each stack maps were corrected for the mean magnification of the sources in the stack subsample. The final flux was therefore computed as $F_{\text{final}} = \frac{F_{\text{stack}}}{\bar{\mu}}$, where $F_{\text{final}}$ is the final flux we used to compute the dust mass (see Sect. 3.3), $F_{\text{stack}}$ is the flux obtained from the stacked map, and $\bar{\mu}$ is the average magnification of each source in the subsample. Alternatively, each source could be corrected for
its magnification before stacking, $F_{final} = F_{corrected stack}$, where $F_{corrected stack}$ is the flux obtained from the stack of the individually magnification-corrected maps, $Map_{stack} = \sum_{i=1}^{n} \frac{Stamp_i}{\mu_i}$, where $Map_{stack}$ is the stacked map obtained from individually magnification-corrected stamps (from which $F_{corrected stack}$ is obtained), $n$ is the number of sources in the subsample, $Stamp_i$ is the map-stamp associated with source $i$, and $\mu_i$ is the magnification associated with source $i$. We decided to correct for the average magnification of the sample after stacking because the signal in each stack stamp consists of the flux from the source, which is gravitationally magnified, and some noise, which is independent of the magnification. By correcting each stack stamp for its associated magnification, we corrected the flux from the source and the foreground noise. The efficiency of stacking is based on the fact that the noise is random and approximately similar in every stack stamp. By downscaling the noise in each stamp by a different number (unrelated to the noise properties), we might effectively reduce the efficiency and reliability of the stacking analysis. However, by correcting the final stack with the average magnification, the analysis would be biased toward high-magnification sources. Correcting each stamp individually pre-stacking should avoid this bias, and we therefore repeated our analysis with this correction method. The 1.2 mm fluxes derived in both cases are very consistent, and the values overlap within the error margins.

### 3.3. Measuring $M_{dust}$

To derive dust masses using the 1.2 mm continuum flux extracted from each stacking stamp, we used a single modified blackbody curve under the approximation of an optically thin regime. Following Kovács et al. (2010) and Magnelli et al. (2020), we defined $M_{dust}$ as

$$M_{\text{dust}} = \frac{5.03 \times 10^{-31} (S_{\text{vobs}} / f_{\text{CMB}}) D_{\text{L}}^2}{(1+z)^4 B_{\text{vobs}} (T_{\text{obs}}) \kappa_{\text{v0}}} \left(\frac{\nu_0}{\nu_{\text{rest}}}\right)^{\beta}, \quad (1)$$

where $M_{dust}$ is the dust mass in $M_{\odot}$, $S_{vobs}$ is the flux at the observed frequency in Jy (and corrected for the average magnification, as explained in Sect. 3.2), $\nu_{obs} = \nu_{rest}/(1+z)$ is the observed frequency, $f_{CMB}$ is the correction factor to account for the cosmic microwave background (CMB; see details further in the text), $D_L$ is the luminosity distance in meter at redshift $z$, $B_{vobs}(T_{obs})$ is Planck's blackbody function in Jy sr$^{-1}$ at the observed temperature $T_{obs}$ ($T_{obs} = T_{rest}/(1+z)$), $\kappa_{v0}$ is the photon cross section to mass ratio of dust in m$^2$ kg$^{-1}$, and $\beta$ is the dust emissivity index. $S_{vobs}$ was obtained by integrating over a fixed circular region in the center of the stack stamp (see Section 3.1).

Following Magnelli et al. (2020) and Pozzi et al. (2021), we decided to use a single cold component to account for the dust temperature, and we used a mass-weighted dust temperature of $T_{rest} = 25$ K. While the total dust in galaxies is thought to consist of a warmer ($20 < T < 60$ K) and a colder ($T < 30$ K) component, studies of local galaxies have shown that the cold component causes most of the dust budget (from 96% to 99%, see Orellana et al. 2017) and the majority of the Rayleigh-Jeans emission (see also Shivaei et al. 2022, for a study of the contribution of the cold and warm dust components to the dust continuum emission at $z \sim 2$). Scoville et al. (2014, 2016, 2017) argued that this approximation is probably still valid in higher-redshift galaxies. We decided to adopt $T_{rest} = 25$ K at all redshifts probed by our analyses, but we note that dust masses evolve as $T^{-1}$, meaning that our dust-mass results are effectively

highly dependent on the estimated dust temperature, and that the masses derived in our analyses might be off by factors of a few if the assumed dust temperature is not correct (see Sect. 5).

While it might seem more appropriate to choose a dust temperature that evolves with redshift, we decided to maintain a fixed value to facilitate the comparison with similar previous studies. We used $\beta = 1.8$, the Galactic dust emissivity index measurement from the Planck data (Planck Collaboration XXI 2011), which also correlates well with values observed in high-redshift galaxies (e.g. Chapin et al. 2009; Magnelli et al. 2012; Faisst et al. 2020). Theoretical studies predicted values for $\beta$ that range from 1.5 to 2.0 (e.g. Draine 2011), and similarly, observational studies such as Faisst et al. (2020) reported $\beta$ values between 1.6 and 2.4 (with a median value of 2.0). The effect of $\beta$ is weaker than the choice of the dust temperature, however, and we found that choosing $\beta = 1.5$ (or 2.0) would only modify our results by ~9% at $z = 1.5$ and ~15% at $z = 5.5$. We used $\kappa_{v0} = 0.0431 \text{ m}^2 \text{ kg}^{-1}$ with $\nu_0 = 352.6 \text{ GHz}$ (Li & Draine 2001; Magnelli et al. 2020). We refer to Sect. 5 for a more complete discussion of the choice of parameters for the dust-mass computation.

Following da Cunha et al. (2013), we corrected our flux measurements for the effect of the CMB. First, we corrected our dust temperature for the additional heating through the CMB (following Equation (12) of da Cunha et al. 2013). This effect is minor because the dust temperature is only increased by ~3% at $z = 6$ (the effect being even less important at lower redshift). Second, and most importantly, the CMB acts as a bright observing background, leading to an underestimation of the total flux. This was corrected for by following Equation (18) of da Cunha et al. (2013). This effect is much more important, as it yields an upward measurement to our results that ranges from ~1.03 at $z \sim 1$ to ~1.31 at $z \sim 6$.

### 3.4. Computing the uncertainty

To compute the overall uncertainties associated with our measurements, we combined different sources of uncertainties. The first source was the direct RMS from each stack analysis (computed from source-free stacks of each subsample; see Sect. 3.1). In addition, the intrinsic flux distribution of the sources in the sample needs to be accounted for in stacking analyses. To do this, we performed a bootstrapping analysis for each subsample. In each subsample, the dust mass was recomputed 1000 times, each time with a different sample, randomized from the original subsample without replacement (see Jolly et al. (2020) for a detailed description of the bootstrapping routine included in LINESTACKER). The distribution of the dust masses computed in this manner is shown in the appendix (available on Zenodo), where we also show the dust mass obtained from the original subsample (vertical black line). The peaks and shapes of the distributions are overall consistent with the values derived from the original samples. However, the peak of the distributions in most of the subsamples is slightly lower than in the original stacks. Similarly, some distributions present a faint tail toward higher dust masses. These two pieces of information combined indicate a skewed distribution, which can easily be explained because the samples contain individual detections. These effects are small, however, and should be well represented by the associated uncertainties on the dust-mass measurements.

From the distributions derived from each bootstrapping analyses, we extracted the 16th to 84th percentile (see the appendix available on Zenodo), which we summed quadratically
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** 1.2 mm continuum stack maps. (left panel) $9.76 \times 9.76$ arcsec$^2$ (61 × 61 pixels) mean stacking stamps, split into bins of stellar mass and redshift. Each map is normalized by the corresponding standard deviation, computed in associated empty stacks (see Sect. 3.1). The number of sources stacked (N) is indicated for each bin. (right panel) Similar to the left panel, but using SFR bins instead of stellar mass ones.</figcaption></figure>

with the stack RMS<fnref n="1" />. In addition, because the redshift is used multiple times in equation 1, its 16th to 84th percentile in each subsample was used to propagate the associated uncertainty on the dust-mass computation. Finally, following Sun et al. (2022) and Fujimoto et al. (2023), we used a magnification uncertainty of 20% of the magnification associated with each subsample, and we propagated it to the aforementioned errors.

The 16th and 84th percentile of the (magnification-corrected) SFRs and stellar masses in each subsample were used to compute the error associated with the SFRs and stellar masses.

The combination of these measurement uncertainties was used to plot the error bars in the different figures shown in this work. However, only the RMS associated with each cube was used to qualify a detection or an upper limit (as stated in Sect. 3.1).

# 4. Results

## 4.1. Mean stacking results and dust continuum detection

The standard deviation in stack stamps reaches levels as low as ~4.3 × 10⁻³ mJy per beam in the $z = 1–2$ bin of the lowest stellar mass and SFR subsamples (both containing ~1000 sources). The total flux in the central circular aperture with a radius of 2 arcsec varies between ~1.78 × 10⁻² mJy and ~1.04 mJy for stack detections (flux above 3σ, not corrected for magnification). Securely extracted dust masses (with an S/N above 3σ) range from ~2.8 × 10⁶ M☉ to ~2.56 × 10⁸ M☉ (after correction for magnification). Figure 4 shows the stacked maps we obtained for each of the subsamples, from which we extracted the results presented below.

The evolution of the dust mass with stellar mass and redshift is shown in Figure 5 and tabulated in Table 3. While detections in the lowest-redshift bin ($1 < z < 2$) indicate a clear linear relation between log($M_{\text{dust}}$) and log($M_*$), the other redshift bins are less regular (see the left panel of Figure 5). When compared to other works (Santini et al. 2014; da Cunha et al. 2015; Kirkpatrick et al. 2017; Shivaei et al. 2022, also plotted on the figures), our data points are typically below other measurements (individual detections and stacks). This behavior is strongly reduced when quiescent galaxies are excluded from the stacks, however (see the appendix, which is available on Zenodo). This highlights that quiescent galaxies might have lower dust masses at the same stellar mass.

The log($M_{\text{dust}}$) − log(SFR) relation seems to follow a similar trend of increasing dust mass with increasing SFR, but with a clearer linear relation across redshift (although there is a higher number of nondetections; see the left panel of Figure 6 and the stacking results tabulated in Table 4), which indicates that the SFR might be better suited than the stellar mass to trace the dust mass. When compared to previous works (Santini et al. 2014; da Cunha et al. 2015; Kirkpatrick et al. 2017; Shivaei et al. 2022), our measurements follow a similar trend as was already observed, but they reach so far mostly unexplored regimes.

The evolution of the dust mass with redshift (see the right panels of Figures 5 and 6 and Figure 7) shows a relatively clear linear trend that indicates a decrease in the dust mass with increasing redshift.

By fitting quadratic functions to the data shown in Figures 5 and 6, we derived scaling relations between $M_{\text{dust}}$ and $M_*$ at fixed $z$; $M_{\text{dust}}$ and $z$ at fixed $M_*$; $M_{\text{dust}}$ and SFR at fixed $z$; and $M_{\text{dust}}$ and $z$ at fixed SFR. The fits were performed using detections only (ignoring upper limits), and only when the number of data points was ≥ 3. The results from the fits are summarized in Table 5, and they are plotted in the appendix (available on Zenodo).

<footnote n="1">¹ The bootstrap analysis is also affected by the RMS of the stacked images. The two errors are therefore correlated and should in principle not be added quadratically. This is an approximation that probably slightly overestimates the error computation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Average dust mass evolution with stellar mass and redshift. (left) Average dust mass as a function of the stellar mass in each redshift bin. The circles represent detections (above $3\sigma$), and down-pointing arrows represent $3\sigma$ upper limits. The error bars on the dust-mass detections were computed using a combination of different errors (see Sect. 3.4). The error bar along the x-axis represents the 16th to 84th percentile of parameter distribution in the stacked sample. (right) Similar to the left panel, but plotted as a function of redshift in each stellar mass bin. We also show for comparison individual dust-mass measurements detections (da Cunha et al. 2015; Kirkpatrick et al. 2017) and stacking measurements from Santini et al. (2014), Shivaei et al. (2022).</figcaption></figure>

<table><caption>**Table 3.** Stacking results in the stellar mass subsamples.</caption><thead><tr><th>⟨z⟩</th><th>log(⟨|M<sub>*</sub>|⟩/M<sub>☉</sub>)</th><th>⟨Dust mass⟩<br>(M<sub>☉</sub>)</th><th>Flux<br>(mJy)</th><th>RMS<br>(mJy)</th><th>⟨μ⟩<sup>(a)</sup></th></tr></thead><tbody><tr><td>1.42</td><td>8.48</td><td>&lt;2.49 × 10<sup>6</sup></td><td>—</td><td>4.28 × 10<sup>-3</sup></td><td>3.43</td></tr><tr><td>2.51</td><td>8.53</td><td>&lt;1.82 × 10<sup>6</sup></td><td>—</td><td>5.64 × 10<sup>-3</sup></td><td>5.41</td></tr><tr><td>3.45</td><td>8.54</td><td>&lt;1.27 × 10<sup>6</sup></td><td>—</td><td>6.51 × 10<sup>-3</sup></td><td>7.83</td></tr><tr><td>4.47</td><td>8.56</td><td>&lt;9.43 × 10<sup>5</sup></td><td>—</td><td>6.98 × 10<sup>-3</sup></td><td>10.55</td></tr><tr><td>1.42</td><td>9.43</td><td>1.61 × 10<sup>7</sup> ± 4.33 × 10<sup>6</sup></td><td>7.84 × 10<sup>-2</sup></td><td>6.67 × 10<sup>-3</sup></td><td>3.24</td></tr><tr><td>2.54</td><td>9.39</td><td>3.89 × 10<sup>6</sup> ± 2.54 × 10<sup>6</sup></td><td>2.54 × 10<sup>-2</sup></td><td>8.25 × 10<sup>-3</sup></td><td>3.78</td></tr><tr><td>3.41</td><td>9.41</td><td>4.96 × 10<sup>6</sup> ± 2.61 × 10<sup>6</sup></td><td>4.49 × 10<sup>-2</sup></td><td>9.12 × 10<sup>-3</sup></td><td>4.64</td></tr><tr><td>4.44</td><td>9.42</td><td>&lt;2.19 × 10<sup>6</sup></td><td>—</td><td>9.12 × 10<sup>-3</sup></td><td>5.95</td></tr><tr><td>1.41</td><td>10.42</td><td>6.95 × 10<sup>7</sup> ± 1.84 × 10<sup>7</sup></td><td>0.28</td><td>1.27 × 10<sup>-2</sup></td><td>2.67</td></tr><tr><td>2.46</td><td>10.33</td><td>5.57 × 10<sup>7</sup> ± 1.78 × 10<sup>7</sup></td><td>0.3</td><td>2.37 × 10<sup>-2</sup></td><td>3.16</td></tr><tr><td>3.43</td><td>10.35</td><td>3.57 × 10<sup>7</sup> ± 1.32 × 10<sup>7</sup></td><td>0.33</td><td>1.96 × 10<sup>-2</sup></td><td>4.71</td></tr><tr><td>4.46</td><td>10.37</td><td>1.7 × 10<sup>7</sup> ± 1.01 × 10<sup>7</sup></td><td>0.2</td><td>2.02 × 10<sup>-2</sup></td><td>5.67</td></tr><tr><td>1.41</td><td>11.2</td><td>1.02 × 10<sup>8</sup> ± 4.17 × 10<sup>7</sup></td><td>0.39</td><td>2.88 × 10<sup>-2</sup></td><td>2.53</td></tr><tr><td>2.56</td><td>11.22</td><td>&lt;3.58 × 10<sup>7</sup></td><td>—</td><td>9.51 × 10<sup>-2</sup></td><td>4.6</td></tr><tr><td>3.55</td><td>11.38</td><td>&lt;3.21 × 10<sup>7</sup></td><td>—</td><td>4.93 × 10<sup>-2</sup></td><td>2.33</td></tr><tr><td>4.46</td><td>11.39</td><td>1.55 × 10<sup>7</sup> ± 1.11 × 10<sup>7</sup></td><td>0.16</td><td>4.03 × 10<sup>-2</sup></td><td>4.87</td></tr></tbody></table>

**Notes.** When the central flux is lower than three times the RMS of the corresponding stack map, the $3\sigma$ upper limit is shown in place of the dust mass. The errors on the computed dust mass correspond to the quadratic sum of the standard deviation computed from empty stacks with the uncertainties from the bootstrap analysis, to which we finally propagated the magnification error (see Sect. 3.4). Flux and RMS are extracted as indicated in Sect. 3.1: The flux is integrated in the central circular region of the stack stamp, with a radius $R_{\text{integ}} = 2''$, and it is then converted into Jy from Jy/beam. The RMS is computed from the standard deviation in the empty stacks. $^{(a)}$⟨μ⟩ is the average magnification of the sources in the subsample.

## 4.2. Median stacking analysis

To better assess the effect of the population distribution on the derived dust mass, we also performed median stacks in the same way as the average stacks presented in Sect. 3.1. The overall median dust masses derived through our stacking analyses are systematically slightly lower than the mean dust masses (see Figure 8 for an illustration and the appendix available on Zenodo for detailed plots). While the median stellar mass and SFR of each subsample are also systematically lower than the mean (see the figures in the appendix), the difference is not large enough to explain the offset dust mass observed between the mean and

median stacks. This again indicates a possibly skewed distribution of the dust mass in each subsample that results in offsets between the mean and median dust masses. However, we note that the dust-mass trends observed in this analysis (with redshift, stellar mass, and SFR) remain the same in the median stacking analyses.

## 5. Discussion

As our analysis focused on galaxies spanning a wide redshift range and our observation remained at a fixed observed wavelength, it is important to question the assumption made here
<figure><img src="image_6.png" /><figcaption>Fig. 6. Similar to Figure 5, but in SFR bins instead of stellar mass.</figcaption></figure>

<table><caption>Table 4. Similar to Table 3, but for the SFR subsamples.</caption><thead><tr><th>⟨z⟩</th><th>⟨SFR⟩<br>(M<sub>☉</sub>/yr)</th><th>⟨Dust mass⟩<br>(M<sub>☉</sub>)</th><th>Flux<br>(mJy)</th><th>RMS<br>(mJy)</th><th>⟨μ⟩</th></tr></thead><tbody><tr><td>1.39</td><td>0.15</td><td>3.72 × 10<sup>6</sup> ± 1.82 × 10<sup>6</sup></td><td>1.95 × 10<sup>-2</sup></td><td>4.28 × 10<sup>-3</sup></td><td>3.51</td></tr><tr><td>2.51</td><td>0.81</td><td>&lt;1.94 × 10<sup>6</sup></td><td>—</td><td>6.17 × 10<sup>-3</sup></td><td>5.55</td></tr><tr><td>3.45</td><td>0.99</td><td>&lt;1.14 × 10<sup>6</sup></td><td>—</td><td>7.17 × 10<sup>-3</sup></td><td>9.62</td></tr><tr><td>4.47</td><td>0.53</td><td>&lt;6.72 × 10<sup>5</sup></td><td>—</td><td>7.05 × 10<sup>-3</sup></td><td>14.94</td></tr><tr><td>1.45</td><td>2.84</td><td>1.78 × 10<sup>7</sup> ± 4.96 × 10<sup>6</sup></td><td>8.13 × 10<sup>-2</sup></td><td>6.47 × 10<sup>-3</sup></td><td>3.04</td></tr><tr><td>2.53</td><td>1.79</td><td>&lt;3.17 × 10<sup>6</sup></td><td>—</td><td>7.19 × 10<sup>-3</sup></td><td>3.95</td></tr><tr><td>3.43</td><td>2.16</td><td>6.43 × 10<sup>6</sup> ± 2.94 × 10<sup>6</sup></td><td>5.51 × 10<sup>-2</sup></td><td>7.36 × 10<sup>-3</sup></td><td>4.38</td></tr><tr><td>4.45</td><td>4.68</td><td>&lt;2.21 × 10<sup>6</sup></td><td>—</td><td>7.72 × 10<sup>-3</sup></td><td>4.97</td></tr><tr><td>1.47</td><td>24.56</td><td>7.06 × 10<sup>7</sup> ± 1.83 × 10<sup>7</sup></td><td>0.26</td><td>1.38 × 10<sup>-2</sup></td><td>2.5</td></tr><tr><td>2.52</td><td>19.66</td><td>3.4 × 10<sup>7</sup> ± 1.23 × 10<sup>7</sup></td><td>0.2</td><td>2.29 × 10<sup>-2</sup></td><td>3.49</td></tr><tr><td>3.4</td><td>31.67</td><td>2.64 × 10<sup>7</sup> ± 1.05 × 10<sup>7</sup></td><td>0.22</td><td>2.37 × 10<sup>-2</sup></td><td>4.23</td></tr><tr><td>4.46</td><td>10.55</td><td>2.72 × 10<sup>7</sup> ± 1.78 × 10<sup>7</sup></td><td>0.19</td><td>2.29 × 10<sup>-2</sup></td><td>3.4</td></tr><tr><td>1.76</td><td>86.71</td><td>2.56 × 10<sup>8</sup> ± 1.26 × 10<sup>8</sup></td><td>1.04</td><td>7.88 × 10<sup>-2</sup></td><td>2.64</td></tr><tr><td>2.09</td><td>60.65</td><td>1.8 × 10<sup>8</sup> ± 4.46 × 10<sup>7</sup></td><td>0.65</td><td>9.54 × 10<sup>-2</sup></td><td>2.25</td></tr><tr><td>3.87</td><td>81.75</td><td>&lt;6.56 × 10<sup>7</sup></td><td>—</td><td>0.1</td><td>2.35</td></tr><tr><td>4.58</td><td>53.34</td><td>&lt;3.55 × 10<sup>7</sup></td><td>—</td><td>6.11 × 10<sup>-2</sup></td><td>2.45</td></tr><tr><td>2.84</td><td>1.33 × 10<sup>2</sup></td><td>2.54 × 10<sup>8</sup> ± 7.27 × 10<sup>7</sup></td><td>1.0</td><td>0.2</td><td>2.17</td></tr><tr><td>3.5</td><td>1.51 × 10<sup>2</sup></td><td>&lt;4.71 × 10<sup>7</sup></td><td>—</td><td>6.67 × 10<sup>-2</sup></td><td>2.16</td></tr><tr><td>4.44</td><td>3.52 × 10<sup>2</sup></td><td>&lt;4.89 × 10<sup>7</sup></td><td>—</td><td>0.1</td><td>2.93</td></tr></tbody></table>

that the continuum flux traces the same physical process at a rest wavelength $\lambda_{\text{rest}} \sim 0.5$ mm (corresponding to $z \sim 1$) and $\lambda_{\text{rest}} \sim 0.18$ mm (corresponding to $z \sim 5$). More specifically, the validity of the dust-mass equation (see Sect. 3.3) comes from the assumption that the dust is optically thin and that the flux measurement probes the Rayleigh-Jeans dust emission. This means that the part of the SED that is dominated by the emission of the dust at $T = 25$ K dominates the dust-mass budget. This is true at long wavelengths, where we effectively probe the Rayleigh-Jeans tail of the dust emission (e.g. Magdis et al. 2012; Scoville et al. 2016). However, the higher the redshift observed, the closer the region to the peak of the dust emission, which is more sensitive to temperature and total luminosity. This might explain the low detection rate in the high-$z$ subsamples, even when the number of stack sources is high, and for which the corresponding RMS in the stack maps is low.

While numerous studies showed an evolution of the dust temperature with redshift, it is important to distinguish between the peak temperature, which is thought from theoretical and observational works to increase with redshift (e.g. Magdis et al. 2012; Magnelli et al. 2014; Béthermin et al. 2015; Ferrara et al. 2017; Narayanan et al. 2018; Schreiber et al. 2018; Liang et al. 2019; Ma et al. 2019; Sommovigo et al. 2020), and the mass-weighted temperature that probes the Rayleigh-Jeans tail, which can be directly used to measure the dust mass, as stated above. Unlike the peak temperature, the mass-weighted temperature is thought to be mostly constant with redshift (see Scoville et al. 2016; Liang et al. 2019). It should be noted, however, that the mass-weighted temperature still exhibits a small range in possible values, typically from 15 to 45 K (e.g. Liang et al. 2019), and a small increase with redshift may exist. On the other hand, Dudzevičiūtė et al. (2020) argued that the previously observed
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Ratio of the average dust mass to average stellar mass as a function of redshift. Similar to the right panel of Figure 5, but normalized by the stellar mass of the galaxies in the sample. The error bars combine the error on the dust-mass computation and stellar mass distribution.</figcaption></figure>

evolution of the dust temperature is likely due to the luminosity evolution in the samples employed, and they suggested a dust temperature of 30.4 ± 0.3, constant with redshift. In any event, to facilitate comparison with studies similar to ours and to limit complexity, we decided to keep a fixed temperature at all redshifts for the dust-mass computations. Consequently, for a T = 40 K (instead of 25) in the dust-mass computation, our results would change by ~50% in the highest-redshift bin.

When comparing our results to the semi-analytical modeling from Popping et al. (2017), Lagos et al. (2019), Vijayan et al. (2019); Triani et al. (2020, see Figures 9 and 10)<fnref n="2" /> we can see an overall decent agreement in dust-mass predictions as a function of stellar mass at z ~ 1. However, all models but the one presented in Vijayan et al. (2019) predict an overall increase in the dust mass with redshift while our analysis shows the opposite. In addition, the dust mass derived in the high-mass end of our subsamples is typically lower than predicted by the models (although this flattening might be due to the presence of quiescent galaxies in our sample, as shown in the appendix available on Zenodo). In addition, the method that is used to compute the average dust mass in the SHARK model (Lagos et al. 2019) differs from the one used in this work. To highlight this difference, we show in Figures 9 and 10 the average dust masses predicted by the SHARK model and the dust masses computed using equation 1 from the band 6 continuum (S6) predicted by SHARK. It is interesting to note that at z ~ 3 and at z ~ 5, the dust mass derived using the predicted S6 differs by almost one dex from the dust mass predicted otherwise. As mentioned above, this might be a direct consequence of the evolution of the rest wavelength with redshift, which moves the observation closer to the peak of the dust emission. Surprisingly, however, our results seem to agree better with the direct dust-mass measurements compared to those obtained from S6.

As the analysis relies on photometric redshifts for a large fraction of the catalog, it is important to consider the potential impact this might have on the results. In particular, we note that the bins at z > 3 appear to be overpopulated. This is shown in Figure 3, where we compare the SMF of our sample to the SMF of the COSMOS2020 sample (Weaver et al. 2023, 2022). As a consequence, the dust content at z > 3 is likely overestimated.

This is probably due to the misclassification of some lower-redshift sources to higher redshift, as suggested in Kokorev et al. (2022) and mentioned in Sect. 2.2.

Moreover, because the catalog we used was generated using *HST* and *Spitzer* photometry, some of the most dusty galaxies may be missing. This biases the analysis toward lower averaged dust masses. As pointed out in Kokorev et al. (2022), only 145 of the 180 sources that were individually detected with an S/N > 4 in the ALCS data are identified in their catalog. This indicates that some of the most dusty galaxies may be missing from our analysis. This bias could be strongest in the high-redshift bins, where only the UV-brightest galaxies (i.e., typically less obscured by dust) might be observed. The use of catalogs based on data from the James Webb Space Telescope (JWST) might help us to include more of the missing dusty galaxies, especially at high z, and to improve the SED-fitting routines.

As shown in the appendix (available on Zenodo), the inclusion of quiescent galaxies in the stacks biases the analysis toward lower dust masses, especially in the samples that are split according to stellar mass. This highlights that some sources that are included in the high-mass subsamples have a low SFR and also low dust masses. The higher average dust mass observed in the mass-selected samples without sources with $\Delta$MS < -0.5 and the relatively stable average dust mass in SFR-selected samples implies that the SFR probably is a better dust-mass tracer. It shows in addition that quiescent galaxies typically have a lower dust content than their star-forming counterparts at similar stellar mass.

As in Magnelli et al. (2020), we used the dust opacities from Li & Draine (2001) to derive the dust masses. These values were constrained from observations of the diffuse interstellar medium (ISM), and these properties are therefore representative of dust in the diffuse ISM. However, some studies in nearby galaxies (e.g. Galliano et al. 2018) showed that the opacity from Li & Draine (2001) might have been underestimated by a factor of 2–3. Higher opacity values in galaxies do not mean that the values of Li & Draine (2001) are incorrect. Their values were constrained for the diffuse interstellar medium and are therefore well adapted for galaxies only if the bulk of the dust emission emerges from the diffuse interstellar medium of galaxies. The underestimation of the dust opacity might indicate that the majority of the dust emission does not emerge from diffuse regions, but from denser regions in which dust grains are likely to be larger.

In dense regions of the ISM, dust grains grow through accretion and coagulation, which triggers an increase in the dust opacity by a factor 2–3 (e.g. Köhler et al. 2015). The underestimation of the dust opacities in some nearby galaxies might therefore mean that the bulk of the dust emission arises from dense regions in which dust grains are larger than in the diffuse ISM. This seems to be the case in high-redshift galaxies as well, as Magnelli et al. (2020) claimed that the bulk of the dust emission in their galaxy sample emerged from the molecular phase, and therefore, from dense regions in which dust differs from the diffuse ISM. They also stated that there might be a significant increase in the dust emissivity from diffuse to denser regions, which can be explained in terms of grain growth (e.g. Köhler et al. 2015).

There is currently no consensus on the dust properties that must be considered in nearby or in high-redshift galaxies. Moreover, dust grains are probably not the same within galaxies and from one galaxy to the next, which adds to the complexity. Although diffuse ISM-like dust grains might not cause the majority of the dust emission in galaxies, these dust models have the advantage of being thoroughly developed and extensively used. We therefore use these dust properties from the diffuse

<footnote n="2">² The SHARK model (Lagos et al. 2019) does not directly track the buildup and destruction of dust. Instead, dust masses are estimated assuming an empirical z = 0 relation between the dust-to-gas ratio and gas-phase metallicity.</footnote>
<table><caption>**Table 5.** Fitting results.</caption><thead><tr><th colspan="4">M<sub>dust</sub>(SFR) (fixed z)</th></tr><tr><th>z range</th><th>a</th><th>b</th><th>Range of applicability</th></tr></thead><tbody><tr><td>1 ≤ z &lt; 2</td><td>7.2 ± 1.2 × 10<sup>6</sup></td><td>7.8 ± 0.7 × 10<sup>-1</sup></td><td>-1 &lt; log(SFR) &lt; 2</td></tr><tr><td>2 ≤ z &lt; 3</td><td>9.2 ± 5.7 × 10<sup>5</sup></td><td>12.3 ± 1.5 × 10<sup>-1</sup></td><td>1 &lt; log(SFR) &lt; 3</td></tr><tr><td>M<sub>dust</sub>(z) (fixed SFR)</td><td></td><td></td><td></td></tr><tr><th>SFR range (M<sub>☉</sub> per year)</th><th>a</th><th>b</th><th>Range of applicability</th></tr><tr><td>10 &lt; SFR &lt; 50</td><td>10.5 ± 1.6 × 10<sup>7</sup></td><td>-1.1 ± 0.2</td><td>1 ≤ z &lt; 5</td></tr><tr><td>M<sub>dust</sub>(M<sub>*</sub>) (fixed z)</td><td></td><td></td><td></td></tr><tr><th>z range</th><th>a</th><th>b</th><th>Range of applicability</th></tr><tr><td>1 ≤ z &lt; 2</td><td>5.2 ± 13.0 × 10<sup>3</sup></td><td>4.8 ± 1.1 × 10<sup>-1</sup></td><td>8 &lt; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 12</td></tr><tr><td>M<sub>dust</sub>(z) (fixed M<sub>*</sub>)</td><td></td><td></td><td></td></tr><tr><th>M<sub>*</sub> range (M<sub>☉</sub>)</th><th>a</th><th>b</th><th>Range of applicability</th></tr><tr><td>10<sup>9</sup> &lt; M<sub>*</sub> &lt; 10<sup>10</sup></td><td>3.0 ± 1.3 × 10<sup>7</sup></td><td>-1.8 ± 0.7</td><td>1 ≤ z &lt; 4</td></tr><tr><td>10<sup>10</sup> &lt; M<sub>*</sub> &lt; 10<sup>11</sup></td><td>10.5 ± 2.7 × 10<sup>7</sup></td><td>-1.0 ± 0.3</td><td>1 ≤ z &lt; 4</td></tr></tbody></table>

**Notes.** Fits are performed with a function of the form $M_{dust}(x) = ax^b$, where *a* and *b* are free parameters, and *x* is the SFR, *z*, or *M*<sub>*</sub>. The range of applicability refers to the parameter range over which subsamples yielded detections (see the appendix available on Zenodo to visualize the fits).

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Median stacking stamps, split into stellar mass and redshift (left panel) and SFR and redshift (right panel). Similar to Figure 4.</figcaption></figure>

ISM until we have more constraints on the dust grains that cause the emission in these galaxies.

## 6. Summary

We used deep ALMA band 6 data over 33 lensing clusters to perform a set of stacking analyses on a large sample of 4103 lensed sources, whose positions and physical parameters were extracted from the HST/IRAC catalog presented in Kokorev et al. (2022). The sample spanned redshifts from $z = 1$ to $z = 5$, SFRs from $0.001 \text{ M}_\odot \text{ yr}^{-1}$ to $1000 \text{ M}_\odot \text{ yr}^{-1}$, and stellar masses from $10^8 \text{ M}_\odot$ to $10^{12} \text{ M}_\odot$. To study the evolution of the dust mass through cosmic time, we binned the sources into five different redshift bins, and we also grouped sources according to their SFR and stellar mass. We performed a continuum stacking on each of the subsamples using LINESTACKER (Jolly et al. 2020).

Using the mean continuum-stacked flux, we computed average dust masses. From the nondetections, we derived $3\sigma$ upper

limits. Using detections and upper limits, we studied the evolution trend of the dust mass with redshift. We found clear indications for an average decrease in the dust mass with redshift. Similarly, we studied the evolution of the dust mass with stellar mass and SFR. In both cases, we found a positive correlation. From our detections, we derived the relations between the dust mass and SFR or stellar mass in different redshift ranges. The results from our study are mostly consistent with the results from modeling (at least at $z \sim 1$) and from other similar studies. Our analysis allowed us to probe regimes of the stellar mass and SFR that were unexplored so far. The highest-redshift bin shows mostly nondetections, and the corresponding upper limits indicate low average dust masses (even though the overall RMS is low) when compared to other studies. This could be due to the different rest-wavelength probed at high $z$, lying closer to the peak of the dust emission. Alternatively, the reason might be the tendency of individual high-$z$ dust measurements to be biased toward very dust-bright objects.
<figure><img src="image_9.png" /></figure>

**Fig. 9.** Log average dust mass recovered as a function of the (log) stellar mass in each redshift bin. The circles represent detections (above $3\sigma$), and the $3\sigma$ upper limits are represented by down-pointing arrows. We overplot the $z=1$ and $z=6$ dust mass-stellar mass relation from Popping et al. (2017) and the $z \sim 1$ and $z \sim 5$ relations from Lagos et al. (2019) (SHARK model), Triani et al. (2020), and Vijayan et al. (2019). We show the average dust mass directly predicted by the SHARK model and the dust mass computed from the predicted band 6 flux following Eq. (1).

<figure><img src="image_10.png" /></figure>

**Fig. 10.** Similar to Figure 9, but for the SFR instead of the stellar mass. The predictions from Lagos et al. (2019, SHARK model) at $z \sim 1, 3$, and 5 are shown in terms of direct dust-mass predictions and dust mass inferred from predicted band 6 flux, following Eq. (1).

The linear trends observed between stellar mass and dust mass confirm that as galaxies evolve and form more stars, they also accumulate more dust. These dust grains play a vital role in catalyzing the formation of H$_2$. In this way, they affect the overall gas reservoir that is available for future star formation. Dust becomes an integral component of a self-perpetuating cycle, accumulating throughout the stellar life cycle and subsequently bolstering the SFR. In addition, the inverse trend observed between the average dust mass and redshift implies that the buildup of dust in galaxies over cosmic time is a gradual process that is aligned with the overall evolution of the SFR density (e.g. Magnelli et al. 2020). This further indicates a gradual increase in the average metallicity in galaxies with cosmic time, as dust grains are known to be composed of heavy elements that are produced in stars.

## Data availability

Appendices are available on https://zenodo.org/records/14224076

**Acknowledgements.** The authors thank the anonymous reviewer for their thoughtful and constructive suggestions, which greatly helped improve this manuscript. This paper makes use of the ALMA data: ALMA #2018.1.00035.L, #2013.1.00999.S, and #2015.1.01425.S. ALMA is a partnership of the ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, AUI/NRAO, and NAOJ. JBJ thanks Ian Smail for the discussion. KK acknowledges support from the Swedish Research Council (2015-05580), and the Knut and Alice Wallenberg Foundation (KAW 2020.0081). K. Kohno acknowledges the JSPS KAKENHI Grant Numbers JP17H06130, JP22H04939, JP23K20035, JP24H00004, and the NAOJ ALMA Scientific Research Grant Number 2017-06B. AG acknowledges funding from ANID-Chile NCN2023_002 and FONDECYT Regular 1171506. DE acknowledges support from a Beatriz Galindo senior fellowship (BG20/00224) from the Spanish Ministry of Science and Innovation, projects PID2020-114414GB-100 and PID2020-113689GB-100 financed by MCIN/AEI/10.13039/501100011033, project P20_00334 financed by the Junta de Andalucía, project A-FQM-510-UGR20 of the FEDER/Junta de Andalucía-Consejería de Transformación Económica, Industria, Conocimiento y Universidades. FEB acknowledges support from ANID-Chile BASAL CATA FB210003, FONDECYT Regular 1200495 and 1190818, and Millennium Science Initiative Program – ICN12_009.

## References

Aravena, M., Boogaard, L., González-López, J., et al. 2020, ApJ, 901, 79

Bertoldi, F., Carilli, C. L., Cox, P., et al. 2003, A&A, 406, L55

Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113

Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

Birkin, J. E., Weiss, A., Wardlow, J. L., et al. 2021, MNRAS, 501, 3926

Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112

Casey, C. M. 2012, MNRAS, 425, 3094

Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45

Chapin, E. L., Pope, A., Scott, D., et al. 2009, MNRAS, 398, 1793

Coe, D., Salmon, B., Bradač, M., et al. 2019, ApJ, 884, 85

Combes, F. 2018, A&ARv, 26, 5

da Cunha, E., Eminian, C., Charlot, S., & Blaizot, J. 2010, MNRAS, 403, 1894

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2021, MNRAS, 500, 942

Engelbracht, C. W., Rieke, G. H., Gordon, K. D., et al. 2008, ApJ, 678, 804

Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020, MNRAS, 498, 4192

Ferrara, A., Hirashita, H., Ouchi, M., & Fujimoto, S. 2017, MNRAS, 471, 5018

Fujimoto, S., Kohno, K., Ouchi, M., et al. 2023, arXiv e-prints [arXiv:2303.01658]

Galametz, M., Madden, S. C., Galliano, F., et al. 2011, A&A, 532, A56

Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673

González-López, J., Bauer, F. E., Romero-Cañizales, C., et al. 2017, A&A, 597, A41

González-López, J., Novak, M., Decarli, R., et al. 2020, ApJ, 897, 91

Guerrero, A., Nagar, N., Kohno, K., et al. 2023, MNRAS, 526, 2423

Hodge, J. A., & da Cunha, E. 2020, Roy. Soc. Open Sci., 7, 200556

Hunt, L., Bianchi, S., & Maiolino, R. 2005, A&A, 434, 849

Jolly, J.-B., Knudsen, K. K., & Stanley, F. 2020, MNRAS, 499, 3992

Kirkpatrick, A., Pope, A., Sajina, A., et al. 2017, ApJ, 843, 71

Knudsen, K. K., Richard, J., Kneib, J.-P., et al. 2016, MNRAS, 462, L6

Köhler, M., Ysard, N., & Jones, A. P. 2015, A&A, 579, A15

Kohno, K., Fujimoto, S., Tsujita, A., et al. 2023, in Physics and Chemistry of Star Formation: The Dynamical ISM Across Time and Spatial Scales, eds. V. Ossenkopf-Okada, R. Schaaf, I. Breloy, & J. Stutzki, 16

Kokorev, V., Brammer, G., Fujimoto, S., et al. 2022, ApJS, 263, 38

Kovács, A., Omont, A., Beelen, A., et al. 2010, ApJ, 717, 29

Lagos, C. d. P., Robotham, A. S. G., Trayford, J. W., et al. 2019, MNRAS, 489, 4196

Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ, 837, L21

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425

Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97

Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6
Magnelli, B., Lutz, D., Santini, P., et al. 2012, A&A, 539, A155

Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86

Magnelli, B., Boogaard, L., Decarli, R., et al. 2020, ApJ, 892, 66

Magrini, L., Bianchi, S., Corbelli, E., et al. 2011, A&A, 535, A13

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127

Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718

Orellana, G., Nagar, N. M., Elbaz, D., et al. 2017, A&A, 602, A68

Planck Collaboration XXI. 2011, A&A, 536, A21

Popping, G., & Péroux, C. 2022, MNRAS, 513, 1531

Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152

Popping, G., Shivaei, I., Sanders, R. L., et al. 2023, A&A, 670, A138

Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25

Pozzi, F., Calura, F., Zamorani, G., et al. 2020, MNRAS, 491, 5073

Pozzi, F., Calura, F., Fudamoto, Y., et al. 2021, A&A, 653, A84

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85

Scoville, N. Z. 2013, Evolution of star formation and gas, eds. J. Falcón-Barroso, & J. H. Knapen, 491

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Shapley, A. E., Cullen, F., Dunlop, J. S., et al. 2020, ApJ, 903, L16

Shivaei, I., Popping, G., Rieke, G., et al. 2022, ApJ, 928, 68

Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 497, 956

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Sun, F., Egami, E., Fujimoto, S., et al. 2022, ApJ, 932, 77

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27

Tamura, Y. C., Bakx, T. J. L., Inoue, A. K., et al. 2023, ApJ, 952, 9

Triani, D. P., Sinha, M., Croton, D. J., Pacifici, C., & Dwek, E. 2020, MNRAS, 493, 2490

Valiante, R., Schneider, R., Bianchi, S., & Andersen, A. C. 2009, MNRAS, 397, 1661

Venemans, B. P., McMahon, R. G., Walter, F., et al. 2012, ApJ, 751, L25

Vijayan, A. P., Clay, S. J., Thomas, P. A., et al. 2019, MNRAS, 489, 4072

Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1

Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67

Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

Weaver, J. R., Davidzon, I., Toft, S., et al. 2023, A&A, 677, A184

Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165