Mixed Dark Matter and Galaxy Clustering: The Importance of Relative Perturbations

TL;DR

This paper develops a perturbative model for large-scale structure in cosmologies with mixed cold and warm dark matter, highlighting the importance of relative perturbations for accurate galaxy clustering analysis and constraints.

Contribution

It introduces linear solutions for relative perturbations in mixed dark matter models and demonstrates their impact on galaxy bias and power spectrum analysis.

Findings

Neglecting relative perturbations biases warm dark matter constraints.

Upcoming surveys can detect signatures of mixed dark matter scenarios.

Light WDM masses significantly affect galaxy clustering measurements.

Abstract

We develop a perturbative model to describe large-scale structure in cosmologies where dark matter consists of a mixture of cold (CDM) and warm (WDM) components. In such mixed dark matter (MDM) scenarios, even a subdominant warm component can introduce distinctive signatures via its free-streaming effects, altering the evolution of density and velocity perturbations. We present linear-order solutions for both total and relative perturbations in the two-fluid system, identifying novel contributions to galaxy bias caused by the relative density and velocity modes between the components. Incorporating these effects into the galaxy bias expansion, we compute the linear galaxy power spectrum in both real and redshift space. Using Fisher matrix forecasts, we assess the sensitivity of upcoming surveys such as DESI and PFS to MDM scenarios. Our results demonstrate that neglecting relative perturbations can lead to significant biases in inferred constraints on the warm dark matter fraction, particularly for lighter WDM masses ($\lesssim 150~\mathrm{eV}$ and $\lesssim 80~\mathrm{eV}$) for PFS and DESI, respectively. This framework provides a consistent and generalizable approach for incorporating multi-component dark matter dynamics into galaxy clustering analyses.

Equations

$$k_{\mathrm{fs}}^{-1}(t)=\frac{{\lambda }_{\mathrm{fs}}(t)}{2\pi }\sim \frac{⟨v^2(t){⟩}^{1/2}}{\mathcal{H}(t)},$$

$$f(p,z)=\frac{A}{e^{p/T_{\nu }(z)}+1},$$

$$f_w\equiv \frac{{\Omega }_w}{{\Omega }_c+{\Omega }_w}=\frac{{\Omega }_w}{{\Omega }_{\mathrm{DM}}},$$

$$k_{\mathrm{ph}}^{-1}\equiv {\int }_0^{t_{\mathrm{obs}}}\frac{⟨v^2(t){⟩}^{1/2}}{a(t)}dt,$$

$$\frac{\partial f}{\partial \eta }+\frac{\mathbf{p}}{am}\cdot {\nabla }_{\mathbf{x}}f-am{\nabla }_{\mathbf{x}}\Phi \cdot {\nabla }_{\mathbf{p}}f=0,$$

$$\rho =m\int d^{3}pf(\mathbf{x},\mathbf{p},\eta ),$$

$$\rho v_i=\int d^{3}p\frac{p_i}{am}f,$$

$${\sigma }_{ij}+\rho v_i{v}_j=\int d^{3}p\frac{p_i{p}_j}{a^2{m}^2}f,$$

$$\frac{\partial \delta }{\partial \eta }+\nabla \cdot [(1+\delta )\mathbf{u}]=0,$$

$$\frac{\partial \theta }{\partial \eta }+\mathcal{H}\theta +\nabla \cdot (\mathbf{u}\cdot \nabla \mathbf{u})=-{\nabla }^2\Phi -\frac{1}{\rho }{\nabla }^2(\rho {\sigma }_{ij}).$$

$${\sigma }_{ij}=\mathcal{P}{\delta }_{ij}+{\pi }_{ij},$$

$$\frac{\partial \theta }{\partial \eta }+\mathcal{H}\theta +\nabla \cdot (\mathbf{u}\cdot \nabla \mathbf{u})=-{\nabla }^2\Phi -\frac{1}{\rho }{\nabla }^2\mathcal{P},$$

$${\nabla }^2\Phi =4\pi G{a}^2\bar{\rho }{\delta }_m.$$

$$c_s^2\equiv \frac{\partial \mathcal{P}}{\partial \rho }.$$

$$\frac{\partial }{\partial \eta }{\delta }_c+{\theta }_c=0$$

$$\frac{\partial }{\partial \eta }{\theta }_c+\mathcal{H}{\theta }_c+\frac{3}{2}{\Omega }_m(a){\mathcal{H}}^2{\delta }_m=0$$

$$\frac{\partial }{\partial \eta }{\delta }_w+{\theta }_w=0$$

$$\frac{\partial }{\partial \eta }{\theta }_w+\mathcal{H}{\theta }_w+\frac{3}{2}{\Omega }_m(a){\mathcal{H}}^2{\delta }_m-c_s^2(a)k^2{\delta }_w=0,$$

$${\delta }_m^{\prime }+{\theta }_m=0$$

$${\theta }_m^{\prime }+\mathcal{H}{\theta }_m+\frac{3}{2}{\Omega }_m{\mathcal{H}}^2{\delta }_m-\frac{c_{s_i}^2{k}^2}{a^2}(f_w{\delta }_m-\cancel{f_w{f}_c{\delta }_r})=0$$

$${\delta }_r^{\prime }+{\theta }_r=0$$

$${\theta }_r^{\prime }+\mathcal{H}{\theta }_r+\frac{c_{s_i}^2{k}^2}{a^2}({\delta }_m-\cancel{f_c{\delta }_r})=0.$$

$${\delta }_m^{(1)}(\mathbf{k},a)=$$

$${\theta }_m^{(1)}(\mathbf{k},a)=$$

$$\text{mbox}with{R}_p(\mathbf{k})=$$

$${\delta }_r^{(1)}(\mathbf{k},a)$$

$${\theta }_r^{(1)}(\mathbf{k},a)$$

$${\delta }_g(\mathbf{k},\eta )=\sum _{\mathcal{O}}{b}_{\mathcal{O}}(\eta )\mathcal{O}(\mathbf{k},\eta ).$$

$${\delta }_g^{(1)}(\mathbf{x},\eta )$$

$$P_{gg}^{\mathrm{lin}}(k)=b_1^2{P}_{{\delta }_m{\delta }_m}^{\mathrm{lin}}(k)+2b_1{b}_{\delta }^r{P}_{R_0{\delta }_m}^{\mathrm{lin}}(k)+2b_1{b}_{\theta }^r{P}_{R_p{\delta }_m}^{\mathrm{lin}}(k)+P_{ϵ}.$$

$$b_{\delta }^r=\frac{1}{F_g[0]}{\frac{\partial F_g}{\partial {\delta }_r}|}_0,$$

$${\delta }_g^{(1)}(\mathbf{x},\eta )$$

$$b_{\delta }^r\in [(2f_w-1)b_1,f_w{b}_1].$$

$${\theta }_r\approx H_0^{-1}{\delta }_r,$$

Full text

1 Introduction

The $\Lambda$CDM model is the most widely accepted framework for describing the major components of the Universe. It posits the existence of a cosmological constant ($\Lambda$), responsible for the observed late-time acceleration of the Universe, and a dark matter component that is predominantly “cold”—meaning it has a negligible primordial velocity dispersion. Despite its simplicity, $\Lambda$CDM has been remarkably successful in explaining a wide range of cosmological observables, including the anisotropies of the Cosmic Microwave Background (CMB), which probe the early Universe, and the distribution of Large-Scale Structures (LSS), which trace the growth of structure at later times.

However, the increasing precision of cosmological observations has revealed several tensions that suggest $\Lambda$CDM may be an incomplete description of our Universe. For instance, the cold dark matter (CDM) assumption has been challenged by various observations probing different scales in the universe[1, 2]. On galactic scales, CDM simulations predict cuspy central density profiles in dark matter halos, whereas observations of dwarf and low-surface-brightness galaxies often favor cored profiles—a discrepancy known as the “core-cusp problem.” [3, 4] CDM also over-predicts the abundance of satellite galaxies around Milky Way–like halos, commonly referred to as the “missing satellites problem” [5, 6]. Furthermore, Lyman-$\alpha$ forest data and weak lensing measurements from surveys such as KiDS and DES suggest a possible suppression of small-scale power relative to $\Lambda$CDM predictions [7, 8].

