Bound-state effects on dark matter coannihilation: Pushing the boundaries of conversion-driven freeze-out

Bound-state formation can have a large impact on the dynamics of dark matter freeze-out in the early Universe, in particular for colored coannihilators. We present a general formalism to include an arbitrary number of excited bound states in terms of an effective annihilation cross section, taking bound-state formation, decay and transitions into account, and derive analytic approximations in the limiting cases of no or efficient transitions. Furthermore, we provide explicit expressions for radiative bound-state formation rates for states with arbitrary principal and angular quantum numbers $n,\ell$ for a mediator in the fundamental representation of $SU(3)_c$, as well as electromagnetic transition rates among them in the Coulomb approximation. We then assess the impact of bound states within a model with Majorana dark matter and a colored scalar $t$-channel mediator. We consider the regime of coannihilation as well as conversion-driven freeze-out (or coscattering), where the relic abundance is set by the freeze-out of conversion processes. We find that the region in parameter space where the latter occurs is considerably enhanced into the multi-TeV regime. For conversion-driven freeze-out, dark matter is very weakly coupled, evading direct and indirect detection constraints but leading to prominent signatures of long-lived particles that provide great prospects to be probed by dedicated searches at the upcoming LHC runs.


I. INTRODUCTION
Thermal freeze-out of dark matter has proved to be a successful framework for explaining the measured dark matter abundance in the Universe. However, the sizeable couplings of dark matter to the Standard Model (SM) particles required in the simplest realizations of this mechanism have been put under pressure by experimental null-results at colliders [1], direct [2] and indirect [3] detection experiments. Hence, fulfilling the relic density constraint often requires the exploration of 'exceptional' [4] regions, e.g. the region where coannihilation effects increase the effective annihilation rate [5].
Such effects commonly occur in models with a so-called t-channel mediator, where the mediator is taken to be odd under the Z 2 -parity that stabilizes dark matter and for relatively small mass splittings between the mediator and the dark matter particle. Prominent and wellstudied examples are the sfermion coannihilation regions in the minimal supersymmetric standard model (MSSM), see e.g. [6][7][8]. They have, in turn, motivated a wide range of phenomenological studies of t-channel mediators in the simplified model framework exploring different spin assignments and a wide range of coupling strengths [9][10][11][12][13][14].
While the presence of coannihilating mediators can increase the effective dark matter annihilation rate, toward small mass splittings, mediator pair-annihilation alone can become so efficient that dark matter is rendered under-abundant (seemingly) independent of the dark matter coupling. However, this conclusion is only valid if dark matter remains in chemical equilibrium with the mediator during freeze-out through efficient conversions. Dropping this assumption opens up a cosmologically viable part of the parameter space where the relic density is set by conversion-driven freeze-out [15] (or coscattering [16]). 1 In this scenario, thermal decoupling is initiated by the breakdown of efficient conversions between dark matter and the coannihilating partner(s). The required coupling is several orders of magnitude smaller than the one required to initiate the breakdown of dark matter pair-annihilations. This is due to the significantly larger number density for light standardmodel initial states in conversion processes with respect to the dark matter number density.
The boundary between the coannihilation and conversion-driven freeze-out region marks a significant change in the phenomenology within the parameter space of a given model. While the former is characterized by sizeable couplings that give rise to observable signals in conventional dark matter searches, the latter is largely immune to constraints from (in)direct detection but predicts long-lived particles with typical lifetimes of the order of millimeters to meters to be searched for at the LHC. The conversion-driven freeze-out region was unex-plored terrain for a long time, and often flagged underabundant when displaying the viable parameter space in terms of masses, see e.g. [9]. Recently, it has been studied in various contexts [17][18][19][20][21][22] and often constitutes one of a few regions still allowed within a given model [14].
For electrically and color-charged coannihilators -interacting via massless force carriers -non-perturbative effects such as Sommerfeld enhancement and, in particular, bound-state formation can play an important role in dark matter freeze-out. Radiative bound-state formation has been studied for a variety of dark matter models and for general unbroken Abelian and non-Abelian gauge theories [23][24][25][26][27]. The latter is related to earlier results for quarkonium formation inside the quark-gluon plasma obtained in potential nonrelativistic quantum chromodynamics (pNRQCD), see e.g. [28,29]. Recently, nextto-leading-order (NLO) finite temperature corrections of the general singlet-adjoint dipole interactions have been computed [30].
While it has been shown that bound-state formation effects provide sizeable corrections to the effective annihilation cross section for a coannihilation scenario with a colored mediator [26,31,32], it has widely been overlooked that their effects become considerably more relevant for scenarios with small dark matter couplings such as conversion-driven freeze-out. As a consequence of the small coupling, freeze-out is a prolonged process and the mediator annihilation down to significantly smaller temperatures (i.e. later times) becomes important. This increases the relevance of bound-state effects further prolonging the freeze-out process. Furthermore, studies have focussed on the effect of the ground state, while excited bound states become (increasingly) relevant toward smaller temperatures.
In this work, we extend the study of bound-state effects in several aspects: • First, we revisit the formulation of the underlying Boltzmann equations in the presence of excited bound states and derive a general framework for incorporating their effects in terms of an effective annihilation cross section of the coannihilator. This general form requires not only the knowledge of bound-state formation and decay rates but also the transition between the various excited states. However, we formulate two meaningful limiting cases considering fully efficient or non-efficient transitions, the latter of which is considered as a (conservative) benchmark scenario. This part is model independent and applies to any set of bound states in general.
• Secondly, focussing on the case of a colored mediator in the fundamental representation of SU (3) c , we derive general expressions for the bound-state formation rates of arbitrary n, (the principal and angular momentum quantum numbers of the bound state, respectively) and estimates for the transition in some cases. Furthermore, we discuss the impact of higher-order corrections to the bound-state formation and decay rates.
• Finally, we assess the impact of bound states for a colored t-channel mediator model and study the phenomenological consequences of bound-state effects in the coannihilation and conversion-driven freeze-out region. In particular, we observe a drastic shift in the boundary between the two regimes, greatly enlarging the latter region. These considerations allow us to assess the importance of the various corrections in the prescription of bound-state effects studied here.
The remainder of the paper is structured as follows. In Sec. II we introduce the considered benchmark model and review the Boltzmann equations that describe both the coannihilation and conversion-driven freeze-out scenario. In Sec. III we develop our general formalism to include bound states and discuss various limiting cases analytically. In Sec. IV we compute the involved rates for a colored mediator in the fundamental representation of SU (3) c . Section V is dedicated to the phenomenological application before concluding in Sec. VI. Appendices A and B contain further details of the computation of bound-state formation cross sections and discuss NLO QCD effects, respectively.

