© The Authors 2026

1. Introduction

The Euclid mission (Euclid Collaboration: Mellier et al. 2025) will observe 14 000 deg² of extragalactic sky, detecting billions of galaxies at optical wavelengths (with the VIS instrument at 550–900 nm; Euclid Collaboration: Cropper et al. 2025) and near-infrared wavelengths (with the Near-Infrared Spectrometer and Photometer, NISP at 1–2 μm; Euclid Collaboration: Jahnke et al. 2025). The first Quick Data Release (Q1; Euclid Quick Release Q1 2025) provided single-exposure observations covering three deep fields: Euclid Deep Field Fornax (EDF-F); Euclid Deep Field North (EDF-N); and Euclid Deep Field South (EDF-S). Even at the current depths (about magnitude 24.7 in VIS and 23.2 in NISP) the catalogues generated from these observations contain over 10 million galaxies detected in the Euclid filters (Euclid Collaboration: Romelli et al. 2026).

Spectral energy distributions (SEDs) have been fitted to the Euclid-selected galaxies, providing robust photometric redshifts, stellar masses, and star-formation rates (SFRs) for the majority of the galaxies (Euclid Collaboration: Tucci et al. 2026). This catalogue was recently used to constrain the correlation between the stellar mass (M_*) and SFR of star-forming galaxies (known as the galaxy star-forming main sequence, or MS) out to z  =  3 (Euclid Collaboration: Enia et al. 2026). The correlation is related to universal processes that have been converting cold gas reservoirs into stars since at least z  =  6. Because the bulk of the galaxies in the Universe follow the star-forming MS, determining its evolution is crucial for understanding galaxy evolution in general. For example, it is known that the amplitude of the MS increases with redshift for galaxies of all stellar mass (Speagle et al. 2014; Daddi et al. 2022; Popesso et al. 2023), implying that the specific SFRs (sSFRs) of all galaxies were higher in the early Universe. There is also a deviation from a linear trend at high stellar mass, meaning that there is a maximum average SFR at a given epoch, and this characteristic bending mass also increases with redshift. This has been attributed to a change in environments suppressing cold gas accretion with redshift and quenching (e.g. Dekel & Birnboim 2006; Daddi et al. 2022).

While optical and near-infrared (IR) light observed from the Earth is primarily sensitive to the stellar emission from galaxies at all redshifts (with the longer wavelengths being weighted to higher-redshift galaxies), far-IR light (from tens to around 1000 μm; note that this also includes wavelengths often described as submillimetre at the longer end) is sensitive to the thermal emission from warm dust grains in galaxies at all redshifts. Since these dust grains are primarily heated by hot young stars in star-forming regions, there is a tight correlation between the far-IR luminosity (often defined as the integral of the luminosity density between 8 and 1000 μm) and the SFR (e.g. Kennicutt 1998), suggesting a connection to the galaxy MS. It is therefore of interest to check the consistency between the SFRs derived from optical and near-IR photometry to those derived (independently) from far-IR photometry, and to see if there is evolution in the dust properties that can provide insight into the processes changing the MS.

There have been many large extragalactic surveys carried out by far-IR and submillimetre observatories, particularly by the Herschel Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) at 70–160 μm, the Herschel Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010) at 250, 350, and 500 μm, and the Submillimetre Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013) at 450 and 850 μm. Herschel effectively measures the peak of the thermal SEDs of most Euclid galaxies, providing effective constraints on dust temperatures (which are proportional to the peak frequency), while SCUBA-2 probes the Rayleigh–Jeans tail of the thermal SED, whose slope is related to the dust emissivity index β. Moreover, Herschel has surveyed roughly 1000 deg² of extragalactic sky (and SCUBA-2 about 10 deg²), essentially all of which will ultimately overlap with Euclid by the final data release.

Compared to the exquisite angular resolution of Euclid (approximately 0$\stackrel{″}{.}$2), the angular resolution of Herschel and SCUBA-2 is much coarser, ranging between about 10″ and 30″. Far-IR maps therefore do not individually detect most of the galaxies that Euclid sees, but rather blend these galaxies together into a coherent pattern that traces the large-scale structure. Making progress requires stacking, where the average value of the pixels at the locations of many objects is calculated, rather than trying to measure the flux densities of each object directly from the map. More precisely, this operation calculates the covariance between a catalogue and a map (see e.g. Marsden et al. 2009; Wang et al. 2015). The Herschel properties of optical and near-IR catalogues have been investigated via stacking in other fields such as the UKIDSS Ultra-Deep Survey (UDS; Viero et al. 2013, 1 deg²), the Cosmic Evolution Survey (COSMOS; Simpson et al. 2019; Duivenvoorden et al. 2020, 2 deg²), both UDS and COSMOS (Koprowski et al. 2024), and the Galaxy And Mass Assembly (GAMA) fields and Stripe 82 region (Wang et al. 2016). However, the fields used in these studies are either small and suffer from sample variance to some degree, or the optical catalogues do not go beyond redshift 1. With Euclid we can dramatically expand these results to include millions of star-forming galaxies across tens of square degrees. Moreover, future Euclid data releases will continue to overlap with existing Herschel observations, eventually amounting to a billion galaxies over 1000 deg².

We therefore focus on stacking the entire Q1 catalogue on Herschel and SCUBA-2 maps, making this the largest study yet of this kind. In Sect. 2 we describe the Euclid catalogues and far-IR maps used in the analysis. In Sect. 3 we present our stacking pipeline. In Sect. 4 we show our results, and in Sect. 5 we discuss our findings. The paper concludes in Sect. 6. Throughout this paper we assume the cosmological parameters from Planck Collaboration VI (2020).

2. Data

The Euclid Q1 release is split into three fields: the EDF-F (12.1 deg²); the EDF-N (22.9 deg²); and the EDF-S (28.1 deg²). Of these three fields, EDF-F and EDF-N have overlapping coverage from SPIRE. The SPIRE field overlapping with the EDF-F is known as ‘CDFS-SWIRE’, and the field overlapping with the EDF-N is known as the ‘AKARI-NEP’. Here we describe the multiwavelength data in these two fields.

2.1. Euclid catalogues and masks

The Euclid merging (MER) Q1 catalogue (Euclid Collaboration: Romelli et al. 2026) contains the photometry of all VIS- and NISP-detected galaxies, as well as ground-based photometry in the u, g, r, i, and z bands from various telescopes (see Tereno et al., in prep., for details). Here we make use of the Euclid catalogue described in Euclid Collaboration: Enia et al. (2026), which includes additional photometry from the Infrared Array Camera (IRAC) on board the Spitzer Space Telescope (Fazio et al. 2004) at 3.6 and 4.5 μm. This catalogue also includes refitting of SEDs with the additional IRAC data in order to derive photometric redshifts, stellar masses, and SFRs. The final catalogue contains 2 884 906 objects in the EDF-F and 6 221 146 objects in the EDF-N. We note that this contains somewhat fewer objects than the full Q1 catalogue, since stars and other image artefacts were removed before cross-matching to IRAC.

Since in this study we are interested in the average properties of MS galaxies, we used the same colour cuts to remove quiescent galaxies. Specifically, galaxies were removed with NUV  −  r⁺ > 3(r⁺  −  J)  +  1 and NUV  −  r⁺ > 3.1, leaving 1 318 898 objects in the EDF-F and 2 658 118 objects in the EDF-N.

To ensure an accurate stacking analysis, we must account for the fact that some regions within the Euclid footprint do not contain any extragalactic objects due to contamination from bright stars. This can be done by calculating a mask that is associated with the Euclid catalogue, and propagating that mask through all the calculations.

For each of our far-IR images, we therefore calculated the number of Euclid objects in each pixel, then lightly smoothed the map using a Gaussian kernel. We then defined the Euclid catalogue mask to be the regions where this smoothed map has a value less than a given threshold, determined visually by ensuring that the mask agreed well with the actual galaxy distribution. Defined this way, the Euclid mask removes far-IR pixels that have no extragalactic Euclid objects across a sufficiently large scale.

2.2. Far-IR imaging

The SPIRE maps of the EDF-F (or CDFS-SWIRE) and the EDF-N (or AKARI-NEP) were obtained from the Herschel Extragalactic Legacy Project (HELP) archive (Shirley et al. 2021)¹. The HELP data products include unfiltered maps and matched-filtered maps, along with their corresponding noise (or RMS) maps. Throughout this paper we use the unfiltered products because we do not expect these maps to have significant fluctuations caused by instrumental effects, nor significant contamination from dust in the Milky Way. The CDFS-SWIRE image covers 12.8 deg², while the AKARI-NEP image covers 9.0 deg² (although not all of this area overlaps with Euclid). The raw SPIRE data, RMS maps and masks for each of the two fields are shown in Figs. 1 and 2. In both maps extra observations were taken near the centres of the fields to decrease the noise and create small deep fields and wide shallow fields, and the HELP data products combined all of the observations into a single image.

