Euclid view of dusty star forming galaxies at z>~1.5 detected in wide area submillimetre surveys

We investigate the constraints provided by the Euclid space observatory on the physical properties of dusty star forming galaxies (DSFGs) at z>~1.5 detected in wide area sub millimetre surveys with Herschel. We adopt a physical model for the high z progenitors of spheroidal galaxies, which form the bulk of the DSFGs at z>~1.5. We improve the model by combining the output of the equations of the model with a formalism for the spectral energy distribution(SED). After optimising the SED parameters to reproduce the measured infrared luminosity function and the number counts of DSFGs, we simulated a sample of DSFGs over 100 sq deg and then applied a 5 sigma detection limit of 37 mJy at 250 microns. We estimated the redshifts from the Euclid data and then fitted the Euclid and Herschel photometry with the code CIGALE to extract the physicsl parameters. We found that 100 % of the Herschel galaxies are detected in all 4 Euclid bands above 3 sigma. For 87% of the sources the accuracy on 1+z is better than 15%. The sample comprises mostly massive log(Mstar/Msun)~10.5-12.9, highly star forming, log(SFR/Msun/yr)~1.5-4, dusty, log(Mdust/Msun)~7.5-9.9, galaxies. The measured stellar mass have a dispersion of 0.19 dex around the true value, thus showing that Euclid will provide reliable stellar mass estimates for the majority of the bright DSFGs at z>~1.5 detected by Herschel. We also explored the effect of complementing the Euclid photometry with that from Vera C. Rubin Observatory/LSST.


INTRODUCTION
In the 1990s, surveys using the Submillimeter Common User Bolometer Array (SCUBA) at 850 m on the 15 m James Clerk Maxwell Telescope (JCMT; Smail et al. 1997), for the first time detected a population of high-redshift ( ≳ 1.5) galaxies which are extremely luminous at far-infrared (FIR) and submillimeter (sub-mm) wavelengths and are highly dust-obscured at optical wavelengths.These are now commonly referred to as dusty starforming galaxies (DSFGs) or sub-millimeter galaxies (SMGs) and are the dominant star-forming galaxies at  ≳ 1.5 (Casey et al. 2012;Bothwell et al. 2013;Cai et al. 2013) and most likely to be in the early evolutionary phase of present-day elliptical galaxies (Ivison et al. 2013;Fu et al. 2013;Dokkum et al. 2015;Simpson et al. 2017).These galaxies undergo intense star formation with star formation rates (SFRs) from a few hundreds to thousands solar masses per year (Blain et al. 2002;Magnelli et al. 2012;Riechers et al. 2013;Swinbank et al. 2014;MacKenzie et al. 2017;Michałowski et al. 2017;Gullberg et al. 2019;Castillo et al. 2023).During the cosmic noon (at  ∼ 2) (Madau et al. 1996;Lilly et al. 1996) when the star formation rate density (SFRD) of the Universe peaked, these galaxies contributed heavily to the overall galaxy population ★ E-mail: MitraD@cardiff.ac.uk (Magnelli et al. 2011a;Burgarella et al. 2013;Rowan-Robinson et al. 2016).Therefore, the discovery of these galaxies caused a big revolution in the field of extragalactic astronomy making their study vital for a detailed understanding of the formation and evolution of galaxies which led to the development of more sensitive instruments like the SCUBA-2 (Holland et al. 2013), the Large Apex BOlometer CAmera (LABOCA; Siringo et al. 2009),the AZtronomical Thermal Emission Camera (AzTEC; Wilson et al. 2008) and the Herschel Space Observatory (Herschel; Pilbratt et al. 2010).
Our understanding of the DSFGs has increased drastically thanks to Herschel which, during its four years of operation, from 2009-2013, has mapped ∼ 1300 sq.deg. of the sky at wavelengths in the range 100-500 m leading to the detection of more than a million FIR and sub-mm bright galaxies.In particular, the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) and the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012) together surveyed over ∼ 930 sq.deg. of the sky at 250, 350 and 500 m.Despite such surveys, our understanding of the underlying physical processes responsible for the evolution of galaxies is far from complete, but it is speculated that there will be a clearer understanding of these processes in the next decade with big projects like the Euclid space observatory (Laureĳs et al. 2011), which was launched on 1st July 2023, and the Vera C. Rubin Observatory/LSST (Ivezić et al. 2019).The Euclid Wide Survey (EWS; Euclid Collaboration 2022b) will observe 15 000 sq. deg. of the sky at both optical (Cropper et al. 2018) and near-infrared (NIR; Euclid Collaboration 2022a) wavelengths, thus providing crucial constraints on the stellar mass of galaxies.On the other hand, Rubin's LSST (Legacy Survey of Space and Time) Survey will observe 18 000 sq. deg. of the sky in 6 bands -, , , , , and  -having a wavelength coverage of 0.36 − 1.01 m.
Many sophisticated phenomenological models have been developed to study the cosmological evolution of galaxies (Matteucci et al. 2009;Franceschini et al. 2010;Popesso et al. 2011;Béthermin et al. 2012).These models adopt simple functional forms to describe the luminosity function of multiple galaxy populations, which are assigned different spectral energy distributions (SEDs) and varying evolutionary properties.Active galactic nuclei (AGNs) are also considered in some of the models.In this paper, we adopt instead the physically motivated model by Cai et al. (2013, C13 hereafter) for the formation and evolution of  ≳ 1.5 DSFGs, which links the star formation activity and the growth of a central supermassive black hole in a self-consistent way.The C13 model has proved successful in reproducing data on luminosity functions, number counts and redshift distributions over a broad wavelength range, from the mid-IR to the millimetre.We use the model to investigate how Euclid can detect and study the  ≳ 1.5 proto-spheroidal galaxies discovered by the H-ATLAS survey.For this purpose, we improve on the C13 original model by implementing a more sophisticated formalism for the SED, which allows us to predict the optical/near-IR properties of these galaxies.
The paper is organised in the following manner.In Section 2, we provide a brief description of the model and theory used in this paper.In Section 3, we present a detailed methodology of catalogue creation, deriving photometric flux densities and estimating photometric redshifts and other physical properties using SED fitting.In Section 4, we discuss the surveys for which we make forecasts.The results obtained are analysed and discussed in Section 5 and the final conclusions are summarized in Section 6.Throughout the paper, we adopt a flat Λ-CDM cosmology with present-day matter and baryon density (in units of the critical density), Ω ,0 = 0.3153 and Ω ,0 = 0.0493.We set the value of the Hubble-Lemaître constant to ℎ =  0 /100 = 0.6736, the slope of the spectrum of primordial density perturbations to  = 0.9649 and the normalization of the density fluctuations on a scale of 8ℎ −1 Mpc to  8 = 0.8111 (Planck Collaboration VI, 2020).

ADOPTED MODEL FOR THE PROTO-SPHEROIDAL GALAXIES
In this section, we provide a brief description of the C13 model for the proto-spheroidal galaxies and illustrate the formalism used for predicting the SED of these galaxies.

