Matter Power Spectrum in Hidden Neutrino Interacting Dark Matter Models: A Closer Look at the Collision Term

Dark Matter (DM) models providing possible alternative solutions to the small- scale crisis of standard cosmology are nowadays of growing interest. We consider DM interacting with light hidden fermions via well motivated fundamental operators showing the resultant matter power spectrum is suppressed on subgalactic scales within a plausible parameter region. Our basic description of the evolution of cosmological perturbations relies on a fully consistent first principles derivation of a perturbed Fokker-Planck type equation, generalizing existing literature. The cosmological perturbation of the Fokker-Planck equation is presented for the first time in two different gauges, where the results transform into each other according to the rules of gauge transformation. Furthermore, our focus lies on a derivation of a broadly applicable and easily computable collision term showing important phenomenological differences to other existing approximations. As one of the main results and concerning the small-scale crisis, we show the equal importance of vector and scalar boson mediated interactions between DM and light fermions.


Introduction
Precise measurements of cosmic microwave background (CMB) anisotropies have been building strong evidence for the existence of a new form of matter, called dark matter (DM) [1,2]. However, its nature has not yet been uncovered and one of the most important subjects both in astrophysics and in particle physics. Recently vigorous efforts have been devoted to cosmological probes of interaction strengths between the DM and other long-lived particles [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Interestingly, such probes are not limited to couplings of the DM to standard model (SM) particles (e.g., baryons, photons, and neutrinos). Couplings to hidden particles are equally subject to searches. In this paper, we restrict our discussion to hidden light particles, which we call neutrinos for simplicity. However, the formulation developed and given in this paper is applicable to other models with DM couplings to SM particles.
Interacting DM models are not only within the scope of precise measurements of the large-scale structure of the Universe. They also have their motivation in apparent discrepancies between predictions from DM-only N -body simulations and observations on subgalactic scales. Such discrepancies are called the small-scale crisis collectively: the missing satellite problem [17,18]; the cusp vs core problem [19][20][21]; the too big to fail problem [22,23]. The simulations assume that the DM consists of particles with negligible thermal velocities and faint interactions, called cold dark matter (CDM). The small-scale crisis may imply alternatives to the CDM paradigm, while it has to be clarified by state-of-the-art hydrodynamic simulations what role baryonic processes play in the formation and evolution of subgalactic objects [24,25]. One famous alternative is called the warm dark matter (WDM) model, in -1 -
Although there is a growing interest in interacting DM models, it is still unclear what the evolution equations of cosmological perturbations are in such models. This is because it is difficult to handle and simplify collision terms of the Boltzmann equations. Some works start with the relativistic Navier-Stokes equation for the DM imperfect fluid in particle flow manifest (Eckart's) formulation [37,38]. They determine fluid variables with the help of the Chapman-Enskog method to estimate damping scales in matter power spectra in interacting DM models. Others just add a collision term in the evolution equations of cosmological perturbations by analogy to the well-known Thomson scattering term for baryons and photons [39]. One plausible way is to reduce the collision term to the Fokker-Planck equation by assuming the momentum transfer in each collision is smaller than the typical DM momentum. Such formulation is developed for the traditional bino-like DM in [40]. However, the overall factor, i.e., the reaction rate, in the Fokker-Planck equation is controversial so far. A systematic expansion of the collision term in terms of the momentum transfer leads to a reaction rate proportional to the invariant amplitude at a zero momentum transfer t → 0 [41][42][43]. On the other hand, in [44,45], the reaction rate is given by t-averaging like dt(−t)dσ/dt.
The two formulations introduced above result in different phenomenological consequences. We consider a simple model, where the SM sector is extended by a Dirac DM, a Dirac (hidden) neutrino, and a mediator. A similar scenario is investigated in [28]. When the mediator is a scalar, the reaction rate with a zero momentum transfer is negligible and does not change the matter power spectra on and above subgalactic scales within a plausible range of model parameters. A subgalactic damping scale can be achieved by a vector mediator within this formulation [29]. On the other hand, both vector and scalar mediators can suppress the resultant matter power spectra with the t-averaged reaction rate. We address this point by calculating the resultant matter power spectra in the latter formulation numerically. To this end, we derive the evolution equations of cosmological perturbations in two gauges: the conformal Newtonian gauge and the synchronous gauge [46]. We provide an explicit form of gauge transformations between them. We also show a derivation of the t-averaged reaction rate. It may be useful because a similar derivation is given only in an unpublished thesis [44]. In the recent ETHOS (effective theory of structure formation) papers [34,35], they study the structure formation in interacting DM models based on the former treatment of the collision term, while they also mention the importance of the t-averaging in some models. The ETHOS paper and this paper are complementary to each other in a treatment of the collision term.
The paper is organized as follows. In section 2 we start from first principles and give a detailed derivation of the Fokker-Planck equation with the t-averaged reaction rate. Furthermore, the evolution equations of cosmological perturbations in the synchronous gauge are derived for the most general case of an imperfect fluid. We show explicitly in appendix A that our results transform into the form of the conformal Newtonian gauge according to the rules of gauge transformation. In section 3, we give an introduction of the neutrino interacting DM model first. Then, we summarize our analytic results for scalar and vector mediators: the relic density of the DM; the t-averaged reaction rate; the resultant smallest mass of halos. In appendix C, we present details of our calculations of the chemical decoupling and also summarize the results for models with pseudo scalar and pseudo vector mediators. Finally, we show the resultant matter power spectra to stress that not only a vector, but also a -2 -JCAP11(2016)043 scalar mediator can lead to a sizable suppression of matter power spectra. In appendix B, we provide an estimation of a critical wavenumber below which the perfect fluid approximation appears valid.
Throughout this paper, we use the Planck 2013 cosmological parameters [2]: Ω m = 0.3175, H 0 = 67.11, ln(10 10 A s ) = 3.098, and n s = 0.9624. Updating these input parameters to the Planck 2015 ones would not change our results significantly.

Fokker-Planck equation
In this section, the perturbed Fokker-Planck equation is derived. Our starting point is the classical Boltzmann equation with the collision term. We expand it assuming the momentum transfer per collision is smaller than the typical DM momentum. Within this approximation the collision term satisfies detailed balance and respects number conservation. As a further result of this expansion method, the momentum transfer rate can easily be computed by first taking a t-average and secondly a thermal average of the differential scattering cross section. As an important result of the formulation used, the t-average is a direct consequence of the expansion method. Other methods like in [42] expand the scattering amplitude in terms of a small momentum transfer and keep only the zero order. But this approximation shows a completely different phenomenology for certain DM theories as will be shown as an explicit example in section 3. As part of this section, we develop the evolution equations of linear cosmological perturbations in the synchronous gauge. A comparison to previous works is given. The results are equivalent to those in the conformal Newtonian gauge under the gauge transformation law as we show for the first time in appendix A.

Collision term
In this subsection, we derive a Fokker-Planck equation for the DM phase space distribution function f , partially inspired by the unpublished thesis [44]. Our starting point is the classical Boltzmann equation for the DM, where P µ is the conjugate momentum of the spatial coordinate x µ . When we handle the collision term C[f ], it is convenient to take a local inertial frame X µ , where the metric is diag(−1, +1, +1, +1) and the proper momentum is denoted by p µ = (E, p). We normalize the distribution function such that s d 3 p/(2π) 3 (p µ /E)f = n µ , where s are spin degrees of freedom and n µ is the DM number current. If we assume the DM particles to interact elastically with particles in a thermal bath, i.e., DM(1)+TP(2) ↔ DM(3)+TP(4), the collision term takes the form of
If the elastic scattering is T -inversion invariant, |M| 2 's are identical between the forward and backward scatterings, (2.6) In the presence of four-momentum conservation δ 4 (p 1 + p 2 − p 3 − p 4 ), thermal distributions satisfy From (2.6) and (2.7), we obtain the following relation: Thus, the collision term is We can easily check that the above expression satisfies the so-called detailed balance, i.e., C[f 1 ] = 0 if f 1 = f eq 1 and f 3 = f eq 3 , which follows from the T -inversion invariance. We assume that the momentum transferq = p 3 − p 1 is smaller than the typical DM momentum p 1i and expand the collision term up to the second order, with the velocity of the particle v = p/E. After collecting terms, we obtain The collision term is 14) It should be noted that α i = 0 and thus C[f 1 ] = 0 if f 1 = f eq 1 , which implies that the detailed balance is maintained in the approximation.
We expand S eq (p 3 , p 1 ) in terms ofq, noting that p µ 3 = ( E 2 1 + 2p 1iqi +q 2 , p 1 +q). The scalar function S eq depends onq only through p 1iqi andq 2 . Since we keep the terms only up to the second order in terms ofq, the expansion in terms ofq 2 leads to higher order terms in C[f 1 ], which are to be neglected in our treatment. Therefore we expand S eq only in terms of p 1iqi as follows: 1 S eq (p 3 , p 1 ) ≃ S eq 0 (p 1 ,q 2 ) + S eq 1 (p 1 ,q 2 )p 1iqi , (2.16) where S eq 0 (p 1 ,q 2 ) and S eq 1 (p 1 ,q 2 ) are the expansion coefficients defined by (2.16). Substituting (2.16) and (2.17) into (2.15) and using (2.8) and (2.10), we obtain In the second equality of (2.18), we have dropped the term proportional to S eq 0 (p 1 ,q 2 )q i , since it vanishes after integrated in terms of d 3 p 3 = d 3q . In addition, in the last equality, we have replacedq iqj with (1/3)δ ijq 2 for the same reason. It should be noted that, γ ij =γ ij holds up to the second order inq.
In practice it is more convenient when evaluating (2.19) to replace the perturbative quantities S eq 0 (p 1 ,q 2 ) and −q 2 with their non-perturbative counterparts S eq (p 3 , p 1 ) and respectively. These resulting coefficients only differ through higher order terms and amount to an alternate perturbative expansion. Substituting (2.18) and (2.19) into (2.14), we obtain a Fokker-Planck-type equation for f 1 since the collision term becomes

JCAP11(2016)043
where the momentum transfer rate is Here, dσ/dt is the differential cross section and v is the relative velocity of the initial particles. The center of mass momentum is evaluated by and m is the mass of the particle. This equation satisfies two important requirements. First, it maintains the detailed balance: if f 1 = f eq 1 , then C[f 1 ] = 0. Second, it conserves the DM number, If the DM particles decouple from the thermal bath when they are relativistic, momentum transfer in each collision is as large as the typical momentum of the DM, which may spoil our approximation approach, i.e., the Fokker-Planck equation. It may be useful to give the non-relativistic limit. Then, the Fokker-Planck equation is where the momentum transfer rate is This expression is the same as given in [44,45]. The cross section is essentially independent of p 1 since we consider the case that DM is non-relativistic, or in other words, |p 1 | is much smaller than m 1 . We focus on such a case in the following. Before closing this subsection, let us discuss the relation between [41][42][43] and the present paper. The main difference is the presence of t-averaging in the momentum transfer rate of (2.24). Once we set t → 0 in dσ/dt, we can evaluate the t-integral analytically to reproduce the result in [41][42][43]. The t-averaging originates from the approximation in (2.16) and (2.17) and the replacement of the perturbative quantities with the non-perturbative counterparts after (2.19). In this respect, our formulation is not a systematic expansion in terms of the momentum transfer like that in [42]. However, in some cases, the expansion of invariant amplitudes is not a good approximation since the leading order does not give the dominant contribution. One such example is the scalar operator of DM-neutrino interaction investigated in the present paper. There, the leading order is suppressed by a factor of m 2 ν /(−t) when compared to the next-to-leading order.

Perturbation theory in the synchronous gauge
Now we develop a linear theory in the synchronous gauge: Up to the first order of cosmological perturbations, the Fokker-Planck equation is given bẏ with γ = γ 0 (τ ) + γ 1 (x) and the comoving momentum q = ap. 2 The homogeneous and isotropic part, i.e., the leading order iṡ A solution, is parametrized by the DM temperature T χ0 (τ ) and the DM number density per spin degree of freedomn/g χ with g χ = 2s χ + 1. Its evolution is described by The DM temperature is tightly coupled to the temperature of the thermal bath T χ0 = T 0 ∝ 1/a before the kinetic decoupling γ/H > 1. After they decouple, the DM particles start to stream freely and the temperature decreases adiabatically T χ0 ∝ 1/a 2 . The first order perturbation follows: Here, we define the Fokker-Planck operator by In the Fourier space k i = kk i , these equations are rewritten aṡ Hereafter we consider only the scalar perturbations, defining θ TP , η, and h such that θ TP = ik i u i and h ij =k ikj h + k ikj − 1 3 δ ij 6η (the same notation as in [46]). In order to handle the Fokker-Planck operator, we expand f 1 in terms of eigenfunctions of the Fokker-Planck operator, Here Y ℓ m and L α n denote the spherical harmonics and the Laguerre polynomials, respectively. Noting the rotational symmetry, we can write

JCAP11(2016)043
with the Legendre polynomial P ℓ , and vice versa, After a lengthy but straightforward calculation, we obtain the Boltzmann hierarchy: Here we introduce three new quantities: The first quantity is essentially proportional to the Hubble expansion rate: R = aH/2. Only a few of the second and third quantities are non-zero before the kinetic decoupling (T χ0 = T 0 ): A 0 = 1 and B 1 = 1, while the others vanish. Higher orders of the second quantity become non-zero after the kinetic decoupling (T χ0 ≪ T 0 ): Although we need to solve the full Boltzmann hierarchy to obtain a rigorous result, just taking some small moments of n and ℓ can give the fluid approximation (see discussion in subsection 3.3). The perturbations f nℓ with small n and ℓ can be interpreted as primitive variables of the DM imperfect fluid (i.e., mass density ρ, bulk velocity potential θ, pressure P , and anisotropic inertia σ): The dynamics of the DM imperfect fluid is described by the following equations: The pressure perturbation δP can be decomposed into isentropic c 2 χ δ and entropy π perturbations: The sound speed squared of the DM fluid is The evolution of the DM imperfect fluid can be rewritten aṡ

Neutrino interacting dark matter
The section starts with the introduction of the neutrino interacting DM model via a MeVscale boson. This particle combination leads in a valid parameter region to a possible solution to all three small-scale crisis problems if the mediator is of vector type [29]. We reproduce and confirm these results by using the method that is derived in the previous section to describe the DM kinetic decoupling. The used method has a different expansion of the collision term when compared to the aforementioned reference and to others like [42]. Furthermore, by using this alternative description we explicitly show a suppression of the power spectrum for other types of mediators as well. The suppression is sizable enough -9 -

JCAP11(2016)043
to reduce the abundance of dwarf galaxies but unexpected from the point of view of the above literature. In particular, scalar and vector mediators share an analogue phenomenology within our model set-up and the parameter region is relatively similar concerning the minimal size of the first protohalos. Approximation methods to follow the evolution of cosmological perturbations are also given. Finally, the matter linear power spectrum for scalar and vector interactions are presented, showing a suppression of powers on subgalactic scales.

Simplified neutrino model
A simplified model extends SM by a DM fermion and additional light fermions (denoted by ν). The DM fermion and the additional light fermions are assumed to be of Dirac type, coupled by a MeV-scale boson denoted by φ. In particular, this choice allows us to write down the following set of renormalizable dimension four operators without derivatives: (3.1) Here, we assume parity conservation in the interaction Lagrangian and consider each operator type separately. There are four parameters: the DM mass m χ , the light mediator mass m φ , the DM-mediator coupling g χ , and the light fermion-mediator coupling g ν . Specifically, extensions of the simplified model (3.2) into ultraviolet complete models and their constraints have already been investigated by many authors in connection with the small-scale crisis (for an exemplary list of references, see [47][48][49][50]). For simplicity and by analogy to previous works we call the light fermions hidden neutrinos. In the early Universe, the DM and the hidden neutrinos are assumed to be in thermal equilibrium, where a temperature difference when compared to the SM sector hides the additional light fermions. Further, the light boson is in thermal equilibrium with the neutrino sector during the DM chemical freeze-out. For all operators the parameters are chosen such that the relic density of the DM is dominantly determined through χχ → φφ annihilation and not via direct s-channel neutrino production. This is because for the vector, scalar, and pseudo scalar interactions, we assume g ν ≪ g χ (see [51] for a list of possible natural explanations). In this scenario, the DM relic abundance for all operators is independent of the neutrino coupling g ν . In appendix C.1 we provide for all operators the full calculus of the annihilation cross section and the relic abundance. Due to a more complicated but less illuminating phenomenology, we discuss the results for the pseudo scalar and pseudo vector operators in appendix C.2.

Minimal halo mass
Elastic scattering via a MeV-scale boson keeps the DM for a long time in kinetic equilibrium with the hidden neutrino sector. During kinetic equilibrium, the DM density perturbations do not grow but oscillate. This phenomena is known as acoustic oscillations and has been shown in [3-15, 27-31, 33-36] to be the dominant damping mechanism of the density perturbations in the case of a late kinetic decoupling.
In the cosmological perturbation theory, the mode that enters the horizon at the kinetic decoupling defines a cutoff in the linear matter power spectrum of density fluctuations. Only the DM density modes that enter the horizon thereafter can significantly grow and collapse later into halos. Thus, fluctuations on shorter scales are damped. The minimal mass of first protohalos can be estimated by the mass inside a sphere with radius of Hubble horizon at the time of the kinetic decoupling: where the matter density ρ m and the Hubble expansion rate H are evaluated at the kinetic decoupling. Here, we allow for a different light fermion temperature from the photon temperature to hide the additional neutrinos. The ratio between the two temperatures is defined as r = T kd ν /T kd γ , where the superscript kd means the corresponding value at the DM kinetic decoupling that occurs when the momentum transfer rate γ equals to the Hubble expansion rate H.
In the following, we derive an approximation method to estimate the kinetic decoupling temperature T kd ν in order to calculate the corresponding cutoff mass according to (3.5). The general expression for γ (2.24) is adjusted to describe the scattering of the DM with the light fermions. Dividing it by the Hubble expansion rate and by introducing the following dimensionless variables x = |p ν |/T ν , y = T ν /m χ , and z = m φ /m χ , one ends up with the following form: where we multiply by the number of light fermion species N ν . The phase space distribution function f eq ν (x) is the usual equilibrium Fermi-Dirac distribution, where we neglected the mass of the light fermions: Furthermore, the dimensionless quantity g(xy, z) is defined as the t-averaged scattering amplitude squared: where in this convention, is the invariant amplitude squared that are averaged over the initial and final spin states. Equation (3.6) is the basic formula for the kinetic decoupling description of the neutrino interacting DM. In the following, we derive analytic estimates for the scalar and vector operators, which are valid in a broad range of parameters and derive their corresponding M cut scaling patterns. In the case of the pseudo scalar and pseudo vector operators, this approximation that we call the effective propagator description is only valid in a small parameter space, and thus (3.6) has to be solved numerically at some point. The results are given in Scalar operator: (3.11) In the parameter region we are interested in, it turns out that the mass of the mediator is much larger than the kinetic decoupling temperature. In this case, the Mandelstam t in the boson propagator denominator of the scattering amplitudes can be neglected. We call this approximation the effective propagator description. The propagator denominator can be simplified in such a way because t ∈ [0, −4p 2 ν ] and the neutrino momentum is further limited by the phase space distribution: |p ν | ≃ T ν . So t can be neglected in the denominator of the propagator as long as T kd ν ≪ m φ , which is the case in the parameter region of the scalar and vector operators.
Within the effective propagator framework, g(xy, z) is only a polynomial function in its variables and the integral in (3.6) has even an analytic expression. To the leading order in T ν , we find for the vector operator and for the scalar operator (3.13) with α χ/ν = g 2 χ/ν /(4π). To estimate the kinetic decoupling temperature, we set γ/H = 1 in the last two equations, which are solved for T kd ν . 3 The corresponding minimal halo masses that are derived from the kinetic decoupling temperature according to (3. where we normalize r to the SM neutrino temperature ratio: r 0 = (4/11) 1/3 . To be consistent with constraints on additional radiation components, we use the combined results of Big Bang nucleosynthesis (BBN) and CMB constraints given in [52] to derive an upper bound for our model within the 1σ error bar:  Table 1. Upper bounds on (r/r 0 ) 9/2 × (N ν /6) 3/4 derived from [52] are shown. We separate two extreme cases: the mediator is still relativistic at BBN g pol = {1, 3}; its contribution to the radiation components can be neglected (g pol = 0). The factors on the right column reduce the cutoff masses (3.14) and (3.15) by at least one order of magnitude.  Here, we consider the possibility of having a sub-MeV scale mediator contribution to the radiation components at BBN. In table 1 we summarize the upper bounds for two extreme scenario: the mediator does not contribute (g pol = 0); the mediator is still relativistic at BBN and contributes via its internal degrees of freedom (g pol = {1, 3} for the scalar and massive vector mediators, respectively).
First of all, these cutoff masses (3.14) and (3.15) have the same scaling dependence, and thus differ only by a numerical constant and depend mostly on the boson mass. Using the relic density constraint on α χ given by (C.5) and (C.6), we see that M cut is essentially independent of the DM mass. In figure 1, contour lines of a constant M cut are shown for the scalar and vector interactions in the (m φ , α ν )-plane. In order to account for the missing satellite problem and to be consistent with Ly-α forest bounds, the cutoff mass has to be roughly in between 10 7 M ⊙ M cut 5 × 10 10 M ⊙ [29]. We provide the corresponding M cut contour plots for the pseudo scalar and pseudo vector operators and their discussion in appendix C.2.

Matter power spectrum
The minimal halo masses derived in the previous subsection imply that the scalar operator leaves a similar suppression in the resultant matter power spectra to the case of the vector operator. In order to see this explicitly, let us consider our model where the DM scatters light fermions via the scalar operator. The scattering amplitude has a pure t-dependence given by (3.11). In other collision term expansion methods like that in [42], the scattering rate would be declared to be zero at the leading order. But as already shown in the previous subsection, we find that DM models with a scalar interaction can also account for the missing satellite problem.
To emphasize that scalar interactions are as important as vector interactions regarding the small-scale crisis problems, we adjust the free neutrino coupling parameters α ν in (3.14) and (3.15) to give the same cutoff mass and show that their linear matter power spectra are close to each other in figure 2. Here, we modify the public code CAMB [53] suitably to follow the coevolutions of cosmological perturbations of the DM (subsection 2.2) and the other components (e.g., baryons, photons, neutrinos, and gravitational potential). The small effects of the DM-neutrino interactions on the neutrino perturbations are neglected and the perfect fluid approximation (explained below) is used. Clearly, the shape of the power spectrum shows the characteristic features of the dark acoustic oscillations and the power on small scales is suppressed when compared to the CDM prediction.
Additionally, we check the validity of the perfect fluid assumption by comparing the results to the case of an imperfect fluid. To obtain a closed set of equations, we need to develop an approximation for f 03 and f 11 (see (2.50)-(2.53)). One way is setting them to be zero, defining the imperfect fluid approximation. This is valid when T χ /m χ ≪ 1, i.e., the free streaming of the DM particles is negligible after they decouple kinetically γ/H < 1 (see appendix B). Actually, we can also take σ = 0 and π = 0 for the adiabatic perturbations in  We take the same model with the scalar mediator as in figure 2. We take m χ = 1 TeV in both the perfect and imperfect fluid approximations. When the DM mass is lowered to m χ = 1 GeV and γ is kept fixed, the resultant matter power spectrum in the imperfect fluid approximation starts to differ from that in the perfect fluid approximation at wavenumbers larger than k 100 h/Mpc. the same limit, defining the perfect fluid approximation. Before the kinetic decoupling, all the variables f nℓ but f 00 and f 01 remain zero due to the damping term ∼ γ 0 f nℓ in (2.36). The former, corresponding to δ, does not have the damping term. The latter, corresponding to θ, has the source term ∼ γ 0 (θ TP − θ). One non-trivial check is to compare the resultant power spectra in the perfect and imperfect fluid approximations.
When the results from the perfect and imperfect fluid approximations deviate from each other, it does not necessarily mean that the imperfect fluid approximation gives a better description, but it just indicates that the perfect fluid approximation is not valid. To check if the imperfect fluid approximation gives a valid description or not, we need to compare the result from the treatment incorporating the full Boltzmann hierarchy, which is beyond the scope of this paper. Let us stress that the above deviation does not correspond to a limitation of the Fokker-Planck equation, which is valid as long as the momentum transfer in each collision is smaller than the typical DM momentum.
For a smaller DM mass with γ being fixed, we find differences in their power spectrum above a certain critical wavenumber as shown in figure 3. This is because the free streaming is sizable after the kinetic decoupling for the lighter DM. The results from the perfect fluid approximation are reliable below the critical wavenumber. On smaller scales, however, we may need to solve the full Boltzmann hierarchy (2.36). In appendix B, we give a more detailed discussion on the impact of higher order terms in the Boltzmann hierarchy and give a rough estimate of the critical wavenumber, where the results from the perfect and imperfect fluid approximations start to deviate.

Summary and outlook
In summary, we presented a consistent formulation that allows one to start from an underlying DM model and calculate its linear matter power spectrum. Regarding the small-scale crisis, -15 -

JCAP11(2016)043
the method is broadly applicable to essentially generic radiation interacting DM models that lead to a power spectrum suppression when compared to the standard cosmology on subgalactic scales.
In this paper, we focused on the case where the DM is in kinetic equilibrium with light and hidden fermions for a long time and the decoupling process was investigated for mediators of fundamentally different type. The new message is that not only a vector mediator at the MeV scale may solve all three small-scale problems at the same time [29], but we find that new classes of interactions may also solve at least the missing satellite problem. This result was unexpected from the point of view of previous literature [41][42][43], where the leading contribution to the momentum transfer rate is assumed to come from the scattering amplitude evaluated at Mandelstam t = 0. We explicitly derived an expansion method of the collision term where the scattering amplitude is t-averaged in the final form of the momentum transfer rate. This results in a different phenomenology from that in the previous literature for scattering amplitudes proportional to Mandelstam t, e.g., in the scalar, pseudo scalar, and pseudo vector interactions between the DM and the hidden fermions.
With this new insight, the classification of possible DM-radiation interactions, which are suppressing the abundance of dwarf galaxies, has to be revisited. During the preparation of this work, we have been informed that Bringmann et al. [54] have independently derived similar results concerning the possibility of kinetic decoupling at late times with in new classes of interactions. As a consequence, our work and studies by the latter authors may extend the list of realistic WIMP-like DM theories accounting for small-scale discrepancies.
As an important subtlety, we also discussed the validity of the perfect fluid approximation for the calculation of the power spectrum. We derive the consistent equations needed to be solved for an imperfect fluid treatment and compare the power spectra obtained from the perfect and imperfect fluid approximations. As indicated from figure 3, the perfect fluid approximation is limited by free streaming effects on the smallest scales. This may infer that we need to solve the full Boltzmann hierarchy to have reliable results for some models where the DM mass is small.
Our formulation, as a fundamental building block, in combination with N -body simulations would allow one to map DM models into the observational non-linear small-scale structure. We plan to combine baryonic feedback and DM-induced small-scale suppression to investigate the observational outcome. At present or in close future, this kind of sophisticated simulations are expected to shed more light on whether the small-scale crisis will be related to fundamental properties of DM or not. Even if the DM-radiation interaction does not resolve the small-scale crisis, our work and others can help to constraint DM models from a new perspective.

JCAP11(2016)043
x f /r on the model parameters and drop the last factor in (C.5)-(C.8) if the relic density constraint is used to reduce one of the parameters.
Furthermore, we remark that due to the presence of a light mediator and its long range property one has to include the Sommerfeld effect for DM annihilation in principle. This may lead to an O(1) correction of the DM coupling in order to produce the correct relic abundance, but including the effect is beyond the scope of this paper.

C.2 Minimal halo mass of the pseudo scalar and pseudo vector operators
In the case of the pseudo scalar and pseudo vector operators, the parameter space of interest spoils partially the effective propagator description, and thus γ/H does not have a simple power law dependence on the neutrino temperature like in the scalar and vector cases. Nevertheless, we derive analytically the scaling pattern of the cutoff mass from the effective propagator description, and compare it to the cutoff mass derived from the exact numerical evaluation of (3.6). The DM-neutrino scattering amplitudes for the pseudo scalar and pseudo vector operators are given by: Pseudo scalar operator: Pseudo vector operator: Pseudo scalar operator. The DM-neutrino scattering amplitude (C.9) via a pseudo scalar mediator has a pure t 2 -dependence. Within the effective propagator framework, γ/H depends therefore on a different power of T ν when compared to the scalar and vector operators: (C.12) Note that the mass of the mediator is close to the temperature ∼ 1 keV for subgalacitc cutoff masses. This spoils our effective propagator description as can be seen by comparing the exact numerical result with the effective description in figure 4.
Pseudo vector operator. The DM annihilation cross section via a pseudo vector mediator shows a z −4 enhancement in (C.3). At a first look, the limit z → 0 in the cross section seems to diverge and give rise to unitarity violation [56]. By embedding the model into a local U(1) gauge theory where both the mass of the DM and the gauge boson mass arise due to the spontaneous symmetry breaking via an additional scalar field, we show explicitly that this is not the case and the parameter region that we use to produce subgalactic cutoffs is in the perturbative regime. We denote the additional scalar by Φ and the local U(1) gauge invariant action reads  where / D + = / ∂ + ig χ / φγ 5 , D µ,−2 = ∂ µ − i2g χ φ µ , and V (Φ) = −µ 2 Φ ⋆ Φ + λ 2 (Φ ⋆ Φ) 2 and the fields transform such that (C.14) The vacuum expectation value of the field Φ in this potential is given by v = µ 2 /λ. We expand the scalar field around its minimum Φ(x) = v + (h(x) + iΦ 2 (x)) / √ 2 and get the following relevant quantities after symmetry breaking: m χ = λ Y v, m 2 φ /2 = 4g 2 χ v 2 , scalar mass m h = √ 2λv 2 = √ 2µ, Yukawa interaction −λ Y / √ 2χhχ = −m χ /(v √ 2)χhχ and scalargauge boson interaction +4 √ 2g 2 χ v hφ µ φ µ . The invariant amplitude of DM annihilation into two gauge bosons contains three terms: where y = m h /m χ . Now, the limit of m φ → 0 (z → 0), effectively meaning g χ → 0, results in a finite value of the annihilation cross section that is proportional to λ 4 Y . In the following, we  show that all parameters are in the perturbative regime and the scalar contribution (C. 15) can be ignored in the low energy expansion, so that (C.16) reduces to (C.3). When we choose z ∼ 10 −3 , due to the z −4 enhancement in the annihilation cross section (C.4), the DM coupling is forced to be tiny in order to satisfy the relic abundance constraint. A choice of m χ = 100 MeV leads to g χ = 1.0 × 10 −5 . With these choices, we derive λ Y = 0.03, v = where the relic density constraint on α χ (C.7) is inserted into (3.6). Note that the cutoff mass depends now on the DM mass unlike in the scalar and vector operator cases. In figure 5, the exact numerical solution of γ/H is compared to the cutoff derived from (C. 19), showing the valid range of the parameter space for the effective propagator description.