The Dark Matter equation of state through cosmic history
Michael Kopp Constantinos Skordis Daniel B. Thomas Stéphane Ilić





kopp@fzu.cz
CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2, Praha 8 Czech Republic
Department of Physics, University of Cyprus, 1, Panepistimiou Street, 2109, Aglantzia, Cyprus

skordis@fzu.cz
CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2, Praha 8 Czech Republic
Department of Physics, University of Cyprus, 1, Panepistimiou Street, 2109, Aglantzia, Cyprus

daniel.thomas-2@manchester.ac.uk
Jodrell Bank Centre for Astrophysics, School of Physics & Astronomy, The University of Manchester, Manchester M13 9PL, UK 
Department of Physics, University of Cyprus, 1, Panepistimiou Street, 2109, Aglantzia, Cyprus

ilic@fzu.cz
CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2, Praha 8
Czech Republic




Cold Dark Matter (CDM) is a crucial constituent of the current concordance cosmological model. 
Having a vanishing equation of state (EoS), its energy density scales with the inverse cosmic volume and  is thus uniquely described by a single number, its present abundance. 
We test the inverse cosmic volume law for Dark Matter (DM) by allowing its EoS to vary independently in eight redshift bins in the range z=10^5 and z=0.
We use the latest measurements of the Cosmic Microwave Background radiation from the Planck satellite and supplement them with
 Baryon Acoustic Oscillation (BAO) data from the 6dF and SDSS-III BOSS surveys, and with the Hubble Space Telescope (HST) key project data.
We find no evidence for nonzero EoS in any of the eight redshift bins. With Planck data alone, the DM abundance is most strongly constrained around matter-radiation 
equality ω^ eq_g =   0.1193^+0.0036_-0.0035 (95% c.l.), whereas its present day value is more weakly 
constrained = 0.16^+0.12_-0.10 (95% c.l.). Adding BAO or HST data does not significantly change the ω^ eq_g constraint, while
  tightens to 0.160^+0.069_-0.065 (95% c.l.) and 0.124^+0.081_-0.067 (95% c.l.)  respectively.
Our results constrain for the first time the level of “coldness” required of the DM across various cosmological epochs and show that the DM abundance is strictly positive at all times.


The Dark Matter equation of state through cosmic history
    Stéphane Ilić
    December 30, 2023
========================================================





  
Introduction
Cosmological observations indicate that there is insufficient baryonic matter in the Universe
for the correct description of physical processes, if gravitational laws are dictated by General Relativity. 
A natural explanation is that most of the matter fields interact negligibly with light,
and are thus called Dark Matter, but  can still be seen through their gravitational effect.



Dark Matter (DM) is generally thought to be a  stable particle (or particles) not part of the standard model, however
it has so far remained elusive <cit.>.
Cosmologically, it is usually modeled as Cold Dark Matter (CDM), which is part of the successful model 
that is consistent with observations of the Cosmic Microwave Background (CMB) (e.g. <cit.>), 
cosmic shear surveys (e.g. <cit.>), measurements of the background expansion such as BAO probes <cit.>,
  supernovae distance measurements <cit.> and the observed abundance of light elements <cit.>.

The CDM model is defined by a phase space distribution function satisfying the collisionless Boltzmann equation with 
initially vanishing velocity dispersion and curl. This leads to a background CDM density _c (a) ∝  a^-3 (a being the scale factor of the Universe) and 
equation of state (EoS) w=0 while the linearized density and velocity perturbations satisfy the continuity and pressureless Euler equations.[Note that we use w to indicate the EOS of DM and not the EOS of Dark Energy which we assume to be -1 as in ΛCDM.] 
The resulting model arises naturally in the Weakly Interacting Massive Particle (WIMP) paradigm: the candidate particles are effectively 
collisionless and typically have an EoS  w ∼ 10^-24 a^-2 <cit.>, thus 
well described by CDM. The QCD axion is another CDM candidate <cit.>.


Not all DM candidates fit into the CDM paradigm, for instance, warm DM <cit.>, 
ultra light axions <cit.>, collisionless massive neutrinos <cit.>, 
 self-interacting massive neutrinos <cit.>, Chaplygin gas <cit.> and
self-interacting DM <cit.>. In addition, DM may interact with other species such as
neutrinos <cit.>, photons <cit.>,
dark radiation <cit.> 
and Dark Energy <cit.>.