There have been several theoretical and phenomenological studies on different alternatives of dark matter to overcome these possible issues with CDM. One class of alternatives are Axion like particles (ALP’s) that erase small-scale structures due to the wave-like behaviour that they have on small scales [9]. Another is warm dark matter (WDM), consisting of particles with non-negligible velocity dispersion that suppresses structure formation below a characteristic scale. The latter include various different theoretical models including thermal freeze-out [10], non-resonant (Dodelson-Widrow) [11] or resonant (Shi-Fuller) [12] sterile neutrino production, or decay of heavier particles [13] which also suppress small scale structure growth due to their non-negligible velocity dispersion.

Warm dark matter (WDM) models are tightly constrained by a range of cosmological and astrophysical observations. Some of the main methods include the Lyman-$\alpha$ forest power spectrum [14], which probes small-scale structure at high redshift, satellite galaxy counts [15], and strong lensing substructure analyses [16]. These probes are sensitive to the free-streaming effects of WDM, which suppress the formation of small-scale structure, thereby placing lower bounds on the WDM particle mass.

However, almost all constraints using the above methods assume that WDM makes up all of the dark matter. If WDM constitutes only a sub-dominant portion of the total dark matter density, many of these constraints become significantly relaxed. In mixed dark matter (MDM) models, a dominant cold component drives efficient structure formation, while a subdominant warm component subtly modifies growth, especially on intermediate and mildly nonlinear scales. Consequently, WDM particle masses that would be ruled out in pure-WDM scenarios can remain viable when a dominant cold component is present. Allowing a subdominant non-cold component can also introduce parameter degeneracies with other sectors; for example, [17] showed that MDM with a photon-interacting subcomponent can mimic the small-scale suppression from massive neutrinos in large-scale structure, complicating neutrino-mass inference in DESI-like analyses.

In this work, we investigate the phenomenology of MDM scenarios, where the dark matter sector contains both cold and warm components. Specifically, we aim to understand how a subdominant warm component, such as a sterile neutrino produced via the Dodelson-Widrow mechanism, affects large-scale structure observables. Our primary focus is on large-scale galaxy clustering, with the linear galaxy power spectrum in both real and redshift space serving as the main probe. This approach can be generalized to other cosmological observables such as Lyman-$\alpha$ forest power spectrum by using the selection dependent bias formalism developed in [18, 19].

MDM scenarios have been already explored as a way to reconcile small-scale tensions in $\Lambda$CDM in the literature. Several studies have investigated cosmological constraints on such models using different observables. Early work by [20] used Lyman-$\alpha$ forest data to place stringent limits on the warm fraction and mass, demonstrating that even a small warm component can significantly alter small-scale structure. The implementation of MDM scenarios in the Boltzmann solver CLASS [21, 22] enabled more recent analyses, such as [23, 24], who estimated non-linear structure formation in MDM models either with an analytic method or N-body simulations, and performed joint analyses with CMB, large-scale structure data and weak lensing measurements. Refs. [23, 24] provided updated bounds on the parameter space of thermal relics and sterile neutrino dark matter respectively. These studies highlight the importance of combining multiple observables to constrain MDM scenarios. However, they have so far typically relied on single-fluid treatments and do not account for the dynamical effects of relative perturbations between the species. Complementary non-linear treatments include [25], who built an emulator for MDM to predict the matter power spectrum and examined impacts on weak lensing, galaxy clustering, and the halo mass function, and [26], who calibrated a halo-model framework for mixed axion–CDM cosmologies against dedicated simulations.

A key feature of our approach is the inclusion of relative perturbations between the cold and warm components and their impact on galaxy bias. By relative perturbations, we refer to the minimal relative density and velocity perturbations arising from the initial difference in velocity dispersion between the two components, which naturally emerge due to the free-streaming of WDM. More generally, for example in the case of initial conditions generated by multi-field inflation, primordial iso-curvature perturbations between the components could also be present, potentially leaving additional imprints on large-scale structure; see [27] for a corresponding study of primordial baryon-CDM iso-curvature perturbations. Although we do not consider these in this work, they could be added straightforwardly.

We focus on the general bias expansion within the effective field theory of large-scale structure (EFTofLSS), which has been widely used for single-species dark matter. Within the EFT, the presence of multiple species with distinct perturbation histories introduces additional bias operators [28]. Throughout, we restrict our analysis to scales larger than the WDM free-streaming scale, $\lambda_{\rm fs}$, because extending the EFT to scales at or below $\lambda_{\rm fs}$ requires additional assumptions and/or an expansion in the WDM mass fraction [29, 30]. This choice confines our treatment to WDM particle masses $\gtrsim 10\,\mathrm{eV}$, without requiring the WDM mass fraction to be much less than unity. As such, our approach is complementary to Ref. [29].

The paper is structured as follows: In Section 2, we present the linear evolution of matter perturbations in the presence of both cold and warm dark matter, including a treatment of relative modes. In Section 3, we construct the linear galaxy power spectrum in both real and redshift space, accounting for the modified bias expansion. Section 4 presents our Fisher forecast analysis, where we assess the constraining power of DESI [31] and the Prime Focus Spectrograph (PFS) survey on MDM models with different warm component masses and fractions. We conclude in Section 5 with a summary and discussion of future directions.

2 Linear-Order Solutions for the Matter Field

2.1 WDM and Structure Formation

In this section, we review the description of WDM in the context of cosmological structure formation. First, it is important to emphasize that any constraint on WDM must be interpreted in the context of its production mechanism. Ultimately, observational constraints, particularly those from large-scale structure (LSS), do not directly probe the particle mass itself, but rather the suppression of power in the matter power spectrum caused by the free-streaming of the WDM component. This suppression depends not only on the mass, but also on the particle’s momentum distribution, which is determined by how it was produced. While we have referred to the scale $\lambda_{\mathrm{fs}}$, it is worth pausing to clarify its physical meaning. The free-streaming length at a given time is the typical distance that a warm particle travels within a Hubble time, and is defined as:

[TABLE]

where $\langle v^{2}(t)\rangle^{1/2}$ is the RMS WDM particle velocity dispersion, obtained by averaging over the phase-space distribution function, and $\mathcal{H}(t)=aH(t)$ is the comoving Hubble rate. This provides an estimate of the free-streaming scale, $\lambda_{\rm fs}$, below which no particle can gravitationally cluster—random motions erase perturbations on all smaller scales.

For example, consider two WDM candidates with the same mass of 100 eV but produced via different mechanisms. A thermal relic that decoupled from the plasma before Standard Model (SM) neutrinos would miss the subsequent entropy injections; see, e.g., Ref. [23], which presents a simple approximate relation between the relic temperature of a thermal WDM particle and its mass. As a result, its relic temperature today would be lower than that of the SM neutrinos, leading to a shorter free-streaming length and a milder suppression of small-scale structure. It also matters whether the relic is fermionic or bosonic, since Fermi–Dirac versus Bose–Einstein statistics change the relation between mass, temperature, and abundance, and thus the free-streaming scale. In contrast, a sterile neutrino produced via a non-thermal mechanism, such as the Dodelson-Widrow scenario, would share a relic temperature similar to that of SM neutrinos. Despite having the same mass (100 eV), its momentum distribution would be warmer, resulting in a significantly larger free-streaming scale and a stronger suppression of small-scale power.

This illustrates that mass bounds derived from LSS observables cannot be universally applied without specifying the underlying production model. The mapping between mass and free-streaming length is model-dependent, and failing to account for this can lead to misleading conclusions about the viability of a given WDM scenario.

This example highlights why non-resonantly produced sterile neutrinos are particularly convenient for our purposes. Unlike thermally produced warm particles—whose relic temperature depends on the specific decoupling history and must be redefined for each mass—sterile neutrinos generated through non-thermal mechanisms typically exhibit a relic temperature close to that of active neutrinos, regardless of their mass. Moreover, their initial momentum distribution is well approximated by a Fermi–Dirac function rescaled by a normalization factor [11],

[TABLE]

where $A$ is determined by the mean energy density of the warm component given its mass, and $T_{\nu}$ denotes the neutrino temperature today. It is convenient to express this energy density in terms of the fraction of the warm component,

[TABLE]

where $\Omega_{w}$ and $\Omega_{c}$ are the present-day density fractions of the warm and cold dark matter components, respectively. Notice that $\Omega_{w}=f_{w}\Omega_{\rm DM}$. In the sterile-neutrino case, specifying the particle mass $m$ and the fraction $f_{w}$ completely fixes the normalization $A$, making the model fully determined by these two parameters. The CLASS Boltzmann solver[22] includes this class of particles in its default configuration, allowing users to specify just the relic temperature, particle mass, and energy density (or $f_{w}$); the code then automatically computes the appropriate normalization for the distribution function.

