© The Authors 2025

1. Introduction

Understanding the nature of dark matter (DM) is one of the priority targets within the communities of cosmology, astroparticle physics, and high-energy physics. Over the past decade, the Large Hadron Collider (LHC) results and the absence of direct or indirect DM detection have shown that the situation concerning the nature of DM is wide open. Weakly interacting massive particles (WIMPs) are only one candidate among many possibilities (Bertone et al. 2005; Feng 2010). Particle-like DM could have a large range of plausible masses, lifetimes, annihilation cross-sections, and scattering cross-sections.

The standard cosmological model makes the working assumption of a purely stable, decoupled, and cold dark matter (CDM) species, which can be modelled as dust in simulations of the evolution of the Universe since very early times – well before photon decoupling. In the CDM limit, the only measurable parameter related to DM is its relic non-relativistic density today, ρ_cdm, which can be expressed in terms of a dimensionless density parameter, ω_cdm := Ω_cdmh², where Ω_cdm is the fractional density of CDM (relative to the critical density) and h := H₀/(100 kms⁻¹ Mpc⁻¹) is the reduced Hubble parameter. However, in non-minimal scenarios, DM could have several other parameters of possible relevance for fitting cosmological observations, such as: a non-negligible velocity dispersion (Bond & Szalay 1983; Bode et al. 2001), a lifetime not considerably larger than the age of the Universe (e.g. Ichiki et al. 2004; Audren et al. 2014), and cross-sections describing either its self-interaction (Spergel & Steinhardt 2000; Feng et al. 2009) or its feeble interaction with other species (Boehm et al. 2001; Cyr-Racine et al. 2016).

From the point of view of a particle physics model-builder, non-minimal DM models are easy to motivate; typically, they do not require more complicated or more fine-tuned ingredients that the particle physics models leading to plain CDM. The logic pursued successfully by high-energy physicists for almost a century consists of postulating additional symmetries (rather than adding individual particles) in order to explain unaccounted experimental results. The current standard model of particle physics is known to be incomplete (Workman et al. 2022) and the assumption of new symmetries usually comes together with a rich dark sector; that is, several new particles with new interactions, with potentially more than one population surviving until today and

contributing to DM or dark radiation. From this point of view, having just one decoupled, stable, and cold relic in our Universe does not sound much more natural than being surrounded by one or more dark species with potentially non-trivial properties. High-energy physicists often suggest that, given the richness of the visible sector, there is no obvious reason for the dark sector to reduce to a single CDM relic particle.

The astrophysics community is sometimes reluctant to investigate the possible consequences of non-minimal particle-physics assumptions in cosmology as long as the minimal ΛCDM model has not been ruled out. The situation is, however, evolving given the accumulation of tensions or unresolved questions in cosmological observations (like the small-scale CDM crisis, Hubble tension, or S₈ tension; see for instance Verde et al. 2019, Abdalla et al. 2022). In this context, it sounds at least reasonable to investigate the possibility of detecting some effects induced by non-minimal DM models. Of course, it is still possible that future observations only provide bounds on these models and leave us with plain CDM as a preferred case. Even in this case, it would be extremely interesting for particle physics model-builders to have such bounds, since constraints from accelerators or astroparticle experiments usually probe a different regime or different model assumptions than cosmological data.

Non-minimal DM properties may affect the growth of structures in the Universe in different ways, at different times, and on different scales. Thus, they can leave several types of signatures in the two-point correlation function of matter fluctuations in Fourier space, called the matter power spectrum. This spectrum can be reconstructed from several types of cosmological observations at different redshifts. The modified growth of structure induced by non-minimal DM models could also affect other statistical probes of structure formation (higher-order correlation functions, halo mass function, peak and void statistics), but in this work we only consider its impact on the matter power spectrum.

Euclid (Euclid Collaboration: Mellier et al. 2025) is a medium-class mission of the European Space Agency, which will map the local Universe to improve our understanding of the expansion history and of the growth of structures. The satellite will observe roughly 15 000 deg² of the sky through two instruments, a visible imager (VIS, Euclid Collaboration: Cropper et al 2025) and a Near-Infrared Spectrometer and Photometer (NISP, Euclid Collaboration: Jahnke et al. 2025), delivering the images of more than one billion galaxies and the spectra of tens of millions of galaxies out to redshift of about 2. The combination of spectroscopy and photometry will allow us to reconstruct the matter power spectrum up to an accuracy of 1%.

Since the matter power spectrum could be strongly affected by the nature of DM, Euclid is a perfect tool for testing non-minimal DM properties. It may either confirm the standard CDM paradigm or discover some new DM features. The goal of this work is precisely to estimate the sensitivity of Euclid to different DM parameters beyond its mere relic density. Given the wide range of possible alternatives to standard CDM, we cannot explore all possibilities. We instead concentrate on four examples of non-minimal scenarios that are still compatible with current data and that could be either constrained or detected by Euclid. Our choice of models is dictated by simple considerations. First, we should select some representative cases. Since in non-minimal models, DM particles are expected to either free-stream (with some velocity dispersion) and/or decay (with some rate) and/or scatter (with some cross-sections), we go through examples in each of these three categories. A well-motivated example of DM with a velocity dispersion is warm DM. Some simple examples of decaying DM consist of particles with a constant decay rate, decaying into either relativistic or non-relativistic daughter particles; and a representative case of scattering DM is that of DM interacting with dark radiation. Second, we are interested in models such that galaxy redshift surveys could provide stronger bounds than other observables, and in particular, than cosmic microwave background (CMB) and/or Lyman-α forest data. For reasons detailed in the next sections, this would not be the case for pure warm DM or pure decaying DM. Thus, going to the next level of complexity, we assume a mixture of either cold and warm DM, or of stable and unstable particles. In conclusion, we focus on four interesting and representative models: a mixture of cold and warm DM, a mixture of stable and unstable particles decaying into either relativistic or non-relativistic particles, and DM interacting with dark relativistic relics.

Euclid will deliver several types of observations. Among these, the weak lensing (WL) survey and the galaxy clustering (GC) photometric survey will be ideal to constrain DM properties, since they will both provide a measurement of the matter power spectrum down to small scales and up to high redshift. These two surveys will return maps in tomographic bins that can be analysed all together (taking into account cross-correlations between WL and galaxy density maps). In addition to this joint data set, called the photometric probe, Euclid will provide other observations. The Euclid spectroscopic galaxy redshift survey will play an essential role in constraining several cosmological models and parameters. Cluster number counts will also convey very useful information. However, these surveys will not provide information on such small scales as WL, and their implementation in sensitivity forecasts relies on a different methodology than for the photometric probe. In particular, they require a different approach to model non-linear effects for each non-minimal DM model. Thus, for simplicity, we choose to concentrate only on the Euclid photometric probe in this work.

In Sect. 2 of this work, we review the four DM models that we investigate, with a brief discussion of their foundations, their free parameters, and their effects on the linear matter power spectrum. In Sect. 3, we explain how to model the effect of these scenarios at the level of the non-linear power spectrum, using emulators trained on dedicated N-body simulations. In Sect. 4, we summarise the assumptions and numerical pipelines used in our parameter sensitivity forecasts for the Euclid photometric probe. We present our results for each model in Sect. 5 and provide final conclusions in Sect. 7.

2. Non-minimal particle dark matter models

Many particle DM properties can be tested with cosmology (Gluscevic et al. 2019). As was mentioned in the introduction, we only focus here on four particular models. On the one hand, these models constitute representative samples of the three most plausible properties of non-minimal particle DM: a non-negligible velocity dispersion, some decay rate, or a non-negligible scattering rate. On the other hand, within their respective category, they account for the simplest scenarios that can be constrained better by WL and galaxy surveys than CMB and Lyman-α data.

2.1. Cold plus warm dark matter

In this model, a fraction, f_wdm, of the total DM fractional density, Ω_dm, is assumed to be warm, so that Ω_dm = Ω_cdm + Ω_wdm = (1 − f_wdm) Ω_dm + f_wdm Ω_dm. The warm dark matter (WDM) component possesses a thermal (root mean square) velocity v_rms that depends on the temperature-to-mass ratio T_wdm/m_wdm. WDM would revert to CDM in the limit v_rms → 0, or equivalently m_wdm → ∞.

This mixed cold plus warm dark matter (CWDM) model has been studied previously, for instance, in Boyarsky et al. (2009a), Schneider (2015), Murgia et al. (2017), Murgia et al. (2018), Parimbelli et al. (2021), or Hooper et al. (2022). It may account either for cosmologies with two distinct DM components, or also, effectively, for cosmologies with a single DM component with a non-thermal distribution, such as resonantly produced sterile neutrinos (Boyarsky et al. 2009b). This model has been often invoked as a possible solution to the small-scale CDM crisis (Anderhalden et al. 2013; Maccio et al. 2013). Current best constraints come from Lyman-α forest surveys (Hooper et al. 2022), Milky Way satellites (Diamanti et al. 2017), and WL surveys (Hervas-Peters et al. 2024; see Sect. 5.1).

The thermal velocity of WDM defines its maximum free-streaming scale, reached at the time of its non-relativistic transition during radiation domination. On larger wavelengths, cosmological fluctuations have the same evolution as in a model in which all the DM would be cold. On smaller scales, the perturbations of the WDM component are negligible and the growth of CDM density fluctuations is suppressed. Thus, at the level of linear perturbations, the overall effect of WDM is to induce a step-like suppression in the matter power spectrum. The amplitude of the step is controlled by f_wdm¹. The shape of the step is universal for all models in which the WDM phase-space distribution has a thermal shape up to a rescaling factor. This covers two well-known limits: on the one hand, thermal WDM, for which the phase-space distribution is thermal (with no rescaling factor) but the WDM temperature, T_wdm, is reduced compared to the active neutrino temperature, due to its earlier decoupling time; and the Dodelson–Widrow (DW) model (Dodelson & Widrow 1994; Colombi et al. 1996), for which the phase-space distribution is identical to that of active neutrinos (with T_wdm = T_ν) up to a rescaling factor χ ≪ 1 accounting for the efficiency of active-sterile neutrino oscillations in the early Universe with a small mixing angle, χ ∼ sin²θ. For this broad category of models, the location of the step-like suppression is controlled by the thermal velocity; that is, by the temperature-to-mass ratio, T_wdm/m_wdm.

It is convenient to parameterise the location of the step in terms of the rescaled mass

$$\begin{array}{c}\hfill x:=m_{\mathrm{wdm}}\frac{T_{\nu }}{T_{\mathrm{wdm}}},\end{array}$$(1)

where T_ν is the current value of the active neutrino temperature computed in the instantaneous decoupling limit; that is, such that T_ν/T_γ = (4/11)^1/3. For the class of models described above, the effect of WDM is entirely described by the two parameters (f_wdm, x), independently of the chosen model (thermal WDM or DW). In the DW case, x coincides with $m_{\mathrm{wdm}}^{\mathrm{DW}}$. In the thermal case, one has

$$\begin{array}{c}\hfill m_{\mathrm{wdm}}^{\mathrm{thermal}}={\left(94.1{\mathrm{\Omega }}_{\mathrm{wdm}}{h}^2\right)}^{1/4}{\left(\frac{x}{1\mathrm{eV}}\right)}^{3/4}\mathrm{eV},\end{array}$$(2)

where we used the fact that for Fermi–Dirac thermal relics X with a temperature T_X = T_ν one gets m_X = 94.1 Ω_Xh² eV².

In this model, the evolution of linear cosmological perturbation can be computed with the public version of CLASS. Then, to account for the thermal WDM case, we pass to the code the parameters Ω_wdmh² = f_wdm Ω_dmh², $m_{\mathrm{wdm}}^{\mathrm{thermal}}$, and finally $T_{\mathrm{wdm}}/T_{\gamma }={(4/11)}^{1/3}(m_{\mathrm{wdm}}^{\mathrm{thermal}}/x)$ with x inferred from Eq. (2)³. In principle one could use a different set of parameters for the equivalent DW model and find the exact same linear power spectra (Lesgourgues 2011; Blas et al. 2011; Lesgourgues & Tram 2011). Figure 1 shows the power spectrum at redshift zero for several CWDM models rescaled by that of a pure ΛCDM model, for various values of (f_wdm, x) but fixed values of the usual ΛCDM parameters (Ω_m, Ω_b, h, A_s, n_s) accounting respectively for the fractional density of total non-relativistic matter (baryonic plus dark), the fractional density of baryonic matter, the reduced Hubble parameter, and the amplitude and spectral index of the primordial spectrum of scalar (curvature) perturbations. In the left panel, we vary x (or equivalently $m_{\mathrm{wdm}}^{\mathrm{DW}}$) with fixed f_wdm to check that x only controls the location of the step. In the right panel we do the opposite to show that f_wdm controls its amplitude.

Fig. 1. Ratio of the linear (solid lines) and non-linear (dashed lines) power spectra of several CWDM models to that of a pure ΛCDM model with the same cosmological parameters, parameterised by the fraction fwdm and the rescaled mass x. The other parameters (Ωdm, Ωb, h, As, ns) are kept fixed. All spectra are computed today (z = 0). These plots cover all the cases in which WDM has a Fermi–Dirac distribution possibly rescaled by a factor χ, including the limits of the thermal WDM (χ = 1) and Dodelson–Widrow (Twdm = Tν) models. In the latter case x coincides with the physical mass. The non-linear spectra are predicted by the emulator introduced in Sect. 3.1 and plotted up to the maximum wavenumber at which this emulator is trusted.