Model Outline
The model by C13 is based on the fact that spheroids in the local universe (i.e.ellipticals and bulges of disc galaxies) are mostly comprised of old stellar populations that formed at redshifts  ≳ 1.5.On the other hand, disc galaxies consist of comparatively younger stellar populations with luminosity-weighted age ≲ 7 Gyr, indicating a formation redshift  ≲ 1.Therefore, the progenitors of local spheroids (also called proto-spheroidal galaxies or proto-spheroids) are the dominant star-forming galaxies at  ≳ 1.5, whereas star formation in discs takes place mostly at  ≲ 1.5.The model provides a physically motivated co-evolution of the SFR of the proto-spheroids along with the supermassive black hole at their centres.The star formation history of these proto-spheroids is calculated using a set of equations considering the evolution of gas phases and of the AGN, including cooling, condensation to form stars, accretion into the supermassive black hole, and feedback from both the stellar component and the AGN (see Figure 1).The model depicts two ways of dark matter halo formation, i.e. a fast collapse of the halo triggering star formation, followed by a slow growth of the halo outskirts, which have little effect on the inner region of the halo.The first phase is the one driving the build-up of most of the stellar mass in these objects.
Both the star formation and the growth of the super-massive black hole at the centre are controlled by feedback mechanisms from the supernovae and the AGN.The AGN feedback is effective in quenching the star formation in the most massive proto-spheroids, while SFR is mainly affected by SN feedback in less massive ones (Granato et al. 2004;Scannapieco et al. 2008;Rosas-Guevara et al. 2016;Rosito et al. 2021).A detailed explanation of the equations governing the evolution of both the stellar and the AGN components is given in Appendix A.
The bolometric luminosity function of the proto-spheroidal galaxies, i.e. the volume density of galaxies per unit interval of luminosity, is calculated by convolving the halo formation rate function with the galaxy luminosity distribution predicated by the equations of the model, i.e.

Φ(log
where  ST is the halo mass function, for which the model adopts the Sheth & Tormen (1999) approximation, while (log , ;  vir ,  vir ) is the luminosity distribution function of galaxies at redshift  with virialisation mass and redshift given by  vir and  vir respectively.The distribution function is a log-normal distribution given by where  is the dispersion around the mean luminosity, L, which is computed using the formalism illustrated in Sec.2.2 for the stellar component.The AGN luminosity is computed using Eq.(A17).The total luminosity of a galaxy is the sum of the stellar luminosity and the AGN luminosity.The dispersion is considered to reflect the uncertainties in the adopted values of the main parameters that are involved in the equations of the model.

SED modelling
The SED of a galaxy is an observable that reflects the properties of the stars, dust and gas in the galaxy.As such, it can be used to extract important information.In fact, while the UV/optical/near-IR part of the SED depends on the star formation history, metallicity and dust extinction, the far-IR/sub-mm/mm portion of the SED reveals the star formation activity occurring in the dust-enshrouded birth clouds and can be used to constrain the properties of the dust itself.C13 adopted a single representative SED for the whole population of high- proto-spheroidal galaxies, i.e. the SED of the well-studied  = 2.3 lensed galaxy SMM J2135-0102.That choice was dictated by the need of speeding up the computation of the relevant statistical properties on the galaxies, such as number counts and redshift In each plot, the left y-axis refers to the mass evolution of the stellar component: infalling gas (solid red line), cold gas (dot-dashed black line), and stellar mass (blue dashed line); while the right y-axis refers the mass evolution of AGN component: gas falling in the reservoir (solid green line) and super-massive black hole (SMBH) mass (dotted magenta line).We observe that the stellar mass increases as more and more cold gas is condensed into stars, thus leading to a corresponding steady decline in the mass of the gas components.At the same time, the central BH mass increases as a fraction of the cold gas loses its angular momentum via the stellar radiation drag and settles down into a reservoir around the central SMBH.The BH accretes matter from the reservoir via viscous dissipation, which powers the nuclear activity.With the evolution of the galactic age, the AGN feedback sets in.Its effect can be observed in the most massive dark matter (DM) halos where the star formation activity and, consequently, the stellar mass growth, are stopped within less than 1 Gyr from the galaxy's birth.
distributions.Also, C13 were mainly focused on observables at far-IR/sub-mm/mm wavelengths where the SED depends on a very limited number of free parameters, e.g. the dust temperature and the dust emissivity index.In fact, with the SED of SMM J2135-0102, C13 did achieve a good qualitative fit to the multi-wavelength (from the mid-IR to millimeter waves) data on number counts, both global and per redshift slices (e.g.Béthermin et al. 2012).Moreover, the model made a successful prediction of the sub-mm (Negrello et al. 2010(Negrello et al. , 2017) ) and mm (Vieira et al. 2010;Cai et al. 2020Cai et al. , 2022) ) number counts of strongly lensed galaxies.
However, the use of a single SED does limit the applicability of the model to the investigation and the forecast of the statistical properties of dusty proto-spheroidal galaxies at optical/near-IR wavelengths, where the SED is much more complex than at longer wavelengths.Therefore, we have decided to improve the C13 model by implementing an SED modelling that reflects such complexity.The SED formalism is based on that put forward by da Cunha et al. (2008), with a few modifications, and it is described in this section.

Stellar component and dust attenuation
The stellar emission spectrum of the proto-spheroidal galaxies is calculated by coupling the star formation history (SFH, hereafter), as obtained from the equations of the model (Eq.(A7)) for any given value of  vir and  vir , with some single stellar population (SSP) synthesis templates.There are several SSP templates available in literature (Bruzual & Charlot 2003;Jimenez et al. 1998;Buzzoni 1993).Here we use the SSP templates by Bruzual & Charlot (2003), which are provided at wavelengths from 91 Å to 160 m and between ages 1 × 10 5 yr and 2 × 10 10 yr, for metallicities  = 0.008, 0.02 and 0.05 and three IMFs, i.e.Salpeter, Korupa and Chabrier.We adopt the SSP template corresponding to  ⊙ and a Chabrier IMF.Dust attenuation is modelled using the attenuation law by Charlot & Fall (2000).It takes into account attenuation by both birth clouds (BC) and the interstellar medium (ISM), represented by different power laws.This also helps in differentiating between the emission from the young stars still embedded in their dust-enshrouded birth clouds and the old stars located in the ISM, respectively, which is important for IR galaxies (Małek et al. 2018;Buat et al. 2018).
The luminosity emitted by a galaxy per unit interval of wavelength at time  is given by (Charlot & Fall 2000) Taking into consideration the AGN feedback which quenches star formation (see Figure 2), we modify and re-define   ( ′ ) as where  AGN is the time at which the AGN feedback has completely stopped the star formation.In practice, when this happens, all the gas and dust in the galaxy have been removed by the AGN feedback so that τ ( ′ ) = 0.
The amount of starlight absorbed by dust present in the BCs and in the ISM is reprocessed and re-radiated at infrared wavelengths.The dust luminosity is thus given by where and denote the amount of starlight absorbed and re-radiated by dust in the BCs and in the ISM, respectively.For later use in the analysis, we also define the fraction of IR luminosity coming from the dust in the ISM as Clearly, this fraction depends on τ  ,  and the SFH.

Infrared emission by dust
Here, we briefly describe the way in which we compute the infrared part of the SED, which is associated with the dust emission.We assume three main dust components: Polycyclic Aromatic Hydrocarbons (PAHs), which produce strong emission in the near-infrared (NIR) at wavelengths of about 3 − 20 m, small dust grains (size < 0.01 m), which are responsible for the hot mid-infrared (MIR) emission, and large dust grains (size ∼ 0.01 − 0.25 m) which are in thermal equilibrium with the radiation field and produce the colder FIR emission.
(i) PAH emission: in star-forming galaxies, strong emission features at 3.3, 6.2, 7.7, 8.6, 11.3 and 12.7 m are observed.They are attributed to emission from PAH molecules.In the model, the PAH emission is represented by the spectrum of the photo-dissociation region (PDR) in the star-forming region M17 SW from ISOCAM observation of Cesarsky et al. (1996) and extracted by Madden et al. (2006).In order to include the 3.3 m PAH emission feature in the template spectra we extend it blueward beyond 5 m using a Lorentzian profile by Verstraete et al. (2001).The SED of the PAHs in the model is calculated as where,  ,17 is the template spectra used here.
(ii) MIR emission: Besides the emission from PAH molecules, the spectra of star-forming galaxies consist of a MIR continuum emission due to the heating of small dust grains to very high temperatures.This part of the spectra is computed using a 'greybody' function as where,   ( dust ) is the Planck function at temperature  dust and   is the dust mass absorption coefficient of the form where  is the dust emissivity index.Studies show that  ≈ 1 for carbonaceous grains (radiating energy in the MIR) while  ≈ 1.5 − 2 for silicate grains (Hildebrand 1983;Draine & Lee 1984), radiating most of the energy in the far-IR/submm wavelengths.In the model, the MIR continuum emission is characterised by the sum of two greybody functions at temperatures 130 K and 250 K respectively, having equal contribution to the IR luminosity.Mathematically, where  ,130 and  ,250 are calculated using Eq. ( 13).The value of  for the MIR emission is set to 1.
(iii) Dust grains in equilibrium: the FIR portion of the SED comprises emission from bigger dust grains which are in thermal equilibrium and have comparatively low temperatures.The model takes into account two types of grains: warm grains which can be present in BCs and in the ISM, with temperatures  , and  ,  respectively; and cold grains which are found in the ISM at temperatures  ,  .The emissions are computed using greybody functions (Eq.( 13)) with  = 1.5 and 2.0 respectively.
The total IR emission from BCs writes satisfying the condition Likewise, the ISM emission can be computed as satisfying the condition The adopted values of the fractions mentioned in Eqs. ( 17) and ( 19) are discussed in Section 3.3.Fig. 3 shows an example of the different components that make up the modelled SED of a proto-spheroidal galaxy.