In describing the impact of a warm component on the linear matter power spectrum $P_{\mathrm{lin}}$, it is useful to think in terms of two parameters—the mass $m_{w}$ and fraction $f_{w}$—and two characteristic scales—the free-streaming scale $k_{\rm fs}$ and the particle-horizon scale $k_{\rm ph}$. The particle-horizon scale quantifies the largest distance warm particles can travel over the age of the Universe and is given by

[TABLE]

where $a(t)$ is the scale factor and the integral runs from the early Universe to the observation time $t_{\mathrm{obs}}$. Figure 1 illustrates the relevant scales and their roles in shaping the linear matter power spectrum, indicating that the suppression starts to emerge at $k_{\mathrm{ph}}$, and reaches its maximum at $k_{\mathrm{fs}}$, where its amplitude is governed by $f_{w}$.

Even though the warm component initially has a significantly higher velocity dispersion compared to CDM, its temperature $T_{\nu}$ decays as $1/a$ as the universe expands. By the time it cools sufficiently (around $z\approx 30$), the evolution of perturbations for each component can be described using the fluid equations—the first two moments of the Boltzmann equation. The fluid approximation is valid on scales larger than the free-streaming scale, $k\ll k_{\rm fs}$, at least for thermal relics [32], which is the regime we will focus on in this paper; we will consequently also neglect the small-scale Poisson-like fluctuations that can be present for warm components [33]. Depending on the WDM mass, this puts an upper limit on the range of scales that we can describe, in addition to the nonlinear scale $k_{\rm nl}$ that limits the fluid approximation for the CDM component. Note that we do not assume that $f_{w}$ is much less than unity. We next derive the fluid equations describing the joint evolution of the CDM and WDM components.

2.2 Fluid Equations for Warm & Cold Dark Matter

The evolution of collisionless dark matter particles is governed by the Vlasov equation (also known as the collisionless Boltzmann equation), which follows from Liouville’s theorem applied to particles following Hamiltonian trajectories in phase space. This equation describes the time evolution of the phase space distribution function $f(\mathbf{x},\mathbf{p},\eta)$ for non-relativistic particles in comoving-Newtonian gauge:

[TABLE]

where $\Phi(\mathbf{x},\eta)$ is the Newtonian gravitational potential, $\mathbf{p}$ and $\mathbf{x}$ are the comoving momentum and comoving coordinates respectively, $\eta$ being the conformal time and $m$ is the mass of the dark matter particle. This form of the Vlasov equation applies to a single dark matter species, such as cold or warm dark matter considered individually. In scenarios with multiple components, each having distinct phase space distribution functions (e.g., different temperatures or momentum distributions), the equation must be applied separately to each component rather than to the total distribution. In this subsection, we consider equations for a single component (warm or cold), but refrain from adding a additional subscripts for clarity. The fluid approximation proceeds by taking moments of Eq. (2.5). We define the moments of the Vlasov equation over momentum space as follows:

• •

Density:

[TABLE]

• •

Mean (bulk) velocity:

[TABLE]

• •

Velocity dispersion tensor:

[TABLE]

where $\mathbf{v}=\mathbf{p}/(am)$ is the comoving peculiar velocity.

If one takes the zeroth moment of the Vlasov equation one obtains the continuity equation:

[TABLE]

where we have introduced the density contrast for the component considered as $\delta=\delta\rho/\bar{\rho}$. Likewise, taking the first moment of the Vlasov equation, and contracting with $\partial_{i}$, yields the divergence of Euler’s equation:

[TABLE]

where we have defined the velocity divergence as $\theta=\nabla\cdot\mathbf{v}$. We have neglected the vorticity component of the velocity field, as it decays with time at any order in perturbation theory, and only becomes relevant on fully nonlinear scales. The velocity dispersion tensor $\sigma_{ij}$ can be decomposed into an isotropic (diagonal) and anisotropic (traceless) part:

[TABLE]

where $\mathcal{P}$ is the effective pressure, and $\pi_{ij}$ is the anisotropic stress tensor, which accounts for the viscous behavior of the fluid. The latter is typically small and therefore often neglected; however, we retain the pressure term when considering the warm component. Substituting this into Eq. (2.10), we obtain the divergence of the Euler equation with pressure:

[TABLE]

where the gravitational potential is determined by the density contrast via the Poisson equation:

[TABLE]

Here, $\delta_{m}=f_{c}\delta_{c}+f_{w}\delta_{w}$ denotes the total matter density contrast, defined as the weighted sum of the density perturbations of two components where $f_{s}$ for $s\in\{c,w\}$ denotes the fractional contribution of each species to the total dark matter density. To close the system of equations, an equation of state for the pressure $\mathcal{P}$ is needed. We introduce an effective sound speed $c_{s}^{2}$ as

[TABLE]

This expression is analogous to the thermodynamic definition of the sound speed and characterizes the response of pressure to small density perturbations. For warm dark matter, the pressure does not arise from interparticle interactions but instead originates from the intrinsic velocity dispersion of the collisionless particles—determined by how they were generated in the early universe. At leading order, the effective sound speed $c_{s}$ decreases proportionally to $1/a$ as the Universe expands, which can be seen directly through Eq. (2.8). Consequently, the extent to which the pressure term suppresses structure formation depends on both the time of particle production and the initial velocity dispersion. These factors are closely tied to the underlying particle generation process discussed earlier. In the case of CDM particles we can neglect the initial velocity dispersion so that we do not need to include the pressure term; the effective pressure induced by nonlinear evolution of small-scale perturbations [34, 35, 36] is of the same order as nonlinear corrections, whose treatment we defer to future work.

2.3 Solutions of the Fluid Equations for a Two-Fluid System

Linearizing Eqs. (2.9)–(2.10) for the two components, CDM and WDM, respectively, we obtain the following set of equations to describe the linear evolution of the two components in our model:

[TABLE]

where we have moved to Fourier space for convenience and reformulated Eq. (2.13) with the help of the Friedmann equations [37]. Since gravity only cares about the total stress-energy, the gravity sourc term is associated with the total matter density fluctuations ($\delta_{m}$) for both fluids, which are hence coupled via gravity.

The evolution of two-fluid systems, such as the baryon–cold dark matter (CDM) mixture, has been studied extensively in the literature [38, 39, 40, 28]. In those cases, both components can be treated as pressureless perfect fluids (since the baryonic pressure is negligible on quasilinear scales), and no pressure term appears in the perturbation equations for the baryons. As a result, by switching from individual fluid perturbations to a basis of total matter and relative perturbations, the system of four coupled first-order differential equations can be reduced to two decoupled second-order equations.

In our case, however, this simplification cannot be applied directly due to the presence of a pressure term in the warm component. Despite this, it remains useful to adopt the same basis transformation. With appropriate approximations, the system can still be effectively described in terms of total and relative perturbations, even though the equations for relative and total perturbations remain partially coupled.

We now define our new basis as follows: the relative density and velocity divergence perturbation is given by $\delta_{r}\mathrel{\mathop{\ordinarycolon}}=\delta_{c}-\delta_{w}$, and $\theta_{r}\mathrel{\mathop{\ordinarycolon}}=\theta_{c}-\theta_{w}$. In the same way we defined the total matter density perturbations after Eq. (2.13) we define $\theta_{m}\mathrel{\mathop{\ordinarycolon}}=f_{c}\theta_{c}+f_{w}\theta_{w}$. After switching to the new basis Eq. (2.3) becomes:

[TABLE]

where $c_{s_{i}}$ is a dimensionless quantity defined as $c_{s}(a)=\frac{c_{s_{i}}}{a}$ which depends only on $k_{fs}$, and therefore implicitly on $m_{w}$.

In the full set of equations describing the evolution of total and relative perturbations in a two-fluid system with cold and warm components, the equations for $\delta_{m}$ and $\theta_{m}$ are coupled to those of $\delta_{r}$ and $\theta_{r}$ through the pressure gradient term. Specifically, the pressure depends on $\delta_{w}$, which can be expressed as a linear combination of $\delta_{m}$ and $\delta_{r}$. However, since the relative density perturbation $\delta_{r}$ is not sourced by gravity, it is not expected to grow significantly (we will confirm this below). In contrast, $\delta_{m}$ is directly sourced by the gravitational potential and dominates the clustering behavior. Therefore, as a physically motivated approximation, we neglect the subdominant contribution proportional to $\delta_{r}$ in the pressure term, as indicated above. This simplifies the system by decoupling the total and relative perturbations. Notice also that, following our assumption, $\delta_{r}$ is only sourced by free-streaming, and hence suppressed by $(k/k_{\rm fs})^{2}$ on large scales, so that the neglected contributions in Eq. (2.16) scale as $\sim k^{4}\delta_{m}$.