In Sect. 3.1, we show how to compute the impact of CWDM on the non-linear matter spectrum. In Sect. 5.1, we perform Euclid forecasts on the parameter of the CWDM model. For that purpose, we use a Bayesian MCMC approach to fit the CWDM model to mock Euclid data, assuming a logarithmic prior on f_wdm ∈ [2 × 10⁻³,  1] and a flat prior on $m_{\mathrm{wdm}}^{\mathrm{thermal}}\in [10\mathrm{eV},1\mathrm{keV}]$.

Such a logarithmic prior on f_wdm allows us to assess precisely the constraining power of Euclid even when f_wdm is very small (e.g. in the range from 10⁻³ to 10⁻¹). This limit is the most interesting in the case of the Euclid probes since, in this case, the data may be compatible with a relatively small WDM mass, and thus a small step located on relatively large wavelengths, in the range probed by WL and GC surveys in the linear and mildly non-linear regime. Large values of f_wdm (e.g. in the range from 0.1 to 1) imply a strong suppression of the power spectrum that is already excluded by Lyman-α forest data unless the mass is really large – a limit in which, from the point of view of Euclid data, CWDM would be indistinguishable from pure CDM.

2.2. Dark matter with one-body decay

If DM particles are unstable, they may decay in different ways into lighter particles. Cosmological observables are not sensitive to all details concerning the nature of the decay products, but they depend on simple considerations like the fact that these decay products could be relativistic or non-relativistic. In the simplest scenario, all decay products are assumed to be ultra-relativistic and can be considered as a single species, dubbed dark radiation (DR). This simple model of decaying dark matter (DDM) is often called one-body decaying DM and abbreviated as 1b-DDM.

In this section, we assume that DM is made up of two cold species: a fraction 1 − f_ddm of stable dark matter (CDM) and a fraction f_ddm of 1b-DDM decaying into DR. For simplicity, we assume a constant decay rate, Γ_ddm = 1/τ_ddm, where τ_ddm is the lifetime of the decaying species. The current value of the fractional dark radiation density, Ω_dr, is not an independent parameter of the model: it can be computed consistently for each value of f_ddm and Γ_ddm.

This model has been studied previously; for instance, in Ichiki et al. (2004), Audren et al. (2014), Berezhiani et al. (2015), Chudaykin et al. (2016), Chudaykin et al. (2018), Oldengott et al. (2016), Poulin et al. (2016), Pandey et al. (2020), Xiao et al. (2020), Nygaard et al. (2021), Schöneberg et al. (2022), Simon et al. (2022), Holm et al. (2023), or Bucko et al. (2023). It has been often invoked as a possible solution to the Hubble and/or S₈ tension. The best constraints at the moment come from CMB plus baryon acoustic oscillation (BAO) data (Nygaard et al. 2021), galaxy surveys (Simon et al. 2022), and WL surveys (Bucko et al. 2023; see Sect. 5.2).

In this model, the evolution of linear cosmological perturbations can be computed with the public version of CLASS⁴. The code accepts two possible definitions of the decaying DM fraction: one can either pass the value of f_ddm today, taking the effect of decay into account, or the value $f_{\mathrm{ddm}}^{\mathrm{ini}}$ evaluated at some initial time τⁱⁿⁱ ≪ τ_ddm, before any significant decay has occurred,

$$\begin{array}{c}\hfill f_{\mathrm{ddm}}^{\mathrm{ini}}=\frac{{\rho }_{\mathrm{ddm}}^{\mathrm{ini}}}{{\rho }_{\mathrm{dm}}^{\mathrm{ini}}}=\frac{{\rho }_{\mathrm{ddm}}^{\mathrm{ini}}}{{\rho }_{\mathrm{cdm}}^{\mathrm{ini}}+{\rho }_{\mathrm{ddm}}^{\mathrm{ini}}}.\end{array}$$(3)

Here we choose to report results on $f_{\mathrm{ddm}}^{\mathrm{ini}}$, for the purpose of easier comparison with previously published bounds. Some related parameters are ${\mathrm{\Omega }}_{\mathrm{ddm}}^{\mathrm{ini}}$ (respectively ${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$), the fractional density that DDM (respectively total DM) would have today if DDM did not decay. The free parameters of the 1b-DDM model are then (Γ_ddm, $f_{\mathrm{ddm}}^{\mathrm{ini}}$, ${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$, Ω_b, h, A_s, n_s), while the cosmological constant parameter Ω_Λ is adjusted to match the budget equation in a flat universe⁵.

If one varies the two decaying DM parameters (Γ_ddm, $f_{\mathrm{ddm}}^{\mathrm{ini}}$) while fixing the other parameters (${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$, Ω_b, h, A_s, n_s), one changes the predicted age of the Universe, which controls the amplitude of the matter power spectrum on all scales, as well as the redshift of radiation-to-matter equality, which determines the scale of the overall peak in the spectrum. These effects cause an enhancement of the matter power spectrum on scales larger than those crossing the Hubble radius around the time of equality, corresponding to comoving wavenumbers k < 2–3 × 10⁻³ h Mpc⁻¹, and a suppression on smaller scales. For wavenumbers k ≥ 10⁻¹ h Mpc⁻¹, the power spectrum is suppressed by a constant factor with respect to the ΛCDM case. A larger fraction, $f_{\mathrm{ddm}}^{\mathrm{ini}}$, or a higher rate, Γ_ddm, both imply a smaller amplitude of the power spectrum on these scales. Actually, the suppression factor is found to depend essentially on the product ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$, as is illustrated in Fig. 2. In the left panel, we vary $f_{\mathrm{ddm}}^{\mathrm{ini}}$ while keeping the product ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$ fixed to 0.005 Gyr⁻¹ (a value representative of the constraints found in the result Sect. 5.2). Then, the power spectrum of the 1b-DDM model is found to be independent of $f_{\mathrm{ddm}}^{\mathrm{ini}}$ up to the order of one per mille. Thus, we anticipate that the parameter $f_{\mathrm{ddm}}^{\mathrm{ini}}$ alone is difficult to constrain with Euclid data. Instead, in the right panel of Fig. 2, we vary the product ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$ while keeping $f_{\mathrm{ddm}}^{\mathrm{ini}}$ fixed. We clearly see a change in the suppression factor for k ≥ 10⁻¹ h Mpc⁻¹ and in the slope of the power spectrum for k ∼ 10⁻² h Mpc⁻¹, potentially detectable using Euclid probes.

Fig. 2. Ratio of the linear (solid lines) and non-linear (dashed lines) power spectra of several 1b-DDM models to that of a pure ΛCDM model with the same cosmological parameters, parameterised by the fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}$ and the decay rate Γddm. We work in the basis $(f_{\mathrm{ddm}}^{\mathrm{ini}},{\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}})$ to show that only the product of the two DDM parameters affects the linear power spectrum. The other parameters (${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$, Ωb, h, As, ns) are kept fixed, and the spectra are computed today (z = 0). The non-linear spectra are predicted by the emulator introduced in Sect. 3.2 and plotted up to the maximum wavenumber at which this emulator is trusted.

In Sect. 3.2, we compute the impact of the 1b-DDM model on the non-linear matter spectrum. In Sect. 5.2, we fit the 1b-DDM model to mock Euclid data. In order to obtain fast-converging MCMC chains, we adopt some flat priors on $f_{\mathrm{ddm}}^{\mathrm{ini}}$ and ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$, with prior edges detailed in Sect. 5.2, but we expect interesting constraints only on the second parameter.

2.3. Dark matter with two-body decay

In the next-to-simplest cosmological model of DDM, a cold DDM particle with a large mass, m_ddm, and a constant decay rate, Γ_ddm, is assumed to decay into a first massless daughter particle and a second massive daughter particle with mass m_daughter. This model is dubbed two-body decaying DM (2b-DDM). The parent particle is assumed to account for a fraction, $f_{\mathrm{ddm}}^{\mathrm{ini}}$, of the initial CDM budget, defined in the same way as for one-body decay (see Eq. (3)), with the remaining fraction $1-f_{\mathrm{ddm}}^{\mathrm{ini}}$ corresponding to ordinary stable CDM. In each decay, the fraction of energy transferred from the parent particle to the first massless daughter particle, ε, can be related to the mass ratio:

$$\begin{array}{c}\hfill \epsilon =\frac{1}{2}(1-\frac{m_{\mathrm{daughter}}^2}{m_{\mathrm{ddm}}^2}).\end{array}$$(4)

In the limit m_daughter → m_ddm, all the energy goes into the second massive daughter, but since this corresponds to the conversion of one CDM particle into another one, the model is indistinguishable from the standard ΛCDM model. In the opposite limit m_daughter → 0, the two daughter particles are ultra-relativistic and share the same amount of energy, which corresponds to ε = 0.5: this limit is equivalent to the 1-body decay model introduced in the previous section. However, in the more interesting range 0 < ε < 0.5, the second daughter particle can behave as WDM. Aoyama et al. (2014) have shown that for the purpose of computing cosmological observables one only needs to specify the three parameters (f_ddm, Γ_ddm, ε) on top of the usual ΛCDM parameters.

This model has been studied previously, for instance, in Aoyama et al. (2014), Vattis et al. (2019), Haridasu & Viel (2020), Franco Abellán et al. (2022, 2021), Schöneberg et al. (2022), Simon et al. (2022), or Bucko et al. (2024). It has also been invoked as a possible solution to the H₀ and/or S₈ tension. The best constraints at the moment come from CMB plus BAO data (Schöneberg et al. 2022), galaxy surveys (Simon et al. 2022), and WL surveys (Bucko et al. 2024; see Sect. 5.3). Finally, Franco Abellán et al. (2024) have shown how to perform efficient sensitivity forecasts based on machine learning techniques for this model.

At the level of linear perturbation theory, this model is implemented in a branch of CLASS developed and publicly released⁶ by the authors of Franco Abellán et al. (2021, 2022). Figure 3 shows the effect on the linear power spectrum of a variation in the parameters (Γ_ddm, ε, $f_{\mathrm{ddm}}^{\mathrm{ini}}$) for fixed ΛCDM parameters. We see that this model leads to a step-like suppression of the matter power spectrum, which is not surprising since, in this case, DM is split between a CDM and a WDM component. The shape of the step is, however, different from the CWDM case, because the warm component gets produced progressively and affects different scales at different times. Figure 3 focuses on cases with ε ≪ 0.5 for which, in each decay, most of the energy is transferred from one non-relativistic DM species to another one. Thus, while the universe expands, the energy density of total matter evolves almost like in the case of stable DM, ρ_m ∝ a⁻³, and the age of the universe is not significantly affected by the DDM parameters. This explains why in Fig. 3 we do not see any effect of the 2b-DDM parameters on the matter power spectrum on very large scales (k ≪ 10⁻¹ h Mpc⁻¹), as it was the case for 1b-DDM. As a side note, one can observe tiny oscillations in the linear power spectrum ratios of Figure 3. This is most likely a numerical artefact caused by the use of a fluid approximation for the perturbations of the warm species within a fixed range of scales. The same figure shows that these spurious oscillations are smoothed out by the emulator introduced in Sect. 3.3. Thus, they cannot affect our results⁷.

Fig. 3. Ratio of the linear (solid lines) and non-linear (dashed lines) power spectra of several two-body DDM models to that of a pure ΛCDM model with the same cosmological parameters, parameterised by the fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}$, the decay rate Γddm, and the fraction of energy ε going into the ultra-relativistic daughter particle at each decay. The parameters (${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$, Ωb, h, As, ns) are kept fixed, and the spectra are computed today (z = 0). The non-linear spectra are predicted by the emulator introduced in Sect. 3.3 and plotted up to the maximum wavenumber at which this emulator is trusted.

The parameter ε controls the velocity of the daughter particle just after the decay, which reads $v_{\mathrm{k}}=c\epsilon /\sqrt{1-2\epsilon }$ in the centre of mass frame (the subscript k refers to ‘kick’, since the daughter particles get a velocity kick). Thus, by analogy with WDM, ε determines the free-streaming scale and the location of the step in the power spectrum. The parameters (Γ_ddm, $f_{\mathrm{ddm}}^{\mathrm{ini}}$) both control the abundance of 2b-DDM as a function of time and thus the linear growth rate of the total DM density fluctuation δ_dm(a). Hence these parameters both control the amplitude of the step. The ΛCDM limit is recovered for ε = 0 and/or $f_{\mathrm{ddm}}^{\mathrm{ini}}=0$ and/or Γ_ddm = 0.

In Sect. 3.3, we show how to compute the impact of the 2b-DDM model on the non-linear matter spectrum. In Sect. 5.3, we fit the 2b-DDM model to mock Euclid data. We perform our sensitivity forecast with flat priors on $\{{log}_{10}{f}_{\mathrm{ddm}}^{\mathrm{ini}},{log}_{10}({\mathrm{\Gamma }}_{\mathrm{ddm}}/{\mathrm{Gyr}}^{-1}),{log}_{10}\epsilon \}$, with prior edges detailed in that section.

2.4. ETHOS n = 0

The ETHOS framework (Cyr-Racine et al. 2016) is a general attempt to parameterise physically plausible interactions in a dark sector featuring at least one type of non-relativistic relics (playing the role of cold interacting dark matter, IDM) and one type of ultra-relativistic relics (playing the role of interacting dark radiation, IDR). The theory provides a mapping between phenomenological parameters describing the relevant interaction rates and fundamental parameters appearing in the Lagrangian of the dark sector. In particular, the ETHOS index n describes to the power-law dependence of the IDR-IDM interaction rate Γ_{idr − idm} on the temperature of the dark sector.

The case n = 0 is of particular interest, because it corresponds to an IDM-IDR momentum exchange rate Γ_{idm − idr} scaling like the Hubble radius during radiation domination (Buen-Abad et al. 2015; Cyr-Racine et al. 2016; Becker et al. 2021). Thus, in this model, the ratio Γ_{idm − idr}/H (where both Γ_{idm − idr} and H depend on time) remains constant during radiation domination and decreases slowly during matter domination. This means that IDM and IDR can remain in a regime of feeble but steady interactions until equality. The IDR-IDM interactions then become gradually irrelevant at the beginning of matter domination and negligible during the formation of non-linear structures. On the other hand, ETHOS models with n > 0 tend to suppress the power spectrum very sharply below some scale. Thus, similar to pure WDM models, they are easier to constrain with Lyman-α data than with galaxy surveys.

This model can be motivated with some concrete and plausible particle physics set up, such as non-Abelian DM (Buen-Abad et al. 2015). It is interesting from the point of view of cosmology phenomenology because it introduces a very smooth suppression in the matter power spectrum (Lesgourgues et al. 2016; Buen-Abad et al. 2018) – instead of oscillatory patterns or an exponential cut-off as would be the case for ETHOS models with n > 0. The power spectrum suppression shape is also very different from the one caused by a hot or warm DM component. This model is often invoked as a solution to the S₈ tension (Lesgourgues et al. 2016; Buen-Abad et al. 2018) – or even to the Hubble tension, but this is no longer the case with recent data (Schöneberg et al. 2022). Current constraints on this model are obtained with CMB data combined with Lyman-α data (Archidiacono et al. 2019; Hooper et al. 2022) or with the full-shape power spectrum of the BOSS galaxy redshift survey (Rubira et al. 2023; see Sect. 5.4).

This model can be parameterised in terms of the IDR-IDM scattering amplitude, Γ_{idr − idm}(z_*), at some arbitrary reference redshift z_*, of the density of DM (through Ω_idmh²), and of the density of DR (through Ω_idrh²). Following the rest of the literature, we choose a reference redshift z_* = 10⁷ and express the effective comoving rate of IDR scattering off IDM as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{idr}-\mathrm{idm}}(z_{\ast })=-{\mathrm{\Omega }}_{\mathrm{idm}}{h}^{2}c{a}_{\mathrm{dark}}.\end{array}$$(5)

Assuming IDR with a thermal spectrum and two fermionic degrees of freedom, we can parameterise the IDR density in terms of the IDR-to-photon temperature ratio, T_idr/T_γ = ξ_idr ≤ 1, such that ${\mathrm{\Omega }}_{\mathrm{idr}}=\frac{7}{8}{\xi }_{\mathrm{idr}}^4{\mathrm{\Omega }}_{\gamma }$. The contribution of IDR to the effective number of neutrinos is then given by ΔN_eff = (T_idr/T_ν)⁴ with T_ν defined in the instantaneous neutrino decoupling limit; that is, ΔN_eff = (11/4)^4/3ξ_idr⁴ ≃ 3.85 ξ_idr⁴. The ratio ξ_idr is a dimensionless parameter. Γ_{idr − idm} is a rate and a_dark is an inverse distance that we express in Mpc⁻¹ (this definition and choice of units has no other purpose than matching the conventions of the CLASS code and of previous work studying this model)⁸. Finally, in the ETHOS framework, one needs to specify the self-interaction rate between IDR particles. The non-Abelian DM model and the CMB+Lyman-α constraints of Lesgourgues et al. (2016), Buen-Abad et al. (2018), Archidiacono et al. (2019), or Hooper et al. (2022) assumed a strongly self-interacting IDR fluid. One may assume instead free-streaming IDR, and Rubira et al. (2023) consider the two cases. These two different assumptions are expected to have a small impact on CMB constraints (due to the effect of IDR fluctuations dragging the photons fluctuations before decoupling) but a negligible impact on constraints from large-scale structure (because IDR self-interactions are irrelevant for the growth rate of IDM). Here we stick to the assumption of free-streaming IDR.

The most important physical effect of this model on the matter power spectrum comes from the fact that the IDR-IDM interactions tend to slow down the growth rate of DM fluctuations on sub-Hubble scales during radiation domination, and to suppress the power spectrum on small scales at all subsequent times (Lesgourgues et al. 2016; Buen-Abad et al. 2015). Actually, as is mentioned in Archidiacono et al. (2019), the power spectrum suppression is mainly sensitive to the effective comoving scattering rate of IDR off IDM, which is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{idm}-\mathrm{idr}}=\frac{4{\rho }_{\mathrm{idr}}}{3{\rho }_{\mathrm{idm}}}{\mathrm{\Gamma }}_{\mathrm{idr}-\mathrm{idm}}.\end{array}$$(6)

