Excited bound states and their role in dark matter production

We explore the impact of highly excited bound states on the evolution of number densities of new physics particles, specifically dark matter, in the early Universe. Focusing on dipole transitions within perturbative, unbroken gauge theories, we develop an efficient method for including around a million bound state formation and bound-to-bound transition processes. This enables us to examine partial-wave unitarity and accurately describe the freeze-out dynamics down to very low temperatures. In the non-Abelian case, we find that highly excited states can prevent the particles from freezing out, supporting a continuous depletion in the regime consistent with perturbativity and unitarity. We apply our formalism to a simplified dark matter model featuring a colored and electrically charged $t$-channel mediator. Our focus is on the regime of superWIMP production which is commonly characterized by a mediator freeze-out followed by its late decay into dark matter. In contrast, we find that excited states render mediator depletion efficient all the way until its decay, introducing a dependence of the dark matter density on the mediator lifetime as a novel feature. The impact of bound states on the viable dark matter mass can amount to an order of magnitude, relaxing constraints from Lyman-$\alpha$ observations.

Understanding the composition of matter in our Universe constitutes a major challenge of today's fundamental physics.Notably, explanations of both the observed dark matter density and matter-antimatter asymmetry necessitate the introduction of physics beyond the Standard Model (SM) and therewith the computation of interactions among new -and presumably heavy -particles in the early Universe.If such new particles interact via a light force carrier, a significant contribution to their depletion may be given by the formation and subsequent decay of bound states, which has intriguing consequences for their thermal history.For instance, in the context of electroweakly charged dark matter [1][2][3] and colored coannihilation scenarios [4][5][6][7][8][9][10][11][12][13][14][15], it has been shown that the inclusion of bound state effects can strongly alter the prediction for the relic density.
Generally, we can classify radiative bound state formation (BSF) processes in terms of its leading multipole contribution: of a charged scalar field can be extremely relevant [16,17].As the emission carries away charge, it changes the initial and final two-particle state, leading to non-orthogonal states and ultimately to a non-vanishing monopole contribution.The BSF cross section via monopole transitions has been worked out for arbitrary excited states.However, partial-wave unitarity can be problematic already for capture into the ground state [16].
3. Quadrupole: One example where quadrupole moments contribute at the leading order, considers the emission of a real scalar field in the radiative bound state formation process [30][31][32][33][34]. Beyond capture into the ground state, little is known about higher excited states and bound-to-bound transitions.
In this work, we will present a more detailed investigation of the second case, where BSF and transitions among bound states in unbroken gauge theories are dominated by the (chromo) electric dipole contribution.While it is indeed the most considered scenario, it still remains unclear by how much the inclusion of highly excited bound states into the chemical network contributes to the depletion of the dark matter relic density.As shown recently in Ref. [14], the impact of higher excitations can be sizeable, in particular, when considering scenarios beyond the paradigm of weakly interacting massive particles (WIMPs), such as conversion-driven freeze-out [35], and for bound states driven by a perturbative, unbroken non-Abelian gauge symmetry.The greatest obstacle is the accurate evaluation of the dipole matrix elements for BSF and for transitions among highly excited states, especially in non-Abelian theories.General formulas to evaluate the dipole matrix elements for all principle and angular momentum quantum numbers, n and ℓ, respectively, have been provided in [14].
Here, we significantly improve on efficiency and numerical stability of their evaluation, which allows us to explore the contribution of up to half a million bound states (all n ≤ 1000, ℓ ≤ n − 1 states) when considering BSF.For bound-to-bound matrix elements, we are able to include all electric dipole transitions up to n ≤ 100, ℓ ≤ n − 1, which are about one million transitions in total for the processes allowed by the selection rules.
These improvements allow us to address several key scientific questions.First, we consider the velocity dependence of the BSF cross section in vacuum and investigate partial wave unitary properties in Abelian and non-Abelian gauge theories.For the case of SU (N c ), the inclusion of a large number of excitations sheds light on the breakdown of our theoretical framework and the necessity of its unitarization.Next, we consider the interplay of BSF, ionization, bound-to-bound transitions and bound state decays in the thermal bath of the early Universe.Following the formalism of [14,36], we describe their effect on the bound state constituents' abundance via an effective thermally averaged cross section.Focusing on the perturbative coupling regime which turns out to be consistent with unitarity, we pose the important question whether the effective cross section grows slower or faster than the inverse temperature, implying freeze-out or a continuous depletion of the abundance, respectively.The latter case is found for non-Abelian interactions and yields important phenomenological implications.
We exemplify these implications in detail in the last part of our work, where we apply our numerical framework to the superWIMP scenario [37,38], considering a simplified model with a colored and electrically charged t-channel mediator [39,40].We showcase the effect of highly excited bound states and find that the combined effect of strong and elecromagnetic interactions conspire to reduce the relic density by more than an order of magnitude compared to the case when including Sommerfeldenhanced annihilations only.Thereby we improve on earlier results within this scenario considering the ground state only [40,41].We find that this has important implications on the viable parameter space, in particular to relief Lyman-α constraints.
The remainder of this paper is organized as follows.In Sec.II, we discuss BSF in vacuum and investigate the velocity-dependence of the BSF cross section in view of partial wave unitarity.In Sec.III, we study the scaling of the effective cross section with temperature and discuss the implications for freeze-out.The setup and relevant quantities for the t-channel model are reviewed in Sec.IV, and our results for the impact of highly excited states on the superWIMP mechanism for dark matter production are discussed in Sec.V. We conclude in Sec.VI.Appendix A details the evaluation of BSF and transition matrix elements and App.B provides expressions for the cross sections and rates used in our analysis.Finally, in App.C, we show implications for the cosmologically viable parameter space of dark QED including highly excited bound states.

A. Matrix elements
Our starting point for computing radiative BSF in gauge theories is potential non-relativistic effective field theory [42][43][44][45].In this framework, the interaction of two non-relativistic particles with SU (N c ) or U (1) gauge vector fields at the ultra-soft scale can effectively be described by a (chromo-)electric dipole operator, g r•E, with gauge coupling g, relative distance r, and electric field E. In the two-particle subspace of the non-relativistic particles, this operator leads to matrix elements of the form The squared absolute value of these matrix elements is directly related to our physical quantities of interest: bound-state formation cross sections and bound-tobound rates for electric dipole transitions.For concrete examples related to dark matter, see Refs.[20,21,24,46].In App.A, we develop an efficient way of evaluating the matrix elements for systems, where the initial state ψ i and the final state ψ f are the solutions of two-body Schrödinger equations with corresponding potentials of Coulomb type: The effective coupling strength α eff of the initial and final state incorporates the details of the underlying particle physics model.We will consider Abelian gauge theories where the effective couplings are equal, and non-Abelian gauge theories where they can be different.We denote the effective couplings by α eff b and α eff s when referring to bound and scattering states, respectively.In the following, we explore the contribution of highly excited bound states in concrete realisations.