Emission associated with the AGN
The SED of the AGN component is modelled using the smooth torus model introduced by Fritz et al. (2006, hereafter F06).We considered a smooth torus model with the ratio of maximum to minimum radii  max / min = 10 to 100.The torus amplitude angle, Θ, has values of 40 • , 100 • and 140 • .The variation of the gas density inside the torus along the radial and angular direction is given by the equation (, ) =    − | cos(  ) | , where  = 0 and  = 0 mean constant density.The model assumes that the luminosity of the central power source of the torus is described by a broken power law as   () =  0   (erg/s), where  = 1.2 if 0.001m <  < 0.03m,  = 0 if 0.03m <  < 0.125m and  = −0.5 if 0.125m <  < 20m.The equatorial optical depth at 9.7 m,  eq (9.7), ranges from 0.1 to 10.The model considers 10 viewing angles ranging from 0 • (type 1 AGN) to 90 • (type 2 AGN).For more details about the model we refer the reader to F06.
In our simulations, as explained in Sec. 3, we associate to each simulated galaxy an AGN SED by randomly sampling the available range of values of the SED parameters.We then normalise the SED so that the resulting bolometric AGN luminosity equals the one obtained from Eq. (A17).Some of the light emitted from the AGN gets absorbed by the dust in the ambient ISM.So, the AGN luminosity emitted per unit wavelength at time  is given by where   is the luminosity per unit wavelength interval from the template.We define τ,AGN The total AGN luminosity absorbed by dust in the ISM is then given by Therefore, the total dust luminosity of the galaxy is redefined as and the fraction of the dust luminosity coming from the ISM is now The absorbed AGN luminosity is reprocessed by dust and emitted at IR wavelengths.Likewise, the total IR emission from the BC and the ISM is now re-defined as and respectively.Fig. 4 shows the F06 AGN template SED as a function of wavelength for different viewing angles (from 0 • to 90 • ) and a constant density profile of the gas in the torus.An example of the SED of a proto-spheroidal galaxy with a visible AGN component is shown in Fig. 5.
Apart from the dusty torus and the dust in the ISM, "polar dust" can also contribute to the emission from the AGN.This dust extends in the polar direction of the dusty torus over parsec scales (Hönig et al. 2013;Yang et al. 2020).Its emission peaks in the mid-IR (Hönig et al. 2013, see Figure 8), and it is negligible at the wavelengths of interest here, i.e.UV/optical/near-IR wavelengths.Hence, we do not include the effect of polar dust in our model.

Dust mass estimation
To estimate the dust mass in galaxies, we follow the prescriptions of da Cunha et al. (2008), which are based on the formula by Hildebrand (1983).For dust grains in thermal equilibrium at the temperature   the dust mass (  ) is given by where    is the infrared luminosity of the dust grains, while   (  ) and   are defined in Eq. ( 13) and ( 14), respectively.After normalising   using  850  = 0.77 g −1 cm 2 , the mass contribution of warm dust in BC and ISM is calculated as and respectively, and that from the cold dust in the ISM is calculated as Taking into account the contribution from stochastically heated dust grains and considering only a small contribution of a few per cent from PAH grains (Draine & Li 2007), the total dust mass of a galaxy is estimated as

BUILDING A SIMULATED SAMPLE OF PROTO-SPHEROIDAL GALAXIES
Here we describe the way in which we implement the formalism presented in the previous section to generate a sample of simulated proto-spheroidal galaxies.The sample will then be used for forecast studies.

Step 1: sampling the halo formation rate function
We start by randomly sampling the halo formation rate function (Eq.( A25)) in both halo virialization mass and virialization redshift, within the intervals  vir ∈ [1.5, 12] and log( vir / ⊙ ) ∈ [11.3, 13.3], respectively.The total number of randomly sampled pairs of  vir and  vir is dictated by the adopted survey area, which, in turn, determines the sampled volume at any given virialization redshift.For each simulated pair of values of  vir and  vir we then solve the C13 equations to compute the SFH as a function of the galaxy age, with the birth of the galaxy set at the corresponding virialization redshift.

Step 2: generating the UV/optical/near-IR part of the SED
At this point we apply the SED formalism illustrated in Sec.2.2.1 to compute the SED of the dust-attenuated stellar component (Eq.( 3)) as well as the SED of the dust-attenuated AGN (Eq.( 20)) at any given age of the galaxy.The main parameters affecting the stellar component are the age of the BCs,  0 , and the total effective Vband absorption optical depth of the dust in the BCs, τ , and in the ISM, τISM  .For them, we adopt the values measured by Rowlands et al. (2014, R14 hereafter), who performed SED fitting on a sample of 29 250 m-selected high-redshift DSFGs from H-ATLAS, using MAGPHYS (da Cunha & Charlot 2011).They obtained a good fit by allowing  0 to vary uniformly in logarithmic space between 10 7 −10 8 years, while for τ  and τISM  they derived mean values of 5.1 ± 0.6 and 1.0 ± 0.1, respectively.In practice, we take Gaussian priors for τBC  and τISM  with mean 5.1 and 1.0, respectively.The standard deviation is taken as 0.6 and 0.1, respectively, for all the simulated objects, and the value of log( 0 /yr) is randomly assigned to each object within the range 7 to 8.
Before moving on, we checked that the infrared luminosity function (LF) of the simulated sample is consistent with current measurements.For a given redshift, we first computed the age of each object at that redshift and then used Eq. ( 23) to compute the corresponding dust luminosity,   .Because the latter represents the amount of energy re-emitted in the far-IR by dust, we took   to be the total infrared luminosity,  IR .We then binned the objects in  IR to generate a simulated infrared luminosity function to be compared with the data.
The results are illustrated in Fig. 6 for  = 2, 3, and 4. The simulated LF (red curve) shows a very good agreement with the C13 model (black curve) at moderate to bright infrared luminosities, i.e. above  IR ≳ 10 11.5  ⊙ .However, it lies significantly below the prediction from the C13 model at low luminosities.This is a consequence of the way C13 computed the infrared luminosity.Indeed, they calculated  IR from the SFR of the objects using the relation by Kennicutt (1998, K98 hereafter).However, that relation only applies during the early dust-obscured phase of the evolution of galaxies, when the far-IR emission is mainly associated with BCs where the dust is heated by young massive stars (R14).It does not hold for more evolved, lowluminosity, galaxies where the ISM contributes a significant fraction of the dust emission.This is illustrated in Fig. 7 where the SFR of our simulated galaxies is shown as a function of their infrared luminosity, computed using Eq. ( 8), and colour-coded according to the value of   .A high value of   corresponds to a dominant contribution of the ISM to the dust luminosity.Below log( IR / ⊙ ) ≃ 11, where   > 0.5, there is a clear deviation from the K98 relation, with the latter significantly underestimating the infrared luminosity for any given SFR.As a result, more objects with lower  IR are expected when using the K98 relation.This is the reason why our simulated LF curve lies below the C13 one at lower luminosities.
The shape of the LFs is highly dependent on the values of the parameters that determine   via Eq.( 8), especially of  0 ,   and  ISM  .Figure 8 shows how the IR LF at  = 2 is affected by changes in the values of the above parameters.The parameters are changed one at a time while keeping the others fixed at their prescribed values.It is evident that a deviation from the R14 values produces LFs that lie below the reference LF by C13 (black curve).
By increasing the dissipation timescale of the BC, stars are allowed to reside inside the BC for a longer time.One would expect that a prolonged dust absorption of UV/optical light from the stars will lead to a higher IR luminosity, hence to an increased infrared LF.But, in reality, that is not the case.The massive stars, which contribute the bulk of the IR luminosity, have significantly evolved by the time the BC dissipates and hence, have reduced luminosity.So, there are fewer UV/optical photons being emitted from them to be absorbed and re-emitted at IR wavelengths.Moreover, by increasing the value of  0 and by getting it closer to  AGN , the BCs are eventually destroyed by the AGN feedback.This effect also leads to a lower IR luminosity.
Similarly, by increasing  ISM  while keeping   constant (which is equivalent to decreasing  BC  ), a lower fraction of the starlight is re-radiated in the far-IR, leading to a decrease in the number of galaxies for any given  IR compared to the reference case.

