Precise dark matter relic abundance in decoupled sectors

Dark matter (DM) as a thermal relic of the primordial plasma is increasingly pressured by direct and indirect searches, while the same production mechanism in a decoupled sector is much less constrained. We extend the standard treatment of the freeze-out process to such scenarios and perform precision calculations of the annihilation cross-section required to match the observed DM abundance. We demonstrate that the difference to the canonical value of this 'thermal cross-section' is generally sizeable, and can reach orders of magnitude. Our results directly impact the interpretation of DM searches in hidden sector scenarios.

Introduction.-Cosmological observations require the existence of a dark matter (DM) component that makes up about 80 % of the matter in our Universe [1] and likely consists of a new type of elementary particle [2,3]. The most often adopted paradigm for DM production is via freeze-out from the primordial plasma of standard model (SM) particles [4]. This roughly requires weak-scale couplings for DM masses at the electroweak scale -which has been argued to be an intriguing coincidence in view of proposed solutions to the hierarchy problem of the SM [5] -but the same mechanism also works for lighter DM and correspondingly weaker couplings [6]. The formalism to calculate the thermal relic abundance in these scenarios [7,8] is well established and successfully used in a plethora of applications, e.g. for benchmarking the reach of experimental searches for non-gravitational DM interactions [1,[9][10][11][12][13]. Based on this standard prescription, several public numerical codes [14][15][16][17] provide precision calculations of the DM abundance, matching the percent level observational accuracy.
More recently, the focus has shifted to models where DM couples more strongly to particles in a 'secluded' dark sector (DS) than to the SM [18][19][20][21][22]. This development is partially motivated by the fact that more traditional DM candidates are increasingly pressured by the absence of undisputed signals in direct searches as well as at colliders [23][24][25], but also from a theoretical perspective there is no need for sizeable inter-sector couplings. Remarkably, thermal freeze-out works equally well also in these models, providing a compelling potential explanation for the observed DM abundance. As a consequence, couplings needed to achieve this goal are often either implicitly fixed or explicitly targeted in various searches for hidden sector particles [26][27][28][29][30][31][32][33].
Despite this development, relic density calculations in such scenarios have not yet reached the same level of refinement as for thermal freeze-out in the visible sector. As a critical first step towards bridging this gap, we per-form here a concise and comprehensive analysis of modelindependent aspects of such calculations, matching both in spirit and precision the widely adopted treatment of 'standard' freeze-out in the visible sector [34].
In the first part of this Letter we discuss (effectively) massless annihilation products. Since in this case the main differences to the standard analysis have largely been identified previously [19,[35][36][37][38][39][40][41][42][43][44], our focus here is on a pedagogic and easily accessible summary, including the presentation of new benchmark 'thermal' cross-sections that can directly be compared to the corresponding visible sector results [34]. In the second, and main, part of this Letter we then perform a detailed analysis of DM annihilating into DS particles similar in mass. In this case comoving conservation of the total number density necessitates an accurate treatment of chemical potentials of all involved particles, both before and during freezeout, which we provide here for the first time.
Standard freeze-out.-We start by briefly revisiting the canonical approach. The number density n i of DM particles i = χ,χ initially in thermal equilibrium with the SM heat bath at temperature T can be described by the Boltzmann equation [7] dn i dt + 3Hn i = σv (n χ,eq nχ ,eq − n χ nχ) , where H is the Hubble rate, and n χ,eq = nχ ,eq = g χ m 3 Here, x ≡ m χ /T , K j are modified Bessel functions of order j, g χ denotes the internal degrees of freedom (d.o.f.) of χ, σ is the total cross-section for DM annihilations, for a centerof-mass energy √ s ≡ 2m χ √s , and v lab is the velocity of one of the DM particles in the rest-frame of the other.
Let us stress two main assumptions that enter in this widely used form of the Boltzmann equation. The first  is that the DM phase-space distribution is of the form f i ∝ f χ,eq = exp(−E χ /T ), i.e. that the freeze-out happens for m χ T and well before kinetic decoupling (see Ref. [45] for a treatment of early kinetic decoupling). The second assumption is that the annihilation products constitute a heat bath, in the sense that none of them builds up significant chemical potentials. Crucially, both assumptions can be violated in decoupled sectors.
In Fig. 1 we indicate with solid lines the value of σv in the standard scenario (assuming a constant value of this quantity around chemical decoupling) that is needed to obtain a relic density matching the observed cosmological DM abundance of Ω DM h 2 = 0.12 [1]. The orange solid lines show the case of Majorana DM (with g χ = 2 and Ω χ = Ωχ = Ω DM ), updating the conventionally quoted 'thermal relic cross-section' in Ref. [34] with a more recent measurement of Ω DM and recent lattice QCD results for the evolution of SM d.o.f. in the early Universe [46]. For comparison, the blue lines indicate the case of Dirac DM (g χ = gχ = 2 and Ω χ = Ωχ = Ω DM /2) to stress the not typically appreciated fact that the required value of σv is not exactly twice as large as in the Majorana case.
A secluded dark sector.-The idea [18][19][20][21][22]26] that DM could be interacting only relatively weakly with the SM, but much more strongly with itself or other particles in a secluded DS, has received significant atten-tion [29,36,[47][48][49][50][51][52]. In such scenarios, both sectors may well have been in thermal contact at high temperatures, until they decoupled at a temperature T dec . The separate conservation of entropy in both sectors then implies a non-trivial evolution of the temperature ratio, where g SM,DS * denotes the effective number of relativistic entropy d.o.f. in the two sectors. Let us stress that this commonly used relation tacitly assumes that DM is in full equilibrium with at least one species S with vanishing chemical potential, µ S = 0 (implying µ χ = −µχ as long as DM is in chemical equilibrium).
For a precise description of the freeze-out process of χ in such a situation the standard Boltzmann equation (1) needs to be adapted at three places: both i) the equilibrium density n eq and ii) the thermal average σv must be evaluated at T χ rather than the SM temperature T , and iii) the Hubble rate must be increased to take into account the energy content of the DS. During radiation domination, in particular, this means that f g f )ξ 4 and the sums runs over the internal d.o.f. of all fully relativistic DS bosons (b) and DS fermions (f ) (in our numerical treatment, we always use the full expression for g eff ). We note that existing relic density calculations for decoupled DSs very often only take into account a subset of these effects, or implement them in a simplified, not fully self-consistent way.
Model setup.-Let us for concreteness consider a simple setup where the DS consists of massive fermions χ, acting as DM, and massless scalars S with µ S = 0, constituting the heat bath. We assume that the DS decoupled from the SM at high temperatures, such that g SM * (T dec ) = 106.75 and g DS * (T dec ) = g S + (7/4)N χ in Eq. (3), where N χ = 1 (2) for Majorana (Dirac) DM. In Fig. 1 we show the 'thermal' annihilation cross-section for χχ → SS in such a scenario, for different values of g S . The fact that this differs significantly from the standard case illustrates the importance of including the effects outlined above in a consistent way. In this sense, the updated procedure for relic density calculations directly impacts a large number of DS models where annihilation also proceeds via an s-wave [20,26,29,47,[53][54][55][56] even though σv is often velocity-dependent in these cases, impeding a literal interpretation of the curves shown in Fig. 1. In order to facilitate the study of such more realistic scenarios, we have updated the general-purpose relic density routines of DarkSUSY [14] to perform precision calculations of DS freeze-out that self-consistently take into account all three effects discussed above. This allows to consider a broad range of relevant models with in principle arbitrary amplitudes, including p-wave an-nihilation (see also Appendix A) and, e.g., Sommerfeld enhancement. 1 To understand the behavior of the curves in Fig. 1, we first note that a constant σv (as in this specific example) is of course not affected by a change in ξ. For g S = 1, furthermore, the change in g eff and hence the Hubble rate has only a subdominant effect (but becomes somewhat more important for g S = 5). The main effect visible in the figure thus originates from changing n χ,eq (x) → n χ,eq (x/ξ). For large DM masses and hence freeze-out temperatures, in particular, the heating in the DS due to χχ → SS, cf. the nominator of Eq. (3), is more efficient than the heating in the SM, leading to ξ > 1 around freeze-out. This leads to a larger DM density, at a given SM temperature T , which has to be compensated for by a larger σv to match the observed relic abundance. Below DM masses of a few GeV, the drop in the SM d.o.f. until freeze-out is more significant than that in the DS (especially during the QCD phase transition), leading to ξ < 1 and hence the need for a smaller value of σv compared to the standard case represented by the solid lines. We finally note that the energy density of S at late times is independent of m χ . Expressing it in terms of an effective number of relativistic neutrino species, this corresponds to ∆N eff = 0.104(0.202) for Majorana DM with g S = 1(5), and ∆N eff = 0.201(0.275) for Dirac DM -which is below current CMB bounds on this quantity, ∆N eff < 0.29 (95% C.L.) [1], but within reach of next-generation CMB experiments [57,58].
Chemical potentials during freeze-out.-The above treatment still assumes that the annihilation products are in chemical equilibrium with themselves during the entire chemical decoupling process of DM. This is consistent for massless DS particles S, where interactions such as χχ → χχS (or number-changing reactions purely within the S sector) will always enforce µ S = 0. For a fully decoupled DS only containing massive degrees of freedom, however, this is no longer necessarily the case. Largely independent of the concrete model realization, in particular, number-changing interactions of massive particles S will cease to be efficient at the latest when the DS temperature drops below their mass, T χ m S . This implies that all DS particles will generally build up chemical potentials before and during the freeze-out process (and not only the DM particles, as in the standard scenario).
Let us for illustration consider the same setup as before, but mostly focus on DS particles χ and S close in mass. In chemical equilibrium, χχ ↔ SS then enforces   (1) assuming thermal equilibrium of S with an additional massless DS heat bath particle. All curves are based on the same Mχχ→SS 2 = const., adjusted to give the correct relic density in the limit mS → 0. furthermore has µ S = 0; this implies µ χ = −µχ = 0 (where the last equality assumes a vanishing asymmetric DM component). When S later develops a non-vanishing chemical potential, on the other hand, that initial condition leads to µ χ = µχ = µ S > 0. Subsequently, DM decouples chemically from S -but will typically remain in kinetic equilibrium at least until the end of the freezeout process. The phase-space densities of all DS particles are then still given by Fermi-Dirac or Bose-Einstein distributions with temperature T χ and chemical potentials µ χ = µχ = µ S [59]. We recall that Eq. (1) describes the evolution of n χ if µ S = 0 and the effect of quantum statistics can be neglected, i.e. f χ ∝ exp(−E χ /T ). In the following we will instead demonstrate how to accurately determine the evolution of n χ without these two assumptions.
Let us start, for simplicity, with the case where χχ ↔ SS is the only relevant S-number changing reaction, implying, e.g., that S is sufficiently long-lived (below, we will modify this assumption). We then consider the Boltzmann equations for the number densities, where C is the integrated collision operator in its standard form [59], as well as energy conservation in the DS during freeze-out, ∇ µ T 0µ DS = 0. The latter takes the forṁ with total energy density ρ DS ≡ N χ ρ χ + ρ S and pressure P DS ≡ N χ P χ + P S . Finally, as long as χ and S stay in kinetic equilibrium, all cosmological quantities Q ∈ {n a , ρ a , P a | a ∈ {χ,χ, S}} can be interpreted as functions of T χ , µ χ and µ S only. We therefore can usė to transform Eqs. (4) and (5) into a set of differential equations for T χ , µ χ and µ S , which we solve numerically (with µ χ = µχ = µ S = 0 as initial condition). Note that Eqs. (4) and (5) replace Eqs. (1) and (3), and generally only imply entropy conservation during chemical equilibrium. For the specific benchmarks below, chemical decoupling occurs only after significant Boltzmann suppression of χ andχ, such that the respective change in total DS entropy can still be neglected. In Fig. 2 we demonstrate the resulting evolution of the particle abundances Y ≡ n/s, with s the total entropy density in the SM and DS. For definiteness we choose a Majorana DM particle with m χ = 100 GeV and a constant annihilation amplitude (for which we provide a closed expression for C in Appendix B) that would result in the correct relic density in the standard treatment (translating to a value of σv Tχ→0 about 10% larger than the orange lines in Fig. 1). The red curves show the case of m S = 0 for which, following the discussion above, we explicitly set µ S = 0. The resulting evolution of χ (red solid line) therefore coincides exactly with the result of the standard treatment of solving Eq. (1). We note that the increase in Y S around T χ ∼ m χ is due to the Boltzmann suppression of χ, analogous to the increase in n γ /s during e + e − annihilation in the SM.
For more degenerate masses (green and purple lines in Fig. 2), we allow all chemical potentials to evolve freely. This leads to a rise in µ S , compensating the would-be Boltzmann suppression of S, and an asymptotic abun- The greater number of S particles then delays the Boltzmann suppression of n χ from around T χ ∼ m χ to when the mean kinetic energy of S drops below δm, roughly around T χ ∼ δm. For reference we also show an application of Eq. (1) (dotted lines) assuming thermal equilibrium of S with additional massless DS heat bath particles such that µ S = 0 and T χ ∝ a −1 with the scale factor a. Comparing the purple lines (δm/m χ = 10 −2 ), e.g., Boltzmann suppression of χ for the solid line occurs at temperatures T χ around two orders of magnitude smaller than for the dotted line, or a one order of magnitude larger (T χ ∝ a −2 at T χ m S for the solid line). Approximating the annihilation rate by σv n χ ∝ a −3 , whereas the dilution by cosmic expansion is H ∝ a −2 , this im- plies that freeze-out happens when χ is less Boltzmannsuppressed and Y χ is enhanced by ∼ a, i.e. around one order of magnitude. In general, the correct treatment of the chemical potentials thus leads to an enhanced DM abundance compared to the 'naïve' assumption of µ S = 0 and T χ ∝ a −1 . Comparing instead to the m S = 0 case, cf. the standard situation depicted in Fig. 1, Y final χ first decreases up to a mass ratio of m S /m χ = 0.4 (green lines) as a larger m S implies a faster decrease of T χ with time around T χ ∼ m χ such that, in fact, the SM temperature T is somewhat larger around freeze-out. Approximating the annihilation rate as above, freeze-out occurs when H/s ∼ σv n χ /s ∝ 1/T (the SM dominates the energy and entropy densities), leading to a slight decrease in Y χ = n χ /s. For even larger m S , this effect is compensated by the delayed Boltzmann suppression of χ and the DM abundance increases as S and χ become more and more degenerate.
For S close in mass to χ, the final DM relic abundance will not only depend on the decoupling process but also on how S decays after freeze-out. If S was stable, in particular, it would simply contribute to the total DM density, by far overshooting the observed value (unless allowing for sufficiently small temperature ratios ξ T →∞ 1, thus relaxing our assumption of initial thermal contact between SM and DS). In the following we assume a lifetime τ S of S such that the decays occur only after freeze-out, implying a negligible impact on the freeze-out process itself. In Fig. 3 we explore two concrete decay scenarios, by showing the 'thermal' annihilation cross-section for the same mass ratios as discussed in Fig. 2. 2 The first scenario is S decaying to effectively massless DS states, or dark radiation (DR), and indicated by solid lines. The additional effective relativistic d.o.f. resulting from the decay of S will in general depend on the lifetime τ S , because the energy densities of matter and radiation red-shift differently. As already for m S = 0 one has ∆N eff = 0.104 (see above), the case m S ∼ m χ is generally in conflict with the CMB limit even if the decay happens shortly after freeze-out. The second example (dash-dotted lines) considers S decays to SM states. In this case, the resulting entropy injection into the SM plasma will lead to a dilution of DM, lowering the required DM annihilation cross-section. This effect has recently been argued to allow for DM masses above the naïve unitarity limit [38,60,61]. Note that the lifetime τ S = 1 s × (1 GeV/m S ) 2 chosen here for illustration is expected to be in conflict with observations of primordial element abundances for τ S > 0.1 s [62], i.e. m S 3 GeV.
To summarize, the solid lines in Fig. 3 show the required DM annihilation cross-section to obtain the observed DM abundance assuming S decays without injecting entropy in the SM and thus diluting the DM abundance. These lines therefore provide an upper limit to scenarios where S decays into the SM after DM freezeout, as exemplary illustrated by the dash-dotted lines. It is evident that the required DM annihilation crosssection can be very different from the canonical value shown in Fig. 1, in particular for small mass differences. In the extreme case of degenerate masses, m S = m χ , no Boltzmann suppression of χ can occur -independently of the DM annihilation cross-section -implying that the observed DM abundance can only be achieved for sufficiently small temperature ratios ξ T →∞ 1 as discussed above for a stable S.
Discussion.-For the choice of parameters discussed above we explicitly checked (see Appendix B) that the assumption of kinetic equilibrium is always satisfied during the freeze-out process, justifying our ansatz for the phase-space distributions f a . Let us stress that this is particularly important for small mass splittings, where µ S m S makes it mandatory to include the full quantum statistics for all particles when aiming for precision calculations of the relic density. The commonly used assumption of a Maxwell-Boltzmann distribution is, in other words, no longer justified. For the green lines in Fig. 3, for example, with m S = 0.4 m χ , ignoring the effect of quantum statistics leads to an overestimate of the final DM abundance by about 3 % for hidden sector decay (solid line) and more than 10 % for decay into SM particles (dash-dotted lines). While we leave a more detailed investigation for future work, let us stress that the effect is typically larger during the freeze-out process, at higher temperatures, and hence potentially more relevant for (semi-)relativistic freeze-out; for hidden sector decay, furthermore, the impact is also larger on the abundance of S, which is very sensitive to ever more stringent constraints from ∆N eff .
So far we have focussed on a fully secluded DS, in which case the most prominent observables to test such models are Ω DM and ∆N eff . It is however worth mentioning that in many models there are additional tiny couplings to the SM that would allow further experimental signatures. A setup where hidden sector freeze-out can naturally occur while still allowing for sufficiently large couplings to the SM to be probed by particle physics experiments, e.g., are scalar or pseudoscalar mediators with Yukawalike coupling structure [30,32,33,[63][64][65]. Also indirect DM searches for secluded dark sectors [66] provide a potentially promising avenue, in particular for the strongly enhanced annihilation rates necessary to accommodate DM degenerate in mass with its annihilation products.
Conclusions.-In this work we have presented a framework for precision calculations of DM freeze-out in a secluded sector, matching the observational accuracy on the one hand, and the increasing demand for consistent interpretations of phenomenological dark sector studies on the other hand. We have provided new benchmark 'thermal' annihilation cross-sections for relativistic heat bath particles, and demonstrated that the difference to the standard treatment can be even larger for non-relativistic DM annihilation products. The latter case is intrinsically strongly model-dependent, and will be studied in more detail elsewhere. Further interesting extensions, not the least in view of the significant model-building activity in these areas, would be to generalize the precision relic calculations presented here to models where the DM particles in the hidden sector do not obey a Z 2 symmetry [67][68][69], are asymmetric [70] or have a relic abundance set by freeze-in rather than freeze-out [71,72].

A. DM annihilation via p-wave
In the case of s-wave annihilation to massless final states, the velocity-weighted annihilation cross-section is constant in the limit of small DM velocities, resulting in σv = σv lab . This simplified ansatz for σv (neglecting higher-order contributions in v, following common practice) has been presented in Fig. 1 in the main text, both for DM annihilating to SM particles and for situations in which the relic abundance is set via freeze-out in a hidden sector.
Here we complement this by considering instead the case of p-wave annihilation, which also has been frequently considered for DS freeze-out production of DM [29,30,32,55,[63][64][65]. To describe such models, we will again take a simplified ansatz for the cross-section by only keeping the leading term in the DM velocities, where we assume b to be constant. For the thermally averaged cross-section entering in the Boltzmann equation, Eq. (1) in the main text, this implies σv = b× 6(x/ξ) −1 − 27(x/ξ) −2 + ... . The value of b resulting in the correct DM relic abundance in this case is shown in Fig. 4, for the same choice of DM models (Dirac and Majorana fermions, respectively) and heat bath components as in Fig. 1 in the main text. In comparison, the main differences in these figures are that i) the value of b resulting in the correct relic density is about one order of magnitude larger than the value of σv required in the case of s-wave annihilation and that ii) this 'thermal' value of b rises faster with m χ than its s-wave counterpart. Both of this can be traced back to the fact that also for p-wave annihilation it is σv around chemical decoupling, and not b, that sets the relic density. In the SM case, e.g., b/ σv ≈ x cd /6, where x cd depends logarithmically on the DM mass and rises from x cd ≈ 18.8 (for m χ = 100 MeV) to x cd ≈ 31.6 (for m χ = 100 TeV). The above estimate should be corrected by another factor of about 2 because decoupling does not happen instantaneously, and dT σv p−wave / dT σv s−wave ≈ 1/2 (as first stressed in Ref. [73]). The same general trend, finally, is also visible for annihilations in the hidden sector, with ξ = 1. Compared to Fig. 1 in the main text, furthermore, the difference between SM and DS results is somewhat larger because ξ enters directly in σv .

B. Collision term including chemical potentials
For general two-body annihilation processes χχ ↔ SS , and assuming CP -invariance, the integrated collision operator from Eq. (4) in the main text takes the form where dΠ a = d 3 p a /(2π) 3 2E a , integration is implied over all (not only physically distinct) momentum configurations, and |M χχ→SS | 2 is the squared matrix element, averaged (summed) over the spins and other internal degrees of freedom of all initial (final) state particles. We assume all involved particles to be in kinetic equilibrium, i.e. the phase-space distributions take the form f a = 1/[e (Ea−µa)/Tχ ± 1], with a ∈ {χ,χ, S, S } and the − (+) sign is used for bosons (fermions). In the special case of a constant matrix element -which is justified for contact-like interactions and which we adopt as benchmark scenario in the main text -Eq. (8) can be simplified to Moreover, α ≡ T χ p log e E /Tχ − e (E −p β+2µ S )/(2Tχ) e E /Tχ − e (E +p β+2µ S )/(2Tχ) , (12) where p ≡ | p χ + pχ| = (p 2 χ + p 2 χ + 2p χ pχ cos θ) 1/2 and E ≡ E χ + Eχ.
For highly non-relativistic DM, the annihilation crosssection for a constant matrix element becomes independent of the center-of-mass energy, and hence σv lab σv . In this limit, annihilation cross-section and amplitude are related as In the simplest models, the same constant matrix element also describes the scattering process χS ↔ χS, in which case the above expression provides a convenient means of estimating the time of kinetic decoupling for a given value of σv . For m χ ∼ m S , e.g., this happens when the scattering rate falls behind the Hubble rate, n S σ χS↔χS v ∼ H, while for m χ m S it is instead the (smaller) momentum exchange rate γ that provides the relevant scale (see, e.g., Refs. [45,74]). Using this condition, we explicitly checked that S and χ remain in kinetic equilibrium during the freeze-out process.