II. MODEL AND CONVERSION-DRIVEN FREEZE-OUT
We consider a simplified t-channel model with a singlet Majorana fermion χ providing the dark matter candidate, and a colored scalar mediatorq that exhibits a Yukawa coupling involving χ and a right-handed SM quark q, The scalar mediatorq transforms as a triplet under SU (3) c , as a singlet under SU (2) L , and has hypercharge that is identical to the one of right-handed SM quarks. It gives rise to a t-channel annihilation diagram for a pair of χ particles, and the corresponding process χχ →qq leads to a relic abundance of χ via thermal freeze-out. If the masses of χ andq are of comparable size, coannihilation processes need to be taken into account as well, in particular mediator pair annihilation, which dominantly proceeds via the processqq † → gg. (Annihilation into a pair of quarks is p-wave suppressed.) Being a pure QCD process, its cross section is entirely determined by the strong coupling α s . Indeed, this contribution can be so large that the χ relic density falls below the measured dark matter abundance, independently of the value of λ χ [9].
However, this conclusion hinges on the assumption that χ andq are in chemical equilibrium during the freeze-out process, i.e. that the corresponding conversion rates are large compared to the Hubble expansion rate H during the freeze-out process. Since the rates of all conversion processes necessarily involve some power of the coupling λ χ , the assumption of chemical equilibrium can be violated if the coupling strength is small enough. In that case, the conversions have to be included along with (co-) annihilation processes in the Boltzmann equations. This scenario is known as conversion-driven freezeout [15], or coscattering [16].
In general, for the minimal t-channel model considered here, the coupled set of Boltzmann equations reads [15] where x = m χ /T and Y i = n i /s, with number density n i and entropy density s, with where M pl 2.4 × 10 18 GeV is the reduced Planck mass. Yq represents the summed contribution of the mediator and its anti-particle, leading to the various factors 1/2. Here, gq = gq † = N c = 3, and fq is the distribution function that is assumed to be identical for particles and antiparticles as well as all colors. The processes in the first line of each equation denote the usual (co-)annihilation processes into SM particles. The conversion terms in the second line of each equation can be split into processes of the formq → χ and qq † → χχ. The former case requires accompanying SM particles, and can be further decomposed into 1 → 2 and 2 → 2 processes, with where Γ is the decay rate at rest, and where n eq i = T /(2π 2 ) g i m 2 i K 2 (m i /T ) and K i denote modified Bessel functions of the second kind. Depending on kinematic constraints further 1 → 3, 1 → 4 or 2 → 3 process can be relevant, especially for a coupling to top quarks,q =t [12]. In the following, we focus on a coupling to bottom quarks,q =b, and include the processes stated in eq. (6).
The set of Boltzmann equations can describe both coannihilations in and out of chemical equilibrium, with well-known simplifications being possible in the former case by summing both equations [5]. Out of chemical equilibrium, the coupled set of equations needs to be solved. However, since the coupling λ χ is small in this case, all terms except for the ones involving σqq † v and Γ conv can be neglected for conversion-driven freezeout. The former process is considerably Sommerfeld enhanced for small relative velocities, due to the attractive potential generated by gluon exchange in the color singlet configuration of theqq † pair [10]. In addition, the same potential leads to the formation of bound states [25,26,32,33]. In this work, we improve the computations of the relic density in the coannihilation and conversion-driven freeze-out scenario by considering bound-state effects, including an exploration of the role of excited states.

III. INCLUDING BOUND STATES
Within the t-channel model, bound states ofqq † pairs in the color singlet configuration exist and can contribute to the freeze-out process. We consider an extension of the Boltzmann equation by including a set of bound states B i , enumerated by an abstract index i, and with g Bi in-ternal degrees of freedom. Within the model considered here, the bound states are characterized by their n and quantum numbers, i ≡ (n, ) and g B n = 2 + 1, but the discussion in this section applies to any set of bound states in general.
We add a Boltzmann equation for the yield Y Bi = n Bi /s for each bound state, taking into account ionization (or equivalently breaking) into an unboundqq † pair via gluon or photon absorption, direct decay of the bound state into SM particles, and transitions between two bound states. In addition, the collision term in the Boltzmann equation of the mediatorq picks up an extra term due to ionization and its inverse process, recombination [or equivalently bound-state formation (BSF)]. The changes in the Boltzmann equations compared to eqs. (2) and (3) are given by The ionization rate Γ i ion is related to the thermally averaged recombination cross section σ BSF,i v via the Milne relation originating from the detailed balance condition in thermal equilibrium. Indeed, the Milne relation ensures that the ionization and recombination terms drop out in the sum d(Yq + 2 i Y B,i )/dx, consistent with the conservation of the total number ofq andq † in the absence of decays. Note that in the non-relativistic limit where E Bi = 2mq − m Bi > 0 is the binding energy, and we used that Yq denotes the yield of the sum ofq andq † . In addition, detailed balance requires Also, here, we can see that transition terms drop out when summing the Boltzmann equations for all bound states, as required. Before discussing explicit expressions for corresponding rates in Sec. IV, we investigate generic features of the coupled set of equations.

A. Single bound state
We first recall the case of a single bound state B. In a typical cosmological setting, the ionization and decay rates (mediated by the strong interaction) are much larger than H. In this case, the density of bound states almost instantaneously adjusts to a quasistationary number (from the point of view of cosmological versus strong interaction timescales) that can be obtained by setting the left-hand side of the Boltzmann equation for B to zero, turning it into an algebraic equation [34]. For the case of a single bound state (dropping the index i and transition terms), one obtains Inserting this relation in eq. (10) yields the same form as eq. (3) but with the substitution This means it is sufficient to solve the Boltzmann equations forq and χ, while the impact of the bound state is captured by replacing theqq † annihilation cross section by the effective cross section.
In the limit H Γ dec Γ ion the ionization and recombination processes establish equilibrium between the bound state and unboundq (ionization equilibrium). The corresponding rates therefore drop out of the effective cross section, which only depends on the decay rate Γ dec , as can be seen using the Milne relation, eq. (11), The effective cross section increases exponentially with falling temperature, due to the energetic preference for bound states in equilibrium. This increase stops once the ionization rate, which itself becomes exponentially suppressed at low temperatures, falls below the decay rate, and ionization equilibrium breaks down. Therefore, at low enough temperatures, the regime H Γ ion Γ dec becomes relevant, for which In that limit, any bound state that forms decays almost immediately, and therefore the effective cross section is only sensitive to the recombination cross section σ BSF v .