Adopting an Einstein-de Sitter (EdS) universe for the moment, the approximated form of the fluid equations immediately provide us with the following set of solutions:

[TABLE]

We fix the normalization of the growing mode by defining the amplitude of the adiabatic total-matter perturbation at an initial redshift $z_{\mathrm{ini}}=10$ (well into matter domination) as $\delta_{m}^{(0)}(\bm{k})\,\equiv\,\delta_{m,\mathrm{ad}}(\bm{k},z_{\mathrm{ini}})/D_{+}(z_{\mathrm{ini}})\,,$ where $D_{+}(z)$ is the linear growth factor normalized to $D_{+}(z=0)=1$, and $\delta_{m,\mathrm{ad}}$ denotes the adiabatic (growing) part of the total-matter perturbations. On the large scales used in our analysis ($k\ll k_{\mathrm{fs}}$), non-growing contributions are already subdominant at $z_{\mathrm{ini}}$, so $\delta_{m,\mathrm{ad}}(k,z_{\mathrm{ini}})$ can be taken from the CLASS total-matter transfer function. A simple consistency check is that $T_{\delta_{m}}(k,z)/D_{+}(z)$ is time-independent across two nearby redshifts around $z_{\mathrm{ini}}$, confirming isolation of the adiabatic growing mode.

Throughout, we work to leading order in the free-streaming (pressure) expansion, retaining only terms proportional to $(k/k_{\mathrm{fs}})^{2}$ and neglecting higher orders \mathcal{O}\!\big{[}(k/k_{\mathrm{fs}})^{4}\big{]}. We also neglect the adiabatic decaying mode in the total-matter perturbations, as it rapidly becomes irrelevant after matter domination and has no observable impact at the redshifts considered. While we display results using the EdS growth law, any small differences between $\Lambda$CDM and EdS growth on the large, perturbative scales and redshifts of interest can be absorbed into the galaxy-bias coefficients introduced in Section 3.

It is important to emphasize that a scale-independent growth factor $D_{+}(a)$, with $D_{+}(a)=a$ in EdS, is not valid across all scales, as the presence of the warm component leads to a suppression of growth on small scales. However, any scale-dependent behavior in the growth factor scales as $(k/k_{\mathrm{fs}})^{2}$, and is thus suppressed on scales larger than $k_{\mathrm{fs}}$. In Figure 2, we see that on perturbative scales the scale dependence of the growth rate amounts only to sub-percent level for MDM scenarios with a warm component of mass $m_{w}\gtrsim 5\ \mathrm{eV}$. The differences between redshifts $z=1$ and $z=2$ are minimal, indicating that the scale-dependent growth can be safely neglected at and above this mass at both redshifts.

After obtaining the solutions for the total matter perturbations in Eq. (2.3), we substitute them into the pressure term that enters the fluid equations for the relative perturbations. Although the full solution for $\delta_{m}$ contains multiple contributions—organized as a series in $(k^{2}/k_{\mathrm{fs}}^{2})$—it is sufficient to retain only the dominant growing mode, $\delta_{m}\approx a\,\delta_{m}^{(0)}$, when evaluating its effect on the relative modes. Including the leading pressure-sourced correction, $R_{p}$, would only generate contributions proportional to $k^{4}$, i.e., $\mathcal{O}\!\left[(k/k_{\mathrm{fs}})^{4}\right]$, which lie beyond our truncation and are therefore neglected. In this approximation, the pressure term acts as an external source in the equations for the relative perturbations. Solving these sourced equations then yields analytic expressions for the relative perturbations:

[TABLE]

We obtain a logarithmically growing relative-density mode which has not been considered in the literature before. One can verify the accuracy of our approximations with the help of numerical results obtained from CLASS. To do so, we extract transfer functions for $\theta_{r}$ at two different redshifts, denoted as $z_{1},z_{2}$. Since Eq. (2.3) includes two independent modes for $\theta_{r}$, namely $R_{p}(\bm{k})$ and $R_{-}(\bm{k})$, evolving respectively as $1/\sqrt{a}$ and $1/a$, this pair of snapshots allows us to independently determine the two modes. We repeat the same procedure for a second redshift pair, $z_{2},z_{3}$. If the time evolution derived in Eq. (2.3) is correct, the $k$-dependent mode amplitudes $R_{p}^{z_{1}z_{2}}$, $R_{p}^{z_{2}z_{3}}$ and $R_{-}^{z_{1}z_{2}}$, $R_{-}^{z_{2}z_{3}}$ extracted from each redshift pair should match each other precisely. Similarly, once $R_{p}$ and $R_{-}$ are known, one can infer $R_{0}$ from the numerical Boltzmann output for $\delta_{r}$, using Eq. (2.3), and verify consistency by comparing $R_{0}^{z_{1}z_{2}}$ and $R_{0}^{z_{2}z_{3}}$ obtained at different redshifts.

In Figure 3 we demonstrate this procedure, finding excellent agreement between the theoretical predictions and numerical results in an Einstein-de Sitter universe for sufficiently large scales. While percent-level deviations appear at small scales—more pronounced in the $R_{p}$ mode—we emphasize that our analysis adopts a conservative choice for $k_{\mathrm{max}}$, as discussed in Sec. 4. This ensures that we remain within the regime where the analytic mode decomposition is reliable and the underlying approximations are fully justified. Although probing mixtures with smaller masses may still be feasible, we restrict our analysis to $m_{w}\geq 10\ \mathrm{eV}$ as the lightest mass considered, which results in deviations of roughly $4\%$ for $R_{p}$ and $2\%$ for $R_{0}$ on perturbative scales.

Note that the growth functions in Eq. (2.3) will change when moving from Einstein-de Sitter to $\Lambda$CDM. However, this change is numerically small, especially as we will consider observations at $z\gtrsim 1$.

We now have the evolution of all the relevant modes necessary for an accurate description of a two-fluid system consisting of one cold and one warm component. However, before continuing it is worth discussing the hierarchy among these modes. As illustrated in Figure 4, the most dominant mode is, as expected, the growing total matter density mode, which serves as the reference here. This is followed by the constant and the logarithmically growing relative density modes, while the $R_{-}$ mode present in the relative perturbations (the green dashed curve) as well as the sub-leading total matter mode (the purple dashed curve) influenced by the pressure term, are much more suppressed. Consequently, we will omit the last two contributions in the following treatment. It is evident that, at any redshift, the sub-leading total matter mode will be supressed by a factor $\frac{1}{3}f_{w}$ compared to the $R_{p}$ mode.

3 Linear Galaxy Bias

Since we cannot directly observe matter density fluctuations, but instead observe tracers such as galaxies, we need a formalism that systematically connects the galaxy density field to the underlying matter density field. On large scales, this connection can be established using the perturbative bias expansion, as reviewed in [41]:

[TABLE]

Here, the $\mathcal{O}$ represent scalar contractions of the second derivative of the gravitational potential, $\partial_{i}\partial_{j}\Phi$, which is the leading local gravitational observable at each perturbative order. The coefficients $b_{\mathcal{O}}$ are free parameters associated with each operator in the expansion and must be determined from observational data.

This framework offers several important advantages. Since the bias expansion includes all local gravitational observables (up to a given order in perturbation theory), it does not require detailed assumptions about the astrophysical processes that govern galaxy formation and evolution. Instead, the complex small-scale physics is absorbed into a finite number of bias parameters, which can be systematically organized and constrained using data or simulations. Moreover, the expansion is consistent with cosmological perturbation theory and, on large scales where perturbation theory remains valid, allows for analytic control over predictions for observables such as the galaxy power spectrum and bispectrum.

Nonetheless, the approach has intrinsic limitations. The perturbative expansion breaks down on small, nonlinear scales where higher-order terms are no longer suppressed and non-perturbative effects dominate. Additionally, the bias coefficients are not directly observable and must be treated as nuisance parameters or fit empirically, which may introduce partial degeneracies and degrade cosmological parameter inference. Despite these caveats, the perturbative bias expansion remains a robust and the most widely used framework for interpreting galaxy clustering measurements in the linear and mildly nonlinear regime.

3.1 Linear Bias Relation

In the $\Lambda$CDM (or more generally quintessence-CDM) case, the matter field is treated as a single comoving fluid encompassing all relevant degrees of freedom. Consequently, all operators depend on a single mode: the adiabatic growing mode. The adiabatic decaying mode, which starts to decay shortly after matter domination, is generally neglected as we have done here [Eq. (2.3)].