Fig. 1. Top: Herschel-SPIRE data covering the CDFS-SWIRE field (overlapping with the EDF-F) at 250, 350, and 500 μm. The blue contour shows the mask applied to the SPIRE images to remove bad edge pixels. The grey contours show the corresponding Euclid catalogue mask, where masked rectangles designate the locations of bright stars in the field contaminating source extraction. Bottom: Same as the top panel, but showing the RMS of the Herschel-SPIRE data. Coordinates are conventional RA and Dec.

Fig. 2. Same as Fig. 1, but for the Herschel-SPIRE AKARI-NEP field (overlapping with the EDF-N).

In addition to SPIRE maps, the HELP archive also provides PACS maps at 100 and 160 μm, wherever the data were taken. For the EDF-F, the full area observed by SPIRE was also observed by PACS (albeit to a much shallower depth), which we make use of here. The AKARI-NEP was also observed, covering 0.6 deg² of the EDF-N. For these data only the output from the standard PACS processing pipeline is available; some filtering was done to the raw data in order to reduce strong 1/f noise and remove artefacts from bright sources; however, the maps were not matched-filtered. The PACS data, RMS maps and masks are shown in Appendix A.

SCUBA-2 was used to map the AKARI-NEP (EDF-N) field as part of the SCUBA-2 Cosmology Legacy Survey (S2CLS; Geach et al. 2017), and this field was later expanded in the North Ecliptic Pole SCUBA-2 survey (NEPSC2; Shim et al. 2020)². The data products include the unfiltered maps and the matched-filtered maps, along with their corresponding RMS maps. As with the SPIRE images, we use the unfiltered maps here. The total area covered by SCUBA-2 in the AKARI-NEP is 2.9 deg², and the entirety of this field has been observed by Euclid. The SCUBA-2 data, RMS map and mask are shown alongside the PACS maps in Appendix A.

Despite the fact that we can down-weight noisy regions using RMS maps, very noisy pixels near the edges of these maps can still cause significant issues, since the uncertainties for regions with few a ‘hits’ can be underestimated. To mitigate this issue, we also created masks to remove problematic edge pixels. To make the mask for the SPIRE CDFS-SWIRE field we smoothed the noise maps using a Gaussian kernel with a standard deviation of 7.5 pixels, then masked regions where the smoothed noise map has a value > 7 mJy beam⁻¹; for reference, typical noise values are < 4 mJy beam⁻¹. For the AKARI-NEP field this strategy does not work because the noise map is too inhomogeneous, so instead we manually defined a rectangular masked region outside of which the noise pattern begins to deviate from the majority of the map. For the PACS maps of the CDFS-SWIRE field the noise is again quite inhomogeneous, so we also manually defined a central region with a consistent noise level. For the PACS AKARI-NEP maps we smoothed the noise maps using a Gaussian kernel with a standard deviation of 5 pixels and masked regions where the smoothed noise map has a value > 40 mJy beam⁻¹ (compared to the typical noise level of about 30 mJy beam⁻¹). For the SCUBA-2 data, we simply masked pixels with noise values > 30 mJy beam⁻¹ (compared to the typical noise level of about 10 mJy beam⁻¹). Lastly, the AKARI-NEP field contains NGC 6543 (the ‘Cat’s Eye Nebula’), which is particularly bright in the PACS and SCUBA-2 images. Since this is a Galactic object, we masked it before stacking our galaxy catalogue. These masks are also shown in Figs. 1, 2, A.1, and A.2 along with the Euclid catalogue masks.

2.3. Overlap between Euclid Q1 products and far-IR images

After accounting for the masked regions, the total overlapping area between the EDF-F and the SPIRE CDFS-SWIRE data is 10.8 deg², while for the EDF-N and the SPIRE AKARI-NEP data the total overlapping area is 6.8 deg², for a total area of 17.6 deg². Within the unmasked EDF-F region, the catalogue from Euclid Collaboration: Enia et al. (2026) contains about 1.5 million MS galaxies, while the EDF-N regions contains about 1.1 million galaxies (note that the exact values beyond the first digit depend somewhat on the SPIRE wavelength, since each SPIRE band has slightly different coverage). This brings the total number of MS galaxies overlapping with SPIRE to 2.6 million. For PACS, the total area available for stacking, after accounting for the Euclid mask, is 8.9 deg² in the EDF-F and 0.4 deg² in the EDF-N, with 1.3 million and 70 000 MS galaxies, respectively. The total PACS area is thus 9.3 deg², containing 1.4 million MS galaxies. For SCUBA-2, the total area available for stacking (only in the EDF-N) is 2.4 deg², and contains 240 000 galaxies.

3. Stacking method

3.1. Stacking algorithm

With a few small adjustments, we employed the stacking method called SimStack (Viero et al. 2013), which stacks on multiple catalogue bins simultaneously. The main advantage of SimStack is that it is not affected by galaxy clustering, since we do not need to assume that the galaxies being stacked are Poisson distributed. Briefly, after the input stacking catalogue was split into N bins, we solved for the N stacking amplitudes defined by the linear equation

$$\begin{array}{c}\hfill y=S_1{X}_1+\cdots +S_N{X}_N,\end{array}$$(1)

where y is the data map being stacked on (in this case a Herschel or SCUBA-2 map at one wavelength), S_b are the average flux densities of the galaxies in bin b, and X_b are the beam-convolved images of the distribution of galaxies in the same bin. We note that in this analysis we subtracted the noise-weighted mean from the data map y as well as the noise-weighted means from each of the beam-convolved images X_b. This ensures that we did not need to add an arbitrary constant to our stacking model (see Viero et al. 2013). To construct the beam-convolved galaxy distribution images X_b we simply looped through the positions of all the galaxies in a given bin, adding a 1 to each pixel in the model image where a galaxy in bin b is located. Then, we convolved the image by the instrumental beam and subtracted the weighted mean.

Equation (1) is a linear system and can thus be solved via the weighted least-squares method (see Section 3.1 of Viero et al. 2013 for an explicit derivation). Moreover, we can use maps from multiple fields to fit for the stacked flux densities simultaneously. For M total pixels across all of the data maps, we defined X as the M  ×  N matrix where column i contains the beam-convolved images X_i and y is the M  ×  1 vector of the data (note that X_i and y must be flattened from a 2D image to a 1D vector in the same way, and that images from multiple fields must be added to the vectors and matrices in the same order). We also incorporated the weights (and masks) in the diagonal M  ×  M matrix W, where W_ii  =  1/σ_ii² if the pixel is not masked, and nan otherwise. The best-fit stacking amplitudes are then just the solution to the weighted linear least-squares problem,

$$\begin{array}{c}\hfill \widehat{}\}S={({\mathsf{X}}^{\mathsf{T}}\mathsf{W}\mathsf{X})}^{-1}{\mathsf{X}}^{\mathsf{T}}\mathsf{W}\mathbf{y}.\end{array}$$(2)

In traditional stacking, creating small equally sized cutouts around each source and averaging these cutouts together is common. Ideally, for far-IR images of unresolved galaxies, these should all resemble the instrumental point-spread function (PSF), and so these 2D stacks are excellent tools for checking for systematic errors. We can easily generalise the 2D stacking method from regular stacking to SimStack by noting that while the central pixel in a regular 2D stack is the covariance between the map and the catalogue, the offset pixels represent the shifted cross-correlations. We can achieve the same result by calculating cross-correlations between our data y and our model S₁X₁ + … + S_NX_N. In practice, we offset the data images and corresponding weights by integer pixel values (Δx, Δy) relative to the model image, then re-solved Eq. (2) for a new set of stacking flux densities $\widehat{\mathbf{S}}$.

3.2. Stacking bins

To perform the stacking we defined a set of N bins in stellar mass and redshift, produced from the Euclid catalogue. Euclid Collaboration: Enia et al. (2026) focused their study of the galaxy MS on galaxies with 0.2 < z < 3.0, with the understanding that most catastrophic outliers have either mistakenly very low or very high redshifts. For a similar reason they also limited their sample to galaxies with log₁₀(M_*/M_⊙)  ≲  11.5, since the galaxies with larger stellar masses are likely to be outliers with very large redshifts. Finally, they estimated the stellar mass completeness of the sample as a function of redshift, finding > 95% completeness at the lowest redshift (z  =  0.2) around log₁₀(M_*/M_⊙)  ≃  8. We therefore split the Euclid catalogue into redshift bins between z  =  0.2 and 3.0, with a spacing of Δz  =  0.2, and into stellar mass bins between log₁₀(M_*/M_⊙)  =  8.3 and 11.5, with a spacing of Δlog₁₀(M_*/M_⊙)  =  0.4 .

We next added an additional layer to the stacking model that includes all remaining catalogued galaxies that do not fall into any redshift or stellar mass bin. We also included the quiescent galaxies in this additional bin. We finally added one more layer to the stacking model to account for the Euclid mask, required for reducing biases in the stack (see Duivenvoorden et al. 2020, for details). For this mask layer we inverted the Euclid mask, convolved it with the instrumental beam, and then restored the masked pixels. This layer is intended to capture far-IR surface brightness leakage from near-IR galaxies that Euclid was unable to detect due to contamination from, for example, bright stars (although in practice this effect is small).

