Distinguishing thermal histories of dark matter from structure formation

It is important to understand the implications of current observational constraints and potential signatures on the thermal history of dark matter. In this paper, we build the connection between the present-day velocities and the production mechanism of dark matter and find that the current observation on structure formation can be imposed to constrain the decoupling temperatures and the phase-space distribution of dark matter. We further explore the potential of distinguishing different possible thermal histories of dark matter with hypothetical future observational data. Using the freeze-in/-out scenarios as templates, we find that future precision data may uniquely identify the allowed parameter spaces for freeze-in and freeze-out, or even completely rule out one of the scenarios. This method can be more generally applied to other scenarios.


Introduction
Dark matter is a long-standing mystery in modern physics whose full picture involves delicate interplay between the physics on the smallest and the largest scales, and the dynamical evolution from the earliest times to the present day.With ample evidence from across many magnitudes in space and time, we now know that it does not interact appreciably with the standard-model (SM) fields, nor with itself, it is non-relativistic today and is non-baryonic, and it has a relic density that is about a quarter of the total energy of the universe.
Our knowledge about dark matter is however still limited, despite the fact that it comprises the majority of the matter component in the universe.We do not know its fundamental nature such as its mass and spin, whether it is elementary or composite, and whether it is single-component or multi-component.Likewise, although dark matter plays an important role in the dynamical evolution of the universe especially in shaping the CMB power spectrum and in seeding the structure formtion, the thermal history of dark matter remains unknown.We do not know how dark matter is produced in the early universe, nor do we know whether or not the dark-matter particles have ever been in thermal equilibrium.
Numerous models have been proposed to provide viable dark-matter candidates with some potential detection signatures, and a large suite of detection techniques, such as direct detection and indirect detection, have been developed to discover or constrain these candidates.While this top-down approach is important and might shed light on the properties of dark matter, it is of equal importance to consider the bottom-up approach, to figure out 1) what information about the early-universe dynamics of dark matter might be encoded in the observational data, and 2) to what extent one can extract this information from the current or future data to learn about the thermal history of dark matter and discriminate different models or scenarios.
Currently, most of the experimental or observational efforts are based on the nongravitational interactions of dark matter with the SM particles.For example, direct-detection experiments assume that dark matter can scatter off the target material (which is made of electrons and nuclei) of a detector with a appreciable rate.Indirect-detection experiments require dark matter to annihilate sufficiently rapidly in astrophysical environment where the dark-matter density is large.In addition, collider searches demand that dark matter can be produced via collisions of SM particles and manifest as missing transverse energy.On the other hand, all pieces of the known evidence for dark matter, such as the velocity dispersion of galaxies in the Coma Cluster [1,2], the galactic rotation curve [3], the cosmic microwave background (CMB) [4], the gravitational lensing [5] by galaxy clusters and the Bullet Cluster [6], originate from the gravitational interaction between dark matter and the SM.In addition, the measurements on CMB and the Bullet cluster also indicate that dark matter has negligible non-gravitational interaction strength.It is therefore important to explore other possible means of probing dark matter which do not rely on interactions with the SM particles other than gravity.
We are primarily interested in the connection between structure formation and the thermal history of dark matter through the primordial phase-space distribution of dark matter. 1n the one hand, simulations of structure formation and observations that depend on the cosmic structure such as measurements on the Lyman-α forest, the 21-cm signal, the halo mass function and the Milky-Way (MW) satellites can provide useful constraints on the dark-matter phase-space distribution.On the other hand, since the phase-space distribution is ultimately related to the production mechanism of dark matter, constraints from observations may have very different implications for different production scenarios.If an actual deviation from the CDM predictions is detected in future observations, the observational data can even be exploited to uniquely identify different possible thermal histories if predictions from different scenarios are sufficiently distinct.
Current studies of this sort have been focused on the warm-dark-matter (WDM) scenario in which dark matter freezes out while it is still relativistic, but becomes non-relativistic before the matter-radiation equality (MRE) to allow a consistent picture with the observed large scale structure.Since the phase-space distribution of WDM can be fully determined by the WDM mass, constraints in the literature are often quoted as bounds on the WDM mass.The most stringent constraints are derived from the observed Lyman-α forest with recent studies giving two lower bounds m WDM ≥ 3.5 keV or 5.3 keV based on different assumptions on the temperature evolution of the intergalactic medium [59].
For more general scenarios in which the phase-space distribution of dark matter is different from that of WDM, in principle numerical simulations are needed to obtain robust constraints.Nevertheless, several effective methods have been developed to recast the WDM bounds to the constraints in the corresponding scenarios [62,[68][69][70].In addition, numerical simulations based parametrized transfer functions have also shown important progress [62,70,71].
In this work, we shall focus on two most common dark-matter production mechanism, namely, thermal freeze-out and freeze-in.These two scenarios can be realized by the same particle-physics process but with different interaction strength.In the former scenario, the interaction between dark-matter particles and the thermal bath where these particles reside is rapid enough such that thermal equilibrium can be established between them, whereas in the latter scenario such interaction is too feeble to allow equilibration.Notice that, while WDM is one subset of the former scenario, our primary interest lies in a different subset, that is, the case in which freeze-out occurs non-relativistically.It turns out that in both the freeze-in and the non-relativistic freeze-out scenarios, dark-matter particles may acquire non-negligible velocities to impact structure formation if the dark-matter mass is sufficiently small, and thus can be potentially constrained by structure formation.Nevertheless, we shall see that the bounds on these two scenarios can be very different since the functional form of the dark-matter phase-space distributions and its dependence on physical parameters such as the dark-matter mass and the decoupling temperature differ significantly between the two scenarios.We shall then further explore this difference and discuss the potential in distinguishing these two distinct types of thermal histories with future observational data.This paper is organized as the following.In Sec. 2, we study the relation between the present-day velocities and physical quantities at the decoupling time via the phase-space distribution of dark matter resulting from different production mechanism.In Sec. 3, we take a vanilla dark-matter model and numerically calculate the phase-space distribution of dark matter in different regimes while focusing on the freeze-in and freeze-out cases.We also place constraints on these scenarios using observables from structure formation and discuss the different implications for the freeze-in and freeze-out cases.Then, in Sec. 4, we hypothesize several forms of future observational data on small scale structure, and evaluate the possibility of differentiating different scenarios from the future data.Finally, we conclude in Sec. 5.

Dark-matter velocities from different production mechanisms
The phase-space distribution f χ (⃗ x, ⃗ p, t) of the dark-matter particle χ is a function of spatial coordinates, momentum, and time.Since the universe is spatially homogeneous and isotropic at the zeroth order, the distribution function can be reduced to a function that only depends on time and the magnitude of momentum f χ (p, t).With f χ (p, t), quantities that are relevant for cosmology such as the number density n χ (t), the energy density ρ χ (t), and the pressure P χ (t) can be calculated via in which E = p 2 + m 2 χ is the energy of the dark-matter particle with mass m χ , and g χ is the internal degrees of freedom.
The dynamical evolution of the phase-space distribution is governed by the Boltzmann equation in which H ≡ ȧ/a is the Hubble parameter with a being the scale factor, and C[f ] is the collision term which contains integrals of the phase-space distribution of all relevant particle species.In general, the collision term describes how particles are added to, subtract from, or reshuffled within the phase space due to all possible interactions, whereas the Hubble term describes how particles redshift from higher momentum to lower momentum due to the expansion of the universe.Therefore, apparently, the phase-space distribution of dark matter depends on how dark matter is produced.In what follows, we discuss the phase space distribution resulting from two types of scenarios -the freeze-out scenario and the freeze-in scenario.

Production mechanism: freeze-out
In the thermal freeze-out scenarios, dark matter is assumed to be in thermal equilibrium with a thermal bath at earlier times during which the phase-space distribution is expected to be Fermi-Dirac or Bose-Einstein, i.e., where T χ is the temperature of the thermal bath that dark matter is in equilibrium with, and µ is the chemical potential.Note that, T χ is not necessarily identical to the temperature T of the SM thermal bath.At later times, dark matter decouples both chemically and kinematically.Let us assume that the decoupling occurs instantaneously at t dec when the temperature of dark matter is T χ (t dec ).Since the momentum p ∼ a −1 , the distribution at a later time t can be mapped to the distribution at the decoupling time t dec via where p dec ≡ pa(t)/a(t dec ).
To recast this distribution into a thermal form as Eq.(2.5) with physical quantities taken at time t, one would require the exponent in the brackets to match: Obviously, if µ is independent of p, the above equation can only be satisfied by a single value of the momentum.In other words, although the distribution f χ is simply related to its what it is at the decoupling time by redshift, it is in principle not evolving towards a thermal distribution with a different temperature and a single-valued chemical potential.However, for a distribution of the type in Eq. (2.6), most of the dark-matter particles in the phase space are concentrated at the average momentum ⟨p⟩ dec at the decoupling time if µ dec is not much larger than T χ (t dec ).In other words, after integrating the angular dependence, the integrated distribution function p 2 dec f χ (p dec , t dec ) only has non-negligible support at momenta that are not too far away from the average momentum.Thus, for the part of p 2 dec f χ (p dec , t dec ) that is non-vanishing, it is possible to define an effective temperature and chemical potential in the cases of relativistic and non-relativistic decoupling.

Relativistic decoupling vs. non-relativistic decoupling
If dark matter decouples relativistically, i.e., ⟨p⟩ dec ∼ T χ (t dec ) ≫ m, the energy of the darkmatter particle is dominated by the kinetic energy By defining it is easy to find that the factor namely, the distribution after relativistic decoupling still takes the approximate form of a relativistic thermal distribution characterized by effective temperature T eff and chemical potential µ eff which both scale like a −1 .Notice that, the shape of the distribution does not depend on the dark-matter mass even after dark-matter particles become non-relativistic at later times.
On the other hand, if dark matter decouples when it is still non-relativistic, i.e., m ≫ T χ (t dec ), we have .11)In this case, we can define and thus Eq. (2.7) can also be satisfied approximately by all momenta that are close to the average momentum Therefore, for non-relativistic decoupling, the resulting distribution is also effectively thermal with T eff χ ∼ a −2 .If m χ ≫ µ dec , the quantum effect from Bose enhancement or Pauli blocking is negligible, and then the distribution is effectively Maxwell-Boltzmann as f χ (p, t) ∼ exp(−p 2 /(2m χ T eff χ )).Knowing that the temperature and chemical potential can only be defined effectively, we shall from now on drop the superscript "eff".