B. Multiple bound states
Let us now generalize the previous findings to a set of bound states. When assuming as before that all relevant ionization, decay and transition rates are much larger than H, we obtain a set of coupled algebraic equations for the yields Y Bi from setting the left-hand sides of the Boltzmann equations (9) to zero. It can be written as where we used eq. (13) and introduced the total width of B i , From the structure of the Boltzmann equation, it is a priori not clear whether the impact of bound states can be captured by an effective cross section when inserting the solution to eq. (18) into the Boltzmann equation (10) forq. However, this turns out to be the case in general.
To see it, we rewrite eq. (18) in the form where we defined y i ≡ Y B,i /Y eq B,i and y ≡ Yq/Y eq q . Introducing the matrix the solution for the bound-state abundances reads Inserting it in the Boltzmann equation (10) forq indeed yields a contribution that has the form of the annihilation term, involving in particular a factor y 2 − 1. Therefore, provided the rates are large compared to the expansion rate, the impact of a set of bound states can in general be captured by an effective cross section, given by with The effective cross section, eq. (23), describes the impact of an arbitrary number of bound states on theq abundance, which can all individually be populated by recombination processes, decay into SM particles, and undergo a network of transitions among them, with the corresponding rates entering in the determination of R i . For a given setup, the R i can be determined numerically. Nevertheless, it is instructive to study two limiting cases analytically.

No transition limit
In the limit Γ i→j trans Γ i dec , Γ i ion we can neglect the transition terms, such that M ij → δ ij becomes the unity matrix, and the total width depends only on ionization and decay rates. The effective cross section becomes In the absence of transitions, each bound state therefore gives a contribution to the effective cross section that is analogous to the case for a single bound state, see eq. (15). In particular, each summand exhibits the limiting cases of ionization equilibrium (Γ i ion Γ i dec ) or instantaneous decay (Γ i ion Γ i dec ) in close analogy to the case of a single bound state.

Efficient transition limit
In the limit Γ i→j trans Γ i dec , Γ i ion , we expect that the transitions establish chemical equilibrium among the bound states, which is indeed a solution to eq. (18) in that limit. The most straightforward way to derive the effective cross section in that limit is to proceed similarly to the case of coannihilations [5], introducing and summing up all Boltzmann equations (9) for the B i , such that the transition terms drop out. Using (26) to write one obtains with effective ionization and decay rates Setting again the left-hand side of the Boltzmann equation (29) to zero, and inserting the resulting algebraic expression together with eq. (28) into eq. (10) yields where The result is similar in form to the case of a single bound state, eq. (15), but with the ionization and decay rates replaced by a thermal average over all bound states and the recombination cross section replaced by the sum. It turns out that obtaining this result directly from the general expression eq. (23) is tedious. The reason is that naively neglecting the ionization and decay rates in the total width would lead to a singular matrix M ij . However, by carefully expanding the abundances around the chemical equilibrium solution y i =const., and treating Γ i ion /Γ i and Γ i dec /Γ i as small, one ultimately arrives at the same expression, eq. (31).
We also note that using the Milne relation, eq. (11), for each bound state, one finds i.e. the summed recombination cross section and the effective ionization rate satisfy a generalized Milne relation. This implies that, in analogy to the case of a single bound state, within the regime of ionization equilibrium (Γ eff ion Γ eff dec ), the effective cross section becomes independent of the recombination cross section, and only depends on the effective decay rate. In the opposite limit Γ eff ion Γ eff dec of almost instantaneous decay, the decay rate drops out, and the effective cross section depends only on σ BSF v sum .

Ionization equilibrium
The limit of ionization equilibrium is somewhat orthogonal to the two limiting cases considered above. When ionization and recombination processes are assumed to be efficient enough to establish ionization equilibrium, the effective cross section approaches the universal form which is a straightforward generalization of eq. (16) and independent of ionization rates Γ i ion as well as transition rates Γ i→j trans . The reason is that efficient ionization and recombination processes establish chemical equilibrium with the unboundq particles in that case for each bound state. This means, in turn, that they are in chemical equilibrium among each other, such that the transition processes play no role for their relative abundances in that limit. This result agrees with the finding in [35], in which a set of bound states in ionization equilibrium is considered.
Indeed, it is easy to see that eq. (33) follows from both the effective cross section in either the limiting case of no transitions or the case of efficient transitions when assuming in addition that Γ i ion Γ i dec . Moreover, the fact that eq. (33) is even valid independently of the size of transition rates can be seen by noticing that the derivation presented in Sec. III B 2 relies only on the assumption of chemical equilibrium among the bound states, which is satisfied in ionization equilibrium.
Therefore, as long as ionization equilibrium holds, the effective cross section is only sensitive to the bound-state decay rates, independently of the size of transition and ionization rates.
In a realistic setup, the limiting assumptions made above may be too restrictive and at best hold only for a subset of bound states and a subset of the corresponding ionization, decay or transition processes. In this case, the effective cross section can be computed using the general result, eq. (23).

IV. RATES
While the discussion in the previous section was generic, we focus on the set of bound states and ionization, decay and transition rates that are relevant for the scalar mediatorq that carries hypercharge and transforms under the fundamental representation of SU (N c ) with N c = 3 in the following.
A heavy (mq Λ QCD ), non-relativisticqq † pair can be described by two wave functions ψ [R] , one for the color octet ( [8]) and one for the color singlet ([1]) configuration. They obey a Schrödinger equation with kinetic energy is the reduced mass, and potential in Coulomb approximation [26] with effective coupling strength Here C [R] 2 denotes the quadratic Casimir of SU (N c ) with 2 = N c = 3, and α s = g 2 s /(4π) is related to the strong coupling constant. Thus, The singlet configuration feels an attractive potential, while it is repulsive for the octet. Therefore, bound states exist for the singlet only. Note that we treat the m quantum number as an internal degree of freedom of the bound state in the Boltzmann equation, and therefore label the bound states by n and only.
In the Coulomb approximation, the bound states are described by hydrogen-like wave functions ψ [1] n m , with the fine-structure constant replaced by α eff [1] and the electron mass by the reduced mass µ. On the other hand, unbound scattering states ψ [R] p rel exist for both the octet and singlet, with wave functions containing the respective effective coupling strength (see App. A).