In our case, however, we consider the evolution of two distinct matter components that are coupled through gravity, which is analogous to the baryon-CDM mixture considered in [38, 39, 40, 28, 42]. As discussed in Section 2, we neglect the decaying relative mode $R_{-}$ and the constant total matter mode as they are subdominant. This leaves us with two additional degrees of freedom in the bias expansion: a constant relative-density mode, $R_{0}$, and a relative-velocity and relative-density mode, $R_{p}$, which arises due to the effective pressure of the warm component.

The linear-order expression for the galaxy density field can then be written as:111Strictly, the relative modes $R_{0}$ and $R_{p}$ should be evaluated at the initial Lagrangian position $\bm{q}$ corresponding to a time at high redshift, consistent with the validity of our fluid description deep in the matter-dominated era. However, the distinction between Lagrangian and Eulerian coordinates is higher order in perturbation theory, so we can evaluate all fields at the Eulerian position $\bm{x}$ here. The distinction will become significant when considering nonlinear evolution.

[TABLE]

Here, $\epsilon(\bm{x},\eta)$ accounts for stochasticity, which we model as Gaussian white noise that is scale independent on large scales, as appropriate in the linear regime. See Ref. [29] for a bias expansion written in a different basis but equivalent at linear order, applied to a MDM scenario with standard dominant CDM and ALPs as the subdominant component.

Given the linear-order expression for the galaxy density field including the necessary relative modes, writing down the linear galaxy power spectrum is straightforward:

[TABLE]

where we have neglected the auto-correlation between the relative modes since their contribution is highly suppressed compared to $b_{1}^{2}P_{\delta_{m}\delta_{m}}^{\mathrm{lin}}$ (although it is easy to include them, and does not introduce any additional free parameters). Figure 5 illustrates the different contributions quantitatively. The power of $P_{\delta_{m}R_{0}}$ is around $15\%$ of $P_{\delta_{m}\delta_{m}}$ at $k\approx 0.3\,h\,{\rm Mpc}^{-1}$, which is non-negligible (this statement depends on the value of the bias coefficient $b_{\delta}^{r}$, which we discuss below). The stochastic field is not correlated with the other fields, therefore it only contributes via the stochastic power spectrum $P_{\epsilon}\sim n_{g}^{-1}$ which is roughly controlled by the number density of the galaxies in a survey.

In order to estimate the magnitude of the bias parameters, it is useful to recall their physical interpretation. Specifically, one can view the galaxy number density as a local functional of the CDM and WDM perturbations, such that $n_{g}(\bm{x},z)=F_{g}[\delta_{m},\delta_{r}]$. Each bias parameter can then be understood as the response of $n_{g}$ to a variation in the corresponding operator. The bias coefficient $b_{\delta}^{r}$ defined by

[TABLE]

quantifies the response of a galaxy population to a long-wavelength perturbation in the relative density field $\delta_{r}$ between two species, such as cold and warm dark matter. Physically, this measures how the abundance of galaxies changes when embedded in a large-scale environment where the relative density between the components is slightly enhanced or suppressed. The response is normalized by the mean abundance $F_{g}[0]$, yielding a dimensionless bias parameter that characterizes the sensitivity of galaxy formation to relative fluctuations.

In practice, this bias can be estimated using the separate universe approach. In this method, a small, constant $\delta_{r}$ perturbation is introduced into the background density in the initial conditions of a cosmological simulation, while keeping the total matter density fixed. Simulations with slightly positive and negative values of $\delta_{r}$ are run, and the resulting change in galaxy abundance is measured. The derivative in Eq. (3.4) is then extracted by comparing these outcomes, providing a direct numerical estimate of the bias of the simulated galaxy sample. This method was applied to measure the bias with respect to the baryon-CDM relative density perturbation in [27].

While it would be interesting to likewise estimate the bias coefficients $b_{\delta}^{r},b_{\theta}^{r}$ from full numerical simulations of MDM, here we restrict ourselves to providing analytical estimates. For this, it is more convenient to change our basis from $\delta_{m},\delta_{r}$ to $\delta_{c},\delta_{w}$. Eq. (3.2) can be equivalently written as:

[TABLE]

It is not trivial to estimate $b_{\delta}^{r}$ by Eq. (3.1) ; however, we can look at two limiting scenarios for the warm component to obtain a rough estimate of the allowed values of $b_{\delta}^{r}$. First, if $m_{w}\approx 0$, the warm component will free-stream on all scales, and will not participate in the formation of nonlinear structure. This implies $b_{w}\approx 0$ and $b_{c}\approx b_{1}$, leading to the approximation $b_{\delta}^{r}\approx f_{w}b_{1}$. On the other hand, if the free-streaming scale is significantly smaller than the size of halos which host the galaxies considered, so that the formation of these halos is not affected by the thermal velocities, the two species become indistinguishable, so we expect $b_{w}\approx b_{c}\approx b_{1}$. Under this condition, we obtain $b_{\delta}^{r}\approx(2f_{w}-1)b_{1}$. Therefore, $b_{\delta}^{r}$ should lie between these two limiting values:

[TABLE]

After setting the physical limits on $b_{\delta}^{r}$, we can relate it to $b_{\theta}^{r}$ using Eq. (2.3). A nonzero $\theta_{r}(\mathbf{x},z)$ is associated with a nonzero relative density perturbation $\delta_{r}$ derived above, which satisfies

[TABLE]

We can then estimate

[TABLE]

Notice that $b_{\theta}^{r}$ has dimensions of length, which is taken care of by the factor of $H_{0}^{-1}$. Assuming $b_{\delta}^{r}$ is taking values allowed in the Eq. (3.6), we therefore expect $b_{\theta}^{r}$ to be of order $H_{0}^{-1}$ as follows:

[TABLE]

3.2 Linear Velocity Bias

To accurately model the distribution of galaxies observed in redshift space, it is essential to account for galaxy velocities, as they affect the observed power spectrum through redshift-space distortions (RSD). In particular, when working within a two-fluid framework consisting of CDM and WDM, we must determine how perturbations in each fluid contribute to the total galaxy velocity.

In the single-fluid approximation, the equivalence principle ensures that the velocity bias between galaxies and matter vanishes at linear order. However, this no longer holds in the presence of two distinct dark matter components, where galaxies may not follow the total matter velocity exactly. This leads to a non-zero velocity bias at linear order due to the initial relative velocity between CDM and WDM. To incorporate this effect, we define the velocity bias $\bm{v}_{\rm bias}$ as the difference between the galaxy velocity $\bm{v}_{g}$ and the total matter velocity $\bm{v}_{m}$:

[TABLE]

Given the above discussion, we need to write the velocity bias at linear order as

[TABLE]

where the higher derivative terms correspond to the velocity bias in the usual single-fluid case, which are expected to be controlled by the physical length scale $R_{*}$ for galaxy formation [41]. One usually assumes that $R_{*}$ is of order the nonlinear scale or smaller, so that higher-derivative terms are compareable to nonlinear corrections, and we do not consider them further here. The coefficient of the first term, $\beta_{r}$, can be bracketed to take values in between $(-1,1)$, through the following argument. If galaxies purely follow the CDM component, then we have $\beta_{r}=f_{w}$, while in the opposite case, where galaxies only follow the warm component, we have $\beta_{r}=-f_{c}$. Given that the fractions can lie between 0 and 1, we obtain the above quoted overall range for $\beta_{r}$. This range could be narrowed more in a concrete analysis by factoring in the allowed region of $f_{w}$.

Now including the first term in Eq. (3.11), we obtain the linear galaxy power spectrum in redshift space as follows:

[TABLE]

where the first term reproduces the well known Kaiser formula [43], $f\equiv d\ln D_{+}/d\ln a$ being the linear growth rate (not to be confused with $f_{w}$) and $\mu$ the cosine of the angle between line of sight direction and wavevector $\bm{k}$. The last term is the contribution due to the linear velocity bias we have discussed above. Eq. (3.3) and Eq. (3.2) give our complete model describing the linear statistics of galaxy clustering both in real and redshift space.

4 Fisher Forecast