B. Abelian case
As a first concrete model, we consider Quantum Electrodynamics (QED) in the non-relativistic regime of the Fermionic particles.The two-particle states of interest then consist of two oppositely charged particles forming gauge singlets.Standard Model examples are hydrogen recombination and positronium formation.In QED, the potential of both the initial and final state is attractive with identical strength.The corresponding BSF cross section, describing the electric dipole transition process of a scattering state into a bound state with quantum numbers n and ℓ, is where α is the fine-structure constant.The difference of the initial and final state energy is the positive quantity ∆E = p 2 2µ + E B nℓ , where p 2 = µ 2 v 2 with v being the relative velocity and µ the reduced mass.Here, the absolute value of the binding energy is given by E B nℓ = µα 2 2n 2 as in QED α = α eff i = α eff f .The BSF cross section as defined in Eq. ( 3) is averaged over initial and summed over final spin degrees of freedom, as well as summed over the magnetic quantum numbers of the bound state (see App. B for details).Since the electric dipole operator is spin conserving, the same equation that applies to the Fermionic case (QED) also applies to, e.g., a complex scalar field charged under a U (1) gauge symmetry [31].For simplicity, we will commonly refer to both cases as U (1) in the following, as we are mainly interested in physics beyond the SM and we would like to cover both the Fermionic and complex scalar case simultaneously in our discussion.We numerically evaluate the U (1) BSF cross section in Eq. ( 3), as detailed in App.A, for various n, ℓ.We present the result in Fig. 1 in a model independent way, i.e. we multiply the cross section by vµ 2 /α 3 to (i) show the remaining dependence on α/v and (ii) to highlight deviations of the velocity dependence from 1/v.Let us begin with the well-known case of capture into the ground state n = 1, ℓ = 0 (blue dotted line).For v ≪ α, the cross section for capture into the ground state scales as 1/v so that the shown product, v × (σv), approaches a constant (see e.g.Ref. [31], also regarding its magnitude relative to annihilation).Similarly, we find the same velocity scaling when fixing the final state ℓ and summing over all n ≤ 1000 (gray lines).Specifically, the ℓ = 0 gray line is larger by a factor 1.268 than the n = 1, ℓ = 0 blue dotted line (capture into the ground state).We note that this factor is smaller than the upper bound 1.6 derived analytically in Refs.[18,47].For ℓ = 0, we additionally checked that our result for each n ≤ 5 coincides with the analytic results available in Ref. [48] and up to n = 10 for all ℓ with those of [14].
Interestingly, when summing the BSF cross section both in n and ℓ (blue lines), the velocity scaling becomes stronger than 1/v.Here, we sum all n ≤ 10 (dot-dashed), 100 (dashed), 1000 (solid) and always ℓ ≤ n − 1 accordingly.We compare our summed result to the well-known Kramer's logarithm [49,50] (red line): The ratio of the Kramer's logarithm and our fully summed numerical result is shown in the bottom panel of Fig. 1.For v ≪ α this ratio is expected to approach unity when including a very large number of excited states.We confirm this trend within the range of our numerical limitations, n ≤ 1000, ℓ ≤ n − 1, which can also be seen as a non-trivial check of our code. 1The Kramer's logarithm 1 For a given α/v, the amount of excited states needed for a convergent sum can be estimated.To this end, we consider the sum that leads to the Kramer's logarithm: From the denominator, one can estimate that for a percentage has also been mentioned in earlier dark matter related works [48,[51][52][53].
Although the logarithm leads to a slope steeper than 1/v, partial-wave unitarity is not violated here as the sum over different ℓ automatically includes different initial state angular momenta.However, we have checked that each individual angular momentum contribution of the initial state does not violate partial-wave unitarity as it does scale as 1/v for v ≪ α.For instance, in this limit, the BSF cross section of the s → np processes summed over all 2 ≤ n ≤ 1000 is larger by a constant factor 3.8 than the s → 2p.Each angular momentum contribution therefore remains below the partial wave unitarity limit for all v, provided the coupling is sufficiently small.In the non-Abelian case, we will observe a qualitatively different behavior, that implies partial wave unitarity violation even for (in principle) arbitrarily small couplings.

C. Non-Abelian case
As our second example, we consider a non-Abelian SU (N c ) gauge theory, specifically SU (3).Interestingly, it provides a qualitatively different phenomenology from QED even though in both cases the leading BSF processes are based on dipole transitions.While in QED the initial and final state potentials are both attractive with the same strength, in Quantum Chromodynamics (QCD), the initial state potential of the adjoint pair is repulsive.As we shall point out and explore in the following, this feature is accompanied by partial-wave unitarity violation in QCD or in a general SU (N c ) for N c ≥ 2.
In particular, we consider pairs of non-relativistic Fermionic particles in the fundamental and antifundamental representation of SU (3), i.e. the twoparticle space is spanned by the direct sum of singlet and adjoint pair states: 3 ⊗ 3 = 1 ⊕ 8.A SM example is heavy quarkonium formation in SM QCD.The BSF cross section describing the chromo electric dipole transition process of an adjoint scattering state into a singlet bound state is given by where an average over color degrees of freedom is performed.Since the final state is necessarily attractive to support bound states and the gluon carries away an octet color charge, the initial scattering state is always in the repulsive adjoint representation.The same equation also accuracy, the maximum principle quantum number n needs to be roughly an order of magnitude larger than a given α/v.This is also what we observe for our summed numerical result in Fig. 1.
For instance, summing all bound state contributions up to our numerical limit n ≤ 1000, ℓ ≤ n − 1, provides a percentage-level accuracy for α/v < ∼ 100 only, while n ≤ 100, ℓ ≤ n − 1 would require α/v < ∼ 10 for the same accuracy.From the truncation of the sum for n ≤ 1, 10, 100, 1000 (blue), one can infer that higher excited states contribute more for decreasing velocity.The red line shows the analytic approximation Eq. ( 4) when summing over all n, ℓ, known as Kramer's logarithm.The gray lines show the contribution to the sum for fixed ℓ, all of which approach a constant value at large α/v, respecting partial-wave unitarity.
Lower: Ratio of Eq. ( 4) and the summed result.
holds for the more general SU (N c ) gauge group, as well as for a complex scalar field (see e.g.Ref. [24,46] for the case of non-fundamental representations).For simplicity, when referring to the case SU (N c ) in the following we mean either pairs of Fermions or complex scalars in the fundamental and anti-fundamental representation.
In the remainder of this section, we consider a constant coupling in the perturbative regime.We do this to investigate the unitarity violation independently of the effects induced by a running coupling.As long as the beta function is negative, running coupling effects can only enlarge the regime where partial wave unitarity is violated.We will return to the impact of running in the FIG.2: Upper: The adjoint-to-singlet bound-state formation cross section in Eq. ( 6) shown for SU (3), summed over all possible final bound-state quantum numbers n and ℓ ≤ n − 1.From the truncation of the sum for n ≤ 1, 10, 100, 1000 (blue), one can infer that higher excited states contribute more for decreasing velocity.Lower: Rescaled ratio of BSF cross section for fixed initial angular momentum ℓ ′ and the corresponding partial-wave unitarity limit.For a given α, unitarity is violated when the rescaled ratio is above 1/α 3 , for which two examples are shown by the horizontal lines.

sections below.
Adopting this setting, we evaluate the SU (3) BSF cross section in Eq. ( 6) for various n, ℓ and present the results in Fig. 2 in analogy to the previous case of U (1).In the upper panel, one can notice that the velocity dependence of the summed BSF cross section may approach a power law when including as many excited states as our numerical limit allows for (n ≤ 1000, ℓ ≤ n − 1) and considering only the region of α/v < ∼ 100 where n ≤ 1000 are expected to be sufficient to capture the full result accurately.The scaling is much stronger than the previous Kramer's logarithm in U (1).In fact, it is even stronger than 1/v 2 .Specifically, when fitting the summed cross section n,ℓ (σv) nℓ ∝ v −γ for v ≪ α with a power law, we obtain γ ≈ 4.0.
Such a velocity dependence raises concerns regarding the partial-wave unitarity.To investigate this issue, a scattering state with fixed initial angular momentum, denoted by ℓ ′ , needs to be considered.From Eq. ( 6) we separate out this contribution by splitting Eq. (A1) into two contributions and denote the corresponding cross section by (σv) ℓ ′ nℓ , where the superscript denotes the fixed initial scattering state angular momentum.For a given ℓ ′ , we then sum (σv) ℓ ′ nℓ over all possible final state quantum numbers n, ℓ, compatible with the selection rules We require that this quantity respects the ℓ ′ th-partial wave unitarity bound.Specifically, as done in Ref. [18,[53][54][55], we consider the partial-wave unitarity cross section for 2 → 2 total inelastic collisions given by [56] (σv) In the lower panel of Fig. 2, we show the introduced summed (σv) ℓ ′ divided by the corresponding ℓ ′ th-partial wave unitarity cross section from Eq. ( 8).This ratio is multiplied by α −3 , leaving only a dependence on α/v.In this way, the partial-wave unitarity limits for a given α correspond to horizontal lines with values equal to α −3 .We show ℓ ′ = 0 (s-wave) to ℓ ′ = 3 (f-wave) BSF cross section summed in n ≤ 1000 and one particular high initial angular momentum with ℓ ′ = 100.The maximum principle quantum number is taken as our numerical limit (n ≤ 1000), which is sufficient for the first four partialwave summed cross sections to converge for α/v < ∼ 100.From these results, one can infer that each summed (σv) ℓ ′ grows faster than 1/v, i.e. the unitarity limit will be exceeded at a finite velocity, respectively, for any ℓ ′ we can resolve.This is independent of the chosen coupling strength.We point out that among all partial waves the s-wave unitarity bound is always violated at the largest velocities, though all curves approach a common behavior for decreasing velocity.We explicitly checked for the s-wave case that partial wave unitarity is even violated when including only a single, large n bound state, implying that it is not the summation over all possible final states for a given ℓ ′ which is problematic.
Moreover, we observe the same situation for various SU (N c ).To make an even more general statement, let us consider (σv) ℓ ′ as a function of the ratio α eff s /α eff b .Now, SU (N c ) is a special case which lies in the region α eff s /α eff b = [−1/3, 0[, where the lower limit corresponds to N c = 2 and the upper to the large N c limit.The U (1) case corresponds to α eff s /α eff b = 1.Our numerical results suggest that in the range α eff s /α eff b < 1, (σv) ℓ ′ scales stronger than 1/v, while for α eff s /α eff b ≥ 1 no evidence of partial wave-unitarity violation is found within our numerical boundaries.
Note that the mechanism behind the unitarization of BSF for α eff s /α eff b < 1 is still an open problem.While we point out this problem here for BSF via dipole transitions, it is worth noting that a similar question has been recently raised for monopole transitions in Ref. [16] where partial-wave unitarity can be violated already for capture into the ground state level n = 1.
For non-numerical evidence of unitarity violation in non-Abelian gauge theories, a simple analytic expression would be warranted.We managed to get an approximate analytic result by taking two limits in Eq. ( 6): (i) α eff s → 0 and subsequently (ii) ζb → ∞ where ζb = α eff b /(nv).Taking these limits, we obtain the result for the s-wave case in The two limits are justified for relative velocities which fulfill the condition2 In this velocity regime, we compared our direct numerical evaluation of Eq. ( 6) to the analytical result in Eq. ( 9) for a variety of N c and n values and find very good agreement.The fact that the s-wave BSF cross section reaches a constant value for the above velocity regime is another non-trivial check of our numerical implementation also for very large n.
However, the velocity regime may be too restricted to analytically proof unitarity violation for contributions of a single n.Namely, while for SU (N c ) the s-wave BSF cross section approaches the unitarity limit for increasing n, the velocity regime where the analytic expression is valid becomes smaller and eventually -(very) close to the unitarity bound -the condition in Eq. ( 10) cannot be met.Nevertheless, if there exists a theory with α eff s = 0, then there is no lower bound on v and violation of swave unitarity can be shown with the above formula.Notice that α eff s = 0 corresponds to the large N c limit of SU (N c ), which is, however, not justified for all velocities for a finite N c .
In the following, we explore phenomenological consequences focusing on the regime compatible with perturbativity and partial wave unitarity bounds.

III. SUPER CRITICAL BEHAVIOR
The impact of a set of bound states on the freeze-out dynamics of some particle species, j, can under very general conditions be described by the Boltzmann equation where n j is the number density and H the Hubble expansion rate.The effective cross section, ⟨σv⟩ eff , includes all the effects of pair annihilation as well as scatteringbound [4] and bound-bound transitions [1,14,36].Here, we investigate whether the inclusion of an increasing number of excited states can lead to an effective cross section that grows sufficiently fast to maintain efficient depletion of the (comoving) particle number density and, hence, prevent the particle species (e.g.dark matter) from freezing out.We call this condition a super critical behavior.
To obtain the threshold for such a super critical behavior, let us consider a typical scenario where a particle species with mass m is initially in thermal equilibrium with a heat bath with temperature T and entropy density s.We assume s ∝ T3 , H ∝ T 2 , i.e. no (significant) change in the relativistic degrees of freedom of the bath.Introducing the yield as Y j ≡ n j /s and parametrizing time by x ≡ m/T in Eq. ( 11), one can estimate the yield evolution as a function of x as follows.For times where the yield Y j (x) starts to deviate significantly from its equilibrium value, Y j (x) ≫ Y eq j (x), also known as the time of chemical decoupling, x cd , one can neglect the impact of Y eq j (x) in the Boltzmann equation.This allows for an analytic solution for the yield evolution after chemical decoupling (see e.g.Ref. [57]), which up to constants, can be estimated to scale as The integral converges for x 0 → ∞ only if ⟨σv⟩ eff (x) grows slower than x while for ⟨σv⟩ eff ∝ x γ with γ ≥ 1 the integral diverges.Accordingly, the particle species only freezes out for γ < 1 (typical WIMP) while the particle continues to deplete for γ ≥ 1.The critical value γ = 1 leads to logarithmic depletion and sets the threshold for what we define a super critical behavior.Above this threshold, the evolution of the yield approaches the scaling Y j ∝ x 1−γ for x ≫ x cd .In this case, the effective annihilation rate Γ eff ≡ n j ⟨σv⟩ eff is dynamically driven to be proportional to the Hubble rate Γ eff ∝ H.
In the presence of bound states, the effective cross section introduced above can be written as [14,36] where the first term is the usual pair annihilation cross section, thermally averaged.In all cases considered in this work, it includes the Sommerfeld effect [58,59].
The second term contains the thermal average of the BSF cross sections, denoted as σv nℓ .The summation over all bound-state quantum numbers contains a dimensionless, temperature dependent quantity, which obeys 0 ≤ R nℓ ≤ 1. 3 Thus, the presence of bound states al-ways increases the value of the effective cross section and could eventually lead to a super critical behavior.Introducing a simpler index to label a specific combination of quantum numbers, i = (nℓ), R i can explicitly be written as [14,36] The last line defines the total width of a particular bound state.It consists of the ionization rate, the rate of decay (via annihilation of the bound state's constituents), and bound-to-bound transition rates, respectively.The latter contains bound state excitation and de-excitation rates.
In practice, we use the Milne relation (cf.App.B 1) to obtain the excitation rate from the de-excitation rate, and Γ nℓ ion from σv nℓ .Note that the inclusion of boundto-bound transition rates increases n,ℓ σv nℓ R nℓ [36].

A. Dark QED
We now investigate the behavior of the effective cross section in Eq. (13) for our first example of a concrete model.In particular, we consider dark matter as a Dirac Fermion charged under a U (1) gauge group, which has been studied e.g. in Refs.[18,19,21].We shall call it dark QED in the following, where dark photons set the thermal environment with temperature T .
Dark QED has only two parameters, which are the dark matter mass, m, and the dark fine structure constant, α.For our analysis, we consider Sommerfeld enhanced annihilation and s-wave spin-singlet bound state decay into two dark photons.Relevant expressions are listed in App.B 2. We briefly comment on the influence of spin-triplet states below.The electric dipole interaction allows for transitions among the excited states in dark QED.The de-excitation rate is given by Taking into account all leading processes, we show our results for the dark QED effective cross section in Fig. 3. Let us start by neglecting all bound state contributions, considering the case of Sommerfeld enhanced dark matter pair annihilation into two dark photons only (gray line).ferent spin are not directly coupled to each other.We thus leave the spin sum implicit, see App.B for details.
FIG. 3: Effective cross section for heavy Fermions charged under a U (1) for constant coupling α = 0.1, including the contribution from bound states for all n, ℓ up to n = 1, 10, and 100 (blue dotted, dashed, and solid, respectively) and excluding BSF (gray solid).The effective cross section includes Sommerfeld enhanced annihilation, BSF, ionization, and all possible bound-to-bound transitions arising from the electric dipole interaction, as well as spin-singlet s-wave bound state decay.The red long-dashed line displays the slope ∝ x 1 for comparison.
As is well known, the Sommerfeld effect in this case introduces a 1/v dependence of the annihilation cross section for v ≪ α, leading to ⟨σv⟩ eff ∝ x 1/2 for sufficiently low temperatures.Next, we add the contribution of the spin-singlet ground state (blue dotted).Similarly to the Sommerfeld effect, the cross section for capture into the ground state also scales as 1/v for v ≪ α (as seen in Fig. 1).For T much lower than the binding energy, this leads to σv 10 ∝ x 1/2 .In this regime, the spin-singlet decay rate is much faster than the ionization rate due to Boltzmann suppression and consequently R 10 → 1, resulting again in an overall ⟨σv⟩ eff ∝ x 1/2 scaling.Compared to the Sommerfeld enhanced pair annihilation only, the effective cross section is larger by a constant factor in the low temperature regime, as expected [18].
Let us finally add many excited spin-singlet states to the system.We include, according to the selection rules, all possible electric dipole transitions among them via Eq.( 17) to evaluate the effective cross section in Eq. ( 13).In Fig. 3, it can be seen that for n ≤ 10 and ℓ ≤ n − 1 (dashed blue line) the effective cross section increases strongly until around x ∼ 10 3 , although within this regime also even higher excited states become important as seen by the solid and dashed lines separating.One could therefore not deduce that dark QED does not exceed the critical power scaling γ = 1 (indicated by the dashed red line) when including only n ≤ 10.However, the result for n ≤ 100 and ℓ ≤ n − 1 (blue solid line), which includes about 5000 bound states and about 10 6 transitions among them, clearly shows that the effec-FIG.4: Effective cross section for heavy Fermion triplets under SU (3), assuming α(m) = 0.025.The critical scaling ⟨σv⟩ eff ∝ x 1 (red dashed) is exceeded when respecting running couplings (darker blue), i.e. no freeze out occurs, as opposed to using constant coupling strength (lighter blue).The effective cross section includes Sommerfeld enhanced annihilation, BSF and ionization via chromo electric dipole interactions, as well as spin-singlet s-wave bound state decay.Note that no bound-to-bound transitions occur in dark QCD in dipole approximation.The gray line shows the case without including bound states.
tive cross section does not continue a strongly increasing trend but rather converges to a smaller power scaling.When fitting a power law in the regime 104 < ∼ x < ∼ 10 5 , where n ≤ 100 is trustworthy, we get a scaling of about ⟨σv⟩ eff ∝ x 0.6 .Note that from the Kramer's logarithm in Sec.II, it is not clear that the temperature dependence actually follows a power law.
The no-transition limit follows closely, but lies slightly above, the n = 1 line as we explicitly checked.From this we conclude that transitions among the bound states are an important effect, which needs to be taken into account for predicting the relic abundance in dark QED precisely.
We verified, by varying the included n between the two shown cases, that this less steep scaling, ⟨σv⟩ eff ∝ x 0.6 , is already found for including only n ≤ 80, from which we deduce that the inclusion of even higher excited states would not change the scaling.The scaling power is trivially unaffected by the value of the dark matter mass, as the mass cancels out in the shown product ⟨σv⟩ eff × m 2 .Moreover, one can show analytically that the scaling power is even unaffected when changing the value of α, which we have confirmed numerically in a wide range of α.We also explicitly checked that the inclusion of spintriplet bound states leaves the scaling power of the effective cross section at low temperatures unaffected.This numerical observation can be understood since they only differ from the spin-singlet bound states by a smaller decay rate which suppresses their contribution at small x but does not alter the late time scaling of the effective cross section (although it does increase the overall magnitude at late times).
From all this, we conclude that dark QED does not reach a critical scaling of the effective cross section in the low temperature regime within the electric dipole approximation.In other words, dark QED indeed freezes out within the approximations made. 4We show the impact of excited states on the relic abundance within dark QED in App.C, refining earlier results on this subject [18,53].

B. Dark QCD
As our second example of a concrete model, we consider dark matter as a Dirac Fermion in the fundamental representation of a new SU (3) gauge group, see e.g.[21,23], which is often called dark QCD.In the following analysis of dark QCD, we include standard expressions for the Sommerfeld enhanced pair annihilation cross section and the decay rate of the color singlet swave states, listed in App.B 3 neglecting spin-triplet states.Further, for the octet-to-singlet BSF cross section in the chromo electric dipole approximation we use Eq.(6).No singlet-to-singlet transitions can be mediated via the chromo electric dipole operator.Therefore, the effective cross section in Eq. ( 13) reduces to the notransition limit Γ i→j trans = 0 (cf.App.B 3) when allowing only for the leading chromo electric dipole interactions.This simplification allows us to focus exclusively on swave bound states which are the only ones with a nonvanishing decay rate in our approximations.
While in dark QED the coupling is frozen, in dark QCD, we take into account the one-loop running effect induced by gluon self-interactions in all considered quantities, as detailed in App.B 3. Yet, m remains the only dimensionful scale in the theory, implying that ⟨σv⟩ eff ×m 2 is independent of the choice of m.
In Fig. 4, we show ⟨σv⟩ eff × m 2 for the specific choice α(m) ≡ 0.025 across the regime α(m/x) < ∼ 1.We checked that the BSF cross sections are compatible with partial wave unitarity bounds for the velocities that give a sizeable contribution to the thermal average, within the perturbative regime.The scaling power in the absence of bound states (Sommerfeld only) is, unsurprisingly, the same as in the dark QED case at low temperatures, i.e. ⟨σv⟩ eff ∝ x 1/2 .The inclusion of s-wave bound states, however, leads to a much steeper scaling of the effective cross section.Even when omitting the running of the coupling strength, i.e. α = 0.025 at all scales (light blue line, using n ≤ 1000), we get a scaling power of about x 0.9 for x > 105 , which is close to the critical line.When fully including one-loop running (darker blue lines), the scaling of the effective cross section for dark QCD becomes super critical with a scaling of around x 1.1 .
We also find super critical behavior for other choices than α(m) ≡ 0.025.From this we conclude that dark QCD does not freeze-out even in the perturbative regime.However, since the scaling exceeds the critical one only slightly, we find a moderate effect on the abundance.For instance, assuming the dark sector is in thermal equilibrium with the SM bath, the addition of excitations 2 ≤ n ≤ 1000 leads to a reduction of the dark matter abundance by around 50% at x = 10 9 for m ∼ 10 6 GeV and a slightly larger percentage for smaller masses.
In the following sections, we consider a model featuring dark matter and an accompanying particle that is charged under SM QED and QCD, and hence subject to bound state effects.The electric charge of that particle will allow for color singlet-to-singlet transitions, implying that the inclusion of ℓ > 0 states pushes the scaling of the effective cross section further inside the super critical regime (cf.Fig. 6).As a consequence, the corrections to the relic abundance in the perturbative regime will be much larger.

IV. COLORED T-CHANNEL MEDIATOR MODEL
We consider a singlet Majorana Fermion χ being the dark matter candidate, and a scalar mediator q with gauge quantum numbers identical to those of either an up-or down-type right-handed SM quark q R (we focus on the latter case in our numerical results for concreteness).The dark matter field χ interacts with the SM only via the Yukawa interaction while the mediator q has additional interactions with the SU (3) c and U (1) Y SM gauge fields given by the usual kinetic term with covariant derivatives.We assume a mass m q > m χ such that the mediator can decay into the dark matter particle rendering only the latter stable on cosmological time scales.This model belongs to the class of so-called t-channel mediator models, see e.g.[61], that are being actively considered in the context of LHC dark matter searches, see e.g.[62] for a recent account on the subject. 5 Dark matter production is governed by the interaction of χ with the mediator field q, controlled by λ χ , as well as the dynamics of the mediator itself, largely driven by its gauge interactions.In particular, as the q and q † particles are color and electrically charged, they can form bound states, that have an important impact on the freeze-out [1, 9-11, 14, 15, 40, 41].Here, we are particularly interested in including excited bound states, following [9,14].
All of them can be described by the following set of coupled Boltzmann equations for the yields Y j , where x = m q /T , s is the entropy density, H the Hubble rate, and all of which depend on x.Here g j denotes the number of internal degrees of freedom, with g q ≡ 2N c = 6 denoting the sum of q and q † densities, g χ = 2, and g B nℓ = 2ℓ + 1 for bound states with angular momentum ℓ, capturing the degenerate magnetic quantum number.Note that the factor of 1/2 in Eq. ( 19) is due to our convention of including both the q and q † density in Y q .Eqs. ( 19) and ( 20) contain the following collision terms: 1.The effective cross section σ q q † v eff includes direct annihilation (including Sommerfeld enhancement following [35,69]) as well as the impact of bound states, as given by Eq. ( 13).We discuss the relevant BSF, transition and decay processes within the simplified model below, following and extending [14].
2. The rate Γ q→χ conv describes the conversion rate of q into χ particles.It is controlled by the Yukawa coupling, Γ q→χ conv ∝ λ 2 χ , and its size determines whether the freeze-out happens in the coannihilation, conversion-driven or superWIMP regime (see below).At high temperatures, it is dominated by scatterings X q → Y χ with appropriate SM particles X, Y , while at low temperatures the decay process q → qχ dominates.Accordingly, in the low temperature limit -relevant for the superWIMP mechanism considered below -it reads where is the vacuum decay rate of a single mediator particle in the limit m q → 0 (not to be confused with the bound state decay rates, that are dominated by the strong interaction).
In addition, the Boltzmann equations could be complemented by collision terms for the conversion process q q † → χχ, which are, however, negligible within the conversion-driven and superWIMP regimes as they are, again, proportional to λ 4 χ , and irrelevant in the coannihilation limit, and therefore not displayed here.
Let us now discuss in more detail the various possible regimes for dark matter genesis.As mentioned above, which regime is realized depends on the size of the conversion rate.More precisely, the most relevant quantity is the conversion rate for χ into q particles, where the last expression is the low temperature limit.Dark matter genesis is qualitatively different depending on the size of this rate relative to the Hubble rate for temperatures around the mediator mass (during the time when the mediator starts to chemically decouple from the SM bath), (a) In the coannihilation regime the q and χ populations are in mutual chemical equilibrium.The actual size of the conversion rate is irrelevant as long as it is strong enough to maintain chemical equilibrium [63,64].Within the coannihilation regime, the dark matter abundance is determined by the cross sections σ q q † v , σ χχ v and σ χq v (and in addition, as for all cases, the bound state dynamics), and generically here λ χ ∼ O(1) [9,14,15,61].
(b) In the conversion-driven case, the freeze-out of chemical equilibrium among χ and q drives the dynamics and the size of the conversion rate Γ q→χ conv largely influences the dark matter abundance [35,65].In addition, the efficiency by which the q (and q † ) abundance is depleted is relevant [35], controlled by σ q q † v and bound state effects [14].The conversion-driven case occurs for small couplings, typically λ χ ∼ O(10 −6 ), for which the χχ and χq terms in Eqs. ( 19) and ( 20) can be safely neglected.
(c) Finally, in the superWIMP scenario [37,38], the mediator has an even smaller decay rate, and can usually be considered as stable while the freeze-out of q q † annihilation occurs.The population of remaining q and q † particles then decays into χ at a temperature T for which H(T ) ∼ Γ χ→q conv , thereby generating the dark matter abundance [39].Technically, this means that the terms corresponding to inverse decays in Eqs. ( 19) and ( 20) can be neglected, in addition to those for χχ and χq annihilation, while the size of σ q q † v and the bound state dynamics are most important [39][40][41].To the extent that decay and freeze-out occur on different time-scales, the dark matter abundance is also insensitive to the size of the conversion (or equivalently decay) rate in that limit, since eventually each q (and q † ) produces one dark matter particle.Note that in addition to the superWIMP contribution a contribution from freeze-in [66][67][68] has to be considered which stems from inefficient decays (or scatterings) around x ∼ 1, i.e. when the mediator is in thermal equilibrium with the SM bath.The relative importance of superWIMP versus freeze-in contributions depend on the couplings and masses, see e.g.[39,40].However, in our analysis we are particularly interested in regions with a dominant superWIMP contribution.
In this work, we re-evaluate the superWIMP regime when including excited bound state effects.In particular, since the mediator is relatively long lived within this regime, its abundance crucially depends on how much of the mediator is depleted due to bound state dynamics.

B. Bound state rates and processes
The impact of bound states on Y q is captured by the effective cross section defined in Eq. ( 13) entering in the Boltzmann equation, Eq. ( 19).It depends on the set of bound states that are included as well as their formation, transition and decay rates, discussed in the following.
In the considered model, the scalar mediator particle q interacts both via electromagnetic interaction with bottom-like charge Q = −1/3 and strong interactions in the fundamental representation.This leads to differences compared to the case of pure Abelian or non-Abelian interactions discussed in Sec.III.In particular, the potentials determining the bound state spectrum and wavefunctions as well as BSF and decay are driven by QCD, while QED is relevant for transitions among the various energy levels [14].In the following, we briefly review the bound state processes included in our analysis, and then comment on the relevance of further extensions.
Bound state formation is dominated by the chromo electric dipole transition, (q q † ) [8] → B [1] nℓ + g , (25) going from an octet scattering state to a singlet bound state, and emitting an (ultrasoft) gluon.The effective interaction potentials V s(b) = −α eff s(b) /r for the scattering (s) and bound states (b) are with running strong coupling 6 evaluated at MS-scale µ s = m q v rel /2 for the scattering state, and at the Bohr momentum scale µ b = m q α eff b /2/n for the bound state B nℓ .Note that the latter definition is implicit, but can be easily solved either iteratively or numerically for each level n.The BSF cross section is then given by Eq. ( 6), with the effective running coupling entering the respective initial and final state wave-functions.Furthermore, we evaluate the coupling in the prefactor of Eq. ( 3), associated to the gluon emission, at the ultrasoft scale of the gluon energy µ BSF = m q /4 v 2 + (α eff b /n) 2 .Bound state formation is also possible via an electromagnetic dipole transition, but is negligible due to the smaller interaction strength [14] (see also Fig. 5 below).Computational details can be found in App.B 3.
Transitions among bound states are not possible via a single insertion of the QCD dipole interaction due to color conservation.Therefore, we consider electromagnetic dipole interactions as the leading bound-to-bound transitions: with rates computed as detailed in App.B 4. We use the wave-functions evaluated with the corresponding effective QCD couplings at their respective Bohr momentum scales for level n and n ′ and the electromagnetic fine structure constant α EM = 1/128.9 in the coupling prefactor of Eq. ( 17).
For bound state decay we include the process B n,ℓ=0 → gg, see App.B 3. The next-to-leading order correction to this decay channel within QCD has been shown to have only a minor effect in [14], and we therefore omit it here.
Let us briefly comment on possible further transition and decay processes.By simple power counting, electric quadrupole and magnetic dipole transitions are sup- 6 In all numerical computations involving running αs we used RunDec 3 [70] to evaluate the SM strong coupling (employing 5-loop running).
pressed relative to electric dipole transitions.Nevertheless, they could potentially have an impact by allowing new transition channels due to the modified selection rules.We checked (up to n ≤ 6) that electric quadrupole transitions have a negligible impact on the effective cross section.Furthermore, the decay of ℓ = 1 bound states is suppressed by at least a factor α 3 s relative to its de-excitation rate when assuming a power counting α EM ∼ O(α 2 s ).This is due to the higher derivatives of the radial wavefunction entering for higher ℓ, as well as the fact that decays into two gluons vanish for ℓ = 1 states at tree-level, leading to additional suppression due to a three-gluon decay or two-gluon decay at one-loop [71].Lastly, two-gluon transitions in SU (N c ) are expected to be suppressed by phase space factors and the repulsive potential of the necessarily adjoint intermediate state.Nevertheless, it would be interesting to investigate their impact in future work.

V. RESULTS FOR SUPERWIMP SCENARIO
We now discuss our results and phenomenological implications considering superWIMP production within the model introduced in Sec.IV.In the superWIMP scenario, the standard assumption has been that the freeze-out and decay of q are well separated in time, such that the late time χ abundance depends on the mediator freeze-out abundance only.The resulting dark matter density has thus been considered to be independent of the mediator lifetime τ q = 1/Γ q→qχ .In our model, we will show that the superWIMP mechanism proceeds in a qualitatively different way as the consequence of the super critical behavior of the effective cross section in non-Abelian gauge theories found in Sec.III B.
We assume that the mediator decays well before the QCD phase transition, such that the entire freeze-out dynamics takes place within the unconfined phase and involves α s in the perturbative regime (T > 1GeV) only.Numerical results presented throughout this section are shown for the representative benchmark point m q = 4 × 10 6 GeV.We discuss the mass dependence of our findings in Sec.V C.

A. Effective cross section
The abundance of the mediator q is governed by the effective cross section Eq. ( 13), that encapsulates the impact of bound states.In Fig. 5, we show the rates that enter this quantity for an exemplary subset with n ≤ 4 and for all ℓ ≤ n − 1, for a wide range of x = m q /T corresponding to m q /10 ≥ T ≥ 4 GeV.Since the mediator is charged under both SU (3) c and U (1) Y it combines features of the Abelian and non-Abelian cases discussed in Sec.III.
As expected, ionization (and correspondingly BSF) is cleary dominated by the QCD-mediated process, as can be seen by comparing the long-dashed and dotted lines in Fig. 5. Since q q † bound states exist only for the attractive color singlet configuration, color conservation dictates that bound-to-bound transitions are only contributing via an electric dipole interaction mediated by QED.The total transition rate from a given bound state (n, ℓ) into any higher or lower state is shown by the blue lines in Fig. 5.For the ground state (1, 0) only excitation oc-curs such that the rate becomes exponentially Boltzmann suppressed once the temperature drops below the corresponding difference of binding energies.The same is true for (2, 0) due to the selection rule ∆ℓ = ±1 for dipole transitions.For all other levels, the total transition rate approaches a finite value for low temperature (i.e.large x), corresponding to the rate of de-excitation into lower levels.In addition, we include the direct decay of ℓ = 0 FIG.6: Effective cross section for the colored and electrically charged mediator q.It includes the contribution from bound state levels (n, ℓ) for all ℓ and up to n ≤ 1, 10, 100 (blue lines) accounting for all dipole transitions among them.Also shown is the no-transition limit for n ≤ 1000 (green).Both grow more steeply than ∝ x (indicated by a thin black line) when sufficiently high n are included, implying mediator depletion without freeze-out.Here mq = 4 × 10 6 GeV.
states into a pair of gluons, which is the analog of mediator pair annihilation.This rate is practically temperature independent for x ≫ 1.Overall, the total width of any given level is dominated by the QCD-mediated ionization rate at temperatures T above or around the binding energy.For much lower temperatures, decay dominates for ℓ = 0, and QED-mediated transitions (to lower levels) for ℓ ≥ 1.
The resulting effective cross section, which combines QCD mediated BSF and QED transitions among bound states, is shown in Fig. 6, where we take excited states B nℓ with 0 ≤ ℓ ≤ n − 1 and up to a given maximal n into account.For the various blue lines, we include all possible electric dipole transitions.For any given temperature, the effective cross section converges when including a sufficient number of excited states.The reason for this convergence is that each excited state contributes in a limited temperature range only.While the underlying velocity dependence of the BSF cross section becomes increasingly complex at large n, a given bound state n, ℓ only starts to contribute significantly once the temperature drops down to roughly its respective binding energy, T ∼ E B nℓ .In fact, this is important for the validity of the dipole approximation, as it ensures that the temperature is well below the typical momenta of bound states that contribute significantly i.e. below their respective Bohr momenta.For T ≪ E B nℓ in contrast, its contribution is negligible due to the repulsive potential in the scattering state [9,14].
Introducing x n ≡ m q /E B nℓ , we find x 1 ≃ 7 × 10 2 , x 10 ≃ 5 × 10 4 , x 100 ≃ 3 × 10 6 for our benchmark m q = 4 × 10 6 GeV, which correspond to the x values at which the respective excited levels are expected to start contributing significantly to the effective cross section.Overall, the lower the temperature (i.e. the higher the x), the more relevant higher excited levels become for achieving a converged effective cross section.We include states up to n = 100, taking all transitions among them into account (blue solid line).We checked that this suffices to reach converged results within the perturbative regime.In particular, the difference between n ≤ 50 and n ≤ 100 is less than 0.2% for T > 1 GeV.
As visible in Fig. 6, the effective cross section (blue solid line) clearly shows a super critical behavior, i.e. the power scaling is significantly larger than ∝ x (red dashed).We stress that the interplay of bound states formed by the non-Abelian QCD interaction with transitions mediated by QED leads to a significant enhancement of the effective cross section compared to the limit of inefficient transitions, see green line in Fig. 6.The latter shows the result when omitting transition processes.It is similar to the case of dark QCD discussed in Sec.III B. This shows that excited states play an even more prominent role for a mediator charged under both QCD and QED considered here.Nevertheless, the effective cross section increases more steeply than ∝ x even in the no-transition approximation due to running.All the same, when including transitions but neglecting running, the slope is still steeper than ∝ x, i.e. the presence of bound-to-bound transitions causes a super-critical behavior in our model even without running coupling effects.

B. Relic abundance
The evolution of the yields Y q (x) and Y χ (x) for the mediator and the dark matter particle as obtained from solving the coupled Boltzmann equations (19) and (20) are shown in Fig. 7. Let us first discuss the case when including either only mediator annihilation (blue dotted line) or in addition the ground state (blue dot-dashed line, see e.g.[40,41]).In these cases, the yield Y q (x) freezes out at x ∼ O(10 2 −10 3 ).After freeze-out it remains constant until the age of the Universe becomes comparable to the mediator lifetime, i.e. for H ∼ Γ q→qχ .At this point, the mediator particles decay into dark matter, such that the final yield Y χ is identical to the freeze-out value of Y q which is set at much earlier times already.Accordingly, in a wide range of lifetimes the dark matter abundance does not depend on the dark matter coupling.This qualitative picture of the superWIMP mechanism has widely been adopted throughout the literature in the past.
Intriguingly, when including excited bound states, the super critical effective cross section leads to a continuous depletion such that Y q never freezes out.This can be seen in the solid line in Fig. 7, where bound states up to n ≤ 100 are included.The depletion of the number density is dominated by the effective cross section, until the decay of the mediator, q → qχ, becomes efficient (here x > ∼ 10 6 ).In fact, the effective annihilation rate FIG.7: Upper panel: Abundance evolution of χ (red) and q (blue) for mq = 4 × 10 6 GeV and Γq→qχ = 10 −17 GeV.When including no (SE only) or one (n = 1) bound state, the mediator yield Yq freezes out and then subsequently transfers its abundance to Yχ via q → qχ.Taking excited states and transitions among them into account leads to a continuous decrease of Yq (dashed and solid blue line for n ≤ 10 and n ≤ 100, respectively) that is only terminated by mediator decay once Γq→qχ > ∼ H. Lower panel: Ratio of χ yield when including bound states up to n ≤ 1, 10, 100 over the result without bound state contributions.FIG.8: Evolution of mediator and dark matter abundances when including bound states and transitions among them up to n = 100 (solid lines) or no bound states (dotted) for three different decay rates Γq→qχ = 2.5 × 10 −14 , 5 × 10 −16 , 10 −17 GeV.The long dashed line shows the limit Γq→qχ = 0, which decreases due to the large contribution of excited states to the effective cross section.The final value of Yχ therefore depends on Γq→qχ even in the superWIMP regime where the χ abundance generated via freeze-in at low x is negligible.
Γ eff = n q ⟨σv⟩ eff is kept on the edge of being efficient, i.e.Γ eff ∼ H, over the entire period of bound state induced depletion.The dynamics is qualitatively different to the standard picture as there is no temperature regime where the yield of the mediator has frozen out.
The quantitative impact of bound states up to n ≤ 1, 10, 100 on the dark matter abundance is explicitly shown in lower panel normalized by the result including only Sommerfeld enhancement and no bound states.Taking into account the ground state only yields a reduction by a factor ∼ 2 in the final abundance.When considering the full bound state effects (n ≤ 100), we find that the dark matter relic abundance is lowered by more than an order of magnitude.
Note that the mediator decays before the QCD transition such that dark matter production takes place in the unconfined phase involving α s in the perturbative regime.We explicitly verified that our results are insensitive to the behavior of the strong coupling at scales below 1 GeV.To check this we implemented different numerical prescriptions for treating the strong coupling at these low scales, and find that the final abundance is highly insensitive as long as the mediator lifetime ensures a decay before the QCD transition.Furthermore, we checked that our results are not influenced by contributions for which partial-wave unitarity is questionable, as the corresponding velocities are not relevant for the thermally averaged effective cross section in the regime T > 1 GeV.
In Fig. 8, we show the abundance evolution for three different values of the mediator lifetime.(The additional long-dashed line depicts the result when excluding mediator decays, i.e.Γ q→qχ = 0.) For all decay rates, super-WIMP production provides the dominant contribution to the dark matter density.While for the largest decay rate, Γ q→qχ = 2.5 × 10 −14 GeV, freeze-in still contributes almost 10%, it is fully negligible for the smaller rates chosen.Due to the continuous depletion of Y q in the presence of excited states, the time of decay does, indeed, have an impact on the final dark matter abundance, as can be seen from the three different values of Y χ (red lines) for x → ∞.In contrast, when neglecting excited states (dotted lines), the final yield is identical for all three mediator decay rates.A similar behavior can be found when including only a small number of bound states.

C. Implications
In Fig. 9, we finally show the dark matter mass, m χ (left axis labeling), for which the final χ abundance (displayed using the right axis labeling) yields the observed dark matter relic density, Ω χ h 2 ≃ 0.12 [72], as a function of the mediator decay rate, Γ q→qχ .For Γ q→qχ > ∼ 10 −12 GeV, dark matter production is dominated by freeze-in, while for lower decay rates, which we focus on in the following, the superWIMP mechanism sets the abundance.As discussed above, within FIG. 9: Dark matter mass mχ and decay rate Γq→qχ of the colored t-channel mediator for which the final yield matches the observed dark matter density Ωχh 2 = 0.12, using mq = 4 × 10 6 GeV.The abundance is set by the su-perWIMP mechanism for Γq→qχ < ∼ 10 −13 GeV, while freezein dominates for larger decay rates.When taking only Sommerfeld-enhanced q q † annihilation into account (dotted line) the relic density is independent of Γq→qχ within the su-perWIMP regime, i.e. the dotted line is horizontal.The same is true when taking the ground state n = 1 into account (dotdashed line).When including excited states, the mediator depletion due to BSF, transitions and decay of bound states continues until eventually H ∼ Γq→qχ, such that the relic density does depend on the mediator lifetime, corresponding to the black dashed (n ≤ 10) and black solid (n ≤ 100) lines.For Γq→qχ > ∼ 10 −17 GeV the mediator decays before the QCD transition, i.e. within the perturbative regime.The shaded areas bracket the possible values for even lower decay rates, see text for details.The red shaded region is excluded at 95% C.L. by Lyman-α forest observations [40].
this regime, the mediator abundance has so far usually been assumed to freeze out.In that case, each remaining mediator particle subsequently produces one dark matter particle, such that the precise decay rate is irrelevant for superWIMP production.This is indeed the case in Fig. 9, when taking only mediator pair annihilation (dotted) or in addition the ground state (dot-dashed) into account.
When including excited bound states, the mediator abundance continues to deplete until the age of the Universe reaches the mediator lifetime.The remaining mediators then produce dark matter via q → qχ.This implies that the final dark matter abundance does, in fact, depend on Γ q→qχ and, hence, on the dark matter coupling.Therefore also the dark matter mass m χ for which Ω χ h 2 ≃ 0.12 does depend on it.This can be seen in the solid and dashed curves in Fig. 9 corresponding to taking into account all excited states (ℓ ≤ n − 1) with n ≤ 10 and 100, respectively.The smaller dark matter abundance implies a larger dark matter mass, as compared to the case without excited states.
For Γ q→qχ > ∼ 10 −17 GeV, the mediator decay occurs prior to the QCD transition, i.e. within the perturbative regime T > 1 GeV.Note that in this regime the inclusion of bound state up to n = 100 appears to be sufficient.While we find significant contributions from bound states with n > 10 (cf. the difference between the dashed and solid lines in Fig. 9) this contribution is dominated by bound states n < 50.We reiterate that n = 100 is sufficient to find a convergent effective cross section, hence even higher bound states are expected to not alter our results.
For illustration, in Fig. 9, we include decay rates down to around 10 −18 GeV for which a significant fraction of mediators have not decayed at T = 1 GeV.In this region, the gray shaded areas conservatively bracket the uncertainty in the effective annihilation rate arising from the impact of confinement, by assuming that all mediators that are still present at T = 1 GeV either vanish (upper boundary) or fully decay into dark matter particles (lower boundary).However, in this work, we focus on the perturbative regime, Γ q→qχ ≥ 10 −17 GeV, for which we find that the difference between the upper and lower boundary is less than 1%.
Due to the late decay and large mass difference between the mediator and dark matter, the dark matter momentum distribution for the considered scenario can be significantly harder than the one of cold dark matter.The resulting free-streaming effect impacts structure formation on small scales probed by Lyman-α forest observations.As shown in Ref. [40], this results in a lower bound on the dark matter mass, m χ keV > 3.8 × x decay 106.75 g * S (x decay ) that can easily reach into the GeV range.
In Fig. 9, we display the corresponding 95% C.L. exclusion as a red shaded area.Interestingly, excited bound state effects have a significant impact on the implications of this constraint.While with Sommerfeld effect (and n = 1 bound state) only, a decay rate of around 10 −15 (4 × 10 −16 ) GeV would be excluded, the inclusion of excited bound states reveal that the entire region with Γ q→qχ ≥ 10 −17 GeV remains unchallenged by the considered Lyman-α constraint.
So far, we have focused on the benchmark mass m q = 4 × 10 6 GeV.For smaller masses, the effective cross section becomes larger and the mediator abundance decreases.Therefore, for a given x decay , the dark matter mass that leads to Ωh 2 = 0.12 increases and the Lyman-α bound become less constraining.Accordingly, we find that current Lyman-α forest constraints do not challenge the superWIMP scenario for m q < 4 × 10 6 GeV, if we restrict ourselves to decays within the perturbative regime of couplings, i.e. decays that take place before the QCD phase transition, x decay < m q /1 GeV.However, as the so-defined maximal x decay becomes smaller with smaller m q , highly excited states become less relevant somewhat diminishing the large effect of bound states towards small m q in our setup.
For masses larger than 4 × 10 6 GeV, the cosmologically viable dark matter mass decreases and Lyman-α constraints become more restrictive.In particular, starting from masses somewhat above m q = 4×10 6 GeV, they impose an upper bound on x decay that is more restrictive than the above-mentioned perturbativity condition.As a consequence of the tightening Lyman-α constraints on the allowed range of mediator decay rates, higher excitations become less important also towards larger masses.Eventually, for m q > 4 × 10 8 GeV, the entire superWIMP region is excluded (assuming thermalization of the mediator).

VI. CONCLUSION
In this work, we studied the impact of excited bound states on dark matter production.We found that they can be highly relevant, focusing especially on setups where unbroken Abelian and/or non-Abelian gauge interactions are responsible for the bound state dynamics.Considering bound state formation/ionization and (de-)excitation processes described by dipole transitions, we developed an efficient way to compute their rates numerically.It allows us to take into account excitations up to a principal quantum number n = 100 (n = 1000) in the presence (absence) of transitions among them involving more than 5000 individual BSF and O(10 6 ) transitions rates.
With this numerical tool at our disposal, we investigated several theoretical and phenomenological questions.First, considering an Abelian gauge theory like dark QED, we compared the summed BSF cross section to the well-known Kramer's logarithm, confirming its approximate behavior towards small v that increases faster than ∝ 1/v.However, this behavior results from a summation over different initial state partial waves.We checked that each initial angular momentum contribution is compatible with partial-wave unitarity for all velocities, and at sufficiently weak coupling.In contrast, we found that the total BSF cross section in non-Abelian gauge theories generally does violate partial-wave unitarity bounds even for perturbatively small coupling.While a closer investigation of this feature is left for future work, we focused on the phenomenologically relevant regime of velocities in our subsequent results for which unitarity bounds are satisfied.
We exemplified our results for two particle states with constituents transforming as 3 and 3 under SU (3).Due to the repulsive potential of the adjoint scattering state, in radiative BSF via gluon emission, each individual bound state cross section only contributes around a characteristic velocity.This renders very high excitations to be the dominant contribution to the thermally averaged effective cross section at very low temperatures.Consequently, we investigated the important question whether the effective cross section increases slower or faster than x = m/T , in which case the particle's number density would freeze-out or would continue to deplete.
For dark QED, we found a scaling ⟨σv⟩ eff ∝ x 0.6 towards large x.For dark QCD, ⟨σv⟩ eff grows slightly slower (faster) than ∝ x, if we exclude (include) the effect of a running coupling (assuming a negative beta function).However, the dependence of this qualitative difference on the effects due to the running coupling is only present in (dark) QCD, i.e. in the absence of transitions among the states.When color and electric charge is combined, then singlet-singlet transitions are possible and significantly enhance the effective cross section towards small x where higher n dominate, and steepen the effective cross section beyond ∝ x regardless of running effects.
Finally, we studied such a case in detail by considering a scenario where an electrically and color charged mediator -which accompanies dark matter -is subject to the formation of bound states.Such a setup is realized in so-called t-channel mediator models.We focused on the very weak coupling regime for which the dark matter density can be generated through a late decaying mediator particle, i.e. by the superWIMP mechanism.The commonly adopted assumption within this paradigm is that mediator freeze-out and decays can be considered independently, thereby rendering the resulting dark matter density to be independent of the the mediator lifetime (and therewith the dark matter coupling it depends on).
Here we found that considering excited bound state effects, this picture has to be revised.Due to an interplay of QCD bound states and transitions mediated by QED, the resulting effective cross section grows significantly faster than ∝ x, i.e. features a strongly super critical scaling, thereby retaining a sizeable depletion of mediator particles throughout its entire presence before it decays.Therefore, the dark matter relic density does depend on the mediator lifetime.In fact, we found that -restricting ourselves to decays that happen well within the unconfined phase -the yield reduction due to bound state effects can amount to an order of magnitude with respect to the result including Sommerfeld enhanced annihilation only.This has important consequences for the cosmologically viable parameter space.For a given mediator mass, bound state effects shift the dark matter mass required to match Ωh 2 towards larger values by up to an order of magnitude, thereby weakening constraints from structure formation through Lyman-α observations.While we are choosing dark matter physics as an application of the considered bound state effects, we note that they could also play an important role in the context of baryogenesis.For instance, according to [73], the late decay of a color-charged scalar can generate the matterantimatter asymmetry.A significant enhancement of the effective cross section due to bound states could open up parts of the parameter space otherwise excluded by strong bounds from N eff .In addition, an analysis of phenomenological consequences within dark sector models, for which dark matter is charged under an unbroken dark gauge symmetry, would be interesting, also in view of phase transitions and associated gravitational wave signatures [23,25].On the theoretical side, our results motivate an investigation of unitarization of bound state formation processes mediated by unbroken non-Abelian gauge interactions within the regime of perturbatively small couplings.
with z = 4κκ ′ (κ+κ ′ ) 2 < 1 and κ = αµ/n, κ′ = α ′ µ ′ /n ′ and κ ̸ = κ′ .For the special case α ′ µ ′ = αµ (particularly important for QED), see Ref. [21,74].The normalization is given by: For large initial state principle quantum numbers n = O(10), the numerical evaluation of the Hypergeometric functions can be problematic for transitions where z → 1 (n ′ → n).To improve the numerical stability in this regime, we use a lengthy expression for the Hypergeometric function expanded around 1 − z.For the special arguments of the Hypergeometric functions as given above, we can further simplify the expression and finally obtain to all orders in 1 − z: Notice that this substitution is an identity for all sets of {a; b; c; z} arguments given above.In practice, we use the substitution for z > 0.7, which allows us to obtain stable numerical results for all bound-to-bound transitions with n ≤ 100.Increasing the number of digits for z (e.g. via MaxExtraPrecision in Mathematica) allows for even larger n values and also to check the stability, with the cost of loosing efficiency.
Appendix B: Cross sections and rates

Thermal average and Milne relations
In the non-relativistic limit, the thermal average of the bound state formation cross section can be written as B1) where f (∆E) = 1/(e ∆E/T − 1).We note that the thermal averages in non-Abelian gauge theories need to be evaluated carefully due to oscillatory features. 8 As usual, the cross section includes an average of initial state degrees of freedom and a sum over final state degrees of freedom.For the case of bound states with Fermionic spin-1/2 constituents, one has to distinguish between bound states in a spin-singlet or spin-triplet configuration.The main difference between the two cases occurs for the bound state decay rate, see below.Since the 8 For example, for effective coupling α eff s < 0 < α eff b relevant in SU (Nc), (σv) nℓ features n − ℓ − 1 local minima in its velocity dependence.Fewer local minima arise when 0 < α eff s < α eff b .
electric dipole interaction is spin-independent, boundto-bound transitions do not change the spin within our approximations, i.e. transitions occur exclusively among spin-singlet states or spin-triplet states, respectively.The transition rates are identical in both cases.Furthermore, for scattering-to-bound processes, the contribution to the BSF cross section for formation of any single spin degree of freedom of the bound state is the same.Therefore, we account for the formation of spin-singlet or spin-triplet bound states by the prefactor ξ = 1/4 or ξ = 3/4, respectively, in Eqs.(B6) and (B16).Furthermore, the BSF cross sections also apply for bound states composed of charged scalars, with ξ = 1.The reason is that the factor 1/4 = 1/2 × 1/2 from averaging over initial state spin and the factor 4 = 3 + 1 from summing over all final state spin configurations cancel out in the Fermionic case.
For computational efficiency, we use the following Milne relations based on detailed balance to obtain the inverse processes of BSF, i.e. ionization, and bound state de-excitation, i.e. excitation.Assuming leading order non-relativistic expansions for the equilibrium yields Y eq j (the combined particle and anti-particle yield) and Y eq Bi , the ionization rate can be expressed as Note that for bound states with Fermionic constituents, it holds for both spin-singlet and spin-triplet states when using the appropriate BSF cross section as described above and accounting for the factors of spin degrees of freedom contained in the equilibrium yields.The transition rates among the set of bound states are similarly related via Bi . (B3)

Dark QED
In the absence of light Fermions the coupling is not running within dark QED and, hence, the effective coupling strengths are identical for bound and scattering states, i.e. α eff s = α eff b = α here.The annihilation cross section into a pair of massless dark photons, including Sommerfeld enhancement, is given by [75,76] (σv where (B5) The BSF cross section, which includes the summation over the degenerate magnetic quantum number of the final state nℓ and also considers all possible initial angular momenta, is in the case of U (1) interactions given by (σv) nℓ = ξ π α 2 m 2 S BSF can be evaluated numerically in an efficient way by using the radial integral formulas laid out in App.A 1. Further analytic simplifications due to the absence of scale running of the coupling strength and the identical initial and final state gauge representations are possible, though not used in our numerical evaluation.
For spin-singlet bound states, we adopt the decay rate into two dark photons for the s-wave states from Ref. [77,78]: For spin-triplet bound states we adopt the decay rate into three dark photons for the s-wave states from Ref. [79]: Lastly, the de-excitation rate relates to the dipole matrix element via

Dark QCD
In Yang-Mills-theories, gauge Boson self interactions give rise to running of the coupling strength, even in the absence of light Fermions.For our numerical benchmark, we defined the value α(m) ≡ 0.025 and employ oneloop running.Using this choice, the non-perturbative regime α(m/x) = 1 starts at x ≈ 4 × 10 9 , which holds for any mass since m must drop out of dimensionless expressions, being the only mass scale in the theory.The heavy Fermions are from the fundamental representation of SU (3) and thus can form singlet (1) and octet (8) two particle states.We evaluate the wave functions of the initial scattering (bound) states, s (b), at their respective (Bohr-) momentum scale.The effective couplings are thus given by where C F = 4/3 and C A = 3. Eq. ( B11) is an implicit definition easily solved for, either numerically or analytically order by order.The annihilation into two gluons is possible from singlet or octet scattering states and takes the form We neglect spin-triplet bound states for dark QCD.Dark QCD automatically corresponds to the limiting case of no transitions, i.e.Γ i→j trans = 0, therefore the effective cross section Eq. ( 13) simplifies to [14,36] This result can be seen as a straightforward generalization of the single bound state case [4], extended to a sum over individual bound states that do not impact each other.

superWIMP scenario
The expressions for BSF of a colored t-channel scalar mediator is identical to that for dark QCD in Eq. (B16), now using m = m q and the SM strong coupling strength α s to 5-loop accuracy as implemented in RunDec 3 [70] in place of α, as well as ξ = 1.The annihilation cross section reads (σv) ann = 14 27 πα(2m q ) 2 m 2 q 2 7 S and the decay rate is given by [9,14] Γ nℓ→gg dec = δ ℓ,0 m q C F 8n 3 α(m q ) 2 α eff b 3 . (B20)

FIG. 1 :
FIG.1: Upper: Bound-state formation cross section Eq. (3) for U (1), summed over all possible final bound-state quantum numbers n and ℓ ≤ n − 1.From the truncation of the sum for n ≤ 1, 10, 100, 1000 (blue), one can infer that higher excited states contribute more for decreasing velocity.The red line shows the analytic approximation Eq. (4) when summing over all n, ℓ, known as Kramer's logarithm.The gray lines show the contribution to the sum for fixed ℓ, all of which approach a constant value at large α/v, respecting partial-wave unitarity.Lower: Ratio of Eq. (4) and the summed result.

6 FIG. 5 :
FIG.5: Ionization, decay and transition rates for all bound state levels (n, ℓ) up to n = 4 for mq = × 10 6 GeV.The gray dashed vertical line indicates the temperature that corresponds to the binding energy EB of the respective bound state.Note that the transition rates contain all possible excitations and de-excitations (here summed up to n ′ = 20).