Step 3: generating the far-IR/sub-mm/mm part of the SED
For each simulated object we compute the far-IR/sub-mm/mm part of the SED by implementing the formalism illustrated in Sec.2.2.2, and setting the values of the relevant parameters to those proposed by da Cunha et al. (2015).The indices  for warm and cold dust grains were fixed at 1.5 and 2 respectively.The temperature of the warm dust grains in the BCs,  , , is allowed to vary uniformly in the interval 30 K to 80 K, while that of the ISM,  ,  , is fixed at 45 K.The temperature of the cold grains in the ISM,  ,  , is taken from a uniform distribution between 20 K to 40 K.In the BCs,  PAH BC and  MIR BC are allowed to vary uniformly in the ranges 0 − 0.1 and 0.1 − 0.2, respectively.For each choice of both  PAH BC and  MIR BC ,  BC  follows from Eq. ( 17).For the ISM, we extract  ISM  from a uniform distribution between 0.5 and 1.The other parameters ) so that Eq. ( 19) is satisfied.With this parameter set, we generated differential number counts at 850 m, 870 m and 1100 m.We did that by first computing the objects' flux density at several redshifts from  = 1 to  = 8 in steps of 0.1.Each object was randomly assigned a redshift within the corresponding bin assuming a uniform distribution.We then computed the number of objects per unit interval of redshift and per unit interval of flux density at each redshift step.Finally, we integrated the redshift-dependent number counts to produce the differential number counts.The result is shown by the red curve in Fig. 9, where it is compared with the C13 number counts model curve (black line) of the proto-spheroids and with existing data (see also Table 1).
We notice that the simulated number counts significantly underestimate the observed number counts and the C13 curve at low flux densities, i.e. ≲ 6 mJy, while producing, at the same time, an excess of bright objects.This result indicates that the assumed values for the parameters that describe the far-IR/sub-mm/mm SED of the protospheroids need to be revised.In principle, the best-fit values of the SED parameters could be found using a minimum  2 approach on the number counts.However, that method, in this specific case, is very time consuming.Therefore, we opted instead, for a trial and error approach, which led to the following changes in the value of some of the SED parameters.We assigned a higher value to the dust emissivity index of the warm dust grains, raising it from  = 1.5 to  = 2, while we kept the one of the cold dust grains to its original value of  = 2.The temperatures  , and  ,  were allowed to vary uniformly within the intervals 30-45 K and 20-30 K, respectively.The value of  ,  was left unchanged.By using these updated values for the SED parameters we achieved a very good agreement with the observed number counts at the faint end, as illustrated by the Correlation between SFR and total IR luminosity,  IR , for a sample of simulated galaxies at  = 2.  IR is calculated using Eq. ( 8).Data points are colour-coded in accordance with the value of   .The black line shows the relation between SFR and total IR luminosity (integrated from 8 to 1000 m) from K98 with a 30% dispersion shown by the shaded blue area.Galaxies with a low value of   lie on the K98 relation.However, for galaxies where the ISM contributes significantly to the dust luminosity (i.e.those with high values of   ) the K98 relation leads to an underestimate of  IR .solid cyan curve in Fig. 9.However, we were still left with an excess of objects at bright flux densities.
Figure 10 shows that, as expected, the bright tail of the number counts is made up by the intrinsically brightest objects, i.e. those with total infrared luminosity  IR ≳ 10 13  ⊙ .Therefore, we assigned a separate set of SED parameter values to these very bright objects.Because the dust temperature increases with increasing infrared luminosity (Hwang et al. 2010), we let  , and  ,  to vary uniformly in the range 30-70 K and 20-40 K, respectively.At the same time, we also decreased the dust emissivity index of the warm dust grains to  = 1.8.With these modifications, we obtained the number counts shown by the red solid curve in Fig. 11.The agreement with the measured number counts is now satisfactory across the whole range of flux densities probed by the data.

Step 4: generating the final sample
With the adopted values of the SED parameters, which reproduce both the infrared luminosity function and number counts of protospheroidal galaxies, we generate the final catalogue of simulated objects used for forecast studies.We simulate proto-spheroidal galaxies over an area of 100 deg 2 , which is of the same order of the wide area surveys conducted with Herschel.∈  (0.9, 1.1).Top panel: varying  0 uniformly in logarithmic scale in the ranges 8-9 (red), 8.5-9.5 (orange) and 9-10 (green), respectively.Middle panel: varying   uniformly in the ranges 0-1 (orange), 1-2 (green), and 1-1.5 (red), respectively.Bottom panel: varying  ISM  uniformly in the ranges 0-8 (orange), 0-12 (green), and 2-10 (red), respectively.

The Herschel-ATLAS
The Herschel-ATLAS (H-ATLAS; Eales et al. 2010) is a survey of 600 deg 2 conducted in 5 photometric bands (100 m, 160 m, 250 m, 350 m and 500 m), using the Photoconductor Array Camera and Spectrometer (PACS) and the Spectral and Photometric Imaging Receiver (SPIRE) cameras.The H-ATLAS observed fields in the northern and southern hemispheres and on the celestial equator.The fields were chosen to minimise the dust emission from the Milky Way.The observed fields are as follows: (i) A field with an area of 150 deg 2 close to the North Galactic Pole (NGP).
(ii) Three equatorial fields each of approximate area of 56 deg 2 coinciding with the previously surveyed fields in the Galactic and Mass Assembly (GAMA) spectroscopic survey.
(iii) Two fields having a total area of 250 deg 2 , close to the South Galactic Pole (SGP).Because of the higher sensitivity to high-z sources of SPIRE compared to PACS, hereafter we will only focus on the SPIRE wavelengths.The 1 noise for source detection at 250 m, 350 m and 500 m, is 7.4 mJy, 9.4 mJy and 10.2 mJy, respectively.

Euclid wide area survey
The Euclid Wide Survey (EWS; Euclid Collaboration 2022b) will map 15000 sq.deg. of the sky from the Sun-Earth Lagrange point L2, using all four of its filters -visible (  ) and NIR (  ,   and   ).The survey will have a 5  depth of 26.2, 24.3, 24.5 and 24.4 mag for   ,   ,   , and   , respectively.

Rubin's Legacy Survey of Space and Time (LSST)
The Vera C. Rubin Obervatory (previously known as Large Synoptic Survey Telescope; Ivezić et al. 2019), which is under construction in Cerro Pachón in Northern Chile, is planned to conduct a survey named Legacy Survey of Space and Time (LSST) covering approximately 18 000 sq. deg. of the sky in the optical in 6 bands -, , , , , .The survey will have a 5  depth of 23.8, 24.5, 24.03, 23.41, 22.74 and 22.96 mag respectively.
A list of the Euclid and Rubin filters along with their central wavelengths and bandwidths is provided in Table 2.