Since ρ_idr is proportional to ξ_idr⁴ while ρ_idm is normalised by the measurement of Ω_idmh², this rate is controlled mainly by the parameter combination a_dark ξ_idr⁴. Therefore, we expect the amplitude of the suppression in the linear matter power spectrum to depend strongly on a_dark ξ_idr⁴ and weakly on the orthogonal combination, except in the case of sufficiently large ξ_idr⁴, in which the effect of additional radiation with a given ΔN_eff also comes into play. Indeed, an enhancement of ΔN_eff has some well-known effects on the matter and CMB power spectra, explained for instance in Lesgourgues & Verde (2022), and we expect Euclid to be sensitive to this effect (Euclid Collaboration: Archidiacono et al. 2025).

The ETHOS formalism is implemented in the public version of CLASS⁹. We show in Fig. 4 the effect of varying the parameters ξ_idr or a_dark ξ_idr⁴ with fixed values of all other cosmological parameters.

Fig. 4. Ratio of the linear (solid lines) and non-linear (dashed lines) power spectra of several free-streaming ETHOS n = 0 models to that of a pure ΛCDM model with the same cosmological parameters, parameterised by the dark-radiation-to-photon temperature ratio ξidr and interaction strength adark. The effects are displayed in the basis (ξidr, adark ξidr4) to show that the combination adark ξidr4, which gives the scattering rate of IDR off IDM, controls the amplitude of the small-scale suppression of the linear matter power spectrum. The other parameters (${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$, Ωb, h, As, ns) are kept fixed, and the spectra are computed today (z = 0). The non-linear spectra are predicted by the emulator introduced in Sect. 3.4 and plotted up to the maximum wavenumber at which this emulator is trusted.

In the left panel, the scattering rate controlled by a_dark ξ_idr⁴ is fixed, which explains the constant suppression of the linear power spectrum in the large-k limit. When log₁₀(ξ_idr) varies from −1.2 to −0.6, ΔN_eff increases from 6.1 × 10⁻⁶ to 0.015, which is too small to directly affect the matter power spectrum. However, these different values of ξ_idr and thus a_dark have an impact on intermediate scales: they control the maximum scale at which IDM feels the interaction, and thus the wavenumber at which the matter power spectrum starts to be suppressed. When log₁₀(ξ_idr) reaches −0.4, the radiation density gets enhanced by a non-negligible amount, ΔN_eff = 0.097. This results in an additional suppression of the linear power spectrum on small scales.

In the right panel, the amount of IDR is fixed to a small value but the effective scattering rate is increased, leading to more and more suppression. This suppression has a different shape to the case of WDM: it behaves like a transition to a smaller spectral index rather than an exponential cut-off.

In Sect. 3.4, we show how to compute the impact of the ETHOS n = 0 model on the non-linear matter spectrum. In Sect. 5.4, we fit this model to mock Euclid data. We perform our sensitivity forecast with flat priors on {log₁₀(a_darkξ_idr⁴/Mpc⁻¹), log₁₀ξ_idr}, with prior edges detailed in that section.

3. Emulating the non-linear evolution

To predict observable WL and galaxy correlation spectra, one needs to know the non-linear matter power spectrum for each model. Since N-body simulations are computationally too expensive for being run at each point in MCMC chains, it is customary to use a restricted set of N-body simulations to build emulators of the non-linear matter power spectrum. These emulators should be accurate compared to the sensitivity of the experiment within the range of model parameters covered by our priors, and fast to evaluate within MCMC runs. In this section, we describe the emulators used in our MCMC forecasts for each of the four non-minimal DM models described in Sect. 2.

Instead of directly emulating the non-linear power spectrum of the extended cosmological model, $P_{\mathrm{model}}^{\mathrm{nl}}(k,z)$, it is customary to emulate the ratio

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{model}}(k,z)=\frac{P_{\mathrm{model}}^{\mathrm{nl}}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nl}}(k,z)},\end{array}$$(7)

and to compute the final observable spectra using

$$\begin{array}{c}\hfill P_{\mathrm{model}}^{\mathrm{nl}}(k,z)=P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nl}}(k,z){\mathcal{S}}_{\mathrm{model}}(k,z).\end{array}$$(8)

This strategy offers two main advantages. Firstly, it is easier to generate accurate training data for the ratio 𝒮_model(k,z) than for the final spectrum, since several N-body simulation artefacts tend to cancel out in the ratio (e.g. resolution effects at small scale, or residual noise from cosmic variance and mesh assignment on large scale). Secondly, the final spectrum depends on all cosmological and DM parameters, but the ratio 𝒮_model(k,z) does not in some cases. This ratio depends on course on the DM parameters, but not necessarily on each single parameter of the ΛCDM model. In each model, one can perform some explicit tests to investigate this dependence and build the emulator on a reduced parameter space.

In this work, we need to decide which tool we should use for predicting the first factor in Eq. (8); that is, the spectrum $P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nl}}(k,z)$. In principle, we could use fitting functions like Halofit (Smith et al. 2003) or HMcode 2020 (Mead et al. 2021), emulators like the EuclidEmulator2 (Knabenhans et al. 2021) or BACCOemulator (Angulo et al. 2021), etc. In the future, when analysing real data, we shall use the best tool available at that time in order to get unbiased results. But for the purpose of the present work, which is to forecast the sensitivity to DM parameters, one could use essentially any of these tools without changing the results on the DM parameter sensitivity, provided that the same tool is used when generating fiducial data and when fitting theoretical predictions. Our choice shall be specified in the next sections.

In the context of this work, having accurate predictions for the ratio 𝒮_model(k,z) is more important. With a noisy emulator, one could get slightly wrong predictions for the effect of DM parameters on the final observable spectra, and potentially underestimate degeneracies between these parameters and cosmological or baryonic feedback parameters. In the forecasts presented here, we use emulators designed to achieve per-cent level accuracy up to k ∼ 𝒪(10) h Mpc⁻¹ and z ∼ 2.5 (in the next section we provide further details on the accuracy of each emulator). Given the sensitivity of Euclid, this is sufficient for robust forecasts. There are some plans to keep training these emulators and improving their accuracy in order to be sure that, when analysing real data, the error coming from the emulator is clearly subdominant in the total systematic error budget.

3.1. Cold plus warm dark matter

To predict the non-linear suppression in the matter power spectrum in CWDM scenarios, we use an improved version of the emulator already described in Parimbelli et al. (2021). Such an emulator is trained on a large set of N-body simulations, covering a large parameter space, for a total of 100 models with different WDM fractions f_wdm and WDM masses. The simulations explicitly assume thermal WDM, but this assumption is not relevant in the final analysis: as long as one performs the mass conversion described in Sect. 2.1 before calling the emulator, the latter still applies to all models in which WDM has a Fermi–Dirac distribution possibly rescaled by a factor χ. The simulations cover masses down to $m_{\mathrm{wdm}}^{\mathrm{thermal}}=0.03\mathrm{keV}$, but we have checked that the emulator provides a consistent extrapolation down to $m_{\mathrm{wdm}}^{\mathrm{thermal}}=0.01\mathrm{keV}$ for small fractions $f_{\mathrm{ddm}}^{\mathrm{ini}}$ (see Hervas-Peters et al. 2024). For each model, four realisations are run with fixed amplitudes: two with different random phases and two with the opposite phases. The box size is set to 120 h⁻¹ Mpc in order to reconnect with the linear regime at large scales for all redshifts and without any significant discontinuity and to obtain percent-level convergence up to k ≈ 10 h Mpc⁻¹. The (fixed) cosmological parameters are Ω_m = 0.315, Ω_b = 0.049, h = 0.674, n_s = 0.965, and a value of A_s that would give σ₈ = 0.811 in the pure ΛCDM limit (where σ₈ is the square root of the variance of matter fluctuations in spheres of radius 8 h⁻¹ Mpc).

Initial conditions are set at z = 99 with a modified version of the N-GenIC code (Springel et al. 2005), using a linear power spectrum obtained from CLASS (Blas et al. 2011). The simulations are run with the tree-particle mesh (TreePM) code GADGET-III (Springel et al. 2005) and follow the gravitational evolution of 512³ particles. Snapshots are taken starting from z = 3.5 down to z = 0, linearly spaced with Δz = 0.5. Once the power spectra from these snapshots are measured, we take their ratio with respect to the corresponding ΛCDM spectrum and build the emulator following the exact same procedure as in Parimbelli et al. (2021). This new tool emulates the first 20 principal components of the power spectrum suppression using Gaussian processes. It is trained on the redshift range z ∈ [0 − 3.5] and in the range of scales k ∈ [0.07 − 25] h Mpc⁻¹. The performances are found to be comparable to the ones stated in Parimbelli et al. (2021); that is, the difference between the emulated and the simulated suppressions never exceeds ∼1.5%. All in all, the non-linear matter power spectrum in the presence of CWDM is given by

$$\begin{array}{c}\hfill P_{\mathrm{\Lambda }\mathrm{CWDM}}^{\mathrm{nl}}(k,z)=P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nl}}(k,z){\mathcal{S}}_{\mathrm{CWDM}}(k,z),\end{array}$$(9)