In this section we forecast the constraining power of ongoing galaxy redshift surveys on the MDM scenario considered in this work, parametrized by the mass $m_{w}$ and mass fraction $f_{w}$ of the WDM component. Fisher forecasting provides an efficient way to estimate expected parameter constraints for a given survey setup. A particular focus is on the effect of the relative density ($R_{0}$) and velocity ($R_{p}$) modes in the two-fluid cosmology on the constraints, obtained by marginalizing over the corresponding bias coefficients $b_{\delta}^{r},\ b_{\theta}^{r}$. We first consider the real-space galaxy power spectrum, before moving on to the redshift-space case.

We focus on the Subaru Prime Focus Spectrograph (PFS) [44] and DESI [31] surveys. PFS is smaller volume but probes higher redshifts than DESI. This comparison provides insights into how the impact of relative modes depends on survey redshift coverage and observational characteristics, and allows us to identify which survey is most sensitive to different mass ranges of warm particles.

4.1 Real Space Analysis

In a Fisher analysis, we estimate how well a set of model parameters can be constrained by examining the expected curvature of the likelihood function around its maximum. Under the assumption that the likelihood is Gaussian and that the data covariance matrix is diagonal and does not depend on the parameters being constrained, the Fisher information matrix is given by:

[TABLE]

where $\{\theta_{\alpha}\}$ is the set of model parameters, $P_{g}(k_{i})$ is the predicted galaxy power spectrum for a wavenumber bin $i$ centered at $k_{i}$, and $\mathrm{Var}[P_{g}(k_{i})]$ is the variance of the measured power spectrum in a given $k$-bin.

This expression relies on several key assumptions. First, we assume that the likelihood of the observables is Gaussian with diagonal covariance, which is a good approximation on linear scales where fluctuations are nearly Gaussian. To remain in this regime, on which our linear model for $P_{g}(k)$ also relies, we restrict our analysis to wavenumbers $k\lesssim 0.3\,k_{\mathrm{nl}}(z)$, where nonlinear corrections are still small and perturbation theory remains reliable. Here, $k_{\mathrm{nl}}(z)$ denotes the nonlinear scale at redshift $z$, defined implicitly by the condition

[TABLE]

where $P_{m}(k,z)$ is the linear matter power spectrum and $\Delta^{2}(k,z)$ is the variance of (linear) density perturbations. This criterion defines $k_{\mathrm{max}}$ in our analysis of the two surveys.

We focus on the redshift bin around $z=2$ for PFS and $z=1$ for DESI, correspondingly adopting $k_{\mathrm{max}}=0.30\,h\,\mathrm{Mpc}^{-1}$ for PFS and $k_{\mathrm{max}}=0.16\,h\,\mathrm{Mpc}^{-1}$ for DESI, given that $k_{\mathrm{nl}}(z=2)=1\,h\,\mathrm{Mpc}^{-1}$ and $k_{\mathrm{nl}}(z=1)=0.55\,h\,\mathrm{Mpc}^{-1}$. While these $k_{\mathrm{max}}$ values may appear somewhat optimistic, it is worth noting that contributions from the relative modes to the linear galaxy power spectrum scale as $(k/k_{\mathrm{fs}})^{2}$, which is already at the order of one-loop corrections. These scales can be safely modeled using (EFTofLSS) at one-loop order; however, going to one-loop order would introduce many additional free parameters in the galaxy power spectrum, which would ultimately degrade the constraining power of our analysis.

In addition to assuming a Gaussian likelihood for the observables, the Fisher matrix also approximates the posterior distribution of the parameters as Gaussian. This is a less accurate assumption, as posteriors can exhibit non-Gaussian features such as skewness or curved degeneracies, especially when parameters are not tightly constrained. As a result, this approximation can lead to underestimated uncertainties and should be interpreted with the appropriate caution.

Finally, the forecast depends on survey-specific inputs such as volume, redshift range, and galaxy number density. To see this dependence let us have a closer look at the variance [45]:

[TABLE]

where $N_{k_{i}}$ is the number of Fourier modes within each bin for which we evaluate the power spectrum in Eq.(3.3). Here we are not bothered with the details of the survey geometry and assume $N_{k_{i}}$ is the same as in case of a cubic survey with with a volume $V_{\mathrm{Survey}}$. In this simplified case one directly obtains

[TABLE]

where, $k_{F}=2\pi/V_{\mathrm{survey}}^{1/3}$ denotes the fundamental wavenumber of the survey volume. The factor of $1/2$ accounts for the reality of the density field, which implies that only half of the Fourier modes within a given shell are independent. For our analysis, we adopt $V_{\mathrm{survey}}=5.3\,(\mathrm{Gpc}/h)^{3}$ and a mean galaxy number density of $\bar{n}_{g}=3\times 10^{-4}\,(\mathrm{Mpc}/h)^{3}$ by combining the two largest redshift bins of the PFS survey [44]. Similarly, for DESI, we focus on four redshift bins centered around $z_{\mathrm{survey}}=1.1$. In this case, we take $V_{\mathrm{survey}}=40\,(\mathrm{Gpc}/h)^{3}$, $\bar{n}_{g}=3\times 10^{-4}\,(\mathrm{Mpc}/h)^{3}$.

In our analysis, we consider the following parameter vector:

[TABLE]

where the parameters—–$A_{s}$, $n_{s}$, $H_{0}$, $\Omega_{\mathrm{DM}}$, and $\Omega_{b}$—–are standard cosmological parameters describing the primordial power spectrum, background expansion and total dark matter/baryon energy densities respectively. $b_{1}$, $b_{\delta}^{r},b_{\theta}^{r}$ and $\beta_{r}$ are the bias coefficients which we treat as free parameters, and $f_{w}$ is the fraction of the warm component which we will constrain for individual WDM masses.

To improve parameter constraints and break degeneracies, we combine the Fisher matrix from the galaxy survey with a prior Fisher matrix from the Planck CMB experiment. The total Fisher matrix is

[TABLE]

where $F^{\mathrm{CMB}}_{\alpha\beta}$ is derived from the analysis of Planck data [46, 47]. This combination incorporates the additional constraining power of the early universe on the standard $\Lambda$CDM parameters.

In our analysis, the CMB information is used solely to constrain the $\Lambda$CDM sector, not the MDM parameters. This is justified because Planck probes scales much larger than $k_{\mathrm{fs}}$ at recombination, so the CMB effectively measures a single CDM component. While, in principle, CMB data could provide some sensitivity to $f_{w}$ through changes in $N_{\mathrm{eff}}$, this effect is negligible in our case: the sterile neutrinos are produced non-thermally and are already non-relativistic by recombination, making their contribution to $N_{\mathrm{eff}}$ minimal.

While the derivatives with respect to bias parameters in Eq. (4.1) can be performed analytically, in general these derivatives need to be evaluated numerically. In particular, the derivative with respect to $f_{w}$ must be computed numerically, since $f_{w}$ also modifies the transfer functions of the individual fluid components in a nontrivial way. This arises due to the nontrivial back-reaction between the CDM and the additional fluid component, which modifies the evolution of perturbations in both sectors.222 A recent study [48] introduced a modified Boltzmann solver called CLASSIER, an extension of the CLASS code we also use. CLASSIER is designed to handle non-cold relics with improved accuracy and may be a useful tool for cross-checking our results, especially in models with nontrivial multi-component dynamics.

In Fig. 6, we present the constraining power of the linear galaxy power spectrum for the PFS survey, comparing two models. The first (black contours) corresponds to a simplified scenario including only the total matter power spectrum contribution, i.e. it assumes that the bias coefficients for the relative modes vanish: $b_{\delta}^{r}=0=b_{\theta}^{r}$. The second (blue and red contours) incorporates the full model described in Eq. (3.3). For clarity, we omit the cosmological parameters listed in Eq. (4.5) from the plot; however, the full posterior distribution, including all parameter degeneracies, is provided in Appendix B. Fiducial cosmological parameters are taken from [46], with $b_{1}^{\mathrm{fid}}$ adopted from [44], and the warm dark matter fraction set to $f_{w}^{\mathrm{fid}}=0.05$ following [24]. We impose a Gaussian prior on $b_{\delta}^{r}$ with a standard deviation $\sigma_{b_{\delta}^{r}}=b_{1}^{\mathrm{fid}}$, motivated by the expectation that $b_{\delta}^{r}$ should not exceed the total bias by a large factor. Similarly, since $b_{\theta}^{r}$ is related to $b_{\delta}^{r}$ by a factor of $H_{0}^{-1}$, we apply a Gaussian prior with $\sigma_{b_{\theta}^{r}}=H_{0}^{-1}b_{1}^{\mathrm{fid}}$.