Simulation set-up
A simulated source is said to be detected by Herschel if it has a flux density above the 5  limit of 37 mJy at 250 m, and to be detected in one of the Euclid bands if it has a flux density higher than the 3  limit at that band, i.e. 26.75, 24.85, 25.05, 24.95 mag for   ,   ,   and   , respectively.Moreover, we say that an Euclid source is detected by LSST if it has a flux density > 3  in at least 4 out of the 6 bands.
To speed up the calculation, we simulated objects with a minimum virialization halo mass of log( vir / ⊙ ) = 12.5, which corresponds to an infrared luminosity of log( IR / ⊙ ) = 12, as illustrated in  The simulated catalogue comprised 458 994 galaxies prior to the application of any selection criterion.Approximately 33% of them (i.e.151 720 objects) have  250 > 37 mJy.All of these Herschel detected galaxies are also detected in all the four Euclid bands above 3 .Among these sources, 140 524 have a 350 m flux density above 3  and 30 337 are detected above 3  at both 350 m and 500 m.Approximately 75% of the Euclid detected Herschel galaxies are also detected by Rubin/LSST above 3  in at least 4 out of the 6 bands.

Photometric redshifts
The first thing needed to extract the intrinsic properties of galaxies is an estimate of the redshift.When spectra are not available, an estimate of the redshift can be obtained from the photometry.Because we are dealing with Herschel sources detected by Euclid, we can use either the optical/near-IR photometry from Euclid or the sub-mm Herschel photometry, or a combination of both to estimate the redshift.

Redshifts from the optical/near-IR Euclid photometry
To estimate the photo- from the Euclid data we used the photometric redshift estimation code EAZY (Brammer et al. 2008), which is designed to work with UV/optical/near-IR data.EAZY fits the observed SED by linearly combining a number of template SEDs.It thus calculates the best-fit coefficients of the linear combination together with the photometric redshift.We used the EAZY version "v1.3" which comprises 9 (5+1+1+1+1) templates, shown in Figure 13.The original "v1.0" version of EAZY comprises 5 "principal component" templates obtained from Grazian et al. (2006) using the Blanton & Roweis (2007) algorithm.One dusty starburst template ( = 50 Myr,   = 2.75) was added to account for the extremely dusty galaxies.
The new adopted version has these 6 templates modified to include line emissions.It also includes one dust template, a template taken from Erb et al. (2010) and the evolved SSP by Maraston (2005) to account for massive old galaxies at  < 1. Dust absorption from the intergalactic medium (IGM), in accordance with Madau (1995), is also included by EAZY during its calculations.We selected a redshift range from  = 1 to 8 in steps of Δ = 0.01.

Redshifts from the sub-mm photometry
Another redshift estimate was derived from the Herschel/SPIRE photometry at 250 m, 350 m and 500 m.We fitted it with the empirical SED of Pearson et al. (2013), which is the sum of two modified blackbody spectra: where,   is the flux density at the restframe frequency ,  is the normalisation factor,   is the Planck function,  is the dust emissivity index,  1 is the hot dust temperature,  2 is the cold dust temperature and  is the ratio of the mass of the cold dust to that of the hot dust.We adopt  1 = 46.9K,  2 = 23.9K,  = 2 and  = 30.1 (Pearson et al. 2013;Negrello et al. 2017).The template SED was redshifted between 0 ≤  ≤ 8 in steps of 0.01 and a  2 minimisation was performed to estimate the photometric redshift.

Estimating the galaxies properties
Once the photo- is determined, we estimate the physical properties of the simulated galaxies using the SED fitting code CIGALE (Code Investigating GALaxy Emission; Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019).CIGALE is a code written in Python  which is used to efficiently model the multi-wavelength (UV to radio) spectrum of galaxies and to estimate the physical properties of the galaxies.It is based on the principle of "energy balance", where the stellar energy absorbed by dust in the UV-optical is entirely reemitted in the far-IR/sub-mm/mm.CIGALE creates model SEDs based on the modules and parameter values given as input by the user.Once an SFH is chosen, CIGALE couples it with the stellar population synthesis models to create the stellar spectrum of the galaxies.Several SFH functions are available depending on the physics of star formation and on the morphology of the galaxies.In selecting the modules of CIGALE, we took help from Traina et al. (2023), who studied 1620 sub-mm galaxies from 0.5 ≤  ≤ 6 by using the data from the largest available Atacama Large Millimeter Array (ALMA)1 survey in the COSMOS field (Scoville et al. 2007), known as  3 COSMOS (Liu et al. 2019).We chose a delayed SFH with an additional late burst of star formation.For modelling the stellar spectrum, we adopted the Chabrier IMF (Chabrier 2003) and chose the SSP models by BC03.The metallicity was set at the solar value ( ⊙ = 0.02).The dust attenuation was modelled using the attenuation law by Charlot & Fall (2000).
To model the dust emission, we adopted the templates from Draine et al. (2014).Finally, the AGN component was modelled using the templates from Fritz et al. (2006).The chosen modules and parameter values, which are summarized in Table 3, produced a total of 39 029 760 SED templates.CIGALE considered all these templates to perform SED fitting by implementing a  2 minimisation process (Noll et al. 2009).

RESULTS
This section presents the results obtained for the estimation of photometric redshift, stellar mass and SFR for the simulated catalogue for Euclid, using the SED fitting technique.

Derived photometric redshift
In Figure 14, we show the photometric redshift estimated with EASY versus the input redshift for the detected galaxies.There is a good level of agreement between the photometric redshift and the input redshifts for  input ≤ 3. The fraction of outliers, i.e. of objects with |Δ|/(1 +  input ) > 0.15 where |Δ| = | input −  EAZY phot |, as commonly defined in the literature (e.g., Laigle et al. 2016), is  outlier = 0.128.Most of the disastrous estimates of the photometric redshift are associated with galaxies at  input ≳ 3.5, which are wrongly assigned a redshift in the range  EAZY phot ∼ 1 − 1.5.Also, there are some galaxies with input redshift  input ∼ 1 − 2 which are estimated to be at much higher redshift, i.e.  EAZY ∼ 5 − 6. Figure 15 shows some of the best-fit SEDs produced by EAZY along with the corresponding estimated redshift.In order to achieve a good photometric redshift, EAZY mainly exploits the 4000 Å break.At low redshifts (i.e. input ≲ 3), such feature falls within the wavelength range sampled by Euclid, thus leading to photometric redshifts that agree with the input values.However, at higher redshifts, the 4000 Å break is missed by Euclid thus causing catastrophic redshift estimates.Moreover, due to the limited wavelength coverage, the 4000 Å break is wrongly identified as the Lyman− break at 912 Å, which also leads to a wrong redshift estimate.The addition of filters blueward of Euclid can help to identify the Lyman− break, thus leading to a more accurate redshift estimate.This is observed when the Euclid photometry is complemented with the LSST data.Figure 16 shows a plot of the input redshift vs the photometric redshift obtained by EAZY from the Euclid +LSST photometry.Out of the 113 897 Euclid sources that are also detected by LSST in at least four bands, 112 831 (approximately 99%) have |Δ|/(1 +  input ) ≤ 0.15.The outlier fraction  outlier is thus down to just 1%.
Figure 17 shows the contour plot of the derived photometric redshift of the same galaxies using the Herschel/SPIRE photometry alone.The fraction of outlier is  outlier = 0.37, more than twice the one obtained from the Euclid photometry.Moreover, there are 47 618 sources for which the Pearson template gives |Δ|/(1 +  input ) > 0.15 while EAZY produces |Δ|/(1 +  input ) ≤ 0.15.At the same time, there are 10 623 sources with |Δ|/(1 +  input ) ≤ 0.15 when using the Pearson template while EAZY produces |Δ|/(1 +  input ) > 0.15.Therefore, in general, the near-IR data seem to be more successful at estimating the redshift of the simulated galaxies compared to the Herschel photometry.However, it is worth pointing out that for the objects which are also detected above 3  at both 350 m and 500 m (only 20% of the sample), the sub-mm photometry provides a better estimate of the photometric redshift as illustrated in Figure 18.Therefore, the inaccuracy of the sub-mm photometry in estimating the redshift for the bulk of the objects can be attributed to the absence of robust flux density measurements at 350 m and 500 m.
In the following analysis, we will rely on the photometric redshift estimated by EAZY from the Euclid photometry.