3.3. Instrumental beams

The final input required for stacking is a model of the Herschel and SCUBA-2 beams. For PACS, we used the model empirical beams provided by the HELP archive (Shirley et al. 2021). These were measured by stacking WISE-selected galaxies on the PACS maps and fitting an elliptical 2D Gaussian profile to the stack. The ellipticity is due to the fact that the PACS beam is not circularly symmetric, and so for these data we used the best-fit beams for each PACS map separately.

For SPIRE, since the input maps are not filtered, we used the instrumental beams approximated as 2D Gaussians with full width half maximum (FWHM) values of 18$\stackrel{″}{.}$15 at 250 μm, 25$\stackrel{″}{.}$15 at 350 μm, and 36$\stackrel{″}{.}$3 at 500 μm (see Griffin et al. 2010). For the SCUBA-2 image, we used the updated beam profile from Mairs et al. (2021), which is the sum of two Gaussians. The first Gaussian (the main beam) has $\mathrm{FWHM}=11\stackrel{″}{.}0$ and a relative amplitude of 0.98, while the second Gaussian (the error beam) has $\mathrm{FWHM}=49\stackrel{″}{.}1$ and a relative amplitude of 0.02.

3.4. Estimating uncertainties

We estimated the uncertainties in the best-fit stack flux densities following Viero et al. (2013) by both propagating the weight matrix and by performing bootstrap resampling to determine the overlap between neighbouring bins. The statistical covariance matrix from solving the weighted linear least-squares system is analytic and can be calculated as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\widehat{\mathsf{S}}}={({\mathsf{X}}^{\mathsf{T}}\mathsf{W}\mathsf{X})}^{-1},\end{array}$$(3)

and so the statistical uncertainty in ${\widehat{S}}_i$ is just $\sqrt{{\mathrm{\Sigma }}_{\widehat{S},ii}}$. We added this in quadrature with the uncertainty from bootstrap resampling (which dominates the error budget); for this contribution, we generated 100 random catalogues by drawing stellar masses and redshifts for each galaxy from their probability distributions. In principle the stellar masses and redshifts are correlated with each other and would require a full 2D posterior distribution for each galaxy; however, we simplified the procedure by assuming that the distributions are Gaussian with a standard deviation equal to half the 68% confidence interval. We then re-binned each random catalogue, re-solved Eq. (2), and calculated the standard deviation of the 100 random estimates of $\widehat{\mathbf{S}}$.

4. Results

4.1. Average far-IR flux densities of main-sequence Euclid galaxies

We ran our stacking algorithm (essentially SimStack generalised to provide 2D cross-correlations) on the Herschel and SCUBA-2 maps described in Sect. 2. We set the 2D stacking cutout size for each of the SPIRE wavelengths to be 200″ (meaning that we solved Eq. (2) for offsets in a 200″  ×  200″ grid), while for SCUBA-2 we set the cutout size to be 70″, and for PACS we set the cutout size to be 40″. The full 2D cross-correlations are shown in Appendix B. We find significant detections in most bins above log₁₀(M_*/M_⊙)  =  10 , and the detections are generally consistent with the instrumental PSFs.

More precisely, for perfect point sources the 2D profile produced by our algorithm is the cross-correlation of the beam with itself; for a Gaussian, this increases the FWHM by a factor of $\sqrt{2}$. We tested this by computing the averaged 1D radial profiles of our 2D cross-correlations and comparing these to the expected PSF profiles. We find good agreement between the stacked signals and the PSFs for most bins with log₁₀(M_*/M_⊙) > 9.9. Below this stellar mass we find that the 2D profiles can be more extended than the beam, with the largest effect seen at 500 μm where the PSF is largest. This could be caused by a combination of the beam size and incompleteness in the Euclid catalogue. For example, dust galaxies bright at far-IR wavelengths are likely to be undetected by Euclid, and these galaxies could have different clustering properties. However, as detailed in the following sections, we do not expect this extended emission to affect our results and so we do not attempt to correct for it.

As discussed previously, the central pixel values shown in the 2D stacks are the best-fit mean flux densities of the galaxies in the given bin (i.e., the covariances between the catalogue and the maps). For the high S/N detections, we checked that the central pixel in our 2D stacks was consistently the brightest pixel, meaning that there were no astrometric offsets between the Euclid catalogue and our stacking images. We provide these central values and the corresponding uncertainties in Appendix B.

4.2. Average SEDs of main-sequence Euclid galaxies

Our best-fit stacked flux densities contain information about the average far-IR SEDs of the Euclid-selected galaxies in each stellar mass and redshift bin. We modelled the lower frequencies of these SEDs using a standard modified blackbody for ν(1 + z) < ν_α (explained below):

$$\begin{array}{c}\hfill S_{\nu }[\mathrm{mJy}]=A{\left[\frac{\nu (1+z)}{{\nu }_0}\right]}^{\beta }(1+z)B_{\nu (1+z)}(T_{\mathrm{d}})\left[\frac{{10}^{29}\mathrm{mJy}}{\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{Hz}}^{-1}}\right],\end{array}$$(4)

where B_ν(T_d) is the blackbody function

$$\begin{array}{c}\hfill B_{\nu }(T_{\mathrm{d}})=\frac{2h{\nu }^3}{c^2}\frac{1}{{\mathrm{e}}^{h\nu /k{T}_{\mathrm{d}}}-1}·\end{array}$$(5)

Since we dealt with galaxies at z < 3, we ignored corrections related to the cosmic microwave background (CMB; see da Cunha et al. 2013). In these equations ν is the observed frequency, z is the redshift, h is the Planck constant, k is the Boltzmann constant, c is the speed of light, T_d is the dust temperature, and the factor of 10²⁹ converts the observed flux density units from SI to mJy. The quantity ν₀ is not a free parameter of the model but sets the dust mass normalisation (as explained in Sect. 4.3), and we used ν₀  =  353 GHz. We note that this simple parametrisation assumes the dust is optically thin, so we did not need to include the characteristic frequency where the optical depth is unity (e.g. Draine 2006; Drew & Casey 2022). This left the amplitude A, dust temperature T, and dust emissivity index β as the three free parameters.

At rest-frame mid-IR frequencies ν(1 + z) > ν_α, i.e. around a few terahertz, the Wien side of the thermal SED falls off less steeply than an exponential and is typically modelled as a power law with a slope of −α (e.g. Blain et al. 2003; Roseboom et al. 2013; Casey et al. 2014; Reuter et al. 2020). Since our PACS data cover this region of the SED, we needed to include this phenomenological feature. To do so, we modelled the SED for ν(1 + z) > ν_α as

$$\begin{array}{c}\hfill S_{\nu }[\mathrm{mJy}]=A_{\alpha }{\left[\frac{\nu (1+z)}{{\nu }_{\alpha }}\right]}^{-\alpha }(1+z)\left[\frac{{10}^{29}\mathrm{mJy}}{\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{Hz}}^{-1}}\right].\end{array}$$(6)

Here, A_α and ν_α are determined by matching the amplitude and slope of Eqs. (4) and (6). This gives A_α  =  A(ν_α/ν₀)^βB_να(T), and ν_α solves the equation 3  +  β  +  α  = xe^x/(e^x − 1) where x  =  hν_α/kT.

Since we did not have enough photometry information to constrain β on the Rayleigh–Jeans sides or α on the Wien sides of the SEDs, we fixed these values to β  =  1.96 and α  =  2.3 according to the mean values found for local far-IR-selected galaxies (Drew & Casey 2022). Fixing β and α to these values, we performed our SED fits on all redshift and stellar mass bins where all three SPIRE flux densities were detected in the stack with S/N> 3, and where there is also a PACS detection at 100 or 160 μm with S/N> 3. The reason for these constraints is simply that we needed the observed photometry to bracket the peak of the SED to properly constrain the dust temperature (which is directly proportional to the peak frequency). We included a 4% absolute calibration uncertainty and a 1.5% relative calibration uncertainty in the SPIRE bands, a 5% absolute calibration uncertainty in the PACS bands, and a 15% absolute calibration uncertainty in the SCUBA-2 band. The resulting best-fit SEDs are shown in Fig. 3, and the best-fit parameters are given in Appendix C. In this figure we do not show the lowest stellar mass bins as we do not have enough far-IR photometry to fit SEDs. For two of the highest redshift and stellar mass bins, z  =  2.7, log₁₀(M_*/M_⊙)  =  10.9 and 11.3, the fits provide unphysically high dust temperatures (> 60 K), potentially due to contamination in the Euclid catalogue. We discard these two bins for the remainder of this work. As shown, it is mainly the log₁₀(M_*/M_⊙) > 9.9 bins where we have enough far-IR photometry to fit SEDs, and these bins show good agreement between the 2D cross-correlation profiles and the instrumental beams.