Rather than taking the CDM description for granted we consider it timely to examine whether 
the data itself supports any deviation from the CDM paradigm, and thus to further determine or constrain DM properties. 
For our purpose we use the Generalized Dark Matter (GDM) model, first proposed by W. Hu <cit.>. 
The phenomenology  of the GDM model has been thoroughly investigated in <cit.>, where a connection was found with more fundamental theories, including
those of a rich self-interacting dark sector.
In addition, the recent work on the Effective Field Theory of Large Scale Structure (EFTofLSS) <cit.> suggests that, 
even for an initially pressureless perfect fluid, the non-linearities that develop on small scales affect  
the cosmological background and large scale linear perturbations, creating an effective pressure and viscosity such as those found in GDM.

The GDM model has been used to constrain DM properties with either constant or 
specific time dependences of the parameters <cit.>. 
Here, we allow the DM EoS to vary more freely in time than  all previous studies.



  
The model
We consider a flat Universe with only scalar perturbations, see <cit.> for more details and notation.
The background density _g and pressure _g of the DM evolve according to the conservation law

    _g  = - 3 H (1+w) _g ,     _g=w _g ,

where H =ȧ/a is the Hubble parameter and the overdot denotes derivatives with respect to cosmic time t.
The parametric function w(t) is freely specifiable with the case w=0 corresponding to a CDM background (_g = _c). 
The GDM model has two further free functions, the speed of sound, c_s^2, and the (shear) viscosity, c_ vis^2. 
The EoS w is uncorrelated with the two perturbative parameters c_s^2 and c_ vis^2, as shown in <cit.>, 
thus in this work we set these to zero and denote this class of GDM models by wDM. Consequently, replacing  CDM  by wDM  in the ΛCDM model
 leads to ΛwDM.

With this choice, the perturbed wDM fluid equations for the density contrast δ_g and velocity perturbation θ_g are given by

    δ̇_g/1+w   = 3   H ( w δ_g/1+w  +3 a H  c_a^2  θ_g) -  (  1/2ḣ  -  1/a^2θ_g ) 
    θ̇_g     = -  H  θ_g  ,      c_a^2 = _g/_g = w - ẇ/3 H (1+w).

Here, c_a^2 is the adiabatic speed of sound and h is a metric perturbation in synchronous gauge <cit.>.
The Euler equation θ̇_g  = -  H  θ_g is identical to that of CDM, which implies the solution θ_g=0.
An example of wDM is the combination of CDM and Λ interpreted as a single fluid with w=-(1+_c/_Λ)^-1. 
A large degeneracy between Ω_Λ and w is thus expected at late times (see also <cit.>).







  
Methodology
The wDM fluid equations (<ref>), (<ref>) were implemented in the Boltzmann code CLASS <cit.> 
as in <cit.>.  A sufficiently general time-dependence of w was achieved by binning its evolution into N=8 scale factor bins, 
whose edges are ã_i = 10^{0,-1,-1.5,-2,-2.5,-3,-3.5,-4}. 
The bins were smoothly connected using w(a) =  w_i-w_i+1/2( ln (a/ã_i+1)/σ_a) +  w_i+w_i+1/2 
for a_i+1<a<a_i, with bin centers a_i=√(ã_iã_i+1) for 1≤ i≤ N-2 while a_0=1 and a_N-1 = 0.
Because of the aforementioned degeneracy of wDM with CDM and Λ, we chose a wider bin in the late Universe. 
 
The σ_a parameter controls the transition width between bins; it was set to 1/20 so that the transition is small compared to the bin width.
We tested that this choice does not affect our conclusions. 

We define a dimensionless scaled wDM density 

    ω_g ≡ a^3 _g 8 π G/3× (100 km/s/Mpc)^2.

When w=0 through cosmic history, ω_g is a constant equal to the conventional dimensionless CDM density ω_c.
In general however, ω_g varies over time and is fully determined by the N+1 parameters , w_i. 
We use the notation  ω_g^(i)=ω_g(a_i)  and similarly for other functions with subscripts, so that the present day DM abundance 
is = ω_g(a_0).  For functions without a subscript we instead write H_i = H(a_i) and w_i=w(a_i).

Our parameter constraints were obtained as in <cit.> and we present only brief details here. We used the Markov chain Monte Carlo 
code MontePython <cit.> and established convergence of the chains using the Gelman-Rubin criterion <cit.>.
Our total parameter set 
 
    (ω_b, , H_0,n_s, τ, ln 10^10 A_s, w_i)

consists of 6 ΛCDM parameters and the 8 values w_i. We denote the Λ wDM model with 8 bins as “” and 
 the previously studied model <cit.> with w=const as “”.
We assumed adiabatic initial conditions.