Derived physical properties
In order to extract the physical properties of the simulated galaxies, such as the stellar mass,  ★ , and the star formation rate, SFR, we fitted the SED of the simulated Herschel galaxies detected with Euclid by using CIGALE and adopting the photometric redshift produced by EAZY.In the following sections, we discuss the results obtained using the Euclid photometry alone and the Euclid plus Herschel/SPIRE photometry combined.We also analyse the improvement of the above when including the LSST photometry.Overall, the retrieved stellar mass is in good agreement with the input values.The median value of the stellar mass is found to be log( ★ / ⊙ ) = 11.29 ± 0.181.The 1  dispersion in Δ log( ★ ) are 0.16 and 0.19 from the Euclid and Euclid+Herschel photometry, respectively.These dispersions are mainly a reflection of the uncertainties in the derived photometric redshifts.As expected, because the evolved stellar populations − which contribute most of the stellar mass in a galaxy − dominate the rest-frame optical/near-IR part of the spectrum, the stellar mass estimates produced by CIGALE do not improve when the SED fitting also includes the sub-mm Herschel data.However, the addition of the LSST data does improve the stellar mass estimate, by reducing the 1  dispersion in Δ log( ★ ) to 0.14.

Stellar Mass
We observe that, on average, CIGALE overestimates the stellar mass, the effect being particularly relevant for the higher masses.This offset can be attributed to the choice of the SFH and its parameterisation.Here, we choose a multi-component SFH, which, in general, has a tendency to yield a higher mass-to-light ratio, hence higher stellar mass estimates, as pointed out by Michałowski et al. (2012).The discrepancy in stellar mass estimates could also be explained by the difference in the value of the dust attenuation law in birth clouds used in CIGALE versus the one adopted here, i.e. -0.7 versus -1.3, respectively, and therefore to the extension of the extinction curve to the near-IR, where it directly impacts the stellar mass.

SFR
The comparison between the input values of the SFR and the estimated ones is shown in Figure 20.In this case, the SED fit has been performed using both the Euclid and the Herschel photometry.In fact, for dusty objects, the NIR wavelengths are not good tracers of the star formation (Pforr et al. 2012(Pforr et al. , 2013;;Collaboration et al. 2023) and the use of sub-mm data becomes crucial to reliably infer the SFR.The difference in logarithm between the input and the output values of the SFR has a dispersion of 0.26.As the LSST filters cover the UV/optical wavelengths, adding the LSST data to the NIR and FIR photometry can help to trace the SFR more accurately as can be seen in Figure 20 (right panels).In fact, when using the Eu-clid+Herschel+LSST photometry, the 1  dispersion in Δ log(  ★ ) drops to 0.18.However, we observe a systematic overestimation of the SFR by CIGALE when adding the LSST photometric data.As for the stellar mass estimates, the discrepancy in the SFR measurements can also be explained in terms of differences in the adopted values of the slope of the dust attenuation law between CIGALE and our formalism.Overall, estimating the SFR seems to be more challenging  than inferring the stellar mass.This may be due to the more sensitive dependence of the latter on the SFH compared to the former.In fact, the SFR is an instantaneous quantity while  ★ is an integrated quantity.To test this, we compared the values of the SFR averaged over a defined time scale.We introduce two quantities,  ★,10 and  ★,100 , which represent the SFR averaged over the 10 Myr and the 100 Myr prior to the time of observation, respectively, i.   dispersion to what is derived for the instantaneous SFR.Therefore, we can conclude that the use of an instantaneous SFR cannot explain the larger dispersion in the recovered values of SFR compared to what was observed for the stellar mass.However, the galaxies we are dealing with are dust-obscured; therefore, the accuracy in the estimate of the SFR is also dependent on the availability of data at far-IR/sub-mm/mm wavelengths.Only 20% of the simulated DSFGs detected by Herschel above 5  at 250 m are also detected above 3  at 350 m and at 500 m.Indeed, if we restrict the calculation of the SFR to that subsample we find a 1  dispersion of 0.23 and 0.14 for Euclid+Herschel photometry and Euclid+Herschel+LSST photometry, respectively, which clearly shows an improvement in the estimation of the SFR.The median SFR of the sources in our sample is log(  ★ / ⊙ yr −1 ) = 2.77 ± 0.162.

SFR-𝑀 ★ relation
Studies of large galaxy samples have shown that most of the starforming galaxies in the Universe follow an SFR- ★ correlation, called the galaxy main sequence (MS), while a fraction of galaxies lie above that correlation and are regarded as starburst (SB) galaxies (Speagle et al. 2014;Mancuso et al. 2016;Popesso et al. 2022).Figure 22 shows the SFR- ★ plot of our simulated galaxies for several bins of redshifts, along with the redshift-dependent galaxy-MS relation derived by Speagle et al. (2014).Our sample comprises both MS and starburst galaxies.The bulk of the galaxies, approximately 57%, with  ∼ 1 − 2, are above the MS relation and they are in their starburst stage, forming stars at an average rate of ∼ 300 − 3000  ⊙ per year.There are very few cases (11) where the estimated SFR is > 10 4  ⊙ yr −1 .This is due to the wrong redshift estimation by EAZY, leading to the wrong estimation of SFR by CIGALE.Around 15% of the galaxies at  ∼ 2−3 are above the MS relation and in their starburst stage.On the other hand, at high redshifts  ∼ 3.5 − 5.0, we find that the bulk of the sources are MS galaxies.

Dust Luminosity and Dust Mass
The comparison between the measured values of the dust luminosity and the input values is shown in Figure 23.We considered again the results obtained by using the Euclid+Herschel photometry (lefthand panels) and the Euclid+Herschel+LSST photometry (right-hand panels).The histogram of the differences between the input values and the measured values is shown in the bottom panels.We found that, when fitting the Euclid and the Herschel data together, we could estimate the dust luminosity ( dust ) accurately, with a 1  dispersion in Δ log( dust ) of about 0.22.On adding the LSST photometry to the fit, the constraint on  dust becomes tighter and the 1  dispersion in Δ log( dust ) reduces to 0.17.In fact, by tracing the UV/optical wavelengths, the LSST data help to better constrain the extinction of starlight by dust, thus leading to a more precise estimate of  dust .It is worth noticing, though, that in this case CIGALE systematically overestimates the dust luminosity by a factor of about 1.07.This is due to a higher standard value of the slope of the dust-attenuation law adopted by CIGALE.Overall, with a median dust luminosity of log( dust / ⊙ ) = 12.8 ± 0.12, these galaxies sample the bright tail of the population of ultra-luminous infrared galaxies (ULIRGs).
In Figure 24 we show the results for the dust mass estimates.The dust mass has a median value of log( dust / ⊙ ) = 8.9 ± 0.10.This is consistent with the values of  dust of other ULIRGs reported by Dudzevičiūtė et al. (2020) and Pantoni et al. (2021), and shows that these DSFGs have higher dust content than local ( < 1) star forming galaxies selected with Herschel (e.g.: R14).The dust mass is affected by larger uncertainties compared to the dust luminosity and the use of the LSST photometry does not help to reduce that uncertainty.We find, indeed, that the 1  dispersion in Δ log( dust ) is 0.35 and 0.32 for the Euclid+Herschel photometry and the Euclid+Herschel+LSST photometry, respectively.We also observe that CIGALE systematically overestimates  dust by a factor of 1.2.Interestingly, Liao et al. (2024) fit the SED of a sample of 18 bright ( 870 = 12.4 − 19.4 mJy)  2021) when comparing dust mass estimates by CIGALE and MAGPHYS for a sample 61 sources selected from ALMA-identified 870 m-selected DS-FGs from AS2COSMOS, AS2UDS (Stach et al. 2019), and ALESS (Hodge et al. 2013) surveys.These discrepancies in  dust can be mainly explained by the difference in the standard values of  0 adopted by the two SED fitting codes (Liao et al. 2024).