Fig. 3. Modified blackbody SEDs (with β = 1.96 and α = 2.3) fitted to the stacked Herschel and SCUBA-2 flux densities. The redshift and stellar mass of each bin are indicated by the top and right axis labels, respectively. The best-fit parameters are given in Table C.1. The SEDs were fitted to bins where all three SPIRE flux densities are detected with S/N> 3 and at least one PACS flux density is detected with S/N> 3. The panels are otherwise blank. Bins that are > 95% complete in stellar mass are highlighted in blue.

4.3. Derived parameters for main-sequence Euclid galaxies

From the best-fit parameters to Eq. (4) we can derive certain average physical properties for the galaxies in the Euclid catalogue. First, the best-fit temperatures T_d are already the rest-frame dust temperatures because we included the redshifts in the fits. Next, the best-fit amplitudes are directly proportional to the dust mass, M_d, and can be calculated using (e.g. Reuter et al. 2020; Eales & Ward 2024; Jolly et al. 2025)

$$\begin{array}{c}\hfill M_{\mathrm{d}}=\frac{D_{\mathrm{L}}^2(z)A}{{\kappa }_0},\end{array}$$(7)

where κ₀ is the calibration factor that scales the specific luminosity at the rest-frame frequency of ν₀ (the same reference frequency in our SED fit, see Eq. (4)) to a dust mass and D_L is the luminosity distance. Here, we used the factor κ₀  =  0.077 m² kg⁻¹, calibrated to the frequency ν₀  =  353 GHz (Dunne et al. 2000; da Cunha et al. 2008; Dunne et al. 2011), and we included an uncertainty of ±0.02 m² kg⁻¹ (from James et al. 2002). This approach uses the best-fit model to estimate the rest-frame specific luminosity at the frequency ν₀, as opposed to using a single measured flux density. There are many other values of κ₀ used throughout the literature, typically resulting in dust mass differences of a factor of a few; however, picking a different κ₀ only changes the absolute value of the dust mass, rather than trends in stellar mass or redshift (although κ₀ could in principle vary with redshift). On the other hand, the dust emissivity index β can have an effect on the best-fit SEDs, and therefore any stellar mass and redshift trends. While we do not expect β to vary by much compared to the fiducial value of about 2, this could be investigated in the future if more far-IR wavelengths could be added.

Lastly, we can estimate far-IR SFRs using the linear relation

$$\begin{array}{c}\hfill \mathrm{SFR}[{\mathrm{M}}_{\odot }y{r}^{-1}]=1.49\times {10}^{-10}{L}_{\mathrm{FIR}}[{\mathrm{L}}_{\odot }]\end{array}$$(8)

(assuming a Kroupa initial mass function as performed with the Euclid catalogue; see Euclid Collaboration: Tucci et al. 2026), where

$$\begin{array}{c}\hfill L_{\mathrm{FIR}}=4\pi D_{\mathrm{L}}^2(z)[{\displaystyle {\int }_{{\nu }_1}^{{\nu }_{\alpha }}}A{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }{B}_{\nu }(T_{\mathrm{d}})\mathrm{d}\nu +{\displaystyle {\int }_{{\nu }_{\alpha }}^{{\nu }_2}}{A}_{\alpha }{\left(\frac{\nu }{{\nu }_{\alpha }}\right)}^{-\alpha }\mathrm{d}\nu ]\end{array}$$(9)

is the area under the far-IR SED between ν₁  =  c/1000 μm and ν₂  =  c/8 μm in the rest frame. Again, we used β  =  1.96 and α  =  2.3. We note that this derived physical parameter SFR is not independent of T_d and M_d, but instead these three quantities are correlated.

In Appendix C we show the resulting physical parameters T_d, M_d, and SFR in each of the stellar mass and redshift bins where we have sufficient stacked photometry to derive these physical parameters, and we also provide the derived physical quantities. We find a clear trend of increasing dust temperature with redshift, and an increase in the dust mass from redshifts to around z  ≃  2. For the SFRs, we see an increase towards high-z, with a dependence on the stellar mass (as expected from the galaxy MS). These trends are discussed further in Sect. 5.

4.4. The mean brightness of the CIB from stacking

We now turn to estimating the fraction of the cosmic infrared background (CIB) resolved by Euclid-selected galaxies in Q1. The average brightness of the extragalactic sky at far-IR wavelengths was measured by the Cosmic Background Explorer (COBE) satellite, which carried two relevant instruments: the Far Infrared Absolute Spectrophotometer (FIRAS; Mather et al. 1993); and the Diffuse Infrared Background Experiment (DIRBE; see Boggess et al. 1992). Analysis of these data provide the best available estimates of the absolute value of the surface brightness of the sky at far-IR wavelengths. Herschel and SCUBA-2, on the other hand, have no sensitivity to the monopole on the sky, but instead measure differences between sources and an unknown background, amounting to fluctuations caused by individual galaxies. The fraction of the CIB can be estimated by taking the sum of detected sources and dividing by the COBE-estimated background.

To do this, for each bin in the stack we multiplied the mean flux density by the total number of SPIRE galaxies in the bin, then summed the contribution from all of the bins (including surface brightness leakage from the mask and the galaxies that do not fall in any redshift or stellar mass bin or are classified as quiescent). Finally, we divided this total flux density by the area of the SPIRE map after applying the masks (these are the areas given in Sect. 2.3). The results are shown in Table 1. We note that > 70% of our measured CIB surface brightness from star-forming galaxies comes from log₁₀(M_*/M_⊙) > 9.9 bins which show good agreement between the 2D cross-correlation signals and the instrumental beams.

To determine the fraction of the CIB resolved by Euclid, we estimated the absolute value of the CIB measured by FIRAS and DIRBE. FIRAS measurements were recalculated using an improved Galactic emission removal procedure at Planck wavelengths by Odegard et al. (2019), which should be more accurate than FIRAS alone at 350, 500, and 850 μm. From DIRBE, we used the results from a re-analysis of the data from the Cosmoglobe project (Watts et al. 2024), which provides measurements of the CIB intensity at 100, 140, and 240 μm. Additionally, Casandjian et al. (2024) used a combination of DIRBE and Planck data to estimate the CIB at 100, 140, 240, 350, 500, and 850 μm. Lastly, we used the result from FIRAS alone at 250, 350, 500, and 850 μm by taking the best-fit spectral shape from (Fixsen et al. 1998). Here, the fit was performed using a modified blackbody function, with β  =  0.64 and T  =  18.5 K. For the uncertainty we used the uncertainty in the best-fit amplitude. To correct the above 140 μm intensities to 160 μm, the 240 μm intensities to 250 μm, and the 550 μm intensities to 500 μm, we used this same best-fit spectral shape – these corrections range from 1 to 24%. All the CIB estimates are given in Table 1.

5. Discussion

5.1. Redshift trends of mean physical properties

In our stacking analysis, we found significant redshift trends for the dust temperatures, dust masses, and SFRs of average Euclid-selected galaxies. To interpret these results, we used the 95% stellar mass completeness limits from Euclid Collaboration: Enia et al. (2026) to select stacking bins that are > 95% complete. These bins are highlighted in blue in Fig. 3.

Focusing only on the stellar mass and redshift bins where the stellar masses are > 95% complete, we first compared our far-IR-derived SFRs to that predicted from the star-forming MS. The Euclid SFRs were estimated by modelling solely optical and near-IR photometry from Euclid and its supporting observations, and the star-forming galaxies were split into four redshift bins between z  =  0.2 and z  =  3.0. In each bin the MS was parametrised as

$$\begin{array}{c}\hfill \mathrm{SFR}=\frac{{\mathrm{SFR}}_{\max }}{1+{(M_0/M_{\ast })}^{\gamma }},\end{array}$$(10)

with γ fixed to 1 and SFR_max and M₀ left as free parameters. Since our redshift bins are much smaller than those used by Euclid Collaboration: Enia et al. (2026), we instead used the fit to Eq. (10) from Popesso et al. (2023), which was found to be in good agreement with the Euclid results. Popesso et al. (2023) also fixed γ  =  1, letting log₁₀(SFR_max)(t)  =  a₀  +  a₁t and log₁₀(M₀)(t)  =  a₂  +  a₃t. The resulting ratio of far-IR SFRs to the prediction from the MS are shown in Fig. 4. We note that below log₁₀(M_*/M_⊙)  =  9.7 we do not have any far-IR detections in the stellar mass-complete bins, so we only plot the results for the five most massive bins. We find good agreement between the two estimates, with a weighted mean ratio of 1.0 and a weighted standard deviation of 0.2.

Fig. 4. Ratio of our SFRs measured from far-IR photometry to a parametrisation of the star-forming MS, shown as a function of redshift and split into different stellar mass stacking bins. Only redshift and stellar mass bins with > 95% completeness are shown. We use the MS parametrisation from Popesso et al. (2023), which is a continuous function of z/t and is in good agreement with Euclid Collaboration: Enia et al. (2026).