A. Ionization and recombination
The leading-order QCD process for bound-state formation is where the initial state corresponds to a scattering state in the octet configuration due to color conservation. The matrix element can be computed within pNRQCD analogously to hydrogen recombination [29,36], with a dipole interaction Hamiltonian of the form g s ωr·E where E = t a E a is the color-electric field, r is the relative coordinate, and is the energy difference of initial and final state, which corresponds to the energy of the emitted gluon in the non-relativistic limit. The thermally averaged ionization (or breaking) rate and recombination (or bound-state formation) cross section are given by [26] which can be checked to satisfy the Milne relation, eq. (11), with f g (ω) = 1/(e ω/T − 1). The recombination cross section can be expressed as [29] σ BSF,n v rel = ω 2πN 2 with the matrix element for the QCD process given by and | ψ [1] n |r|ψ [8] One can also consider the analogous electromagnetic process, that proceed from a color singlet scattering to bound state and matrix element obtained from the electromagnetic dipole interaction, where Qq = 1/3 is the electric charge ofq. We evaluate the strong couplings entering in the scattering (s) and bound (b) state wave function at renormalization scale of the typical momentum transfer related to bound and scattering states, respectively, using the notation For the strong coupling that enters via the interaction Hamiltonian we choose α BSF s = α s (µ MS = ω), evaluated at the gluon momentum scale. In contrast to [26], we use an identical scale choice for couplings entering either via Abelian or non-Abelian vertices. Within pNRQCD, the latter manifest themselves exclusively by the C A contribution to α eff s for the gluonic recombination process. For the binding energy we use Using a partial wave decomposition as well as an integral representation for the hypergeometric function entering the scattering-state wave function, and the generating function of the Laguerre polynomials contained in the bound-state wave function, we arrive at the following expressions for the bound-state formation cross sections via where we defined and Here, corresponds to the partial wave of the scattering state, which is constrained by the usual selection rule, and the radial part of the wave function yields the overlap integral This expression can be easily evaluated numerically. The result has the structure where s BSF n is a polynomial, with explicit expressions for n ≤ 3 given in Tab. I. For the 1s ground state our result agrees with [26], and for the 2s state it agrees with [29,30]. The result for 2p differs from the one given in [29,30] (by a factor 3 for the s-wave contribution with = 0, and a factor 3/2 for the d-wave contribution with = 2) but matches the result for hydrogen when translated to the electromagnetic case [36].
We show the functions S BSF n (ζ s , ζ b ), which are proportional to the bound-state formation cross section, in Fig. 1. For the figure we assume that the strong coupling entering in ζ s and ζ b is evaluated at a common renormalization scale, such that S BSF n depends only on the ratio α s /v rel . Furthermore, we show the sum over all = 0, . . . , n − 1 for a given n (solid lines), as well as the results for the s-orbital with = 0 (dashed lines). For α s /v rel 1 the bound-state formation cross section scales as (α s /v rel ) 4+2 for all n. The limit is given by where the second line arises from the = + 1 contribu- tion, and the first from = − 1 exists only for > 0. The contribution from = 0 orbitals therefore dominates for α s /v rel 1, as can also be seen by the convergence of solid and dashed lines for each n ≥ 2 in Fig. 1 in that limit.
In the opposite limit α s /v rel 1, where Here s BSF n | 4n−2 corresponds to the polynomial obtained when keeping only the terms with maximal combined power in ζ s and ζ b in s BSF n (ζ s , ζ b ), being 4n−2 , such that f BSF n depends only on the ratio ζ s /ζ b = α eff s /α eff b . Up to the different renormalization scale at which the effective couplings are evaluated, f BSF n approaches a constant for α s /v rel 1. The behavior at small relative velocities is therefore governed dominantly by the first factor in eq. (55). It exhibits a qualitatively different behavior depending on the sign of ζ s . For (qq † ) [8] → B [1] n + g, the repulsive potential relevant for the initial state implies ζ s < 0, leading to an exponential suppression for small relative velocities, S BSF n → 2π|ζ s |e −2π|ζs| f BSF n . For the electromagnetic process (qq † ) [1] → B [1] n + γ, both the initialand final-state wave function are sensitive to the attractive color singlet potential, such that in particular ζ s > 0, and S BSF n → 2πζ s f BSF n grows with ζ s ∝ α s /v rel . The different shape of S BSF n for the two processes can clearly be seen in Fig. 1. For the electromagnetic process, the combined contribution from all angular momentum states S BSF n decreases with increasing values of n, for all velocities v rel . On the other hand, for the strong process the exponential suppression at large ζ s leads to a maximum of S BSF n . Its position shifts to higher values of α s /v rel for excited states with increasing n. In addition, the value at the maximum increases with n. This indicates that excited levels become more and more relevant the smaller the relative velocity, i.e. the lower the temperature that is relevant for determining the relic density.

B. Decay
The leading decay process is due to annihilation of the constituents of the bound state into a pair of gluons, B n → gg. Here, we briefly review the derivation of the decay rate following [23], provide an expression for general n (for = 0) and discuss the role of higher-order corrections.
For a generic 1 → N decay process, B n → X 1 X 2 . . . X N the matrix element M n can be related to the usual Feynman matrix element for the process with N q → µ in the nonrelativistic limit, and boundstate wave function ψ n m ≡ ψ [1] n m in momentum space, normalized such that d 3 x|ψ n m (x)| 2 = 1 in position space. Here, K is the four-momentum of the bound state.
The bound-state decay rate is given by where m B n = 2mq − E B n 2mq, and |M n | 2 = 1 2 +1 m g X j |M n m | 2 is averaged over the 2 +1 states with different m, and summed over final-state degrees of freedom. Furthermore, the usual factor 1/S! is included if S particles in the final state are of identical type. For two-body decays at rest, the integration over the Lorentzinvariant phase space (LIPS) reduces to a factor 1/(8π).
At leading order in the small relative momentum q and in the non-relativistic expansion, with on-shell four-momentum K 2 = m 2 B n , and The wave-function at the origin is non-zero for orbitals with = 0 only, while the decay of bound states with orbital angular momentum would require keeping further terms in the expansion of M s in q, leading to a suppression of the decay rate of order q 2 /K 2 ∼ E n /mq. We therefore focus on the decay of = 0 states in the following.
The decay rate can also be obtained from an effective operator that describes the interaction of an = 0 bound state with a pair of gluons, with a scalar field Φ n that describes the B n, =0 bound state, and a form factor F (Q) ≡ Z n00 /Q 2 , where Q 2 ≡ p 1 · p 2 . The coefficient can be obtained by matching the matrix element for the two-body decay in the full and effective description.
Note that the matching does not require the gluons to be on-shell. Accordingly, the effective operator, eq. (64), can also be used to compute the 2 → 2 scattering processes of the form Bq → gq. Implementing the effective operator in MadGraph5_aMC@NLO [37], we checked that 2 → 2 processes can only compete with the boundstate decay for very early times, x < ∼ 5 − 10, for which the mediator is still in thermal equilibrium with the SM plasma. Hence, these processes are negligible for the dynamics of dark matter freeze-out considered here.
In contrast, NLO corrections to the two-body decay rate are potentially relevant since in ionization equilibrium, the impact of bound states on the effective cross section is determined predominantly by their decay rate, see the discussion in Sec. III. Following earlier results in the context of quarkonium [38][39][40], the virtual and real corrections to the B 10 → gg decay rate at O(α s ) have been computed for the comparable case of stoponium [41], resulting in: where n f is the number of light quarks, T F = 1/2, and b 0 = 11/3 C A − (1/3 + 4/3 n f )T F . The parameter δ 4q is either 1 or 0 depending on whether or not the four-point interaction ofq is introduced. In the simplified model considered here, δ 4q = 0, while in the MSSM, δ 4q = 1.
For the scale choice (63), we find that the correction (65) is reduced to a few percent rendering the leadingorder (LO) and NLO predictions fully compatible with each other within scale uncertainties, see App. B for a detailed discussion. For definiteness, we will consider the LO decay rate in the main results in the following.
The decay of bound states with > 0 is suppressed compared to those with = 0. Nevertheless, due to the large number of such states, it would be interesting to include them, which is, however, beyond the scope of this work. We note that the two-body decay into a pair of gluons is forbidden for = 1 states due to the Landau-Yang theorem. A decay channel that is possible for these states is into a pair of electrically charged particles, via an intermediate photon or Z-boson. Note that an analogous process with an intermediate gluon is forbidden by color conservation. Furthermore, the decay rate into gqq via an intermediate off-shell gluon also vanishes for = 1, as can be checked by expanding the matrix element in eq. (57) to first order in the relative momentum q. However, a decay into three gluons could be mediated by the strong interactions.