SUMMARY AND CONCLUSIONS
We have presented predictions for how well the recently launched space telescope Euclid can study the  ≳ 1.5 DSFGs detected by the Herschel-ATLAS survey.To this end, we have built a simulated catalogue of galaxies exploiting the C13 model for the evolution of proto-spheroidal galaxies.We combined the output of model equations with the SED formalism of D08, which is based on the principle of energy balance in the UV-optical and IR.
By setting the values of the main parameters of the model to those derived by R14 from the study of a sample of DSFGs we obtained a good agreement between the predicted infrared luminosity function of DSFGs and the measured one in the redshift range  = 2 − 4.However, compared with the prescriptions of R14, we fine-tuned the values of both the dust emissivity index and the temperature of the warm dust in the BCs and of the cold dust in the ISM to reproduce the observed number counts of at 850 m, 870 m and 1.1 mm.
After these preliminary checks, we simulated a catalogue of 458 994 galaxies with  IR ≳ 10 12  ⊙ within a survey area of 100 deg 2 , which is of the same order of magnitude of the total area surveyed by H-ATLAS, and applied a 5  detection condition at 250 m for Herschel and a 3  detection condition for all the 4 Euclid bands.All the 151 720 sources detected by Herschel at 250 m are also detected in all the Euclid bands.We used EAZY on the Euclid photometry to estimate the photometric redshifts of the above sources.We then exploited CIGALE to extract the physical properties of the galaxies from the SED fitting, assuming the previously derived photometric redshift.
Our main findings are as follows (i) For 87 % of the sample, the Euclid photometry alone provides a good redshift estimate up to  ∼ 3, with a discrepancy ≤ 15% from the true value of 1 + .For input redshifts  ≳ 3, we observe some catastrophic outliers, where the estimated redshift is significantly lower than the input one, i.e.  phot ∼ 1 − 1.5.This discrepancy is attributed to the fact that the 4000 Å break, which is used as the main feature to constrain the redshift, is not sampled by Euclid above  ∼ 3. Also, there are some galaxies at  ∼ 1 − 2 which are estimated to be at  ∼ 4 − 6, due to the Balmer break (4000 Å) being mistaken as the Lyman break (912 Å).Adding the shorter wavelength data from LSST does help to break some degeneracies.In fact, the outlier fraction reduces to just 1% for the 75% of the Euclid detected galaxies that are also detected above 3  in at least 4 LSST bands.
We also investigated the use of the sub-mm Herschel photometry to derive the photometric redshift by fitting those data with the empirical template of Pearson et al. (2013).However, we found a higher fraction of outliers compared to what was obtained with the Euclid photometry, i.e.  outlier = 0.35 compared to  outlier = 0.13.This is mainly due to the lack of sufficient sub-mm data for most of the sources in the sample.In fact, only 20% of the Herschel de- tected galaxies (i.e.those above 5  at 250 m) are also detected at both 350 m and 500 m above 3 .For this sub-sample the outlier fraction is reduced to 13%.In summary, for the vast majority of the simulated DSFGs, we cannot gain any significant improvement on the photo-z from the sub-mm measurements and that is why we have adopted the photometric redshift estimates based on the Euclid data alone.
(ii) The distribution of the stellar masses, as recovered from the Euclid and Euclid+Herschel photometry, has a dispersion of 0.16 dex and 0.19 dex respectively, around the true value.As expected, adding the Herschel photometry to the SED fit does not significantly improve the results, as the stellar mass is sensitive to the rest-frame optical/near-IR portion of the SED.The median value of the stellar mass is found to be log( ★ / ⊙ ) = 11.29 ± 0.181, while the range of values of log( ★ / ⊙ ) is 10.5 − 12.9.
(iv) From the SFR− ⊙ relation, we found that our simulated catalogue comprised both MS galaxies and starburst galaxies.Approximately 55% of the galaxies in our sample lie above the MS relation and are starbursting, with a median SFR of log( / ⊙ yr −1 ) = 2.87.
(vi) The photometric redshift can be improved by complementing the Euclid data with the LSST observations.From the combined Euclid+LSST photometry, EAZY provides an accurate redshift, i.e.Δ/(1 + ) ≤ 0.15, for 99% of the 87% Euclid objects that are also detected in at least 4 of the 6 LSST bands.For this subsample of galaxies that have both Euclid and LSST photometry, the dispersion on the recovered stellar mass and star formation rate reduces to 0.14 dex and 0.18 dex respectively.By covering the UV/optical wavelengths using LSST filters, we could constrain the extinction of starlight by dust in the UV, which reduced the dispersion in the estimation of  dust .However, the dispersion in dust mass shows no significant improvement with the inclusion of LSST photometry.

Figure 1 .
Figure1.Evolution of the properties of the stellar and AGN components of proto-spheroidal galaxies as a function of galactic age for a halo of mass  vir = 10 11.3  ⊙ (left) and a halo of mass  vir = 10 12.3  ⊙ (right).Both have virialisation redshift  vir = 3.In each plot, the left y-axis refers to the mass evolution of the stellar component: infalling gas (solid red line), cold gas (dot-dashed black line), and stellar mass (blue dashed line); while the right y-axis refers the mass evolution of AGN component: gas falling in the reservoir (solid green line) and super-massive black hole (SMBH) mass (dotted magenta line).We observe that the stellar mass increases as more and more cold gas is condensed into stars, thus leading to a corresponding steady decline in the mass of the gas components.At the same time, the central BH mass increases as a fraction of the cold gas loses its angular momentum via the stellar radiation drag and settles down into a reservoir around the central SMBH.The BH accretes matter from the reservoir via viscous dissipation, which powers the nuclear activity.With the evolution of the galactic age, the AGN feedback sets in.Its effect can be observed in the most massive dark matter (DM) halos where the star formation activity and, consequently, the stellar mass growth, are stopped within less than 1 Gyr from the galaxy's birth.

Figure 2 .
Figure 2. SFR as a function of the galaxy age in an object forming in a halo of mass  vir = 10 12.32  ⊙ virialised at redshift  vir = 2.The inset plot shows the evolution of the stellar mass with galaxy age.The black dashed line marks the time ( AGN ) when AGN feedback quenches SF in this galaxy and the SFR drops to zero.For this galaxy,  AGN = 618.0Myr.As there is no more SF going on, the stellar mass remains constant afterwards as illustrated by the stellar mass vs. time plot.

Figure 3 .
Figure 3. Example of the SED of a proto-spheroidal galaxy (red curve) with halo mass of  vir = 10 11.60  ⊙ and formation redshift  vir = 5.9.The SED accounts for the effect of the suppression of the stellar emission in the UV/optical/near-IR due to dust absorption and for the reprocessing at longer wavelengths of the starlight that has been absorbed by dust.The galaxy is observed at  = 5.4 at an age of 106.2 Myr when the SFR is ∼ 99  ⊙ yr −1 .The evolution of the SFR with time is shown in the inset plot where the dashed black line marks the age of the galaxy at the time of observation.The blue line shows the unattenuated stellar spectrum.The IR emission from the BCs is shown in orange while that from the ISM is shown in green.

