© The Authors 2026

1. Introduction

Over the past decade, line-intensity mapping (LIM) has gained attention as a promising strategy to probe both the large-scale distribution of matter and the galactic properties across cosmic time (see Kovetz et al. 2017; Bernal & Kovetz 2022, for a general overview). While the 21 cm line remains the most extensively studied tracer (see e.g. Morales & Wyithe 2010; Pritchard & Loeb 2012; Liu & Shaw 2020), the scope of LIM studies has rapidly broadened to include several other emission lines, with multiple pathfinder surveys now operating or planned across a wide range of wavelengths¹.

Alongside this experimental progress, theoretical work has explored how LIM observables relate to the physics of reionisation and astrophysical properties (see e.g. Lidz et al. 2009; Gong et al. 2012; Visbal et al. 2015; Comaschi & Ferrara 2016; Breysse & Rahman 2017; Breysse & Alexandroff 2019; Sato-Polito et al. 2020; Zhou et al. 2021). For instance, the amplitude and shape of the LIM power spectrum (PS) encode information about the abundance, clustering, and luminosity distribution of the emitting galaxies. Therefore, the PS itself, in combination with other statistics, provides a means to constrain both the moments and parametric form of the luminosity function (LF), thus offering a new route to quantify galaxy properties in the early Universe (see e.g. Li et al. 2016; Ihle et al. 2019; Marcuzzo et al. 2025).

Beyond astrophysical applications, LIM is increasingly being recognised as a tool with significant potential for cosmological tests. Because the large-scale PS traces the underlying matter distribution – modulated by biasing, redshift-space distortions (RSDs), and emission properties – it can be used to extract cosmological information (see e.g. Karkare & Bird 2018; Moradinezhad Dizgah et al. 2019; Bernal et al. 2019; Moradinezhad Dizgah & Keating 2019; Muñoz et al. 2020; Bauer et al. 2021; Bernal et al. 2021; Moradinezhad Dizgah et al. 2022a,b, 2024; Scott et al. 2023; Adi et al. 2023; Fronenberg et al. 2024; Shmueli et al. 2025), similar to galaxy clustering or the Lyman-α forest, but it extends to fainter populations and less explored epochs. This makes LIM particularly appealing for probing non-standard cosmological scenarios (Karkare et al. 2022), despite the additional challenges posed by line interlopers, bright continuum foregrounds, and uncertainties in modelling the signal that remain more severe than in mature observational probes (see e.g. Schaan & White 2021a,b, and references therein).

Continuing the endeavour of the community in identifying possible applications of LIM in cosmological contexts, we investigate its potential to shed new light on the nature of dark matter (DM), one of the biggest unsolved problems in physics (see Cirelli et al. 2024 for a recent review). In the standard model of cosmology, DM is assumed to be ‘cold’, meaning that the free-streaming length associated with its thermal velocity in the matter-dominated era is negligibly small compared to the linear size of the perturbations that lead to the formation of dwarf galaxies. Popular particle candidates for cold DM (CDM) include axions (Weinberg 1978; Preskill et al. 1983) – which are non-thermally produced in the early Universe – and thermal relics (e.g. Hochberg et al. 2014; Arcadi et al. 2018), which decouple from thermal equilibrium with standard-model particles and freeze out during the radiation-dominated era. CDM relics are non-relativistic when they decouple from the primordial plasma and are thus trapped within the primordial density perturbations produced during inflation.

Over time, some tensions have emerged between observations on sub-galactic scales and the predictions of the CDM scenario based on N-body simulations (e.g. Bullock & Boylan-Kolchin 2017, and references therein). Although baryonic physics could mitigate the discrepancies, it is possible to devise a scenario in which the DM particles are relativistic at decoupling and become non-relativistic during the radiation-dominated era. This causes the particles to freely stream out of kiloparsec-scale potential wells, giving rise to a cutoff in the linear matter PS relative to CDM (e.g. Bond & Szalay 1983) and resulting in a deficit of small haloes and delayed galaxy formation in low-mass systems – features that can be probed observationally. This ‘warm’ DM (WDM; e.g. Bode et al. 2001; Viel et al. 2005) scenario can be tested against observations using the flux PS of the Lyman-α forest (Viel et al. 2013; Iršič et al. 2017b, 2024; Villasenor et al. 2023), fluctuations in the density of cold stellar streams (Banik et al. 2021), dwarf-galaxy counts (Nadler et al. 2021), the gravitational lensing of distant quasars (Gilman et al. 2020), and combinations of these strategies (Zelko et al. 2022). Current data put a lower bound on the mass of the WDM thermal relics of a few kiloelectronvolts – at a 95% credibility level (CL).

A probe that is directly sensitive to the abundance of low-mass haloes is therefore particularly valuable since WDM primarily modifies structure formation below a well-defined mass scale. LIM is a natural candidate for these studies. Previous work has explored the impact of WDM on the LIM signal of the 21 cm line during (Sitwell et al. 2014) and after (Carucci et al. 2015) reionisation. Somewhat counterintuitively, Carucci et al. (2015) found enhanced 21 cm fluctuations in WDM relative to CDM because the depletion of low-mass haloes, combined with observational constraints on the cosmic abundance of neutral hydrogen, shifts more H I into massive, highly biased haloes. However, they also concluded that SKA1-LOW would be unable to distinguish WDM from CDM at z = 5 without several thousand hours of integration. This limitation reflects astrophysical considerations: In the post-reionisation Universe, diffuse gas in haloes with M ≲ 10⁹ M_⊙ is expected to photoevaporate rapidly due to insufficient self-shielding against the ultraviolet background (e.g. Rees 1986; Barkana & Loeb 1999). For WDM particle masses consistent with current limits, the associated cutoff in the halo mass function (HMF) lies close to this evaporation scale, implying that little 21 cm emission originates from the low-mass CDM haloes absent in WDM.

These considerations motivate the exploration of tracers of star-formation activity, such as [C II], CO rotational transitions, and Balmer recombination emission. In particular, [C II] and CO trace cold, dense gas associated with the inner star-forming regions of galaxies – gas that is far less susceptible to photoevaporation and therefore provides a cleaner probe of the low-mass structures where WDM and CDM predictions diverge most. This expectation is supported by both theoretical and observational studies showing that reionisation primarily removes diffuse, low-density gas from dwarf-galaxy haloes while leaving the dense, self-shielded component largely intact. As a result, star formation in low-mass systems is expected to taper gradually rather than cease abruptly (Susa & Umemura 2004; Hoeft et al. 2006; Okamoto & Frenk 2009; Sawala et al. 2012; Simpson et al. 2013). Recent simulations further reveal that the survival and re-ignition of star formation depend sensitively on the retention of cold, dense gas and on halo concentration, with some low-mass galaxies capable of initiating or renewing star formation as late as z ≃ 2 in favourable conditions (Moreno et al. 2026).

Observational evidence echoes this diversity. Nearby dwarfs display a wide range of star-formation histories shaped by reionisation, including systems that quenched early, systems that experienced delayed ignition, and ‘two-component’ galaxies where old and young populations overlap (Benítez-Llambay et al. 2015). Deep JWST imaging of isolated low-mass galaxies, such as Leo P, reveals early star formation followed by a long pause after reionisation and a subsequent reignition at later times (McQuinn et al. 2024). Complementary HST and VLA studies of ultra-faint dwarfs near M31 show similarly varied behaviour: Some quenched shortly after reionisation, while others formed stars well into the post-reionisation era despite extremely low present-day gas content (Jones et al. 2026). Even in isolated systems, spatially resolved JWST data reveal radial and azimuthal variations in stellar age gradients consistent with localised star-formation episodes triggered by the interaction of dense gas with the surrounding medium (Cohen et al. 2025).