However, we do see a systematic trend where more massive galaxies have higher SFRs than expected from the MS (and vice versa for less massive galaxies), as well as a trend where we estimate overall higher SFRs than the MS for z  ≳  1.5. It is worth noting that the SFRs used in Popesso et al. (2023) are a compilation from the literature and include estimates using ultraviolet luminosities, Hα, and far-IR luminosities, and the authors do not note any systematic differences between the estimators. We also compared our SFRs to the best-fit MS from Euclid Collaboration: Enia et al. (2026) after interpolating between the redshift bins, and directly to the SFRs from the Euclid catalogue after averaging over the same redshift and stellar mass bins, finding similar results. It is thus more likely that the systematic differences we observe here are due to differences in the assumed far-IR SED shapes (for example the dust emissivity index β, the slope α, whether an evolving dust temperature with redshift was assumed, or if the dust temperature depended on the stellar mass), or biases in the stellar masses and photometric redshifts in the Euclid catalogue. However, a more complete understanding of these issues is beyond the scope of this paper.

Next, in Fig. 5 we show T_d versus time (bottom axis) and z (top axis), split by the different stellar mass bins used in our stacking analysis. We find that average dust temperatures increase with redshift from about 20 to 35 K, without any obvious trend in stellar mass. Looking at the dust temperature trend plotted linearly as a function of time, we see that the data decrease steeply between 2 and 6 Gyr (z  =  3 and 1), then plateau to a constant dust temperature until the present day. To capture this behaviour, we fitted a phenomenological function of the form

$$\begin{array}{c}\hfill T_{\mathrm{d}}(t)=T_2+(T_1-T_2){\mathrm{e}}^{-t/\tau },\end{array}$$(11)

Fig. 5. Best-fit dust temperatures from our SED fitting, Td, as a function of time (bottom axis) and redshift (top axis), considering only the redshift and stellar mass bins that are > 95% complete. We show the dust temperature evolution for five different stellar mass bins, with the stellar mass values of the centres of the bins given in the legend. The solid curve is a fit to the simple form T2  +  (T1  −  T2) e−t/τ; the dotted line is the quadratic-in-redshift fit from Koprowski et al. (2024) and the dashed line is the linear-in-redshift fit from Schreiber et al. (2017). We also show published mean temperature estimates for star-forming galaxies at low redshifts (0.01 < z < 0.05; Lamperti et al. 2019) in blue.

where T₁ is the mean dust temperature of all star-forming galaxies at t  =  0 Gyr, T₂ is the dust temperature all star-forming galaxies approach as t becomes large, and τ is the characteristic timescale. We find best-fit values of T₁  =  (79.7 ± 7.4) K, T₂  =  (23.2 ± 0.1) K, and τ  =  (1.6 ± 0.1) Gyr.

For comparison, in Fig. 5 we also show a fit to the average dust temperature of star-forming galaxies as a function of redshift from Schreiber et al. (2017), which was later used to generate a simulated Euclid catalogue with far-IR photometry (see Sect. 5.3). Here, a similar stacking analysis was used to estimate far-IR photometry in bins of redshift and stellar mass. However, as opposed to fitting a modified blackbody to the photometry (as done here), an empirical template was used to fit the data. They found that the dust temperatures increased linearly as a function of redshift, but as seen in Fig. 5 this approach appears to overestimate the dust temperatures relative to our analysis around z  ≈  1.5. We also include the results from Koprowski et al. (2024), who stacked optically selected catalogues on far-IR images and fitted a modified blackbody to the photometry (with β and α fixed to the same values as used here). They found that the dust temperatures follow a quadratic polynomial in redshift, which agrees well with our exponential decay function. We note that neither of these studies found any trend in stellar mass. In Fig. 5 we also show the average dust temperature derived from a sample of local (0.01 < z < 0.05) star-forming galaxies from Lamperti et al. (2019). Since this sample includes bright galaxies, no stacking was required to estimate far-IR flux densities. The SEDs were fitted to the same modified blackbody function used here, although wavelengths > 100 μm were not used in the fit. Therefore, the transition to a power-law was not needed, and β was kept as a free parameter with best-fit values ranging from 1–2. After averaging over all of the dust temperatures in the sample, we find a mean value of 23.1 ± 0.1 K, consistent with our fit.

We next plotted the dust mass, M_d, and dust-to-stellar mass ratio, M_d/M_*, in the same way for the same stellar-mass-complete bins. The results are shown in Fig. 6. Here, we find an increase in dust mass (and the dust mass ratio) from z  =  0.2 to around z  =  1, followed by a plateau, for all stellar mass bins. We also find that the dust-to-stellar mass fraction decreases with increasing stellar mass at all redshifts.

Fig. 6. Left: Best-fit dust masses from our SED fits, Md, as a function of redshift, considering only the redshift and stellar mass bins that are > 95% complete. Right: Same as the left panel but showing the best-fit dust mass-to-stellar mass ratio, where the stellar masses are from the Euclid catalogue and are used to define the stacking bins. We also show SCUBA-2 stacking results from Millard et al. (2020), scaled to dust mass assuming the same modified blackbody SEDs used here and our best-fit dust temperature as a function of time, and ALMA 1.2 mm stacking results from Jolly et al. (2025) scaled in the same way. In both panels, the solid-coloured curves show the predicted trends for the corresponding stellar masses by combining the star-forming MS from Popesso et al. (2023) with our fit to the dust temperature as a function of time.

These trends can be explained using our fit to Eq. (11) and the fit to Eq. (10) from Popesso et al. (2023), or indeed any fit to the MS that agrees with our far-IR-measured SFRs. By combining these two equations we can solve for the dust mass as a function of t (and therefore z); see Appendix D for details. The resulting curves for the dust mass and dust-to-stellar mass ratio as a function of redshift are shown in Fig. 6 for each stellar mass bin. It is important to emphasise that these are not fits to the data, but predictions based on fits to the dust temperatures, SFRs, and stellar masses. We find an increase in dust mass and dust-to-stellar mass ratio up to z  =  1. Beyond this redshift, the curves predict a decrease in these quantities as opposed to a plateau. However, we note that the actual functional forms are very sensitive to the best-fit parameters to the dust temperature and SFR as a function of t, and we do not have enough stacked photometry beyond z  =  1.5 to properly constrain the model. We also see the same trend of decreasing dust-to-stellar mass fraction with increasing stellar mass. The decrease in the dust-to-stellar mass ratio at redshifts < 1 has been observed and discussed by Béthermin et al. (2015) as a consequence of decreasing gas mass (although increasing metallicity works against this). At high redshift galaxies have high gas-to-stellar mass ratios as material is accreted from the large-scale environment, while at low redshift this gas has been consumed by star formation or ejected and the accretion rate is much lower. Since the dust mass is roughly proportional to the gas mass, we expect the dust mass to be lower in low-redshift galaxies as well. While metallicity competes against this trend (and it is important to keep in mind that high-redshift galaxies at a given stellar mass are not the same objects as low-redshift galaxies with the same stellar mass), the lower gas mass dominates the overall evolution and we observe decreasing dust-to-stellar mass ratios for z < 1. Similarly, more massive galaxies have lower sSFRs (SFR/M_*), so they produce less dust per unit stellar mass, explaining the decrease in the dust-to-stellar mass fraction with stellar mass.

For comparison, we show the results from Millard et al. (2020) and Jolly et al. (2025). Both studies stacked on all optically selected galaxies (i.e. star-forming and quiescent) rather than just star-forming galaxies as done here, and both only stacked one far-IR band (SCUBA-2 850 μm for Millard et al. 2020, and ALMA 1.2 mm for Jolly et al. 2025). To convert the single far-IR photometry points to dust masses, we used the same modified blackbody function fit to our data and our dust temperature model as a function of t to scale the flux densities. The results are shown in the right panel of Fig. 6, where there is general agreement with our findings, although the dust mass-to-stellar mass ratios are lower at z < 1.5, likely due to the presence of quiescent galaxies in the comparison stacking catalogues. Interestingly, at higher redshifts (z > 2) where there is less contamination from quiescent galaxies, the stacked flux densities from Millard et al. (2020) suggest a plateauing mass ratio as opposed to a drop-off. This can be further investigated with future Euclid data releases once the deep fields get deeper, allowing us to probe higher redshifts with better statistical power.

The increasing dust temperatures with redshift have been attributed to the fact that MS galaxies have higher sSFRs at higher redshifts (e.g. Liang et al. 2019; Koprowski et al. 2024) and therefore contain more massive, young hot stars. This can be seen in our stacking results in Fig. 7, where we plot the dust temperature as a function of sSFR, separated by stellar mass in the same way as the previous plots. Galaxies with higher sSFRs clearly tend to have higher dust temperatures, although the trend depends on stellar mass and redshift. To understand why, we note that we can predict the sSFR–T_d relation for each stellar mass and redshift bin by combining the star-forming MS (Eq. (10)) and our fit to the dust temperature–t relation (Eq. (11)); see Appendix D for details. In Fig. 7 we show these predictions, where for each stellar mass bin we calculated the relation between z  =  0.2 (bottom-left of each curve) and z  =  3.0 (top-right of each curve). Again, we emphasise that these are not fits to the data shown in the plot, but predictions based on other fits. At a given redshift, the dust temperatures of all star-forming galaxies are independent of stellar mass. However, higher stellar mass galaxies have lower sSFRs (due to the bending of the MS), which leads to the mass dependence seen in Fig. 7. We note that the apparent convergence towards a constant sSFR at high redshift is simply a selection effect, since we measured the dust temperatures of lower mass galaxies out to higher redshifts than higher mass galaxies (see Fig. 3).