C. Transitions
Since bound states exist only in the color singlet configuration, transitions between energy levels cannot proceed via single gluon emission or absorption. In this work, we do not consider transitions involving two gluons, which can be mediated by the strong interaction. Instead, we provide a lower bound on the size of transition rates by considering the electromagnetic process which is allowed by color and charge conservation. The transition matrix element squared obtained from the electric dipole interaction is given by [36] where ω = |E n −E n | is the photon energy. The matrix element is averaged over m and m , The transition rate from higher to lower energy levels is given by Fermi's golden rule The rate of the inverse process of photoabsorption can be obtained from the detailed balance condition eq. (13). Using the hydrogen-like wave functions and the generating function of the Laguerre polynomials (see App. A) we find | ψ n |r|ψ n | 2 = δ , +1 + δ , +1 (2 + 1)(2 + 1) |I trans where with N n (κ) = κ 3/2 4(n − − 1)! n 4 (n + )! 2κ n . Here with effective strong coupling defined as in eq. (47) and evaluated for the respective energy level as indicated by the subscript. They differ only in the scale choice of the strong coupling constant, related to the typical Bohr momentum of the two energy levels. We checked agreement with various explicit expressions given for specific n and , and all n as well as for n = n in [36], when translating the result to the analogous hydrogen transition rates.

D. Effective cross section
Using the ionization, decay and transition rates discussed above we can compute the effective cross section, eq. (23), that encapsulates the impact of bound states on the freeze-out dynamics. The contribution to the effective cross section, eq. (23), due to bound states, is shown in Fig. 2 for various approximations as a function of x = m χ /T . In the left panel, we include bound states up to n = 6 and for all ≤ n − 1.
The contributions from individual n, states to the effective cross section are indicated by the colored and gray lines in the left panel of Fig. 2. For small x, the = 0 states dominate. The reason is that in this limit, ionization equilibrium holds, and the effective cross section is determined by the decay rate, see eq. (33). For large x, the contribution from each n, level becomes suppressed due to a combination of two effects: (i) the suppression due to the repulsive interaction in the scattering state discussed in Sec. IV A, and (ii) Boltzmann suppression for T E B n . Consequently, each individual contribution features a maximum. Its position shifts to the right for higher n. This implies that excited states dominate the effective cross section for large x. The larger x, the higher n have to be taken into account to obtain a converged result for the total effective cross section.
The line labeled "R i -solution" shows the total result obtained when including all rates as given above using the general expression eq. (23) for the effective cross section. For comparison, we show the limit of efficient transitions, eq. (31), as well as the limit of no transitions, eq. (25). For small x, that is, large enough temperature, all results agree and approach the ionization equilibrium result, eq. (33), that is also shown. The effective cross section can in this limit be written as That is, in ionization equilibrium, excited states lead to a 20% correction to the effective cross section. The factor in front of the sum is the ground-state contribution to eq. (33).
The impact of excited states is much larger for large x, where they give the dominant contribution. The precise value depends in this regime on the recombination as well as transition rates. The efficient transition limit provides an upper bound on the effective cross section (since all orbitals contribute), while the limit of no transitions provides a lower limit (only the bound states with a sizeable decay rate into SM particles contribute, being = 0 orbitals in our approximation). The actual effective cross section is therefore expected to lie in between these two limits. The "R i -solution" result taking into account the electromagnetic transition rates considered in this work can only be considered as illustrative since additional processes mediating further transitions are expected to play an important role. We therefore conservatively adopt the no-transition limit in our numerical analysis in the following.
The effective cross section in the no-transition approximation eq. (25) is shown in the right panel of Fig. 2. We show the result summed up to some maximum n, for n = 1, 3, 6, 10, 15, respectively. While each individual contribution becomes suppressed at large x, the summed result continues to grow with increasing x. The decline at very large x is due to the restriction to n ≤ 15. For x 10 5 , we consider the effective cross section with n ≤ 15 as converged. We leave an exploration of the full result including transitions to future work, and use the no-transition limit with n ≤ 15 as the default choice in the following. For a discussion of the impact of a certain class of higher-order corrections (related to collisional ionization and recombination processes and the associated virtual contributions) computed in [27,30] as  well as to bound-state decay we refer to App. B.