Taken together, these results emphasise that dense star-forming gas can persist in low-mass haloes long after reionisation, making tracers such as [C II] ideal for probing the small-scale structure most affected by WDM. In this work, we investigate whether future ground-based LIM experiments targeting the [C II] fine-structure line can provide competitive constraints on the WDM particle mass. The [C II] transition arises mainly in photon-dominated regions at the boundaries of molecular clouds exposed to far-ultraviolet radiation, making it an excellent probe of ongoing star formation. Although previous estimates suggest that the [C II] LIM signal at z ≃ 5 is dominated by haloes of 10¹¹–10¹² M_⊙ (Yue et al. 2015), contributions from less massive galaxies may remain non-negligible, particularly at lower redshifts. This makes the [C II] LIM PS a promising observable for constraining the abundance of the low-mass haloes suppressed in WDM cosmologies.

To explore the potential of this approach, we adopted a representative observational setup based on the Deep Spectroscopic Survey (DSS), to be carried out with the 6 m Fred Young Submillimeter Telescope (FYST) (CCAT-Prime Collaboration 2023). We selected this ground-based survey, already scheduled for deployment, as it offers a practical near-term route to carrying out LIM observations. Thanks to its high and dry location, which enables efficient sub-millimetre measurements, combined with a wide field of view and fast mapping speed, this setup provides both the sensitivity and the survey volume required for a robust PS measurement. In addition, the DSS is planned to target two fields that overlap with the Extended Cosmic Evolution Survey (E-COSMOS; Aihara et al. 2018; see also Scoville et al. 2007), the Extended Chandra Deep Field South (E-CDFS; Lehmer et al. 2005), and the Hubble Ultra Deep Field (H-UDF; Beckwith et al. 2006) survey areas, providing an almost unique opportunity to obtain a comprehensive multi-wavelength characterisation of the same underlying galaxy population. We assess how survey design choices – such as sky coverage, instrumental sensitivity, and spectral resolution – affect the constraining power of the [C II] LIM signal. Focusing on a redshift window centred at z ≃ 3.6, we evaluate the extent to which the [C II] PS can reveal the small-scale suppression that characterises WDM models, offering a novel and complementary path to probing the nature of DM beyond the CDM paradigm.

To enable comparison with other results from the literature, we also express our results in terms of fuzzy DM (FDM; e.g. Iršič et al. 2017a; Lazare et al. 2024), an alternative model in which DM consists of ultra-light bosons that inhibit structure formation on small scales due to quantum effects. This conversion is possible because both WDM and FDM produce a characteristic cutoff in the linear matter PS, albeit for different physical reasons. By matching the half-mode wavenumbers (i.e. the scale at which this suppression occurs), we can translate a constraint on the WDM particle mass into an effective constraint on the FDM boson mass. We note, however, that this correspondence is approximate and model dependent.

The paper is organised as follows. Sect. 2 introduces WDM and FDM scenarios, detailing how the beyond-CDM PS is derived from its CDM counterpart. In Sect. 3 we briefly summarise the halo-model formalism used to compute the LIM PS. In Sect. 4 we consider the high-redshift [C II] LF and describe the abundance-matching procedure used to connect halo masses and [C II]-line luminosities. Sect. 5 introduces the different survey setups considered throughout this work and presents the predicted [C II] PS with relative uncertainties. In Sect. 6 we outline our Bayesian inference pipeline applied to mock data and report the resulting constraints on the WDM particle mass. We show their conversion to FDM as well. In Sect. 7, we discuss critical aspects of our modelling and inference strategy, with particular attention given to the relative contribution of different halo masses to the LIM signal and, consequently, to the impact of varying the faint-end slope of the [C II] LF. Finally, in Sect. 8 we provide a summary of our findings and draw the main conclusions. The possibility of including other redshift bins and/or emission lines is also taken into account.

We assumed a flat Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology, adopting h = 0.674 for the dimensionless Hubble parameter. The present-day density parameters are set to Ω_m = 0.315 for matter, Ω_b = 0.049 for baryons, and Ω_Λ = 0.685 for the cosmological constant. The primordial PS of density fluctuations is described by a spectral index n_s = 0.965 and a normalisation σ₈ = 0.811. The linear matter PS in the CDM scenario, P_CDM, is obtained using the Code for Anisotropies in the Microwave Background (CAMB;² Lewis et al. 2000).

2. Cold, warm, and fuzzy dark matter

Fermionic DM of mass m_WDM∼ keV/c² produced in thermal equilibrium results in suppressed small-scale clustering with respect to CDM. This is because the DM particles are relativistic at production and free stream out of primordial perturbations within a mass-dependent characteristic scale. We parameterised the damping in the linear matter PS as

$$\begin{array}{c}\hfill P_{\mathrm{WDM}}(k)=T^2(k)P_{\mathrm{CDM}}(k),\end{array}$$(1)

where (e.g. Viel et al. 2005)

$$\begin{array}{c}\hfill T(k)={[1+{\left(\tau k\right)}^{2\nu }]}^{-5/\nu },\end{array}$$(2)

with ν = 1.12 and

$$\begin{array}{c}\hfill \tau =0.049{\left(\frac{m_{\mathrm{WDM}}}{1\mathrm{keV}}\right)}^{-1.11}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{0.25}\right)}^{0.11}{\left(\frac{h}{0.7}\right)}^{1.22}{h}^{-1}\mathrm{Mpc}.\end{array}$$(3)

This was obtained by fitting the output of Boltzmann codes for wavenumbers k < 5 h Mpc⁻¹.

For both CDM and WDM models, we computed the HMF and linear bias by using the fits to numerical simulations (Sheth & Tormen 1999; Sheth et al. 2001). Following Parimbelli et al. (2021), we determined the linear mass variance smoothed on the halo scale with the so-called smooth-k window function, which avoids the excess of low-mass haloes predicted by the standard top-hat filter in models with a small-scale cutoff while providing a reasonable fit to N-body simulation results:

$$\begin{array}{c}\hfill \stackrel{\sim }{W}(k)={[1+{(k{R}_{\mathrm{s}})}^{{\beta }_{\mathrm{s}}}]}^{-1},\end{array}$$(4)

with β_s = 4.8 and $R_{\mathrm{s}}=3.3^{-1}\sqrt[3]{3M/(4\pi {\overline{\rho }}_{\mathrm{m}})}$, where M denotes the halo mass and ${\overline{\rho }}_{\mathrm{m}}$ is the mean comoving matter density. For completeness, we include in Appendix A a graphical illustration of the effects of WDM relative to CDM on the PS, HMF, and halo bias.

We also considered the possibility that DM is made of ultra-light bosons with kiloparsec-scale De Broglie wavelengths, the FDM. In this case, we determined the mass of the FDM particles, m_FDM, through the relation (Hu et al. 2000)

$$\begin{array}{c}\hfill k_{0.5}=4.5{\left(\frac{m_{\mathrm{FDM}}{c}^2}{{10}^{-22}\mathrm{eV}}\right)}^{4/9}{\mathrm{Mpc}}^{-1},\end{array}$$(5)

with T²(k_0.5) = 0.5.

3. The theoretical [C II] power spectrum

We describe the large-scale [C II] intensity fluctuations using the halo-model formalism, which provides a physically motivated framework for linking the statistical properties of matter to those of the underlying halo population (Cooray & Sheth 2002). This approach has been extensively applied to LIM, where the line emission is connected to haloes through empirically or theoretically motivated luminosity–mass relations (e.g. Silva et al. 2015; Gong et al. 2017). A summary of the main aspects of our implementation is presented below, while further details and related prescriptions can be found in Marcuzzo et al. (2025, M25 hereafter).

Under the assumption that line emission occurs in DM haloes and neglecting absorption and scattering, the mean specific intensity observed at redshift z for a transition with rest-frame frequency ν_rf is given by

$$\begin{array}{c}\hfill {\overline{I}}_{\nu }(z)=\frac{c}{4\pi H(z){\nu }_{\mathrm{rf}}}{\overline{\rho }}_{\mathrm{L}}(z),\end{array}$$(6)