Fig. 7. Dust temperature as a function of sSFR from our best-fit far-IR SEDs, colour-coded to show the measured stellar mass dependence. Here, we only show physical properties for bins that are > 95% complete. The coloured curves show the predicted trends for the corresponding stellar masses by combining the star-forming MS from Popesso et al. (2023) with our fit to the dust temperature as a function of time, ranging from z  =  0.2 (starting at the bottom-left of each curve) to z  =  3.0 (ending at the top-right of each curve).

One feature worth noting is that according to the star-forming MS, the SFRs (and sSFRs) of galaxies continuously decrease to zero as a function of time, yet the dust temperatures do not. Instead, the average dust temperatures appear to be converging to a constant value of about 23 K, independent of stellar mass. Indeed, it seems reasonable to expect that the dust temperatures of MS galaxies are not strictly proportional to the SFR or sSFR, since this would imply that galaxies with no star formation would have zero dust temperature, which is impossible. Instead, Eq. (11) implies that as the SFRs of galaxies fall below a certain threshold, the dust is no longer heated by hot young stars but by the existing cooler and older stellar population, which changes on timescales much longer than the current age of the Universe. Indeed, Chapman et al. (2003) showed that IRAS-selected galaxies are bi-modal in luminosity-versus-far-IR colour (a proxy for temperature) with a break around 10 M_⊙ yr⁻¹, and discussed a possible transition from cirrus-dominated SEDs to SFR-dominated SEDs. Similarly, Groves et al. (2012) found that the mean dust temperature within M31 increases from about 17 K in the disc to 35 K in the bulge, despite the lack of any strong star-forming regions throughout the galaxy, suggesting that the heating is driven primarily by the density of the old stellar population.

5.2. Euclid’s contribution to the CIB

Looking at the CIB estimates in Table 1, the latest estimates cannot be said to be in close agreement; hence these CIB determinations are still dominated by systematic effects. This highlights the difficulty in subtracting zodiacal light and emission from the Milky Way, i.e. assessing the level of appropriate zero points when trying to determine the level of the extragalactic background. We can therefore only roughly estimate the fraction of the CIB resolved by the current Q1 Euclid catalogue: about 30–80% at 100 μm; 40–70% at 160 μm; 70–120% at 250 μm; 70–120% at 350 μm; 60–80% at 500 μm; and 30–40% at 850 μm.

A similar calculation was previously carried out using catalogues from the COSMOS field stacked on SPIRE images using the same stacking algorithm used here (Duivenvoorden et al. 2020). The authors found that r-band catalogues down to about magnitude 26 or K-band catalogues down to about magnitude 24 can recover essentially all of the CIB measured by FIRAS at all three SPIRE wavelengths. This provides a good benchmark for Euclid – the depth of the VIS images is about 24.7 in VIS and 23.2 in NISP, which approaches the depth of the COSMOS catalogues, and with a substantially larger numbers of catalogued objects.

In particular, at the SPIRE wavelengths we recovered > 60% of the CIB. These results are in line with Duivenvoorden et al. (2020), considering the r and K bands they used to create near-IR catalogues in their stacking study do not match perfectly with Euclid’s VIS and NISP instruments. It is reasonable to conclude that while we have not yet resolved the entire CIB with Euclid in Q1, > 60% is consistent. As the Euclid Deep Fields get deeper, we can therefore expect to approach a more complete recovery of the CIB through Euclid-selected galaxies.

5.3. Comparison to the MAMBO simulation

A number of simulations able to reproduce Euclid observations have been investigated, including the MAMBO (Mocks with Abundance Matching in BOlogna; Girelli et al. 2020) mock catalogue. The MAMBO catalogue is based on the Millennium Simulation (Springel et al. 2005) with physical properties prescribed using Empirical Galaxy Generator (EGG) code (Schreiber et al. 2017) across 3.14 deg². In particular, magnitudes in the Euclid bands, as well as far-IR flux densities at the Herschel and SCUBA-2 bands used here, were calculated from the simulation by Euclid Collaboration: Parmar et al. (2026) and used to construct a mock Euclid catalogue at the depth of the Q1 catalogue (equal to the depth of the Euclid Wide Survey). Simulated stacked far-IR flux densities of star-forming MS galaxies were then calculated within the same redshift and stellar mass bins used here, from which we fitted the same modified blackbody SEDs.

From this simulation we found slightly higher far-IR SFRs (although by a factor < 2) and dust masses (also larger by a factor < 2), with no significantly different redshift trends compared to our measurements. Interestingly, the simulated dust temperatures as a function of time are well-recovered by our pipeline, as show in Fig. 8. The simulation does not show the same plateauing behaviour as with our measurements, but instead decrease monotonically with time across the entire redshift range investigated here. This difference can be attributed to details of EGG, which assigns dust temperatures using an equation linear in redshift (Eq. (14) in Schreiber et al. 2017) as opposed to quadratic in redshift, which better matches our observations.

Fig. 8. Same as Fig. 5, but with the data points derived from fits to simulated stacked far-IR photometry in the MAMBO simulation (Euclid Collaboration: Parmar et al. 2026).

5.4. Star-formation rate density

Lastly, we calculated the far-IR-derived SFR density (SFRD) as a function of redshift coming from the Euclid catalogue. We multiplied the number of SPIRE stacking galaxies by the mean SFR, then summed the contributions from each stellar mass bin at a given redshift, and divided them by the volume of the redshift slice in the maps. Lastly, we averaged over every second redshift bin. In Fig. 9 we show the results for star-forming galaxies in purple, compared to several published estimates (Behroozi et al. 2013; Madau & Dickinson 2014; Koprowski et al. 2017). Since these literature curves include all galaxies (not just star-forming ones), we re-ran our stacking pipeline on the full Euclid catalogue using the same stellar mass and redshift bins, then re-fitted the same SEDs to derive the mean SFRs of all galaxies in each bin. We show the resulting SFRD for all Euclid galaxies as the green points.

Fig. 9. SFRD as a function of redshift solely for star-forming galaxies (purple) and the full Euclid catalogue (green). The black curves show published SFRD fits from Koprowski et al. (2017), Behroozi et al. (2013), and Madau & Dickinson (2014), with the latter scaled by 0.63 to convert from Salpeter to Chabrier initial mass functions. Beyond redshift 1.5 our estimates of the SFRD are incomplete because we are not able to recover enough far-IR photometry to fit SEDs.

We find that at z < 1.5 our total SFRD agrees well with the literature, but past this redshift our results are incomplete. This makes sense looking at Fig. 3, which shows that we are not able to recover enough far-IR photometry for high-stellar-mass galaxies at high redshifts to fit SEDs and derive SFRs. Future Euclid data releases will include more of these galaxies overlapping with more SPIRE data, allowing us to measure the complete SFRD past z  =  1.5. These SFRD values should be regarded as lower bounds, since stacking only recovers the mean flux of cataloged galaxies and can miss contributions from galaxies that were not detected or heavily obscured galaxies.

6. Conclusions

We stacked over 2 million star-forming galaxies from the Euclid Q1 catalogue across 17.6 deg² of far-IR imaging, providing robust statistics for their mean far-IR flux densities. We performed our stacking on Herschel-PACS 100- and 160-μm maps, Herschel-SPIRE 250-, 350-, and 500-μm maps, and SCUBA-2 850-μm maps. To avoid biases related to clustering, we used the SimStack algorithm, which simultaneously fits flux densities to far-IR beam-convolved model images of galaxy distributions in different bins.

Given the large number of galaxies available for stacking, we split the Euclid star-forming catalogue into eight stellar mass bins from log₁₀(M_*/M_⊙)  =  8.3 to log₁₀(M_*/M_⊙)  =  11.5, and 14 redshift bins from z  =  0.2 to z  =  3.0. Beyond these ranges, the Euclid stellar masses and redshifts are no longer reliable. We find significant stacked detections in most bins at all wavelengths.

Using these average flux densities, we modelled the average far-IR SEDs of Euclid galaxies using a modified blackbody function transitioning to a power law at high frequencies, finding good fits where we are able to measure the stacked flux densities. From our fits we derived average dust temperatures, dust masses, and far-IR-derived SFRs, and we find significant redshift evolution in all of these parameters. In particular, we investigate the difference between our far-IR-derived SFRs and the SFRs predicted from the star-forming MS. We find consistent values between the two estimates, with no significant trend in redshift. The average dust temperature decreases as a function of time following a functional form T₂  +  (T₁  −  T₂) e^−t/τ, with no stellar mass dependence. We argue that the dust temperatures of MS galaxies below z  =  1 have converged to a constant value (T  ≃  23 K) because the dust is now primarily heated by existing cooler and older stellar populations as opposed to hot young stars in star-forming regions. We also find that the average dust-to-stellar mass ratio increases for galaxies of all stellar mass up to z  ≃  1, and decreases with increasing stellar mass. This is a consequence of higher gas mass-to-stellar mass ratios at higher redshifts, although decreasing metallicity competes against this trend. Similarly, more massive galaxies have lower dust-to-stellar mass ratios due to their lower sSFRs compared to low mass galaxies. Lastly, we showed that the correlation between dust temperature and SFR (and therefore sSFR) is stellar mass-dependent since dust temperatures are stellar mass-independent.