where the last term is precisely what the emulator predicts and $P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nl}}(k,z)$ is computed with the version of Halofit revisited by Takahashi et al. (2012) and Bird et al. (2012).

We plot a few examples of predictions for the non-linear spectrum at z = 0 (compared to the linear predictions of CLASS) in Fig. 1. We can clearly see that the suppression of power induced by the WDM component on small scales is much smaller in the non-linear (rather than linear) power spectrum. This is a well-known effect of mode-mode coupling when perturbations become non-linear. The smaller is the redshift, the less pronounced is the power spectrum suppression on scales smaller than the maximum free-streaming scale.

A few considerations about the simulations must be made here. For the sake of computational efficiency, all the particles in all the realisations are initialised as cold particles, even in the runs containing WDM. This assumption has a twofold implication. First, we assume that the differences between a CWDM model and ΛCDM reside in the initial conditions and in their linear power spectra; second, we are neglecting WDM thermal velocities. We tested the impact of these two assumptions by running a further realisation, with f_wdm = 0.2 and $m_{\mathrm{wdm}}^{\mathrm{thermal}}=0.13\mathrm{keV}$, in which we initialise 512³ CDM particles as well as 512³ more particles as Type2, with the correct thermal velocities¹⁰. This value of f_wdm has been chosen because, below this fraction, current data are compatible with any value for m_wdm; the value of the mass has been chosen in order to have a ∼50% suppression in the linear power spectrum at k ∼ 5 h Mpc⁻¹. We show the results of this test in Fig. 5. In the left plot, we compare the matter power spectrum suppression at various redshifts when neglecting thermal velocities (solid orange lines) and when fully considering them (dashed violet lines). As can be noted, differences between the two treatments are only relevant at z ≳ 2 and for k ≳ 5 h Mpc⁻¹. The right plot shows instead the ratios between the angular power spectra of cosmic shear (or equivalently WL, orange), the cross-correlation between GC and galaxy lensing (red), and GC (purple), computed according to the prescriptions described in Sect. 4, using each of the two sets of power spectra in the left plot. We use a single bin here for simplicity, ranging from z = 0 to z = 3.5, and neglect intrinsic alignment. Differences are well below percent level; for comparison, at ℓ = 10⁴, the Euclid sample variance is expected to be ∼1.6%. We can conclude that our assumptions do not introduce any systematic effects in the analysis.

Fig. 5. Left: Effect of neglecting the WDM thermal velocities in CWDM simulations with fwdm = 0.2 and $m_{\mathrm{wdm}}^{\mathrm{thermal}}=0.13\mathrm{keV}$. In each panel, solid orange lines represent the suppression in the non-linear matter power spectrum when neglecting WDM thermal velocities; dashed violet lines do the same when implementing WDM as a second fluid in the simulation, with its own thermal velocity field. We plot as vertical lines the mean interparticle separation in blue and the Nyquist frequency in red. Right: Ratio of angular power spectra C(ℓ) for cosmic shear (orange), the cross-correlation of GC and galaxy lensing (red), and GC (purple), defined in Eq. (30), and computed using either the power spectra that neglect or consider thermal velocities. These C(ℓ) are computed for simplicity in a single redshift bin ranging from z = 0 to 3.5 with the galaxy distribution of Eq. (34). We show the 0.25% and 0.5% regions as dark and light shaded areas. The maximum ℓ value corresponding to the optimistic and pessimistic settings for Euclid are drawn as vertical lines for each probe (cosmic shear or equivalently WL in dotted yellow, GC and cross-correlation in dotted violet).

3.2. Dark matter with one-body decay

We employed the fitting functions found by Hubert et al. (2021) to model the non-linear matter power spectrum in the presence of one-body decay. These fits are inspired by fitting functions published in Enqvist et al. (2015) and built upon N-body simulations implementing DDM into the PKDGRAV3 code (Potter et al. 2017).

We have seen in Sect. 2.2 that 1b-DDM induces a suppression in the linear matter power spectrum that is asymptotically constant on intermediate and small scales, with a suppression factor proportional to ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$, or to $f_{\mathrm{ddm}}^{\mathrm{ini}}/{\tau }_{\mathrm{ddm}}$. The amplitude and redshift dependence of this suppression factor is given by

$$\begin{array}{c}\hfill {\epsilon }_{\mathrm{lin}}(z)=\alpha f_{\mathrm{ddm}}^{\mathrm{ini}}{\left(\frac{\mathrm{Gyr}}{{\tau }_{\mathrm{ddm}}}\right)}^{\beta }{\left(\frac{1}{0.105z+1}\right)}^{\gamma },\end{array}$$(10)

where α, β, γ are functions of ω_b := Ω_bh², h, and ω_m := Ω_bh² + Ω_dmh². We refer to Hubert et al. (2021) and Bucko et al. (2023) for their detailed form. We note that the suppression functions ε_lin(z) and ε_nonlin(k, z) introduced respectively in Eqs. (10, 11) should not be confused with the parameter ε of the 2b-DDM model. The non-linear evolution imprints an additional suppression that can be inferred from N-body simulations. Enqvist et al. (2015) provided a fit to the non-linear suppression function ε_nonlin(k, z) in the case $f_{\mathrm{ddm}}^{\mathrm{ini}}=1$ that Hubert et al. (2021) generalised to arbitrary values of the DDM fraction. The suppression function is estimated from N-body simulations for a fixed cosmology. Since only late-time DM decays are of interest, the initial conditions of such N-body simulations are identical to those in a ΛCDM scenario. However, to account for the 1b-DDM, the particle masses are being gradually decreased in the simulation as a function of the rate Γ_ddm, the fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}$ and the simulation time, mimicking the decay process (for more details, see Hubert et al. 2021). The suite of N-body simulations used to construct the fitting functions was run with a box size of 500 h⁻¹ Mpc evolving 1024³ particles. The cosmological parameters were fixed to fiducial values Ω_m = 0.307, Ω_b = 0.048, 10⁹A_s = 2.43, h = 0.678, and n_s = 0.96. The convergence of the 1b-DDM N-body simulations was studied in Hubert et al. (2021) with the conclusion that the implementation of the model is trustworthy at least up to k ≃ 6.4 h Mpc⁻¹. Finally, Hubert et al. (2021) argue that ε_nonlin(k, z) is nearly cosmology-independent and can be extrapolated to cosmologies well beyond those probed in our work.

The fitting function provides the suppression of the matter power spectrum with respect to the fiducial ΛCDM cosmology, P_{1b − DDM}(k, z)/P_ΛCDM(k, z) = 1 − ε_nonlin(k, z), with

$$\begin{array}{c}\hfill {\epsilon }_{\mathrm{nonlin}}(k,z)=\frac{1+a_1{(k/{\mathrm{Mpc}}^{-1})}^{p_1}}{1+a_2{(k/{\mathrm{Mpc}}^{-1})}^{p_2}}{\epsilon }_{\mathrm{lin}}(z).\end{array}$$(11)

The suppression function interpolates from the linear behaviour on intermediate scales, ε_nonlin(k, z)→ε_lin(z), to a power-law suppression on small scales with ε_nonlin(k, z)∝k^{p₁ − p₂}. The factors a₁, a₂, p₁, and p₂ are given for each lifetime τ_ddm and redshift z by

$$\begin{array}{cc}\hfill a_1& =0.7208+2.027\left(\frac{\mathrm{Gyr}}{{\tau }_{\mathrm{ddm}}}\right)+\frac{3.031}{1+1.1z}-0.18,\hfill \\ \hfill a_2& =0.0120+2.786\left(\frac{\mathrm{Gyr}}{{\tau }_{\mathrm{ddm}}}\right)+\frac{0.6699}{1+1.1z}-0.09,\hfill \\ \hfill p_1& =1.045+1.225\left(\frac{\mathrm{Gyr}}{{\tau }_{\mathrm{ddm}}}\right)+\frac{0.2207}{1+1.1z}-0.099,\hfill \\ \hfill p_2& =0.992+1.735\left(\frac{\mathrm{Gyr}}{{\tau }_{\mathrm{ddm}}}\right)+\frac{0.2154}{1+1.1z}-0.056.\hfill \end{array}$$(12)

These fitting functions are publicly available as a part of the DMemu package,¹¹ and designed to reproduce the results of N-body simulations with a precision better than 1% up to k = 13 h Mpc⁻¹. We note that, at a given redshift, the fitting functions of Eq. (12) depend only on τ_ddm (or Γ_ddm), while ε_lin depends only on ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$. Thus, the non-linear evolution lifts the degeneracy between Γ_ddm and $f_{\mathrm{ddm}}^{\mathrm{ini}}$ observed at the level of the linear power spectrum.

In order to match the linear predictions of CLASS on large and intermediate scales with those of the fitting functions on intermediate and small scales without introducing any discontinuity, we use the following ansatz to calculate the non-linear matter power spectrum of the 1b-DDM model:

$$\begin{array}{cc}\hfill P_{1\mathrm{b}-\mathrm{DDM}}(k,z)& =P_{1\mathrm{b}-\mathrm{DDM},\mathrm{lin}}(k,z)\frac{P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM},\mathrm{lin}}(k,z)}\hfill \\ \hfill & \times \frac{1-{\epsilon }_{\mathrm{nonlin}}(k,z)}{1-{\epsilon }_{\mathrm{lin}}(z)},\hfill \end{array}$$(13)

with the non-linear ΛCDM spectrum evaluated with the version of Halofit revisited by Takahashi et al. (2020) and Bird et al. (2012). Then, firstly, on intermediate (linear) scales, the second and third factor in the right-hand side of Eq. (13) go to one, and one recovers P_{1b − DDM}(k, z)→P_{1b − DDM, lin}(k, z). Secondly, on smaller (non-linear) scales, after noticing that we can rewrite Eq. (13) as

$$\begin{array}{cc}\hfill P_{1\mathrm{b}-\mathrm{DDM}}(k,z)& =\frac{P_{1\mathrm{b}-\mathrm{DDM},\mathrm{lin}}(k,z)}{1-{\epsilon }_{\mathrm{lin}}(z)}\frac{P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM},\mathrm{lin}}(k,z)}\hfill \\ \hfill & \times [1-{\epsilon }_{\mathrm{nonlin}}(k,z)],\hfill \end{array}$$(14)

and that the first fraction tends towards P_ΛCDM,lin(k,z), we get

$$\begin{array}{c}\hfill P_{1\mathrm{b}-\mathrm{DDM}}(k,z)⟶P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)[1-{\epsilon }_{\mathrm{nonlin}}(k,z)],\end{array}$$(15)

that is, the approximation to the non-linear 1b-DDM power spectrum provided by the emulator. Equation (13) is designed to provide a smooth transition between these two limits. According to this ansatz, the ratio P_{1b − DDM}(k, z)/P_ΛCDM(k, z) is given by the boost factor

$$\begin{array}{c}\hfill {\mathcal{S}}_{1\mathrm{b}}(k,z)=\frac{P_{1\mathrm{b}-\mathrm{DDM},\mathrm{lin}}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM},\mathrm{lin}}(k,z)}\frac{1-{\epsilon }_{\mathrm{nonlin}}(k,z)}{1-{\epsilon }_{\mathrm{lin}}(z)}.\end{array}$$(16)

We already saw in Fig. 2 the ratio of 1b-DDM-to-ΛCDM linear power spectra, as well as the ratio of non-linear spectra given by Eq. (16). Figure 6 is similar to Fig. 2 but shows additionally the raw result of the emulator; that is, the ratio P_{1b − DDM}(k, z)/P_ΛCDM(k, z)≃1 − ε_nonlin(k, z) above k ≳ 0.05 h Mpc⁻¹ (dashed lines). The linear prediction (solid lines) and the raw emulator (dashed lines) match each other quite well around k = 0.05 h Mpc⁻¹, but switching abruptly from one to the other at a given wavenumber would introduce a small discontinuity in the spectrum. Dotted lines show the boost factor defined in Eq. (16) and used in our pipeline. This factor provides a very smooth interpolation from the prediction of CLASS to that of the emulator.

Fig. 6. Effect of 1b-DDM parameters on the linear (solid) and non-linear (dashed) matter power spectrum. Left: Effect of varying the decay rate Γddm with a fixed fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}=1$. Right: Effect of varying the fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}$ with a fixed decay rate Γddm = (1/13.5) Gyr−1. The other parameters (${\mathrm{\Omega }}_{\mathrm{dm}}^{\mathrm{ini}}$, Ωb, h, As, ns) are kept fixed, and the spectra are computed today (z = 0). Dashed lines show the predictions of the emulator of Hubert et al. (2021). Our pipeline relies on the prescription of Eq. (16), shown as a dotted line, which smoothly interpolates from the linear to non-linear behaviour.

3.3. Dark matter with two-body decay

To model the two-body decays up to non-linear scales, we use the emulator published in Bucko et al. (2024), which can provide the 2b-DDM-to-ΛCDM non-linear power spectrum ratio

$$\begin{array}{c}\hfill {\mathcal{S}}_{2\mathrm{b}}(k,z)=\frac{P_{2\mathrm{b}-\mathrm{DDM}}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)}\end{array}$$(17)