Production mechanism: freeze-in
If the interaction between dark matter and the SM particles is too feeble to establish thermal equilibrium, dark-matter particles may still be produced through the freeze-in mechanism (see Refs. [72][73][74][75] and references therein).In principle, it is possible for these particles to establish thermal equilibrium within a dark sector sourced by the SM sector [75][76][77][78][79][80][81][82][83][84][85].However, in this type of scenarios, as dark matter equilibrates with the dark thermal bath, the distribution function of dark matter is simply thermal, which is no different from what is studied in the previous subsection.In other words, the equilibration process erases the "memory" before the dark freeze-out.Given these considerations, we shall only take into account scenarios in which dark-matter particles, once produced, never equilibrates.
The freeze-in production can proceed through various channels, and the analytical approximation of the final distribution function in many different cases has been studied in literature [69,[86][87][88][89]. Usually, such analytical expressions are written as a function of the comoving quantity p/(M a dec /a), where M is the mass of the heaviest particle in the production channel and thus marks the temperature scale at which the dark-matter production is complete.At lower temperatures, the production of dark matter is suppressed due to an insufficiency in either the kinetic energy or the number density of the source particles depending on whether M is one of the masses of the source particles or the products.Therefore, one can define a quantity T χ (t) ≡ M a dec /a(t) which plays the role of an effective "temperature" after dark matter is produced at t dec .The decoupling temperature in the dark sector is then simply T χ (t dec ) = M .For concreteness, we shall use the following analytical formula which is suitable for freeze-in through two-body decay or 2 → 2 scattering: where C is a constant that depends on the internal degrees of freedom, the decay width/cross section, and the masses of the particles involved in the production channel.For simplicity, we shall treat it as a normalization factor to match the relic abundance of dark matter.

Relation between mass and velocity
Since the phase-space distribution of dark matter after production takes different functional forms in different scenarios, the impact on structure formation and its dependence on other physical parameters, especially the dark-matter mass, is expected to be different among these Table 1: Values of A n and A p in each case, where FD, BE and MB stand for Fermi-Dirac, Bose-Einstein and Maxwell-Boltzmann distributions, respectively.scenarios.It is therefore of great interest to find the relation between the dark-matter mass and velocity in each case.We first notice that the distribution function after decoupling allows us to calculate the number density and average momentum of dark matter straightforwardly.Assuming that µ dec is negligible for simplicity, we have and where we use "RFO", "NRFO" and "FI" for relativistic freeze-out, non-relativistic freeze-out and freeze-in, respectively, the values of A n and A p in each case are listed in Table 1, and we have approximated both the Bose-Einstein and the Fermi-Dirac distributions by the Boltzmann distribution in cases of non-relativistic decoupling.Notice that we have kept µ in the second line of Eq. (2.15) since the chemical potential increases and asymptotes m χ for non-relativistically decoupled dark matter.
On the other hand, Eq. (2.16) alone enables us to relate the average velocity of dark matter today to relevant physical quantities at the decoupling time.Using the scaling relation between T χ and a in each case and the conservation of entropy in the SM sector, it is easy to find that where g ⋆,s is the relativistic degrees of freedom for entropy density in the SM thermal bath, and we have defined Note that, the use of entropy conservation implies our assumption that the SM sector is not transferring energy with any other sector after dark matter decouples.The relations in Eqs.(2.18) and (2.19) both suggest that, given one functional form of the phase-space distribution, a constraint on dark-matter velocities can be recast as a constraint Notice that, in the case of non-relativistic freeze-out, the curves cannot extend to arbitrarily small mass without violating the constraint on the dark-matter relic abundance.on the dark-matter mass.For typical values of the parameters, the average velocity of darkmatter particles can be larger by roughly one order of magnitude in the case of non-relativistic decoupling than in the case of freeze-in or relativistic freeze-out.This means, for a constraint on the average velocity ⟨v⟩ 0 , the corresponding constraints on the dark-matter mass m χ is typically stronger in the case of non-relativistic decoupling.
To allow a better comparison, we plot the relation between ⟨v⟩ 0 and m χ with the condition from the relic abundance imposed in FIG. 1.To be specific, we directly use the relations in Eq. (2.18) for freeze-in and relativistic freeze-out since the relation can be uniquely determined with at most one additional parameter C. On the other hand, we use both Eq.(2.18) and Eq.(2.19) for the non-relativistic freeze-out scenario in order to eliminate x dec .The current Lyman-α constraints on WDM mass, m WDM ≥ 3.5 keV and m WDM ≥ 5.3 keV [59], can then be converted to preliminary constraints on the present-day average velocity via the relation in the relativistic freeze-out scenario, which gives ⟨v⟩ 0 ≲ 2.1 × 10 −8 or ⟨v⟩ 0 ≲ 1.2 × 10 −8 , as indicated by the black horizontal lines. 3We can then apply these constraints on the average velocity to the cases of non-relativistic freeze-out and freeze-in.
We first notice that the constraints on m χ in the freeze-in scenario is very similar to that in the WDM scenario (which is exactly the relativistic freeze-out scenario) if C = 1.However, these lower bounds can drastically increase or decrease if C becomes smaller or bigger which is consistent with Eq. (2.18).In the non-relativistic freeze-out scenario, the constraints on m χ are overall much stronger than that in the WDM scenario, where the bounds can be m χ ≳ 5 keV or m χ ≳ 120 keV depending on different choices of η dec .
It is worth pointing out that the orange curves in the case of non-relativistic freeze-out do not extend arbitrarily to the left.This is because while the mass-velocity relations in Eq. (2.18) are constrained by the abundance Ω χ , the relations in Eq. (2.19) are derived directly from the relation between momentum and velocity and thus do not necessarily respect the constraint on the dark-matter relic abundance.Therefore, parameters in Eq. (2.19), {m χ , T dec , η dec , x dec }, in principle cannot vary independently from each other without violating the observational bound on the dark-matter relic abundance.Indeed, while holding η dec fixed and vary m χ in Eq. (2.19), one would need to adjust x dec (and thus T dec ) in order to satisfy the constraint in Eq. (2.18) at the same time.This is reasonable since x dec is essentially determined by comparing the interaction rate and the Hubble expansion rate when actually solving the Boltzmann equation, and thus depends on m χ through the cross-section and number density.The dependence on m χ in x dec can also be seen by noticing that the slopes of the orange curves differ distinctively from that of the red and blue curves, despite that all relations in Eq. (2.18) contains a factor of m It is also worth mentioning that, for the same decoupling temperatures, one in general does not expect ⟨v⟩ 0 from freeze-in to be larger than that from non-relativistic freeze-out.This can be seen from Eq. (2.19).To make the present-day average velocity from freezeout smaller than that from freeze-in, one need x dec ≲ 3 -a value close to the regime of relativistic freeze-out which makes our estimate unreliable.Therefore, the only possibility for any of the blue curves to surpass the η dec = 1 orange curve is to have η dec > 1 in the freeze-in case.While the information about η dec is encoded in the parameter C, and having η dec > 1 does not explicitly violate any assumption in the calculation that we have presented so far, the existence of a dark thermal bath hotter than the SM thermal bath would be severely constrained by current observations.We shall therefore neglect any part of the blue curve above the solid orange curve.
Other than constraining the dark-matter mass, the constraint on the average velocity also have non-trivial implication on physical quantities at the decoupling time which are extremely important for understanding the thermal history of dark matter.Equating Eq. (2.18) with Eq. (2.19), we find the relation between the dark-matter mass and the temperatures at decoupling . (2.21) The above equations suggest that, for a given dark-matter mass m χ , one cannot uniquely determine the temperatures T dec and T χ (t dec ) at the moment of decoupling.Nevertheless, if a bound from observations can be place on the dark-matter mass, regardless of how this bound is obtained, such information can then be used to constrain the parameter space spanned by T dec and T χ (t dec ) in a given production scenario which reveals information about the thermal history of dark matter.This can be seen in FIG. 2, where we show such constraints in each case in different panels.
For the relativistic freeze-out (or WDM) scenario in the upper left panel, the contours of m χ are also contours of the present-day velocity of dark matter ⟨v⟩ 0 .Therefore, the current constraint on m χ from structure formation, as indicated by the dashed and solid red contours, In addition, the blue contours in the top right panel shows the value of x dec .As described by the text, the gray areas indicate that, when using Eq.(2.21), either the solution of m χ is inconsistent with the assumption of the calculation (top left and bottom panels), or no solution can be found in this region (right panel).
follows the same trend as the dotted black contours.One immediate conclusion that can be drawn from such constraint is that WDM cannot decouple directly from the SM thermal bath.Instead, dark matter can only decouple relativistically from a dark thermal bath with a lower temperature.Depending on T dec , the temperature of the SM sector at the decoupling time, T χ (t dec ) has to be smaller by a factor of ∼ 0.3 or 0.1 at least.
In the case of non-relativistic freeze-out, we follow the procedure in FIG. 1 to turn the constraint on m χ in the case of relativistic freeze-out into the constraint on the present-day average velocity, i.e., ⟨v⟩ 0 ≲ 2.1 × 10 −8 or ⟨v⟩ 0 ≲ 1.2 × 10 −8 , which is again indicated by the dashed and solid red curves.Similar to examples shown in FIG. 1, we also observe that the constraints on dark-matter mass are now depending on the decoupling temperatures.The constraint for relativistic freeze-out, m WDM ≥ 3.5 keV (or m WDM ≥ 5.3 keV), now covers a broad range from roughly 5 to 70 keV (or from 8 to 120 keV) in the case of non-relativistic freeze-out.However, the exact constraints on mass cannot be unambiguously determined until the decoupling temperatures are specified.Moreover, since η dec can approach unity above the solid or dashed red curves, the current WDM bound from structure formation cannot completely rule out the possibility that dark matter might decouple from the SM thermal bath if such decoupling is non-relativistic.In addition, we also notice that the second line of Eq. (2.21) does not always have a solution on the entire η dec -T dec plane.This simply means that, for certain combination of decoupling temperatures, it is not possible to obtain the correct relic abundance without introducing additional assumptions.The corresponding region of parameter space is therefore shaded in gray.
For the freeze-in scenario, we first notice that there is an additional degree of freedom from the normalization factor C. Indeed, comparing the freeze-in case in Eq. (2.18) with that in Eq. ( 2. 19), we see that a change in C in the former equation amounts to a change in the product η dec /g ⋆,s (T dec ) in the latter equation.Physically, a variation in C (while holding the masses of the particles in the corresponding process fixed) would result in a change in the number density of dark matter at the decoupling temperature T χ (t dec ) = M , and thus one needs to adjust the mapping between T χ (t dec ) and T dec in order to obtain the correct relic abundance.Therefore, instead of plotting against η dec , we choose C 1/3 η dec as the horizontal axis such that we can unambiguously determine contours of m χ on the plane.For the average velocity, the red contours correspond to the choice of fixing η dec = 1, which can be seen as having dark matter produced directly from the SM thermal bath.In this situation, the contours of ⟨v⟩ 0 differs slightly from that of m χ due to the dependence on T dec in Eq. (2.19), and the current bound on average velocity covers a range from 10 to 30 keV (or from 17 to 50 keV).Equivalently, one can also fix C = 1, and then the contours of m χ can be directly turned into contours of ⟨v⟩ 0 using Eq.(2.18).In this situation, the present constraint is very similar to the case of relativistic freeze-out, which is m χ ≳ 3.3 keV (or m χ ≳ 5.0 keV).

Comments on ∆N eff
In the previous discussion, we have only re-interpreted and applied the WDM constraint from studies of the large scale structure.When exploring the implications of such constraint on the production of dark matter, other factors might also play a role and can be even more restrictive.
The major constraint relevant for the production of dark matter comes from the big bang nucleosynthesis (BBN) which occurs when the temperature of the SM sector drops below ∼ 2 MeV.First, dark matter (and other dark-sector species, if there is any) provides energy density in addition to the SM thermal bath which can affect the expansion rate of the universe during BBN.Second, if there is any energy transfer between the SM thermal bath and dark matter (or the dark sector) after neutrino decoupling, the temperature ratio between photons and neutrinos would be changed.Both effects are going to modify the prediction for the abundance of light elements.Moreover, if there is any dark species that stay relativistic after the matter radiation equality, the anisotropy power spectrum of the cosmic microwave background (CMB) will also be affected.Nevertheless, a systematic study on these effects is model-dependent and beyond the scope of this work.Here, we shall only briefly comment on the possible effects on N eff , the effective number of relativistic neutrino species.
In the case of relativistic freeze-out, we have seen from FIG. 2 that dark matter can only decouple from a dark thermal bath with a smaller temperature.However, this implies the existence of extra relativistic species which can potentially be constrained by N eff , the effective number of relativistic neutrino species.We can estimate this effect by noticing that the ratio between the energy densities of the dark sector and the SM neutrinos after neutrino decoupling is ∼ (11/4) 4/3 η 4 , up to corrections from the number of relativistic freedoms in the two sectors and the quantum statistics.Suppose there is no appreciable change in η from the beginning of BBN to the CMB time, the allowed range of η dec in FIG. 2 (η dec ≲ 0.3 or 0.1 depending on T dec ) indicates that the energy density of the dark-sector radiation is at most ∼ 3% that of the SM neutrinos, which can contribute to N eff at most comparable to the uncertainty of current measurements, and thus can easily evade the current bound by slightly decreasing η dec .Such a rough estimate can also be applied to the case of non-relativistic freeze-out or freeze-in.
For the case in which dark matter decouples directly from the SM thermal bath, which is not excluded by structure formation in the case of non-relativistic freeze-out and freeze-in, N eff will be modified only if dark matter decouples after neutrinos decouple.In the nonrelativistic freeze-out scenario, if dark matter is always in thermal equilibrium with the SM thermal bath during BBN but before decoupling, for m χ below a few MeV, dark matter can contribute appreciably to the total energy density of the universe when it is relativistic or close to being relativistic.As temperature drops, this energy density will be transferred to the SM sector.This will increase N eff if dark matter is primarily converted into neutrinos, or decrease N eff if it is preferentially converted into photons.In this situation, we would expect the η dec = 1 edge in the top right panel in FIG. 2 to be severely constrained by N eff below O(MeV).Nevertheless, if, during the BBN epoch, dark matter is decoupled at first, but then equilibrate with the SM while being non-relativistic before its eventual freeze-out, the impact on N eff can then be significantly relaxed.This type of scenario has been explored in [90], and its impact on BBN has been discussed in [91].We therefore do not consider the η dec = 1 region of sub-MeV dark matter in the non-relativistic freeze-out scenario as ruled out.
In the freeze-in scenario, the constraint from BBN is much easier to evade as the ratio ρ χ /ρ R ∼ n χ ⟨E⟩ /(g ⋆ (T )T 4 ) do not increase as we go back in time, if we neglect the possibility of equilibration within the dark sector which might exponentially suppress n χ as dark matter becomes non-relativistic.Therefore, during the entire BBN epoch, the fraction of the total energy density of the universe occupied by dark matter is maximum by the end of BBN, i.e., when T ∼ O(10) keV.Suppose dark matter decouples before the end of BBN (which is consistent with the allowed region in FIG. 2), a rough estimate shows that the ratio nm χ /ρ R ∼ a/a eq ∼ T eq /T , is at most O(10 −4 ) during the BBN epoch, where T eq ∼ O(1) eV is the temperature at the matter-radiation equality.Thus, unless the kinetic energy of dark matter is orders of magnitude larger than its mass energy by the end of BBN, which might again run into trouble with constraints from structure formation, we do not expect any observable effect on BBN in the freeze-in scenario.

Constraints on phase-space distribution
The analysis in the previous section provides a useful estimate in constraining the mass and decoupling-temperatures of dark matter in several different scenarios.However, such constraints are preliminary as they are obtained by simply converting the current Lymanα constraints on WDM mass into constraints on the average velocity of dark matter, and the latter provides very limited information other than a characteristic length scale below which the formation of structure is suppressed.On the other hand, the actual astronomical observables such as the Lyman-α forest often concern the structure of the universe over a wide range of length scales.Therefore, to more rigorously constrain the scenarios considered in this work, it is important to go beyond the estimate from a single average velocity and study structure formation in more detail.
For this purpose, in what follows, we shall first demonstrate how to obtain the phase space distribution of dark matter by solving the Boltzmann equation at the level of the phasespace distribution f χ (p, t).The phase-space distribution is then used as an input to compute the density perturbations and the matter power spectrum via the public code CLASS [92,93].Eventually, we shall use the matter power spectrum to derive more rigorous constraints in our scenarios.

Phase-space distributions from freeze-in and freeze-out
As we have shown in Sec. 2, the general form of the Boltzmann equation for the phase-space distribution f χ (p, t) reads where C[f ] is the collision term, and for a general process a where ±" describes the Bose-enhancement/Pauli-blocking effects from the quantum statistics, and the two amplitudes correspond to the forward and inverse processes respectively.In principle, the collision term in Eq. (3.1) contains all the relevant processes such as the annihilations, scatterings or decays that dark-matter particles participate.
To proceed further, we consider a concrete example in which the dark-matter particle χ is predominantly produced through a 2 → 2 annihilation ψ + ψ ↔ χ + χ, where we assume that χ and χ are distinct particles and set g χ = 2.For the particle ψ, we assume that it is either a massless species in the SM sector or a massless species that thermalizes in the dark sector.In the latter case, we assume that the total energy density of the dark sector is negligible compared with that of the SM thermal bath, so that the Hubble parameter is dominated by the contribution from SM-sector energy density.And this can be easily satisfied as long as η is not too close to unity.Note that, by choosing ψ to be massless, dark matter will decouple non-relativistically in either the freeze-in scenario or the freeze-out scenario since the decoupling condition is T χ ≲ m χ for the former, and n χ being Boltzmann suppressed, i.e., e −mχ/Tχ ≪ 1, for the latter.Nevertheless, we shall simply call the latter scenario freeze-out for brevity when there is no confusion.For simplicity, we shall also neglect elastic scattering processes like ψ + χ ←→ ψ + χ in this work.The effects of elastic scatterings on the phase-space distribution have been discussed in Ref. [94].Therefore, in our example, the collision term for the 2 → 2 annihilation process can be simplified as where s and t are the Mandelstam variables, p * χ (s) is the momentum of χ in the center-ofmass frame, |M| 2 is the squared amplitude for the annihilation χ χ → ψψ, the superscript "eq" denotes the value of f χ in thermal equilibrium, and we have approximated f eq χ by the Maxwell-Boltzmann distribution, i.e., f eq χ (p χ ) = e −Eχ/Tχ , and ignored the quantum-statistics factors.The explicit form of E max/min χ (s) and t max/min (s) can be found in Ref. [94] and are not repeated here for brevity.
Note that, the above equation is general in the sense that the freeze-out and freeze-in scenarios are unified under one single framework -it allows one to switch from the freezeout scenario to the freeze-in scenarios by simply tuning the size of the squared amplitude |M χ χ→ψψ | 2 .In particular, the collision term always tends to bring f χ to its equilibrium value f eq χ as can be seen from the part in the squared brackets in Eq. (3.3).In the freeze-out scenario, the squared amplitude is large enough such that the annihilation rate of χ can be sufficiently rapid to set f χ = f eq χ until χ decouples.On the contrary, in the freeze-in scenario, the squared amplitude is too small such that f χ ≪ f eq χ and the annihilation of χ has negligible effect on the final dark-matter density.As is pointed out in Ref. [94], the transition between the two regimes can be conveniently quantified by introducing a quantity to specify the interaction rate of the annihilation process, where ⟨σ χχ→ψψ v⟩ is the thermally averaged cross (see Appendix A for definition).Thus, the freeze-out scenario correspond to γ χ ≫ 1, whereas the freeze-in scenario corresponds to γ χ ≪ 1.In addition, there exists a transition regime when γ χ ∼ 1, where the inverse process χ χ → ψψ has significant impact on the production of χ, even though the chemical equilibrium between χ and ψ is never established.In Fig. 3, we present contours of fixed relic abundance, Ω χ = 0.25, in the m χ -|M| 2 plane for cases with η dec = 0.1, 0.2, 0.5 and 1.In all the cases, we assume that |M| 2 is a constant for simplicity, and its value is adjusted to obtain Ω χ = 0.25 after solving the Boltzmann equation.As indicated by the color, the interaction rate characterized by γ χ increases from γ χ ≪ 1 to γ χ ∼ 1, and eventually to γ χ ≫ 1 along each contour as one goes from the lower right part to the upper right part in a clockwise direction.For each contour, there is a lower limit m min χ on the dark-matter mass below which the correct relic abundance cannot be obtained, and this smallest mass increases when η dec decreases.Above this lower limit, there are always two values of the squared amplitude that are able to generate the correct relic abundance at a fixed mass, which give rise to two possible thermal histories.
Motivated by the results in Ref. [94], we can roughly identify the freeze-in, the transition and the freeze-out regimes with γ χ < 0.1, 0.1 ≤ γ χ ≤ 10, and γ χ > 10.It is easy to notice that, for each fixed η dec , the transition regime only takes a small portion in the entire mass As indicated by the color bar, the interaction rate γ χ varies from γ χ ≪ 1 to γ χ ≫ 1 along each contour.Thus, each contour corresponds to a path that goes through the freeze-in, transition, and freeze-out regimes.The diamond markers are placed at γ χ = 0.1 and 10 to roughly identify these three regimes.
range, whereas most of the available parameter space corresponds to the freeze-in and the freeze-out regimes as there is no upper limit for the m χ other than the potential unitarity bound [95].Based on this observation, we shall therefore focus on the freeze-in and the freezeout regimes in what follows.We shall assume η dec = 1 in these two regimes, and show in Sec.3.3 how to extend our results in the parameter space with η dec < 1.
In the left panel of Fig. 4, we show the evolution of dark-matter comoving mass density N χ m χ , where N ≡ n χ a 3 is the comoving number density.We have normalized the present-day value of the scale factor a(t 0 ) = 1 such that N χ (t 0 )m χ ≃ ρ χ (t 0 ).The examples shown here are selected from the data points on the η dec = 1 contour in Fig. 3.The solid and dashed bunches of curves correspond to the freeze-out and the freeze-in scenarios, respectively, whereas the dotted bunch of curves shows the value of N eq χ m χ -the comoving mass density if χ is in thermal equilibrium.Notice that, the dotted curves, and thus the solid curves which overlap with the dotted curve before freezing out, are lowered as m χ decreases.On the other hand, the dashed curves appear to have the opposite behavior before freezing in.The freeze-out and the freeze-in bunches become closer when m χ is dialed smaller.Nevertheless, the distinct separation of them suggests that for each m χ within this mass range there are always two different thermal histories that give rise to the same relic-abundance of dark matter, even though the dark-matter particles are created via the same 2 → 2 process.
In the right panel of Fig. 4, we plot the resulting phase-space distribution g χ (p) ≡ p 3 f χ (p, t 0 ) (normalized by the comoving number density N χ ) which is the appropriate distribution with respect to the momentum on a logarithmic scale [64][65][66].We find that, despite having the same relic abundance, the two production mechanisms result in two distinct sets of momentum distributions, among which the momentum distributions from the freeze-out production always possess an overall higher momentum and a smaller width.Another no- ticeable feature is that the momentum distributions from the same production mechanism tend to be extremely similar which is consistent with our discussion in Sec. 2. The variation in the dark-matter mass m χ only causes a very slight horizontal shift.In the freezein scenario, this slight shift to the left as m χ increases is only caused by the change of g ⋆ (T ) after dark matter is produced since the effective temperature of dark matter today T χ (t 0 ) = M a dec /a 0 = T 0 [g ⋆,s (T 0 )/g ⋆,s (M )] 1/3 .In the freeze-out scenario, the same effect from g ⋆ (T ) also exists and is even stronger since the effective temperature T χ ∼ a −2 .However, as is shown in the left panel, the requirement on producing the correct relic abundance forces dark matter to decouple "later" (at a smaller m χ /T ) when m χ increases.Since this additional effect allows dark matter to thermalize with the radiation bath for a longer time, it tends to work in the opposite direction by shifting the distribution to higher momentum.It turns out that the second effect is always stronger than the first one within the mass range that we are showing since g ⋆ (T ) does not vary appreciably below O(100) keV.Thus, the distribution for freeze-out always shifts to the right when m χ increases.Analytically, such behavior can be understood from the second line of Eq. (2.19) where the competition between the dependence on m χ and the dependence on x dec is explicit.
Although both the thermal history and the final phase-space distribution in the freeze-in and freeze-out regimes are extremely different, in principle, one can always dial down the mass to enter the transition regime and find a continuous and smooth interpolation between the two types of distribution functions.This is similar to the results shown in Ref. [94], where, for fixed m χ , the production of χ goes through the three regimes continuously by simply increasing γ χ , and the resulting phase-space distributions after dark-matter production also interpolate smoothly between the two limits set by the freeze-in and freeze-out mechanism.However, for the η dec = 1 case here, the transition regime correspond to m χ ≲ O(1) eV, which is well excluded by current constraints.We thus do not plot the results in the transition regime for η dec = 1.Instead, for completeness, we present in Fig. 5 the results in the transition regime for η dec = 0.1 which is more likely to be allowed.Just like the previous figure, for each m χ > m min χ , there are always two solutions with different γ χ that give rise to the same relic abundance -one with a larger γ χ and thus is closer to establishing thermal equilibrium, and the other with a smaller γ χ and is thus further away from equilibrium, even though the annihilation of dark-matter particles are not negligible in both solutions.Such behavior is also reflected in the final phase-space distributions shown in the right panel, where the solutions from having a larger γ χ always results in a momentum distribution with overall higher momenta but smaller width.As m χ decreases towards m min χ , the two solutions approach each other, which is consistent with our observation in Fig. 3.Eventually, when m χ → m min χ ≃ 12.75 keV, the two solutions would merge.This is evident in both panels as the curves from the two solutions for m χ = 12.75 keV almost overlap with each other.
The reason that there is no solution for m χ < m min χ can also be glimpsed from the left panel of Fig. 5.As we have mentioned before, the particle interactions always tend to bring the χ particles into thermal equilibrium.Thus, we see in all cases that the m χ N χ curve goes up and evolves towards m χ N eq χ before reaching its maximum.However, since the maximum is bounded by the value of m χ N eq χ in the relativistic regime, this also means that it is impossible to obtain the correct relic abundance if m χ N eq χ < ρ χ (t 0 ).In fact, as can be inferred by the behavior when m χ → m min χ , further lowering m χ would bring down m χ N eq χ and render it impossible to generate enough dark-matter particles even if the plateau of m χ N eq χ is still higher than ρ χ (t 0 ).

Density perturbations
Beyond the relic abundance of dark matter, the thermal history of dark matter also manifests itself in the process of structure formation via the phase-space distribution.To better understand how this works, let us briefly review how the phase-space distribution affects structure formation.Following Ref. [96], the perturbed distribution function can be written as where τ ≡ dt/a(t) is the conformal time, q ≡ ap is the comoving momentum, n ≡ ⃗ p/p is the direction of the momentum, and f χ,0 (q) is the unperturbed zeroth-order distribution, i.e.,the distribution obtained in last section.In the Synchronous gauge, the line element is given by where τ is the conformal time, and the perturbation of the metric h ij can be decomposed as , where h ≡ h ii is the trace part of h ij (not to be confused with the little h in the Hubble constant), and h ∥ ij , h ⊥ ij and h T ij are the traceless parts, satisfying By definition, h ∥ ij and h ⊥ ij can be written in terms of a scalar filed µ and a divergenceless vector ⃗ A, respectively: Therefore, the scalar mode of the metric perturbations are characterised by the two scalar fields h and µ, while the vector and the tensor modes are characterised by ⃗ A and h T ij , respectively.
In this work we will only consider the scalar modes of h ij , which can be Fourier transformed as where h( ⃗ k, τ ) and η( ⃗ k, τ ) are the trace part and the traceless part of the scalar mode in the Fourier space which are related to the scalar fields h and µ in the real space.Since we have assumed that χ is collisionless after production, the evolution of Ψ obeys the first-order collsionless Boltzmann equation in the Fourier-space where ϵ ≡ aE = q 2 + a 2 m 2 is the comoving energy.Since Eq. (3.10) is independent of the azimuth angle, we can expand Ψ in a Legendre series Inserting Eq. (3.11) into Eq.(3.10) and using the orthonormality of the Legendre polynomials, we can obtain the Boltzmann hierarchies for Ψ l : (3.12) Notice that, even though we assume that dark-matter particles are collisionless, they are still coupled with other species such as photons and baryons through metric perturbations encoded in h and η.
With the moments Ψ l , the perturbation in physical quantities, such as the energy density, pressure, energy flux and shear stress, can be computed straightforwardly by integrating different moments.For our purpose, the fractional density perturbation can be expressed as where ρχ is the zeroth-order homogeneous background energy density defined in Eq. (2.2).
The effects from the phase-space distribution at different length scales can then be reflected by the evolution of the Fourier modes δ k which we obtain via the CLASS code.In particular, we have truncated the Boltzmann hierarchy at ℓ = 50, and have found that our general results are not affected significantly as long as ℓ > 8.We present in Fig. 6 the evolution of δ k for the CDM scenario and the freeze-in/freeze-out scenarios with three dark-matter masses m χ = 10, 20 and 30 keV at three wavenumbers k = 7.1, 14.2 and 42.8 h/Mpc.Overall, the suppression effect from having non-negligible velocities become weaker at larger length scales which correspond to smaller wavenumbers k.For small enough k whose corresponding length scale is larger than the maximum free-streaming distance of the darkmatter particles, no suppression in δ k is expected.This is evident in the left panel, where k = 4.3 h/Mpc, as all the curves overlap with the black curve which shows the evolution in the CDM scenario.As k increases to larger values, the deviation from CDM occurs.Thus, in the middle panel with k = 14.2 h/Mpc, we see the red, green and blue curves all peel off from the black curve.For the same production scenario (freeze-in or freeze-out), a smaller mass corresponds to a stronger suppression.On the other hand, for the same m χ , the suppression in the freeze-out scenario is stronger than that in the freeze-in scenario.Both behavior can be understood from the right panel of Fig. 4 -1) although the momentum distribution is nearly the same for the same production mechanism, the velocity is larger for smaller darkmatter mass; 2) for the same m χ , phase space distribution from freeze-out scenario possesses overall larger momenta than that from the freeze-in scenario.As k further increases, the deviation from the CDM scenario become more distinct in all the cases while the differences among all the cases also become very significant.Thus, in the right panel, we see that at k = 42.8 h/Mpc all the curves are well separated, and the suppression in the freeze-out scenario is always stronger than that in the freeze-in scenario.The examples in Fig. 6 demonstrate three important points.First, for a particular type of production mechanism, the suppression on the perturbation modes δ k increases with the wavenumber k.Second, such suppression is stronger for smaller dark-matter mass m χ .Third, dependence of the suppression on the wavenumber k and the dark-matter mass m χ is significantly different in the freeze-in and the freeze-out scenarios.The differences in these perturbation modes δ k therefore offer us an avenue to examine the predictions of different thermal histories, which, in turn, can be used to constrain or even distinguish these thermal histories.

