First Cosmological Constraint on the Effective Theory of Dark Matter-Proton Interactions
Kimberly K. Boddy Vera Gluscevic
astro-ph.CO hep-ph





Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

Department of Physics & Astronomy, University of Hawaii, Honolulu, HI 96822, USA


School of Natural Sciences, Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA

Department of Physics, University of Florida, Gainesville, Florida 32611, USA

Department of Physics, Princeton University, Princeton, NJ 08544



  We obtain the first cosmological constraints on interactions between dark matter and protons within the formalism of nonrelativistic effective field theory developed for direct detection.
For each interaction operator in the effective theory, parameterized by different powers of the relative velocity of the incoming particles, we use the Planck 2015 cosmic microwave background (CMB) temperature, polarization, and lensing anisotropy to set upper limits on the scattering cross section for all dark matter masses above 15 keV.
We find that for interactions associated with a stronger dependence on velocity, dark matter and baryons stay thermally coupled for longer, but the interaction strengths are suppressed at the low temperatures relevant for Planck observations and are thus less constrained.
At the same time, cross sections with stronger velocity dependencies are more constrained in the limit of small dark matter mass.
In all cases, the effect of dark matter–proton scattering is most prominent on small scales in the CMB power spectra and in the matter power spectrum, and we thus expect substantial improvement over the current limits with data from ground-based CMB experiments and galaxy surveys.

First Cosmological Constraint on the
 Effective Theory of Dark Matter–Proton Interactions
    Vera Gluscevic
    
=========================================================================================






§ INTRODUCTION


There are various astrophysical, cosmological, and collider-based probes engaged in an effort to detect interactions between dark matter (DM) and Standard Model particles.
The most sensitive low-energy constraints on DM interactions typically come from direct searches for rare collisions between DM particles from the local Galactic halo and atomic nuclei in low-background underground detectors <cit.>.
Nuclear-recoil measurements from these experiments exclude large portions of the relevant parameter space, especially for DM particle masses above a few  <cit.>.
Data from a new generation of direct-detection experiments are forthcoming from a number of projects worldwide and are expected to deliver another order of magnitude in sensitivity beyond the current detection limits <cit.>.

However, information about DM gained from direct-detection measurements is subject to several caveats.
First, the conversion of the observed nuclear-recoil rate into a limit on the DM–baryon interaction cross section relies on detailed knowledge of DM astrophysical parameters, such as the local energy density and velocity distribution.
N-body simulations show departures from the usual assumptions of a Maxwell-Boltzmann phase space distribution on small scales <cit.>, but the only way to directly measure these parameters would be through direct-detection experiments.
Second, nuclear recoil-based measurements are not sensitive to DM masses below a few GeV; sub-GeV DM would need to travel at speeds larger than the escape velocity of the Galactic halo in order to produce a detectable nuclear recoil <cit.>.
While electronic recoil measurements can probe the parameter space of sub-GeV DM <cit.>, new technologies are necessary to address scattering with nucleons in this mass regime <cit.>.
Finally, there is a “ceiling” on the cross sections that underground experiments are able to probe, due to their extensive shielding: DM particles that interact too strongly would lose momentum before reaching the detector <cit.>.

Alternative probes of low-energy DM physics can offer complementary information and may sidestep some of the aforementioned caveats of direct detection.
For example, direct-detection surface runs <cit.>, molecular spectroscopy <cit.>, and balloon-borne experiments <cit.> have a higher a “ceiling” and place bounds on large interactions between DM particles and baryons.
There are also constraints that arise from the Galactic structure <cit.>, observations of galaxy clusters <cit.>, cosmic rays <cit.>, and other astrophysical observations <cit.>.
Additionally, a recent proliferation of precision cosmological measurements provides a testing ground for the very same interaction physics, but in the early Universe.
Previous studies used the cosmic microwave background (CMB) spectral distortion limit from COBE FIRAS to constrain the DM–baryon scattering for DM masses less than 100  <cit.>.
References <cit.> used CMB temperature anisotropy measurements to constraint DM–proton interactions; Ref. <cit.> considered the limit in which DM is much heavier than the proton, and our previous work <cit.> covered masses between a keV and a TeV, but only for velocity-independent contact interactions.

In this study, we expand on previous work to search for cosmological evidence of any nonrelativistic DM–proton effective interaction (including velocity-dependent interactions) for particle masses above 15 .[Warm DM limits exclude masses below a few  <cit.>, and power spectra computations become progressively more difficult to perform at high accuracy for scattering with stronger velocity dependencies and for DM with lower masses. For these reasons, we focus on masses greater than 15 .]
For this purpose, we use the latest temperature, polarization, and lensing anisotropy measurements from the Planck 2015 data release <cit.>.
We adopt the effective theory formalism, originally developed for DM direct detection, which renders our results directly comparable to those from laboratory experiments and provides a framework to systematically investigate all possible low-energy DM interactions with protons.
We include DM scattering with helium nuclei, accounting for the nuclear responses triggered by different interaction operators.
We find no evidence for either velocity-independent or velocity-dependent DM scattering and thus present state-of-the-art cosmological constraints on DM–proton interactions; the key results are summarized in Figure <ref>.

In Sec. <ref>, we review the nonrelativistic effective theory of DM interactions with baryons.
In Sec. <ref>, we embed this formalism into the Boltzmann equations that describe the evolution of cosmological perturbations, allowing for the presence of DM–proton scattering in the early Universe.
In Sec. <ref>, we describe the data and present our results.
In Sec. <ref>, we discuss our results and future avenues of investigation.