up to z ≃ 2.3 and k ≃ 6 h Mpc⁻¹. The emulator was trained on approximately 100 PKDGRAV3N-body simulations directly implementing the late-time DM decays, while starting from ΛCDM-like initial conditions at z_ini = 49. Bucko et al. (2024) set L_box = 125, 250, 512 h⁻¹ Mpc and N = 256³, 512³, 1024³ depending of each specific DDM configuration, in such way to achieve converged simulations up to k_max = 6h Mpc⁻¹. Bucko et al. (2024) argue that the suppression 𝒮_2b(k, z) is approximately independent of cosmology and fix the standard cosmological parameters to Ω_m = 0.307, Ω_b = 0.048, 10⁹A_s = 2.43, h = 0.678, and n_s = 0.96 in the simulations. At each simulation time step, a number of DM particles is randomly selected for decay. The decay into a lighter daughter particle is accounted for through a velocity kick with amplitude v_k and random direction. In the limit ε ≪ 0.5 considered here, v_k is approximately given by c ε. These kicks lead to suppression in the matter power spectrum below the free-streaming length of the massive daughter particles controlled by v_k ∼ c ε.

The emulator predicts 𝒮_2b(k, z) using a combination of a ‘Principal Component Analysis’ (PCA) with feed-forward ‘sinusoidal representation networks’ (SIRENs; see Sitzmann et al. (2020)). Within the emulation process, the PCA is used to compress the power spectrum ratios 𝒮_2b(k, z), taking into account 5 principal components. Then, the SIREN architecture is trained in a supervised fashion to predict these principal components given the input parameters of 2b-DDM model and the redshift of interest. The loss function of the network is the square distance of the input and output 2b-DDM-to-ΛCDM ratio, reconstructed from the PCA components predicted by the network. The emulator covers the case of an arbitrary fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}\in [0,1]$ of long-lived DDM particles with ${\tau }_{\mathrm{ddm}}:={\mathrm{\Gamma }}_{\mathrm{ddm}}^{-1}\ge 13.5\mathrm{Gyr}$ decaying into non-relativistic daughters with v_k ≲ 5000 km s⁻¹, corresponding to ε < 0.017. The emulator can predict ratios of 2b-DDM and ΛCDM nonlinear matter power spectra up to z = 2.3 and k ≃ 6 h Mpc⁻¹, with a precision better than 1% at the 68% CL. It is implemented inside the publicly available DMemu package introduced after Eq. (12). We already compared the emulator result to the linear 2b-DDM-to-ΛCDM linear power spectrum ratio in Fig. 3.

Like in other cases, the final non-linear power spectrum of the 2b-DDM model is obtained by mutiplying the non-linear power spectrum of the ΛCDM model (computed using Halofit) with the emulated ratio 𝒮_2b(k, z).

3.4. ETHOS n = 0

The non-linear matter power spectrum of the ETHOS n = 0 model is predicted by a dedicated emulator that will be presented in Bucko et al. (in prep.). Like for the 1b-DDM and 2b-DDM cases, this emulator will be released within the DMemu package. It assumes the particular case in which IDR consists of two free-streaming fermionic degrees of freedom. It predicts the ETHOS-to-ΛCDM non-linear power spectrum ratio

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{ETHOS}n=0}(k,z)=\frac{P_{\mathrm{ETHOS}\mathrm{n}=0}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)}\end{array}$$(18)

up to z = 3 and k ≃ 5 h Mpc⁻¹. The architecture used to train the ETHOS emulator is similar to the one used in the 2b-DDM scenario, described in Sect. 3.3, with slight modifications. First of all, only 4 PCA components are used to compress the input ETHOS-to-ΛCDM matter power spectra, while the SIREN architecture involves two dense hidden layers with 256 neurons each. The emulator provides below 1% errors at the aforementioned scales and redshifts, within the range of ETHOS models defined by the N-body simulations discussed in the next paragraphs.

The emulator is built upon a suite of N-body simulations which have been run using PKDGRAV3 with L_box = 325 h⁻¹ Mpc and N = 512³, assuming a fiducial cosmology with ω_idm := Ω_idmh² = 0.1202, ω_b = 0.02236, h = 0.6727, n_s = 0.9649, and 10⁹A_s = 2.101. One massive neutrinos species with m_ν = 0.06 eV was included. Instead of using a typical back-scaling approach to generate the initial conditions, Bucko et al. (in prep.) follow an alternative method described in Tram et al. (2019). The “true” initial conditions are generated using the C0NCEPT code (Dakin et al. 2022). In combination with CLASS, C0NCEPT also computes the linear evolution of all species (photons, metric, neutrinos, IDR, IDM). This information is used to calculate the gravitational potential at each time step in the PKDGRAV3 simulation. In this way, the DM particles of the simulation feel their own gravity, taken into account at the non-linear level, plus the gravity from other species, modelled at the linear level.

The simulations used to train the emulator implement the ETHOS n = 0 only through modified initial conditions at z_ini = 49. The effect of IDM-IDR scattering at z < z_ini is neglected. This assumption is valid only for model parameters such that the scattering rate Γ_{idm − idr} is negligible compared to the Hubble rate at z = z_ini. Since the rate Γ_{idm − idr} is computed with respect to conformal time, it should be compared to the conformal Hubble rate ℋ = aH. Rubira et al. (2023) provide an analytical approximation for the (redshift-dependent) interaction rate to Hubble rate ratio,

$$\begin{array}{cc}\hfill \frac{{\mathrm{\Gamma }}_{\mathrm{idm}-\mathrm{idr}}}{\mathcal{H}}& \simeq 0.0152\left(\frac{a_{\mathrm{dark}}}{1000{\mathrm{Mpc}}^{-1}}\right){\left(\frac{{\xi }_{\mathrm{idr}}}{0.1}\right)}^4\hfill \\ \hfill & \times \frac{{(1+z)}^2}{{[{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\gamma }{(1+z)}^4(1+{\xi }_{\mathrm{idr}}^4)+{\mathrm{\Omega }}_{\mathrm{\Lambda }}]}^{1/2}}.\hfill \end{array}$$(19)

Assuming Ω_m = 0.27 and ξ_idr⁴ ≪ 1, this gives approximately Γ_{idm − idr}/ℋ ≃ a_darkξ_idr⁴/(0.48 Mpc⁻¹) at z = 49. Bucko et al. (in prep.) suggest to trust the simulations and the emulator as long as this ratio is smaller than 0.1. In first approximation, this is the case for a_darkξ_idr⁴ < 0.05 Mpc⁻¹. We shall see in Sect. 5.4 that this region is appropriate to study the sensitivity of Euclid to ETHOS parameters, at least when the fiducial model is assumed to be ΛCDM (or close to it).

We plot the magnitude of the ratio given by Eq. (19) at z = 49 as a function of (ξ_idr, a_dark) in the left panel of Fig. 7. The region where the emulator is to be trusted lays below the solid black line. Dashed lines in the left panel of Fig. 7 correspond to models with either log₁₀(a_dark/Mpc⁻¹) = 1.0 (dashed) or ξ_idr = 0.25 (dash-dotted), for which we show the emulator predictions in the right panel of the same figure. Namely, the dashed curves show the power spectrum suppression as a function of ξ_idr, ranging from the ΛCDM limit for ξ_idr = 0.0 up to ξ_idr = 0.4, while fixing log₁₀(a_dark/Mpc⁻¹) to 1.0. Similarly, the dash-dotted lines show the emulator output for ξ_idr = 0.25 varying the coupling strength a_dark from 10⁻³ Mpc⁻¹ to 10^9.5 Mpc⁻¹. We note that some of the cases shown in the right panel stand outside of the region Γ_{idm − idr}/ℋ ≤ 0.1 in which the emulator is to be fully trusted.

Fig. 7. Left: Ratio of the interaction rate between IDM and IDR (Γidm − dr) and comoving Hubble rate (ℋ) as a function of the dark-radiation-to-photon temperature ratio ξidr and interaction strength adark, computed at the redshift zini = 49 at which the N-body simulations used to construct the ETHOS emulator are initialised. We display the contours of equal ratio as solid white lines and highlight the threshold value of 0.1 in black. We further depict the region with adark = 10 Mpc−1 (dashed grey) and ξidr = 0.25 (dash-dotted grey). Right: Power spectrum suppression 𝒮ETHOS(k, z) predicted by the emulator for parameters chosen along each of the two grey lines of the left panel.

The effect of the ETHOS parameters on the non-linear power spectrum is also shown in Fig. 4 in the (ξ_idr, a_darkξ_idr⁴) basis. While at the linear level the suppression of the matter power spectrum is mainly controlled by the combination a_darkξ_idr⁴, we see that at the non-linear level ξ_idr also plays a significant role for fixed a_darkξ_idr⁴.

3.5. Baryonic feedback

Baryonic feedback processes can alter the gas distribution around DM halos, causing a deviation between the total matter distribution and the distribution of DM (e.g. Chisari et al. 2018; van Daalen et al. 2020). These processes induce a suppression of the total matter power spectrum P_m on small scales that may be somewhat similar to the effect of non-minimal DM and induce degeneracies between baryonic and DM parameters (see e.g. Hubert et al. 2021). Hence, for our purposes, it is crucial to incorporate these processes into our modelling framework. In this study, we use the BCemu framework¹² (Giri & Schneider 2023) to address this concern. The BCemu framework serves as an emulator for the suppression 𝒮_bf(k) caused by baryonic feedback. Consequently, the total non-linear matter power spectrum in a given cosmological model, P_m, can be expressed as

$$\begin{array}{c}\hfill P_{\mathrm{m}}(k,z)={\mathcal{S}}_{\mathrm{bf}}(k,z)P_{\mathrm{m},\mathrm{no}\mathrm{bf}}(k,z),\end{array}$$(20)

where P_m, no bf(k, z) is the total matter power spectrum neglecting baryonic feedback effects at redshift z.

BCemu has been used in several recent WL studies (e.g. Schneider et al. 2022; Grandis et al. 2024). It is based on the baryonic correction modelling framework of Schneider & Teyssier (2015), Schneider et al. (2019), and Giri & Schneider (2021). This framework parameterises the stellar and gas profiles at a given redshift with seven baryonic parameters. Giri & Schneider (2021) analysed these baryonic parameters and found that three parameters are enough to model the suppression seen in hydrodynamical simulations at scales k ≲ 10 h Mpc⁻¹ and at a given redshift. We use this three-parameter model in this work.

Two of these parameters describe the gas profile in halos of given virial radius, r_vir, and virial mass, M_vir, modelled as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{gas}}(r)\propto \frac{{\mathrm{\Omega }}_{\mathrm{b}}/{\mathrm{\Omega }}_{\mathrm{m}}-f_{\mathrm{star}}(M_{\mathrm{vir}})}{{[1+10\frac{r}{r_{\mathrm{vir}}}]}^{\beta (M_{\mathrm{vir}})}{[1+\frac{r}{{\theta }_{\mathrm{ej}}{r}_{\mathrm{vir}}}]}^{\frac{2}{5}[7-\beta (M_{\mathrm{vir}})]}},\end{array}$$(21)

with a total stellar fraction, f_star(M_vir), and a mass-dependent index,

$$\begin{array}{c}\hfill \beta (M_{\mathrm{vir}})=\frac{3M_{\mathrm{vir}}/M{\prime }_{\mathrm{c}}}{1+M_{\mathrm{vir}}/M{\prime }_{\mathrm{c}}}.\end{array}$$(22)

The former function is assumed to be known,

$$\begin{array}{c}\hfill f_{\mathrm{star}}(M_{\mathrm{vir}})=0.055{\left(\frac{{10}^{11.3}{h}^{-1}{M}_{\odot }}{M_{\mathrm{vir}}}\right)}^{0.2}.\end{array}$$(23)

Thus, in this model, the gas profile only depends on two free parameters: a critical mass M′_c such that small halos with M_vir ≪ M_c have a gas profile shallower than the Navarro–Frenk–White profile, and an ejection factor θ_ej giving the ratio of the gas ejection radius to the virial radius. The BCemu model also involves assumptions concerning the stellar profile of the central galaxy. The fraction of stars in the central galaxy, f_cga(M_vir), is given by a relation similar to f_star(M_vir), but with a different exponent,

$$\begin{array}{c}\hfill f_{\mathrm{cga}}(M_{\mathrm{vir}})=0.055{\left(\frac{{10}^{11.3}{h}^{-1}{M}_{\odot }}{M_{\mathrm{vir}}}\right)}^{0.2+{\eta }_{\delta }},\end{array}$$(24)

where the index η_δ is an additional free parameter.

In summary, the minimal BCemu model relies on three free parameters (M′_c, θ_ej, η_δ) impacting, respectively, the overall suppression induced by baryonic feedback, the maximum scales affected by the suppression, and the upturn of 𝒮_bf(k,z) at large k. In order to deal only with dimensionless parameters, Schneider et al. (2022) redefine the first one as M_c := M′_c/(1 h⁻¹ M_⊙). Figure 2 in Schneider et al. (2019) shows the impact of these parameters on the matter power spectrum. In the BCemu model, the only cosmology dependence of 𝒮_bf(k,z) comes through the baryon fraction Ω_b/Ω_m.

In particular, we assume no explicit dependence of the baryonic feedback suppression function 𝒮_bf(k,z) on the parameters describing non-standard DM models. This assumption was shown to be valid at least for k < 5 h Mpc⁻¹ in the CWDM scenario (see Sect. 3.4 in Parimbelli et al. 2021). This conclusion is expected to apply also to the other DM scenarios studied here in which, like in the CWDM case, DM particles are decoupled at low redshift and behave either as CDM or WDM. In our analysis, smaller scales with k > 5 h Mpc⁻¹ only have a small contribution to the spectra $C_{\mathit{ij}}^{\mathrm{XY}}(\ell )$ involving the first two WL redshift bins. Thus, the findings of Parimbelli et al. (2021) suggest that we can safely neglect the impact of non-standard DM on baryonic feedback. More generally, we can think of the effect of non-standard DM and of BF on the non-linear matter power spectrum as two leading-order effects, of a few percent each within the range of scales relevant in our analysis, and that of non-standard DM on BF as a next-to-leading order effect of a few percent squared; that is, a few per mille. It is thus reasonable to neglect this correction in a first analysis.