where ${\overline{\rho }}_{\mathrm{L}}(z)$ is the mean comoving luminosity density, i.e. the first moment of the LF.

To relate galaxy luminosities to halo properties, we introduce the conditional luminosity function (CLF), ϕ(L|M, z), which quantifies the average number of galaxies with luminosity L inside haloes of mass M at redshift z. This allows one to express the LF as

$$\begin{array}{c}\hfill \mathrm{\Phi }(L,z)={\displaystyle {\int }_0^{\infty }}\varphi (L|M,z)\frac{\mathrm{d}{\overline{n}}_{\mathrm{h}}}{\mathrm{d}M}(M,z)\mathrm{d}M.\end{array}$$(7)

We further define the n-th moments of the CLF as

$$\begin{array}{c}\hfill {\eta }_n(M,z)={\displaystyle {\int }_0^{\infty }}{L}^n\varphi (L|M,z)\mathrm{d}L,\end{array}$$(8)

which correspond to the total Lⁿ-weighted luminosity in haloes of mass M. For instance, η₁ gives the mean total [C II] luminosity at a fixed halo mass, and thus the mean emissivity becomes

$$\begin{array}{c}\hfill {\overline{\rho }}_L(z)={\displaystyle {\int }_0^{\infty }}{\eta }_1(M,z)\frac{\mathrm{d}{\overline{n}}_{\mathrm{h}}}{\mathrm{d}M}(M,z)\mathrm{d}M.\end{array}$$(9)

These moments capture the distribution of line luminosities within haloes of varying mass and play a central role in the following determination of the LIM PS model.

Under the assumptions of linear biasing, Poisson statistics, and the distant-observer approximation, the redshift-space PS of intensity fluctuations can be decomposed as

$$\begin{array}{c}\hfill P(k,\mu ,z)=P_{\mathrm{clust}}(k,\mu ,z)+P_{\mathrm{shot}}(z),\end{array}$$(10)

where P_clust describes the large-scale clustering and P_shot accounts for the shot noise due to the discrete nature of sources. Moreover, $\mu =\widehat{\mathit{k}}·\widehat{\mathit{n}}$ is the cosine of the angle between the wavevector and the line of sight.

On large scales, the clustering term is sourced by biased, linearly evolving tracers of the matter field and can be modelled as follows:

$$\begin{array}{c}\hfill P_{\mathrm{clust}}(k,\mu ,z)={\overline{I}}_{\nu }^2(z){[b(z)+f(z){\mu }^2]}^2\mathcal{D}(k,\mu ,z)P_{\mathrm{m}}(k,z),\end{array}$$(11)

where b(z) denotes the luminosity-weighted linear bias, f(z) the linear growth rate, and P_m(k, z) the linear matter PS, taken as P_CDM or P_WDM depending on the underlying cosmological model. The damping function 𝒟, which models RSDs, can be parameterised using various functional forms. In this work, we adopt the squared Lorentzian form,

$$\begin{array}{c}\hfill \mathcal{D}(k,\mu ,z)={[1+\frac{{(k\mu \sigma )}^2}{2}]}^{-2},\end{array}$$(12)

although our conclusions do not significantly change if a different commonly used shape is employed. The parameter σ indicates a characteristic comoving displacement that encodes the strength of the damping generated by the RSDs, which is expected to be roughly comparable to the pairwise velocity dispersion scaled by (aH)⁻¹. Finally, the luminosity-weighted bias, b(z), can be expressed in terms of ${\overline{\rho }}_{\mathrm{L}}$ as

$$\begin{array}{c}\hfill b(z)=\frac{1}{{\overline{\rho }}_L(z)}{\displaystyle {\int }_0^{\infty }}{\eta }_1(M,z)b_{\mathrm{h}}(M,z)\frac{\mathrm{d}{\overline{n}}_{\mathrm{h}}}{\mathrm{d}M}(M,z)\mathrm{d}M.\end{array}$$(13)

On the other hand, the shot-noise component is given by

$$\begin{array}{c}\hfill P_{\mathrm{shot}}(z)=\frac{{\overline{I}}_{\nu }^2(z)}{{\overline{n}}_{\mathrm{eff}}(z)},\end{array}$$(14)

where ${\overline{n}}_{\mathrm{eff}}$ is the effective number density of emitters and can be defined as

$$\begin{array}{c}\hfill {\overline{n}}_{\mathrm{eff}}(z)={\left[\frac{1}{{\overline{\rho }}_L^2(z)}{\displaystyle {\int }_0^{\infty }}{\eta }_2(M,z)\frac{d{\overline{n}}_{\mathrm{h}}}{\mathit{dM}}(M,z)dM\right]}^{-1}.\end{array}$$(15)

4. Mass–luminosity relation for [C II] emitters

Building on the halo–galaxy connection for [C II] emitters identified in the MARIGOLD simulations (Khatri et al. 2025), and following the methodology of M25, we assigned a [C II] luminosity, L, to each halo via abundance matching, assuming L increases monotonically with halo mass. This procedure implicitly assumes that [C II] emission is dominated by a single source per halo and neglects scatter in luminosity at fixed mass.

To implement this procedure, one requires the LF of [C II] emitters at the relevant redshift. Despite recent observational progress, this function remains highly uncertain owing to the limited number of detections currently available. We therefore relied on model predictions calibrated against existing measurements to construct a plausible LF, which serves as the basis for the abundance-matching relation.

Fig. 1. [C II] luminosity–halo mass relation inferred from abundance matching. In both panels, the solid lines correspond to the L(M) relations derived from our Schechter fits to the ALPINE data with fixed values for α. The dashed red line is based on the fit by Yan et al. (2020, Y20). The left and right panels assume a CDM and a 0.5  keV HMF, respectively (see Appendix A for a direct comparison between the two scenarios). In the left panel, for comparison, the dot-dashed black line indicates the L(M) relation obtained by Silva et al. (2015, S15), while the dotted gold and dark gold lines show results from the MARIGOLD simulations (Khatri et al. 2025, K25) at z = 5 and z = 4, respectively.

We modelled the [C II] LF using the Schechter form

$$\begin{array}{c}\hfill \mathrm{\Psi }(L)\equiv \frac{\mathrm{d}n}{\mathrm{d}{log}_{10}L}={\mathrm{\Psi }}_{\ast }{\left(\frac{L}{L_{\ast }}\right)}^{1+\alpha }exp(-\frac{L}{L_{\ast }}),\end{array}$$(16)

where Ψ_*, L_*, and α denote the normalisation, characteristic luminosity, and faint-end slope, respectively. To bracket the current observational uncertainties, we adopted two parallel LF parameterisations (see Fig. 1 in M25). The first is a high-normalisation, or ‘optimistic’, model based on the multi-wavelength and multi-redshift fit of Yan et al. (2020). The second is a low-normalisation, or ‘pessimistic’³, model obtained by fitting only the targeted sample from the ALPINE [C II] survey (Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020). For the optimistic and pessimistic cases, we used (Ψ_*, L_*, α) = (10^−3.08 Mpc⁻³ dex⁻¹, 10^9.5 L_⊙, −1.1) and (10^−3.02 Mpc⁻³ dex⁻¹, 10^8.73 L_⊙, −1.1), respectively. These two parameterisations define the fiducial models analysed throughout this work. To assess the sensitivity of our results to the assumed faint-end shape, we additionally performed a series of fits to the ALPINE data with fixed slopes in the range α ∈ [ − 1.9, −0.5].