§ DARK MATTER EFFECTIVE FIELD THEORY



The effective field theory (EFT) for DM interactions with nucleons enables a systematic description of processes relevant for probes of low-energy DM physics <cit.>.
In Sec. <ref>, we summarize the general EFT approach following Refs. <cit.> and express the scattering amplitude for each interaction operator in a form that is useful in a cosmological setting.
In Sec. <ref>, we derive the associated momentum-transfer cross sections, relevant for investigating the effect of DM scattering on cosmological observables.




 §.§ Formalism



We begin by considering the nonrelativistic elastic scattering between a DM particle and a nucleon.
The complete set of Hermitian observables that describes the scattering process is as follows <cit.>: iq⃗/m_N is the momentum transfer per nucleon mass; v⃗^⊥≡v⃗ + q⃗/(2μ_χ N) is the relative velocity in a direction perpendicular to the momentum exchange (and μ_χ N is the reduced mass of the DM–nucleon system); and S⃗_χ and S⃗_N are the DM spin and the nucleon spin, respectively.
The momentum transfer incorporates the angular information of the scattering process via |q⃗|^2=2μ_χ N^2 v^2 (1-cosθ), where v≡|v⃗| and θ is the scattering angle in the center-of-mass frame.
The maximum possible momentum transfer is |q⃗|_max=2μ_χ N v.
Working to second order in momenta and velocities, various combinations of these four quantities give rise to the following 14 operators,[Following Ref. <cit.>, we fix c_2^τ=0 in order to omit 𝒪_2=|v⃗^⊥|^2, as it does not arise at leading order from the nonrelativistic reduction of a relativistic operator; conversely, we keep 𝒪_15, which is third order in momenta and velocities, because it can arise at leading order.] derived in Ref. <cit.>:

    𝒪_1       = 1_χ 1_N    𝒪_9       = S⃗_χ·(S⃗_N ×iq⃗/m_N) 
    𝒪_3       = S⃗_N ·(iq⃗/m_N×v⃗^⊥)    𝒪_10   = S⃗_N ·iq⃗/m_N
    𝒪_4       = S⃗_χ·S⃗_N    𝒪_11   = S⃗_χ·iq⃗/m_N
    𝒪_5       = S⃗_χ·(iq⃗/m_N×v⃗^⊥)    𝒪_12   = S⃗_χ·(S⃗_N ×v⃗^⊥) 
    𝒪_6       = -(S⃗_χ·iq⃗/m_N) (S⃗_N ·iq⃗/m_N)    𝒪_13   = (S⃗_χ·v⃗^⊥) (S⃗_N ·iq⃗/m_N) 
    𝒪_7       = S⃗_N ·v⃗^⊥   𝒪_14   = (S⃗_χ·iq⃗/m_N) (S⃗_N ·v⃗^⊥) 
    𝒪_8       = S⃗_χ·v⃗^⊥   𝒪_15   = (S⃗_χ·iq⃗/m_N) [(S⃗_N ×v⃗^⊥) ·iq⃗/m_N]  .

Since DM may interact with both protons and neutrons, it is convenient to work in an isospin basis, in which the interaction Hamiltonian has the form

    ℋ = ∑_τ=0^1 ∑_i=1^15 c_i^τ𝒪_i t^τ ,

where i labels the interaction and τ labels isospin.
The isospin operators are t^0≡1 and t^1≡τ_3, and the isospin coupling coefficients are c_i^0 and c_i^1.
The coupling coefficients to protons c_i^(p)=(c_i^0 + c_i^1)/2 and neutrons c_i^(n)=(c_i^0 - c_i^1)/2 set the strength of the DM interactions with the corresponding nucleon.
For a given relativistic theory, reducing the DM–nucleon interaction into its nonrelativistic counterpart generically yields a linear combination of operators.
For example, 𝒪_1 through 𝒪_11 are associated with interactions that occur through the exchange of a heavy spin-0 or spin-1 mediator in a relativistic theory, and the coupling coefficients as written may depend on factors of |q⃗|^2.
In this study, we avoid choosing specific underlying theories to maintain generality and thus work directly with the operators in Eq. (<ref>), treating the coupling coefficients as constants.
As detailed in Sec. <ref>, by considering individual operators, we can study the cosmological effects of DM–nucleon interactions in a systematic manner and place conservative upper limits on each individual coupling coefficient (neglecting operator interference, which we discuss at the end of this section).
We previously constrained interactions via 𝒪_1 and 𝒪_4 (referred to as the standard “spin-independent” and “spin-dependent” interactions), corresponding to the simplest case of n=0 <cit.>.
We expand upon that work by investigating the remaining 12 operators.
Constraining full relativistic theories is left for future work.

Thus far, we have focused on DM scattering with individual nucleons.
For the purposes of computing signals in direct-detection experiments, these nucleons are embedded in atomic nuclei, and the full response of the nucleus must be treated appropriately.
In a cosmological setting, we are interested in hydrogen and helium nuclei, which dominate the energy density of baryons in the Universe.
Scattering with helium was either neglected or treated too simplistically in previous literature <cit.>; here, we incorporate scattering with helium, accounting for the velocity dependence that arises from its composite nature.
We use the label B∈{p,He} to denote the baryonic species—either the proton or the helium nucleus.
The Hermitian observable S⃗_N becomes S⃗_B and v⃗^⊥ becomes v⃗^⊥_B (which depends on the reduced mass of the DM–B system, μ_χ B).