Matter power spectrum and transfer function
The observations from the previous section shows that it is reasonable to examine the whole spectrum of δ k when studying the impact of dark-matter velocities on structure formation.Such physical quantity is called the matter power spectrum: Since we are always interested in studying the deviation from the CDM scenario, it is convenient to focus on the squared transfer function defined as where P CDM (k) is the matter power spectrum obtained for cold dark matter.In Fig. 7, we show a few examples of the squared transfer function associated with the cases studied in Fig. 4. The results from the freeze-out and the freeze-in scenarios are presented in left and right panels, respectively.As expected, we observe that in each panel the suppression on T 2 (k) is stronger as the dark-matter mass decreases.Comparing across the panels, we also find that for fixed m χ , the small-scale suppression on T 2 (k) is always greater in the freeze-out scenario than in the freeze-in scenario.
Similar to the WDM scenario, where the transfer function can be fitted by [50,52]  with ν = 1.12 and 11 it proves convenient to analytically fit the transfer function resulting from the freeze-in and the non-relativistic freeze-out scenarios.We find that, although the WDM fitting formula can fit well the transfer functions in the non-relativistic freeze-out scenario with a different α, the fit for the freeze-in scenario is not as good.A better fit for the freeze-in scenario can be done using the {α, β, γ} formula proposed in Refs.[62,71]: We shall therefore adopt Eq. (3.18) for both the freeze-in and freeze-out cases.In analogous to Eq. (3.17), we choose α to be a function the dark-matter mass.In order to allow the parameter α to have a different mass dependence in different scenarios, we further decompose α as where we have inserted the factor (g * ,s (T 0 )/g * ,s (m χ )) 1/3 to include the shift effect due to the change of g ⋆,s after dark matter is produced.With α specified, the parameters β and γ can be fixed based on the goodness of fit.We find that the transfer functions can be fitted to high accuracy by choosing {α 1 , α 2 , β, γ} = {0.32,−0.80, 2.24, −4.46} for freeze-out, and {α 1 , α 2 , β, γ} = {0.30,−0.87, 2.28, −1.51} for freeze-in.Note that the values of β and γ for freeze-out are chosen to match the corresponding exponents 2ν and −10/ν for the WDM fitting functions.Therefore, the overall shape of the fitted transfer function is the same in the freeze-out and the WDM scenarios.In addition, Eq. (3.16) or Eq.(3.18) suggests a symmetry transformation: which holds approximately true in the actual transfer functions and is exactly true in the fitting formula.Therefore, when plotted against k on logarithmic scale, a variation in α caused by a change in m χ amounts to sliding T 2 (k) rigidly in the horizontal direction.Indeed, the reason for this behavior lies in the similarity between the momentum distributions within either the freeze-out or the freeze-in scenario, which we have already seen in Fig. 4. We shall prove this in the Appendix B.
The above fitting formulas allow us to conveniently compare the transfer functions in our cases with the standard WDM scenario over a large range of dark-matter mass.To estimate how the Lyman-α constraints in the WDM scenario are projected to scenarios studied in this work, namely, freeze-in and non-relativistic freeze-out, we adopt the widely used "half-mode analysis" [61] and the "δA analysis" [62].
In the half-mode analysis, the half-mode k 1/2 is defined via The Lyman-α constraint for a particular model is then set by demanding that where T 2 WDM (k) is obtained from Eq. (3.16) by assuming a certain reference WDM mass.On the other hand, the δA criterion is based on comparing the quantity A defined as where k min/max = 0.5 and 20 h/Mpc respectively, corresponding to the MIKE/HIRES and XQ-100 combined dataset used in [59], and P 1D (k) is the "one-dimensional" matter power spectrum The quantity δA is then defined as the fractional deviation from the CDM scenario The Lyman-α constraint for a particular scenario is then estimated by requiring that δA ≤ δA ref , where δA ref is a reference value in a WDM scenario.In particular, we find that δA ref ≈ 0.50 for m WDM = 3.5 keV, and δA ref ≈ 0.36 for m WDM = 5.3 keV.With these two methods, the current constraints on WDM, m WDM ≥ 3.5 (or 5.3) keV, can be mapped easily to the scenarios of freeze-in and non-relativistic freeze-out.We find that, in the freeze-in scenario, the constraints are given by Mappings between the constraint on m WDM and that on m χ based on the half-mode analysis (solid curves) or the δA analysis (dashed curves) in the freeze-in (blue) and the nonrelativistic freeze-out (red) scenarios.Notice that the two curves overlap completely in the latter scenario.As indicated by the black lines which show the current bounds, a Lyman-α constraint on the WDM mass can be directly mapped to a constraint on m χ in a particular scenario.
which is expected as the transfer functions possess the same shape as the WDM ones.In each case, the larger (smaller) value in each case corresponds to the projection of the more stringent (conservative) bound on WDM.Overall, the results obtained from the two methods are extremely similar with the relative difference in the freeze-in case no more than 5%.
The results above shows that the dark-matter mass in both the freeze-in and the nonrelativistic freeze-out scenario are constrained at much larger values than those in the WDM scenario.This is consistent with our discussion in Sec.2.3, which is based on the present-day average velocity ⟨v⟩ 0 .The reason roots in the fact that the WDM can only be generated from a dark sector that is much colder than the SM thermal bath, whereas the cases we are comparing here all have η dec = 1.It is therefore also foreseeable that such constraints can be weaker if η dec < 1 which is indeed the case as we can see by the end of this section.In addition, we have also find that the constraints obtained here are considerably weaker than the preliminary results in Sec.2.3, where the constraints are estimated to be m χ ≳ 30 (or 50) keV for freeze-in, and m χ ≳ 70 (or 120) keV.This therefore suggests that while a single free-streaming velocity can provide a good order of magnitude estimate on the constraints from structure formation, the constraints it provides is in practice not very accurate.
Beyond the current bounds at 3.5 or 5.3 keV, the two methods can be more generally used to recast the Lyman-α constraint for any reference WDM mass, as is shown in Ref. [66].Indeed, one can create a mapping between a WDM mass and m χ in a specific scenario when their Lyman-α constraints are the same using one particular method.In Fig. 8, we show such mappings between m χ and m WDM in the freeze-in and the non-relativistic freeze-out scenarios obtained by allowing m WDM to vary in the two recasting procedures.Thus, if the Lyman-α constraints on WDM is updated in the future, one can easily find the equivalent constraints for the freeze-in ore non-relativistic freeze-out scenarios.