V. VIABLE PARAMETER SPACE
To determine the relic abundance, we solve the coupled set of Boltzmann equations (2) for Y χ and (3) for Yq. We compute the involved conversion and annihilation cross sections, σq k→χl (s) and σ χχ (s), σ χq (s), σqq †(s), respectively, with MadGraph5_aMC@NLO [37]. We take into account the leading conversions in α s and regularize m χ =1TeV, Δm=20GeV the soft divergence occurring in the processqg → χb (see the discussion in [15]) by introducing a thermal mass for the gluon [42]. To include the impact of bound states, we replace the annihilation cross section ofqq † pairs by the effective cross section, eq. (23). In addition, we include Sommerfeld enhancement in the contribution from direct mediator annihilation as described in [15]. Figure 3 exemplifies the effective cross section. The long-dashed curve ('pert.') shows the perturbative direct annihilation cross section while the short-dashed ('Som'.), dot-dashed ('BS, n = 1') and solid ('BS, n ≤ 15') curves display the effective cross section after successively including Sommerfeld enhancement, bound-state formation effects of the ground state and excited bound states up to n = 15 (in the no transition limit), respectively. In the following, we choose the latter for our main results. We also show the effective cross section under the assumption of ionization equilibrium in the limit of large n as the gray dotted curve ('ion-eq'). For two benchmark points in the conversion-driven freeze-out scenario, the evolution of the abundances is shown in Fig. 4. Because of the small coupling λ χ , the χ particle cannot annihilate efficiently by itself, and its abundance is reduced only due to conversions intoq. While the colored mediatorq starts to depart from thermal equilibrium at x > ∼ 25, the χ abundance already significantly exceeds the equilibrium value at this time. Subsequently, for x > ∼ 25, conversion processes -which are on the edge of being efficient -gradually transform χ intoq particles. This leads to a prolonged duration of the freeze-out dynamics, which can last until x ∼ O(10 2 −10 3 ). The mediatorq continues to annihilate and is in addition depleted due to bound-state formation.
The duration is further enhanced for a small relative mass splitting ∆m/m χ , which implies that the equilibrium abundances of the mediator and χ are comparable until x ∼ m χ /∆m even if the conversion processes were fully efficient, i.e. in the usual coannihilation scenario. Eventually, the mediators decay viaq → bχ, thereby transferring their remaining abundance to the population of χ particles. For the chosen value of the coupling λ χ for the two benchmark points shown in Fig. 4, the amount of conversions is sufficient to reduce the χ abundance to a final value that matches the observed relic density, Ωh 2 = 0.12 [43].
In Fig. 5, we show the coupling λ χ that is required to achieve Ωh 2 = 0.12 as a function of the mass splitting, ∆m, for fixed mass m χ = 1 TeV (left panel) and as a function of the dark matter mass, m χ , for fixed ∆m = 5 GeV (right panel). The drastic change in the coupling at ∆m 35 GeV and m χ 2850 GeV, respectively, is due to the transition between the conversiondriven freeze-out (to the left) and coannihilation regime (to the right), see below for details.
The gray lines in Fig. 5 show the impact on the relic density for various levels of approximation, relative to our fiducial choice with Sommerfeld enhancement and bound states up to n = 15. The relic density differs up to a factor of order 10 relative to the perturbative leading-order approximation, and for small ∆m/m χ . Relative to the case when including Sommerfeld enhancement, we find differences of up to a factor of order five. The gray line labeled 'BS, n = 1' corresponds to the case when including the ground state. As apparent from the relatively small deviation of this curve from one, excited states with n ≤ 15 only yield a comparably small correction for most of the shown parameter space.
The kink at the transition between the conversiondriven freeze-out and coannihilation regime that can be seen in most curves arises due to the sudden increase of the coupling at this point that causes χχ and χq annihilation processes to become relevant. Accordingly, in the latter regime, not only does the relative importance of non-perturbative effects on the effective mediator annihilation change but also the importance of the effective mediator annihilation with respect to χχ and χq annihilation changes. This causes the quicker decrease of seen in the left panel of Fig. 5. Here, we can also observe that all curves approach unity toward large mass splittings as both the Boltzmann suppression of the mediator abundance during freeze-out and the larger coupling, λ χ , diminishes the relative importance of the mediator annihilation.
In the parameter slice with m χ = 1 TeV, chosen in the left panel, freeze-out mainly occurs while the system is still close to ionization equilibrium. This can also be seen from the gray dotted curve showing the result assuming ionization equilibrium (for all n). It only deviates significantly for low ∆m where freeze-out extends to large x. For even smaller relative mass splittings, ∆m/m χ , considered in the right panel, this effect is even more pronounced as freeze-out extends to larger x (even in the coannihilation region). Here, the result for ionization equilibrium differs by orders of magnitude from the one of our fiducial choice reaching Ωh 2 /0.12 < ∼ 10 −4 in the considered range of m χ (outside the displayed range in Fig. 5).

A. Boundary between coannihilation and conversion-driven regime
In this section, we determine the part of parameter space of the model for which conversion-driven freezeout is relevant. For small mass splitting ∆m ≡ mq − m χ and mass m χ , theqq † annihilation process becomes very efficient, and would deplete the relic abundance below the observed dark matter density, if χ andq were in chemical equilibrium during freeze-out. Within the region of parameter space where this happens, the correct dark matter abundance can only be explained if the as- sumption of chemical equilibrium does not hold. The dynamics are described by conversion-driven freeze-out in this regime, and one obtains a viable relic density for couplings λ χ 1. On the other hand, for points in parameter space where theqq † annihilation cross section is small enough, the standard scenario of coannihilation yields the observed dark matter abundance, with λ χ ∼ O (1). The division between these regimes can effectively be obtained with high precision by solving the Boltzmann equation using the conventional coannihilation approximation in the limit λ χ 1. The relic density obtained in this limit matches the observed dark matter abundance along a line in the two-dimensional parameter space (m χ , ∆m), which we refer to as the boundary line.
Altogether, the correct relic density can be reproduced for any point within the two-dimensional parameter space for a suitable value of λ χ , via conversion-driven freezeout below the boundary line, and via conventional coannihilation above the boundary line. Note that the effectiveqq † cross section, eq. (23), including bound-state effects is relevant both in the coannihilation as well as the conversion-driven regimes, and therefore also for determining the boundary between them.
In Fig. 6, we show the boundary line in the (m χ , ∆m) plane obtained for various approximations, which successively include a number of effects. When using the perturbative tree-levelqq † annihilation cross section only, one obtains the line labeled 'pert.'. This is the result one would obtain when using standard tools for the relic density computation [44][45][46] without further modification. The line labeled 'Som.' is obtained when including Sommerfeld enhancement ofqq † annihilation, and this approximation has been used in previous works in the context of conversion-driven freeze-out with colored mediators [12,15]. 3 The regime of conversion-driven freeze-out extends significantly when including the bound-state effects considered in this work. The line labeled 'BS, n = 1' in Fig. 6 corresponds to including the contribution from the ground state only. Finally, adding excited states up to n = 15 within the default approximation discussed in Sec. IV D yields the thick solid line. We observe that the conversion-driven freeze-out region reaches to significantly higher values of m χ and also ∆m due to the impact of bound states.
Let us briefly comment on the role of excited states. For m χ /∆m < ∼ O(10 2 ), freeze-out dominantly takes place in the regime of ionization equilibrium. In that case, excited states lead to a correction of the effective cross section of order 20%, due to the additional available decay channels, see eq. (75). For m χ /∆m > ∼ O(10 2 ), the freezeout extends to lower temperatures. In this regime, a combination of two effects leads to a significant enhancement of the impact of excited states. First, since ionization equilibrium breaks down for the ground state, its contribution to the effective cross section drops. Secondly, the bound-state formation rate for excited states exceeds the one of the ground state by many orders of magnitude at low temperatures. Hence, excitations remain in ionization equilibrium toward smaller temperatures and dominate the effective cross section. 4 Potentially, the region of conversion-driven freeze-out could even become larger when including transitions between the bound states, which is beyond the scope of this work. To provide a maximal upper bound we show the result that would be obtained when assuming ionization equilibrium to hold during the entire freeze-out and including all n using eq. (75), indicated by the gray dotted line. The full result when including transitions is expected to lie significantly below this line, and above the solid line, cf. the respective results for σqq †v BS eff in the left panel of Fig. 2. For the regime where the gray and thick solid lines differ from each other, ionization equilibrium breaks down during the freeze-out. The boundary therefore becomes insensitive to uncertainties from transitions among bound states where both lines converge, i.e. for m χ < ∼ 2 TeV.