We compared our results to a recent mock Euclid catalogue with derived far-IR photometry, finding good agreement for the SFRs and dust masses. However, we showed that the simulated catalogue predicts consistently decreasing dust temperatures below z  =  1, in disagreement with our observation. We attribute this discrepancy to the model used to produce the dust temperatures, which assigns dust temperatures using a monotonically decreasing function of redshift as opposed to a functional form which converges to a constant value at low z.

In the future, Euclid will observe more area overlapping with far-IR surveys and will obtain deeper VIS and NISP imaging in the Euclid Deep Fields, where some of the best existing far-IR and submillimetre imaging lies. These advances will provide even better statistics than are currently available, allowing our stacking analysis to extend to higher redshifts and lower stellar masses. In addition, upcoming far-IR facilities such as the Cerro Chajnantor Atacama Telescope (CCAT-Prime Collaboration 2023) will play an important role in better constraining the Rayleigh–Jeans tail of the average far-IR SEDs, improving the SED constraints and derived physical parameters.

Acknowledgments

This work has made use of the Euclid Quick Release Q1 data from the Euclid/mission of the European Space Agency (ESA), 2025, https://doi.org/10.57780/esa-2853f3b. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org). This work was supported by the Natural Sciences and Engineering Research Council of Canada and the Canadian Space Agency. This research used the Canadian Advanced Network For Astronomy Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and The Digital Research Alliance of Canada with support from the National Research Council of Canada, the Canadian Space Agency, CANARIE, and the Canadian Foundation for Innovation. The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen), Thales Alenia Space (Turin), and Astrium (Toulouse), with in excess of a hundred subcontractors. PACS was developed by a consortium of institutes led by MPE (Germany) and including: UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE was developed by a consortium of institutes led by Cardiff University (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). We acknowledge use of data obtained with the JCMT, which is operated by the EAO on behalf of NAOJ, ASIAA, KASI, and CAMS as well as the National Key R&D Program of China (No. 2017YFA0402700). Additional funding support is provided by the STFC and participating universities in the UK and Canada. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. ELSA: Euclid Legacy Science Advanced analysis tools (Grant Agreement no. 101135203) is funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or Innovate UK. Neither the European Union nor the granting authority can be held responsible for them. UK participation is funded through the UK Horizon guarantee scheme under Innovate UK grant 10093177.

References

Appendix A: PACS and SCUBA-2 field overviews

In Fig. A.1 we show the Herschel-PACS 100- and 160-μm data and RMS maps overlapping with the EDF-F and EDF-N. We also show the SCUBA-2 data and RMS map overlapping with the EDF-N in Fig. A.2. In both figures the Euclid mask is shown as the grey contours. The bright resolved source in the EDF-N is NGC 6543 (the ‘Cat’s Eye Nebula’), which we mask prior to the stacking.

Fig. A.1. Top: Same as Fig. 1, but for the PACS observations of the CDFS-SWIRE (EDF-F) field. Bottom: Same as the above panel, but for the PACS observations of the AKARI-NEP (EDF-N) field.

Fig. A.2. Same as Fig. 1 but for the SCUBA-2 observations of the AKARI-NEP (EDF-N) field.

Appendix B: Stacking cutouts and flux densities

Here we show the 2D cross-correlations from our stacking algorithm and provide a table of stacked flux densities. The 2D stacking results for the star-forming galaxies in the Euclid catalogue are shown in Fig. B.1 (for SPIRE), Fig. B.2 (for PACS), and Fig. B.3 (for SCUBA-2). We show both the signal (in units of mJy) in the left column and the S/N in the right column, and the results for the mask and the remaining galaxies in the Euclid catalogue (i.e., non-star-forming and outside the main redshift and stellar ranges we considered here) are shown in the top row. The corresponding stacked flux densities (the value of the central pixel in the 2D cross-correlations) are provided in Table B.1. Bins that are >  95% complete in stellar mass (Euclid Collaboration: Enia et al. 2026) are highlighted in blue. The 2D cross-correlation profiles are expected to follow the autocorrelations of the instrumental beams. We tested this by computing the averaged 1D radial profiles and compared these to the expected PSF profiles. We found good agreement between the stacked signals and the PSFs for most bins with log₁₀(M_*/M_⊙) > 9.9. Below this stellar mass we found that the 2D profiles can be more extended than the beam, with the largest effect seen at 500 μm where the PSF is largest. However, this does not affect the calculations and results in this work so we do not attempt to correct for it.

Fig. B.1. Herschel-SPIRE results from stacking star-forming MS galaxies from the Euclid catalogue (Euclid Collaboration: Enia et al. 2026) with redshifts between 0.2 and 3.0, and stellar masses between log10(M*/M⊙)  =  8.3 and 11.5. Each 2D cutout is 200″  ×  200″. Bins that are >  95% complete in stellar mass (Euclid Collaboration: Enia et al. 2026) are highlighted in blue. Left column: Stacking signal in units of mJy. Right column: The S/N of the stacked flux densities. The ‘Mask’ and ‘Other gals.’ panels differ from the colour bar and range from S/N  =   − 1 to S/N  =  60. Top row: SPIRE 250-μm. Middle row: SPIRE 350-μm. Bottom row: SPIRE 500-μm.

Fig. B.2. Same as Fig. B.1 but for PACS 100 (top row) and 160 μm (bottom row). Here, each 2D cutout measures 40″  ×  40″. The ‘Mask’ and ‘Other gals.’ panels differ from the colour bar and range from S/N  =   − 1 to S/N  =  20.

Fig. B.3. Same as Fig. B.1, but for SCUBA-2 850 μm. Here, each 2D cutout measures 70″  ×  70″.

Appendix C: Best-fit far-IR SED parameters and derived quantities

Here we show the resulting physical parameters T_d, M_d, and SFR in each of the stellar mass and redshift bins where we have sufficient stacked photometry to derive these physical parameters, and we provide the best-fit SED parameters and derived physical quantities. Fig. C.1 shows the physical parameters for each redshift and stellar mass bin in our stacking analysis (where bins that are >  95% complete in stellar mass are highlighted in blue – see Euclid Collaboration: Enia et al. 2026.), and best-fit far-IR SED parameters and derived parameters are provided in Table C.1. We omit showing the lowest stellar mass bins since they contain no data.

Fig. C.1. Physical parameters derived from the best-fit modified blackbody SEDs in Fig. 3. Bins that are >  95% complete in stellar mass (Euclid Collaboration: Enia et al. 2026) are highlighted in blue. Top: Best-fit dust temperatures, Td. Middle: Dust mass (Md), calculated by scaling the best-fit amplitude (see e.g. Reuter et al. 2020; Eales & Ward 2024; Jolly et al. 2025). Bottom: SFRs calculated from LIR (the integral of the best-fit SED from 8 to 1000 μm) multiplied by a factor of 1.49  ×  10−10 M⊙ yr−1 L${}_{\odot }^{-1}$.

Appendix D: Mathematical models combining the star-forming main sequence with dust properties

In this appendix we outline the mathematical details for combining the star-forming MS with the dust temperature, T_d, and the dust mass, M_d. We start by re-defining the MS, explicitly showing the dependence on time and stellar mass; here we focus on the parametrisation from Popesso et al. (2023), but it is easy to adjust our derivations for other parametrisations. The MS follows

$$\begin{array}{c}\hfill \mathrm{SFR}(t,M_{\ast })=\frac{{\mathrm{SFR}}_{\max }(t)}{1+{(M_0(t)/M_{\ast })}^{\gamma }},\end{array}$$(D.1)

where we use γ  =  1, log₁₀(SFR_max)(t)  =  a₀  +  a₁t, log₁₀(M₀)(t)  =  a₂  +  a₃t, a₀  =  2.71, a₁  =   − 0.186, a₂  =  10.86, and a₃  =   − 0.0729. Our equation for the dust temperature is only a function of time and not stellar mass and follows

$$\begin{array}{c}\hfill T_{\mathrm{d}}(t)=T_2+(T_1-T_2){\mathrm{e}}^{-t/\tau },\end{array}$$(D.2)

where T₁  =  79.7 K, T₂  =  23.2 K, and τ  =  1.6 Gyr. The SFR is related to the time-evolving far-IR SED through

$$\begin{array}{c}\hfill \mathrm{SFR}(t,M_{\ast })=(1.49\times {10}^{-10})4\pi D_{\mathrm{L}}^2(t){\displaystyle {\int }_{{\nu }_1}^{{\nu }_2}}{S}_{\nu }(A(t,M_{\ast }),T_{\mathrm{d}}(t))\mathrm{d}\nu ,\end{array}$$(D.3)