In our abundance-matching procedure, the observed [C II] LF was used as a constraint that both cosmological scenarios must satisfy. Since the LF is an empirical property of the galaxy population, its bright end must be reproduced in both CDM and WDM. Only the faint end is allowed to differ between the two models, as it is not yet well constrained observationally and naturally reflects the differing halo abundances predicted by each cosmology. The resulting abundance-matching relations are shown in Fig. 1 for the two DM scenarios (left: CDM; right: WDM). The WDM case adopts an extreme particle mass of 0.5 keV to maximise the contrast with CDM. Unless the faint-end slope of the LF is very shallow, the scarcity of low-mass haloes in WDM implies the existence of a minimum [C II] luminosity: There are simply not enough haloes to host galaxies below this threshold. As a result, the abundance matching cannot be fully satisfied, and the input LF must be truncated in the WDM case (see Fig. 2). This minimum luminosity corresponds to the horizontal asymptotes of the L(M) relation at low halo masses. This feature partly reflects our assumption of one source per halo, but it arises primarily because the WDM HMF turns over and rapidly declines at low masses (see Appendix A). The resulting plateau in the L(M) relation is therefore an artefact of the abundance-matching procedure, signalling the drop in halo abundance rather than a physical property of the luminosity–mass relation. In the halo model of Sect. 3, the effective PS components depend on the product $L(M)\mathrm{d}{\overline{n}}_{\mathrm{h}}/\mathrm{d}M$, so the fall-off of the HMF further suppresses the contribution from these low-mass systems. Finally, while combining L(M) with the CDM HMF reproduces the LF across all luminosities, in WDM scenarios this reconstruction is only possible down to the minimum luminosity associated with the L(M) plateau. Below this point the predicted abundance drops sharply, as shown in Fig. 2, reflecting the underlying paucity of low-mass haloes.

Fig. 2. Reconstructed [C II] LFs obtained by combining the abundance-matched L(M) relations with the HMFs of the CDM and WDM cosmologies. The figure uses the LF of our optimistic case and adopts a WDM particle mass of 0.5 keV. While the CDM model reproduces the input LF across all luminosities, the WDM case shows a sharp downturn at the faint end, reflecting the dearth of low-mass haloes and the corresponding minimum luminosity imposed by the abundance-matching procedure.

Fig. 3. Predicted [C II] PS. Left: Δ2(k, z ≃ 3.6) computed in the CDM scenario for our optimistic (solid red line) and pessimistic (solid blue line) cases for α = −1.1. The shaded areas represent the associated statistical uncertainty as for our survey reference setup. The dotted line is a graphical representation of the white-noise level, PWN, while the dashed and dot-dashed lines refer to the clustering and shot-noise components, respectively. Right: ratio between the WDM and CDM power spectra. The dashed and dotted black lines correspond to WDM models with mWDM = 10 keV and 0.5 keV, respectively. The shaded regions reflect the CDM uncertainties from the left panel. In both panels, the markers illustrate the adopted binning scheme.

5. The ‘observed’ [C II] power spectrum

In this section, we present predictions for the [C II] LIM PS that will be observed with the FYST DSS at z ≃ 3.6 and explore how it is affected by changes in the nature of DM. All technical aspects concerning observational setup, instrumental effects, foreground removal, and PS error estimation follow the methodology outlined in M25; here we only summarise the key elements and refer the reader to that work for further details. In our baseline configuration, a 16 deg² patch of sky is mapped with a total integration time of 8000 hours, a spectral resolution of R = 100, and thus a resulting white-noise level of P_WN ≃ 2.4 × 10¹⁰ h⁻³ Mpc³ Jy² sr⁻². In the remainder of this paper, we refer to this configuration as ‘reference setup’.

To investigate the dependence of the results on survey design, we considered variations in sky coverage, sensitivity and spectral resolution (R). Specifically, we analysed configurations with areas of 160 deg², 1600 deg², and half of the sky. In parallel, we examined a scenario where the instrumental sensitivity is enhanced by a factor of $\sqrt{10}$, leading to a tenfold reduction in P_WN. Finally, we assessed the performance of a futuristic instrument with improved resolution, increasing R to 500 and rescaling the total survey time accordingly.

To obtain a realistic estimate of the observed [C II] PS, we started from its theoretical definition given in Eq. (10) and related equations, incorporating the instrumental effects through the following expression:

$$\begin{array}{c}\hfill P_{\mathrm{obs}}(k,\mu ,z)=P(k,\mu ,z)W_{\perp }(k,\mu )W_{\Vert }(k,\mu ),\end{array}$$(17)

where W_⊥(k, μ) accounts for the Gaussian smoothing due to the telescope beam, and W_∥(k, μ) describes the Lorentzian damping introduced by the finite spectral resolution. We then defined the direction-averaged PS (i.e. the monopole moment) as

$$\begin{array}{c}\hfill P_0(k,z)=\frac{{\displaystyle {\int }_{{\mu }_{\min }}^{{\mu }_{\max }}}{P}_{\mathrm{obs}}(k,\mu ,z)\mathrm{d}\mu }{{\displaystyle {\int }_{{\mu }_{\min }}^{{\mu }_{\max }}}\mathrm{d}\mu },\end{array}$$(18)

where μ_min = k_f^∥/k and μ_max = min(1, k_max^∥/k), with k_f^∥ and k_max^∥ being the fundamental and maximum wavenumbers along the line of sight, respectively. The integration limits were selected from the full range [0, 1] to approximately mitigate continuum foregrounds⁴ and to account for the effects of finite spectral resolution.

The associated uncertainty is given by

$$\begin{array}{c}\hfill {\sigma }_{P_0}(k)=\frac{P_0(k)+P_{\mathrm{WN}}}{\sqrt{N_{\mathrm{m}}(k)}},\end{array}$$(19)

where N_m denotes the number of independent Fourier modes contributing to each k-bin and can be defined as

$$\begin{array}{c}\hfill N_{\mathrm{m}}(k)=\frac{\min (k,k_{\max }^{\Vert })k\mathrm{\Delta }k{V}_{\mathrm{surv}}}{4{\pi }^2},\end{array}$$(20)

with V_surv being the survey volume. For practical purposes, we evaluated P₀ and σ_P0 in discrete k-bins with fixed width Δk. Throughout the analysis, we adopted Δk = 5 k_f^∥, but we have verified that our results are not sensitive to this specific choice.

Fig. 3 (left panel) displays the predicted [C II] PS for the reference setup. We plot the quantity Δ², which is connected to P₀ in Eq. (18) via the following relation:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^2(k,z)=\frac{k^3}{2{\pi }^2}{P}_0(k,z),\end{array}$$(21)

Solid lines represent the theoretical model, while shaded regions denote the statistical uncertainty derived from Eq. (19). Red and blue curves correspond to our optimistic and pessimistic LF assumptions, respectively, both adopting a faint-end slope of α = −1.1. The discrete points overlaid on each curve indicate the binned measurements resulting from the actual k-binning procedure.

The right panel of Fig. 3 shows the ratio of the WDM and CDM power spectra, focusing on the two most extreme WDM scenarios considered in our Bayesian analysis (see Sect. 6). Specifically, these examples were selected to show how models associated with very different WDM particle masses compare to the CDM prediction. We underline that the higher signal amplitude found in WDM scenarios results from an L(M) relation that exceeds the CDM counterpart at halo masses around the HMF cutoff, where the reduced abundance of WDM haloes is compensated by their higher luminosities. This follows from the requirement that both models reproduce the observed bright end of the [C II] LF, thereby raising the WDM luminosities at intermediate halo masses (i.e. 10¹¹ − 10¹² h⁻¹ M_⊙). Because these masses also provide the bulk of the [C II] PS contribution (see Fig. 8), this compensation directly affects the overall amplitude. As already mentioned in the Introduction, an analogous feature has also been reported in H I studies (Carucci et al. 2015). This behaviour is further examined in Sect. 7.1, where we discuss the modelling choices and assumptions that lead to this outcome.

Fig. 8. Differential contribution of different halo masses to ${\overline{I}}_{\nu }$, b, and ${\overline{n}}_{\mathrm{eff}}^{-1}$ (from top to bottom) based on our halo model for CDM at z = 3.6.

6. Bayesian inference

