Dark Freeze-out Cogenesis

We propose a new mechanism where a multi-component dark sector generates the observed dark matter abundance and baryon asymmetry and thus addresses the coincidence between the two. The thermal freeze-out of dark matter annihilating into meta-stable dark partners sets the dark matter relic abundance while providing the out-of-equilibrium condition for baryogenesis. The meta-stable state triggers baryon asymmetry production by its decay well after the freeze-out and potentially induces a period of early matter domination before its decay. The dark matter and baryon abundances are related through number conservation within the dark sector (cogenesis). The"coincidence"is a natural outcome with GeV- to TeV-scale symmetric dark matter and the dark sector's interactions with the Standard Model quarks. We present a UV-complete model and explore its phenomenological predictions, including dark matter direct detection signals, LHC signatures of new massive particles with color charges and long-lived particles with displaced vertices, dark matter-induced nucleon conversions, (exotic) dark matter indirect detection signals, and effects on the cosmological matter power spectrum. As a side result, we provide a novel analytical treatment for dark sector freeze-out, which may prove useful in the study of related scenarios.


Introduction
The nature and dynamical origins of the baryon asymmetry and the observed dark matter (DM) are two long-standing puzzles in particle cosmology. Meanwhile, the observation that their abundances are strikingly similar, Ω DM /Ω B ≈ 5 [1], presents a coincidence problem which is suggestive of a potential connection between the origins of DM and baryons in the early Universe. Although there exists a cornucopia of theoretical explanations which address each of these pieces separately, this cosmic coincidence has inspired new directions in DM model-building (see e.g., [2,3] for reviews).
Weakly interacting massive particles (WIMPs) as DM candidates have been leading the DM model-building efforts since the 1990s [4,5]. This paradigm is motivated largely by the observation that DM with weak-scale interactions and masses can produce the correct DM abundance through thermal freeze-out (WIMP "miracle"). However, conventional WIMPs have become increasingly constrained by indirect/direct detection, and collider experiments, see, e.g., [6,7] and references therein. This has led to the exploration of alternative DM candidates beyond the WIMP paradigm. For example, asymmetric dark matter (ADM) [8][9][10][11][12][13][14] is an alternative inspired by the DM-baryon coincidence. In this framework, the DM particle is distinct from its antiparticle, and an asymmetry in their respective population densities is generated in the early universe. The core idea of ADM is based on relating and generating DM and baryon/lepton asymmetries through shared interactions. With the exception of a few mechanisms (e.g. [15][16][17][18]) the observed coincidence is achieved with ADM masses of O(GeV).
Recently, there have been attempts at unifying WIMP DM and ADM mechanisms [3,[19][20][21][22][23][24][25][26][27]. Among the existing WIMP-related proposals, [20] is highly sensitive to various initial conditions, while both [21] and WIMP DM annihilation triggered "WIMPy baryogensis" [19] have sensitivity to "washout", i.e., a reduction in the baryon asymmetry due to inverse decays or L-violating scatterings. The mechanism of "Baryogenesis from metastable WIMPs" [22] was then proposed as an alternative where the prediction is more robust against model details: the baryon asymmetry is generated by a long-lived WIMP that undergoes CP-and B-violating decays after its thermal freeze-out. Such models also provide a strong cosmological motivation for long-lived particle searches at collider experiments and have become a benchmark for related studies [28][29][30]. However, in its original form, it does not involve the specifics of DM, except for assuming it is WIMP-like. The more recently proposed "WIMP Cogenesis" [27] is a new realization of ADM with DM specifics explicitly incorporated, and fully inherits the desirable qualities of both WIMP and ADM production mechanisms. However, the permitted DM mass which produces the observed DM-baryon coincidence is limited to a few GeV.
From the model building perspective it is desirable to further develop a framework which incorporates the merits of WIMP-like and ADM-like mechanisms. Specifically, a mechanism which maintains a symmetric DM candidate permitting a wider range of DM masses, O(GeV − TeV), while supplying a predictive connection between Ω DM and Ω B presents a worthwhile theory target. In this work, we consider a realization of this goal in the framework of an isolated dark sector which may have a thermal temperature different from the Standard Model (SM) photon bath (see e.g., [31][32][33][34][35][36][37]). In particular, the dark sector is composed of a stable state as DM and a meta-stable dark partner which is the "parent" of the baryon asymmetry. Both states are initially in equilibrium with each other through the efficient annihilation of DM into dark partners, before freezing out. The dark partners are long-lived, and their decays into SM particles trigger the generation of baryon asymmetry. The stable and meta-stable states, and consequently, Ω DM and Ω B , are connected through the overall number-conservation in the dark sector as manifest in the annihilation process.
Independent of the baryogenesis motivation, we note that the particular dark freezeout scenario we consider in this work is new to the literature. While there are existing studies on freeze-out in a decoupled dark sector [35,[37][38][39][40], some also involving meta-stable annihilation final states, the present scenario distinguishes itself from known possibilities. For instance, in the Co-Decaying dark matter case [35], the meta-stable states decay during the freeze-out process which exponentially depletes the DM abundance, unlike in our case. We provide both numerical and analytical approaches for tracking the dynamical evolution during the freeze-out going beyond the instantaneous freeze-out approximation.
The basic idea of Dark freeze-out Cogenesis is illustrated in the schematic of Fig. 1. As we will elaborate later in the paper, a simple realization consists of three dark states χ 1,2,3 , where χ 1 is taken to be the stable DM candidate and χ 2 is the meta-stable component whose baryon number-and CP-violating decays trigger baryogenesis. Departure from equilibrium is achieved when χ 1 and χ 2 freeze-out in the isolated sector. Finally, χ 3 is required to induce CP-violation. When the meta-stable χ 2 component undergoes CP-and baryon number-violating decays, its abundance is converted into the baryon asymmetry. The latter is directly tied to the stable DM component's abundance through number conservation in the dark sector. We demonstrate our general ideas by introducing a UV complete model, i.e., a particle content and their interactions solely based on renormalizable interactions.
The remainder of the paper is organized as follows: in Sec. 2 we outline two scenarios of dark sector freeze-out and describe the dynamics of the cosmological evolution of matter abundances around the dark freeze-out. In Sec. 3, we discuss the dynamics related to the subsequent χ 2 decays and the ensuing production of the baryon asymmetry. A UV-complete particle physics model is introduced in Sec. 4. Various phenomenological implications and observational signatures predicted by the model are detailed in Sec. 5 before concluding in Sec. 6. Further details on the freeze-out calculation are relegated to App. A.

Dark Sector Freeze-out
In this section, we consider a dark sector with states χ 1 and χ 2 that evolve independently from the thermal SM bath. This is commensurate with initial conditions where χ 1,2 decouple from SM while still being relativistic. The state χ 1 is stabilized by a Z 2 symmetry, so that it plays the role of the DM candidate. In contrast, a lighter state χ 2 is meta-stable, and decays out-of-equilibrium after χ 1 freezes out, leading to the observed baryon asymmetry (once interference effects with a state χ 3 , to be introduced later, are taken into account). The dark sector is assumed to be self-thermalized before freeze-out, characterised by a common dark temperature T , which is set to be comparable to or lower than the photon temperature T , T ≤ T , so that the SM sector always dominates the radiation energy of the Universe.
Depending on the mass hierarchy between χ 1 and χ 2 , there are various options for the details of DM freeze-out. In the following, we focus on two principal possibilities, for each of which we present analytical solutions: m χ 1 m χ 2 ("hierarchical scenario") and m χ 1 = m χ 2 /(1 − δ) with δ 1 ("nearly degenerate scenario"). 1 Throughout the paper, we use the subscripts, "i", "f.o.", and "f " to represent the initial value, the value at χ 1 freeze-out, and the value well after χ 1 freeze-out, respectively. The subscript "0" is used for present values (except for σ 0 below), while the superscript denotes dark sector quantities that are measured in units of T .

Hierarchical scenario:
In the hierarchical scenario, we assume that DM χ 1 is significantly heavier than χ 2 , and reaches its final abundance via non-relativistic freeze-out through χ 1 χ 1 → χ 2 χ 2 annihilation while χ 2 remains relativistic. For the parameters of our interest, this generally requires m χ 2 m χ 1 /20. We introduce a time-dependent function, ξ, to characterize the temperature ratio between the dark sector (T ) and the SM (T ), with its initial value smaller than unity, ξ i 1. Given that the number of particle species in the SM is much more than in the dark sector we consider, this ensures that the Universe's expansion is solely driven by the visible sector particle content during the radiation-dominated epoch.
While the specific dark freeze-out scenario we consider was not studied in the literature, the analysis of the Boltzmann equation evolution is analogous to some of the existing work on freeze-out in a decoupled dark sector which are useful references for our case, e.g., [31][32][33][34]. In particular, here we use the approach in [33] to calculate the final abundance in terms of the yield variable, Y j ≡ n j /s, where n j is the number density of species i and s = (2π 2 /45)g * S T 3 is the total (dark and visible sector) entropy density, which, by assumption, essentially coincides with the SM one; g * S (g * ) are the effective degrees of freedom in entropy (energy) and we take g * = g * S before neutrino decoupling.
The thermally averaged annihilation cross section σ χ 1 v is conventionally parametrized as where x ( ) = m χ 1 /T ( ) and σ 0 is a reduced "n-wave" annihilation cross section of χ 1 with n = 0, 1, 2, . . . corresponding to s, p, d, . . . -wave annihilation with a respective relative velocity scaling as v 2n in (σ χ 1 v) before taking the non-relativistic thermal average. Concretely, we obtain being the Hubble expansion rate evaluated at a photon temperature T = m χ 1 ; m Pl = 1.2 × 10 19 GeV is the Planck mass. The freeze-out point, x f.o. , measured in terms of photon temperature, is given by Here, c = 0.3 yields good agreement with the exact numerical solution for ξ ≥ 0.01. It is worth mentioning that although the temperature ratio has been treated as a constant function above, one should in practice take its value at freeze-out (referred to as ξ f.o. ), which is slightly larger than an initially attained temperature ratio ξ i : χ 1 annihilation heats up the dark sector, and increases this ratio by a factor 2 1/3 as can be obtained from the conservation of entropy in the self-thermalized dark sector. In practice we adopt a step function for the temperature ratio, ξ, w.r.t. x as: ξ(x ) = ξ i for x ≤ 4 and 2 1/3 ξ i for x ≥ 4. Non-relativistic freeze-out mostly happens at x ≥ 4, where most of the entropy density in dark sector is stored in χ 2 [42].
, one may re-write the solution to the DM abundance as (2.4) in analogy to the standard thermal freeze-out solution: Y χ 1 ,f = (n + 1) x n+1 f.o. /λ. In addition, since the freeze-out point x f.o. , measured in units of the dark temperature T and broadly in the range 3-30, is not sensitive to mild changes of parameters, an s-wave cross section σ χ 1 v /(2 1/3 ξ i ) ∼ O(1) pb is needed to achieve the observed DM abundance in an otherwise standard cosmological thermal history. In case of later co-moving entropy increase by a factor ∆ (to be discussed later), the annihilation cross section should be comparatively smaller, satisfying σ χ 1 v ∆/(2 1/3 ξ i ) ∼ O(1) pb, in order to generate the observed DM relic density. Compared to the standard case T = T , dark freeze-out with (2 1/3 ξ i )/∆ < 1 needs to happen at a relatively higher SM temperature T in order to realize the correct relic abundance.
For example, considering the interaction term ig jχj γ 5 χ j A (j = 1, 2) induced by a heavy pseudoscalar state A (this is one of the interactions realized in our UV complete model presented in Sec. 4), the annihilation cross section is dominantly s-wave and given by where the limit relates to non-relativistic relative initial motion with m A 2m χ 1 . Therefore, in order to generate the observed DM abundance, this points to Since χ 2 particles remain relativistic during the freeze-out (in the case of m χ 2 /m χ 1 x f.o. ), their abundance is fixed in terms of the temperature ratio as Y χ 2 ,f = 0.42g χ ξ 3 i /g * − Y χ 1 ,f . This is because in our setup, by assumption, all interactions in the dark sector preserve the total number of χ 1 and χ 2 , such that Y χ 1 + Y χ 2 = const., up to changes in co-moving entropy.
In the end, the observed baryon asymmetry Y B produced from χ 2 decays is proportional to the "would-be" χ 2 abundance Y χ 2 ,f if it had not decayed: Y B ∝ CP Y χ 2 ,f , where CP is the CP asymmetry of the decay and τ is the χ 2 -lifetime at rest. A later entropy increase, which may be induced by the late decay of χ 2 , would dilute the χ 2 abundance, requiring a larger CP to yield the observed baryon asymmetry today. More explanation about the final yield of matter abundances will be given in Sec. 4.

Nearly degenerate scenario:
In the second, nearly degenerate scenario, we assume that χ 1 is marginally heavier than χ 2 and χ 1 is still stabilized by a Z 2 symmetry. The mass degeneracy is characterised by a dimensionless parameter, δ ≡ (m χ 1 − m χ 2 )/m χ 1 1. The final abundances of DM candidate χ 1 and of meta-stable χ 2 are dominantly set by non-relativistic freeze-out of χ 1 χ 1 ↔ χ 2 χ 2 . The crucial difference from the previous scenario is that χ 2 is now nonrelativistic during the freeze-out of χ 1 . This, in turn, leads to a faster decrease of dark temperature. While dark freeze-out scenarios involving two near-degenerate states have been discussed in literature (see e.g., [35,[43][44][45][46][47]), in our scenario the decay of the metastable state occurs well after the freeze-out era, resulting in different dynamics. Here we provide an analytical solution to this problem, which to our knowledge is novel.
Due to the mass degeneracy, the difference from the earlier discussed hierarchical scenario is that, both χ 1 and χ 2 are non-relativistic at DM freeze-out. At the same time, the total co-moving number of χ 1 and χ 2 particles is still conserved, and the ratio of n χ 1 and n χ 2 in equilibrium depends on their mass splitting and the dark temperature. In this scenario, the temperature ratio T /T can deviate significantly from its initial value during Universe expansion, so we must follow both number and energy density in the dark sector. They are described by the respective Boltzmann equations, as well asρ where the dot signifies the derivative with respect to cosmic time and ρ (p ) is the energy density (pressure) of the dark sector. For a dark sector in kinetic equilibrium, as assumed here, we are to solve the set of three Boltzmann equations to obtain the dark temperature, T , and two number densities, n χ 1,2 (or equivalently, two chemical potentials µ χ 1,2 ). Note, however, that σ χ 1 v and σ χ 2 v are related through the principle of detailed balance; see Eq. (A.2) in the appendix. For illustration, we again take the pseudoscalar mediated model with unaltered annihilation cross section given by Eq. (2.5), and solve the set of equations above for several benchmark parameters. The results are given in Fig. 2. The right panel shows that freezeout is associated with sizable x since the Boltzmann suppression of Y χ 1 with respect to Y χ 2 is determined by the mass splitting (orange dot-dashed lines), 2 instead of the DM mass (red dotted lines), when both DM annihilation and its inverse process are in equilibrium. This is dictated through overall dark number conservation, and by the fact that the inverse process, χ 2 χ 2 → χ 1 χ 1 , only becomes suppressed once T δm. Nevertheless, when the same solutions are shown as a function of x (left panel of Fig. 2) one observes that freezeout happens with x f.o. ∼ 3, smaller than the typical value, x f.o. ∼ 20, in the standard case (ξ ≡ 1) of TeV DM freeze-out, to get the observed DM abundance; the reason for it is explained below Eq. (2.4).
While these equations can only be solved numerically, which is how we obtained our final results, here we also provide analytical explanation/approximation to better demonstrate the underlying physics. An empirical relation between T and T can be obtained using the observation that a decoupled non-relativistic species cools down adiabatically with the scaling T ∝ T 2 ; for a relativistic species T ∝ T instead (away from epochs of entropy transfer). The transition between these two scaling laws happens approximately at T ∼ m χ 2 /4. Using the insights we gained from our numerical study, a good approximation which is an approximation of the exact quasi-static equilibrium solution, Y QSE χ 1 (purple dashed lines); see Eq. (A.4) for more details. Both show that the suppression of Y χ 1 is solely due to the mass splitting, independent of the absolute mass of χ 1 . The observed DM abundance is shown as gray dot-dashed lines.
at any temperature T , where ξ i is the initial temperature ratio, 2 1/3 is the reheating factor discussed earlier. That is, we obtain the relation x = βx 2 at T m χ 2 , with the dimensionless constant β ≡ (1 − δ)/(4 × 2 2/3 ξ 2 i ). See Fig. 12 in the Appendix for a comparison with full numerical results.
Adopting Eq. (2.9) allows us to solve Eq. (2.7) analytically as detailed in App. A. Following the notation of Eq. (2.2), we obtain an approximate solution of a non-relativistic freeze-out, where the freeze-out time is given by (2.11) The solution is different from Eq. (2.2) of the previous scenario, mostly due to T ∝ T 2 at freeze-out here. Nevertheless, in the s-wave limit, we have   Figure 3: Contours of x f.o. (left panel) and Y χ 1 , f (right panel) obtained from the analytical solution of Sec. 2.2 as a function of mass splitting, δ, and the initial temperature ratio between the dark and SM sectors, ξ i . In both panels, we choose m χ 1 = 1 TeV, as well as The dashed black lines show the observed DM abundance. The red (gray) stars depict the parameter choices used in Fig. 2 (Fig. 13), where the Boltzmann equations are solved numerically for δ = 0.1 and ξ i = 0.1.
in both scenarios. It in turn suggests that Ω DM is generally not sensitive to the DM mass, as it cancels in the ratio x f.o. /m χ 1 and for as long as the s-wave annihilation cross section σ 0 is held constant. Hence, the yield decreases for increasing σ 0 .
The analytical solutions are illustrated in Fig. 3 by scanning the parameter space for m χ 1 = 1 TeV and g 1 g 2 /m 2 A = 1/(4 TeV) 2 . For ξ = 0.1 and δ = 0.1, it gives the freeze-out point x f.o. ≈ 3, as well as the observed DM abundance, in good agreement with the exact numerical results of Fig. 2. For additional numerical results of n χ 1,2 evolution with various parameter sets, see Fig. 13 in the appendix. For instance, reducing the mass splitting by a factor of 10 to δ = 0.01 gives x f.o. ≈ 8 and a final abundance about ten times larger (middle gray star in Fig. 3). This can be explained as follows: estimating the freeze-out condition That is, smaller mass splitting leads to a lower freeze-out temperature, and thus larger value of the final χ 1 yield. In contrast, reducing the temperature ratio gives an less populated dark sector, and thus a smaller final abundance of χ 1 .
Finally, we turn to the χ 2 abundance. Number conservation of the sum of χ 1 and χ 2 directly translates into Y χ 2 ,f ≈ 0.42ξ 3 i g χ /g * − Y χ 1 ,f at the end of freeze-out. Eventually, we demand χ 2 to decay into the visible sector, in analogy to the hierarchical scenario before. 3. The left panel shows the evolving ratio of matter-to-radiation energy densities, ρ M /ρ R (dashed line), and ∆-the ratio of entropy density defined by a 3 (ρ M + ρ R ) 3/4 to its initial value (solid line); see Eq. (2.14). In the right panel, the purple solid (black dashed) line shows the evolution of the SM temperature T in our model (in standard cosmology). In both panels, the vertical gray lines depict the beginning (end) of this epoch at a α (a β ), while horizontal gray lines correspond to unity in left panel, and T = m χ 2 in right panel.

Potential matter-domination caused by χ 2
In both the hierarchical and near-degenerate scenarios described above, at the end of DM freeze-out the χ 2 abundance is given by As we focus on ξ i ≥ 10 −2 and the DM abundance should satisfy Y χ 1 ,f ≈ 10 −10 (4 GeV/m φ ) < 10 −10 for m χ 1 well above GeV, the equation above gives Y χ 2 ,f Y χ 1 ,f , thus the second term on the R.H.S. can be neglected in practice. Note that the insensitivity of Y χ 2 ,f to the annihilation cross section is a result of the dark number conservation and the Boltzmann suppression of Y χ 1 ,f relative to Y χ 2 ,f around freeze-out time due to the mass splitting. Interestingly, depending on the lifetime of χ 2 , the χ 2 abundance, now mainly decided by ξ i , may lead to an epoch of early matter-domination (EMD). This would modify the prediction for current-day particle abundances, due to entropy injection at the end of EMD.
If χ 2 decays while still being relativistic, the Universe's expansion, by assumption, is dominated by SM radiation, and the final baryon asymmetry remains unaltered from Eq. (2.17). In contrast, the decay of a non-relativistic species transfers additional entropy to the SM sector. The dilution factor, measured in terms of the ratio of total co-moving entropies before (S α ) and after (S β ) decay may be estimated as where the subscripts, α and β, mark the beginning and end points when the Universe exits and re-enters the radiation-domination, respectively. Note that the radiative entropy (S) does not include the non-relativistic matter contribution; the last quantity is thereby used during the EMD epoch. Using the energy conservation lawρ = −3H(ρ + p) and the equation of state p = w ρ, one obtains where ω = 1/3 (0) for a radiation-dominated (matter-dominated) Universe. During the transition, ω varies and needs to be calculated numerically. The results for one benchmark parameter set is shown in Fig. 4. We may also obtain a simplified analytic estimate for the dilution factor ∆ by assuming a sudden transition between matter-and radiation-dominated epochs. For this, one may take ρ ∝ a −3 during matter domination, thus ∆ ≈ (a β /a α ) 3/4 from Eq. (2.14), i.e., ∆ is mostly decided by the duration of this period. Matter-domination starts at the scale factor a = a α , when ( The matter-domination epoch ends when the Hubble rate equals the χ 2 decay rate, and ρ M (a β ) = 3Γ 2 Evidently, when ∆ ≈ 1, the matter-domination epoch was either very brief or never happened. The expression above can be tested against the numerical solution of Fig. 4. 3 For example, while Eq. (2.16) yields ∆ ≈ 41, Fig. 4 points to ∆ ≈ 43. For better guidance, we point out that ∆ 10 3 for m χ 2 smaller than 1 TeV, due to g χ 2 g * and the Big Bang Nucleosynthesis (BBN) bound on the lifetime of χ 2 .
In our setup, the EMD epoch occurs after DM freeze-out. Therefore, while the dilution reduces the values of Y χ 1 and Y χ 2 separately by a common factor ∆, it does not change their ratio. That is, in absence of other baryon asymmetry-depleting processes we have, 17) and the predicted relation between the DM abundance and baryon asymmetry is not affected.
3 χ 2 Decays, Baryogenesis, and the DM-Baryon Coincidence As outlined earlier, the χ 2 abundance following the dark freeze-out is meta-stable, and its CP-and B-violating, out-of-equilibrium decays trigger baryogenesis. We require that χ 2 decays after freeze-out of both χ 1 and χ 2 , but before primordial nucleosynthesis: . This way, freeze-out and baryogenesis can be treated separately, and the wash-out effect from inverse decay of χ 2 is ineffective during the baryogenesis stage. 3 Note that the decoupling of heavy SM particles contributes to the evolution of T in the right panel of Washout processes such as udd →ūdd are suppressed by the large mediator mass (required by collider constraints, see Secs. 4 and 5) in the effective operator for the process. A more detailed discussion on suppressing various washout processes can be found in [22] which is in analogy to the scenario considered here. The initial condition for Y χ 2 of this stage of evolution (i.e., χ 2 decay) is set by the wouldbe abundance of χ 2 after the freeze-out of χ 1 , which is an essential factor for predicting the baryon asymmetry. Recall from Sec. 2.3, via the overall number conservation within the dark sector, the χ 2 abundance after the freeze-out of i . Assuming that it violates baryon number by one unit (which is realized in the model shown in the next section), the comoving density of baryon asymmetry where CP is the CP asymmetry in χ 2 decays, and Γ W is the rate of washout processes. The ellipses stand for any possible additional sources of the baryon asymmetry prior to χ 2 -decay, 4 which we simply assume to be zero in the following. As mentioned earlier, the rate of washout processes that can potentially reduce the baryon asymmetry is weak by construction and the exponential factors can be dropped. This gives a simple approximate solution for the co-moving baryon asymmetry in either, hierarchical or near-degenerate, scenario, where Y χ 1 ,f is given in Eqs. (2.4) and (2.10). The result in Eq. (3.2) directly connects the co-moving baryon asymmetry to the relic abundance of χ 1 . Assuming that χ 1 composes all of DM observed today, we obtain the ratio of the DM and baryon abundances observed today as where m p is the proton mass; the degrees of freedom are to be evaluated at the time of DM freeze-out g * ≈ g * S ∼ 100. 5 In the case of s-wave annihilation (n = 0), neglecting the second term allows to further simplify the expression to to the baryon asymmetry can be accounted for by the addition of a term Y initial where the exponential suppression is due to "wash-out." 5 In the numerical evaluation we carry along the temperature dependence of g * and g * S, whereas in the analytical estimates we take them as constants.
The result in Eq. (3.3) for Ω B takes a similar form as the conventional WIMP miracle, augmented by the CP-asymmetry factor CP and proton-to-χ 1 mass ratio m p /m χ 1 . Taking similar values for the annihilation cross section, splitting δ, and dark sector temperature ξ i to those in Fig. 2 and Fig. 13, the observed coincidence between dark and baryonic matter abundances can be achieved with m χ 1 ∼ O(10 − 10 4 ) GeV and CP ≈ O(10 −8 − 10 −2 ). In the next section, we implement the general dynamics of χ 1 -χ 2 freeze-out and χ 2 decay-triggered baryogenesis in a specific model which achieves exactly this with dedicated numerical analysis.

Exemplary UV-complete Model
In this section we present a concrete model with fermionic dark states, which include interactions realizing the χ 1 χ 1 → χ 2 χ 2 freeze-out and generation of the baryon asymmetry. For the latter, the three Sakharov conditions [48] are met by the CP-and baryon numberviolating out-of-equilibrium decay of χ 2 to SM quarks.
We start by considering dark sector interactions necessary for freeze-out. There are a number of options to achieve this and our core mechanism is insensitive to the detailed realization of it. Here, we consider interactions mediated by a pseudoscalar, A, or a scalar, S. The associated Lagrangian then reads, where the SM singlet Majorana fermions χ j (j = 1, 2, 3) are the members of the dark sector. The interactions mediated by A allow for s-wave annihilation, whereas the leading contribution in a velocity expansion of the annihilation cross section mediated by S is p-wave. In order to avoid additional complications of a dynamical role that χ 3 may play during dark freeze-out, we may take g ( ) 3 1. This allows us to neglect χ 1,2 χ 1,2 ↔ χ 3 χ 3 processes and χ 3 -induced energy transfer between dark and SM sectors even for β = O(1) (β is to be introduced in Eq. 4.3). In addition, for m A(S) ≥ 2m χ 1 we may also neglect the processes χ 1,2 χ 1,2 ↔ AA(SS). Taken together, we then map onto the analysis of dark sector freeze-out of the preceding section. In passing, we note that a hierarchical coupling structure among g ( ) j can be addressed by an underlying dark flavor physics model, but which is beyond the scope of this work.
In order to predict the yields of DM and baryon asymmetry from this model, we first analyze the freeze-out process due to the annihilation χ 1 χ 1 ↔ χ 2 χ 2 mediated by either the pseudoscalar A (s-wave) or the scalar S (p-wave). The non-relativistic thermally averaged annihilation cross-section (away from the resonances) reads where there is no interference between the two mediators, and the proportionality to √ δ indicates the available limited final state phase space. Apparently, in the case of g j = g j and m S = m A , the s-wave component mediated by A dominates the DM annihilation in the non-relativistic regime. Figure 5: Tree-level and one-loop diagrams contributing to CP violation in χ 2 decays.
We now discuss the part that enables baryogenesis. The associated Lagrangian is given by where the chiral projectors ensure that only the right-handed quarks u j and d j of generation j participate, kl is an antisymmetric symbol for flavor indices k, l. As previewed in general discussions, χ 1 is charged under a Z 2 symmetry. Hence, it does not couple to quarks and plays the role of the DM. χ 2 is the meta-stable state whose decays trigger baryogenesis. We include a third member in the dark sector, χ 3 , which enables the CP asymmetry-inducing interference between tree and loop-level diagrams. A complex scalar, or "di-quark" φ is also introduced. This state transforms as an anti-triplet under SU (3) C , has the same SM charge as the right-handed up-type quarks, and couples to the right-handed down-type quarks of the SM. All Yukawa couplings α j , β j , η jk are generic, complex numbers. Because of stringent LHC constraints on color-charged particles φ and to allow for a tractable exposition, we shall consider the hierarchy m φ m χ j . In addition, we take m χ 2 > m χ 3 so that contributions to the baryon asymmetry solely arises from χ 2 decays as shown in Fig. 5. The analogous diagram for χ 3 decay with χ 2 in the loop does not contribute to an analogous CP asymmetry from χ 3 , since here the kinematic cut is forbidden as m χ 2 > m χ 3 . 6 Finally, an ordering |α| |β| < |η| will be assumed to ensure the branching ratio Br φ→d k d l ≈ 1, and a sufficiently early decoupling of dark particles from the SM. For instance and as shown below, the smallness of |α| together with a heavy φ guarantees that χ 1,2 decouple from both SM particles and χ 3 before DM freeze-out. The hierarchy |α| |β| further forces χ 3 to have a lifetime shorter than χ 2 , decaying into SM quarks well before BBN.
The Lagrangian (4.3) enables χ 2 decay that triggers baryogenesis. With the assumed mass hierarchy m φ m χ j , we may integrate out the heavier di-quark scalar, yielding the 6 If mχ 2 < mχ 3 , it would be χ3, instead of χ2, whose decay contributes to the CP-asymmetry by switching the subscripts 2 ↔ 3. Alternatively, if mχ j m φ , the kinematic cut can be made through φ and up-quark propagators, where both χ2 and χ3 decays would have non-vanishing CP , and, in principle, contribute to the baryon asymmetry.

effective operators
They induce the three-body decay χ 2 → udd with |∆B| = 1. The associated decay rate at T < m χ 2 is given by where the factor of three in Eq. (4.5) accounts for color multiplicity. The CP-asymmetry produced in decays arises through the interference between tree-level and 1-loop diagrams shown in Fig. 5. The asymmetry as a function of m χ 2 , m φ , and Yukawa couplings is given by [25] in the process of χ 2 → u j d k d l (andū jdkdl ) with an intermediate up-type quark u m . Since we are interested in a broad range of m χ 1 spanning over O(10 − 10 4 ) GeV, all generations may contribute when m χ 2 m t . Nevertheless, the constraints from (di-)nucleon decay and neutron-antineutron oscillation on the effective operators in Eq. (4.4) are very strong if (u j d k d l ) that couples to χ 3 contains only light quarks, such as (uds); see e.g., [49]. Thus, here we focus on the combination of (cds), or (udb), where sizable coupling-combination βη is allowed, and χ 2,3 decays into them are kinematically accessible. Finally, we note that the model is exempt from constraints from the null measurements of the neutron electric dipole moment (EDM). The reason is that the interference diagrams in Fig. 5 leading to CP violation only involve right-handed quarks [50,51].
Having established all the ingredients and hierarchies of a UV complete model, we now turn to the identification of prospective parameter ranges. As discussed above, DM decoupling with ξ i ≤ 1 prefers an earlier freeze-out in comparison with the standard case. Under the assumption that χ 1 freezes out for T m χ 2 /10, the requirement that χ 2 decays after χ 1 freeze-out but before BBN yields the preferred range of allowed decay couplings as (4.7) Note that the last inequality is also approximately the condition that χ 2 decouples from the SM bath during χ 1 freeze-out: if the bound were violated, the process χ 2 ↔ u j d k d l becomes efficient and dark and SM sectors thermalize. As a result, for m φ = 10 TeV, the hierarchy |α| |η| introduced for baryogenesis leads to |α| 10 −2 for m χ 2 = 1 TeV, and |α| 0.14 for m χ 2 = 30 GeV. Similarly, requiring χ 3 to decay fast, say at T ∼ m χ 3 , suggests |β||η| 0.1(m φ /10 TeV) 4 (10 GeV/m χ 3 ) 3 , which can be easily satisfied in the parameter region of our interest here. So the existence of χ 3 neither leads to additional early matter domination, nor affects the standard BBN processes. Putting all together and plugging the specific results into Eq. (3.3), we obtain the relic abundance of DM and the baryon asymmetry. Taking Eq. (3.4), their ratio in the exemplary UV model for s-wave dominated freeze-out becomes 8) for which, and for numerical results below, we have assumed the Yukawa couplings between χ 2 and the up-type quarks to be real such that the CP-asymmetry in Eq. (4.6) can be written in terms of a complex phase θ β of β : Im[(ααβ * β * ) 2 ] → |α| 2 |β| 2 sin(2θ β ). We observe, that the final ratio of Ω B to Ω χ 1 is sensitive to the mass splitting δ and m A /m χ 1 , instead of m χ 1 alone. In Fig. 6 we show contours of Ω B /Ω χ 1 = 0.2 as a function of the χ 1 mass and DM annihilation coupling g ≡ √ g 1 g 2 in the hierarchical (left) and the nearly degenerate (right) cases of our exemplary model. Note that the contributions from the scalar S is subleading and is thus not taken into account. To obtain the numerical results, we plug Eqs. (4.2) and (4.6) into the solution for the χ 1 freeze-out abundance given by Eqs. (2.4) and (2.10), which is then used together with Y χ 2 ,f ≈ 0.42ξ 3 i g χ /g * − Y χ 1 ,f to find the abundance of χ 2 just after DM freeze-out. This, in turn, gives rise to the observed baryon asymmetry, Y B ∝ CP Y χ 2 ,f . The contours have fixed mass splitting in the range of 0.97 < δ < 0.997 (or equivalently, 0.003m χ 1 < m χ 2 < 0.03m χ 1 ) in the hierarchical case and 0.05 < δ < 0.5 in the degenerate case. The pseudoscalar mass is taken to be m A = 2.1 TeV in Fig. 6. In the limit of m A m S and/or g g, the annihilation channel via the scalar S dominates the DM freeze-out, which case is shown in Fig. 7. In our analysis, we fix θ β = π/4 to obtain the maximal CP-asymmetry.

Experimental Signatures
The concrete model outlined in the previous section lends itself to observational signatures and experimental tests, in particular in relation to the stable DM candidate χ 1 and its unstable dark sector partners χ 2,3 . These include conventional and novel direct and indirect detection signals, as well as modifications to primordial density fluctuations. We will sketch the prospects of these searches below, leaving a more-in-detail analysis for future work.

New particle searches at colliders
To produce the baryon asymmetry with renormalizable operators, the color anti-triplet, di-quark scalar φ has been introduced. As a result of the color charge it carries, LHC bounds on its mass are quite strong. As outlined in Sec. 4, the di-quark scalar φ decays to two down-type quarks. At the LHC, φφ * are pair-produced dominantly via gluon fusion, and subsequently decay to 4 jets: pp → φφ * → 4j, as shown in Fig. 8. LHC searches for diquark scalars constrain the mass of φ. For instance, the recent CMS search [52] constrains the production cross section of di-quark scalar resonances to be below fb at √ s = 13 TeV, excluding m φ 7 TeV. The χ j (j = 1, 2, 3) particles can also potentially be produced at colliders. The production of χ 1 is shown on the left panel of Fig. 9, and would lead to a missing transverse energy (MET) signal. However, due to the combination of loop factor, the required large mass of the di-quark scalar φ in the loop diagram, and the additional initial state radiation (ISR) photon or jet required for tagging the event, the χ 1 production rate is found to be too small to be detectable. The production of χ 2 and χ 3 can proceed via t-channel φ-exchange as shown in the right panel of Fig. 9. The subsequent χ 2,3 decays can lead to prompt jets, displaced vertices or MET, depending on their decay length. Again, m φ 7 TeV suppresses the production rates of χ 2 and χ 3 . Since successful co-genesis requires |α| 10 −2 for m χ 2 = 1 TeV, the production cross section of χ 2 falls significantly below fb, well beyond the LHC reach. The prospect for χ 3 production is better, as the coupling |β| can be O(1), leading to fb-sized LHC production cross sections for √ s = 13 TeV and m φ = 10 TeV [53].
The prospect of detecting χ 3 depends on its decay length following production. It may vary over a wide range, where γ = E χ 3 /m χ 3 is the Lorentz factor; the same expression also applies to χ 2 after replacing m χ 3 and β with m χ 2 and α. χ 3 decay can hence lead to prompt multijet events, MET, or displaced jet signals. For the size of the production cross section, the displaced vertex channel is the most promising one, with the potential of detection at the High-Luminosity LHC [28], due to the generally low background [54][55][56][57].
We would like to note that although with the LHC it is generally challenging to detect χ j particles in the model as outlined above, the detection prospects are expected to improve notably with the proposed future higher energy colliders [58][59][60][61][62]. Furthermore, model-variations can increase the production rate of χ j . For instance, interaction terms, such as a trilinear A-A-H, may open up new channels at higher energy colliders for both χ 1 and χ 2 production with potentially appreciable rate, independently of m φ . In addition, while we consider φ dominantly decays into a pair of down-type quarks, a sizable branching ratio of φ → χ 3 u is possible, depending on the ratio |β/η| 2 , which provides a new, potentially efficient channel for χ 3 production as shown in the right panel of Fig. 8; for χ 2 the corresponding ratio is much smaller and less prospective.
In summary, the di-quark and χ j 's in the co-genesis framework that utilizes quarkcouplings in its UV-complete representation could leave observable signals in various search channels at high luminosity run of the LHC, yet are potentially challenging. On the other hand, the proposed future colliders at the high energy frontier have more promising capacity to reveal detectable signatures from this framework. Figure 10: Diagrams contributing to DM-nucleus scattering χ 1 N → χ 1 N . The diagram on the left is the dominant contribution in the minimal model while the diagram on the right may dominate if the S-Higgs mixing is present.

Direct detection of χ
In the example model detailed in Sec. 4, there are no interactions between the stable DM candidate χ 1 and quarks at tree-level. However, DM-nucleus interactions can be induced at 1-loop as shown in the left panel of Fig. 10. Integrating out the heavier mediating states generates low-energy effective operators which induce interactions between χ 1 and nucleons N = n, p. They are either spin-independent, L ∝N N S, or spin-dependent, L ∝N (iγ 5 )N A. The effective interactions lead to DM-nucleon elastic scattering cross sections which can be estimated using dimensional analysis: which applies in the limit m φ,S m χ 2 m N [63]. Here the first component, mediated by S, is a spin-independent cross section, while the latter, mediated by A, is both velocity suppressed and spin-dependent; the coefficient f N ≈ 0.3, is counting for the valence quark content of the nucleon [64]. For instance, under the assumption that g 1 g 2 /m 2 A = g 1 g 2 /m 2 S = 1/(4 TeV 2 ), m χ 1 = 3 TeV, and m χ 2 = 0.1 TeV, which gives the observed DM abundance with ξ ∼ 1, setting m φ = 10 TeV together with Eq. (4.7) require |α| 1, resulting in the spin-independent part σ SI χ 1 N 10 −56 cm 2 ; the spin-dependent component is even smaller due to the velocity suppression. We conclude that in the current setup, there are no direct detection prospects. Finally, we point out that there can be contributions to direct detection at tree-level when the scalar S mixes with the SM Higgs. In order for the dark sector and SM to remain decoupled as per our working hypotheses, the S-Higgs mixing angle must be small, θ S 10 −7 . In the case of m S ≥ m h , the scattering cross section between χ 1 and proton through mixing can be estimated as 10 −42 (g 1 θ S ) 2 cm 2 [65], and thus will not be detectable in the foreseeable future.
An interesting possibility arises through t-channel scattering of protons (or neutrons) with χ j at one loop via off-shell A (S), χ 2 , and φ. An example of this process with a proton initial state is shown in Fig. 11. Therefore, it is important to note that the relative kinetic energy of the galactic χ 1 DM flux at earth interacting with nuclei falls short to kinematically access this channel, but rather requires a "boosted" χ 1 component. A particularly interesting environment may involve neutron stars in regions of high χ 1 density accelerating χ 1 to sufficiently large energy [66], such that this nucleon conversion process occurs and leads to observable effects such as anomalous heating. Note that the energy spectrum of the incoming and outgoing states can differ significantly from the elastic scattering usually considered [67,68]. We leave a detailed investigation of nucleon conversion via capture in neutron stars, and other stellar objects, to future work.
In addition, DM can be captured by celestial objects, and consequently annihilate inside. Given the short lifetime of χ 2 , one particularly interesting signature is that the boosted χ 2 produced from DM annihilation in the centre of Earth decays inside a largevolume terrestrial detector, for which we take IceCube [76] as an example. The signal requires the in-flight decay length of χ 2 , l χ 2 ≡ τ c m 2 χ 1 /m 2 χ 2 − 1, to be on the order of or longer than the Earth's radius, r Earth . If satisfied, the probability for χ 2 to decay inside its detector of fiducial volume 1 km 3 is given by Assuming a sensitivity threshold of several events per year in the detector same as [77] and l χ 2 /c = 0.02 sec, the probability given by Eq. (5.3) implies that the DM annihilation rate, Γ A , in the Earth's centre should reach 10 6 particles per second, in order to be detectable in IceCube, which is in agreement with [78]. We may consider the prospects of such signature in the concrete UV model. As discussed in Sec. 5.2, the non-relativistic scattering between χ 1 and nucleons is dominated by the S-mediated interaction, leading to a spin-independent cross section σ SI χ 1 N . In the parameter region of interest, where the annihilation cross section is σ χ 1 v ≤ 1 pb, equilibrium between capture and annihilation in Earth is not reached [76,79]. That is, the actual DM annihilation rate in Earth scales as (σ SI χ 1 N ) 2 σ χ 1 v , up to subleading effects [76]. Moreover, a comparison of Figs. 6 and 8 in [76] allows us to estimate the coefficient of this scaling for m χ 1 ≥ 100 GeV, leading to Γ A ∼ 10 9 s −1 10 TeV As an example, for 100 GeV DM with σ χ 1 v = 0.1 pb, detecting O(1) events per year requires Γ A ∼ 10 6 s −1 , implying a sensitivity of σ SI χ 1 N ∼ 10 −45 − 10 −46 cm 2 . This is comparable to current direct detection constraints, but still needs to be improved by several orders of magnitude to eventually probe our model. Such improvement may be achieved if the DM annihilation is significantly enhanced at low velocities (e.g., through Sommerfeld enhancement with a light mediator), bringing the Earth capture and annihilation into equilibrium. Note that this simple re-scaling does not apply for m χ 1 ≤ 100 GeV, where the bound gradually gets weaker, except for DM masses that trigger resonant capture of abundant elements in Earth [76]. A more detailed investigation of captured DM abundances is deferred to future work.

Modifications to primordial density fluctuations
New particles and their associated dynamics in our framework may also leave an imprint on cosmological observables. One aspect is structure formation. For instance, the metastable χ 2 can induce an EMD epoch, during which the subhorizon density perturbations grow linearly. In contrast, the motion of DM particles leads to damping of density perturbations at small scales. As studied in the literature (e.g., [80] and most recently in [81]), the total effect on density perturbations is mainly decided by two scales: the horizon size at the end of the EMD epoch and the scale below which primordial fluctuations are suppressed by DM streaming out of over-and under-dense regions.
The scale associated with perturbation growth is decided by the horizon size at cosmic time t = τ , where χ 2 decays and the EMD epoch ends. Its value can be expressed in terms of co-moving distance as where s(T ) is the radiation entropy density at temperature T , T end is the reheating temperature post the EMD epoch; the photon temperature at present is T 0 = 2.35 × 10 −4 eV [1]. Perturbations at scales below L end thus grow linearly w.r.t. the scale factor a during this EMD epoch [80,82]. On the other hand, density perturbations of DM can be suppressed due to random energy transfers either via sound waves or traveling particles, referred to as acoustic or collision(-less) damping [83][84][85]. Since χ 1 does not significantly scatter after its freeze-out, collisionless damping dominates in the suppression of structure. This free-streaming scale is estimated by where we have neglected the logarithmic dependence on the exact time when DM kinetically decouples. Finally, we note that the last equality of Eq. (5.6) assumes a standard cosmology.
In presence of a matter-dominated epoch, the ratio of a 0 /a(t) becomes larger than in standard cosmology at t < τ (compare solid and dashed lines in right panel of Fig. 4), thus a "stretching factor" ∆ 1/3 needs to be multiplied in this case.
Although there exist two competing effects as outlined above (enhancement and suppression), observationally, the most relevant scale is the larger one of L end and L fs . In our framework, the scale governing additional perturbation growth, L end , is larger than the suppression scale, L fs , as the former is the horizon size when DM is already non-relativistic. As a result, we generally expect a peak in density perturbations between the two scales, leading to a potentially observable enhancement in the local subhalo abundance. The ensuing altered predictions of the ionization history as well as of the halo-mass function in the late Universe can be probed with various strategies [81,[86][87][88][89][90]. In particular, future Pulsar Timing Arrays (PTA) will be able to constrain the abundance of local subhalos down to a halo mass of 10 −10 M , and thus have the sensitivity to probe effects on primordial perturbations at co-moving scales as small as 10 −7 Mpc [89,91]. Such a scale corresponds to τ ∼ 10 −6 s, or, conversely, a reheating temperature of T end ∼ 300 MeV following the end of the EMD era.

Conclusions
In this work, we identify a new mechanism for the joint generation of the baryon and DM abundances. DM χ 1 makes a thermal freeze-out via annihilation into a lighter metastable dark partner χ 2 through an overall dark number conserving process χ 1 χ 1 → χ 2 χ 2 . The lighter state χ 2 subsequently decays to SM quarks. Its chemical decoupling and CP-and B-violating interactions with SM ensure the fulfillment of the Sakharov conditions, while its interactions with χ 1 ensure fulfillment of the relic density requirement by observations. The lifetime of χ 2 is assumed such that its decay happens after χ 1 -χ 2 freeze-out yet before primordial nucleosynthesis. By itself, this only requires small couplings to SM, and, in general the dark and observable sector temperatures, T and T , may differ.
We present a novel analytical treatment for a two-state dark sector freeze-out and subsequent baryogenesis assuming T ≤ T and for which T = T is contained as a special case. There are then two principal options. First, in the hierarchical scenario there is a mass hierarchy m χ 1 /m χ 2 10 such that χ 2 is relativistic during χ 1 freeze-out. The relic abundance of DM, χ 1 , is established following the dark freeze-out, and depends on the annihilation cross-section and the ratio T /T , resembling the prediction familiar from WIMPs. The χ 2 abundance prior to its decay is also fixed by the freeze-out, yet is insensitive to the annihilation cross-section. Consequently, the ensuing prediction of the baryon asymmetry is predominantly determined by the CP asymmetry CP and the initial value of T /T . In the nearly degenerate scenario (δ 1) both χ 1 and χ 2 freeze out non-relativistically. However, because of overall χ 1 + χ 2 number conservation in the dark sector, its freeze-out yield again essentially coincides with that of radiation such as in the previous case. Therefore, the baryon asymmetry likewise only depends on CP and the temperature ratio.
In either scenario, χ 2 may dominate the energy budget of the Universe before its decay, leading to an early matter dominated era. The associated entropy injection dilutes both baryon asymmetry and DM abundance. However, it leaves their relative proportion unchanged. In the case of s-wave annihilation, the prediction for Ω B /Ω DM only depends on the DM annihilation cross section-which may take on its usual thermal value O(pb)-and on CP , once the initial temperature ratio is fixed; see Eq. (3.4). We verify our analytical estimates by numerically solving the Boltzmann equations, and find excellent agreement for benchmark parameter points.
We then realize our general ideas by introducing a UV complete model. Here, a heavy pseudoscalar A mediates the dominant, s-wave annihilation of the fermionic DM state χ 1 , while a massive CP-even scalar mediator S may play a more favorable role for detection prospects in DM direct detection experiments. The connection to the SM is made through Yukawa interactions between fermionic states χ 2,3 , an SU (3) c charged scalar φ, and SM quarks. These interactions mediate the tree-level B-violating decay of χ 2 which-through its interference with the loop-induced decay by the intermediate state χ 3 -becomes CPviolating. The colored state φ can, e.g., be pair-produced through gluons and decay to four jets, and is currently constrained by the LHC data to be at least of multi-TeV mass. Pair production of χ 1,2,3 would lead to missing energy signals and/or in the form of displaced vertices. Some of these collider signatures can be within reach of the High-Luminosity LHC, while future high energy colliders are more promising for detection. The conventional DM direct detection signal from this model is strongly suppressed due to the decoupled nature of the dark sector. Model-specific indirect signals of baryogenesis from DM capture and annihilation in the Earth and other stellar objects, exotic signatures at large volume neutrino experiments, and departures of a standard matter power spectrum warrant further study.
In summary, we identify a new generic mechanism connecting the DM and baryon abundances. We term it "Dark freeze-out Cogenesis" on account of a separate thermal dark sector evolution in the early Universe, and that both DM abundance and baryon asymmetry are seeded during the same freeze-out dynamics in a dark sector. Our mechanism works over a vast range of dark sector masses. The price to pay is restrictions on the size of the interactions with the SM that ensure certain decoupling conditions to hold. This renders direct and indirect tests more challenging, while these restrictions may be lifted with future model-building efforts in this direction. for χ 1 and χ 2 , the evolution of χ 1 followṡ The principle of detailed balance dictates that the right hand side of Eq. (A.1) vanishes in equilibrium, which gives the thermally averaged cross section for χ 2 annihilation, σ χ 2 v , in terms of σ χ 1 v , where δ ≡ (m χ 1 −m χ 2 )/m χ 1 is the dimensionless mass splitting, g χ 1 ,χ 2 are the number of internal degrees of freedom in χ 1 /χ 2 , and x = m χ 1 /T . The exponential factor arises from the equilibrium number density with a vanishing chemical potential, n eq = g mT 2π 3/2 e −m/T . Plugging Eq. (A.2) into Eq. (A.1) and replacing the number densities n χ 1,2 with the comoving densities Y χ 1,2 and cosmic time t in favor of x = m χ 1 /T = x ξ, the χ 1 abundance, Y χ 1 ≡ n χ 1 (T )/s(T ), evolves as with other quantities explained in the main text. At T m χ 1 , number conservation requires Y χ 1 + Y χ 2 ≡ Y χ 1 ,i + Y χ 2 ,i ≈ 0.42g χ /g * S ξ 3 i , if setting initially the temperature ratio ξ = ξ i at T ≥ m χ 1 . When the dark sector is in thermal equilibrium, the quasi-static condition, dY χ 1 /dx ≈ 0, needs to be satisfied for Eq. (A.3), leading to Y χ 1 ≈ (1 − δ) −3/2 e −δx Y χ 2 = (1 − δ) −3/2 e −δx (Y χ 1 ,i + Y χ 2 ,i − Y χ 1 ).
(A. 4) In general, the denominator is approximately one for cold freeze-out, i.e., if x f.o. 1/δ. We emphasize that our QSE solution is significantly different from that of a WIMP freeze-out. In the WIMP case, DM pair annihilate to final states that remain in chemical equilibrium with the thermal bath so that the associated Boltzmann equation is directly dY WIMP /dx = −λ/x 2+n [Y 2 WIMP − (Y eq WIMP (µ = 0)) 2 ], leading to a QSE solution Y QSE = Y eq WIMP (µ = 0), when all relevant processes are sufficient. In our framework, non-vanishing chemical potentials arise for both dark particles due to the total number conservation of χ 1 and χ 2 , which is why now its QSE solution has a dependence on the initial total dark abundance.
An implicit assumption here is that both particles are non-relativistic when freeze-out happens, i.e., our results below are not valid for δ → 1. The latter case is instead studied in the main text in the hierarchical scenario of Sec. 2.1.
To solve the equation (A.3) at later time when the annihilation process gradually becomes out of equilibrium, we recast it in terms of ∆ neglecting sub-leading terms under the assumption δx 1. Note that ξ = T /T is not a constant due to the fact that for a free non-relativistic particle T ∝ 1/a 2 while T ∝ 1/a for the radiative thermal bath. That is, one can set ξ = 1/(βx), and thus x = βx 2 in the non-relativistic limit T m χ 2 . As argued in the main text, a reasonable choice for the coefficient is β = (1 − δ)/(4 × 2 2/3 ξ 2 i ), as confirmed by the numerical results in Fig. 12. Using this, we may write (A.6) Following standard procedure, we assume the L.H.S. of Eq. (A.5) to vanish at freezeout x f.o. , and take ∆(x f.o. ) = cY QSE (x f.o. ), to solve for x f.o. . As fiducial value, we choose c = 0.3 with a rather mild dependence on the variation of that number. It then follows that at x = x f.o. we may write, Substituting Y χ 1 ,i +Y χ 2 ,i ≈ 0.42g χ ξ 3 i /g * S and ξ = 1/(βx), the exponential above is expressed as where the last term is obtained from Eq. (A.8). We may verify our analytical results above by comparing to the numerical results shown in Fig. 13, where one finds good agreement. As expected, the evolution of T is insensitive to the detailed processes. This is because it is essentially governed by the dominant number abundance Y χ 2 , and, thus, by the lighter state m χ 2 and ξ i .