In Refs. <cit.>, DM interactions with composite nuclei are approximated as the sum of interactions with the individual nucleons.
Additionally, the nuclear wave functions are taken to have the standard shell-model form, and the underlying single-particle basis is the harmonic oscillator with a parametric size

    a_B = √(41.467/(45 A^-1/3 - 25 A^-2/3))  ,

where A>1 is the atomic number.
The proton is pointlike and thus a_p=0.
The effect of compositeness when scattering with a nonpointlike nucleus is encoded in the nuclear response function W_B,k^ττ^'(y), where y≡(|q⃗|a_B/2)^2.
The index k labels the type of response (not to be confused with the wave number defined in Sec. <ref>), as in the standard treatment of semileptonic weak interactions <cit.>.
There is also a DM response function R_k^ττ^' (v⃗_B^⊥ 2,|q⃗|^2/m_N^2), which incorporates the coupling coefficients of the various operators.
The resulting spin-averaged amplitude squared for the scattering between DM and the baryon B is <cit.>

    _B = 1/m_v^44π/2S_B + 1∑_τ, τ^'∑_k
      (|q⃗|^2/m_N^2)^ξ_k
      R_k^ττ^'(v⃗_B^⊥ 2,|q⃗|^2/m_N^2)
      W_B,k^ττ^'(y)  ,

