Saha equilibrium for metastable bound states and dark matter freeze-out

The formation and decay of metastable bound states can significantly decrease the thermal-relic dark matter density, particularly for dark matter masses around and above the TeV scale. Incorporating bound-state effects in the dark matter thermal decoupling requires in principle a set of coupled Boltzmann equations for the bound and unbound species. However, decaying bound states attain and remain in a quasi-steady state. Here we prove in generality that this reduces the coupled system into a single Boltzmann equation of the standard form, with an effective cross-section that describes the interplay among bound-state formation, ionisation, transitions and decays. We derive a closed-form expression for the effective cross-section for an arbitrary number of bound states, and show that bound-to-bound transitions can only increase it. Excited bound levels may thus decrease the dark matter density more significantly than otherwise estimated. Our results generalise the Saha ionisation equilibrium to metastable bound states, potentially with applications beyond the dark matter thermal decoupling.


Introduction
The thermal decoupling of dark matter (DM) from the primordial plasma is among the most widely considered mechanisms for its production. In thermal-relic DM scenarios, and for DM masses around and above the TeV scale, the interactions of DM or its co-annihilating partners typically manifest as long-range due to the hierarchy between the large DM mass and the smaller mass of the force mediators. Multi-TeV scenarios that do not feature such a hierarchy typically fail to provide for sufficient DM annihilation in the early universe and are thus not viable in this context, as suggested by model-independent unitarity arguments [1,2].
Long-range interactions imply the existence of bound states. The formation of meta-stable bound states and their subsequent decays into radiation open new annihilation channels for DM that can significantly reduce its predicted relic density [3]. Besides being supported by unitarity arguments, the emergence of this type of inelasticity in the multi-TeV regime has been shown by explicit calculations in a variety of theories [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. In addition to DM freeze-out, its implications extend to DM indirect detection and collider signatures, as discussed in the original works [22][23][24] and more recently in [25][26][27][28][29][30][31]. 1 If metastable bound states exist in the spectrum, then DM freeze-out in the early universe is governed by a system of coupled Boltzmann equations for the unbound species and the bound states [3]. Besides the standard (co-)annihilations of the unbound particles directly into radiation [49,50], this system includes the bound-state formation (BSF), ionisation and decay processes, as well as bound-to-bound transitions. The latter have been frequently (albeit not always) neglected in the past, with most previous computations focusing on the effect of the ground states on the DM density [2][3][4][5][6]. The reason has been two-fold: the capture into the ground state is indeed dominant in some models [51][52][53], while excited states tend to decay into radiation more slowly and be ionised rapidly until lower temperatures than the ground state, when the DM density is smaller, which limits their efficacy in depleting the DM abundance. However, several studies have now shown that this does not hold in general; capture into excited states can be faster than into the ground state, while de-excitations increase (decrease) the effective decay (ionisation) rate of the excited states [15,[19][20][21]. Considering excited levels is hence crucial for predicting the DM density accurately [20,21].
into radiation or via transitions to other bound levels, attain and remain always in a quasi-steady state in an expanding universe [4]. At high temperatures, this is ensured by the rapid BSF and ionisation processes, while at low temperatures the (effective) decay rate typically surpasses the expansion rate of the universe. As we shall show, the quasi-steady state that emerges depends on the interplay of the bound-state ionisation, transition and decay processes, and generalises the Saha ionisation equilibrium to metastable bound states. In this quasi-steady state, the equations that govern the DM freeze-out simplify significantly. More generally, this extension of Saha equilibrium pertains to all systems with metastable bound states in any thermodynamic environment.
Considering the above, in this work, we develop a general formalism for the DM thermal decoupling in the presence of an arbitrary number of metastable bound states, allowing for transitions among them. Under the assumption that the bound levels are in a quasi-steady state, we prove that the coupled system of equations reduces to a single Boltzmann equation that has the standard form of a Boltzmann equation for DM freeze-out with an attractor solution and an effective DM depletion cross-section. We show that the attractor solution remains the standard equilibrium solution for the free species in the absence of bound states and derive the effective DM depletion crosssection in closed form, in terms of the bound-state ionisation, transition and decay rates, for an arbitrary number of bound states. This single object encapsulates the system's entire dynamics, and provides a simple way to incorporate BSF in DM freeze-out computations and public DM codes.
The same object also specifies how Saha equilibrium varies between the two known regimes of stable and rapidly decaying bound states, and allows insight into the effect of bound-to-bound transitions. We prove that, despite the complex interplay of the various processes involved, boundto-bound transitions can only increase the effective DM depletion rate. This is an important result that renders excited bound levels more important than previously anticipated. In addition, in the limit of very rapid transitions, we show that the effective DM depletion rate becomes independent of the bound-to-bound transition rates.
We note that our formalism is independent of the microphysics that underlies the relevant processes. Depending on the model, bound states may form, get ionised and transition between each other radiatively [2-6, 14, 15, 19-21, 51, 52], via scattering on the relativistic plasma [16][17][18][19][20] or other re-arrangement processes [54,55]. Our parametrisation encompasses all possibilities, and our analysis pertains to all scenarios that feature metastable bound states. This paper is organised as follows. In section 2, we introduce the coupled system of Boltzmann equations that describe the DM freeze-out in the presence of bound states, reduce it to a single equation with an effective DM depletion cross-section that we derive in closed form, and discuss the generalised Saha ionisation equilibrium that emerges. In section 3, we examine the effect of bound-to-bound tran-sitions and show that they can only increase the effective DM depletion cross-section. We illustrate our results in a toy model in section 4, and conclude in section 5.

Coupled system
For generality, in the following, we consider the possibility of several co-annihilating species that share the same quantum number that stabilises DM. Bound states can form from any combination of these species. Let Y j ≡ n j /s and Y B ≡ n B /s be the number-density-to-entropy-density ratios of the free species j and the bound state B respectively. As is standard, we use the time parameter where m can be taken to be the mass of the lightest coannihilating species. The energy and entropy densities of the universe are ρ = (π 2 /30)g * ρ T 4 and s = (2π 2 /45)g * S T 3 , where g * ρ and g * S stand for the corresponding degrees of freedom (dof) of the universe respectively. We define [50] g The evolution of Y j and Y B is governed by the coupled Boltzmann equations where λ ≡ π 45 m Pl m g The equilibrium densities in the non-relativistic regime are where E B < 0 is the binding energy of the B bound state. The effective dof depend on the temperature as follows with g i and g B being the actual dof of the free species and bound states. We defined the fractional mass differences We are ultimately interested in the total yield of all coannihilating species whose equilibrium value is with Clearly, Y eq j = (g j,eff /g eff )Y eq . Note that the bound states are metastable and their abundance becomes eventually negligible, so we do not include them in eq. (8).
In eqs. (3), Γ dec B , Γ ion B→ij and Γ trans B→B are, respectively, the rates of B decay into radiation, B ionisation (a.k.a. dissociation) into ij, and B transition into another bound level B . The rates Γ j→i describe the transitions between free particles, due to decays, inverse decays and/or scatterings on the thermal bath; overall, these processes do not change the DM number density, but retain kinetic and chemical equilibrium among the dark species. Note that in eqs. (3), we must use the thermally averaged rates, Γ . The thermal average introduces Lorentz dilation factors for decay processes, which however are insignificant in the non-relativistic regime, as well as Bose-enhancement factors in the case of bound-to-bound transitions and capture or ionisation processes. The thermally-averaged rates and cross-sections of inverse processes are related via detailed balance that we have already employed in writing eqs. (3), Typically the various co-annihilating species either down-scatter or decay into the lightest one, that then makes up all of the DM today. Assuming that m has been chosen to be its mass, the fractional DM relic density is 2 where Y ∞ is the final value of the total yield (8) and f = 1 or 2 if the lightest species is self-conjugate or non-selfconjugate respectively. 3 In eq. (12), s 0 2839.5 cm −3 and ρ c 4.78 · 10 −6 GeV cm −3 are the entropy density and critical energy density of the universe today [56].

Effective equation
The system of Boltzmann eqs. (3) is computationally expensive to solve. However, it is possible under reasonable assumptions that remain valid throughout the DM chemical decoupling from the primordial plasma, to reduce it to a single effective Boltzmann equation, as we shall now see. In particular: (i) We assume that the i ↔ j conversion processes remain sufficiently fast to keep the various free species in chemical equilibrium; this implies where w is the same for all i, and in general depends on x; it essentially reparametrises the chemical potential of the free species, w = exp(µ free /T ). The assumption (13a) appeared first in the context of DM in ref. [49]. It is satisfied in many particle physics models; however, if species carrying the quantum number that stabilises DM are not in chemical equilibrium with the latter, then they have to be treated separately (see, e.g., [57,58]).
(ii) At high temperatures, BSF and ionization processes are very efficient, while later on decays into radiation, directly or via bound-to-bound transitions, dominate. In both cases, the rates that appear in the right-hand side of eq. (3b) remain very large such that the bound-state densities attain a quasi-steady state where This assumption was first proposed in ref. [4], and yields a system of linear equations for Y B whose solution can be re-employed in eq. (3a). We will discuss the validity of this approximation in section 2.4, using the solution that we obtain in the following.
Under the assumption (13a), the Boltzmann eqs. (3a) imply that the total yield (8) obeys the equation where we note that the i ↔ j conversion terms have canceled each other, as they do not change the total density. Using the definitions of the previous section, it is straightforward to show that and thus bring the first term of eq. (14) in the canonical form, with Y eq being the attractor solution and i,j (g i,eff g j,eff /g 2 eff ) σ ann ij v rel being the effective annihilation cross-section [49].
Under the assumption (13b), it is also straightforward to express the bound-state term of eq. (14) in terms of the actual densities of the free species only and the equilibrium densities of both the free and bound states. However, the attractor solution and the form of the effective cross-section that arise from this term are not immediately evident. Using both the assumption (13b) and detailed balance, in the following we prove that the attractor solution for Y remains Y eq and determine the bound-state contribution to the effective cross-section in closed form.
We begin with defining the total ionisation and boundto-bound transition rates, as well as the total interaction rate for the bound level B, Equations (11a) and (16a) allow to define the total capture cross-section into a level B, Then, we define the column vectors of the bound-state yields Y and Y eq , 4 4 Y and Y eq should not be confused with the yields Y and Y eq of the unbound species; we never use matrix notation for the latter in the present work.
as well as the square matrices Γ dec , Γ ion , Γ trans , Γ tot , T, We note the opposite order of the indices on the left and right sides of eq. (18e). No summation over repeated indices is implied in any of eqs. (18). Clearly, Γ dec , Γ ion , Γ trans and Γ tot are diagonal, while T has only off-diagonal non-vanishing elements. Using this matrix notation and under the assumption (13b), the Boltzmann eq. (3b) for the bound species yields the linear system We must now solve eq. (19) for Y. First, we consider the detailed balance eq. (11b) for bound-to-bound transitions, which in matrix form reads Summing over B , or equivalently We also note that since Γ tot is diagonal, Γ tot − T and its inverse also possess the property (20a), namely By virtue of the uniqueness of the inverse, P =P, and eq. (20e) follows.
Combining eqs. (19) and (20c) to eliminate Γ ion Y eq , and solving for Y, we find Considering the B element of eq. (21), we obtain where in the second step we took into account eq. (20e) and that Γ dec is diagonal. Equation (22) can be employed directly into the Boltzmann eq. (14), which with the help of (15), takes now the canonical form with the effective cross-section where g i,eff and g eff are given by eqs. (6a) and (10), and we define the efficiency factors Equation (25), in the context of eqs. (23) and (24), is one of the main results of the present work. It furnishes a general proof that the system of coupled Boltzmann eqs. (3) can be reduced into a single one, under the assumptions (13); this had previously been shown only in specific cases with a small number of bound states, while neglecting [4] or considering [20,21] bound-to-bound transitions. The closed form of eq. (25) renders it suitable for incorporating a large number of bound levels to thermal decoupling calculations.
The efficiency factors are bounded by where the limiting values are discussed in more detail in section 2.3. However, the inequality (26) is not immediately evident from the definition (25). To prove it, we first note that eq. (20c), expressed as Then, considering eqs. (25) and (27), it suffices to show that (Γ tot − T) −1 0 in order for (26) to be true. This i.e. all its off-diagonal elements, −T, are less than or equal to zero, and (ii) there exists a vector v > 0 such that (Γ tot − T) v > 0; indeed, v = Y eq satisfies this requirement as per eq. (20c) [59, theorem 1, condition K 33 ].

Generalisation of Saha ionisation equilibrium to
metastable bound states To obtain insight in the dynamics of systems with metastable bound states, we now examine more carefully the evolution of the bound-state densities, using the results obtained in the above. This will ultimately also allow us in the following subsection to determine the range of validity of our approximation. Equation (22) can be re-expressed using (25), as where for generality we used the number densities n = sY , and set w = n free /n eq free = exp(µ free /T ), with 'free' denoting any of the unbound co-annihilating species, presumed in chemical equilibrium with each other [cf. eq. (13a)]. Using eqs. (16d) and (27), eq. (28) can be also expressed in the following form that will be useful below, Equation (28), together with the expression (25) for r B in terms of the inherent parameters of the theory, generalise Saha equilibrium to metastable bound states. In terms of chemical potentials, it reads We discern two limits of eqs. (28) and (29).
• At temperatures much larger than the binding energies, T |E B |, the ionisation rates are large, consequently r B 1. Thus n B /n eq B (n free /n eq free ) 2 is independent of B. This implies that the bound states are in chemical equilibrium with the free species, µ B = 2µ free , as well as among themselves, however the equilibrium is maintained via the scattering-tobound transitions and does not depend on the boundto-bound transitions or bound-state decays. In terms of the bound-state densities, this equilibrium is similar to that of stable bound states. Nevertheless, the DM depletion due to BSF may be sizeable even during this period of r B 1.
• In the limit r B → 1, typically attained at T |E B |, the bound-state densities formally approach their equilibrium values with zero chemical potential, However, if the free species have developed a large chemical potential such that 1−r B exp(−2µ free /T ), it is possible that the BSF processes prevent n B from reaching n eq B . Neglecting transitions for simplicity, eq. (28a) yields n B n eq B + n 2 free σ BSF B v rel /Γ dec B . Considering that Y free = n free /s decreases only slowly after DM freeze-out, and taking into account the form of the Coulomb Sommerfeld factor for the crosssections that defines their asymptotic behaviour at late times, S = 2πζ/(1 − e −2πζ ) with ζ = α/v rel , we estimate that the second contribution to n B decreases roughly as  (29) between the chemical potentials of metastable bound states and free species, µB and µ free respectively, and rB, given by eq. (25) in terms of the bound-state decay, ionisation and transition rates. rB = 0 reproduces the standard Saha equilibrium for stable bound states, while rB 1 corresponds to bound states that decay much more rapidly than they get ionised.
between the incoming free particles in the BSF processes is attractive (α > 0) or repulsive (α < 0) respectively. In either case, the second term in n B continues to dominate over the first at low temperatures. Equation (28) relates the densities of free and bound species in the entire continuum between the above limits. The relation (29) between the chemical potentials is shown in fig. 1. While eqs. (28) and (29) were derived in the context of an expanding universe, they hold in any (homogeneous and isotropic) thermodynamic environment where metastable bound states attain a (quasi-)steady state.

Validity of the steady-state approximation
We now return to the steady-state approximation (13b) in order to determine its regime of validity. Using the solution (28) for the Y B /Y eq B ratio, the Boltzmann eq. (3b) for the bound level B can be re-written after trivial rearrangement of the terms, as where From the definition (25) of the r B factors, we find that where in the second step we separated the sum into B = B and B = B. With this, w B simplifies to i.e. the attractor solution of eq. (30) is indeed (28), as expected by construction. For the validity of our approximation, we must now require that for the attractor solution (28), the left-hand side of eq. (30) is much smaller than each term in the right-hand side, i.e. |d(w where H is the Hubble parameter. 5 Let us now estimate d ln(w B Y eq B )/dx. This of course depends on the evolution of w, which is determined by the Boltzmann eq. (23). Before DM begins to freeze-out, i.e. at x 30, w 1 hence w B 1 and d ln(w B Y eq B )/dx −2. When freeze-out begins, w starts growing exponentially with x, roughly w ∝ x −3/2 e x (though bound states slow down this growth since they prolong the decoupling.) We discern two regimes, following the discussion in section 2.3. As long as the temperature is much higher than the binding energies, T |E B |, then r B 1 and w B w 2 , hence d ln(w B Y eq B )/dx ∼ |E B |/m − 3/(2x) −3/(2x). At T |E B |, r B → 1 and we recall the discussion in section 2.3 stemming from eq. (28a). Collecting everything, we estimate that At T |E B |, the ionisation rates are typically very large, such that the conditions (33) are easily satisfied. 6 If the capture into a particular bound level, hence also the ionisation rate of that level, happen to be suppressed, for example due to a symmetry, then the bound-to-bound transitions may suffice to populate that level and satisfy the condition (33). This includes up-scattering processes, which are not Boltzmann suppressed at high temperatures.
At T |E B |, the ionisation rates become exponentially suppressed. However, by then typically Γ dec B H, unless a particular bound state is (nearly) forbidden to decay directly into radiation due to a symmetry. In this case, it may still effectively decay via transitions to other levels, thereby satisfying (33). This dynamics is very different from that encountered in the case of stable bound states, e.g. Hydrogen recombination in the early universe, where the steady-state conditions cease to be satisfied at T |E B | and the bound-state densities freeze-out. 7 Formally, the above considerations imply that the solutions (28) are good approximations if the conditions (33) are satisfied for all bound levels, since (28) was employed in deriving the Boltzmann eq. (30) for B, to account for the densities of all levels B = B. However, if there are bound levels that do not satisfy (33), then these levels decouple from the rest of the system. Hence, the densities of all other levels, as well as the depletion of the DM density, must not be affected significantly by them. We shall now show that the prescription for computing the DM density given by eqs. (23) to (25) can still be safely used, even though the number densities of the decoupled levels may deviate from (28).
We first consider the contribution to the effective DM depletion cross-section (24) due to capture into a levelB that does not satisfy the condition (33). Even though incorrectly estimated, this contribution would be limited by where we took into account that rB 1 and n eq B /n eq ∼ e −x 1. The above implies that the estimated effect of capture into theB level on the relic density would be indeed negligible, as it should.
Next we examine the effect of the transitions to and from levelB on σ eff v rel . Using the derivative of σ eff v rel with respect to the transition rates that will be computed in section 3.4, eq. (51), we find that the change in σ eff v rel due to a non-zero Γ trans 7 Of course, as is well known, the Hydrogen recombination in the early universe is further delayed due to the very small baryon-tophoton ratio.
where we employed similar considerations as before. Thus, the estimated effect on the relic density is again negligible.
We conclude overall that the effective method consisting of eqs. (23) to (25) may be used for the computation of the DM density without much concern about whether the steady-state condition always holds. On the other hand, the validity of the generalised Saha eq. (28) relies on the condition (33), which should be ensured for general applications.

No transitions
If the transitions between bound states are negligible or outright forbidden, Γ trans where we took into account that since Γ dec and Γ ion are diagonal, then ( This result was first derived in ref. [4] by considering a single bound level, and the accuracy of the approximation with respect to the solution of the coupled Boltzmann equations has been checked in [2,[29][30][31].

Rapid transitions
We now consider the possibility that the transitions between bound states are very rapid, in particular Γ trans In Appendix A, we show that in this limit, i.e. (Γ tot − T) −1 is independent of the column index. This renders the efficiency factors r B of eq. (25) independent of B, Equation (39) signifies that when the bound states are connected, only the averaged decay and ionisation rates, and consequently the averaged BSF cross-section, are important. This is consistent with and generalises the findings of [20], where transitions between bound states in a twolevel system were considered. Looking back at eq. (28), it is evident that as long as transitions are rapid, the ratios Y B /Y eq B become independent of B. This means that the bound states are in chemical equilibrium with each other.
Equation (39) constitutes another major result of the present work that can be readily employed to study systems with a large number of bound levels where bound-tobound transitions are rapid.

Rapid versus no transitions
Comparing eqs. (37) and (39), we find after elementary manipulations that bound-to-bound transitions increase the contribution to the effective cross-section (24) due to the capture into a bound level B, if Clearly, this condition is satisfied by some bound states and not by others. For higher angular momentum states, Γ dec B tends to be suppressed by higher powers in the coupling, and for excited states more generally Γ ion B tends to remain significant until later. The condition (40) thus tends to be more easily satisfied for the excited states, which tend to have lower decay-to-ionisation ratios.
The above may seem to suggest that the net effect of bound-to-bound transitions on the effective cross-section (24) depends on the relative strengths of the various BSF cross-sections and decay rates, since transitions between two levels affect their corresponding r B factors in opposite directions. However, the overall effect of transitions is in fact to enhance the DM depletion. Indeed, the total BSF contributions to the effective cross-section in the limits of sections 3.1 and 3.2 are With the crosssections expressed in this form, in Appendix B we show, using the Cauchy-Bunyakovsky-Schwarz inequality, that σ eff v rel no trans σ eff v rel rapid trans .
The equality holds when all bound states have the same decay-to-ionisation rate [cf. condition (40)], which is of course an unrealistic case. However, eqs. (41a) and (41b) approach each other in the limits, Γ dec This reduces to where the total cross-section for capture into a bound level B, σ BSF B v rel has been defined in eq. (16d). The former is realised at high temperatures, T |E B |. The effective cross-section is independent of the actual BSF crosssections, but grows exponentially until it saturates to the latter at low temperatures, T |E B | [13]. In the passage between these two regimes, bound-to-bound transitions matter. For systems that involve many bound levels, and the binding energies and decay rates span a large range of values, this regime may have a significant duration, thereby rendering bound-to-bound transitions very important. This is particularly evident if one or more bound levels have sizeable BSF cross-sections but vanishing or very suppressed decay rates, typically due to a symmetry. Neglecting transitions, such levels do not contribute to the depletion of the DM abundance, while considering transitions, their contribution may be very significant. Such an example can be found in refs. [19,20].

The effect of bound-to-bound transitions globally
While the discussion above offers important insight in the effect of bound-to-bound transitions on the depletion of the DM abundance, it does not unambiguously prove that this is monotonic. To examine the impact of transitions beyond the two limiting cases of sections 3.1 and 3.2, here we consider the variation of the efficiency factors r B and the effective cross-section σ eff v rel , with respect to the bound-to-bound transition rates. All indices below refer to bound states, and no summation over repeated indices is implied.
Starting from eq. (25) for r B , The derivative of a matrix can be expressed in terms of the derivative of its inverse. In particular, Based on the definitions (18) of Γ tot and T, Equations (46) and (47) and the symmetry property (20e) of (Γ tot − T) −1 imply Note that both sides of eq. (48) are symmetric in the interchange k ↔ . With this, and recalling the definition (25) of r B , eq. (45) becomes From eq. (49) we deduce that dr B /dT k may be either positive or negative, which is consistent with discussion around the condition (40). Next, we examine the variation of the effective crosssection (24), Recalling eq. (27), the above becomes This unambiguously proves that bound-to-bound transitions can only enhance the depletion of the DM abundance.

Toy model
While bound-state processes in realistic particle physics models involve considerable complexity that is beyond the scope of this work, in order to illustrate some of our results here we introduce a toy model. We consider in particular a three-level system with one free species χ of mass m and dof g χ = 1 and two bound levels 1 and 2, with masses m 1 = 2m − |E 1 |, m 2 = 2m − |E 2 | and dof g 1 = g 2 = 1; E 1 < E 2 < 0 are the binding energies. To avoid reference to dimensionful quantities, we define for B = 1, 2, where λ is defined in eq. (4a). Clearly, ε 1 > ε 2 > 0. The equilibrium densities of the three levels are The detailed balance eqs. (11a) and (11b) imply The effective Boltzmann eq. (23) becomes with For our system, where we set 1 − ε 1 1 − ε 2 1 for simplicity, since ε B 1 for weakly bound states. In the limit of very rapid transitions, γ trans 2→1 γ dec B , γ ion B , the above become For the numerical implementation of the model, we fix g * S = 10 2 for simplicity. To imitate a system that freezesout at x ∼ 30 in the presence of direct annihilations only, we set τ ann = 10 15 , and pick τ BSF 2 = 10τ BSF 1 = 10 2 τ ann to ensure that BSF is significant with respect to direct annihilations and that capture to excited states is significant with respect to capture to the ground state. We set ε 1 = 10 ε 2 = 10 −3 , to ensure two well separated and weakly coupled bound states, and consider three different cases for the bound-state decay rates, summarised in table 1. Note that the values ε ∼ 10 −3 and γ dec ∼ 10 11 correspond roughly to those of the ground state of a Coulomb potential with α ∼ 0.1, for m ∼ TeV.  For the models of table 1, we present the evolution of τ eff and the χ abundance in fig. 2. Comparing the cases γ 2→1 = 0 and γ 2→1 larger than all other rates, we see that transitions enhance τ eff in all cases. The effect is more pronounced when the excited state decays via down-scattering to the ground state (case A). When the ground state decays via up-scattering to the excited state, the effect of transitions is small albeit non-zero (case B). If the direct decay rates of both bound levels are significant, bound-tobound transitions have a moderate effect (case C). Figure 3 shows the monotonic decrease of the final density with increasing transition rate, in agreement with the result of section 3.4. We find that Y no trans final /Y max final 2.82, 1.035, 1.25 for cases A, B and C respectively.

Conclusion
The emergence of long-range effects around and above the TeV mass regime within the thermal-relic DM scenario impel a thorough reconsideration of DM freeze-out. It is now well established that the formation and decay of metastable bound states can alter the relic density of stable species even by orders of magnitude, thereby changing the predicted DM mass and couplings to other particles and affecting all experimental signatures. Designing targeted experimental searches for DM and accurately inter-preting the experimental results therefore necessitates a robust calculation of the DM freeze-out that will account for all bound-state effects. This is particularly important as experiments are now at the onset of exploring the multi-TeV energy scale.
In this work, we have presented a framework for the computation of the DM relic density when metastable bound states exist in the spectrum of the dark sector. This requires in general solving a coupled system of Boltzmann equations that describes the interplay of direct annihilation, bound-state formation, ionisation, transition and decay processes. Here we have proven in full generality that this system can be reduced to one Boltzmann equation of the standard form for DM freeze-out, with the entire dark sector dynamics encoded in a single object, the effective DM depletion cross-section, that we computed in closed form [cf. eqs. (23) to (25).] This result allows to easily incorporate bound-state effects in freeze-out computations, including in public codes.
Of particular interest is the importance of bound-tobound transitions for the DM depletion via BSF. We have shown that overall, bound-to-bound transitions can only increase the DM depletion rate [eq. (51).] This suggests that excited states may play a more important role, thus DM may be heavier than previously anticipated. In fact, because the capture into excited states renders higher partial waves relevant for the DM freeze-out [2], this result supports the conclusion that partial-wave unitarity does not constrain the mass of thermal relic DM as previously considered [2].
We finish by noting that our results generalise the Saha ionisation equilibrium to the case of metastable bound states. We have derived in closed form the relation between the densities of the bound states and the unbound species, in terms of the inherent parameters of the theory under consideration: the bound-state ionisation, transition and decay rates, combined into a single parameter for each bound state, the efficiency factor [cf. eqs. (28) and (29)]. This result can be adapted to and employed in studying any system that features metastable bound states.

Appendices Appendix A. The limit of rapid transitions
We want to determine the inverse of the matrix M ≡ Γ tot − T, in the regime where the transition rates between bound states are much larger than the decay and ionization rates. We write M as Then, adding to the j-th column all other columns, k = j, each multiplied by Y eq k /Y eq j , the j element for The above operations leave the determinant invariant, as well as all entries not contained in the i-th row or j-th column. Then, expanding with respect to either the i-th row or the j-th column, one can notice that to leading order in (and up to a sign), the determinant is given by the

Appendix B. Rapid versus no transitions
We consider two series a n and b n , and the expressions . (B.1) For our purposes, S n and L n correspond to eqs. (41a) and (41b) respectively, with a B ≡ Γ dec The difference between the two expressions is Rearranging further the above, we obtain .

(B.3)
Defining u j ≡ a j / a 2 j + b 2 j and v j ≡ a 2 j + b 2 j , this becomes