We investigate the potential of [C II] PS measurements to constrain the nature of DM through a Bayesian analysis on mock observations. The focus is on two key parameters: the WDM particle mass, m_WDM, and the line-of-sight displacement, σ, which partially encapsulates the effect of RSDs. The statistical framework follows standard practice and employs a Gaussian likelihood function:

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{\theta }|\mathit{D})\propto exp\{-\frac{1}{2}{\displaystyle \sum _i}\frac{{[D_i-M_i(\mathit{\theta })]}^2}{{\sigma }_i^2}\},\end{array}$$(22)

with θ ≡ {σ, m_WDM} being the model parameters, D ≡ {D_i} the mock data⁵ (i.e. the [C II] PS monopole evaluated in k_i-bins), M_i the corresponding model predictions for a given parameter set, and σ_i the statistical uncertainties, when the covariance matrix is assumed to be diagonal.

In Section 6.1, we assume CDM as the true underlying model and assess which lower bounds on m_WDM can be placed at a 95% CL by fitting the mock data with WDM-based models, thus ruling out excessively warm scenarios. In Section 6.2, we consider data generated with m_WDM = 3 keV and examine whether the true mass can be successfully recovered with both upper and lower limits.

When producing the mock data, a fiducial value of σ = 3 h⁻¹ Mpc was adopted to describe the damping caused by RSDs along the line of sight, following the rationale provided in Sect. 4.6 of M25. Throughout the analysis, the LF parameters were assumed to be precisely known from independent observations. Results are presented for both our optimistic and pessimistic LFs with a faint-end slope α = −1.1. We assess how the constraints vary under improvements in survey area, sensitivity, and spectral resolution as discussed in Sect. 5.

6.1. CDM-generated data

In this section, we consider mock data generated under a CDM cosmology, aiming to determine a lower limit for the WDM particle mass across different survey configurations and LF assumptions.

A flat prior was imposed on the displacement parameter σ within the interval [0.2, 5] h⁻¹ Mpc. The prior range for σ spans values typical of central galaxies in low-mass haloes up to those representative of satellite systems, including contributions from internal gas dynamics and large-scale relative motions. As our main focus is on m_WDM, we marginalised the posterior distributions over σ in order to get the final constraints. For the WDM particle mass, m_WDM, we considered values above 0.5 keV because lower masses are strongly disfavoured by existing constraints from galaxy clustering, Lyman-α forest, and reionisation (e.g. Menci et al. 2012; Dayal et al. 2017; Iršič et al. 2017b).

Performing statistical inference on m_WDM in a CDM universe requires particular care as the CDM case formally corresponds to the limit m_WDM → +∞. In the literature, this is commonly addressed by reparameterising the problem in terms of the inverse particle mass w = (m_WDM/keV)⁻¹ (e.g. Markovic et al. 2011; Smith & Markovic 2011; Iršič et al. 2017b, 2024; Enzi et al. 2021; Rudakovskyi et al. 2021; Villasenor et al. 2023). A flat prior in w then corresponds to an informative prior in m_WDM that scales as ${\pi }_{m_{\mathrm{WDM}}}\propto m_{\mathrm{WDM}}^{-2}$, thereby assigning higher a priori probability to lower particle masses⁶. Since current observations can only provide lower limits on m_WDM, this choice can be regarded as conservative.

Several authors (e.g. Rudakovskyi et al. 2021; Villasenor et al. 2023) have examined the dependence of their results on the prior by instead adopting a uniform distribution in m_WDM, truncated at a finite maximum value $m_{\mathrm{WDM}}^{\max }\simeq 10\mathrm{keV}$. This range necessarily excludes the formal CDM limit, as including it would require an improper prior. In practice, the marginalised likelihood from current datasets is insensitive to larger m_WDM values, since they yield predictions indistinguishable from CDM within measurement uncertainties. However, when setting a lower limit on m_WDM at a fixed CL (e.g. 95%), the result inevitably depends on $m_{\mathrm{WDM}}^{\max }$ through the posterior normalisation. Consequently, the inferred lower bound can be artificially shifted by changing $m_{\mathrm{WDM}}^{\max }$ (see Appendix B for an illustrative example).

To assess the sensitivity of our analysis to prior assumptions, we explored several alternative forms. Motivated by the absence of a physically preferred parameterisation, we adopted a family of priors ${\pi }_{\mathit{w}}\propto {\mathit{w}}^{\beta }$ with a uniform hyperprior on β ∈ [ − 0.9, 0], where the upper limit β = 0 corresponds to the commonly used flat prior in w. It is worth noting that the case β → −1 approaches a flat prior in the logarithm of m_WDM/keV (γ = −1), which cannot be realised exactly because it would be improper; for this reason, we restricted the hyperprior range to values slightly above −1. Over the selected domain w ∈ (0, 2], the normalised prior distribution with a fixed value of β can be expressed as

$$\begin{array}{c}\hfill {\pi }_w=\frac{1+\beta }{2^{1+\beta }}{w}^{\beta }.\end{array}$$(23)

In the remainder of the analysis, all results are presented after marginalising over the hyperprior, which is equivalent to using

$$\begin{array}{c}\hfill {\pi }_w=\frac{(y-1)-{(2/w)}^{0.9}(-0.9y+y-1)}{1.8y^2},\end{array}$$(24)

with y = ln(w/2). This prior probability distribution has its 95^th percentile at w ≃ 1.821, corresponding to m_WDM ≃ 0.549 keV.

We also tested a (regularised) two-dimensional Jeffreys prior, proportional to the square root of the determinant of the Fisher information matrix, which ensures reparameterisation invariance but is not necessarily non-informative. In our setup, the Fisher information decreases rapidly for m_WDM ≳ 10 keV, as the predicted power spectra become increasingly insensitive to the particle mass. As a result, the Jeffreys prior suppresses the region corresponding to the CDM limit (w → 0) used to generate the mock data, effectively penalising the true model. This suppression shifts the posterior peak away from both the likelihood maximum and the fiducial parameters. Since this shift reflects the structure of the prior rather than information contained in the data, we do not show results obtained with the Jeffreys prior. For completeness, we note that the resulting lower bounds on m_WDM are weaker by a factor of 1.2–1.3, depending on the configuration considered.

Fig. 4. Likelihood ratio ℓ computed for the optimistic case. Left: ℓ in the reference setup. Centre: Same as the left panel, but with R = 500. Right: Graphical explanation of the banana-shaped ℓ shown in the left panel. Each pair of values (w, σ) in the legend indicates the inverse particle mass in keV−1 and the RSD displacement parameter in units of h−1 Mpc, respectively. The stars of corresponding colours in the left panel serve as identifiers for each combination in the parameter space, while the red circle denotes the true values used to generate the data.

Fig. 5. Marginalised posterior (over σ) computed for the whole set of sky areas in the optimistic scenario, with R = 100 and PWN ≃ 2.4 × 1010 h−3 Mpc3 Jy2 sr−2. Left: Posterior for a fixed β = 0 (uniform prior in w). Right: Result marginalised over β. The dashed lines of corresponding colours indicate the mWDM threshold up to which CDM and WDM can be distinguished (95% CL). The dotted lines show the assumed prior distribution for w.