Halo mass function
The distinction between the freeze-in scenario and the non-relativistic freeze-out scenario exists not only in the linear regime, but also in the non-linear regime.One of the important physical quantity for assessing the non-linear evolution of the density perturbation is the halo mass function dn/d ln M , defined as the comoving number density of dark-matter halos at a logarithmic interval of the halo mass M .The halo mass function is in principle obtained from N-body or hydrodynamic simulations, which are computationally expensive.Nevertheless, approaches have been developed to analytically calculate or parametrize the halo mass function.For example, the Press-Schechter formalism [97] enables one to assess the halo mass function using the linear matter power spectrum P (k).Recent study based on stochastic differential equations [98], and that based on the mass and energy cascade theory in dark matter flow [99], as well as the work towards a general parametrization [100,101] have shown good agreement with results from simulations.In this paper, we follow the extended Press-Schechter framework [97,[102][103][104][105][106] to derive the halo mass function.In particular, we use a sharp-k window function which has the advantage that the result is sensitive to the shape of the linear matter power spectrum rather than to the shape of the window function [56].The halo mass function at the present-day thus takes the form where ρ = Ω m ρ crit is the average energy density of matter, δ c ≃ 1.686 is the critical overdensity, the mapping between R and M is given by with c W ≈ 2.5 [58], and ν(M ) ≡ δ 2 c /σ 2 (t 0 , M ), with σ 2 (t 0 , M ) being the spatial average of the variance of δ(⃗ x, t): where A ≃ 0.3222 and α = 0.3.A detailed derivation can be found in Refs.[56,58,65,107].We present in Fig. 9 the halo mass functions for a range of dark-matter masses from 10 keV to 100 keV, where the results for the non-relativistic freeze-out and the freeze-in scenarios are shown in left and right panels, respectively.The halo mass function in the CDM scenario is shown in both panels by the black curve.Apparently, the halo mass function in the CDM scenario (dn/d ln M ) CDM increases monotonically as the halo mass M decreases.On the other hand, the number of dark-matter halos with smaller masses is clearly suppressed within the range of masses considered in both the freeze-in and the non-relativistic freezeout scenarios.For instance, if m χ = 10 keV, the number of halos starts to be suppressed below M ≲ 10 11 M ⊙ /h.This suppression effect essentially originates from the suppression in the matter power spectrum P (k) which enters the evaluation of dn/d ln M by itself and through the variance σ 2 (t, R).Since the sharp-k window function cuts off all contribution from modes whose wavenumber k > 1/R(M ), it is not surprising that the halo mass function is only suppressed at a halo mass M when the matter power spectrum is suppressed at the wavenumber k < 1/R(M ).Therefore, we find that the halo mass function generally follows the trend in the matter power spectrum -the halo mass functions in our scenarios coincide exactly with the (dn/d ln M ) CDM at the larger-mass end, but gradually peel off from the CDM curve at smaller masses.Similar to what we have seen in the squared transfer function, a smaller m χ leads to stronger suppression at larger mass scales in each scenario, while for the same m χ the suppression in dn/d ln M is always stronger in the non-relativistic freeze-out scenario.
It is also interesting to notice that for the same dark-matter mass, as the halo mass M decreases, while the initial peel-off from (dn/d ln M ) CDM appear at locations that are close in M in both scenarios, the decrease at further smaller M is always more drastic in the non-relativistic freeze-out scenario than in the freeze-in scenario before entering the regime of rapid acoustic oscillations where dark-matter particles are all free-streaming.The reason behind this phenomenon can be glimpsed from the right panel of Fig. 4, which shows that while the distribution function resulting from the freeze-in scenario is overall colder, it is also wider.The tail of the freeze-in distribution function at larger momentum stretches close to the higher-momentum end of the distribution in the non-relativistic freeze-out scenario which gives rise to a peel-off location close to the one in the freeze-out case.On the other hand, the lower-momentum tail of the freeze-in distribution function extends far into the low-momentum region which prevents a complete suppression in the halo formation on the corresponding scales.
Indeed, the shape of the halo mass function encodes the information of the phase-space distribution.It is not difficult to conjecture that a correlation between the halo mass function and the phase-space distribution exists via the following mapping where v(k) stands for a dark-matter velocity whose corresponding free-streaming length d FSH ∼ 1/k.Recent developments on this correlation can be found in Refs.[64,65] where FIG. 10.The expected number of subhalos N SH with masses M > 10 8 h −1 M ⊙ within a MW-sized galaxy as a function of the dark-matter mass.The results from the WDM, freezein, and non-relativistic freeze-out scenarios are colored in green, blue and red, respectively, whereas the two horizontal lines indicate the result from a pure CDM N-body simulation [55], and the estimate on the number of observed MW satellites.
the information from structure formation can be exploited to reproduce the primordial phasespace distribution of dark matter.

Milky-Way satellite counts
Another important physical observable relevant for the non-linear evolution of density perturbations is the number of satellite galaxies N SH residing within a host galaxy of mass M 0 .Considering the observability, only the satellites that are more massive than a threshold M min are taken into account.Thus, the expected number of detectable subhalos is given by where dN SH /dM is the subhalo mass function.Following Refs.[58,65,108], for the shape-k window function that we have been using, the subhalo mass function takes the form where N SH is a normalization factor.It is of particular interest to apply the above formula in the case of a MW-like galaxy of mass M 0 = 1.77 × 10 12 M ⊙ /h, for which the results from the N-body simulation in the Aquarius project [55] predicted that N SH = 157 for the number of subhalos above M min = 10 8 M ⊙ /h in a pure CDM scenario. 5Based on this number, we obtain numerically that N SH ≈ 46.9.
On the observational side, we follow Refs.[58,65,68,109,110] to estimate the observed number of MW satellites N (obs) SH above M min by multiplying the 15 ultra-faint satellites discovered by the Sloan Digital Sky Survey (SDSS) by a factor of 3.5 to account for the limited sky coverage and then add to it 11 classical satellites which were discovered pre-SDSS.This method yields N (obs) SH ≈ 63, and thus for a viable dark-matter candidate, we expect that number of theoretically predicted MW subhalos N SH > N (obs) SH .In Fig. 10, we show the expected number of MW satellites as a function of the darkmatter mass m χ in the freeze-in and non-relativistic freeze-out scenarios, where we have used the {α, β, γ} fitting formula in Eq. (3.18) together with the decomposition in Eq. (3.19) to obtain the matter power spectra for arbitrary m χ .As we see, while the results in the freezein and non-relativistic freeze-out scenario both start from zero as m χ → 0 and asymptotes to N CDM SH = 157 at larger enough m χ , their difference in the intermediate range, m χ ∼ O(10) -O(100) keV, is clearly distinct.Similar to the discussion on the halo mass function, this difference in N SH originates from the fact that for the same mass m χ the distribution resulting from the non-relativistic freeze-out scenario is always overall hotter than that from the freeze-in scenario.Therefore, the observed number of the MW satellites imposes a stronger constraint on the non-relativistic freeze-out scenario.With N Comparing to the Lyman-α constraints in Eqs.(3.26) and (3.27), it is noticeable that the constraint from the MW satellite count is slightly weaker for the WDM scenario, but significantly weaker in either the freeze-in scenario or the non-relativistic freeze-out scenario.

Improved constraints on decoupling temperatures
Thus far in this section, we have assumed η dec = 1 when deriving the constrains on the dark matter mass.To extend our results to the regime with η dec < 1 one in principle needs to solve Eq. (3.1) with a T χ that is distinct from the temperature of the standard-model thermal bath.However, we notice that the shape of the distribution in the non-relativistic freezeout or the freeze-in scenario is very well described by the Maxwell-Boltzmann distribution or the distribution function in Eq. (2.14).Therefore, instead of doing a computationally expensive scan with η dec < 1, it is possible to identify the result from one combination {m χ , η dec } to that from anther combination {m dec } provided that these two combinations produce approximately the same present-day velocity distribution.Indeed, this identification amounts to requiring m in the non-relativistic freeze-out scenario, or dark matter can be obtained form any other source, there is still a huge overlap in the allowed mass range in which it is essentially impossible to uniquely identify the thermal history of dark matter.Therefore, in this section, we explore the extent to which one can unambiguously distinguish these two scenarios from future observations on the cosmic structure.We emphasize that, while we use the comparison between the freeze-in and the non-relativistic freeze-out scenarios as an example, an analysis of this kind can be more generally applied to comparisons between many other models or scenarios.
Since the form of future data is somewhat arbitrary, for generality, let us assume that future measurements on a physical quantity X(y) at the location y in the parameter space can be expressed as where X ± (y) is the upper/lower limits of X(y), X ref (y) is the best-fit reference value obtained from the data, and σ ± X (y) is the relative uncertainty for the upper/lower bound.In what follows, we shall explore two types of assumptions on future data, with the first type (Type I) assumes a upper and lower limit on the mass of dark matter from a particular reference model; and the second type (Type II) assumes a series of error bars at a collection of wavenumbers {k i } on a reference matter power spectrum P (k) (for example, see Refs.[111][112][113][114][115]).For convenience, we shall take the WDM model as the reference model in both cases.However, we emphasize that the WDM mass associated with the WDM model should simply be treated as a parameter that characterizes the deviation from CDM, and, in principle, any other dark-matter model can be used a reference.

Type I: analysis with a WDM mass range
Let us assume that a deviation from CDM is detected in future observations, and this deviation is best fitted by a reference WDM model with a WDM mass m WDM .We shall also assume that the relative uncertainty of the upper and lower limits in m WDM is the same such that It is tempting to apply both the δA and half-model analyses that we have shown in Sec.3.2.2 to find the upper and lower limits on m χ in our freeze-in and non-relativistic freezeout scenarios.However, while both methods are valid for recasting the lower bound, it is not clear how to treat the upper bound on an equal footing in the half-mode analysis.If one requires the squared transfer function associated with the upper bound to have a stronger power at any k < k 1/2 , the half-mode analysis can be overly conservative when σ m WDM is small for cases in which the test transfer function has a shape that differs significantly from the WDM transfer functions.Therefore, we shall only apply the δA method which is less sensitive to the detailed shape of the transfer function to recast the WDM bounds, i.e., the upper/lower limit on m WDM is mapped to the upper/lower limit on m χ in our scenarios by demanding them to produce the same δA.
We show in Fig. 12 the mapping between the allowed mass range for the freeze-in and the non-relativistic freeze-out scenarios as a function of m WDM with different choices of σ m WDM .Clearly, as long as σ WDM is not too large (σ WDM ≲ 0.4) the allowed regions for the nonrelativistic freeze-out and the freeze-in scenarios can be well separated as long as the reference WDM mass is not too small.For example, the two regions can be separated if m ref WDM ≳ 30 keV for σ m wdm = 0.4, and for σ m wdm = 0.2 or 0.1, they are completely separated for the entire range shown in the figure.Given the separation of the allowed regions, if additional information on the allowed range of dark-matter mass can be obtained from other sources, such as direct and indirect detection experiments, the constraints given by structure formation can be combined to distinguish between the freeze-in scenario and the non-relativistic freeze-out scenario.Such a comparison can in principle be done using different reference models and can be made between other scenarios.
It is worth pointing out that the δA analysis is designed to impose Lyman-α constraints.If the future constraints are derived primarily from other sources, such as the large-scale 21-cm data (see, e.g., Ref. [115]), one would need a different proxy other than δA to reliably recast the constraints.It therefore motivates us to explore the more complicated Type-II assumptions on future data which deal directly with data points on the matter power spectrum.

Type II: analysis with future constraints on P (k)
Future observations can be represented as constraints on the matter power spectrum (or transfer function) at a collection of wavenumbers {k i } which are consistent with a best-fitted matter power spectrum from a WDM model with a reference mass m ref WDM .Thus for any k ∈ {k i }, we have Obviously, we notice that the relative errors in the squared transfer function are always the same with that in the matter power spectrum since T 2 (k) = P (k)/P CDM (k).To specify the upper and lower limits, and the dependence on the scales k, we shall test three different possibilities on σ ± P (k): 1. Constant symmetric relative errors on P (k) In this case, we demand the relative errors of the upper and lower bounds are constant and symmetric with respect to P ref (k) on a linear scale.Thus, we can write σ ± P (k) = σ P to reduce the relative errors at different k to a constant.This method has the advantage that the error bars on either P (k) or T 2 (k) do not exceed below zero when the value of the WDM reference quantity is small.However, as we can see from Eq. ( 4.3), it tends to give more constraining power to error bars at the places where the reference P WDM (k) or T 2 WDM (k) is small and thus might be sensitive to largest wavenumbers at which measurements are made.

Constant symmetric absolute errors on log P (k)
Since the matter power spectrum is often presented on a logarithmic scale in literature, we also include the possibility that the error bars on P (k), when plotted on logarithmic scale, appear to be symmetric and constant.In this case, one can define a constant quantity with which it is easy to find that the relative errors on log P (k) and P (k), Notice that while σ log P (k) is symmetric with respect to log P WDM (k) but k-dependent, σ ± P is k-independent but asymmetric with respect to P WDM (k).This approach is convenient when σ P (k) can exceed O(1), especially when the upper/lower limit is a few times larger/smaller than the best-fit value.Nevertheless, the upper and lower limits might have unbalanced constraining power.It also tends to give more constraining power to the data points at which P WDM (k) or T 2 WDM is small.

Constant symmetric absolute errors on T 2 (k)
In this case, we demand that the error bars on the transfer function are constant and symmetric with respect to T 2 WDM (k).Similar to the second case, we can define a constant quantity Notice that the relative errors of P (k) are symmetric but k-dependentσ P (k) approaches a constant for small k, but grows drastically at large k.The advantages of this approach are that the constraining power is not dominated by the hypothetical data points at the largest k, and that it seems to treat the upper and lower bounds equally when comparing transfer functions as long as T 2 WDM (k) is not sufficiently small such that the lower limit T 2 − (k) becomes formally negative.In the case where the lower limit does turn negative beyond certain k, the lower bounds will lose their constraining power completely.
In addition to the form of the uncertainties in each case, we assume that the hypothetical future data consist of N data points evenly distributed along the log k-axis within the interval [k min , k max ].We then set the criterion that a test power spectrum is allowed only if P + (k i ) > P test (k i ) > P − (k i ) for any k i within the collection of data points at {k i }.For our work, the test power spectra are those obtained from the freeze-in and the non-relativistic freeze-out scenarios.In particular, since any data on P (k) can be turned into data on T 2 (k), with the help of Eq. ( 4.3), we shall base our analysis on the comparison between the squared transfer functions directly, using the fitting formulae in Eq. (3.16) and Eq.(3.18).
We present in Fig. 13 the allowed mass ranges for the freeze-in and the non-relativistic freeze-out scenarios for the three different assumptions on the uncertainties of the hypothetical future data.We fix the number of data points N = 20, and we have tested that further increasing N do not bring appreciable change in the allowed regions.Similar to the previous analysis on Fig. 12, we look for clear separation of the allowed mass ranges associate at the same reference WDM mass such that one is able to distinguish the two types of thermal history if additional information on dark-matter mass is available from other sources.
In the left panels, we fix k min = 1 h/Mpc and k min = 500 h/Mpc.We first notice in all cases that, for sufficiently large m ref WMD , while one can obtain the lower limits for the dark-matter mass m χ in corresponding scenarios, it is not possible to obtain the upper limits as the contours become vertical.This is not surprising since the resulting bounds on the nonrelativistic freeze-out and the freeze-in scenarios depend on both the observational precision and the scale at which P WDM (k) deviates from P CDM (k).In the case where the upper error bars on T 2 (k) all exceed one, i.e., the data set is consistent with the CDM predictions within the range of k, obtaining an upper limit for our scenarios is impossible since the free-streaming effects can only suppress the matter power spectrum.In other words, an upper limit on the dark-matter mass can only be found if a deviation from CDM predictions is confirmed in this analysis.
As the reference WDM mass decreases, the upper limits on m χ do appear as the corresponding contours slowly turn in the horizontal direction which indicates that the hypothetical data set starts to show deviation from CDM.If we decrease m ref WDM further, a distinct separation of the allowed mass ranges would eventually occur when both the uncertainty and are sufficiently small.
For the cases in the top and middle panels, this separation increases as the allowed mass range shrinks drastically when m ref WDM decreases.This is simply because for smaller m ref WDM , the part of the region where the transfer function is suppressed becomes larger within the fixed interval [k min , k max ], and thus the entire data set becomes more constraining.In these two cases, we find that the freeze-in scenario can even be ruled out completely if m ref WDM ≲ 15 keV, leaving the non-relativistic freeze-out scenario as the only viable possibility between the two scenarios.This happens when the uncertainties of the data set cannot accommodate any transfer function from the freeze-in scenario, and it is worth pointing out that this constraining power is not possessed by the Type-I data (which essentially also include the type of analysis based on the present-day average velocity of dark matter).However, it is more difficult to completely rule out the non-relativistic freeze-out scenario given that its corresponding transfer function has a shape that is extremely similar to that from a WDM scenario.Nevertheless, an exclusion of the freeze-out scenario is in principle possible if the reference scenario is modified.On the other hand, in the bottom left panel, the shrinking of the allowed region is not very noticeable at small m ref WDM .This is due to the fact that the absolute errors on T 2 WDM (k) is set to constant in this case, and thus the constraining power does not vary significantly when m ref WDM decreases.While varying m ref WDM would change the constraining power of the data set within a fixed range of k, we are also interested in understanding the how the range of k at which the data set is presented affects our analysis.In the right panels of Fig. 13, we show the mass ranges as a function of k max while holding m ref WDM = 10 keV fixed.We notice in all cases that for a given uncertainty there exists a critical value of k max at which the contours for the freeze-in and the non-relativistic freeze-out scenarios both become vertical, which means the upper limits for m χ cannot be obtained if k max is below this value.Similar to the situation in the left panels, the appearance of this critical value simply suggests that the data set cannot exclude CDM.We shall therefore call this critical value k CDM lim , and it can be determined by solving the equation Beyond k CDM lim , we also notice in the top and middle panels that for a given uncertainty there exist another critical value of k max , which we call k FI lim , such that the allowed mass region of the freeze-in scenario vanishes.As we have explained for the similar behavior in left panels, the appearance of k FI lim suggests that no transfer function from the freeze-in scenario can be fit within the error bars.Apparently, a data set that covers a range beyond k FI lim is able to rule out freeze-in completely, and if the uncertainty of the data decreases, k FI lim would also decrease making it easier to rule out the freeze-in scenario.Similar to k CDM lim , k FI lim also depends on m ref WDM since a change in it would shift the entire reference WDM transfer function.However, for both k CDM lim and k FI lim , the dependence on m ref WDM is somewhat degenerate with the dependence on the data uncertainty.
While the critical wavenumbers k CDM lim and k FI lim vary with the data uncertainty, the value of the transfer function at these critical points also depend on how precise these data are.Since we have been using the fitting functions Eq. (3.16) and Eq.(3.18) in our analysis, we need to make sure that the uncertainties of the fitting functions do not surpass that of the data.To examine this, we plot in the left panel of Fig. 14  For the relation between T 2 WDM (k FI lim ) and σ + P (k FI lim ), we first find that it does not depend on the reference WDM mass due to the translational symmetry in the transfer functions that we are considering (see Eq. (3.20)), and the fact that the most constraining data points for excluding the freeze-in scenario are always those close to k max .However, the relation at k FI lim is sensitive to the relative uncertainty σ + P (k CDM lim ) as T 2 WDM (k FI lim ) quickly decreases when the interaction only, it is therefore pressing to understand what one can infer from the effects that only rely on the gravitational interaction of dark matter.
In this paper, we explore the connection between structure formation and the darkmatter thermal history using the phase-space distribution after dark matter is produced.We first revisited the calculation to relate the present-day average velocity of dark matter to relevant physical quantities at the time when dark matter just decouples.We find that the decoupling temperatures of the dark and the SM sectors, as well as the mass of dark matter are correlated with the present-day velocities.This correlation can vary significantly across different types of dark-matter production mechanisms.We therefore re-interpret the Lyman-α constraints in the well studied WDM scenario as constraints on the present-day average velocity of dark matter and then apply them to other production scenarios such as freeze-in and non-relativistic freeze-out.We find that the constraints on the dark-matter mass can differ a lot among the two scenarios and the WDM scenario.After mapping out the correlation between the physical quantities at the early and the late times, we are able to project these preliminary bounds on dark-matter mass or average velocity to constrain the decoupling temperatures which concern the physics in the early universe.
While the present-day average velocity can provide a good estimate on the scale below which structure formation is suppressed, it cannot capture all the effects from the entire phasespace distribution in general.To provide more reliable constraints, we numerically solve the Boltzmann equation to obtain the phase-space distribution f χ after the production of dark matter.By requiring the results to match the present-day energy density of dark matter, we find that, above certain threshold mass m min χ , there are always two amplitudes for the same particle-physics processes that are able to generate the correct relic abundance.These solutions are associated with either the freeze-in regime if the interaction strength is feeble, the freeze-out regime if the particle interactions are rapid, or the transition regime if the interaction strength is intermediate.We then focus on pure freeze-in and the non-relativistic freeze-out scenarios, and use the corresponding phase-space distributions as input for the CLASS code to obtain the matter power spectra and the transfer functions.We find that the transfer functions in each scenario all possess the same shape and can be related to each other via a translational symmetry.In addition, all of them can be well fitted using the {α, β, γ} fitting formula for which we have introduced the dependence on the dark-matter mass m χ in the parameter α.The new fitting parameters then enable us to scan over a large range of the parameter space to constrain our scenarios.In particular, we have used the well-known half-mode analysis and the δA analysis to recast the current Lyman-α constraints on WDM as constraints in the freeze-in and the non-relativistic freeze-out scenarios.These constraints are weaker than the preliminary ones from the average velocity but still show significant large differences when comparing the results in the WDM, the freeze-in, and the non-relativistic freeze-out scenarios.
With the matter power spectrum, we have also discussed the effects on the physical quantities in the non-linear regime such as the halo mass function and the MW satellite counts using an extended Press-Schechter approach.The halos mass function is difficult to constrain directly from astronomical observations since dark-matter halos are not visible.Efforts on constraining the halo mass function thus rely on the visible content that resides in the halos and the comparison with simulations.For example, observables such as the galaxy cluster number counts, the galaxy cluster power spectrum, the lensing of Type Ia supernovae were exploited in Ref. [116], the possibility of using the magnification bias effect on highredshift submillimeter galaxies were discussed in Ref. [117].Ref. [118] used the cumulative luminosity function of galaxies at redshift z = 6 and placed a lower limit on WDM mass at 2.71 keV at 2σ level (or 2.27 keV for Sheth-Tormen [104] and 1.96 keV for Press-Schechter mass function), though this limit is sensitive on the minimum halo mass.Ref. [119] showed that the halo-mass-function parameters can be measured at a few percent level by weak lensing in future survey.For the MW satellite count, we found the constraints directly from it are weaker than that from the Lyman-α observations.On the other hand, recent studies (see e.g., Refs.[120][121][122][123]) on MW satellite count are shown to be able to place more stringent constraint on WDM mass.In particular, effects such as the tidal stripping from the host halo, the disruption of subhalo due to the Galactic disk, and reionization that were included in these works are able to lower the subhalo abundance and thus tighten the constraint.Nevertheless, we stick to the simple approach with Eq. (3.34) since it allows us to recast the WDM bounds in other scenarios analytically in a straightforward manner.Future observations with better sensitivity and larger sky coverage may further improve this bound.
Since no deviation from the CDM predictions has been confirmed thus far from current astrophysical or cosmological data, observational constraints can only provide lower bounds on the dark-matter mass in typical dark matter models.Although such bounds can be quite different in different dark-matter scenarios which have already enabled a discrimination between different dark-matter thermal histories to some extent, there is still a huge overlap in the allowed parameter space for these scenarios since no upper limit is present.We therefore seek to explore the extent to which future data can better improve the ability to distinguish different thermal histories.
As the form of future data is unknown to us, we examine two general types of hypothetical future data among which the Type-I data refer to upper and lower limits on the dark-matter mass of a reference model, and the Type-II data refer to a series of data points with error bars superposed on a reference matter power spectrum or transfer function.Obviously, the Type-I data should be seen as a fully processed result, whereas the Type-II data are closer to the raw data themselves though they are not the data on the actual physical observables.Using WDM as a reference scenario, and with several different assumptions on the data uncertainties, we find that in all cases one can significantly narrow down the allowed mass ranges for the freeze-in and the non-relativistic freeze-out scenarios -with reasonably small data uncertainties, the allowed regions for the two scenarios in consideration can be completely separated.In the presence of Type-II data, we find that it is even possible to fully exclude one of the scenarios (in our analysis the freeze-in scenario) while maintaining the viability of the other scenario.This ability, while requiring measurements with sufficiently small uncertainties to be made beyond certain critical wavenumber, is not possessed by the analysis with Type-I data which also includes the type of analysis based on the average velocity only.Our results therefore suggest that there is a great opportunity in distinguishing different thermal histories of dark matter from future observational data.
We first notice that the velocity distributions g v (v) ≡ g p (p) |d log p/d log v| ≈ g p (p) are identical to each other (after normalization by the area under the curve) up to some rigid translation along the log v-axis when the distributions only have non-vanishing support for non-relativistic velocities.In particular, changing m χ to m ′ χ amounts to translating the velocity distribution as On the other hand, the velocity distribution can be mapped into a distribution g k (k) in the space of wavenumbers by mapping the velocity v to the wavenumber k using the freestreaming horizon k hor [64,65]: In the limit v/a MRE ≪ 1, one can show numerically that since log(k ′ /k) is approximately a constant as v changes, the Jacobian J is also approximately a constant.Therefore, g ′ k (k ′ ) ≈ J(v)g v (v) = g k (k), i.e., the transformation in the k-space is also an approximately rigid translation.Approximating J as a true constant, we can write which means k ′ = ζk with ζ = (m ′ χ /m χ ) 1/J .The next step is to figure out how the squared transfer function transform.In Ref. [64], it has been shown that the logarithmic derivative of the transfer function at a wavenumber k is correlated with the fraction of dark-matter particles that are able to freestream a distance larger than ∼ 1/k.In particular, it is shown that the relation holds to high accuracy, where is called the hot-fraction function.Integrating Eq. (B.7), it is straightforward to show that where we have used the fact that lim k→0 log T 2 (k) = 0. Since F (k) is the integral of g k (k), it is easy to verify that the k-space transformation implies Therefore, the squared transfer function transforms as where we have used k = ζ −1 k ′ in the third line.The entire transformation then goes like Noticing that J ≈ 1.1 for k ∈ [1, 10 3 ] h/Mpc from numerical result, it is therefore reasonable to identify ζ = (m ′ χ /m χ ) 1/J with α(m χ )/α(m ′ χ ).Thus, we have approximately proved that the underlying reason for the translational symmetry is the invariance of g p (p) (which amounts to a rigid shift of g v (v) in the log v space) under the change of m χ as we have seen in Fig. 4.

FIG. 1 .
FIG.1.Relation between the dark-matter mass and present-day average velocity for different decoupling scenario, where we have chosen Ω χ = 0.25.The black horizontal lines indicate the critical velocities associated with the two WDM constraints m WDM ≥ 3.5 keV and m WDM ≥ 5.3 keV, and thus the part of the curve in the shaded region can be considered as excluded.Notice that, in the case of non-relativistic freeze-out, the curves cannot extend to arbitrarily small mass without violating the constraint on the dark-matter relic abundance.

FIG. 2 .
FIG.2.Constraints on decoupling temperatures in the cases of relativistic freeze-out (top left), non-relativistic freeze-out (top right), and freeze-in (bottom).The color map and the black dotted contours indicate the masses of dark matter.They are also contours of the present-day average velocity ⟨v⟩ 0 in the top left panel.On the other hand, ⟨v⟩ 0 is shown by the red dotted contours in the top right and bottom panels, where we have assumed that η dec = 1 in the freeze-in case.The red dashed and solid contours thus represent the current constraints on ⟨v⟩ 0 , with the arrow pointing to the constrained region of the parameter space.In addition, the blue contours in the top right panel shows the value of x dec .As described by the text, the gray areas indicate that, when using Eq.(2.21), either the solution of m χ is inconsistent with the assumption of the calculation (top left and bottom panels), or no solution can be found in this region (right panel).

FIG. 3 .
FIG.3.Contours of Ωχ = 0.25 in the m χ -|M| 2 plane for different values of η dec .For given a m χ , the (constant) squared amplitude |M| 2 is determined by solving the Boltzmann equation to obtain Ωχ = 0.25.As indicated by the color bar, the interaction rate γ χ varies from γ χ ≪ 1 to γ χ ≫ 1 along each contour.Thus, each contour corresponds to a path that goes through the freeze-in, transition, and freeze-out regimes.The diamond markers are placed at γ χ = 0.1 and 10 to roughly identify these three regimes.

FIG. 4 .
FIG. 4. Left panel:The evolution of the comoving mass density in the freeze-out (solid) and the freeze-in (dashed) scenarios, and when dark matter is in thermal equilibrium with the SM sector (dotted).Note that the mass m χ is different for each curve, and so is the corresponding temperature along the horizontal axis.Right panel: Normalized momentum distributions of dark matter today associated with the examples shown on the left.The freeze-out and the freeze-in scenarios are represented by the solid and the dashed curves, respectively.

FIG. 5 .
FIG. 5. Left panel:The evolution of the comoving mass density in the transition regime where, at the same dark-matter mass, the solution with a larger/smaller amplitude is represented by the solid/dashed curve.The evolution of m χ N χ when χ is kept in equilibrium with the thermal bath is shown by the dotted curves.Right panel: The resulting phase-space distributions of dark matter associated with the examples in the left panel.All the distributions are normalized by the their comoving number density N χ .

FIG. 6 .
FIG. 6. Evolution of the density perturbation modes |δ k | with k = 4.3, 14.2 and 42.8 h/Mpc for the CDM, freeze-out and freeze-in scenarios.For the latter two scenarios, three darkmatter mass, m χ = 10, 20 and 30 keV are chosen.The vertical dashed line marks the matter-radiation equality. )