Figure 5 .
Figure 5. Example of the SED (stellar + AGN components) of a proto-spheroidal galaxy (black curve) with halo mass of  vir = 10 11.65  ⊙ and formation redshift  vir = 5.8.The galaxy is observed at  = 5.1 at an age of 175.4 Myr when the SFR is ∼ 86  ⊙ yr −1 .The blue curve is the unattenuated stellar spectrum, while the red curve shows the starlight attenuated by dust in the UV/optical/near-IR and reprocessed by the same dust at longer wavelengths.The IR emission from the BCs is shown in orange while that from the ISM is shown in green.The unattenuated SED of the AGN component is shown in brown while the attenuated AGN SED is shown in cyan.The grey curve depicts the IR emission from the ISM due to the absorption of AGN luminosity by ISM dust.The purple curve shows the total AGN SED.The AGN SED is typical of a Type 2 AGN ( = 90 • ) with density profile ( ) =  −0.5 ,  eq (9.7) = 10,  max / min = 60 and Θ = 100 • .

Figure 7 .
Figure7.Correlation between SFR and total IR luminosity,  IR , for a sample of simulated galaxies at  = 2.  IR is calculated using Eq.(8).Data points are colour-coded in accordance with the value of   .The black line shows the relation between SFR and total IR luminosity (integrated from 8 to 1000 m) from K98 with a 30% dispersion shown by the shaded blue area.Galaxies with a low value of   lie on the K98 relation.However, for galaxies where the ISM contributes significantly to the dust luminosity (i.e.those with high values of   ) the K98 relation leads to an underestimate of  IR .

Figure 8 .
Figure 8.Effect of changing the values of  0 , τ  and τISM  on the shape of the total IR LF at  = 2.In all panels the black curve represents the theoretical curve obtained by C13 and the blue curve represents the simulated LF obtained by using the optimal values of the parameters: log  0 ∈  (7, 8),   ∈  (4.5, 5.7) and  ISM

Figure 9 .
Figure 9.Comparison between the simulated differential number counts (red curve) with the C13 model (dashed black dots) and the observed data at 850 m, 870 m and 1100 m.The cyan curve shows the number counts obtained by modifying the parameter values by trial and errpr.The new curve fits the faint end but produces an excess of objects at the brightest luminosities.The data points are taken from Simpson et al. (2019), Shim et al. (2020),Garratt et al. (2023), Stach et al. (2018), Simpson et al. (2020), Scott et al. (2010), and Scott et al. (2012).

Figure 10 .
Figure 10.Simulated differential number counts at 850 m, with contributions from intervals of total infrared luminosity.The dot-dashed black curve represents the model number counts by C13.

Fig. 12 .
Fig. 12.In fact, simulated sources with  250 = 37 mJy have a minimum halo mass of 10 12.5  ⊙ and a corresponding minimum IR luminosity of 10 12  ⊙ .The simulated catalogue comprised 458 994 galaxies prior to the application of any selection criterion.Approximately 33% of them (i.e.151 720 objects) have  250 > 37 mJy.All of these Herschel detected galaxies are also detected in all the four Euclid bands above 3 .Among these sources, 140 524 have a 350 m flux density above 3  and 30 337 are detected above 3  at both 350 m and 500 m.Approximately 75% of the Euclid detected Herschel galaxies are also detected by Rubin/LSST above 3  in at least 4 out of the 6 bands.

Figure 12 .
Figure 12.Distribution of infrared luminosity,  IR , and virial mass,  vir , of the simulated proto-spheroidal galaxies as a function of their 250 m flux density,  250 .The vertical black line denotes the Herschel-ATLAS 5  flux density limit of 37 mJy at 250 m.

Figure 13 .
Figure 13.The set of EAZY SED templates used in this work.

Figure 14 .
Figure 14.Derived photometric redshift using EAZY vs the input redshift of the galaxies detected by Herschel which are also detected by all the four bands of Euclid at ≥ 3 .The solid red line denotes  input =  EAZY phot while the dashed red lines define the region where |Δ | ≤ 0.15(1 +  input ).

Figure 19
Figure19shows the comparison between the derived stellar masses,  CIGALE ★

Figure 15 .
Figure 15.Examples of simulated near-IR SEDs fitted with EAZY, illustrating the cases of both good and catastrophic redshift estimates.The blue curve denotes the simulated SED while the orange curve is the best-fit SED produced by EAZY.The green points denote the "observed" flux density values at the four Euclid bands.

Figure 16 .
Figure 16.Derived photometric redshift using EAZY vs the input redshift of the galaxies detected by Herschel which are also detected by all the four bands of Euclid and any four out of six bands of LSST at ≥ 3 .The solid red line denotes  input =  EAZY phot while the dashed red lines define the region where |Δ | ≤ 0.15(1 +  input ).
Figure 21 shows the distribution of the input values versus the measured values of  ★,10 and  ★,100 along with the difference in logarithm between the true and the estimated values of the SFR from the Euclid+Herschel and the Euclid+Herschel+LSST photometry.Both Δ log(  ★,10 ) and log(Δ  ★,100 ) have a 1  dispersion of 0.26 when estimated from the Euclid+Herschel photometry.Upon adding the LSST data to the former, the 1  dispersion reduces to 0.18.Even though  ★,10 and  ★,100 are integrated properties, we get similar

Figure 17 .
Figure 17.Contour plot of the derived photometric redshifts using Pearson templates on Herschel photometry vs the input redshift of the Euclid detected Herschel galaxies.The red solid line denotes  input =  pearson phot while the red dashed lines bound the region where |Δ | ≤ 0.15(1 +  input ).

Figure 18 .
Figure 18.Contour plot of the derived photometric redshifts using Pearson templates on Herschel photometry vs the input redshift of the Herschel/Euclid galaxies detected above 3  at both 350 m and 500 m.The solid red line denotes  input =  pearson phot while the dashed red lines bound the region where |Δ | ≤ 0.15(1 +  input ).

Figure 19 .
Figure 19.Comparison between the input and the estimated values of the stellar masses obtained with CIGALE from Euclid photometry alone (left panels), Euclid + Herschel photometry (middle panels) and Euclid + Herschel+ LSST photometry (right panels).The top panels show the contour of the estimated vs the true stellar masses for both cases.The black straight line marks the 1:1 relation.The black dashed line marks the boundaries of the 1  dispersion around the mean.The bottom panels show the distribution of differences between the logarithms of the input stellar mass and of the estimated stellar mass.The black solid line denotes the mean value while the dashed lines mark the 1  dispersion around the mean.

Figure 20 .
Figure 20.Same as Figure 19, but for the SFR, with the SED fit being performed on Euclid+Herschel/SPIRE photometry (left) and Euclid+Herschel/SPIRE+LSST photometry (right).

Figure 21 .
Figure 21.Same as Figure 19, but for the SFR averaged over 10 Myrs (first two top and bottom figures) and 100 Myrs (last two top and bottom figures) respectively, with the SED fit being performed on Euclid+Herschel/SPIRE photometry and Euclid+Herschel/SPIRE+LSST photometry.

Figure 22 .
Figure 22.Star formation versus stellar mass for the simulated sample of galaxies detected by Herschel and Euclid.The top left panel shows the correlation for the entire catalogue of galaxies.The points are colour-coded according to their photometric redshift,  EAZY phot .The remaining plots show the SFR versus  ★ for different redshift bins as follows: 1 ≤  ≤ 2 (top right), 2 <  ≤ 3 (middle left), 3 <  ≤ 4 (middle right), 4 <  ≤ 5 (bottom left), and 5 ≤  ≤ 6 (bottom right).The black line in each of the plots represents the galaxy-MS relation from calculated at the central redshift of the interval.The blue-filled region denotes the 0.3 dex scatter around that relation.The dashed black line shows the MS relation scaled up by a factor of four, the conventional boundary between MS and starburst galaxies.

Table 1 .
Information about the measured differential number counts at 850 m, 870 m and 1100 m shown in Figs 9-11.

Table 2 .
Euclid filters along with their central wavelengths and bandwidths.

Table 3 .
Parameter values given as input to CIGALE for SED fitting (Fritz et al. 2006)tz et al. 2006)