As seen in the figure, both models (no relative-mode bias vs. the full model) exhibit a positive correlation between $f_{w}$ and $b_{1}$; however, the degeneracy is more pronounced when relative-mode bias is allowed. This reflects a greater degree of compensation between the two parameters in the full model. The constraints on the warm component are generally strong for a 10 eV sterile neutrino in both cases, but introducing the additional nuisance parameters $b_{\delta}^{r}$ and $b_{\theta}^{r}$ leads to a significant weakening of the constraints on $f_{w}$—by roughly a factor of 2.5—compared to the simpler model with only the adiabatic growing mode. This degradation is primarily driven by a significant degeneracy between $f_{w}$ and $b_{\delta}^{r}$. The strong correlation between the parameters also suggests that neglecting the relative-mode bias in our model would induce a shift in the inferred value of $f_{w}$; a detailed calculation of this systematic shift will be presented explicitly in the next subsection 4.2. Given the amplitude hierarchy between $P_{\delta_{m}\delta_{m}}$ and $P_{\delta_{m}R_{0}}$, including additional observables may help break this degeneracy and better exploit the signal from relative perturbations—particularly for lighter secondary components.

We also observe that the constraint on $b_{\theta}^{r}$ is fully prior-dominated, while $b_{\delta}^{r}$ remains constrained by the data in the 10 eV case. Different fiducial choices for $b_{\delta}^{r}$ slightly affect the posterior for $f_{w}$, though this sensitivity is modest compared to the dominant effect of including the relative modes in the bias expansion. Lastly, we note that the constraints on $f_{w}$ are not sensitive to the choice of $f_{w}^{\mathrm{fid}}$, but they do depend on $b_{1}^{\mathrm{fid}}$, which significantly impacts the signal-to-noise ratio.

On the other hand, as the mass of the warm component is raised, the relative perturbations between the cold and warm particles become sufficiently small so that the data can no longer effectively constrain the nuisance parameter $b_{\delta}^{r}$ either. In Fig. 7, which shows the same as Fig. 6 but for $m_{w}=80\>\rm eV$, we observe that the uncertainty on $b_{\delta}^{r}$ is entirely dominated by the prior, similar to $b_{\theta}^{r}$. Moreover, for WDM with this mass, the inclusion of relative modes in the model becomes significantly less impactful compared to the case of lighter warm components. The change in $\sigma_{f_{w}}$ is at the level of $13\%$ for 80 eV particles, and this difference will diminish further for higher masses.

4.2 Redshift Space Analysis

We now extend our Fisher analysis to redshift space, where peculiar velocities of galaxies introduce anisotropies in the observed clustering pattern. This redshift space distortion (RSD) effect modifies the galaxy power spectrum by introducing a dependence not only on the wavenumber $k$, but also on the cosine of the angle between the wavevector and the line of sight, $\mu\equiv\cos\theta$. Consequently, the Fisher matrix formalism must be updated to reflect this anisotropic nature of the signal.

In redshift space, the Fisher matrix becomes:

[TABLE]

where $P_{g}(k_{i},\mu_{j})$ is the redshift-space galaxy power spectrum and $\mathrm{Var}[P_{g}(k_{i},\mu_{j})]$ is the variance in a given $(k,\mu)$-bin, which takes the form:

[TABLE]

$N_{k_{i}\mu_{j}}$ is the number of independent Fourier modes in the bin defined by $(k,\mu)$, given by:

[TABLE]

Here, we have considered joint binning in $k$ and $\mu$ (also known as “wedges”), while one could equivalently choose Legendre multipoles $P_{g,\ell}(k)$. The advantage of the chosen binning is that the covariance remains diagonal. We adopt 10 uniform bins in $\mu\in[0,1]$, which we find sufficient to capture the relevant angular structure; increasing the number of bins further yields negligible gain in constraining power.

In the redshift-space analysis, we introduce the velocity bias parameter $\beta_{r}$ as free parameter [Eq. (3.11)], which accounts for a possible linear velocity bias between galaxies and the underlying dark matter distribution. We choose a fiducial value for $\beta_{r}$ at the midpoint of the physically allowed range $[-f_{c},f_{w}]$, and impose a Gaussian prior such that the $1\sigma$ uncertainty spans this interval.

For the linear growth rate, we adopt the commonly used approximation $f\approx\Omega_{m}^{0.55}$, which is valid within a $\Lambda$CDM cosmology. Given that we do not expect scale-dependent growth in the linear perturbative regime of interest, this approximation remains valid for MDM scenarios considered here. While $f$ formally depends on the total matter density $\Omega_{m}=\Omega_{\rm{DM}}+\Omega_{b}$, we do not treat $\Omega_{m}$ as a free parameter in our analysis. Instead, we vary $\Omega_{\rm{DM}}$ and $\Omega_{b}$ separately. Since $\Omega_{b}$ is tightly constrained by CMB observations and contributes subdominantly to $\Omega_{m}$, we allow $f$ to vary only with $\Omega_{\rm DM}$ in our Fisher forecasts.

The main result of the analysis can be seen in Fig. 8. Focusing first on the constraints on $f_{w}$ from $P_{g}(k,\mu)$, we observe a slight improvement over the real-space case without RSD in both the full model—which includes relative modes and velocity bias—and the simplified model with only adiabatic modes. The trend of loosening the constraints on $f_{w}$ by including relative modes is also present in redshift space but the impact is not as strong as in the real space case.

The parameter $\beta_{r}$ remains poorly constrained due to its weak observational imprint. As shown in the figure, its posterior is largely dominated by the prior, indicating limited sensitivity in the data. Moreover, there are no significant degeneracies between $\beta_{r}$ and the other model parameters. As a result, treating $\beta_{r}$ as a free parameter does not degrade the constraints on the remaining parameters, despite its own large uncertainty. The information on $b_{\theta}^{r}$ is also prior dominated as in the real-space case. Another notable result is the improved constraint on $\Omega_{\mathrm{DM}}$, which arises thanks to the linear growth rate constraint from the large-scale anisotropic galaxy power spectrum.

Given the significant impact of the relative-mode bias on the forecasted error bar of $f_{w}$, it is worth asking what systematic bias in $f_{w}$ would result in a redshift-space power spectrum analysis when neglecting the relative-mode-bias contributions. To estimate this, we use the “Fisher bias” calculation, which yields a parameter shift of [37]

[TABLE]

where $\left(F^{-1}\right)_{f_{w}\beta}$ is the $(f_{w},\beta)$ element of the inverse Fisher matrix, and $B_{\beta}$ is the bias vector entry for parameter $\beta$:

[TABLE]

which is proportional to the difference between the full model and the no-relative-mode-bias case in each bin. Once $B_{\beta}$ is computed for all parameters in the no relative-bias model, the shift $\Delta f_{w}$ is obtained from Eq. (4.10) using the inverse Fisher matrix.

For the minimum-mass case, we find that the parameter shift in $f_{w}$ lies in the range $[0.2,\,3.7]\,\sigma_{f_{w}}$ for DESI, depending on the choice of fiducial values within the interval defined in Eq. (3.4). For PFS, the shift in $f_{w}$ ranges from $[0.3,\,5.56]\,\sigma_{f_{w}}$, again determined by the chosen fiducial value of $b_{\delta}^{r}$. This shows that neglecting the relative-mode bias contributions can lead to substantial mis-estimates of the WDM mass fraction $f_{w}$ for current surveys, so that these contributions are essential in real analyses. As expected, the significance of the parameter shift decreases for higher masses, as the relative modes become less significant.

So far, we have considered constraints at a fixed given WDM mass $m_{w}$. In an actual analysis, one would derive constraints on the mass fraction $f_{w}$ while scanning over $m_{w}$. Fig. 9 shows the expected constraints, $\sigma_{f_{w}}$, as a function of $m_{w}$ for the two different surveys, PFS and DESI. As expected, the constraints monotonically worsen as the mass is increased from the minimum considered here, $m_{w}=10\,\mathrm{eV}$. The mass fraction is no longer effectively constrained for $m_{w}\gtrsim 600\,\mathrm{eV}$ in the case of PFS, and for $m_{w}\gtrsim 400\,\mathrm{eV}$ for DESI. Although the comoving particle horizon scale $k_{\mathrm{ph}}$ does not change significantly between redshifts 2 and 1, the nonlinear scale $k_{\mathrm{nl}}$—and therefore the $k_{\mathrm{max}}$ that we use in our analysis—decreases at lower redshifts. This limits the ability of lower-redshift surveys to probe high-mass warm components, whose effects manifest on smaller scales.