As anticipated, we performed the analysis considering both the optimistic and pessimistic LF assumptions, each characterised by a faint-end slope of α = −1.1. In all cases, we find that WDM masses above 10 keV are virtually indistinguishable from CDM at the scales probed by LIM. The left panel of Fig. 4 illustrates a representative two-dimensional likelihood ratio, ℓ = ℒ(θ|D)/ℒ(θ_best-fit|D), in the optimistic scenario, assuming the reference survey configuration. The likelihood contours reveal that WDM models with low m_WDM values can still provide excellent fits to the CDM mock data, provided the displacement parameter σ associated with RSDs is sufficiently large. This reflects a partial degeneracy between the impact of the warm nature of DM and that of the line-of-sight velocity dispersion: Low particle masses enhance power, while large values of σ suppress it, allowing the two effects to partially cancel out. This correlation is exacerbated by the limited spectral resolution, which smooths out small-scale features and makes it difficult to disentangle the scale-dependent damping induced by RSDs from the effects of varying m_WDM. Indeed, the central panel of Fig. 4 demonstrates how the degeneracy m_WDM–σ is mitigated when adopting a higher spectral resolving power, R = 500, while keeping the same survey area and instrumental sensitivity. As discussed in M25, increasing R diminishes the effect of spectral smoothing to a level comparable with that of RSDs. This leads to a more accurate estimation of σ and significantly weakens its degeneracy with m_WDM.

The compensation mechanism introduced above is further illustrated in the right panel of Fig. 4, where we compare the fractional difference in the PS signal relative to the fiducial CDM PS (dotted red line) and related uncertainty (shaded region) for three representative cases: (i) a CDM model with σ = 5 h⁻¹ Mpc (solid grey line), corresponding to the upper bound of the prior; (ii) a 1 keV model with σ = 3 h⁻¹ Mpc (solid dark green line), matching the fiducial value for σ used to generate the mock data; and (iii) a 1 keV model with σ = 5 h⁻¹ Mpc (solid purple line). The resulting ℓ values (0.32, 0.26, and 0.94, respectively) highlight how the interplay between the two parameters of the m_WDM–σ plane can yield to good fits even in extreme WDM cosmologies.

Fig. 5 shows representative examples of the posterior distribution for w, marginalised over σ. The left panel corresponds to a prior with fixed β = 0, while the right panel uses the effective prior obtained after marginalising over the hyperparameter β. In both cases, the posterior closely follows the shape of the prior at small values of w, where the data provide essentially no constraining power. With β = 0, the posterior approaches a finite plateau as w → 0, whereas the β-marginalised prior leads to a posterior that diverges in this limit, although it remains normalisable. Importantly, this difference has no practical impact on our analysis: The high-w tail, which drives the 95% CL, is only marginally affected by the choice of prior. This also confirms that our results do not depend on the arbitrary upper bound of w = 2 (m_WDM = 0.5  keV). In this region, the likelihood strongly suppresses low masses, ensuring that the posterior constraints are dominated by the data rather than by prior assumptions.

The resulting constraints, obtained for R = 100 and marginalised over β, for all four considered survey areas, and under both sensitivity settings, are summarised in the top panel of Fig. 6. The bottom panel displays how our limits on m_WDM map onto FDM constraints via Eq. (5). For comparison, Fig. 6 also shows the current limits from Lyman-α forest analyses⁷ (Iršič et al. 2017a, 2024), which are broadly comparable to those obtained for our 1600 deg² configuration with improved sensitivity. We note, however, that the constraining power of the [C II] LIM PS can be further enhanced by incorporating multiple redshift bins and/or additional emission lines (see Sect. 8).

Fig. 6. Lower bounds on mWDM (top) and mFDM (bottom) for R = 100. The pessimistic case is represented by the blue squares, while the optimistic case is denoted by the red diamonds. Filled markers correspond to the reference sensitivity; empty markers indicate a sensitivity increased by a factor of $\sqrt{10}$ (i.e. a white-noise PS decreased by a factor of ten). Shown for reference are the dashed green and dotted orange lines, marking the lower bounds inferred by Iršič et al. (2024) and Iršič et al. (2017a), respectively.

6.2. 3 keV WDM-generated data

We generated mock data using a WDM model with m_WDM = 3 keV to determine whether this particle mass can be reliably recovered within our Bayesian framework, thereby setting both lower and upper bounds. We applied the priors on σ and m_WDM as described in Sect. 6.1.

For most survey configurations, ℓ remains significantly non-zero even for the CDM scenario (m_WDM → ∞). As a result, we cannot set meaningful upper limits at the 95% CL, as CDM remains consistent with the 3 keV mock data within the errors, in line with the findings discussed in the previous subsection. Under favourable survey setups, however, the likelihood constrains both parameters tightly, allowing us to determine both upper and lower bounds on m_WDM. An illustrative example is presented in Fig. 7, where we show the 2D ℓ and the posterior marginalised over σ as computed in the optimistic scenario for a survey with R = 500, P_WN ≃ 2.4 × 10⁹ h⁻³ Mpc³ Jy² sr⁻², and a sky coverage of 1600 deg². In this case, uniform priors on m_WDM can be safely adopted, without any additional treatment of the CDM limit. The resulting posterior and corresponding constraints at 95% CL (i.e. 2.62 < m_WDM [keV] < 4.54) are displayed as dark-green lines in the right-hand panel. For comparison, we overplot in orange those obtained by assuming a uniform prior on w = 1/m_WDM. Both approaches produce similar results, although the latter slightly favours lower particle masses as expected from the discussion in Sect. 6.1.

Fig. 7. Likelihood ratio ℓ (left) and marginalised posterior (right) obtained from 3 keV mock data in the optimistic scenario for a survey covering 1600 deg2 with R = 500 and $P_{\mathrm{WN}}\simeq 2.4\times {10}^9{h}^{-3}{{\mathrm{Mpc}}^3{\mathrm{Jy}}^2\mathrm{sr}}^{-2}$. In the right panel, the dark-green and orange vertical lines denote the 95% credible regions, computed assuming uniform priors in mWDM and $m_{\mathrm{WDM}}^{-1}$, respectively.

7. Impact of the faint-end slope

In this section, we critically assess the robustness of our results by examining how key modelling choices influence the inferred constraints on m_WDM, focusing on the CDM-generated data. We analyse the relative contribution of different halo mass scales to the LIM components, and investigate the role of the faint-end slope of the LF, exploring how variations in α and correlated Schechter parameters impact the signal and the subsequent statistical inference in the two cosmological scenarios.

7.1. Differential contributions by halo mass

In Fig. 8, we show the differential contribution of halo masses to ${\overline{I}}_{\nu }$, b, and ${\overline{n}}_{\mathrm{eff}}^{-1}$, in a CDM universe both for the optimistic LF scenario and across the range of fits obtained with α ∈ [ − 1.9, −0.5] in the pessimistic case. The clustering signal is strongly dominated by DM haloes with mass M ≃ 10¹¹–10¹² h⁻¹ M_⊙ and receives non-vanishing (but still minor) contributions from masses below 10¹⁰ h⁻¹ M_⊙ only when α < −1.4. The high-mass cutoff is provided by the HMF while the low-mass one comes from L(M). The situation is even more extreme for the effective volume per emitter which basically receives no contributions from haloes with M < 10¹¹ h⁻¹ M_⊙. It is worth noting that, although Fig. 8 refers to the CDM case, the relative weighting of halo masses remains essentially unchanged in WDM scenarios. The only difference arises in more extreme WDM models combined with a steep faint-end slope of the LF, where the scarcity of low-mass haloes leads L(M) to flatten and both the low- and high-mass cutoffs to be determined by the HMF.

As anticipated in Sect. 5, the enhanced amplitude of the PS in WDM scenarios when α = −1.1 originates from this dominance of intermediate-mass haloes in combination with our assumption of a perfectly known LF parameterisation. Specifically, since the luminosity–mass relation in our model (Fig. 1) was derived via abundance matching (Sect. 4) applied to the same LF but using either a CDM or WDM HMF, the resulting L(M) curves overlap at high masses but become systematically higher in the WDM case at masses below the scale where the WDM HMF begins to diverge from its CDM counterpart (see Appendix A). This feature, built into the model by construction, contributes to the boost in power observed for WDM relative to CDM.