where ξ_k=0 for k∈{M, Σ^'', Σ^'} and ξ_k=1 for k∈{Φ^'',Φ^''M,Φ^',Δ,ΔΣ^'}.
We have inserted the mass of the weak scale m_v≡(√(2)G_F)^-1/2≈246 to render the coupling coefficients dimensionless.[This mass scale is an arbitrary normalization, and it does not impact the numerical value for the constraints on the cross sections reported in Sec. <ref>.]
We refer the reader to Ref. <cit.> for the full expression for R_k^ττ^' and to Refs. <cit.> for the form of W_B,k^ττ^' for various elements.

In the context of cosmology, we need only consider scattering on protons (with response types M, Σ^'', and Σ^') and helium nuclei (with response type M).
It will also be useful to isolate the velocity and angular dependencies in Eq. (<ref>).
For this purpose, we first define the dimensionless quantity x≡|q⃗|^2/|q⃗|^2_max and make the following substitutions:

    |q⃗|^2/m_N^2 = x v^2 (2μ_χ B/m_N)^2  ,     
      |v⃗_B^⊥|^2 = v^2(1-x)  ,     
      y = x v^2(μ_χ B a_B)^2  .

Next, we express the DM response function as a product of a numerical coefficient ℛ_k,ij, the coupling coefficients of 𝒪_i and 𝒪_j, and powers of |v⃗|_B^⊥ 2 and (|q⃗|/m_N)^2, such that

    R_k^ττ^'
      = ∑_i,jc_i^τ c_j^τ^'ℛ_k,ij(|v⃗|_B^⊥ 2)^α_ij(|q⃗|^2/m_N^2)^β_ij
      = ∑_i,j c_i^τ c_j^τ^'ℛ_k,ij(2μ_χ B/m_N)^2β_ij
      v^2(α_ij+β_ij) (1-x)^α_ij x^β_ij .

The values of α_ij, β_ij, and ℛ_k,ij are listed in Table <ref>.[For the most generic form of R_k^ττ^', the powers of |v⃗|_B^⊥ 2 and (|q⃗|/m_N)^2 are determined not only by the operators, but also by the response type k involved; thus, they should be written as α_k,ij and β_k,ij. As evident in Table <ref>, however, the values of α_ij and β_ij are the same for fixed i and j across all k (with nonzero ℛ_k,ij) under consideration. Hence, we may drop the k index for the subset of response types that are relevant for cosmology.]
We factorize the nuclear response function in a similar manner to obtain

    W_B,k^ττ^' = 2S_B+1/4π𝒲_B,k^ττ' e^-2y
      = 2S_B+1/4π𝒲_B,k^ττ' e^-2x v^2 (μ_χ Ba_B)^2 ,

where the values of the numerical factors 𝒲_B,k^ττ' are listed in Table <ref>.
This expression is valid for hydrogen and helium,[For arbitrary nuclei, the expression for the nuclear response function is a polynomial in y multiplied by e^-2y. However, the nuclear response is constant for the proton, and the polynomial in y is simply a constant for helium.] noting that the velocity dependence in the exponential is removed for hydrogen when setting a_p=0.
We may now recast Eq. (<ref>) as

    _B = 1/m_v^4∑_i,j∑_τ, τ^'∑_k
      c_i^τ c_j^τ^'ℛ_k,ij𝒲_B,k^ττ'(2μ_χ B/m_N)^2β_ij v^2(α_ij+β_ij)
      (1-x)^α_ij x^β_ij e^-2x v^2 (μ_χ Ba_B)^2 .

It is possible to have nonzero terms for which i≠j, indicating interference between operators; specifically, as seen in Table <ref>, there is interference between 𝒪_4 and 𝒪_6 and between 𝒪_12 and 𝒪_15.
However, when handling a single operator at a time, the interference between operators plays no role.
There is constructive interference between 𝒪_4 and 𝒪_6; hence, we expect the upper limits on the coupling coefficients in Sec. <ref> to be conservative: if another operator contributes to a signal, then the upper limit can, in principle, be made stronger.
For 𝒪_12 and 𝒪_15, there is destructive interference, and we do not expect to obtain the most conservative upper limits on the coupling coefficients for the individual analyses of these two operators.




 §.§ Cross Sections



The differential cross section, written as a function of x, for DM scattering with a baryon B is

    dσ_B/dx = μ_χ B^2/π_B  .

In Sec. <ref>, we express our cosmological constraints in terms of the total cross section σ_B(v) for scattering with protons, found by integrating Eq. (<ref>) over x:

    σ_B(v)
      = ∫_0^1 dσ_B/dx  dx
         = ∑_i,jσ_B^(ij) v^2(α_ij+β_ij)_1F_1 (1+β_ij, 2+α_ij+β_ij;-2μ_χ B^2 a_B^2 v^2),

where

    σ_B^(ij)≡μ_χ B^2/π m_v^4(2μ_χ B/m_N)^2β_ij∑_ττ^'∑_k c_i^τ c_j^τ^'ℛ_k,ij𝒲_B,k^ττ^' ,

and _1F_1 is the confluent hypergeometric function of the first kind.
However, as mentioned in Sec. <ref>, the key quantity that affects cosmological observables is the momentum-transfer cross section σ_MT,B(v), which weights the differential cross section by a factor of (1-cosθ)=2x to preferentially pick out scattering processes with large momentum transfer:

    σ_MT,B(v)
      = 2 ∫_0^1 x dσ_B/dx  dx
         = ∑_i,jσ_B^(ij) v^2(α_ij+β_ij)2(1+β_ij)/2+α_ij+β_ij_1F_1 (2+β_ij, 3+α_ij+β_ij;-2μ_χ B^2 a_B^2 v^2)  .

The velocity dependence of σ_B(v) and σ_MT,B(v) may arise either from the structure of the interaction itself [i.e., from the DM response function, captured in the power-law index 2(α_ij+β_ij)] or from the composite nature of the nucleus (i.e., from the nuclear response function, captured in the _1F_1 function).
The latter is nontrivial only for helium in this study, while for hydrogen it evaluates to 1, leaving the cross section and momentum-transfer cross section with a simple power-law velocity dependence.

It is important to note the various velocity dependencies for discerning observational signatures of different interactions; as described in Sec. <ref>, the velocity dependence influences the time evolution of the momentum transfer between DM and baryon fluids, as well as the time of their thermal decoupling.
The rate of momentum transfer in turn controls the relative size of the effect of scattering on different density modes in the early Universe, producing different signatures in the CMB and the matter power spectrum. We discuss the effect of scattering on cosmological observables in detail in Sec. <ref>.

In this study, we focus on DM interactions with protons and neglect interactions with neutrons embedded in helium nuclei; we can thus replace all coupling coefficients labeled with isospin by the coupling coefficients with protons, which we simply denote as c_i, dropping the superscript “(p)” henceforth.
Since we also focus on a single operator at a time, Eqs. (<ref>), (<ref>), and (<ref>) simplify to

    σ_B^(i)(v)
         = σ_B^(i) v^n
      _1F_1 (1+β, 2+α+β;-2μ_χ B^2 a_B^2 v^2) 
    σ_MT,B^(i)(v)
         = σ_B^(i) v^n
      2(1+β)/2 + α + β_1F_1 (2+β, 3 + α + β;-2μ_χ B^2 a_B^2 v^2) 
    σ_B^(i)   = μ_χ B^2 c_i^2/π m_v^4(2μ_χ B/m_N)^2β∑_ττ^'∑_k
      ℛ_k,ij𝒲_B,k^ττ^' ,

where we eliminated the indices on α and β, with the understanding that they correspond to α_ii and β_ii, respectively.
Additionally, we define the quantity n≡2(α + β) to emphasize that the cross section has a single power-law index of the relative velocity for a given operator.
As shown in Table <ref>, the operators cover n∈{0,2,4,6}.
However, the entire phenomenology of these different interactions does not reduce to the choice of n; namely, we find in Sec. <ref> that different operators with the same n but a different set of (α,β) yield different mass dependencies for the cosmological constraints on the coupling coefficients.




§ COSMOLOGICAL SETTING



In this section, we embed the EFT formalism into a cosmological setting.
We start in Sec. <ref> by reviewing the modified system of Boltzmann equations that describe the evolution of cosmological perturbations in the presence of DM–proton interactions.
We make the necessary modifications to all relevant background quantities, such as the DM temperature and the heat and momentum transfer rates between the DM and baryon fluids.
In Sec. <ref>, using our implementation of these equations in the code  <cit.>, we illustrate the effect of interactions on cosmological observables.




 §.§ Boltzmann equations



We work in Fourier space to express the evolution of the DM and baryon density fluctuations, δ_χ and δ_b, and their velocity divergences, θ_χ and θ_b, respectively, as

    δ̇_χ   = -θ_χ - ḣ/2   θ̇_χ   = -ȧ/aθ_χ + c_χ^2 k^2 δ_χ
      + R_χ (θ_b - θ_χ)
      
    δ̇_b       = -θ_b    - ḣ/2   θ̇_b       = -ȧ/aθ_b   + c_b^2 k^2 δ_b
      + R_γ (θ_γ - θ_b)
      + ρ_χ/ρ_b R_χ (θ_χ - θ_b)

in the synchronous gauge, where k is the wave number of a given Fourier mode (not to be confused with the index denoting DM response function types of Sec. <ref>); a is the scale factor; h is the trace of the scalar metric perturbation <cit.>; c_b and c_χ are the speeds of sound in the two fluids <cit.>; and ρ_b and ρ_χ are their respective energy densities.
The overdot represents a derivative with respect to conformal time, and the subscript γ pertains to photons.
The rate coefficients R_γ and R_χ arise from scattering processes that change particle momenta, resulting in a drag force that affects the evolution of θ_b and θ_χ.
The standard term R_γ is associated with Compton scattering <cit.>; the term R_χ arises from DM–baryon scattering <cit.>.

In a single collision, the momentum of the DM particle changes by <cit.>

    |Δp⃗_χ| = m_χ m_B/m_χ + m_B |v⃗_χ - v⃗_B|
      (n̂ - v⃗_χ - v⃗_B/|v⃗_χ - v⃗_B|)  ,

where m_B is the mass of baryon B, n̂ is the direction of the scattered DM particle in the center-of-mass frame, v⃗_χ is the DM velocity, and v⃗_B is the velocity of baryon B.
The relative velocity is v⃗=v⃗_χ - v⃗_B.
In line with Sec. <ref>, quantities labeled with “B” indicate a particular baryonic species and those labeled with “b” refer to the baryon fluid as a whole.
We assume phase space distribution functions of the form

    f_χ(v⃗_χ)    = n_χ/(2π)^3/2v̅_χ^3/2exp[-(v⃗_χ-V⃗_χ)^2/(2v̅_χ^2)] 
    
      f_B(v⃗_B)    = n_B/(2π)^3/2v̅_B^3/2exp[-(v⃗_B-V⃗_b)^2/(2v̅_B^2)]  ,

where V⃗_χ and V⃗_b are the peculiar velocities of the DM and baryon fluids, and v̅_χ^2=T_χ / m_χ and v̅_B^2=T_b / m_B are the DM and baryon velocity dispersions, respectively.
The resulting drag force per unit mass, or drag acceleration, on the DM fluid is given by

    dV⃗_χ/dt   = -1/m_χ∑_B ∫ d^3v_χ d^3v_B f_χ(v⃗_χ) f_B(v⃗_B)
      ∫dσ/dΩ |v⃗_χ - v⃗_B|
      |Δp⃗_χ| 
       = -∑_B ρ_B/m_χ+m_B1/(2π)^3/21/(v̅_B^2 + v̅_χ^2)^3/2∫ d^3v v⃗[σ_MT,B^(i)(v) v]
      exp{-[v⃗-(V⃗_χ-V⃗_b)]^2/2(v̅_B^2 + v̅_χ^2)} .

In the limit where the DM–baryon relative bulk velocity is small compared to thermal velocity,[The condition (V⃗_χ - V⃗_b)^2 ≪ (v̅_χ^2 + v̅_B^2) is satisfied for interactions that couple DM to protons at early times only; in all the cases we consider, DM–proton decoupling occurs at z>10^4, as shown in Figure <ref>.] the final expression for the rate coefficient at a given redshift is

    R_χ = aρ_b ∑_B
         Y_B/m_χ+m_Bσ_B^(i)2(1+β)/2+α+β𝒩_n
      (T_b/m_B+T_χ/m_χ)^1/2+α + β
       ×[1+(2μ_χ B a_B)^2 (T_b/m_B+T_χ/m_χ)]^-(2+β) ,

where Y_B is the mass fraction of the baryon B; 𝒩_n≡2^(5+n)/2Γ(3+n/2) / (3√(π)) is a numerical factor; and T_b and T_χ are the baryon and DM temperatures, respectively.

Finally, since we are interested in sub-GeV DM, we cannot neglect terms with T_χ in the above equations, as was the approach in Ref. <cit.>.
We thus track the coupled evolution of the DM and baryon temperatures, given by

    Ṫ_χ   = -2ȧ/aT_χ
      + 2 R^'_χ (T_b - T_χ) 
    Ṫ_b    = -2ȧ/aT_b
      + 2μ_b/m_χρ_χ/ρ_bR^'_χ (T_χ - T_b)
      + 2μ_b/m_e R_γ (T_γ - T_b)  ,

where m_e is the electron mass; μ_b≈m_H (n_H+4n_He)/(n_H+n_He+n_e) is the mean molecular weight of the baryons; n_H and n_He are the number densities of protons and helium nuclei, respectively; and the heat-exchange rate coefficient R_χ^' is given by Eq. (<ref>), but with an additional factor of m_χ/(m_χ + m_B) multiplying each summand.




 §.§ Effect of Scattering on Cosmological Observables







We have modified the linear Boltzmann solver  to implement the Boltzmann equations given in Sec. <ref>, incorporating the effect of DM–proton scattering on the evolution of cosmological perturbations, for all operators in Eq. (<ref>).
Figure <ref> shows the relevant background quantities, T_χ and R_χ, as functions of redshift for a subset of operators that correspond to interactions whose cross sections scale with different powers of velocity (labeled as v^n in the legends of the figure).
To illustrate the evolution of these quantities, we fix S_χ=1/2 and m_χ=1, set the coupling coefficients to their respective 2σ upper limits reported in Table <ref>, and keep other cosmological parameters at their best-fit Planck 2015 values <cit.>.
A stronger velocity dependence (larger n) leads to more momentum transfer at early times and to a later thermal decoupling time of the DM and baryon fluids.
This difference in the evolution of R_χ determines the relative size of the effect of DM–baryon scattering on different density perturbation modes, traceable through cosmological observables such as the CMB power spectra C_ℓ and the three-dimensional matter power spectrum P(k).

In Figure <ref>, we illustrate the effect of scattering on these observables for 𝒪_1 (with n=0) and 𝒪_8 (with n=2) by comparing to the standard “CDM” case with no DM–proton interactions.
For DM–proton scattering, we fix S_χ=1/2 and set the coupling coefficients to their 2σ upper limits reported in Table <ref>, keeping all other cosmological parameters at their best-fit Planck 2015 values, as in Figure <ref>.

For the left panel of Figure <ref>, we set m_χ=1.
To get a sense for the range of ℓ multipoles in which the scattering signal is most prominent, we show the percent difference between the “CDM” case and the cases with DM–proton scattering.
As was previously noted in Refs. <cit.>, on large angular scales in the CMB, the coupling of DM to baryons presents similarly to baryon loading: the tight coupling with DM particles at early times effectively increases the total mass of the baryons, enhancing power at low multipoles (c.f. negative percent differences).
On small scales, the drag force between baryons and DM dissipates the momentum of the baryon–photon fluid, damping baryon acoustic oscillations and suppressing power more for modes that enter the cosmological horizon earlier (roughly corresponding to larger multipoles).
Furthermore, acoustic peaks shift to smaller angular scales as a result of the decrease in the speed of sound in primordial plasma.
Since the EE polarization spectra display sharper oscillatory features in multipole space compared to TT, the shift produces higher-amplitude oscillations seen in this difference plot.

In the right panel of Figure <ref>, we show the matter power spectrum today, comparing the “CDM” case to the cases with DM–proton scattering for m_χ=1 and m_χ=1.
When DM couples to protons, P(k) exhibits oscillatory suppressions at large k due to DM tightly following the behavior of the baryon–photon fluid and undergoing “dark oscillations” at early times, when the corresponding modes enter the horizon.
For heavier DM particles, the oscillations are shifted toward larger k, indicating earlier thermal decoupling from baryons; at a fixed DM energy density, increasing the DM mass results in a lower particle number density and thus a reduced interaction rate to maintain thermal equilibrium.

With the ability to compute the power spectra for a given cosmology, we proceed to search for signals consistent with DM–baryon interactions within the data from the Planck satellite.








§ NUMERICAL RESULTS



We perform a joint likelihood analysis of Planck 2015 temperature, polarization, and lensing anisotropy measurements, using the Planck Likelihood Code v2.0 () <cit.> to search for evidence of DM–proton scattering in the early Universe.
For high multipoles,  we use the nuisance-marginalized joint TT, TE, EE likelihood ( lite), with ℓ=30–2508 for TT and ℓ=30–1996 for TE and EE.
Using  lite over the full likelihood requires substantially fewer computational resources, and we find no significant effect on the derived values of cosmological parameters and constraints.[Planck high-multipole polarization could have systematic issues <cit.> that may affect cosmological parameter estimation; however, excluding ℓ≥30 polarization would degrade our reported constraints on DM interactions by less than 8%.]
The lensing likelihood we use contains SMICA map-based lensing reconstruction for multipoles in the range ℓ=40–400.

To sample the cosmological parameter space, we employ  <cit.> with the  <cit.> implementation of the nested sampling algorithm <cit.>.
We have verified that our pipeline recovers the published ΛCDM parameter values and constraints <cit.> to within 0.14σ, in the limit of vanishing DM coupling coefficients.
Since the convergence criteria for nested likelihood sampling are not clearly defined (in contrast to the case of MCMC methods), we check for convergence in our samples by varying the number of live points passed to  and ensuring that the inferred constraints do not change by more than an order of percent in numerical value between different sampling runs.
We assume that baryons, photons, and DM are tightly coupled and in thermal equilibrium at the start of the integration of the Boltzmann equations; this condition is satisfied for all interactions at their 1σ and 2σ exclusion limits in Table <ref>.
Additionally, we consider only a flat geometry.

We perform a likelihood analysis for each individual interaction operator, repeating the fitting procedure for the following fixed values of DM mass: m_χ∈{15, 1, 1, 1 }.
Thus, in addition to the six standard ΛCDM parameters (baryon density Ω_bh^2, DM density Ω_χ h^2, Hubble parameter h, reionization optical depth τ, amplitude of the scalar perturbations A_s, and scalar spectral index n_s), each sampling run also includes the coupling coefficient c_i for 𝒪_i as a free fitting parameter; we assume wide flat priors for all parameters.
We obtain the 68% and 95% confidence-level exclusion limits for the coupling coefficients c_i from their posterior probability distributions, and we convert these limits to ones for the corresponding interaction cross section using Eq. (<ref>).





We find no evidence of DM–proton scattering for m_χ>15.
We thus report a complete set of cosmological upper limits on the coupling coefficients and cross sections associated with the 14 different operators in Eq. (<ref>).
Our inferred upper limits on the cross sections for proton scattering, σ_p^(i), with spin-1/2 DM are listed in Table <ref>.[The lowest-mass data points for 𝒪_15 are missing. As we alluded to in a previous footnote, our code is not sufficiently optimized to produce accurate power spectra for the extremely strong velocity dependence of n=6 at low DM mass. We leave a detailed treatment of this regime for future work.]
These cross sections are evaluated at v=(220 /)/c, the relative velocity relevant for direct detection.
In the Appendix, Tables <ref> and <ref> list the upper limits in terms of the quantities σ_p^(i) and c_i^2, respectively, to allow for a more straightforward comparison to other works.
Note that since 𝒪_1 is velocity independent, σ_p^(1)=σ_p^(1), consistent with the notation in our previous study <cit.>.
The limits for DM spins other than 1/2 are easily obtained by rescaling the cross section by S_χ (S_χ + 1) / (3/4) (in accordance with Table <ref>) for any operator except 𝒪_1, 𝒪_3, 𝒪_7, and 𝒪_10, which are independent of DM spin.

We choose the subset of operators {𝒪_1, 𝒪_5, 𝒪_8, 𝒪_15} to represent the various types of velocity-dependent interactions with n∈{0,2,4,6}, and we choose 𝒪_4 as a representative of spin-dependent interactions.
For these operators, we sample the likelihood on a finer grid of logarithmically-spaced masses: m_χ={15, 32, 1, 32, 1, 32, 1 }.
We show the resulting 95% confidence-level exclusion curves in Figure <ref>.
Finally, in Figures <ref>, <ref>, <ref>, and <ref>, we show the posterior probability distributions of ΛCDM parameters and coupling coefficients for 𝒪_1, 𝒪_8, 𝒪_5, and 𝒪_15, respectively, for fixed m_χ=1 and S_χ=1/2.








§ DISCUSSION AND CONCLUSIONS



We analyzed Planck 2015 TT, EE, TE, and lensing power spectra in search for evidence of DM–proton interactions in the early Universe.
Our results are consistent with absence of interactions through operators dependent or independent of spin and velocity, and we thus report the first cosmological upper limits on the full nonrelativistic effective field theory of DM–proton scattering for masses above 15 .
The main result is summarized in Figure <ref>.
In this section, we discuss in more detail our results and their implications for future studies.

We first note that comparing the cross-section constraints for 𝒪_1 and 𝒪_4 (shown in the right panel of Figure <ref>) illustrates the impact of including scattering with helium in our analyses.
These two operators represent the standard spin-independent and spin-dependent scattering with no dependence on relative velocity, and the main difference in their upper limits comes from the fact that DM cannot interact through 𝒪_4 with helium, which has zero nuclear spin, while 𝒪_1 interacts with both helium and hydrogen.
The coefficient R_χ, defined in Eq. (<ref>), quantifies the rate of momentum transfer between DM and all baryonic species that participate in a given interaction.
If there is only scattering with protons (as is the case for 𝒪_4), the maximum of R_χ occurs for DM masses near a GeV, above which R_χ rapidly decreases (asymptoting to ∼1/m_χ behavior at m_χ≫1 ), resulting in a loss of constraining power.
Including scattering on helium (for 𝒪_1) shifts its maximum above a GeV (since helium is roughly four times as massive as the proton), thereby substantially improving the constraint at high mass.
As previously reported in Ref. <cit.>, although helium contributes only about a quarter of the total baryonic mass, the modest shift of the maximum in R_χ amounts to an improvement of about a factor of 6 near 1 in sensitivity of CMB probes to spin-independent scattering, while it has no impact on constraints for spin-dependent scattering.

We now examine the constraints on the coupling coefficients in the left panel of Figure <ref> for operators whose cross sections have different scalings of relative velocity.
There is a hierarchy of constraints: interactions with the stronger velocity dependencies are more suppressed in the regime of very low velocities (of interest to cosmological studies), so CMB observables have less constraining power on the coupling coefficients of such interactions.
It is also evident that CMB observations are most sensitive to DM masses around a , when considering interactions with protons.
Roughly speaking, the sensitivity drops at higher masses, because the amount of momentum transferred in the scattering process saturates (μ_χ B approaches the mass of the proton) but the drag force per unit mass drops as 1/m_χ; the sensitivity drops at lower masses, because the momentum transfer scales as m_χ in that regime (μ_χ B approaches m_χ).
Understanding the detailed mass dependence in either panel of Figure <ref> for a given operator is nontrivial and does not immediately follow from the analytic expressions of Secs. <ref> and <ref>.
One reason is that the CMB observables are controlled by two quantities: the coefficients for the rate of momentum transfer R_χ and of heat transfer R_χ^'.
The former appears in Eq. (<ref>) and controls the impact DM interactions have on density and velocity fluctuations at a given redshift; the latter appears in Eq. (<ref>) and controls the time of thermal decoupling between the DM and baryon fluids.
Their scaling with m_χ is different, and the resulting mass dependence observed in Figure <ref> follows from a combination of the two.

When comparing the left and right panels of Figure <ref> or examining the entries of Table <ref>, it is important to keep in mind that we treat the coupling coefficients as the primary parameter to constrain from the data, for a fixed DM mass.
On the other hand, the scattering cross section is a derived quantity, obtained by the relation in Eq. (<ref>); it scales quadratically with the coupling coefficient, but its dependence on m_χ is determined by the operator at hand.
We emphasize that the low-mass limiting behavior of the exclusion curve for the cross section is not solely determined by n=2(α+β), the power of the cross section scaling with relative velocity.
From Eq. (<ref>), we see that velocity dependence in the cross section arises from the perpendicular component of the relative velocity (associated with the power α) and the momentum transfer (associated with the power β) that appear in the DM response function; however, the momentum transfer is a function of μ_χ B, which introduces a mass scaling for the cross section that involves β but not α.
This dependence on β explains the different low-mass behaviors of the exclusion curves associated with the coupling coefficients of different operators, shown in the left panel of Figure <ref>.

In the high-mass limit of m_χ≫1, the situation is greatly simplified.
The coefficient for the rate of heat transfer R_χ^' approaches that for the rate of momentum transfer R_χ, and the thermal term T_χ / m_χ for DM in Eq. (<ref>) is negligible, regardless of operator.
Thus, cosmological observables are controlled solely by σ_B^(i) / m_χ.
Moreover, σ_B^(i) in Eq. (<ref>) scales as c_i^2 and has no m_χ dependence.
All exclusion curves in Figure <ref> exhibit the behavior σ_B^(i)∼c_i^2∼m_χ.
This relation allows extrapolation of our results in Table <ref> to arbitrarily high masses: the upper limit on the DM–proton cross section for m_χ> is obtained by scaling the upper limit at 1 by m_χ/.

Let us now take a look at the possible degeneracies between the DM scattering signal and other cosmological parameters.
The addition of high-ℓ polarization and lensing improves our reported constraint by about 30%; thus, most of Planck's constraining power comes from the temperature anisotropy at the smallest angular scales measured with a high signal-to-noise ratio (where the Planck error bars are their smallest in the left panel of Figure <ref>).
In this regime, where TT measurements at ℓ∼1200 dominate the constraint, the effect of DM–proton scattering is not strongly degenerate with any of the standard cosmological parameters.
There is, however, a mild degeneracy with the scalar spectral index n_s, as seen in Figure <ref> for 𝒪_1.
The origin of this degeneracy can be understood from the left panel of Figure <ref>: increasing the strength of the DM–proton coupling progressively suppresses power on smaller scales, and a larger value of n_s is needed to restore C_ℓ to within the measurement error, giving rise to the mild positive correlation between c_1 and n_s.
Analogous plots of the posterior probabilities for the cases of velocity-dependent scattering are shown in Figures <ref>, <ref>, and <ref> for 𝒪_8, 𝒪_5, and 𝒪_15, respectively.
The degeneracy with n_s reverses sign when the leading effect of scattering on relevant angular scales (ℓ∼1200) gives rise to an enhancement of power rather than a suppression (such as the case for scattering through 𝒪_8, for example, shown in Figure <ref>).

Finally, we briefly discuss future prospects for cosmological probes of DM interactions.
Given that the effect of DM–baryon interactions is stronger on smaller angular scales, high-resolution ground-based CMB observations from existing experiments (such as the Atacama Cosmology Telescope <cit.> and the South Pole Telescope <cit.>) and from future experiments (such as the Simons Observatory <cit.> and CMB-Stage 4 <cit.>) could improve upon our limits, particularly in the regime of sub-GeV DM.
As shown in the right panel of Figure <ref>, DM scattering is also imprinted on P(k) and is progressively more prominent at larger values of k.
At the level of current constraints, DM scattering becomes inefficient at redshifts z<10^4 due to Hubble expansion, so its effect on P(K) is similar to a k-dependent change in the initial conditions, for all effective operators considered in this work.
We thus expect that galaxy-survey and large-scale structure (LSS) measurements from BOSS <cit.>, DES <cit.>, LSST <cit.>, DESI <cit.>, etc. can help improve constraints on DM interactions; see, for example, results from the Lyman-α data analysis in Ref. <cit.>.
Moreover, there is evidence that including the Lyman-α forest power spectrum leads to smaller values of the scalar spectral index (when also allowing the spectral index to run) <cit.>, indicating that Lyman-α data may aid in constraining DM interactions on multiple fronts.
However, unlike the case of CMB, LSS analyses require modeling of complicated baryonic effects as well as treatments of nonlinearities and other systematic effects that arise on small scales; a calculation of the nonlinear effects depends on the cosmological context and has not yet been fully addressed within a cosmology with DM–baryon scattering, in which P(k) exhibits oscillatory features shown in Figure <ref>.
For these reasons, the CMB is currently the most robust cosmological probe of DM–baryon interactions; we have thus only focused on CMB measurements and leave the extension to LSS analyses for future work.

In summary, this study paves a path for comprehensive and model-independent cosmological studies of low-energy DM physics.
Cosmological probes are particularly attuned to testing for the presence of sub-GeV dark matter, whose interactions may not be detectable in standard laboratory searches for weakly-interacting massive particles, but may leave an imprint in the early Universe.
Many forthcoming observations of the CMB and large-scale structure promise to further deepen the coverage of the relevant parameter space for DM candidates.

























VG gratefully acknowledges the support of the Eric Schmidt fellowship at the Institute for Advanced Study and the hospitality and support of Juna Kollmeier through the Visitor Program at Carnegie Observatories, where this work started.
The authors thank John Beacom, Joanna Dunkley, Marc Kamionkowski, and Timothy Morton for useful discussions.
Posterior probability distributions in this study were visualized using  <cit.>.




§ APPENDIX


In this Appendix, we present the 68% and 95% confidence-level upper limits on two quantities related to the scattering cross section σ_p^(i) of Table <ref>, for a more straightforward comparison with other literature.
Table <ref> lists the upper limits on σ_p^(i), corresponding to quantity σ_0 in Ref. <cit.>, for example.
Table <ref> lists the upper limits on c_i^2, which corresponds to the limits in Figure 6 of Ref. <cit.>, for example.