In both surveys, we observe weakened constraints on $f_{w}$ for low-mass particles when bias with respect to relative modes is included (full vs open symbols). However, this trend persists up to $m_{w}\approx 150\,\mathrm{eV}$ for PFS, while it disappears at $m_{w}\approx 80\,\mathrm{eV}$ for DESI. This difference arises because relative modes become increasingly suppressed at late times, making them less relevant for low-redshift surveys like DESI. In light of our Fisher bias results, we expect systematic biases in $f_{w}$ to occur over the mass range in Fig. 9 where the relative-mode bias has an impact on constraints. Thus, in an actual analysis, these relative modes should always be included.

5 Discussion and Conclusions

In this work, we have derived predictions for large-scale galaxy clustering in MDM scenarios consisting of both cold and warm components. By incorporating relative perturbations between the cold and warm species—–specifically, the relative density and velocity modes—–we extended the standard galaxy bias expansion to account for two-fluid dynamics. Together with [29], this represents the first application of a full two-fluid perturbative bias formalism to MDM scenarios. Previous work on related MDM scenarios has neglected the contributions of relative perturbations to galaxy statistics [23, 24].

We derived analytical solutions for the total and relative perturbations using a fluid description, valid outside the free-streaming scale of WDM, $k\ll k_{\rm fs}$ (but with no restriction on the WDM mass fraction $f_{w}$), and showed that, in addition to a constant relative-density mode $R_{0}$, a logarithmically growing relative mode, $R_{p}$, emerges due to the effective pressure of warm dark matter. These modes imprint additional contributions to the galaxy power spectrum, characterized by new bias parameters $b_{\delta}^{r}$ and $b_{\theta}^{r}$. We incorporated these effects into both real- and redshift-space power spectra, enabling a consistent comparison with future observational data on large scales.

Using a Fisher matrix analysis, we forecast the constraining power of upcoming galaxy surveys such as PFS and DESI on the warm dark matter fraction $f_{w}$ for a range of particle masses. We found that including relative modes introduces parameter degeneracies—–particularly between $f_{w}$ and $b_{\delta}^{r}$—–which degrade constraints on $f_{w}$ by up to a factor of 2.5 for light warm components (e.g., $m_{w}=10\,\mathrm{eV}$). However, for heavier warm components ($m_{w}\gtrsim 80\,\mathrm{eV}$), the relative modes become suppressed, and their impact on parameter constraints diminishes significantly. We also examined the potential bias in $f_{w}$ that arises from incorrect assumptions about the relative bias $b_{\delta}^{r}$. In the minimum-mass case, these shifts can reach up to $3.7\sigma$ for DESI and $5.6\sigma$ for PFS, depending on the choice of fiducial values. These results highlight the sensitivity of $f_{w}$ to the inclusion of relative modes $R_{0},R_{p}$ particularly at low masses; as expected, the parameter bias becomes insignificant at high masses.

Our results show broad consistency with previous constraints in the literature, although our methodology differs in several aspects. In [20], constraints on MDM models were derived using Lyman-$\alpha$ forest data, which are particularly sensitive to small-scale power suppression induced by warm components. Their analysis demonstrated that even a modest warm fraction can lead to a significant suppression of structure on small scales. However, their treatment was based on a single-fluid approximation, which does not account for relative perturbations between the cold and warm components. It would thus be interesting to incoporate the relative perturbations in the model for Lyman-$\alpha$ forest data. Moreover, while Ref. [20] obtained tighter constraints for heavier warm particles, the constraints are weaker for lighter masses ($m_{w}\lesssim 100\,\mathrm{eV}$), primarily due to the limited sensitivity of Lyman-$\alpha$ forest observables to larger scales, where low-mass WDM leaves its main imprint. In contrast, our analysis based on large-scale galaxy clustering retains sensitivity to lighter warm species. This highlights the complementary nature of different large-scale structure probes. It is also worth highlighting that DESI and PFS are competitive in different WDM mass ranges (Fig. 9), providing a nice complementarity between the two surveys. This is due to their different redshift and number densities.

One of the key advantages of the general EFT bias framework adopted here is that it allows for a well-defined expansion at higher order in perturbation theory, and a joint prediction of other observables such as the galaxy bispectrum (three-point function) or cross-correlations with other tracers. We plan to pursue this in upcoming work. Our formalism can likewise be applied to other, non-WDM mixed scenarios (see [29] for the ALP case). Beyond extending to higher-order statistics, fast non-linear forward models for mixed dark matter developed in [25, 26] could be used to estimate the bias coefficients associated with relative perturbations by seperate universe approach. In practice, one can impose long-wavelength relative density/velocity backgrounds and use these emulation/halo-model tools to measure the response of halo abundance and clustering, thereby calibrating the relative-mode bias terms.

Acknowledgement

We thank Julien Lesgourgues and Julia Stadler for help with the CLASS numerics, and Emmanuele Castorina, Diego Redigolo, Ennio Salvioni, and Francesco Verdiani for discussions. ŞÇ thanks Beatriz Tucci Schiewaldt and Alejandro Pérez Fernández for their help with the Fisher analysis, and Eiichiro Komatsu, Toshiki Kurita, and Ivana Nikolac for insightful discussions. We thank the Galileo Galilei Institute for Theoretical Physics for the hospitality and the INFN for partial support during the completion of this work.

Appendix A Details of the CLASS transfer functions

In this section, we describe the key steps required to extract the relevant modes needed to build our model using CLASS. Once the code is successfully installed, working with mixtures composed of ordinary cold dark matter (CDM) and non-resonantly produced sterile neutrinos becomes straightforward. The first step is to define the cosmological parameters as follows:

⬇

cosmo_values = {

’omega_b’: fiducial_values[’omega_b’],

’N_ncdm’: 1,

’m_ncdm’: 10, # in units of eV

’omega_cdm’: fiducial_values[’omega_d’] * (1 - fiducial_values[’f_w’]),

’omega_ncdm’: fiducial_values[’omega_d’] * fiducial_values[’f_w’],

’A_s’: fiducial_values[’A_s’],

’n_s’: fiducial_values[’n_s’],

’H0’: fiducial_values[’H0’],

’output’: ’tCl,mTk,vTk,mPk’,

’P_k_max_1/Mpc’: 10.0,

’z_pk’: 1000.0,

’gauge’: ’sync’,

’ncdm_fluid_trigger_tau_over_tau_k’: 51,

’ncdm_fluid_approximation’: 3

}

It is also possible to include multiple non-cold relic species by simply increasing the value of N_ncdm. Detailed information on the implementation of non-cold relics can be found in [21].

After running the code with the desired cosmological parameters, transfer functions for individual components can be accessed using the get_transfer(z) function. The total matter density and velocity divergence transfer functions are directly provided, but to isolate the relative modes, one must combine the individual transfer functions $T_{\delta_{c}},T_{\delta_{w}},T_{\theta_{c}},T_{\theta_{w}}$, which correspond to the density and velocity divergence perturbations of the cold and warm components, respectively. For example, the relative velocity divergence transfer function is given by:

[TABLE]

Since these two modes scale differently with time, one can compute the transfer functions at two different redshifts and solve algebraically to isolate $R_{p}(\bm{k})$ and $R_{-}(\bm{k})$.

CLASS also outputs the sound speed $c_{s}^{2}(a)$ of the warm component, which is useful for estimating pressure effects. Because $R_{p}(\bm{k})$ follows the time evolution of the total matter growing mode, with an additional dependence on $k^{2}$ and $c_{s_{i}}^{2}$, one can cross-check the extracted values at different redshifts as a consistency check.

Similarly, the transfer function for the relative density perturbation can be written as:

[TABLE]

Since $R_{p}(\bm{k})$ and $R_{-}(\bm{k})$ have already been determined, one can now easily solve for $R_{0}(\bm{k})$.

Finally, it is worth mentioning that the total matter sub-leading mode appearing in Eq. (2.3) sourced by pressure cannot be directly extracted using transfer functions. It is overshadowed by residual radiation contributions, making it practically impossible to isolate from the individual component transfer functions numerically.

Appendix B Supplementary material for detailed results

This section presents the Fisher parameter posteriors for all parameters listed in Eq. (4.5) for completeness. We only show the analysis of MDM with $m_{w}=10\>\rm eV$, in both real space (Fig. 10) and redshift space (Fig. 11), since this choice is the most informative in terms of parameter degeneracies. As can be seen from these figures, we do not find any significant change in the final constraints or in the structure of the parameter degeneracies for the standard cosmological parameters, whether the relative-mode bias is included or excluded in our model.