B. Coannihilation regime
While the main focus of this work is on the impact of bound states on conversion-driven freeze-out, we also assess the relevance in the coannihilation regime. As is already apparent from Fig. 6, bound states and Sommerfeld enhancement have a significant impact on the boundary, and therefore on coannihilations as well. In Fig. 7, we show the contours in the (m χ , ∆m) plane for which freeze-out in the coannihilation regime yields the correct dark matter relic abundance for three values of the coupling, λ χ = 0.169, 0.5, 1, respectively. The former choice is motivated by supersymmetry, for which the χ particle can be viewed as the bino and the mediator as the righthanded sbottom quark within the MSSM. In this case, the coupling is fixed by the bottom hypercharge. We note that for large λ χ > ∼ O(1), additional annihilation diagrams forqq † → bb as well asqq → bb contribute, which are modified by bound-state formation. In this work, we are mainly interested in the case of small λ χ , and therefore do not take these contributions into account, since their cross section scales as λ 4 χ and is subleading compared to the QCD contributions toqq † annihilation. The red lines in Fig. 7 correspond to the case with perturbative leading-order annihilation, and the blue lines correspond to our fiducial approximation that includes Sommerfeld enhancement and bound states up to n = 15. It is apparent that the blue contours allow for significantly larger masses m χ for a given λ χ . For example, for λ χ = 0.5 and ∆m = 20 GeV, the mass for which the relic density matches the observed value shifts from m χ 1.2 TeV to 2 TeV when including the aforementioned corrections. For the MSSM value, λ χ = 0.169, the contour almost coincides with the boundary, and the mass shifts from m χ 0.9 TeV to 1.8 TeV (for ∆m = 20 GeV). In addition, for a very small mass splitting, including bound states allows for mediator masses in the multi-TeV regime, around m χ = 3 TeV for ∆m = 5 GeV. This shift can be expected to be of major relevance for experimental searches for colored t-channel mediators within the coannihilation regime. It re-opens part of the parameter space that is constraint by conventional dark matter searches.

C. Conversion-driven regime and collider limits
In Fig. 8 we show the viable parameter space within the regime of conversion-driven freeze-out. The value of the coupling that is required to obtain the measured dark matter abundance is of order 10 −6 − 10 −7 in that case. We show several contours for λ χ /10 −7 = 2, 3, 5, 7. The smallness of the coupling implies that this production mechanism is compatible with null results from direct and indirect dark matter detection experiments, while still providing an explanation of the abundance of dark matter that is insensitive to the initial conditions.
The decay length cτ of the mediator, where τ is its lifetime, is shown by the gray contour lines in Fig. 8. It is of the order of a few centimeters to 1 m within most of the parameter space, going down to 1 mm close to the boundary. For the freeze-out computation, we limit ourselves to the parameter space where ∆m > m b , such that the two-body decayq → χb is kinematically allowed. For even smaller mass splitting, conversions proceed via scatterings, and the mediator would be stable on detector timescales.
The primary signal of conversion-driven dark matter production with a colored mediator are searches for heavy, (meta-)stable colored particles at the LHC. For ∆m < m b , the colored mediator becomes detector stable as its decay is four-body suppressed. We can directly apply the limit from the 13 TeV ATLAS search [48] derived for an R-hadron containing a b-squark. It excludes masses below 1250 GeV. The resulting limit is shown in Fig. 8 as a solid blue curve (and blue shaded exclusion region). For larger ∆m the decay length is in the range 1 mm ∼ 1 m such that a sizeable fraction of decays take place inside the inner detector. To estimate the reach of the same search for this case, we employ the reported cross section upper limits for the muon-system-agnostic analysis for a b-squark R-hadron. We rescale them by the relative suppression of the cross section upper limits toward small lifetimes reported in the similar ATLAS analysis [49] where the case of a gluino R-hadron has been considered. Note that this introduces a certain level of approximation. A recasting of the search is, however, beyond the scope of this work. We use the cross-section predictions from [50]. The resulting limit is displayed as the blue, dashed curve in Fig. 8. Furthermore, we display the limit from the recasting of the CMS 13 TeV R-hadron search [51] performed in [15] as the blue, dotdashed curve.
Being only sensitive to the fraction of R-hadrons traversing a significant part of the detector, the sensitivity of these searches is exponentially suppressed for small lifetimes. Dedicated analyses exploiting the displaced nature of the decay are, hence, expected to greatly improve the sensitivity to this scenario. While several such analyses have been performed by the collaborations, their target model differs considerably from the one considered here, significantly reducing their reach or raising questions about their applicability as pointed out in [52] (contribution 7). For instance, the sensitivity of the displaced jets search [53] considerably suffers from the imposed cut on the invariant mass of the displaced tracks. While the respective choice was optimized for the scenario considered in the search, it reduces the signal of the one considered here by around two orders of magnitude [52]. This is due to its relatively small mass splittings ∆m of order tens of GeV in our scenario, resulting in softer tracks. The search has been targeted to mass splittings of the order of hundreds of GeV. Another example of a potentially sensitive search is the one for disappearing tracks. The existing searches are targeted to charginos whose long lifetime arises due to a tiny mass splitting, O(100 MeV), to the dark matter particle. Accordingly, in the decay, an ultra-soft pion is emitted facilitating the use of a disappearance condition. In our scenario, the emitted b-jet is considerably harder than in the targeted model. However, the search is estimated to still provide sensitivity to the model considered here, as shown in the approximate recasting of [54] performed in [52]. In this recasting, the probability of the R-hadron to cause a charged track was also taken into account. We overlay the respective limit as the purple dotted curve in Fig. 8.
We conclude that, after including the impact of bound states, a wide part of the parameter space for conversiondriven freeze-out is still viable, and provides a clear target for long-lived particle searches at future LHC runs.