We refer interested readers to Giri & Schneider (2021) for a more detailed description of the BCemu parameters.

The redshift evolution of 𝒮_bf(k,z) is modelled by making each of the three baryonic parameters b redshift dependent as

$$\begin{array}{c}\hfill b(z)=b(0){(1+z)}^{-{\nu }_b},b\in \{{\log }_{10}{M}_c,{\theta }_{\mathrm{ej}},{\eta }_{\delta }\},\end{array}$$(25)

where ν_b is a free parameter. This leads to a total of six parameters to model the baryonic feedback. In our choice of fiducial values and priors, we restrict the values of ν_b such that the baryonic parameters {log₁₀M_c, θ_ej, η_δ} remain within the range of the parameter space where BCemu is trained.

While it is known that the WL signal is modified at small scales by baryonic feedback effects, the situation is much less clear regarding the GC signal. Since galaxies act as tracers of the underlying DM distribution, they are not directly affected by the ejection of gas via feedback processes. We rather expect an indirect effect caused by the relaxation of the DM potential reacting to the ejection of gas. Since we do not know the true amplitude of this indirect effect, we consider two extreme cases where the GC is either unchanged by baryonic feedback or it is affected in the same way as the WL. We expect the truth to lie somewhere between these two cases.

4. Forecast methodology

4.1. Likelihood

We use a standard formalism to describe the Euclid photometric likelihood already presented, for instance, in Audren et al. (2013b), Euclid Collaboration: Blanchard et al. (2020), Euclid Collaboration: Archidiacono et al. (2025), Casas et al. (2024). The galaxy images of the WL survey and the galaxy positions of the GC photometric survey are binned into N redshift bins. In each bin, the raw data can be processed into two-dimensional spherical maps of either the lensing potential field density in the WL case or the galaxy density field in the GCph case. The maps are decomposed into spherical harmonics with coefficients a_i(ℓ, m) for each redshift bin i. Each a_i is assumed to obey a Gaussian distribution with covariance matrix $C_{\mathit{ij}}(\ell )=\frac{1}{2\ell +1}{\displaystyle \sum _m}{a}_i(\ell ,m){[a_i(\ell ,m)]}^{\ast }$. The C_ij(ℓ) are the observed power spectra of WL or GCph in harmonic space and can be compared to theoretical predictions.

In our forecasts, we assume that the power spectra observed by Euclid coincides with the theoretical predictions of a given fiducial cosmology with spectra $C_{\mathit{ij}}^{\mathrm{fid}}(\ell )$ arising from multipoles a_i^fid such that $C_{\mathit{ij}}^{\mathrm{fid}}(\ell )=\frac{1}{2\ell +1}{\displaystyle \sum _m}{a}_i^{\mathrm{fid}}(\ell ,m){[a_i^{\mathrm{fid}}(\ell ,m)]}^{\ast }$. The likelihood ℒ of the observed data given a theoretical model with spectrum $C_{\mathit{ij}}^{\mathrm{th}}(\ell )$ is then given by

$$\begin{array}{cc}\hfill \mathcal{L}& =\mathcal{N}{\displaystyle \prod _{\ell ,m}}{[det{\mathsf{C}}^{\mathrm{th}}(\ell )]}^{-1/2}\hfill \\ \hfill & \times exp\{-f_{\mathrm{sky}}\frac{1}{2}{\displaystyle \sum _{\mathit{ij}}}{a}_i^{\mathrm{fid}}(\ell ,m){(C_{\mathit{ij}}^{\mathrm{th}})}^{-1}(\ell ){[a_j^{\mathrm{fid}}(\ell ,m)]}^{\ast }\},\hfill \end{array}$$(26)

where 𝒩 is a normalisation factor, and partial sky coverage is approximately accounted for through multiplication with the sky fraction f_sky. This can be rewritten as (Audren et al. 2013b)

$$\begin{array}{cc}\hfill {\chi }^2& :=-2ln\frac{\mathcal{L}}{{\mathcal{L}}_{\max }}\hfill \\ \hfill & =f_{\mathrm{sky}}{\displaystyle \sum _{\ell }}(2\ell +1)\{\mathrm{Tr}[{({\mathsf{C}}^{\mathrm{th}})}^{-1}(\ell ){\mathsf{C}}^{\mathrm{fid}}(\ell )]\hfill \\ \hfill & +ln\frac{det{\mathsf{C}}^{\mathrm{th}}(\ell )}{det{\mathsf{C}}^{\mathrm{fid}}(\ell )}-N\}\hfill \\ \hfill & =f_{\mathrm{sky}}{\displaystyle \sum _{\ell }}(2\ell +1)\{\frac{d_{\ell }^{\mathrm{mix}}}{d_{\ell }^{\mathrm{th}}}+ln\frac{d_{\ell }^{\mathrm{th}}}{d_{\ell }^{\mathrm{fid}}}-N\},\hfill \end{array}$$(27)

where N is the size of the matrices C^th(ℓ) and C^fid(ℓ), while

$$\begin{array}{c}\hfill \mathsf{C}(\ell ):=\left[\begin{array}{cc}{C}_{\mathit{ij}}^{\mathrm{LL}}(\ell )& C_{\mathit{ij}}^{\mathrm{GL}}(\ell )\\ C_{\mathit{ij}}^{\mathrm{LG}}(\ell )& C_{\mathit{ij}}^{\mathrm{GG}}(\ell )\end{array}\right],d_{\ell }:=det\mathsf{C}(\ell )\end{array}$$(28)

for each of the theoretical and fiducial spectra. Finally, the mixed determinant is defined as