FIG. 7 .
FIG.7.The squared transfer functions resulting from the freeze-out (left panel) and freezein (right panel) scenarios for a series of dark-matter masses from m χ = 10 keV to 100 keV.

m χ ≳ 21 . 9 (
FIG. 8. Mappings between the constraint on m WDM and that on m χ based on the half-mode analysis (solid curves) or the δA analysis (dashed curves) in the freeze-in (blue) and the nonrelativistic freeze-out (red) scenarios.Notice that the two curves overlap completely in the latter scenario.As indicated by the black lines which show the current bounds, a Lyman-α constraint on the WDM mass can be directly mapped to a constraint on m χ in a particular scenario.

FIG. 9 .
FIG.9.The halo mass functions for non-relativistic freeze-out (left panel) and freeze-in (right panel) with dark-matter mass m χ ranging from 10 to 100 keV.

1 FIG. 12 .
FIG.12.Mapping between the resulting allowed range of m WDM to the allowed range of m χ in the freeze-in (shaded in red) and the non-relativistic freeze-out (shaded in blue) scenarios based on the δA analysis.The solid, dashed and dotted contours correspond to the choices σ WDM = 0.4, 0.2 and 0.1, respectively.

1 FIG. 13 .
FIG.13.Mappings from the hypothetical constraints on the WDM matter power spectra to the allowed mass ranges for the freeze-in (blue shaded) and the non-relativistic freeze-out (red shaded) scenarios.The results in the top, middle and bottom panels are obtained by assuming constant symmetric relative errors on P (k), constant symmetric absolute errors on log P (k) and constant symmetric absolute errors on T 2 (k), respectively.Contours with different line styles correspond to specific choices of σ P , ∆ log P or ∆ T 2 which characterize the errors on hypothetical data.In the left panels, we fix k min = 1 h/Mpc and k max = 500 h/Mpc and assume N = 20 hypothetical data points evenly distributed on log scale.In the right panels, we fix m ref WDM = 10 keV and vary k max .

2 γv γ 2 v 2 + a 2 − 1 , (B. 3 )
where ξ is an O(1) constant factor.One can then define a distribution function in the log k spaceg k (k) ≡ J(v)g v (v) , (B.4)in which J(v) ≡ |d log v/d log k| is the Jacobian.With these, the translation in the v-space can be mapped to a translation in the k-space which goes like log k ′ − log k = (log v − log v ′ )/J = log m ′ the relation between T 2 WDM (k lim ) and σ + P (k lim ), where k lim can be either the one for CDM or for freeze-in.The relation at k CDM lim can be obtained directly from Eq. (4.8) which clearly shows that it does not depend on which kind of data set we use, nor does it depend on the reference WDM mass.It is also clear that T 2 WDM (k CDM lim ) ≳ 0.5 throughout the entire region and does not decrease significantly as long as σ + P (k CDM lim ) is not much bigger than O(1).Therefore, the location ofk CDM lim can be reliably identified as the fitting function has a high accuracy in this regime with the relative error less than 5%.