All this implies that the large differences in the low-mass end of the HMF between CDM and WDM do not imprint a strong signature in the [C II] PS which is most sensitive to substantially higher halo masses. In other words, the constraining power of LIM experiments on the mass of WDM particles derives from the (smaller) variations of the HMF and linear bias parameter at M ≃ 10¹¹–10¹² h⁻¹ M_⊙. Consequently, we are unable to fully exploit the potential of [C II] as a tracer to determine the nature of DM. As noted in Sect. 1, unlike H I, which has been employed for previous studies involving LIM PS, [C II] is present in low-mass haloes that are abundant in CDM but largely absent in WDM. However, the [C II] PS is not particularly sensitive to these small haloes. Enhancing their contribution would require an ‘unconventional’ tracer whose emission is stronger in low-mass haloes and weaker in high-mass ones – potentially as a result of physical processes that inhibit its development in more massive environments.

7.2. Results for a steeper faint-end slope

The results presented in Sects. 5 and 6 refer to the assumption of α = −1.1, and therefore correspond to a scenario in which low-mass haloes (M < 10¹⁰ h⁻¹ M_⊙) are negligible contributors to the overall signal (see Sect. 7.1). As a consequence, the expected suppression of power in WDM scenarios due to the absence of this small-halo population becomes irrelevant for our purposes.

To maximise and exploit the potential contribution of [C II] emission from low-mass haloes discussed above, we now turn to steeper faint-end slopes of the LF. In particular, we repeated the analysis for different values of α, with special attention to the case α = −1.9, which amplifies the role of faint, low-mass haloes in the total emission (see Fig. 8). In this regime, the degeneracy in the m_WDM–σ parameter space discussed in Sect. 6.1 changes direction, as both parameters now act to suppress the power (left panel of Fig. 9). Consequently, to reproduce the CDM mock data with σ = 3 h⁻¹ Mpc under increasingly severe WDM scenarios, σ must decrease in order to counterbalance the damping induced by RSDs and preserve the overall signal amplitude (right panel of Fig. 9).

Fig. 9. Likelihood ratio ℓ computed for the pessimistic case with α = −1.9 in a convenient setup (R = 100, 160 deg2, and PWN ≃ 2.4 × 109 h−3 Mpc3 Jy2 sr−2) to highlight its reversed shape in the parameter space. The left and right panels are related in the same way as in Fig. 4.

We further find that the improvement in the m_WDM constraints remains modest relative to the results obtained for the pessimistic LF with α = −1.1 in the reference setup, owing to the large uncertainties associated with the mock data. The advantage of a steeper faint-end slope becomes more pronounced when we consider more favourable survey setups. For example, in a plausible near-future configuration with R = 100, 160 deg², and P_WN ≃ 2.4 × 10¹⁰ h⁻³ Mpc³ Jy² sr⁻², we obtain a 95% CL of 1.35 keV. This limit represents an improvement of approximately 75% relative to the pessimistic case with α = −1.1, exceeding even that derived under the optimistic LF assumption of ∼8% for the corresponding survey setup.

8. Conclusions

The nature of DM remains one of the central open questions in modern cosmology. While the CDM paradigm is the standard framework, alternative scenarios such as WDM composed of thermal relics provide viable solutions to some of the known challenges to ΛCDM. In this work we have built upon the model of M25, extending it to WDM cosmologies and developing a Bayesian framework to assess the constraining power of the [C II] PS on m_WDM. We presented forecasts for the FYST DSS at z ≃ 3.6 and explored the potential of future surveys with larger sky coverage, higher sensitivity, and/or increased spectral resolution. Our analysis led to the following main results:

• Under the assumption of a CDM background, the DSS in its reference setup is expected to set a 95% CL lower limit of 1.10 keV and 0.58 keV on the WDM particle mass in the optimistic and pessimistic (α = −1.1) scenarios, respectively. These bounds improve substantially with more ambitious setups, reaching 5.82 keV and 1.90 keV for half-sky coverage and an increase in sensitivity by a factor of $\sqrt{10}$ (see Fig. 6, with the lower panel showing the corresponding constraints expressed in terms of m_FDM).

• Enhancing the spectral resolution from R = 100 to R = 500 boosts sensitivity to the damping scale σ, thereby reducing its degeneracy with m_WDM. In favourable survey conditions and under the optimistic assumption for the LF, this translates into a tightening of constraints by a factor of up to ∼ 1.8. In the reference setup, however, the improvement is negligible, as the likelihood is anyway truncated by the prior imposed on σ (Fig. 4).

• By fitting mock data generated with m_WDM = 3 keV, we find that the true value can be recovered with both upper and lower bounds only in the optimistic LF scenario and under highly ideal survey configurations. For R = 100, this requires both half-sky coverage and enhanced sensitivity, while for R = 500, improved sensitivity combined with a 1600 deg² footprint is sufficient. In this latter case, for example, we obtain a 95% CL that lies in [2.62, 4.54] keV with uniform priors in m_WDM, and in [2.58, 4.12] keV with uniform priors in $m_{\mathrm{WDM}}^{-1}$, which demonstrates how the latter choice slightly favours the regime of the small WDM particle masses (see Fig. 7).

In deriving the results presented above, we explored several prior choices: (i) uniform in m_WDM, (ii) uniform in $\mathit{w}=m_{\mathrm{WDM}}^{-1}$, (iii) proportional to w^β, marginalising over the hyperparameter β ∈ [ − 0.9, 0], and (iv) a Jeffreys prior. The latter three yield consistent constraints, whereas the uniform prior in m_WDM leads to results that depend strongly on the choice of its upper bound (see Appendix B).

Further insights and additional considerations arising from the analysis are summarised below:

• The elongated likelihood contours obtained for α = −1.1 (Fig. 4) arise from the compensating effects of σ and m_WDM on the PS amplitude. This degeneracy stems from the form of the L(M) relation derived by matching an extrapolated LF – anchored to existing data but extended to faint, unobserved regimes – to increasingly extreme WDM HMFs, which predict progressively fewer low-mass haloes near and below the cutoff scale (see Sect. 4 and Appendix A).

• As discussed in Sect. 7.1, Fig. 8 shows that the LIM signal in terms of the [C II] PS is dominated by haloes with masses around 10¹¹–10¹² h⁻¹ M_⊙, while low-mass haloes, which mark the sharpest difference between CDM and WDM, contribute only marginally. This limits the sensitivity of current experiments and explains why very ambitious surveys are required to obtain competitive constraints. Nonetheless, adopting a steeper faint-end slope (α = −1.9, as discussed in Sect. 7.2) increases the relative contribution of faint sources, reversing the m_WDM–σ degeneracy (Fig. 9) and yielding tighter bounds. The effect is modest in the reference case but grows with survey improvements.

In summary, our forecasts underline both the potential and the current limitations of [C II] LIM as a probe of DM physics. While the DSS alone is unlikely to deliver stringent constraints, future surveys combining wider sky coverage, greater sensitivity, and higher resolution will significantly enhance its constraining power.

In this context, it is important to note that this work focuses on a single line and a narrow redshift range⁸ and therefore covers only a limited portion of the available observational domain. The resulting constraints should be regarded as somewhat conservative, as they do not yet exploit the full potential of forthcoming LIM surveys. Future LIM experiments will naturally provide access to multiple emission lines and extend coverage across a broad range of redshifts. Combining information from different epochs enhances the sensitivity to departures from CDM and simultaneously enables a more detailed characterisation of the evolution of the [C II] luminosity–mass relation. Moreover, by jointly analysing signals from multiple tracers, such as [C II], CO, Ly-α, or Hα, and combining measurements across several redshift intervals, it becomes possible to probe a wider range of halo masses and astrophysical environments, thereby obtaining tighter constraints on m_WDM. Although the exact quantitative gain depends on the degree of correlation between lines and the overlap in redshift coverage, multi-line and multi-epoch analyses are expected to substantially reduce uncertainties and yield more robust constraints on the nature of DM. These considerations highlight the strong potential of coordinated LIM surveys to overcome current limitations and maximise their cosmological return.

Acknowledgments