where ν₁  =  c/1000 μm, ν₂  =  c/8 μm, and S_ν is the rest-frame SED. The dust mass is calculated using similar quantities, with

$$\begin{array}{c}\hfill M_{\mathrm{d}}(t,M_{\ast })=\frac{D_{\mathrm{L}}^2(t)A(t,M_{\ast })}{{\kappa }_0}.\end{array}$$(D.4)

We note that we are assuming only two free parameters in the time-evolving far-IR SED, the overall multiplicative amplitude A (a function of both time and stellar mass) and the dust temperature T_d (a function of only time) – see Sect. 4.2.

We start by solving for M_d(t). This effectively amounts to calculating the ratio between Eqs. D.1 and D.3:

$$\begin{array}{c}\hfill A(t,M_{\ast })=\frac{{\mathrm{SFR}}_{\max }(t)/(1+{(M_0(t)/M_{\ast })}^{\gamma })}{(1.49\times {10}^{-10})4\pi D_{\mathrm{L}}^2(t){\displaystyle {\int }_{{\nu }_1}^{{\nu }_2}}{S}_{\nu }(1,T_{\mathrm{d}}(t))\mathrm{d}\nu }.\end{array}$$(D.5)

This is then inserted into Eq. D.4 to find M_d(t), and the dust-to-stellar mass ratio can be calculated by dividing by the stellar mass.

We can also calculate time-dependent correlations between quantities. As an example we show how to derive the equations relating SFR and T_d, and we note that the same procedure can produce equations relating any other combination of far-IR properties. Eq. D.1 says that SFRs increase monotonically with time for galaxies of all stellar mass; therefore, given an SFR and a stellar mass we can calculate the appropriate age of the Universe. While in the current mathematical form this cannot be done analytically, a simple root-finding algorithm will provide a function f that provides the age for a given SFR and stellar mass, t  =  f(SFR, M_*). The dust temperature can then be written as

$$\begin{array}{c}\hfill T_{\mathrm{d}}(\mathrm{SFR},M_{\ast })=T_2+(T_1-T_2){\mathrm{e}}^{-f(\mathrm{SFR},M_{\ast })/\tau },\end{array}$$(D.6)

where the SFR range being plotted is understood to correspond to a time range. This can be converted to T_d as a function of sSFR (SFR/M_*) by dividing Eq. D.1 by M_* and writing a new root-finder to obtain a function g that provides the age for a given sSFR and stellar mass, t  =  g(sSFR, M_*).

All Tables

All Figures

Fig. 1. Top: Herschel-SPIRE data covering the CDFS-SWIRE field (overlapping with the EDF-F) at 250, 350, and 500 μm. The blue contour shows the mask applied to the SPIRE images to remove bad edge pixels. The grey contours show the corresponding Euclid catalogue mask, where masked rectangles designate the locations of bright stars in the field contaminating source extraction. Bottom: Same as the top panel, but showing the RMS of the Herschel-SPIRE data. Coordinates are conventional RA and Dec.

In the text

	Fig. 2. Same as Fig. 1, but for the Herschel-SPIRE AKARI-NEP field (overlapping with the EDF-N).
In the text

Fig. 3. Modified blackbody SEDs (with β = 1.96 and α = 2.3) fitted to the stacked Herschel and SCUBA-2 flux densities. The redshift and stellar mass of each bin are indicated by the top and right axis labels, respectively. The best-fit parameters are given in Table C.1. The SEDs were fitted to bins where all three SPIRE flux densities are detected with S/N> 3 and at least one PACS flux density is detected with S/N> 3. The panels are otherwise blank. Bins that are > 95% complete in stellar mass are highlighted in blue.

In the text

Fig. 4. Ratio of our SFRs measured from far-IR photometry to a parametrisation of the star-forming MS, shown as a function of redshift and split into different stellar mass stacking bins. Only redshift and stellar mass bins with > 95% completeness are shown. We use the MS parametrisation from Popesso et al. (2023), which is a continuous function of z/t and is in good agreement with Euclid Collaboration: Enia et al. (2026).

In the text

Fig. 5. Best-fit dust temperatures from our SED fitting, Td, as a function of time (bottom axis) and redshift (top axis), considering only the redshift and stellar mass bins that are > 95% complete. We show the dust temperature evolution for five different stellar mass bins, with the stellar mass values of the centres of the bins given in the legend. The solid curve is a fit to the simple form T2  +  (T1  −  T2) e−t/τ; the dotted line is the quadratic-in-redshift fit from Koprowski et al. (2024) and the dashed line is the linear-in-redshift fit from Schreiber et al. (2017). We also show published mean temperature estimates for star-forming galaxies at low redshifts (0.01 < z < 0.05; Lamperti et al. 2019) in blue.

In the text

Fig. 6. Left: Best-fit dust masses from our SED fits, Md, as a function of redshift, considering only the redshift and stellar mass bins that are > 95% complete. Right: Same as the left panel but showing the best-fit dust mass-to-stellar mass ratio, where the stellar masses are from the Euclid catalogue and are used to define the stacking bins. We also show SCUBA-2 stacking results from Millard et al. (2020), scaled to dust mass assuming the same modified blackbody SEDs used here and our best-fit dust temperature as a function of time, and ALMA 1.2 mm stacking results from Jolly et al. (2025) scaled in the same way. In both panels, the solid-coloured curves show the predicted trends for the corresponding stellar masses by combining the star-forming MS from Popesso et al. (2023) with our fit to the dust temperature as a function of time.

In the text

Fig. 7. Dust temperature as a function of sSFR from our best-fit far-IR SEDs, colour-coded to show the measured stellar mass dependence. Here, we only show physical properties for bins that are > 95% complete. The coloured curves show the predicted trends for the corresponding stellar masses by combining the star-forming MS from Popesso et al. (2023) with our fit to the dust temperature as a function of time, ranging from z  =  0.2 (starting at the bottom-left of each curve) to z  =  3.0 (ending at the top-right of each curve).

In the text

	Fig. 8. Same as Fig. 5, but with the data points derived from fits to simulated stacked far-IR photometry in the MAMBO simulation (Euclid Collaboration: Parmar et al. 2026).
In the text

Fig. 9. SFRD as a function of redshift solely for star-forming galaxies (purple) and the full Euclid catalogue (green). The black curves show published SFRD fits from Koprowski et al. (2017), Behroozi et al. (2013), and Madau & Dickinson (2014), with the latter scaled by 0.63 to convert from Salpeter to Chabrier initial mass functions. Beyond redshift 1.5 our estimates of the SFRD are incomplete because we are not able to recover enough far-IR photometry to fit SEDs.

In the text

	Fig. A.1. Top: Same as Fig. 1, but for the PACS observations of the CDFS-SWIRE (EDF-F) field. Bottom: Same as the above panel, but for the PACS observations of the AKARI-NEP (EDF-N) field.
In the text

	Fig. A.2. Same as Fig. 1 but for the SCUBA-2 observations of the AKARI-NEP (EDF-N) field.
In the text

Fig. B.1. Herschel-SPIRE results from stacking star-forming MS galaxies from the Euclid catalogue (Euclid Collaboration: Enia et al. 2026) with redshifts between 0.2 and 3.0, and stellar masses between log10(M*/M⊙)  =  8.3 and 11.5. Each 2D cutout is 200″  ×  200″. Bins that are >  95% complete in stellar mass (Euclid Collaboration: Enia et al. 2026) are highlighted in blue. Left column: Stacking signal in units of mJy. Right column: The S/N of the stacked flux densities. The ‘Mask’ and ‘Other gals.’ panels differ from the colour bar and range from S/N  =   − 1 to S/N  =  60. Top row: SPIRE 250-μm. Middle row: SPIRE 350-μm. Bottom row: SPIRE 500-μm.

In the text

	Fig. B.2. Same as Fig. B.1 but for PACS 100 (top row) and 160 μm (bottom row). Here, each 2D cutout measures 40″  ×  40″. The ‘Mask’ and ‘Other gals.’ panels differ from the colour bar and range from S/N  =   − 1 to S/N  =  20.
In the text

	Fig. B.3. Same as Fig. B.1, but for SCUBA-2 850 μm. Here, each 2D cutout measures 70″  ×  70″.
In the text

Fig. C.1. Physical parameters derived from the best-fit modified blackbody SEDs in Fig. 3. Bins that are >  95% complete in stellar mass (Euclid Collaboration: Enia et al. 2026) are highlighted in blue. Top: Best-fit dust temperatures, Td. Middle: Dust mass (Md), calculated by scaling the best-fit amplitude (see e.g. Reuter et al. 2020; Eales & Ward 2024; Jolly et al. 2025). Bottom: SFRs calculated from LIR (the integral of the best-fit SED from 8 to 1000 μm) multiplied by a factor of 1.49  ×  10−10 M⊙ yr−1 L${}_{\odot }^{-1}$.

In the text