VI. CONCLUSION
In this work, we revisited the computation of the relic density in the presence of bound-state effects during dark matter freeze-out. With respect to previous work, we improved the calculations in various aspects and demonstrated the respective phenomenological implications on the cosmologically viable parameter space in the coannihilation and conversion-driven freeze-out scenario.
In the first part of this work, we reformulated the Boltzmann equations including arbitrary excitations of bound states and derived a general framework for incorporating their effects in terms of an effective annihilation cross section. While a full treatment of these effects requires the knowledge of all involved bound-state formation, decay, and transition rates, we introduced meaningful limiting cases when assuming fully efficient or nonefficient transitions. We provided simple analytical expressions for the effective cross section in these limits, as well as a general result. Furthermore, we showed that for an arbitrary set of bound states in ionization equilibrium, the effective cross section is independent of bound-state formation and transition rates, and only depends on a weighted sum of bound-state decay rates.
For the case of a colored coannihilator, we computed the radiative bound-state formation rates for arbitrary excitations with quantum numbers n, , and estimate the lowest order transition rates. Furthermore, we investigated the impact of NLO corrections to bound-state decays. We further discuss the relevance of NLO effects on bound-state formation and decay in App. B.
We then solved the coupled Boltzmann equation for the mediator and the dark matter particle in a t-channel model and assessed the impact of bound states for coannihilations as well as conversion-driven freeze-out. On the one hand, in ionization equilibrium, the effective mediator annihilation cross section is insensitive to the boundstate formation but directly proportional to the boundstate decay rates. Including excited states increases the effective cross section by about 20% in that case. On the other hand, after the breakdown of ionization equilibrium of the ground state, higher excitations become increasingly important. At the same time, a large boundstate formation rate extends the duration of ionization equilibrium down to smaller temperatures. Nevertheless, we found that freeze-out significantly extends beyond the period of ionization equilibrium for small relative mass splittings between the mediator and dark matter, phenomenologically most relevant in the region of high masses, m χ > ∼ 2 TeV. In this region of parameter space, our fiducial approximation that neglects boundstate transitions is expected to underestimate the effects of excited bound states, motivating further studies. In addition, we demonstrated that NLO corrections to the bound-state formation rate itself play only a moderate role in the setup considered here.
Evaluating the cosmologically viable parameter space, we found that the region for which conversion-driven freeze-out is relevant extends significantly when including bound-state effects, ranging up to the multi-TeV region. In addition, our findings imply that significantly higher dark matter masses are viable also within the coannihilation region. This has immediate consequences for dark matter searches. For instance, considering a mass splitting of 20 GeV and a coupling of ∼ 0.169, as predicted in the MSSM, the dark matter mass that matches the relic density is shifted from around 900 GeV to 1.8 TeV by the inclusion of the discussed effects. On the other hand, when keeping the masses fixed at m χ = 900 GeV and ∆m = 20 GeV, the coupling would change from 0.169 to around 5 × 10 −7 as it lies in the conversion-driven freezeout regime.
Dark matter produced via conversion-driven freeze-out is compatible with (in)direct detection limits due to a very weak coupling but yields signatures of long-lived particles at the LHC. We discussed the applicability of existing searches for R-hadrons, disappearing tracks and displaced jets, which exclude masses below about 0.6 − 1.2 TeV. Because of the increase of the viable parameter space for conversion-driven freeze-out, extending into the multi-TeV region, the scenario provides great prospects for long-lived particle searches at future LHC runs.
The computations considered here can be improved in future work in several ways, regarding the description of transitions among bound states, the decay of excited states with angular momentum, as well as the inclusion of thermal corrections.
with the radial overlap integral where f n (ρ) = F n (r)/p 3/2 rel is the dimensionless radial wave function of the bound state.
To compute the radial integral we use an integral representation of the hypergeometric function that appears in the scattering wave function, Note that by substituting s → 1 − s one finds that F (ρ) is real, such that we can drop the complex conjugate in I R . In addition, we use the generating function of the Laguerre polynomials for the bound-state wave function, to write The ρ integration can be performed using the definition of the Γ function, and we obtain with some rational coefficients c r , that will be unimportant in the following. We can generate the required integral by differentiating + +3 with respect to b. Because of the selection rule, + + 3 ≥ 2 + 2 > 2 , such that the sum over r in the square bracket drops out, as announced. Using setting z ≡ ζ b n 1+t 1−t = −ib and using sin(aπ) = i sinh(ζ s π) yields 1 0 ds s −iζs (1 − s) +iζs (z + i(2s − 1)) + +4 = π sinh(πζ s ) (1 + z 2 ) e −2ζsarccot(z) , (A16) which allows us to evaluate the radial integral. Using finally gives the result, eq. (52), for the radial integral. section, eq. (23), in the no-transition limit, eq. (25), is shown in Fig. 9. In [27,30] it was pointed out that the correction to the bound-state formation cross section becomes very large for small enough x, corresponding to T > ∼ E B n . Nevertheless, for these temperatures, ionization equilibrium holds to a large extent. In ionization equilibrium, the effective cross section becomes insensitive to the bound-state formation cross section. Therefore, the effect of the NLO corrections considered in [30] on the effective cross section is almost negligible for small x (left part of the left panel in Fig. 9). For large x, on the other hand, the temperature is so small that the finite-temperature contribution of the NLO corrections gives a negligible contribution. In this region, the zerotemperature correction dominates. This is the reason why the difference between LO and the NLO correction considered in [30] is moderate in the right part of the left panel in Fig. 9. However, it becomes more relevant for excited states, due to the larger effective strong coupling, given our scale choice eq. (47).
In order to further assess the impact of NLO corrections, we show the dependence of the effective cross section when changing all scales at which the strong coupling is evaluated by a factor of two or a half, respectively, by the colored bands in Fig. 9. Within the perturbative uncertainty, both results are consistent with each other. We observe that including the NLO corrections considered in [30] leads only to a small reduction of the scale uncertainty (right panel of Fig. 9). This indicates that further sources of higher-order corrections, including those listed above, would have to be taken into account for a complete NLO analysis.
The effect of the NLO corrections on the boundary between the coannihilation and conversion-driven regime is shown in Fig. 10, and compared to the impact of taking excited states into account. We find that the latter is  [30] on the boundary between the coannihilation and conversion-driven regime (red: with BSF NLO correction, blue: without). The impact on the boundary is smaller than the difference that arises when including excited states (n ≤ 15) as opposed to the ground state only (n = 1), which is shown for comparison for both cases, respectively. significantly more important.

NLO corrections to bound-state decay
Real and virtual correction to the decay B 10 → gg have been computed in [41]. The relative correction at NLO in the limit of massless quarks is given by eq. (65) in the main text. Note that collinear singularities in the real correction cancel when including the virtual piece [41], analogously to heavy quarkonium decay [38][39][40].
In Fig. 11   scale. The central line corresponds to µ MS = mq while the lower and upper boundaries of the red shaded band correspond to the choices µ MS /mq = 1/2 and 2, respectively. We adopted µ MS = mq in the main text, while µ MS = 2mq is used e.g. in [41]. We observe that the NLO correction is significantly smaller for µ MS = mq, at the level of a few percent. This justifies using the LO decay rate in our main analysis for this scale choice. In Fig. 11 (right panel) we show the dependence of the decay rate on µ MS at LO and NLO, respectively. As expected, the NLO result is significantly less sensitive to the scale choice. Note that for these figures we have set n f = 5 and neglected the contribution from the top quark since the use of the massless approximation is in general not well justified in that case. Using the expressions for the real corrections for massive quarks obtained in [41] confirms that the top quark contribution would amount to a small change of the already small NLO correction. Note that in Fig. 11 we only vary α ann s while keeping α eff b fixed.
In Fig. 12 we show the impact on the boundary line between conversion-driven freeze-out and coannihilation when taking into account NLO corrections to the decay. As expected, their impact is very small both for the ground state only and when taking into account excitations. Note that to obtain the NLO line when taking excited states into account we have assumed that Γ NLO dec /Γ LO dec is identical for all states with arbitrary n and = 0.