The authors gratefully acknowledge the Collaborative Research Center 1601 (SFB 1601 sub-project C6) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 500700252. They also acknowledge the International Max Planck Research School for Astronomy and Astrophysics (IMPRS A&A) at the Universities of Bonn and Cologne for supporting EM through a research contract. EM is a member of the IMPRS A&A, the Bonn Cologne Graduate School (BCGS), and guest researcher at the Max Planck Institute for Radio Astronomy (MPIfR) in Bonn. PK held the same affiliations (IMPRS A&A, BCGS, and MPIfR) at the time when most of this research was carried out. CP is grateful to SISSA, the University of Trieste, and IFPU, where part of this work was carried out, for hospitality and support. MV is partially supported by the Fondazione ICSC, Spoke 3 “Astrophysics and Cosmos Observations”, Piano Nazionale di Ripresa e Resilienza Project ID CN00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S (M4C2-19)” – Next Generation EU (NGEU). MV is also supported by the INAF Theory Grant “Cosmological Investigation of the Cosmic Web” and by the INFN INDARK grant. The research of AMD is supported by the Agence Nationale de la Recherche (ANR), grant ANR-23-CPJ1-0160-01.

References

Appendix A: Halo abundance and clustering in CDM and WDM models

We compare CDM and WDM models at z ≃ 3.6, highlighting their effects on the matter PS, HMF, and linear bias parameter (Fig. A.1).

Fig. A.1. Comparison between CDM and WDM models at z ≃ 3.6. Top panels: matter PS (left), HMF (centre), halo bias (right). The dotted black line corresponds to CDM, while the coloured lines represent WDM models with three representative particle masses. Bottom panels: Ratio of WDM to CDM, illustrating the deviation of each WDM model from the CDM baseline.

Appendix B: Uniform prior in m_WDM

For illustrative purposes, in this Appendix we combine the likelihood function obtained for the optimistic LF in our reference survey setup (see the left panel of Fig. 4) with a uniform prior in m_WDM defined over two different ranges. Because the likelihood becomes flat at large m_WDM, the inferred upper bound on the WDM particle mass depends strongly on the chosen maximum value $m_{\mathrm{WDM}}^{\max }$, thereby highlighting the prior dependence of this inference problem.

Fig. B.1. Marginalised posterior distributions (differential in mWDM) for the reference survey in the optimistic scenario. The posteriors are obtained assuming uniform priors in mWDM over two different mass ranges. Shaded regions indicate the 5% tails, defining the 95% CL for mWDM, corresponding to 1.57 keV and 6.16 keV for $m_{\mathrm{WDM}}^{\max }=10$ and 100 keV, respectively.

All Figures

Fig. 1. [C II] luminosity–halo mass relation inferred from abundance matching. In both panels, the solid lines correspond to the L(M) relations derived from our Schechter fits to the ALPINE data with fixed values for α. The dashed red line is based on the fit by Yan et al. (2020, Y20). The left and right panels assume a CDM and a 0.5  keV HMF, respectively (see Appendix A for a direct comparison between the two scenarios). In the left panel, for comparison, the dot-dashed black line indicates the L(M) relation obtained by Silva et al. (2015, S15), while the dotted gold and dark gold lines show results from the MARIGOLD simulations (Khatri et al. 2025, K25) at z = 5 and z = 4, respectively.

In the text

Fig. 2. Reconstructed [C II] LFs obtained by combining the abundance-matched L(M) relations with the HMFs of the CDM and WDM cosmologies. The figure uses the LF of our optimistic case and adopts a WDM particle mass of 0.5 keV. While the CDM model reproduces the input LF across all luminosities, the WDM case shows a sharp downturn at the faint end, reflecting the dearth of low-mass haloes and the corresponding minimum luminosity imposed by the abundance-matching procedure.

In the text

Fig. 3. Predicted [C II] PS. Left: Δ2(k, z ≃ 3.6) computed in the CDM scenario for our optimistic (solid red line) and pessimistic (solid blue line) cases for α = −1.1. The shaded areas represent the associated statistical uncertainty as for our survey reference setup. The dotted line is a graphical representation of the white-noise level, PWN, while the dashed and dot-dashed lines refer to the clustering and shot-noise components, respectively. Right: ratio between the WDM and CDM power spectra. The dashed and dotted black lines correspond to WDM models with mWDM = 10 keV and 0.5 keV, respectively. The shaded regions reflect the CDM uncertainties from the left panel. In both panels, the markers illustrate the adopted binning scheme.

In the text

	Fig. 8. Differential contribution of different halo masses to ${\overline{I}}_{\nu }$, b, and ${\overline{n}}_{\mathrm{eff}}^{-1}$ (from top to bottom) based on our halo model for CDM at z = 3.6.
In the text

Fig. 4. Likelihood ratio ℓ computed for the optimistic case. Left: ℓ in the reference setup. Centre: Same as the left panel, but with R = 500. Right: Graphical explanation of the banana-shaped ℓ shown in the left panel. Each pair of values (w, σ) in the legend indicates the inverse particle mass in keV−1 and the RSD displacement parameter in units of h−1 Mpc, respectively. The stars of corresponding colours in the left panel serve as identifiers for each combination in the parameter space, while the red circle denotes the true values used to generate the data.

In the text

Fig. 5. Marginalised posterior (over σ) computed for the whole set of sky areas in the optimistic scenario, with R = 100 and PWN ≃ 2.4 × 1010 h−3 Mpc3 Jy2 sr−2. Left: Posterior for a fixed β = 0 (uniform prior in w). Right: Result marginalised over β. The dashed lines of corresponding colours indicate the mWDM threshold up to which CDM and WDM can be distinguished (95% CL). The dotted lines show the assumed prior distribution for w.

In the text

Fig. 6. Lower bounds on mWDM (top) and mFDM (bottom) for R = 100. The pessimistic case is represented by the blue squares, while the optimistic case is denoted by the red diamonds. Filled markers correspond to the reference sensitivity; empty markers indicate a sensitivity increased by a factor of $\sqrt{10}$ (i.e. a white-noise PS decreased by a factor of ten). Shown for reference are the dashed green and dotted orange lines, marking the lower bounds inferred by Iršič et al. (2024) and Iršič et al. (2017a), respectively.

In the text

Fig. 7. Likelihood ratio ℓ (left) and marginalised posterior (right) obtained from 3 keV mock data in the optimistic scenario for a survey covering 1600 deg2 with R = 500 and $P_{\mathrm{WN}}\simeq 2.4\times {10}^9{h}^{-3}{{\mathrm{Mpc}}^3{\mathrm{Jy}}^2\mathrm{sr}}^{-2}$. In the right panel, the dark-green and orange vertical lines denote the 95% credible regions, computed assuming uniform priors in mWDM and $m_{\mathrm{WDM}}^{-1}$, respectively.

In the text

	Fig. 9. Likelihood ratio ℓ computed for the pessimistic case with α = −1.9 in a convenient setup (R = 100, 160 deg2, and PWN ≃ 2.4 × 109 h−3 Mpc3 Jy2 sr−2) to highlight its reversed shape in the parameter space. The left and right panels are related in the same way as in Fig. 4.
In the text

	Fig. A.1. Comparison between CDM and WDM models at z ≃ 3.6. Top panels: matter PS (left), HMF (centre), halo bias (right). The dotted black line corresponds to CDM, while the coloured lines represent WDM models with three representative particle masses. Bottom panels: Ratio of WDM to CDM, illustrating the deviation of each WDM model from the CDM baseline.
In the text

	Fig. B.1. Marginalised posterior distributions (differential in mWDM) for the reference survey in the optimistic scenario. The posteriors are obtained assuming uniform priors in mWDM over two different mass ranges. Shaded regions indicate the 5% tails, defining the 95% CL for mWDM, corresponding to 1.57 keV and 6.16 keV for $m_{\mathrm{WDM}}^{\max }=10$ and 100 keV, respectively.
In the text