We used the Planck 2015 data release <cit.> of the CMB anisotropies power spectra, composed of the low-l T/E/B likelihood and the 
full TT/TE/EE high-l likelihood with the complete “not-lite” set of nuisance parameters.
 [For full details, see the Planck papers and wiki http://wiki.cosmos.esa.int/planckpla2015/index.php/.] 
These likelihoods combined are referred to as Planck Power Spectra (PPS). 
We also added selectively the HST key project prior on H_0 <cit.>,
BAO from the 6dF Galaxy Survey <cit.> and the Baryon Oscillation Spectroscopic Survey Sloan Digital Sky Survey <cit.>, and the 
Planck CMB lensing likelihood (respectively referred to as HST, BAO and Lens thereafter).

We set uniform priors on τ and H_0 such that 0.01<τ and 45 ≤ H_0 ≤ 90 respectively. We used the same priors on Planck nuisance parameters
  and the same neutrino treatment as in <cit.>. The helium fraction was set to Y_ He=0.24667 <cit.>. 




  
Results
Our main results are constraints on the time dependence of DM EoS w(a) and  abundance ω_g(a) shown in Figs. <ref> and  <ref>.
For comparison, we also show the constraints on the  model already discussed in <cit.>. 
We list the 95% confidence regions of all parameters in Table <ref>.

In Fig. <ref> we observe that lies in the 99% confidence region of the  model, which in turn lies in the 99% confidence region of the  model, 
such that the constraints are nested like the models themselves.
There is no evidence for significant deviations of the DM EoS from 0 at any time. Consequently, any model selection criteria will favor .




 The constraints on w are the strongest between a_6 and a_5 enclosing the matter-radiation equality a_ eq≃ 3×10^-4,
and are about a factor 2 weaker compared to the  model.
In other bins the constraints on w weaken significantly.   Adding 
the BAO or HST dataset has only a minor effect on  constraints and only 
tightens limits in the rightmost bin. 
 As was the case for the  model, <cit.>, 
adding CMB lensing does not significantly  improve the constraints.

Let us now compare in more detail the DM abundance ω_g(a) of the  and  models focussing only on the two dataset combinations PPS and PPS+BAO+Lens.
In Fig.<ref> we see that, like w(a), ω_g is most tightly constrained between a_6=10^-3.75 and a_5=10^-3.25, in fact almost as tightly as for the  model (see inset in
 Fig.<ref>).  Around a=0.4 there is a squeeze in the constraints of ω_g from PPS, which extends to a∼ (0.08,0.4) when BAO  or HST are included. 
 At all times a vanishing DM abundance (ω_g=0) is inconsistent with the data.
More quantitatively, we find for the  model at 95% c.l. ω^ eq_g =   0.1193^+0.0036_-0.0035 and =0.16^+0.12_-0.10 
with PPS only, whereas for PPS+BAO+Lens we get ω^ eq_g =   0.1189^+0.0032_-0.0033  and = 0.169^+0.067_-0.065. 
For -PPS+BAO+Lens we find ω^ eq_g=0.1193^+0.0026_-0.0026, whereas the result is ω_c=ω^ eq_g=0.1184^+0.0022_-0.0022, 
see <cit.>.



Consider the tightly constrained region around a_ eq, a_6<a<a_5, as shown in Fig.<ref>  (see also the inset in the same figure). 
As discussed in <cit.> the GDM abundance ω_g(a) and expansion rate H(a) in the early Universe determine
the time of matter radiation equality and thereby the amount of potential decay until recombination. This in turn sets the relative 
heights of the first few peaks of the CMB temperature angular power spectrum. Both the  and  models constrain ω_g 
around  a_ eq at a similar level (see above). The degeneracy between H_0 and ω_g in the  model <cit.> 
translates into a degeneracy between H_6 and  in the  model as seen in the left panel of Fig. <ref>. 
Indeed, the H_6- contours reveal how well the CMB constrains a combination of the expansion rate and  the abundance of DM  around a_ eq.
The degeneracy between w-ω_g <cit.> in the  model due to the same effect is also seen as
a degeneracy between  w_6 and  (right panel of Fig. <ref>), however, in the  model it is weakened as w_6 has only an indirect effect on
a_ eq, contrary to .  Similar correlations exist for H_5 - and - w_5 but in the opposite direction.







The origin of the squeeze around a∼ (0.08-0.4) is of an entirely different nature.
 The angular diameter distance d^_A to the last scattering surface is given by 
 
    d^_A  = a_∫_a_^1 d ln a  (aH)^-1 = a_ (η_0 - η_) 
     .

Here, the second equality has been written in terms of the conformal time (η = ∫ dt/a) today, η_0, and at the last scattering surface, η_.
The largest contribution to  d_A^ comes from the first ln a bin where η grows from ∼ 4000Mpc to ∼ 14000Mpc
and constitutes ∼  70% of the total.
 This contribution can be strongly constrained by geometric probes. 
Within the first scale factor bin, we have  H =  √( (a^-3 (1+w_0)-1)+H_0^2) (with H and H_0 in units of 100 km/s/Mpc).
Hence for a one-parameter family of  and w_0 the combination  (aH)^-1 is
 approximately constant. As this is the largest contributor to d^_A we expect w_0 and ω^(0)_g to be anticorrelated, as is indeed
observed in Fig.<ref> (lower left panel).
The inclusion of BAO or HST data significantly improves and shifts the constraints on η(a) and is in turn reflected in the ω_g and w constraints.




For the  model, PPS alone allows for very low H_0, as low as 45km/s/Mpc, corresponding to our hard prior on H_0, see the blue contours in the top left and bottom right panels of  
Fig.<ref>.  Adding BAO (red dotted lines) or HST (black dashed lines) shrinks the posteriors and also moves the mean of H_0 back towards the range consistent with the
  model <cit.>. BAO (and also HST) data leads to a degeneracy between H_0 and w_0 and between H_0 and , as in the  model. However, 
as the present day values of H and ω_g are no longer anchored to their early Universe values,
the degeneracy axis is rotated and the contours are not as flattened. 

In the middle panels of Fig.<ref> we display the 2D marginalized posteriors of the H_0-Ω_Λ and -Ω_Λ planes.
As ω_b  is well constrained, wDM and Λ are the only relevant species in the late (flat) Universe and
are expected to have their abundances anticorrelated.  Indeed, the parameter  is anticorrelated
with  Ω^(0)_Λ for all data sets. 
The combination of CDM and Λ may be modeled by wDM; in that model, however,  w changes steeply only within the w_0 bin so that this behavior is unaffected.
When BAO or HST data are included the slope and size of the contours change strongly as the late Universe behavior 
dissociates from the early Universe in the  model. The negative values of Ω^(0)_Λ are correlated with the low values of H_0,
and, while allowed by PPS, they disappear when H_0 is better constrained after including BAO or HST data.








  
Implications
In the wDM model the DM abundance ω_g may deviate from its expected (constant) CDM value throughout cosmic history,
causing only minimal changes to the clustering properties of DM. 
Hence, the constraints on w and ω_g are conservative. 
One could also conservatively allow for general c_s^2(a,k) and c_ vis^2(a,k) and marginalize over them. 
However, as w is almost uncorrelated 
 with c_s^2 and c_ vis^2, we expect such procedure to give constraints similar to those here. 
In the cases of warm DM and EFTofLSS the parameters w, c_s^2 and c_ vis^2 are interrelated so that the w constraints will be driven
 by c_s^2 and c_ vis^2, and hence, tightened further <cit.>.
Adding spatial curvature and/or neutrino mass would likely widen the ω_g constraints on the squeeze at a∼{0.08-0.4} <cit.> 
and in the latter case on the tightly constrained region around a_ eq as well.


When applying our constraints to generic theories of Dark Matter, including those coming from modifications of gravity, one must keep in mind our underlying assumption of adiabaticity.
As models of modified gravity will typically have additional fields leading to more types of isocurvature modes, we expect our constraints to be less applicable in those cases.
However, within our adiabatic assumption we expect our constraints to be valid for any theory of Dark Matter or modified gravity. The cosmological background in any such theory 
will have to evolve as in ΛCDM (see for example <cit.>), around matter radiation equality and before decoupling. 
Typical examples include DM-DE coupled models <cit.>. Explicit realizations where a CDM-like
 background decays into DE are given by the quasidilaton models of massive gravity <cit.> and by 
axion models <cit.>.










  
Conclusion
We have constrained the EoS w and abundance ω_g of Dark Matter, in 8 temporal bins covering 5 decades in cosmic scale factor, 
using the CMB data from the Planck satellite, and separately including BAO and HST data. 
We found that w is consistent with zero and the DM abundance is strictly positive at all cosmological epochs considered here, 
see Fig. <ref> and Fig. <ref>, and thus the  concordance model remains unchallenged.
This is the first time that the level of DM “coldness” across cosmic time has been explicitly constrained. 



The research leading to these results
has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC Grant Agreement n. 617656 “Theories
 and Models of the Dark Sector: Dark Matter, Dark Energy and Gravity”. The Primary Investigator is C. Skordis.



AJ
A&A
ApJ
A&A Rev.
ApJ
MNRAS
ARA&A
AJ
QJRAS
Phys. Rep.
Nature
A&A Supp.
Ap&SS      
ApJS
Phys. Rev. D
JCAP
New Astron. Rev




apsrev