$$\begin{array}{c}\hfill d_{\ell }^{\mathrm{mix}}:={\displaystyle \sum _{k=1}^N}det\left[\{\begin{array}{cc}{C}_{\mathit{ij}}^{\mathrm{th}}(\ell )\hfill & \mathrm{for}j\ne k\hfill \\ C_{\mathit{ij}}^{\mathrm{fid}}(\ell )\hfill & \mathrm{for}j=k\hfill \end{array}\right],\end{array}$$(29)

such that in each term of the sum, the determinant is evaluated over a matrix in which the k-th column of the theory matrix C^th has been substituted by the k-th column of the fiducial matrix C^fid.

We then perform MCMC forecasts (Audren et al. 2013b; Casas et al. 2024) using this likelihood. The likelihood is incorporated into the MontePython package¹³ (Audren et al. 2013a; Brinckmann & Lesgourgues 2019) for Bayesian parameter inference. The role of MontePython is to fit the fiducial spectra under the assumption of a given theoretical model with a set of free parameters. A few independent Monte Carlo Markov Chains sample the likelihood by exploring the parameter space according to the Metropolis-Hastings algorithm, until some convergence criterium is reached. The best-fit model coincides by construction with the fiducial model, while the marginalised credible interval of each parameter provide an estimate of the sensitivity of Euclid to this parameter.

4.2. Observable power spectra

The model for the spectra $C_{\mathit{ij}}^{\mathit{XY}}$ used in the likelihood, where X = L (respectively X = G) refers to the WL (respectively GC) probe and i = 1, …, N_i to the bin number, is detailed in Euclid Collaboration: Blanchard et al. (2020), Euclid Collaboration: Archidiacono et al. (2025), Casas et al. (2024). The final expression is given by

$$\begin{array}{cc}\hfill C_{\mathit{ij}}^{\mathit{XY}}(\ell )& ={\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\mathrm{d}z\frac{W_i^X(k(\ell ,z),z)W_j^Y(k(\ell ,z),z)}{c^{-1}H(z)r^2(z)}\hfill \\ \hfill & \times P_{\mathrm{m}}(k(\ell ,z),z)\hfill \\ \hfill & +N_{\mathit{ij}}^{\mathit{XY}}(\ell ),\hfill \end{array}$$(30)

where W_i^X(k, z) are the window functions of the X probe, H(z) is the Hubble rate at redshift z, r(z) the comoving distance to an object at redshift z, P_m(k, z) the matter power spectrum evaluated at wavenumber k, and $N_{\mathit{ij}}^{\mathit{XY}}(\ell )$ the noise spectrum. The boundaries z_min and z_max, defined in Table 1, specify the redshift range covered by the survey. The relation k(ℓ, z) is inferred from the Limber approximation (Kaiser 1992; Kilbinger et al. 2017),

$$\begin{array}{c}\hfill k(\ell ,z)=\frac{\ell +1/2}{r(z)},\end{array}$$(31)

which is sufficiently accurate for ℓ > ℓ_min, where ℓ_min is given in Table 1 (see however Tanidis & Camera 2019). Assuming a Poissonian distribution of galaxies, the noise spectra read

$$\begin{array}{c}\hfill N_{\mathit{ij}}^{\mathrm{LL}}=\frac{{\sigma }_{ϵ}^2}{{\overline{n}}_i}{\delta }_{\mathit{ij}},N_{\mathit{ij}}^{\mathrm{GG}}=\frac{1}{{\overline{n}}_i}{\delta }_{\mathit{ij}},N_{\mathit{ij}}^{\mathrm{LG}}=N_{\mathit{ij}}^{\mathrm{GL}}=0,\end{array}$$(32)

where ${\overline{n}}_i$ is the expected average number of galaxies per steradian in the i-th bin, and σ_ϵ² is the variance of the observed ellipticities, also given in Table 1. The galaxy field that GCph measures is assumed to be a linear tracer of the underlying matter field, such that the galaxy power spectrum is given by P_g(k, z) = b²(z)P_m(k, z) with some bias function b(z). Here, for simplicity, we neglect additional effects on the photometric galaxy power spectrum such as lensing magnification or redshift-space distortions (Yoo et al. 2009; Bonvin & Durrer 2011; Challinor & Lewis 2011; Yoo & Zaldarriaga 2014), although these effects are expected to play a non-negligible role in the analysis of real Euclid data (see Lepori et al. 2022 and Tanidis et al. 2024). Sticking to linear bias is conservative as long as we rely on pessimistic assumptions concerning the minimum angular scale or maximal multipole $l_{\max }^{\mathrm{GCph}}$ described in Table 1. In the optimistic case, we should be aware that non-linear biasing may come into play on the smallest scales used in the analysis, and introduce a possible degeneracy with DM parameters that is neglected here.

Then the GC window functions reads

$$\begin{array}{c}\hfill W_i^{\mathrm{G}}(z)=\frac{n_i(z)H(z)b(z)}{c},\end{array}$$(33)

where n_i(z) is the observed galaxy density distribution normalised to unit area in redshift bin i. Since there is no reliable model for b(z), it is modelled as a step-like function given by b(z) = b_i in the redshift range z_i⁻ < z < z_i⁺ of redshift bin i. Each b_i is treated as a free nuisance parameter and marginalised over in the forecast. Taking photometric redshift errors into account, the observed distribution of galaxies n_i(z) in bin i is given by the true galaxy distribution,

$$\begin{array}{c}\hfill n(z)=n_0{\left(\frac{z}{z_0}\right)}^{2}exp[-{\left(\frac{z}{z_0}\right)}^{1.5}],\end{array}$$(34)

with $z_0=z_{\mathrm{mean}}/\sqrt{2}$, and by the redshift error probability distribution,

$$\begin{array}{cc}\hfill p_{\mathrm{ph}}(z_{\mathrm{p}}|z)& =\frac{1-f_{\mathrm{out}}}{\sqrt{2\pi }{\sigma }_{\mathrm{b}}(1+z)}exp\{-\frac{1}{2}{\left[\frac{z-c_{\mathrm{b}}{z}_{\mathrm{p}}-z_{\mathrm{b}}}{{\sigma }_{\mathrm{b}}(1+z)}\right]}^2\}\hfill \\ \hfill & +\frac{f_{\mathrm{out}}}{\sqrt{2\pi }{\sigma }_0(1+z)}exp\{-\frac{1}{2}{\left[\frac{z-c_0{z}_{\mathrm{p}}-z_{\mathrm{b}}}{{\sigma }_0(1+z)}\right]}^2\}.\hfill \end{array}$$(35)

The normalised distribution n_i(z) then reads (Ma et al. 2005; Joachimi & Schneider 2009; Joachimi & Bridle 2010)

$$\begin{array}{c}\hfill n_i(z)=\frac{n(z){\displaystyle {\int }_{z_i^{-}}^{z_i^{+}}}\mathrm{d}{z}_{\mathrm{p}}{p}_{\mathrm{ph}}(z_{\mathrm{p}}|z)}{{\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\mathrm{d}\stackrel{\sim }{z}n(\stackrel{\sim }{z}){\displaystyle {\int }_{z_i^{-}}^{z_i^{+}}}\mathrm{d}{z}_{\mathrm{p}}{p}_{\mathrm{ph}}(z_{\mathrm{p}}|\stackrel{\sim }{z})}.\end{array}$$(36)

The parameters entering this model are listed in Table 1. The WL window functions are given by

$$\begin{array}{c}\hfill W_i^{\mathrm{L}}(k,z)=W_i^{\gamma }(z)-{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m}}\frac{{\mathcal{F}}_{\mathrm{IA}}(z)}{D(k,z)}{W}_i^{\mathrm{IA}}(z),\end{array}$$(37)

where the latter term corrects for intrinsic alignment (IA) effects, W_i^IA(z) = c⁻¹n_i(z)H(z), and W_i^γ(z) is the shear-only window function

$$\begin{array}{cc}\hfill W_i^{\gamma }(z)& =\frac{3}{2}{c}^{-2}{H}_0^2{\mathrm{\Omega }}_{\mathrm{m}}(1+z)r(z)\hfill \\ \hfill & \times {\displaystyle {\int }_z^{z_{\max }}}\mathrm{d}z\prime n_i(z\prime )[1-\frac{r(z)}{r(z\prime )}].\hfill \end{array}$$(38)

D(k, z) is the linear growth factor, defined as D(k, z)≔[P_m, lin(k, z)/P_m, lin(k, z = 0)]^1/2. In linear ΛCDM cosmology, perturbations grow independently of scale and the k dependence exactly cancels out. Instead, the particle DM models considered in this work lead to some scale-dependent linear growth, such that the function D is a function of (k, z). The factor ℱ_IA is modelled as

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{IA}}(z)={(1+z)}^{{\eta }_{\mathrm{IA}}}{[⟨L⟩(z)/L_{\ast }(z)]}^{{\beta }_{\mathrm{IA}}}\end{array}$$(39)

and depends on the mean galaxy luminosity divided by a characteristic luminosity ⟨L⟩(z)/L_*(z), which is read from the file scaledmeanlum_E2SA.dat provided by the authors of Euclid Collaboration: Blanchard et al. (2020). In practice, we vary η_IA and 𝒜_IA as nuisance parameters but fix β_IA and 𝒞_IA due to the strong degeneracies between the former and the latter.

We note that Eq. (38) is derived under the assumption that the non-relativistic matter density scales like a⁻³: the factor (1 + z) in front of the integral actually comes from the product ρ_m(z) a²(z). In the ΛCDM, CWDM, and ETHOS models, the assumption ρ_m ∝ a⁻³ is excellent (as long as massive neutrino effects are neglected). In the 2b-DDM case, it is still excellent since we are only interested in the limit ε ≪ 1 in which the decays convert a negligible fraction of the non-relativistic energy density ρ_m into relativistic energy density ρ_r. However, in the 1b-DDM case, the product ρ_m a³ decreases slightly between the highest and lowest redshift probed by the survey, which spans an interval of proper time Δt. The relative variation in ρ_m a³ over this interval is given by $f_{\mathrm{ddm}}^{\mathrm{ini}}{\mathrm{\Gamma }}_{\mathrm{ddm}}\mathrm{\Delta }t$, and remains below a few percent for the 1b-DDM models studied in the next sections. Thus we neglect this sub-dominant effect and stick to Eq. (38)¹⁴.

In our forecast we rely either on a pessimistic or optimistic assumption concerning the range of scales at which our model is trusted and data is included. In the pessimistic case, we truncate the data at ${\ell }_{\max }^{\mathrm{WL}}=1500$ for WL and ${\ell }_{\max }^{\mathrm{GC}}=750$ for GCps. In the optimistic case, we use ${\ell }_{\max }^{\mathrm{WL}}=5000$ for WL and ${\ell }_{\max }^{\mathrm{GC}}=1500$ for GCps.

The matter power spectrum P_m(k, z) that appears in Eq. (30) is usually computed in four steps. First, we call a Boltzmann code to compute the linear matter power spectrum of a ΛCDM model with the same cosmological parameters as the non-standard DM model we are interested in. Second, we ask the same Boltzmann code to use a standard algorithm to infer the non-linear power spectrum for this ΛCDM model. In the forecasts of this work, for simplicity, we use the version of Halofit revisited by Takahashi et al. (2012) and Bird et al. (2012) as a baseline, or HMcode 2020 (Mead et al. 2021) in the case where neutrinos are assumed to have a mass of 0.06 eV. Third, we use one of the emulators described in Sects. 3.1 to 3.4 to transform this into a non-linear power spectrum for the non-standard DM model of interest. Fourth, when baryonic feedback corrections need to be taken into account, we call the emulator described in Sect. 3.5 to add baryonic corrections. For the WL auto-correlation spectra, $C_{\mathit{ij}}^{\mathrm{LL}}$, the matter power spectrum of Eq. (30) always includes baryonic feedback. For the GCph auto-correlation spectra, $C_{\mathit{ij}}^{\mathrm{GG}}$, we consider the two cases in which the power spectrum incorporates such corrections or not. We note that for the cross-correlation spectra, when baryonic feedback is included in WL but not GCph, we use for $C_{\mathit{ij}}^{\mathrm{LG}}(\ell )=C_{\mathit{ji}}^{\mathrm{GL}}$,

$$\begin{array}{cc}\hfill C_{\mathit{ij}}^{\mathrm{LG}}(\ell )& ={\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\mathrm{d}z\frac{W_i^{\mathrm{L}}(k(\ell ,z),z)W_j^{\mathrm{G}}(k(\ell ,z),z)}{c^{-1}H(z)r^2(z)}\hfill \\ \hfill & \times \sqrt{P_{\mathrm{m}}^{\mathrm{BF}}(k(\ell ,z),z)P_{\mathrm{m}}^{\mathrm{no}\mathrm{BF}}(k(\ell ,z),z)}\hfill \\ \hfill & +N_{\mathit{ij}}^{\mathrm{LG}}(\ell ).\hfill \end{array}$$(40)

4.3. Boltzmann code

We need to call a Boltzmann code for two purposes: first, to compute the comoving distance-redshift relation r(z) and the (scale-dependent) growth factor D(k, z); and second, to compute the non-linear matter power spectrum P_m(k, z). However, the strategy of the emulators described in Sects. 3.1–3.4 implies a calculation of the non-linear matter power spectrum for the equivalent ΛCDM sharing the same value of the standard cosmological parameters {ω_b, ω_cdm, h, A_s, n_s} as the non-standard DM model of interest. Instead, the distance-redshift relation and the scale-dependent growth factor should be computed according to the background and linear theory equations describing the true non-standard DM model. We solve this issue by calling the Boltzmann code twice at each point in parameter space: first for the non-standard DM model with linear output, to infer r(z) and D(k, z); and second for the equivalent ΛCDM with Halofit or HMcode 2020 corrections switched on, to get the desired non-linear matter power spectrum.

We choose to use CLASSv 3.2 (Lesgourgues 2011; Blas et al. 2011) as our Boltzmann code since the CWDM, 1b-DDM and ETHOS models are implemented in the main public branch of the code¹⁵, while 2b-DDM is implemented in a public but separate branch class_decays¹⁶.

When calling CLASS, we should specify a maximum wavenumber k_max. A given Fourier mode k of a field observed at redshift z projects under a given angle θ contributing mainly to a mutipole ℓ, with the relation between k, ℓ, and z given by Eq. (31). Thus the choice of k_max should reflect the maximum multipole ℓ_max and minimum redshift z_min contributing to the power spectra in a given analysis. Using the relation k(ℓ, z) at fixed ℓ, one can express each $C_{\mathit{ij}}^{\mathit{XY}}(\ell )$ as an integral over k rather than z, and plot the cumulative contribution of different k values to $C_{\mathit{ij}}^{\mathit{XY}}(\ell )$. To make a robust choice for k_max, we show in Fig. 8 the contribution of different k values to the power spectra of the first redshift bins, C₀₀^LL(ℓ) and C₀₀^GG(ℓ). The figure shows that in the pessimistic case (${\ell }_{\max }^{\mathrm{WL}}=1500$), choosing k_max = 10 Mpc⁻¹ is sufficient to include 99.5% of the contribution the C₀₀^LL(ℓ), and a fortiori to all other spectra. For the optimistic case (${\ell }_{\max }^{\mathrm{WL}}=5000$), we find that k_max = 30 Mpc⁻¹ is sufficient.

Fig. 8. Cumulative contribution of different k values to $C_{\mathit{ij}}^{\mathit{XY}}(\ell )$ for a given ℓ in the nearest redshift bin (ij = 00). Left: Case of WL, XY = LL. Right: Case of GC, XY = GG. For each ℓ, 99.5% of the contribution stands below the black isocontour. In the pessimistic case, we include all values of (ℓ,k) on the left of and below the dashed grey lines in the calculation of our observables; in the optimistic case, on the left and below the dashed orange line. Thus, we always include at least 99.5% of the contribution to each $C_{\mathit{ij}}^{\mathit{XY}}(\ell )$.

The $C_{\mathit{ij}}^{\mathit{XY}}(\ell )$ are computed with a trapezoidal integral over 200 values of z on a linear grid between z_min and z_max. The integral is performed for discrete values of ℓ, with a logarithmically spaced grid of 100 values of ℓ between ℓ_min and ℓ_max. Finally, a second-order spline interpolation is used to get the spectra for every integer ℓ.

4.4. Parameters and priors

We list in Table 2 the free parameters used in our forecasts (not including the DM parameters specific to each model, which shall be specified in each Sect. 5). The table also provides the assumed fiducial values and priors. We remind that forecast errors are anyway nearly independent of the chosen fiducial values, especially for the ΛCDM parameters, which have nearly Gaussian posteriors. The first five parameters are the cosmological parameter of the standard ΛCDM model. The following six parameters describe the baryonic feedback model implemented in BCemu. The last 12 nuisance parameters account for linear bias in each bin and intrinsic alignment parameters.

We use a flat prior on each of these parameters. For the first three BCemu parameters and the two intrinsic alignment parameter, we pass explicit prior edges to ensure that these parameters remain a range making sense physically. For the other parameters, as long as we only perform parameter inference with the Metropolis-Hastings algorithm, it is strictly equivalent to pass to MontePython some very remote prior edges – such that the chains never reach the prior boundaries – or to require the code to use non-informative priors (abbreviated as n.i. in Table 2), in which case the code achieves the same feature automatically. We note that our chains never reach the prior edges passed for the intrinsic alignement parameters, so we are effectively using non-informative priors also for 𝒜_IA, η_IA.

Finally, theoretical predictions depend on a few additional parameters that are usually kept fixed because the set of free nuisance parameters from Table 2 are sufficient to account for the uncertainty on the model. For intrinsic alignment, we fix the parameter β_IA defined in Eq. (39) to β_IA = 2.17; for baryonic feedback, we fix the BCemu parameters μ = 1, γ = 2.5, δ = 7, η = 0.2.

5. Results and discussion

5.1. Cold plus warm dark matter

Main results. For the CWDM model, we perform forecasts using the free parameters $m_{\mathrm{wdm}}^{\mathrm{thermal}}$ and f_wdm introduced in Sect. 2.1. Our fiducial values and priors are summarised in Table 3 for these parameters and Table 2 for all other free parameters. The fiducial model is chosen to be a pure ΛCDM model. As was already stated in Sect. 2.1, we use a linear prior on the mass and a logarithmic prior on the WDM fraction (that is, a flat priors on its logarithm). The latter choice allows us to explore the constraining power of Euclid for very small WDM fractions (see also Schneider et al. 2020). This limit is particularly interesting to study with Euclid since, in this case, Euclid bounds can be competitive with respect to Lyman-α bounds (Hooper et al. 2022).

Indeed, Lyman-α data probe smaller scales than WL and GC surveys, and can in principle better constrain models with a large mass. However, when the WDM fraction is small, the effect of WDM on the power spectrum is also small. Then, it could be unconstrained by Lyman-α data, and if WDM is light enough its effects can manifest themselves on relatively large scales. In this case, the precise measurement of the power spectrum at larger scales with Euclid remains decisive.

Our main results are summarised in Fig. 9, where we show the marginalised 95% credible intervals for f_wdm and m_wdm. The horizontal axis can be interpreted as the thermal WDM mass (bottom axis) or as the Dodelson–Widrow mass (top axis, see Sect. 2.1 for details). The different colours of the contours mark the different probes that have been used to obtain the constraints. These contours are already marginalised over all other cosmological and nuisance parameters (including baryonic feedback parameters). For each of the six cases shown in Fig. 9, we ran 36 chains summing up to ∼1.4 millions of steps (MS) in each optimistic case or ∼2.6 MS in each pessimistic case. The Gelman-Rubin convergence criterium (Gelman & Rubin 1992) reached about |R − 1|∼0.002 for most parameters, with a worse value of ∼0.008 for a few parameters.

Fig. 9. Left: Edges of the 95% credible interval on the WDM mass mwdm and fraction fwdm for the CWDM model, with pessimistic assumptions and three data combinations: WL alone, WL plus GC from the photometric survey (3 × 2pt), and 3 × 2pt combined with Planck CMB data. For the 3 × 2pt and 3 × 2pt + Planck data sets, baryonic feedback has been assumed to affect the WL power spectrum but not the GC power spectrum. The posterior is marginalised over other cosmological parameters, baryonic feedback parameters, and nuisance parameters (accounting for bias uncertainty and intrinsic alignment). The model is equivalent to pure ΛCDM towards the lower horizontal axis (small fwdm) and right vertical axis (large mwdm). The forecast assumes a flat prior on the mass of thermal WDM (lower axis) and a logarithmic prior on the WDM fraction (left axis), but we show the relation to Dodelson–Widrow masses in the upper axis (see Sect. 2.1 for definitions). Right: Same with optimistic assumptions.

The constraints from WL are considerably looser (by about a factor five) in the pessimistic rather than optimistic case. Indeed, in the pessimistic case, the WL data is only fitted up to ℓ_max = 1500, while models with a thermal mass of a few hundreds of keV only affect larger multipoles. As long as one sticks to pessimistic assumptions, adding information from GC (in the photometric survey) and on clustering-lensing correlations (3 × 2pt) makes a small difference, because in this case the clustering information is taken into account only up to ℓ_max = 750. The addition of CMB data from Planck also has a very small impact, given that CMB is sensitive to the clustering properties of pressureless DM (including WDM) only at second order in perturbations, through CMB lensing effects – as was explained in Voruz et al. (2014). Figure 9 shows that in the pessimistic case, Euclid 3 × 2pt + Planck data have a potential to rule out masses $m_{\mathrm{wdm}}^{\mathrm{thermal}}\lesssim 280\mathrm{eV}$ (95%CL) in the extreme case where these particles make up the totality of DM, or $m_{\mathrm{wdm}}^{\mathrm{thermal}}\lesssim 75\mathrm{eV}$ for f_wdm = 0.1 (95%CL). We note that even with pessimistic assumptions, the WDM mass can be constrained even for WDM fractions slightly below 0.1, while current bounds from high-resolution Lyman-α data cannot distinguish models with f_wdm = 0.1 from the pure ΛCDM limit (Hooper et al. 2022).

The picture drastically improves when one assumes optimistic settings with ℓ_max = 5000 for WL and ℓ_max = 1500 for GC. The data are then able to probe the presence of WDM with a much smaller value of the maximum free-streaming scale; that is, a larger mass. For the same WDM fraction, using WL data alone, the mass bounds become approximately five times tighter in the optimistic case. Despite of its limitation to ℓ_max ≤ 1500, GC data turns out to be very sensitive to the suppression induced by WDM even with a large mass, such that the 3 × 2pt probe is about twice more sensitive than the WL probe alone. However, the combination with Planck data makes no difference also in that case – at least when the mock data is assumed to account for a pure CDM model. The reason is that, for the large WDM masses that remain compatible with the data, the maximum free-streaming scale of WDM is very low, such that even CMB lensing is unaffected by the suppression induced by WDM. In the optimistic case, the Euclid 3 × 2pt probe has a potential to rule out all WDM masses with $m_{\mathrm{wdm}}^{\mathrm{thermal}}\lesssim 930\mathrm{eV}$ for f_wdm = 1 and $m_{\mathrm{wdm}}^{\mathrm{thermal}}\lesssim 230\mathrm{eV}$ for f_wdm = 0.1. It can constrain the mass even when f_wdm is as low as a few times 10⁻²; in other words, when only a few percent of the total DM is warm. This region of parameter space is far from current Lyman-α bounds, and even future Lyman-α surveys are unlikely to probe such small WDM fractions.

It is still unclear whether the final Euclid sensitivity will be closer to that of our pessimistic or optimistic forecast. At least, we expect that these two forecasts are bracketing the true constraining power of the future data. We shall see that, compared to other non-minimal DM models discussed in the next sections, CWDM is particularly sensitive to the choice of a cut-off multipole ℓ_max. This is due to the step-like nature of the effect of WDM on the matter power spectrum: up to a given wavenumber, the ΛCDM and CWDM models are strictly equivalent, and then the power drops. This means that the constraining power of a data set on the CWDM parameters depends more on the minimum scale (and thus maximum multipole and redshift) included in the analysis than on the actual error bars on the power spectrum. As was discussed above, this is particularly true for large values of f_wdm; for tiny WDM fractions, the precision with which the power spectrum is constrained remains crucial.

Importance of baryonic feedback. In Fig. 10, we evaluate the impact of different assumptions concerning baryonic feedback effects. We compare the bounds derived from the 3 × 2pt probe only under three assumptions. The baseline case (orange contours) is the same as in our previous discussion and in Fig. 9: the six nuisance parameters describing baryonic feedback are marginalised over, and baryonic feedback is assumed to affect only the WL probe; that is, the total matter power spectrum. The grey contours are derived assuming instead that baryonic feedback affects the two probes (WL and GC) in the same way: in other words, the same BCemu corrections are applied to the total matter and galaxy power spectra. Finally, the magenta contours were obtained with fixed rather than marginalised baryonic feedback parameters: they account for the unrealistic situation in which baryonic feedback effects would be perfectly known, given some independent measurements.

Fig. 10. Left: Same as Fig. 9, but only for the 3 × 2pt dataset and with different assumptions on baryonic feedback (BF): fixed BF (magenta), BF affecting only the WL power spectrum (orange), or BF affecting both the WL and GC power spectra (grey). The “truth” is expected to lay between the latter two cases (orange and grey), which give anyway very similar results. Right: Same with optimistic assumptions.

In principle, introducing more freedom in baryonic physics may result in looser constraints on WDM parameters due to parameter degeneracies. On the other hand, it is not obvious that such degeneracies are present due to the different shape of the effects imprinted by either WDM or baryons, not only as a function of scale but also as a function of redshift. In particular, the redshift dependence of the DM-induced suppression is reversed compared to the one from baryonic effects; its amplitude gets smaller at smaller redshift, due to non-linear clustering and mode-mode coupling, while overall the opposite is true for baryons, at least in most of the redshift range probed by Euclid.

We first compare the orange and magenta contours. In the pessimistic case (left panel), we find that introducing more freedom in the baryonic model and marginalising over baryonic feedback parameters does degrade a bit the bounds on WDM parameters. However, this is no longer true in the optimistic case, which proves that the information contained in the data at large mutipoles is sufficient to disentangle between the physical effects of WDM and baryonic feedback. This underlines the wealth of cosmological information encoded in the deep non-linear regime of the power spectrum.

We now switch to the comparison between the orange and grey contours. The impact of baryonic feedback on the galaxy power spectrum is not understood and modelled as well as its impact on the total matter powers spectrum. However, as was discussed in Sect. 3.5, we expect baryonic effects to be smaller in the galaxy powers spectrum – or at least not bigger. Thus, the true sensitivity of Euclid should lay somewhere between the forecasts corresponding to the grey and orange contours. However, these contours are very close to each other in both pessimistic and optimistic cases. This is likely due to the reverse redshift dependence of WDM and baryonic effects, which allows GC data to discriminate between them. Such a test validates the bounds on CWDM parameters obtained in the previous paragraphs with baryonic feedback included only in the WL probe.

Comparison with current bounds. It is interesting to compare the expected sensitivity of Euclid to current bounds on CWDM parameters obtained by Hervas-Peters et al. (2024) using an existing WL survey, the Kilo Degree Survey (KiDS – see Asgari et al. 2021). We start by comparing the sensitivity of the Euclid WL probe alone with that of KiDS. A fair comparison requires similar priors. However, the KiDS analysis assumes logarithmic priors on both the WDM mass and fraction, with slightly different prior edges than in our previous analysis. Given that Bayesian credible intervals do depend on priors, especially when constraining some parameters describing a model extension that is not required by the data, we repeat some of our forecasts with different top-hat priors matching exactly the ones of KiDS: log₁₀f_wdm ∈ [log₁₀(0.005), log₁₀(1)] and ${log}_{10}(m_{\mathrm{wdm}}^{\mathrm{thermal}}/\mathrm{keV})\in [{log}_{10}(10),{log}_{10}(1.5)]$. The results are shown in Fig. 11.

Fig. 11. Comparison of our Euclid forecasts for the WL-only probe with current bounds from KiDS. In this particular case, we switch to the same top-hat prior on ${log}_{10}{f}_{\mathrm{wdm}}^{\mathrm{ini}}$ and ${log}_{10}(m_{\mathrm{wdm}}^{\mathrm{thermal}}/\mathrm{eV})$ as in the KiDS analysis of Hervas-Peters et al. (2024). In the pessimistic case, we also adopt the same baryonic feedback recipe as Hervas-Peters et al. (2024) with a marginalisation over three baryonic feedback parameters (instead of six in our baseline treatment).

In the pessimistic case, we find that the Euclid sensitivity is not so different from that of KiDS. This result may sound surprising, given the much larger number of galaxy images expected from Euclid. There are two reasons for this.

The first reason is that, as was explained before, in the case of CWDM, the bounds depend a lot on the minimum scale (or maximum wavenumber k_max) included in the analysis, more than on the error bar on the measured power spectrum. It also depends on the maximum redshift of the data, since the signature of WDM is more clear at high redshift and partially washed out at small redshift. We note also that a given multipole ℓ probes smaller wavenumbers at high redshift, as was shown in Eq. (31). In the KiDS analysis and in our Euclid forecast with pessimistic assumptions, the data is conservatively cut at ℓ_max = 1500. Thus, within the redshift range covered by both experiments, the maximum wavenumber at each redshift k_max(ℓ_max, z) is the same, and the sensitivity to WDM is roughly similar despite of the smaller Euclid error bars. For instance, the highest redshift bin of KiDS peaks around z ∼ 1.1, probing up to k ∼ 1 h Mpc⁻¹. Euclid adds information at higher redshift, with the highest redshift bin peaking around z ∼ 1.7, but at such a high redshift the multipole ℓ_max = 1500 projects only to k ∼ 0.5 h Mpc⁻¹, such that the information gain on the CWDM model is marginal. We note that the arguments presented in this paragraph are only valid in the case of the Euclid pessimistic case and for the CWDM model, which has the same power spectrum as ΛCDM up to some large wavenumber k (given by the free-streaming scale). For instance, we shall see that for the 1b-DDM model Euclid is much more constraining than KiDS even with pessimistic assumption because, in that case, the power spectrum contains information on the DM parameters on larger scales/smaller wavenumbers.

A second reason is that, in our analysis, we marginalise over six baryonic feedback parameters, including the three ν_ℬ parameters accounting for a drift of the main feedback parameters ℬ with redshift. In the KiDS analysis, the three ℬ parameters are instead assumed to be redshift-independent. However, the effect of WDM on the matter power spectrum is redshift-dependent (it decreases when the redshift decreases due to mode-mode coupling). Thus, in our analysis, it is easier to cancel the effect of WDM with a shift in the baryonic feedback parameter, and there is more degeneracy between the WDM and baryonic parameters. This tends to lower our forecasted sensitivity and to compensate the fact that Euclid measurements of the lensing power spectrum are expected to be much more accurate. In Fig. 11, we choose to present the results of the Euclid pessimistic forecast with the same treatment as in KiDS analysis of Hervas-Peters et al. (2024); that is, with only three free baryonic feedback parameters, while in the Euclid optimistic forecast we stick to our more conservative baseline treatment with a marginalisation over six parameters ℬ and ν_ℬ.

Fig. 12. Same as Fig. 9, but for the 1b-DDM model, parameterised by the DDM fraction $f_{\mathrm{ddm}}^{\mathrm{ini}}$ and decay rate Γddm. The forecast assumes flat priors on ($f_{\mathrm{ddm}}^{\mathrm{ini}}$, ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$) because the effect of 1b-DDM on the linear matter power spectrum scales with the product ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$ (see Sect. 2.2). The model is equivalent to pure ΛCDM in the small ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$ limit. The shaded grey area restricts the parameter space to the region where τddm = 1/Γddm ≥ 31.6 Gyr in which the emulator was trained.

In the optimistic case, we still find that the sensitivity of the Euclid WL probe is approximately three times bigger than that of KiDS, despite of our more conservative modelling of baryonic feedback. In addition, we have already seen that with the inclusion of all the information from the 3 × 2pt probe, we can gain a factor three in sensitivity. All in all, under optimistic assumptions, Euclid could be about ten times more constraining, while providing at the same time some bounds that will be more robust against baryonic feedback effects. Furthermore, we note that the results of our optimistic forecast confirm previous findings from Schneider et al. (2020) using more realistic assumptions for the survey characteristics, especially regarding the tomographic galaxy binning.

5.2. Dark matter with one-body decay

Main results. For the 1b-DDM model, we perform forecasts using the parameter combinations $f_{\mathrm{ddm}}^{\mathrm{ini}}$ and ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$ defined in Sect. 2.2. The fiducial model is chosen to be a pure ΛCDM model. We showed in Sect. 3.2 that the small-scale suppression induced by 1b-DDM on the linear power spectrum depends only on the product ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$, while the non-linear evolution adds a bit of sensitivity to Γ_ddm alone. Thus, the data are expected to provide bounds mainly on the product of the two DDM parameters. Our fiducial values and priors are summarised in Table 4 for these parameters and Table 2 for all other free parameters. The fiducial model is chosen to be a pure ΛCDM model. As was already stated in Sect. 2.2, we use linear priors on $f_{\mathrm{ddm}}^{\mathrm{ini}}$ and ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$. We note that the N-body simulations used by Hubert et al. (2021) to build the emulator were limited to models with a DM lifetime much larger than the age of the Universe, namely τ_ddm ≥ 31.6 Gyr. We thus conservatively include in our run a prior Γ_ddm ≤ 0.0316 Gyr⁻¹. This additional prior excludes a small triangle in the parameter space defined by the priors of Table 4.

Figure 12 presents the 95% confidence level (CL) isocontours of the marginalised posterior for the 1b-DDM parameters inferred from our forecasts for the Euclid WL probe alone (WL), the full Euclid photometric probe (3 × 2pt), and Euclid 3 × 2pt combined with Planck, considering both pessimistic (left panel) or optimistic (right panel) assumptions. The fiducial model, ΛCDM, spans the lower horizontal axis (${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}=0$). The shaded grey area restricts the parameter space to the region in which the emulator was trained. The fact that the contour edges remain nearly horizontal is consistent with the fact that 1b-DDM effects depend mainly on the product ${\mathrm{\Gamma }}_{\mathrm{ddm}}{f}_{\mathrm{ddm}}^{\mathrm{ini}}$. The small tilting of the contours comes from the fact that non-linear corrections to the 1b-DDM effects do depend on Γ_ddm alone. For each of the six cases shown in Fig. 12, we ran 96 chains summing up to ∼1.6 MS in each optimistic case or ∼3.8 MS in each pessimistic case. The Gelman-Rubin convergence criterium reached about |R − 1|∼0.01 for most parameters, with a worse value of 0.05 for a few parameters.