#  Excited bound states and their role in dark matter production

Tobias Binder<EMAIL_ADDRESS>Physik Department T31, Technische
Universität München, James-Franck-Straße 1, D-85748 Garching, Germany  Mathias
Garny<EMAIL_ADDRESS>Physik Department T31, Technische Universität
München, James-Franck-Straße 1, D-85748 Garching, Germany  Jan Heisig
<EMAIL_ADDRESS>Institute for Theoretical Particle Physics and
Cosmology, RWTH Aachen University, D-52056 Aachen, Germany  Department of
Physics, University of Virginia, Charlottesville, Virginia 22904-4714, USA
Stefan Lederer<EMAIL_ADDRESS>Physik Department T31, Technische
Universität München, James-Franck-Straße 1, D-85748 Garching, Germany  Kai
Urban<EMAIL_ADDRESS>Physik Department T31, Technische Universität München,
James-Franck-Straße 1, D-85748 Garching, Germany

###### Abstract

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.

††preprint: TUM-HEP 1469/23††preprint: TTK-23-21

###### Contents

  1. I Introduction
  2. II Bound-state formation in vacuum
    1. II.1 Matrix elements
    2. II.2 Abelian case
    3. II.3 Non-Abelian case
  3. III Super critical behavior
    1. III.1 Dark QED
    2. III.2 Dark QCD
  4. IV Colored t-channel mediator model
    1. IV.1 Review of production mechanisms
    2. IV.2 Bound state rates and processes
  5. V Results for superWIMP scenario
    1. V.1 Effective cross section
    2. V.2 Relic abundance
    3. V.3 Implications
  6. VI Conclusion
  7. A Evaluation of large _n_ dipole transitions
    1. A.1 Scattering-to-bound
    2. A.2 Bound-to-bound
  8. B Cross sections and rates
    1. B.1 Thermal average and Milne relations
    2. B.2 Dark QED
    3. B.3 Dark QCD
    4. B.4 superWIMP scenario
  9. C Relic abundance for dark QED

##  I Introduction

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 Mitridate:2017izz ; Bottaro:2021snn ; Bottaro:2022one  and colored
coannihilation scenarios Ellis:2015vaa ; Ellis:2015vna ; Liew:2016hqo ;
Kim:2016zyy ; Harz:2017dlj ; Harz:2018csl ; Biondini:2018pwp ;
Biondini:2018ovz ; Fukuda:2018ufg ; Biondini:2019int ; Garny:2021qsr ;
Becker:2022iso , 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:

  1. 1.

_Monopole:_ Bound-state formation via emission of a _charged scalar field_ can
be extremely relevant Oncala:2019yvj ; Oncala:2021swy . 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 Oncala:2019yvj .

  2. 2.

_Dipole:_ Known examples where radiative BSF is dominated by the contribution
of the (dark electric) dipole moment, are $U(1)$ vonHarling:2014kha ;
Kamada:2019jch ; Binder:2020efn ; Biondini:2023zcz  or $SU(N_{c})$
Harz:2018csl ; DeLuca:2018mzn ; Asadi:2021pwo ; Binder:2021otw ; Asadi:2022vkc
; Becker:2023opo ; Assi:2023cfo ; Biondini:2023zcz  gauge symmetry extensions
of the SM. In these cases, the emitted particle is a massless gauge _vector
field_. Another possibility considers the accompanying dark matter particle to
have SM electroweak Asadi:2016ybp ; Mitridate:2017izz ; Bottaro:2021snn ;
Bottaro:2022one  and/or color charge Ellis:2015vaa ; Ellis:2015vna ;
Liew:2016hqo ; Kim:2016zyy ; Harz:2017dlj ; Harz:2018csl ; Biondini:2018pwp ;
Biondini:2018ovz ; Fukuda:2018ufg ; Biondini:2019int ; Garny:2021qsr ;
Becker:2022iso ; Biondini:2022ggt . A famous example of the latter is squark
coannihilation in the context of the Minimal Supersymmetric extension of the
Standard Model (MSSM), or simplified $t$-channel mediator models inspired by
it.

  3. 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 Wise:2014jva ; Petraki:2015hla ; An:2016kie ;
Biondini:2021ycj ; Biondini:2021ccr . 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.
Garny:2021qsr , 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
Garny:2017rxs , 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 $\ell$,
respectively, have been provided in Garny:2021qsr . 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\leq 1000,\;\ell\leq n-1$ states) when considering BSF. For bound-to-
bound matrix elements, we are able to include all electric dipole transitions
up to $n\leq 100,\;\ell\leq 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 Garny:2021qsr ; Binder:2021vfo , 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 Covi:1999ty ;
Feng:2003uy , considering a simplified model with a colored and electrically
charged $t$-channel mediator Garny:2018ali ; Decant:2021mhj . 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
Sommerfeld-enhanced annihilations only. Thereby we improve on earlier results
within this scenario considering the ground state only Decant:2021mhj ;
Bollig:2021psb . We find that this has important implications on the viable
parameter space, in particular to relief Lyman-$\alpha$ 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.

##  II Bound-state formation in vacuum

###  II.1 Matrix elements

Our starting point for computing radiative BSF in gauge theories is _potential
non-relativistic effective field theory_ Pineda:1997bj ; Beneke:1998jj ;
Beneke:1999zr ; Brambilla:1999xf . 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\,\mathbf{r}\cdot\mathbf{E}$, with gauge coupling $g$,
relative distance $\mathbf{r}$, and electric field $\mathbf{E}$. In the two-
particle subspace of the non-relativistic particles, this operator leads to
matrix elements of the form

$\displaystyle\bra{\psi_{f}}\mathbf{r}\ket{\psi_{i}}=\int\text{d}^{3}r\;\psi_{f}^{\star}(\mathbf{r})\;\mathbf{r}\;\psi_{i}(\mathbf{r}).$
(1)

The squared absolute value of these matrix elements is directly related to our
physical quantities of interest: bound-state formation cross sections and
bound-to-bound rates for electric dipole transitions. For concrete examples
related to dark matter, see Refs. Binder:2020efn ; Binder:2021otw ;
Biondini:2023zcz ; Biondini:2023yxt .

In App. A, we develop an efficient way of evaluating the matrix elements for
systems, where the initial state $\psi_{i}$ and the final state $\psi_{f}$ are
the solutions of two-body Schrödinger equations with corresponding potentials
of Coulomb type:

$\displaystyle V_{i,f}(r)=-\frac{\alpha_{i,f}^{\text{eff}}}{r}\;.$ (2)

The effective coupling strength $\alpha^{\text{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 $\alpha_{b}^{\text{eff}}$ and $\alpha_{s}^{\text{eff}}$
when referring to bound and scattering states, respectively. In the following,
we explore the contribution of highly excited bound states in concrete
realisations.

###  II.2 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 $\ell$, is

$\displaystyle(\sigma v)_{n\ell}=\frac{4\alpha}{3}\Delta
E^{3}|\bra{\psi_{n\ell}}\mathbf{r}\ket{\psi_{\mathbf{p}}}|^{2},$ (3)

where $\alpha$ is the fine-structure constant. The difference of the initial
and final state energy is the positive quantity $\Delta
E=\frac{\mathbf{p}^{2}}{2\mu}+E_{{\cal B}_{n\ell}}$, where
$\mathbf{p}^{2}=\mu^{2}v^{2}$ with $v$ being the relative velocity and $\mu$
the reduced mass. Here, the absolute value of the binding energy is given by
$E_{{\cal B}_{n\ell}}=\frac{\mu\alpha^{2}}{2n^{2}}$ as in QED
$\alpha=\alpha_{i}^{\text{eff}}=\alpha_{f}^{\text{eff}}$. 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 Petraki:2015hla . 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,\ell$. We present the result in Fig. 1 in a model
independent way, _i.e._ we multiply the cross section by $v\mu^{2}/\alpha^{3}$
to (i) show the remaining dependence on $\alpha/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$, $\ell=0$ (blue dotted
line). For $v\ll\alpha$, the cross section for capture into the ground state
scales as $1/v$ so that the shown product, $v\times(\sigma v)$, approaches a
constant (see _e.g._ Ref. Petraki:2015hla , also regarding its magnitude
relative to annihilation). Similarly, we find the same velocity scaling when
fixing the final state $\ell$ and summing over all $n\leq 1000$ (gray lines).
Specifically, the $\ell=0$ gray line is larger by a factor 1.268 than the
$n=1$, $\ell=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.
vonHarling:2014kha ; Bethe:1957ncq . For $\ell=0$, we additionally checked
that our result for each $n\leq 5$ coincides with the analytic results
available in Ref. Petraki:2016cnz  and up to $n=10$ for all $\ell$ with those
of Garny:2021qsr .

Figure 1: _Upper:_ Bound-state formation cross section Eq. (3) for $U(1)$,
summed over all possible final bound-state quantum numbers $n$ and $\ell\leq
n-1$. From the truncation of the sum for $n\leq 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,\,\ell$, known as Kramer‘s logarithm. The gray lines show the contribution
to the sum for fixed $\ell$, all of which approach a constant value at large
$\alpha/v$, respecting partial-wave unitarity. _Lower:_ Ratio of Eq. (4) and
the summed result.

Interestingly, when summing the BSF cross section both in $n$ and $\ell$ (blue
lines), the velocity scaling becomes stronger than $1/v$. Here, we sum all
$n\leq 10$ (dot-dashed), $100$ (dashed), $1000$ (solid) and always $\ell\leq
n-1$ accordingly. We compare our summed result to the well-known Kramer’s
logarithm Kramers1923XCIIIOT ; 1978JETP…48..639K  (red line):

$\displaystyle\sum_{n,\ell}(\sigma
v)_{n\ell}\simeq\frac{32\pi}{3\sqrt{3}}\frac{\alpha^{2}}{\mu^{2}}\frac{\alpha}{v}[\log(\alpha/v)+\gamma_{E}]\,\text{,
for }v\ll\alpha.$ (4)

The ratio of the Kramer’s logarithm and our fully summed numerical result is
shown in the bottom panel of Fig. 1. For $v\ll\alpha$ 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\leq
1000$, $\ell\leq n-1$, which can also be seen as a non-trivial check of our
code.111For a given $\alpha/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:
$\displaystyle\sum_{n}\frac{1}{n[n^{2}+(\alpha/v)^{2}]}\simeq[\log(\alpha/v)+\gamma_{E}]\text{,
for }v\ll\alpha\,.$ (5) From the denominator, one can estimate that for a
percentage accuracy, the maximum principle quantum number $n$ needs to be
roughly an order of magnitude larger than a given $\alpha/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\leq
1000,\ell\leq n-1$, provides a percentage-level accuracy for $\alpha/v\lesssim
100$ only, while $n\leq 100,\ell\leq n-1$ would require $\alpha/v\lesssim 10$
for the same accuracy. The Kramer’s logarithm has also been mentioned in
earlier dark matter related works Belotsky:2015osa ; An:2016gad ;
Petraki:2016cnz ; Smirnov:2019ngs .

Although the logarithm leads to a slope steeper than $1/v$, partial-wave
unitarity is not violated here as the sum over different $\ell$ 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\ll\alpha$. For instance, in this limit, the BSF cross section of the
$\mathrm{s}\to n\mathrm{p}$ processes summed over all $2\leq n\leq 1000$ is
larger by a constant factor $3.8$ than the $\rm{s}\to 2\rm{p}$. 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.

###  II.3 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}\geq 2$.

In particular, we consider pairs of non-relativistic Fermionic particles in
the fundamental and anti-fundamental representation of $SU(3)$, _i.e._ the
two-particle space is spanned by the direct sum of singlet and adjoint pair
states: $3\otimes\bar{3}=1\oplus 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

$\displaystyle(\sigma
v)_{n\ell}=\frac{C_{F}}{N_{c}^{2}}\frac{4\alpha}{3}\Delta
E^{3}|\bra{\psi_{n\ell}^{[{\bf 1}]}}\mathbf{r}\ket{\psi_{\mathbf{p}}^{[{\bf
adj}]}}|^{2},$ (6)

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 holds for the more
general $SU(N_{c})$ gauge group, as well as for a complex scalar field (see
_e.g._ Ref. Binder:2021otw ; Biondini:2023yxt  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 sections below.

Figure 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 $\ell\leq n-1$. From the truncation of the sum for
$n\leq 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 $\ell^{\prime}$ and the
corresponding partial-wave unitarity limit. For a given $\alpha$, unitarity is
violated when the rescaled ratio is above $1/\alpha^{3}$, for which two
examples are shown by the horizontal lines.

Adopting this setting, we evaluate the $SU(3)$ BSF cross section in Eq. (6)
for various $n,\ell$ 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\leq
1000$, $\ell\leq n-1$) and considering only the region of $\alpha/v\lesssim
100$ where $n\leq 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 $\sum_{n,\ell}(\sigma
v)_{n\ell}\propto v^{-\gamma}$ for $v\ll\alpha$ with a power law, we obtain
$\gamma\approx 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 $\ell^{\prime}$, needs to be considered. From Eq.
(6) we separate out this contribution by splitting Eq. (29) into two
contributions and denote the corresponding cross section by $(\sigma
v)^{\ell^{\prime}}_{n\ell}$, where the superscript denotes the fixed initial
scattering state angular momentum. For a given $\ell^{\prime}$, we then sum
$(\sigma v)^{\ell^{\prime}}_{n\ell}$ over all possible final state quantum
numbers $n,\ell$, compatible with the selection rules $\ell=\ell^{\prime}\pm
1$,

$(\sigma v)^{\ell^{\prime}}=\sum_{n,\ell}\,(\sigma
v)^{\ell^{\prime}}_{n\ell}\,.$ (7)

We require that this quantity respects the $\ell^{\prime}$th-partial wave
unitarity bound. Specifically, as done in Ref. vonHarling:2014kha ;
Smirnov:2019ngs ; Baldes:2017gzw ; Bottaro:2023wjv , we consider the partial-
wave unitarity cross section for $2\rightarrow 2$ total inelastic collisions
given by Griest:1989wd

$\displaystyle(\sigma
v)^{\ell^{\prime}}_{\text{uni.}}=\frac{\pi(2\ell^{\prime}+1)}{\mu^{2}v}.$ (8)

In the lower panel of Fig. 2, we show the introduced summed $(\sigma
v)^{\ell^{\prime}}$ divided by the corresponding $\ell^{\prime}$th-partial
wave unitarity cross section from Eq. (8). This ratio is multiplied by
$\alpha^{-3}$, leaving only a dependence on $\alpha/v$. In this way, the
partial-wave unitarity limits for a given $\alpha$ correspond to horizontal
lines with values equal to $\alpha^{-3}$. We show $\ell^{\prime}=0$ (s-wave)
to $\ell^{\prime}=3$ (f-wave) BSF cross section summed in $n\leq 1000$ and one
particular high initial angular momentum with $\ell^{\prime}=100$. The maximum
principle quantum number is taken as our numerical limit ($n\leq 1000$), which
is sufficient for the first four partial-wave summed cross sections to
converge for $\alpha/v\lesssim 100$. From these results, one can infer that
each summed $(\sigma v)^{\ell^{\prime}}$ grows faster than $1/v$, _i.e._ the
unitarity limit will be exceeded at a finite velocity, respectively, for any
$\ell^{\prime}$ 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 $\ell^{\prime}$ which is problematic.

Moreover, we observe the same situation for various $SU(N_{c})$. To make an
even more general statement, let us consider $(\sigma v)^{\ell^{\prime}}$ as a
function of the ratio $\alpha_{s}^{\text{eff}}/\alpha_{b}^{\text{eff}}$. Now,
$SU(N_{c})$ is a special case which lies in the region
$\alpha_{s}^{\text{eff}}/\alpha_{b}^{\text{eff}}=[-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
$\alpha_{s}^{\text{eff}}/\alpha_{b}^{\text{eff}}=1$. Our numerical results
suggest that in the range $\alpha_{s}^{\text{eff}}/\alpha_{b}^{\text{eff}}<1$,
$(\sigma v)^{\ell^{\prime}}$ scales stronger than $1/v$, while for
$\alpha_{s}^{\text{eff}}/\alpha_{b}^{\text{eff}}\geq 1$ no evidence of partial
wave-unitarity violation is found within our numerical boundaries.

Note that the mechanism behind the unitarization of BSF for
$\alpha_{s}^{\text{eff}}/\alpha_{b}^{\text{eff}}<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. Oncala:2019yvj  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)
$\alpha_{s}^{\text{eff}}\rightarrow 0$ and subsequently (ii)
$\tilde{\zeta}_{b}\rightarrow\infty$ where
$\tilde{\zeta}_{b}=\alpha_{b}^{\text{eff}}/(nv)$. Taking these limits, we
obtain the result for the s-wave case in $SU(N_{c})$

$\displaystyle(\sigma
v)^{\ell^{\prime}=0}_{n,\ell=1}\simeq\frac{C_{F}}{N_{c}^{2}}\frac{4\alpha}{3}\frac{32\pi\alpha^{\text{eff}}_{b}}{\mu^{2}}n(n^{2}-1),\text{
for (i) and (ii).}$ (9)

The two limits are justified for relative velocities which fulfill the
condition222For adjoint-to-singlet BSF in $SU(N_{c})$
$\alpha_{s}^{\text{eff}}=-\alpha/(2N_{c})$ and
$\alpha_{b}^{\text{eff}}=C_{F}\alpha$, where $C_{F}=(N_{c}^{2}-1)/(2N_{c})$.

$\displaystyle 2\pi|\alpha_{s}^{\text{eff}}|\ll
v\ll\frac{\alpha_{b}^{\text{eff}}}{2n^{2}}.$ (10)

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 $\alpha_{s}^{\text{eff}}=0$, then there is no lower bound on $v$ and
violation of s-wave unitarity can be shown with the above formula. Notice that
$\alpha_{s}^{\text{eff}}=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

$\displaystyle\dot{n}_{j}+3Hn_{j}=-\langle\sigma
v\rangle_{\text{eff}}[n^{2}_{j}-(n^{\text{eq}}_{j})^{2}]\,,$ (11)

where $n_{j}$ is the number density and $H$ the Hubble expansion rate. The
_effective cross section_ , $\langle\sigma v\rangle_{\text{eff}}$, includes
all the effects of pair annihilation as well as scattering-bound Ellis:2015vaa
and bound-bound transitions Mitridate:2017izz ; Garny:2021qsr ; Binder:2021vfo
. 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\propto T^{3}$, $H\propto T^{2}$, _i.e._ no (significant)
change in the relativistic degrees of freedom of the bath. Introducing the
yield as $Y_{j}\equiv n_{j}/s$ and parametrizing time by $x\equiv 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)\gg Y^{\mathrm{eq}}_{j}(x)$, also known as the
time of chemical decoupling, $x_{\text{cd}}$, one can neglect the impact of
$Y^{\mathrm{eq}}_{j}(x)$ in the Boltzmann equation. This allows for an
analytic solution for the yield evolution after chemical decoupling (see
_e.g._ Ref. Scherrer:1985zt ), which up to constants, can be estimated to
scale as

$Y_{j}(x_{0})\propto\frac{1}{\int_{x_{\text{cd}}}^{x_{0}}\mathrm{d}x\,x^{-2}\,\langle\sigma
v\rangle_{\text{eff}}(x)}.$ (12)

The integral converges for $x_{0}\to\infty$ only if $\langle\sigma
v\rangle_{\text{eff}}(x)$ grows slower than $x$ while for $\langle\sigma
v\rangle_{\text{eff}}\propto x^{\gamma}$ with $\gamma\geq 1$ the integral
diverges. Accordingly, the particle species only freezes out for $\gamma<1$
(typical WIMP) while the particle continues to deplete for $\gamma\geq 1$. The
critical value $\gamma=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}\propto x^{1-\gamma}$
for $x\gg x_{\text{cd}}$. In this case, the effective annihilation rate
$\Gamma_{\text{eff}}\equiv n_{j}\langle\sigma v\rangle_{\text{eff}}$ is
dynamically driven to be proportional to the Hubble rate
$\Gamma_{\text{eff}}\propto H$.

In the presence of bound states, the effective cross section introduced above
can be written as Garny:2021qsr ; Binder:2021vfo

$\big{\langle}\sigma v\big{\rangle}_{\text{eff}}=\big{\langle}\sigma
v\big{\rangle}_{\text{ann}}+\sum_{n,\ell}\big{\langle}\sigma
v\big{\rangle}_{n\ell}R_{n\ell}\,,$ (13)

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 Hisano:2003ec ; Hisano:2004ds . The second term contains the thermal
average of the BSF cross sections, denoted as $\big{\langle}\sigma
v\big{\rangle}_{n\ell}$. The summation over all bound-state quantum numbers
contains a dimensionless, temperature dependent quantity, which obeys $0\leq
R_{n\ell}\leq 1$.333Within the electric dipole approximation, bound states
with different spin are not directly coupled to each other. We thus leave the
spin sum implicit, see App. B for details. Thus, the presence of bound states
always 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\ell)$, $R_{i}$ can explicitly
be written as Garny:2021qsr ; Binder:2021vfo

$\displaystyle R_{i}$ $\displaystyle\equiv$ $\displaystyle
1-\sum_{k}(M^{-1})_{ik}\frac{\Gamma_{\text{ion}}^{k}}{\Gamma^{k}}\,,$ (14)
$\displaystyle M_{ik}$ $\displaystyle\equiv$
$\displaystyle\delta_{ik}-\frac{\Gamma_{\text{trans}}^{i\to
k}}{\Gamma^{i}}\,,$ (15) $\displaystyle\Gamma^{i}$ $\displaystyle\equiv$
$\displaystyle\Gamma_{\text{ion}}^{i}+\Gamma_{\text{dec}}^{i}+\sum_{k\neq
i}\Gamma_{\text{trans}}^{i\to k}\,.$ (16)

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 $\Gamma_{\text{ion}}^{n\ell}$ from
$\big{\langle}\sigma v\big{\rangle}_{n\ell}$. Note that the inclusion of
bound-to-bound transition rates _increases_ $\sum_{n,\ell}\big{\langle}\sigma
v\big{\rangle}_{n\ell}R_{n\ell}$ Binder:2021vfo .

###  III.1 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. vonHarling:2014kha ; Kamada:2019jch ; Biondini:2023zcz . 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, $\alpha$. 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

$\displaystyle\Gamma_{\text{de-ex}}^{n^{\prime}\ell^{\prime}\rightarrow
n\ell}=\frac{4\alpha(2\ell+1)}{3}\Delta
E^{3}|\bra{\psi_{n\ell}}\mathbf{r}\ket{\psi_{n^{\prime}\ell^{\prime}}}|^{2},$
(17)

see App. A.2. The excitation rate is related via detailed balance, see App. B.

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). As is well known,
the Sommerfeld effect in this case introduces a $1/v$ dependence of the
annihilation cross section for $v\ll\alpha$, leading to $\langle\sigma
v\rangle_{\text{eff}}\propto 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\ll\alpha$ (as seen in Fig. 1). For
$T$ much lower than the binding energy, this leads to $\big{\langle}\sigma
v\big{\rangle}_{10}\propto 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}\rightarrow 1$, resulting again in an overall
$\langle\sigma v\rangle_{\text{eff}}\propto 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
vonHarling:2014kha .

Figure 3:  Effective cross section for heavy Fermions charged under a $U(1)$
for constant coupling $\alpha=0.1$, including the contribution from bound
states for all $n,\,\ell$ 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 $\propto x^{1}$ for comparison.  Figure 4:
Effective cross section for heavy Fermion triplets under $SU(3)$, assuming
$\alpha(m)=0.025$. The critical scaling $\langle\sigma
v\rangle_{\text{eff}}\propto 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.

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\leq 10$ and $\ell\leq n-1$ (dashed blue
line) the effective cross section increases strongly until around $x\sim
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
$\gamma=1$ (indicated by the dashed red line) when including only $n\leq 10$.
However, the result for $n\leq 100$ and $\ell\leq n-1$ (blue solid line),
which includes about 5000 bound states and about $10^{6}$ transitions among
them, clearly shows that the effective 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 $10^{4}\lesssim x\lesssim 10^{5}$,
where $n\leq 100$ is trustworthy, we get a scaling of about $\langle\sigma
v\rangle_{\text{eff}}\propto 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, $\langle\sigma v\rangle_{\text{eff}}\propto x^{0.6}$,
is already found for including only $n\leq 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 $\langle\sigma
v\rangle_{\text{eff}}\times m^{2}$. Moreover, one can show analytically that
the scaling power is even unaffected when changing the value of $\alpha$,
which we have confirmed numerically in a wide range of $\alpha$. We also
explicitly checked that the inclusion of spin-triplet 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.444Interestingly, when considering processes
with infinitely many dark photons (classical picture), other works come to a
different conclusion Belotsky:2015osa ; Konstantin:2020rfx . We only included
ultra-soft processes with one dark photon. We show the impact of excited
states on the relic abundance within dark QED in App. C, refining earlier
results on this subject vonHarling:2014kha ; Smirnov:2019ngs .

###  III.2 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._ Asadi:2021pwo ; Biondini:2023zcz , 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 s-wave 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 no-
transition limit $\Gamma_{\text{trans}}^{i\to j}=0$ (_cf._ App. B.3) when
allowing only for the leading chromo electric dipole interactions. This
simplification allows us to focus exclusively on s-wave bound states which are
the only ones with a non-vanishing 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 $\left<\sigma
v\right>_{\text{eff}}\times m^{2}$ is independent of the choice of $m$.

In Fig. 4, we show $\left<\sigma v\right>_{\text{eff}}\times m^{2}$ for the
specific choice $\alpha(m)\equiv 0.025$ across the regime $\alpha(m/x)\lesssim
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._ $\left<\sigma
v\right>_{\text{eff}}\propto 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._ $\alpha=0.025$ at
all scales (light blue line, using $n\leq 1000$), we get a scaling power of
about $x^{0.9}$ for $x>10^{5}$, 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 $\alpha(m)\equiv
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\leq n\leq 1000$ leads to a reduction of the dark matter
abundance by around $50\%$ at $x=10^{9}$ for $m\sim 10^{6}\,{\rm{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 $\ell>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 $\chi$ being the dark matter candidate,
and a scalar mediator $\tilde{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 $\chi$ interacts with the SM only via the Yukawa interaction

${\cal L}_{\text{int}}=\lambda_{\chi}\tilde{q}\bar{q}_{R}\chi+\text{h.c.}\,,$
(18)

while the mediator $\tilde{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_{\tilde{q}}>m_{\chi}$ 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._ Garny:2015wea , that are
being actively considered in the context of LHC dark matter searches, see
_e.g._ Arina:2023msd  for a recent account on the subject.555The term
_mediator_ refers to the dark matter-SM interaction, not to be confused with
the long-range force carrier which in this case is the gluon (and photon).
Occasionally, this class of models has alternatively been dubbed _charged
parent particle model_.

Dark matter production is governed by the interaction of $\chi$ with the
mediator field $\tilde{q}$, controlled by $\lambda_{\chi}$, as well as the
dynamics of the mediator itself, largely driven by its gauge interactions. In
particular, as the $\tilde{q}$ and $\tilde{q}^{\dagger}$ particles are color
and electrically charged, they can form bound states, that have an important
impact on the freeze-out Mitridate:2017izz ; Biondini:2018pwp ; Harz:2018csl ;
Biondini:2018ovz ; Decant:2021mhj ; Bollig:2021psb ; Garny:2021qsr ;
Becker:2022iso . Here, we are particularly interested in including excited
bound states, following Harz:2018csl ; Garny:2021qsr .

###  IV.1 Review of production mechanisms

There are three distinct possibilities for how the freeze-out dynamics occurs,
known as coannihilation Griest:1990kh ; Edsjo:1997bg , conversion-driven
freeze-out Garny:2017rxs ; DAgnolo:2017dbv  and superWIMP production
Covi:1999ty ; Feng:2003uy , respectively. In addition, the model can also
feature freeze-in production McDonald:2001vt ; Asaka:2005cn ; Hall:2009bx .

All of them can be described by the following set of coupled Boltzmann
equations for the yields $Y_{j}$,

$\displaystyle\frac{\mathrm{d}Y_{\tilde{q}}}{\mathrm{d}x}=\frac{1}{3H}\frac{\mathrm{d}s}{\mathrm{d}x}\Bigg{[}\frac{1}{2}\big{\langle}\sigma_{\tilde{q}\tilde{q}^{\dagger}}v\big{\rangle}_{\text{eff}}\left(Y_{\tilde{q}}^{2}-Y_{\tilde{q}}^{\mathrm{eq}\,2}\right)$
(19)
$\displaystyle\;+\big{\langle}\sigma_{\chi\tilde{q}}v\big{\rangle}\left(Y_{\chi}Y_{\tilde{q}}-Y_{\chi}^{\mathrm{eq}}Y_{\tilde{q}}^{\mathrm{eq}}\right)+\frac{\Gamma_{\text{conv}}^{\tilde{q}\to\chi}}{s}\left(Y_{\tilde{q}}-Y_{\chi}\frac{Y_{\tilde{q}}^{\mathrm{eq}}}{Y_{\chi}^{\mathrm{eq}}}\right)\Bigg{]}\,,$
$\displaystyle\frac{\mathrm{d}Y_{\chi}}{\mathrm{d}x}=\frac{1}{3H}\frac{\mathrm{d}s}{\mathrm{d}x}\Bigg{[}\big{\langle}\sigma_{\chi\chi}v\big{\rangle}\left(Y_{\chi}^{2}-Y_{\chi}^{\mathrm{eq}\,2}\right)$
(20)
$\displaystyle\;+\big{\langle}\sigma_{\chi\tilde{q}}v\big{\rangle}\left(Y_{\chi}Y_{\tilde{q}}-Y_{\chi}^{\mathrm{eq}}Y_{\tilde{q}}^{\mathrm{eq}}\right)-\frac{\Gamma_{\text{conv}}^{\tilde{q}\to\chi}}{s}\left(Y_{\tilde{q}}-Y_{\chi}\frac{Y_{\tilde{q}}^{\mathrm{eq}}}{Y_{\chi}^{\mathrm{eq}}}\right)\Bigg{]}\,,$

where $x=m_{\tilde{q}}/T$, $s$ is the entropy density, $H$ the Hubble rate,
and

$Y_{j}^{\mathrm{eq}}=\frac{g_{j}}{s}\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}\mathrm{e}^{-\sqrt{m_{j}^{2}+p^{2}}/T},$
(21)

all of which depend on $x$. Here $g_{j}$ denotes the number of internal
degrees of freedom, with $g_{\tilde{q}}\equiv 2N_{c}=6$ denoting the sum of
$\tilde{q}$ and $\tilde{q}^{\dagger}$ densities, $g_{\chi}=2$, and $g_{{\cal
B}_{n\ell}}=2\ell+1$ for bound states with angular momentum $\ell$, 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 $\tilde{q}$ and
$\tilde{q}^{\dagger}$ density in $Y_{\tilde{q}}$.

Eqs. (19) and (20) contain the following collision terms:

  1. 1.

The effective cross section
$\big{\langle}\sigma_{\tilde{q}\tilde{q}^{\dagger}}v\big{\rangle}_{\text{eff}}$
includes direct annihilation (including Sommerfeld enhancement following
Ibarra:2015nca ; Garny:2017rxs ) 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 Garny:2021qsr .

  2. 2.

The rate $\Gamma_{\text{conv}}^{\tilde{q}\to\chi}$ describes the conversion
rate of $\tilde{q}$ into $\chi$ particles. It is controlled by the Yukawa
coupling, $\Gamma_{\text{conv}}^{\tilde{q}\to\chi}\propto\lambda_{\chi}^{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\tilde{q}\to Y\chi$ with appropriate SM particles
$X,Y$, while at low temperatures the decay process $\tilde{q}\to q\chi$
dominates. Accordingly, in the low temperature limit – relevant for the
superWIMP mechanism considered below – it reads

$\Gamma_{\text{conv}}^{\tilde{q}\to\chi}=\Gamma_{\tilde{q}\to
q\chi}\frac{K_{1}\left(m_{\tilde{q}}/T\right)}{K_{2}\left(m_{\tilde{q}}/T\right)}\,,$
(22)

where $\Gamma_{\tilde{q}\to
q\chi}=\lambda_{\chi}^{2}/(16\pi)\,m_{\tilde{q}}\left(1-m_{\chi}^{2}/m_{\tilde{q}}^{2}\right)^{2}$
is the vacuum decay rate of a single mediator particle in the limit $m_{q}\to
0$ (not to be confused with the bound state decay rates, that are dominated by
the strong interaction).

  3. 3.

The dark matter pair annihilation rate
$\big{\langle}\sigma_{\chi\chi}v\big{\rangle}\propto\lambda_{\chi}^{4}$ and
coannihilation rate
$\big{\langle}\sigma_{\chi\tilde{q}}v\big{\rangle}\propto\lambda_{\chi}^{4}$
are strongly suppressed for $\lambda_{\chi}\ll 1$ and practically irrelevant
within the conversion-driven and superWIMP regimes.

In addition, the Boltzmann equations could be complemented by collision terms
for the conversion process $\tilde{q}\tilde{q}^{\dagger}\to\chi\chi$, which
are, however, negligible within the conversion-driven and superWIMP regimes as
they are, again, proportional to $\lambda_{\chi}^{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 $\chi$ into $\tilde{q}$ particles,

$\Gamma_{\text{conv}}^{\chi\to\tilde{q}}=\Gamma_{\text{conv}}^{\tilde{q}\to\chi}\frac{Y_{\tilde{q}}^{\mathrm{eq}}}{Y_{\chi}^{\mathrm{eq}}}\to\Gamma_{\tilde{q}\to
q\chi}\frac{g_{\tilde{q}}m_{\tilde{q}}^{3/2}}{g_{\chi}m_{\chi}^{3/2}}e^{-\frac{m_{\tilde{q}}-m_{\chi}}{T}}\,,$
(23)

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),

$\displaystyle\Gamma_{\text{conv}}^{\chi\to\tilde{q}}\gg H(m_{\tilde{q}})$
$\displaystyle\text{coannihilation}\,,$ (24a)
$\displaystyle\Gamma_{\text{conv}}^{\chi\to\tilde{q}}\sim H(m_{\tilde{q}})$
$\displaystyle\text{conversion-driven}\,,$ (24b)
$\displaystyle\Gamma_{\text{conv}}^{\chi\to\tilde{q}}\ll H(m_{\tilde{q}})$
$\displaystyle\text{superWIMP/freeze-in}\,.$ (24c)

  1. (a)

In the coannihilation regime the ${\tilde{q}}$ and $\chi$ 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
Griest:1990kh ; Edsjo:1997bg . Within the coannihilation regime, the dark
matter abundance is determined by the cross sections
$\big{\langle}\sigma_{\tilde{q}\tilde{q}^{\dagger}}v\big{\rangle}$,
$\big{\langle}\sigma_{\chi\chi}v\big{\rangle}$ and
$\big{\langle}\sigma_{\chi\tilde{q}}v\big{\rangle}$ (and in addition, as for
all cases, the bound state dynamics), and generically here
$\lambda_{\chi}\sim{\cal O}(1)$ Garny:2015wea ; Harz:2018csl ; Garny:2021qsr ;
Becker:2022iso .

  2. (b)

In the conversion-driven case, the freeze-out of chemical equilibrium among
$\chi$ and $\tilde{q}$ drives the dynamics and the size of the conversion rate
$\Gamma_{\text{conv}}^{\tilde{q}\to\chi}$ largely influences the dark matter
abundance Garny:2017rxs ; DAgnolo:2017dbv . In addition, the efficiency by
which the $\tilde{q}$ (and $\tilde{q}^{\dagger}$) abundance is depleted is
relevant Garny:2017rxs , controlled by
$\big{\langle}\sigma_{\tilde{q}\tilde{q}^{\dagger}}v\big{\rangle}$ and bound
state effects Garny:2021qsr . The conversion-driven case occurs for small
couplings, typically $\lambda_{\chi}\sim{\cal O}(10^{-6})$, for which the
$\chi\chi$ and $\chi\tilde{q}$ terms in Eqs. (19) and (20) can be safely
neglected.

  3. (c)

Finally, in the superWIMP scenario Covi:1999ty ; Feng:2003uy , the mediator
has an even smaller decay rate, and can usually be considered as stable while
the freeze-out of $\tilde{q}\tilde{q}^{\dagger}$ annihilation occurs. The
population of remaining $\tilde{q}$ and $\tilde{q}^{\dagger}$ particles then
decays into $\chi$ at a temperature $T$ for which
$H(T)\sim\Gamma_{\text{conv}}^{\chi\to\tilde{q}}$, thereby generating the dark
matter abundance Garny:2018ali . Technically, this means that the terms
corresponding to inverse decays in Eqs. (19) and (20) can be neglected, in
addition to those for $\chi\chi$ and $\chi\tilde{q}$ annihilation, while the
size of $\big{\langle}\sigma_{\tilde{q}\tilde{q}^{\dagger}}v\big{\rangle}$ and
the bound state dynamics are most important Garny:2018ali ; Decant:2021mhj ;
Bollig:2021psb . 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
$\tilde{q}$ (and $\tilde{q}^{\dagger}$) produces one dark matter particle.
Note that in addition to the superWIMP contribution a contribution from
freeze-in McDonald:2001vt ; Asaka:2005cn ; Hall:2009bx  has to be considered
which stems from inefficient decays (or scatterings) around $x\sim 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._ Garny:2018ali ; Decant:2021mhj . 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.

###  IV.2 Bound state rates and processes

The impact of bound states on $Y_{\tilde{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 $\tilde{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 wave-functions as well as BSF and decay are driven by QCD,
while QED is relevant for transitions among the various energy levels
Garny:2021qsr . 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,

$(\tilde{q}\tilde{q}^{\dagger})^{[8]}\to{\cal B}_{n\ell}^{[1]}+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)}=-\alpha_{s(b)}^{\text{eff}}/r$ for the scattering ($s$) and bound
states ($b$) are

$\alpha_{s}^{\text{eff}}=-\frac{1}{6}\alpha_{s}(\mu_{s}),\quad\alpha_{b}^{\text{eff}}=\frac{4}{3}\alpha_{s}(\mu_{b})\,,$
(26)

with running strong coupling666In all numerical computations involving running
$\alpha_{s}$ we used RunDec 3 Herren:2017osy  to evaluate the SM strong
coupling (employing 5-loop running). evaluated at $\overline{\text{MS}}$-scale
$\mu_{s}=m_{\tilde{q}}v_{\text{rel}}/2$ for the scattering state, and at the
Bohr momentum scale $\mu_{b}=m_{\tilde{q}}\alpha_{b}^{\text{eff}}/2/n$ for the
bound state ${\cal B}_{n\ell}$. 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
$\mu_{\text{BSF}}=m_{\tilde{q}}/{4}\,\left(v^{2}+(\alpha_{b}^{\text{eff}}/n)^{2}\right)$.
Bound state formation is also possible via an electromagnetic dipole
transition, but is negligible due to the smaller interaction strength
Garny:2021qsr  (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:

${\cal B}_{n\ell}^{[1]}\leftrightarrow{\cal B}_{n^{\prime}\ell\pm
1}^{[1]}+\gamma\,,$ (27)

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^{\prime}$ and the electromagnetic
fine structure constant $\alpha_{\text{EM}}=1/128.9$ in the coupling prefactor
of Eq. (17).

For bound state decay we include the process ${\cal B}_{n,\ell=0}\to 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 Garny:2021qsr , 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
suppressed 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\leq 6$) that electric
quadrupole transitions have a negligible impact on the effective cross
section. Furthermore, the decay of $\ell=1$ bound states is suppressed by at
least a factor $\alpha_{s}^{3}$ relative to its de-excitation rate when
assuming a power counting $\alpha_{\text{EM}}\sim\mathcal{O}(\alpha_{s}^{2})$.
This is due to the higher derivatives of the radial wavefunction entering for
higher $\ell$, as well as the fact that decays into two gluons vanish for
$\ell=1$ states at tree-level, leading to additional suppression due to a
three-gluon decay or two-gluon decay at one-loop Cacciari:2015ela . 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.

Figure 5: Ionization, decay and transition rates for all bound state levels
$(n,\ell)$ up to $n=4$ for $m_{\tilde{q}}=4\times 10^{6}\,$GeV. The gray
dashed vertical line indicates the temperature that corresponds to the binding
energy $E_{\mathcal{B}}$ of the respective bound state. Note that the
transition rates contain all possible excitations and de-excitations (here
summed up to $n^{\prime}=20$).  Figure 6:  Effective cross section for the
colored and electrically charged mediator $\tilde{q}$. It includes the
contribution from bound state levels $(n,\ell)$ for all $\ell$ and up to
$n\leq 1,10,100$ (blue lines) accounting for all dipole transitions among
them. Also shown is the no-transition limit for $n\leq 1000$ (green). Both
grow more steeply than $\propto x$ (indicated by a thin black line) when
sufficiently high $n$ are included, implying mediator depletion without
freeze-out. Here $m_{\tilde{q}}=4\times 10^{6}$GeV.

##  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
$\tilde{q}$ are well separated in time, such that the late time $\chi$
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 $\tau_{\tilde{q}}=1/\Gamma_{\tilde{q}\to q\chi}$. 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.2.

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 $\alpha_{s}$ in the perturbative regime ($T>1$GeV) only.
Numerical results presented throughout this section are shown for the
representative benchmark point $m_{\tilde{q}}=4\times 10^{6}$ GeV. We discuss
the mass dependence of our findings in Sec. V.3.

###  V.1 Effective cross section

The abundance of the mediator $\tilde{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\leq 4$
and for all $\ell\leq n-1$, for a wide range of $x=m_{\tilde{q}}/T$
corresponding to $m_{\tilde{q}}/10\geq T\geq 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 $\tilde{q}\tilde{q}^{\dagger}$ 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,\ell)$ into any higher or lower state is shown by the blue lines in
Fig. 5. For the ground state $(1,0)$ only excitation occurs 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 $\Delta\ell=\pm 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 $\ell=0$ states
into a pair of gluons, which is the analog of mediator pair annihilation. This
rate is practically temperature independent for $x\gg 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 $\ell=0$, and QED-mediated transitions (to
lower levels) for $\ell\geq 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 ${\cal B}_{n\ell}$ with $0\leq\ell\leq 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,\,\ell$ only starts to contribute significantly once the temperature drops
down to roughly its respective binding energy, $T\sim E_{{\cal B}_{n\ell}}$.
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\ll E_{{\cal B}_{n\ell}}$ in contrast, its contribution is negligible due to
the repulsive potential in the scattering state Harz:2018csl ; Garny:2021qsr .

Introducing $x_{n}\equiv m_{\tilde{q}}/E_{{\cal B}_{n\ell}}$, we find
$x_{1}\simeq 7\times 10^{2}$, $x_{10}\simeq 5\times 10^{4}$, $x_{100}\simeq
3\times 10^{6}$ for our benchmark $m_{\tilde{q}}=4\times 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\leq 50$ and $n\leq 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 $\propto 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.2. 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 $\propto 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 $\propto x$, _i.e._ the presence of bound-to-
bound transitions causes a super-critical behavior in our model even without
running coupling effects.

###  V.2 Relic abundance

Figure 7:  _Upper panel:_ Abundance evolution of $\chi$ (red) and $\tilde{q}$
(blue) for $m_{\tilde{q}}=4\times 10^{6}$ GeV and $\Gamma_{\tilde{q}\to
q\chi}=10^{-17}$ GeV. When including no (SE only) or one ($n=1$) bound state,
the mediator yield $Y_{\tilde{q}}$ freezes out and then subsequently transfers
its abundance to $Y_{\chi}$ via $\tilde{q}\to q\chi$. Taking excited states
and transitions among them into account leads to a continuous decrease of
$Y_{\tilde{q}}$ (dashed and solid blue line for $n\leq 10$ and $n\leq 100$,
respectively) that is only terminated by mediator decay once
$\Gamma_{\tilde{q}\to q\chi}\gtrsim H$. _Lower panel:_ Ratio of $\chi$ yield
when including bound states up to $n\leq 1,10,100$ over the result without
bound state contributions.  Figure 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 $\Gamma_{\tilde{q}\to q\chi}=2.5\times 10^{-14}$, $5\times 10^{-16}$,
$10^{-17}$ GeV. The long dashed line shows the limit $\Gamma_{\tilde{q}\to
q\chi}=0$, which decreases due to the large contribution of excited states to
the effective cross section. The final value of $Y_{\chi}$ therefore depends
on $\Gamma_{\tilde{q}\to q\chi}$ even in the superWIMP regime where the $\chi$
abundance generated via freeze-in at low $x$ is negligible.

The evolution of the yields $Y_{\tilde{q}}(x)$ and $Y_{\chi}(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._
Decant:2021mhj ; Bollig:2021psb ). In these cases, the yield
$Y_{\tilde{q}}(x)$ freezes out at $x\sim{\cal 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\sim\Gamma_{\tilde{q}\to
q\chi}$. At this point, the mediator particles decay into dark matter, such
that the final yield $Y_{\chi}$ is identical to the freeze-out value of
$Y_{\tilde{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_{\tilde{q}}$ never freezes out. This can be seen in the solid line in Fig.
7, where bound states up to $n\leq 100$ are included. The depletion of the
number density is dominated by the effective cross section, until the decay of
the mediator, $\tilde{q}\to q\chi$, becomes efficient (here $x\gtrsim
10^{6}$). In fact, the effective annihilation rate
$\Gamma_{\text{eff}}=n_{\tilde{q}}\langle\sigma v\rangle_{\text{eff}}$ is kept
on the edge of being efficient, _i.e._ $\Gamma_{\text{eff}}\sim 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\leq 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 $\sim 2$ in the final
abundance. When considering the full bound state effects ($n\leq 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 $\alpha_{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._ $\Gamma_{\tilde{q}\to q\chi}=0$.) For all
decay rates, superWIMP production provides the dominant contribution to the
dark matter density. While for the largest decay rate, $\Gamma_{\tilde{q}\to
q\chi}=2.5\times 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_{\tilde{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_{\chi}$ (red lines) for $x\to\infty$. 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.

###  V.3 Implications

Figure 9:  Dark matter mass $m_{\chi}$ and decay rate $\Gamma_{\tilde{q}\to
q\chi}$ of the colored $t$-channel mediator for which the final yield matches
the observed dark matter density $\Omega_{\chi}h^{2}=0.12$, using
$m_{\tilde{q}}=4\times 10^{6}$ GeV. The abundance is set by the superWIMP
mechanism for $\Gamma_{\tilde{q}\to q\chi}\lesssim 10^{-13}$ GeV, while
freeze-in dominates for larger decay rates. When taking only Sommerfeld-
enhanced $\tilde{q}\tilde{q}^{\dagger}$ annihilation into account (dotted
line) the relic density is independent of $\Gamma_{\tilde{q}\to q\chi}$ within
the superWIMP regime, _i.e._ the dotted line is horizontal. The same is true
when taking the ground state $n=1$ into account (dot-dashed line). When
including excited states, the mediator depletion due to BSF, transitions and
decay of bound states continues until eventually $H\sim\Gamma_{\tilde{q}\to
q\chi}$, such that the relic density does depend on the mediator lifetime,
corresponding to the black dashed ($n\leq 10$) and black solid ($n\leq 100$)
lines. For $\Gamma_{\tilde{q}\to q\chi}\gtrsim 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-$\alpha$
forest observations Decant:2021mhj .

In Fig. 9, we finally show the dark matter mass, $m_{\chi}$ (left axis
labeling), for which the final $\chi$ abundance (displayed using the right
axis labeling) yields the observed dark matter relic density,
$\Omega_{\chi}h^{2}\simeq 0.12$ Planck:2018vyg , as a function of the mediator
decay rate, $\Gamma_{\tilde{q}\to q\chi}$. For $\Gamma_{\tilde{q}\to
q\chi}\gtrsim 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 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 $\tilde{q}\to q\chi$. This
implies that the final dark matter abundance does, in fact, depend on
$\Gamma_{\tilde{q}\to q\chi}$ and, hence, on the dark matter coupling.
Therefore also the dark matter mass $m_{\chi}$ for which
$\Omega_{\chi}h^{2}\simeq 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 ($\ell\leq n-1$) with $n\leq 10$ and $100$, respectively. The
smaller dark matter abundance implies a larger dark matter mass, as compared
to the case without excited states.

For $\Gamma_{\tilde{q}\to q\chi}\gtrsim 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,
$\Gamma_{\tilde{q}\to q\chi}\geq 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-$\alpha$ forest observations. As shown in Ref. Decant:2021mhj , this
results in a lower bound on the dark matter mass,

$\frac{m_{\chi}}{\text{keV}}>3.8\times
x_{\text{decay}}\left(\frac{106.75}{g_{*S}(x_{\text{decay}})}\right)^{\\!1/3\,}\,,$
(28)

that can easily reach into the GeV range. Here,
$x_{\text{decay}}=(\Gamma_{\tilde{q}\to q\chi}/H(m_{\tilde{q}}))^{-1/2}$ is
the temperature parameter at which the decay becomes efficient,
$\Gamma_{\tilde{q}\to q\chi}=H(m_{\tilde{q}}/x_{\text{decay}})$, _cf._ upper
axis labeling in Fig. 9. The formula assumes $\Omega_{\chi}h^{2}=0.12$ and
$m_{q}\ll m_{\tilde{q}}$.

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\times 10^{-16})$ GeV
would be excluded, the inclusion of excited bound states reveal that the
entire region with $\Gamma_{\tilde{q}\to q\chi}\geq 10^{-17}$ GeV remains
unchallenged by the considered Lyman-$\alpha$ constraint.

So far, we have focused on the benchmark mass $m_{\tilde{q}}=4\times 10^{6}$
GeV. For smaller masses, the effective cross section becomes larger and the
mediator abundance decreases. Therefore, for a given $x_{\text{decay}}$, the
dark matter mass that leads to $\Omega h^{2}=0.12$ increases and the
Lyman-$\alpha$ bound become less constraining. Accordingly, we find that
current Lyman-$\alpha$ forest constraints do not challenge the superWIMP
scenario for $m_{\tilde{q}}<4\times 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_{\text{decay}}<m_{\tilde{q}}/1\,$GeV. However, as the so-defined maximal
$x_{\text{decay}}$ becomes smaller with smaller $m_{\tilde{q}}$, highly
excited states become less relevant somewhat diminishing the large effect of
bound states towards small $m_{\tilde{q}}$ in our setup.

For masses larger than $4\times 10^{6}$ GeV, the cosmologically viable dark
matter mass decreases and Lyman-$\alpha$ constraints become more restrictive.
In particular, starting from masses somewhat above $m_{\tilde{q}}=4\times
10^{6}$ GeV, they impose an upper bound on $x_{\text{decay}}$ that is more
restrictive than the above-mentioned perturbativity condition. As a
consequence of the tightening Lyman-$\alpha$ constraints on the allowed range
of mediator decay rates, higher excitations become less important also towards
larger masses. Eventually, for $m_{\tilde{q}}>4\times 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 $\mathcal{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 $\propto 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 $\bar{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 $\langle\sigma v\rangle_{\text{eff}}\propto
x^{0.6}$ towards large $x$. For dark QCD, $\langle\sigma
v\rangle_{\text{eff}}$ grows slightly slower (faster) than $\propto 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 $\propto 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 $\propto 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 $\Omega h^{2}$ towards larger values by up to an order of magnitude,
thereby weakening constraints from structure formation through Lyman-$\alpha$
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 Heeck:2023lxl , the
late decay of a color-charged scalar can generate the matter-antimatter
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_{\text{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
Asadi:2021pwo ; Asadi:2022vkc . 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.

## Acknowledgments

We are thankful to Martin Beneke, Florian Herren, and Gramos Qerimi for
discussions. This work was supported by the DFG Collaborative Research
Institution Neutrinos and Dark Matter in Astro- and Particle Physics (SFB
1258) and the Collaborative Research Center TRR 257. J.H. acknowledges support
by the Alexander von Humboldt foundation via the Feodor Lynen Research
Fellowship for Experienced Researchers and Feodor Lynen Return Fellowship.

##  Appendix A Evaluation of large _n_ dipole transitions

The evaluation of ultra-soft transitions such as bound-state formation and de-
excitation for large $n$ is a demanding task. Here, we briefly summarize the
derivation of the used expressions that enable a numerically stable and
sufficiently fast computation of the scattering-to-bound and bound-to-bound
dipole transitions up to about $n=1000$ and $n=100$, respectively. A
Mathematica notebook containing the implemented expressions is available on
request.

###  A.1 Scattering-to-bound

We turn now to the evaluation of the scattering-to-bound dipole transitions
matrix elements. The angular part can be performed as Garny:2021qsr :

$\displaystyle\sum_{m}|\bra{\psi_{n\ell
m}}\mathbf{r}\ket{\psi_{\mathbf{p}}^{\prime}}|^{2}=\frac{4\pi}{|\mathbf{p}|^{5}}\sum_{\ell^{\prime}}(\ell^{\prime}\delta_{\ell^{\prime},\ell+1}+\ell\delta_{\ell,\ell^{\prime}+1})|I_{R}|^{2},$
(29)

where we have $|{\bf p}|=mv/2$. Our starting point for the evaluation of the
radial integral part is (A.13) in Ref. Garny:2021qsr :

$\displaystyle I_{R,\text{BSF}}=$
$\displaystyle\frac{2\zeta_{b}^{3/2}(\ell+\ell^{\prime}+3)!}{n^{2}\sqrt{(n-\ell-1)!(n+\ell)!}}\frac{2^{\ell^{\prime}}e^{\pi\zeta_{s}/2}}{|\Gamma(1+\ell^{\prime}-i\zeta_{s})|}\times$
$\displaystyle\left(\frac{2\zeta_{b}}{n}\right)^{\ell}\left(\frac{d}{dt}\right)^{n-\ell-1}\frac{1}{(1-t)^{2\ell+2}}\times$
$\displaystyle\int_{0}^{1}ds\frac{s^{\ell^{\prime}-i\zeta_{s}}(1-s)^{\ell^{\prime}+i\zeta_{s}}}{\left(\frac{\zeta_{b}}{n}\frac{1+t}{1-t}+i(2s-1)\right)^{\ell+\ell^{\prime}+4}}\Big{|}_{t=0}\,,$
(30)

where $\zeta_{b}=\alpha^{\text{eff}}_{b}/v$ and
$\zeta_{s}=\alpha^{\text{eff}}_{s}/v$. To arrive at this expression, the
representation of the bound-state wave function in terms of associated
Laguerre polynomials and their generating function was used while for the
scattering state the Hypergeometric function ${}_{1}F_{1}$ in terms of the
integral representation was used. In this work, we directly evaluate the
integral first as

$\displaystyle\int_{0}^{1}ds\frac{s^{\ell^{\prime}-a}(1-s)^{\ell^{\prime}+a}}{\left(-ib+i(2s-1)\right)^{\ell+\ell^{\prime}+4}}=$
$\displaystyle(-i)^{-\ell-\ell^{\prime}}\,\frac{\Gamma(1-a+\ell^{\prime})\Gamma(1+a+\ell^{\prime})}{(1+b)^{4+\ell+\ell^{\prime}}\Gamma(2(1+\ell^{\prime}))}\,\times$
$\displaystyle{}_{2}F_{1}\left(1-a+\ell^{\prime},4+\ell+\ell^{\prime},2(1+\ell^{\prime}),\frac{2}{1+b}\right)\,,$
(31)

where $a=i\zeta_{s},b=i\frac{\zeta_{b}}{n}\frac{1+t}{1-t}$. Applying the $t$
derivatives in Eq. (30) to the Hypergeometric function ${}_{2}F_{1}$, allows
us to obtain the following recursive relations:

$\displaystyle|I_{R,\text{BSF}}^{\ell^{\prime}=\ell+1}|=2n^{3}\sqrt{(\ell^{2}+\zeta_{s}^{2})((\ell+1)^{2}+\zeta_{s}^{2})}P(n,\ell)\times$
$\displaystyle\Big{|}((2+\ell)\tilde{\zeta_{b}}+\zeta_{s})R_{n,\ell}(n-\ell-5)$
$\displaystyle-2((2+\ell)\tilde{\zeta_{b}}+2\zeta_{s})R_{n,\ell}(n-\ell-4)$
$\displaystyle+6\zeta_{s}R_{n,\ell}(n-\ell-3)$
$\displaystyle+2((2+\ell)\tilde{\zeta_{b}}-2\zeta_{s})R_{n,\ell}(n-\ell-2)$
$\displaystyle+(\zeta_{s}-(2+\ell)\tilde{\zeta_{b}})R_{n,\ell}(n-\ell-1)\Big{|}\;,$
(32)

and

$\displaystyle|I_{R,\text{BSF}}^{\ell^{\prime}=\ell-1}|=n^{3}P(n,\ell)\times$
$\displaystyle\Big{|}-\big{[}(\ell(1+\ell)\tilde{\zeta_{b}}(-3+(1+2\ell))\tilde{\zeta_{b}}^{2})$
$\displaystyle+(-1-3\ell+3(1+\ell)(1+2\ell)\tilde{\zeta_{b}}^{2})\zeta_{s}$
$\displaystyle+6(1+\ell)\tilde{\zeta_{b}}\zeta_{s}^{2}+2\zeta_{s}^{3}\big{]}R_{n,\ell}(n-\ell-5)$
$\displaystyle-2\big{[}\ell(1+\ell)\tilde{\zeta_{b}}(3+(1+2\ell)\tilde{\zeta_{b}}^{2})+2\zeta_{s}+6\ell\zeta_{s}$
$\displaystyle-6(1+\ell)\tilde{\zeta_{b}}\zeta_{s}^{2}-4\zeta_{s}^{3}\big{]}R_{n,\ell}(n-\ell-4)$
$\displaystyle+\big{[}6\zeta_{s}(1+3\ell+(1+\ell)(1+2\ell)\tilde{\zeta_{b}}^{2}-2\zeta_{s}^{2})\big{]}R_{n,\ell}(n-\ell-3)$
$\displaystyle+2\big{[}\ell(1+\ell)\tilde{\zeta_{b}}(3+(1+2\ell)\tilde{\zeta_{b}}^{2})-2(1+3\ell)\zeta_{s}$
$\displaystyle-6(1+\ell)\tilde{\zeta_{b}}\zeta_{s}^{2}+4\zeta_{s}^{3}\big{]}R_{n,\ell}(n-\ell-2)$
$\displaystyle+\big{[}\ell(1+\ell)\tilde{\zeta_{b}}(-3+(1+2\ell)\tilde{\zeta_{b}}^{2})+\zeta_{s}+3\ell\zeta_{s}$
$\displaystyle-3(1+\ell)(1+2\ell)\tilde{\zeta_{b}}^{2}\zeta_{s}$
$\displaystyle+6(1+\ell)\tilde{\zeta_{b}}\zeta_{s}^{2}-2\zeta_{s}^{3}\big{]}R_{n,\ell}(n-\ell-1)\Big{|}\;,$
(33)

where $\tilde{\zeta_{b}}=\zeta_{b}/n$. The common prefactor is

$\displaystyle P(n,\ell)$
$\displaystyle=\sqrt{\frac{(n-\ell-1)!}{(n+\ell)!}}\frac{2^{2+2\ell}\tilde{\zeta_{b}}^{\ell+\frac{3}{2}}}{n^{\frac{7}{2}}\zeta_{s}(1+\tilde{\zeta_{b}}^{2})^{3+\ell}}\times$
$\displaystyle\sqrt{\frac{2\pi\zeta_{s}}{1-e^{-2\pi\zeta_{s}}}}e^{-2\zeta_{s}\text{arccot}(\tilde{\zeta_{b}})}\sqrt{\prod_{j=0}^{\ell-1}\left(j^{2}+\zeta_{s}^{2}\right)}$
(34)

and the common recursion is given by

$\displaystyle R_{n,\ell}(x)=\left\\{\begin{array}[]{lr}0&x<0\\\ 1&x=0\\\
\frac{-2(3+\ell)(\tilde{\zeta_{b}}^{2}-1)+4\tilde{\zeta_{b}}\zeta_{s}}{1+\tilde{\zeta_{b}}^{2}}&x=1\\\
\frac{2(2+\ell+x)(1-\tilde{\zeta_{b}}^{2})+4\tilde{\zeta_{b}}\zeta_{s}}{x(1+\tilde{\zeta_{b}}^{2})}R_{n,\ell}(x-1)&\\\
\;\;-\,\frac{(4+2\ell+x)}{x}R_{n,\ell}(x-2)&{\rm else}\,.\end{array}\right.$
(40)

The expressions allow for a fast generation of the matrix elements entering
the BSF cross section, also for large $n,\ell$. They have been cross-checked
against the expressions in Ref. Garny:2021qsr  for $n\leq 10$ and $\ell\leq
n-1$ for the QCD case, as well as for the simpler QED limit
$\zeta_{s}=\zeta_{b}=\zeta$. To achieve the results presented in this work up
to $n=1000$ we split the prefactor (34) into two multiplicative pieces to
avoid numerical underflow.

###  A.2 Bound-to-bound

We turn now to the evaluation of the bound-to-bound dipole transition matrix
elements. The angular part can be performed as Garny:2021qsr :

$\displaystyle\sum_{m,m^{\prime}}|\bra{\psi_{n\ell
m}}\mathbf{r}\ket{\psi_{n^{\prime}\ell^{\prime}m^{\prime}}^{\prime}}|^{2}$
$\displaystyle=$ $\displaystyle(\ell^{\prime}\delta_{\ell^{\prime},\ell+1}$
$\displaystyle+\ell\delta_{\ell,\ell^{\prime}+1})|I_{R,\text{trans}}|^{2}\,,$
(41)

leaving only the radial integral over the initial and final bound state wave-
functions. For convenience, we define
$|\bra{\psi_{n\ell}}\mathbf{r}\ket{\psi_{n^{\prime}\ell^{\prime}}}|^{2}$ as
the squared matrix element average over $m,m^{\prime}$. We start our
evaluation, by representing both wave functions in terms of the Hypergeometric
functions ${}_{1}F_{1}$.777 Another evaluation can be made by representing the
wave functions of the bound states in terms of the associated Laguerre
polynomials:  $\displaystyle
L_{n}^{\alpha}(x)=\sum_{i=0}^{n}(-1)^{i}\left(\begin{array}[]{c}n+\alpha\\\
n-i\end{array}\right)\frac{x^{i}}{i!}.$ (44) Performing the radial integral in
this representation leads to  $\displaystyle I_{R,{\rm
trans}}=N_{n\ell}(\kappa)N_{n^{\prime}\ell^{\prime}}(\kappa^{\prime})(\tilde{\kappa}+\tilde{\kappa}^{\prime})^{-4-\ell-\ell^{\prime}}\times$
(45)
$\displaystyle\sum_{k=0}^{n-\ell-1}\sum_{k^{\prime}=0}^{n^{\prime}-\ell^{\prime}-1}\left(\begin{array}[]{c}n+\ell\\\
n-\ell-k-1\end{array}\right)\left(\begin{array}[]{c}n^{\prime}+\ell^{\prime}\\\
n^{\prime}-\ell^{\prime}-k^{\prime}-1\end{array}\right)\times$ (50)
$\displaystyle\frac{(-2)^{k+k^{\prime}}\Gamma(4+\ell+\ell^{\prime}+k+k^{\prime})}{k!k^{\prime}!}\frac{\tilde{\kappa}^{k}(\tilde{\kappa}^{\prime})^{k^{\prime}}}{(\tilde{\kappa}+\tilde{\kappa}^{\prime})^{k+k^{\prime}}}\,.$
(51) With the help of Mathematica, both sums can be performed, resulting in
similar recursions as we have obtained for the scattering-to-bound case. In
practice, however, it turns out that these bound-to-bound recursions are less
numerically stable within standard digit precision than what we have obtained
by using the Hypergeometric functions in Eqs. (52),(53) and (55).  As pointed
out recently in Ref. Biondini:2023zcz , one can perform the radial integration
and obtain closed expressions in terms of Hypergeometric functions by
following Ref. GordonZurBD . Adopting this procedure, we obtain for the
general case

$\displaystyle|I_{R,\text{trans}}^{\ell^{\prime}=\ell+1}|=N_{n,\ell}(\tilde{\kappa})N_{n^{\prime},\ell+1}(\tilde{\kappa}^{\prime})\times$
(52) $\displaystyle
2^{-1}(\ell+1)(2\ell+3)\Gamma(2\ell+2)(1-z)^{\frac{n^{\prime}-n}{2}}\frac{z^{\ell+1}}{(\tilde{\kappa}^{\prime})^{2}\tilde{\kappa}}\times$
$\displaystyle\bigg{|}\frac{\tilde{\kappa}^{\prime}-\tilde{\kappa}n+\tilde{\kappa}^{\prime}n^{\prime}}{(\tilde{\kappa}+\tilde{\kappa}^{\prime})^{2}}\prescript{}{2}{F}{{}_{1}}(-n+\ell+1;n^{\prime}+\ell+2;2\ell+2;z)$
$\displaystyle+2\frac{\tilde{\kappa}n-\tilde{\kappa}^{\prime}n^{\prime}}{(\tilde{\kappa}^{2}-\tilde{\kappa}^{\prime
2})}\prescript{}{2}{F}{{}_{1}}(-n+\ell+1;n^{\prime}+\ell+1;2\ell+2;z)$
$\displaystyle-\frac{\tilde{\kappa}^{\prime}+\tilde{\kappa}n-\tilde{\kappa}^{\prime}n^{\prime}}{(\tilde{\kappa}-\tilde{\kappa}^{\prime})^{2}}\prescript{}{2}{F}{{}_{1}}(-n+\ell+1;n^{\prime}+\ell;2\ell+2;z)\bigg{|}\,,$

and

$\displaystyle|I_{R,\text{trans}}^{\ell^{\prime}=\ell-1}|=N_{n,\ell}(\tilde{\kappa})N_{n^{\prime},\ell-1}(\tilde{\kappa}^{\prime})\times$
(53) $\displaystyle
2^{-1}\ell(2\ell+1)\Gamma(2\ell)(1-z)^{\frac{n^{\prime}-n}{2}}\frac{z^{\ell}}{\tilde{\kappa}^{2}\tilde{\kappa}^{\prime}}\times$
$\displaystyle\bigg{|}-\frac{\tilde{\kappa}+\tilde{\kappa}n-\tilde{\kappa}^{\prime}n^{\prime}}{(\tilde{\kappa}-\tilde{\kappa}^{\prime})^{2}}\prescript{}{2}{F}{{}_{1}}(-n+\ell-1;n^{\prime}+\ell;2\ell;z)$
$\displaystyle-2\frac{\tilde{\kappa}n-\tilde{\kappa}^{\prime}n^{\prime}}{(\tilde{\kappa}^{2}-\tilde{\kappa}^{\prime
2})}\prescript{}{2}{F}{{}_{1}}(-n+\ell;n^{\prime}+\ell;2\ell;z)$
$\displaystyle+\frac{\tilde{\kappa}-\tilde{\kappa}n+\tilde{\kappa}^{\prime}n^{\prime}}{(\tilde{\kappa}+\tilde{\kappa}^{\prime})^{2}}\prescript{}{2}{F}{{}_{1}}(-n+\ell+1;n^{\prime}+\ell;2\ell;z)\bigg{|}\,,$

with
$z=\frac{4\tilde{\kappa}\tilde{\kappa}^{\prime}}{(\tilde{\kappa}+\tilde{\kappa}^{\prime})^{2}}<1$
and $\tilde{\kappa}=\alpha\mu/n$,
$\tilde{\kappa}^{\prime}=\alpha^{\prime}\mu^{\prime}/n^{\prime}$ and
$\tilde{\kappa}\neq\tilde{\kappa}^{\prime}$. For the special case
$\alpha^{\prime}\mu^{\prime}=\alpha\mu$ (particularly important for QED), see
Ref. Biondini:2023zcz ; GordonZurBD . The normalization is given by:

$\displaystyle
N_{n,\ell}(\tilde{\kappa})=\frac{\tilde{\kappa}^{3/2}}{\sqrt{n}}\frac{2}{(2\ell+1)!}\left(\frac{(n+\ell)!}{(n-\ell-1)!}\right)^{1/2}\,.$
(54)

For large initial state principle quantum numbers $n=\mathcal{O}(10)$, the
numerical evaluation of the Hypergeometric functions can be problematic for
transitions where $z\rightarrow 1$ ($n^{\prime}\rightarrow 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$:

$\displaystyle\prescript{}{2}{F}{{}_{1}}(a;b;c;z)\rightarrow(-1)^{-a}(-1+z)^{-a-b+c}\times$
(55)
$\displaystyle\frac{\Gamma(1-a)\Gamma(c)}{\Gamma(b)\Gamma(1\\!-\\!a\\!-\\!b+c)}\prescript{}{2}{F}{{}_{1}}(c-a;c-b;c-a-b+1;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\leq 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

###  B.1 Thermal average and Milne relations

In the non-relativistic limit, the thermal average of the bound state
formation cross section can be written as

$\big{\langle}\sigma v\big{\rangle}_{i}=\left(\frac{\mu}{2\pi
T}\right)^{3/2}\\!\int\\!\mathrm{d}^{3}v\,\mathrm{e}^{-\frac{\mu
v^{2}}{2T}}\left[1+f(\Delta E)\right]\,(\sigma v)_{i}\,,$ (56)

where $f(\Delta E)=1/(\mathrm{e}^{\Delta E/T}-1)$. We note that the thermal
averages in non-Abelian gauge theories need to be evaluated carefully due to
oscillatory features.888For example, for effective coupling
$\alpha_{s}^{\text{eff}}<0<\alpha_{b}^{\text{eff}}$ relevant in $SU(N_{c})$,
$(\sigma v)_{n\ell}$ features $n-\ell-1$ local minima in its velocity
dependence. Fewer local minima arise when
$0<\alpha_{s}^{\text{eff}}<\alpha_{b}^{\text{eff}}$.

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 electric dipole interaction is spin-independent, bound-to-
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 $\xi=1/4$ or $\xi=3/4$,
respectively, in Eqs. (61) and (71). Furthermore, the BSF cross sections also
apply for bound states composed of charged scalars, with $\xi=1$. The reason
is that the factor $1/4=1/2\times 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_{j}^{\text{eq}}$ (the
combined particle and anti-particle yield) and
$Y_{\mathcal{B}_{i}}^{\text{eq}}$, the ionization rate can be expressed as

$\Gamma_{\text{ion}}^{i}=\frac{s}{4}\,\frac{(Y_{j}^{\mathrm{eq}})^{2}}{Y_{{\cal
B}_{i}}^{\mathrm{eq}}}\big{\langle}\sigma v\big{\rangle}_{i}\,.$ (57)

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

$\Gamma_{\text{trans}}^{i\to j}=\Gamma_{\text{trans}}^{j\to i}\frac{Y_{{\cal
B}_{j}}^{\mathrm{eq}}}{Y_{{\cal B}_{i}}^{\mathrm{eq}}}\,.$ (58)

###  B.2 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._
$\alpha_{s}^{\text{eff}}=\alpha_{b}^{\text{eff}}=\alpha$ here. The
annihilation cross section into a pair of massless dark photons, including
Sommerfeld enhancement, is given by doi:10.1002/andp.19314030302 ;
Sakharov:1948yq

$\displaystyle(\sigma
v)_{\text{ann}}=\frac{\pi\alpha^{2}}{m^{2}}\,S_{0}\left(\frac{\alpha}{v}\right)\,,$
(59)

where

$S_{0}(\zeta)\equiv\frac{2\pi\zeta}{1-e^{-2\pi\zeta}}\,.$ (60)

The BSF cross section, which includes the summation over the degenerate
magnetic quantum number of the final state $n\ell$ and also considers all
possible initial angular momenta, is in the case of $U(1)$ interactions given
by

$\displaystyle\left(\sigma
v\right)_{n\ell}=\xi\frac{\pi\,\alpha^{2}}{m^{2}}\frac{2^{9}}{3}\,S_{\text{BSF}}(n,\ell,\frac{\alpha}{nv},\frac{\alpha}{v})\,,$
(61)

where $\xi$ accounts for spin factors as defined in Sec. B.1, and

$\displaystyle S_{\text{BSF}}(n,\ell,\tilde{\zeta}_{b},$
$\displaystyle\zeta_{s})\equiv\,\frac{1}{2^{6}n\tilde{\zeta}_{b}}\left(1+{\tilde{\zeta}_{b}}^{2}\right)^{3}\times$
$\displaystyle\left[(\ell+1)|I_{R}^{\ell^{\prime}=\ell+1}|^{2}+\ell|I_{R}^{\ell^{\prime}=\ell-1}|^{2}\right].$
(62)

$S_{\text{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. pirenne1946proper ;
doi:10.1111/j.1749-6632.1946.tb31764.x :

$\displaystyle\Gamma_{\text{dec,\,QED}}^{n\ell\to\gamma\gamma}$
$\displaystyle=\delta_{\ell,0}\,\frac{m\,\alpha^{5}}{2n^{3}}\,.$ (63)

For spin-triplet bound states we adopt the decay rate into three dark photons
for the s-wave states from Ref. 1949PhRv…75.1696O :

$\displaystyle\Gamma_{\text{dec,\,QED}}^{n\ell\to 3\gamma}$
$\displaystyle=\frac{4(\pi^{2}-9)\alpha}{9\pi}\times\Gamma_{\text{dec,\,QED}}^{n\ell\to\gamma\gamma}\,.$
(64)

Lastly, the de-excitation rate relates to the dipole matrix element via

$\displaystyle\Gamma_{\text{trans,\,QED}}^{n^{\prime}\ell^{\prime}\to n\ell}$
$\displaystyle=$
$\displaystyle\frac{4\alpha}{3}(2\ell+1)\left(\frac{m\alpha^{2}}{4}\left|\frac{1}{n^{2}}-\frac{1}{n^{\prime
2}}\right|\right)^{3}$ (65)
$\displaystyle{}\times|\bra{\psi_{n^{\prime}\ell^{\prime}}}\mathbf{r}\ket{\psi_{n\ell}}|^{2},$

according to App. A.2.

###  B.3 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 $\alpha(m)\,\equiv\,0.025$ and
employ one-loop running. Using this choice, the non-perturbative regime
$\alpha(m/x)=1$ starts at $x\approx 4\times 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 (${\bf 1}$) and octet
(${\bf 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

$\displaystyle\alpha_{b}^{\text{eff}}$ $\displaystyle\equiv
C_{F}\,\alpha\left(\frac{m}{2}\frac{\alpha_{b}^{\text{eff}}}{n}\right)\,,$
(66) $\displaystyle\alpha_{s}^{\text{eff}}$
$\displaystyle\equiv\frac{2C_{F}-C_{A}}{2}\,\alpha\left(\frac{m}{2}v\right)\,,$
(67)

where $C_{F}=4/3$ and $C_{A}=3$. Eq. (66) 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

$\displaystyle(\sigma
v)_{\text{ann}}=\frac{7}{27}\frac{\pi\alpha(2m)^{2}}{m^{2}}\left(\frac{2}{7}S_{0}^{[\bf{1}]}+\frac{5}{7}S_{0}^{[\bf{8}]}\right)\,.$
(68)

The Sommerfeld factors are given by

$\displaystyle S^{[\bf{1}]}_{0}$
$\displaystyle=\frac{\alpha(mv/2)}{v}\frac{2\pi C_{F}}{1-e^{-2\pi
C_{F}\alpha(m\frac{v}{2})/v}}\,,$ (69) $\displaystyle S^{[\bf{8}]}_{0}$
$\displaystyle=\frac{\alpha(mv/2)}{v}\frac{2\pi(C_{F}-C_{A}/2)}{1-e^{-2\pi(C_{F}-C_{A}/2)\alpha(mv/2)/v}}\,.$
(70)

The BSF cross section in a general $SU(N_{c})$ theory takes the form
Garny:2021qsr

$\displaystyle\left(\sigma
v\right)_{n\ell}=\xi\frac{\pi\,\alpha_{b}^{\text{eff}}\,\alpha_{\text{BSF}}}{m^{2}}\frac{2^{9}C_{F}}{3N_{c}^{2}}\,S_{\text{BSF}}(n,\ell,\frac{\alpha_{b}^{\text{eff}}}{nv},\frac{\alpha_{s}^{\text{eff}}}{v}),$
(71)

with $\alpha_{\text{BSF}}=\alpha\left(mv^{2}/4+E_{{\cal B}_{n\ell}}\right)$
and $S_{\text{BSF}}$ defined in Eq. (62). Here $\xi$ accounts for spin factors
as defined in Sec. B.1. The spin-singlet bound state decay rate into two
gluons is given by Harz:2018csl ; Garny:2021qsr ; Biondini:2023zcz

$\displaystyle\Gamma_{\text{dec,\,QCD}}^{n\ell\to gg}$
$\displaystyle=\delta_{\ell,0}\,\frac{m\,C_{F}}{4n^{3}}\alpha(m)^{2}\left(\alpha_{b}^{\text{eff}}\right)^{3}\,.$
(72)

We neglect spin-triplet bound states for dark QCD.

Dark QCD automatically corresponds to the limiting case of no transitions,
_i.e._ $\Gamma_{\text{trans}}^{i\to j}=0$, therefore the effective cross
section Eq. (13) simplifies to Binder:2021vfo ; Garny:2021qsr

$\big{\langle}\sigma v\big{\rangle}_{\text{eff}}\to\big{\langle}\sigma
v\big{\rangle}_{\text{ann}}+\sum_{i}\big{\langle}\sigma
v\big{\rangle}_{i}\,\frac{\Gamma^{i}_{\text{dec}}}{\Gamma^{i}_{\text{ion}}+\Gamma^{i}_{\text{dec}}}\,.$
(73)

This result can be seen as a straightforward generalization of the single
bound state case Ellis:2015vaa , extended to a sum over individual bound
states that do not impact each other.

###  B.4 superWIMP scenario

The expressions for BSF of a colored $t$-channel scalar mediator is identical
to that for dark QCD in Eq. (71), now using $m=m_{\tilde{q}}$ and the SM
strong coupling strength $\alpha_{s}$ to 5-loop accuracy as implemented in
RunDec 3 Herren:2017osy  in place of $\alpha$, as well as $\xi=1$. The
annihilation cross section reads

$\displaystyle(\sigma
v)_{\text{ann}}=\frac{14}{27}\frac{\pi\alpha(2m_{\tilde{q}})^{2}}{m_{\tilde{q}}^{2}}\left(\frac{2}{7}S_{0}^{[\bf{1}]}+\frac{5}{7}S_{0}^{[\bf{8}]}\right)\,,$
(74)

and the decay rate is given by Harz:2018csl ; Garny:2021qsr

$\displaystyle\Gamma_{\text{dec}}^{n\ell\to gg}$
$\displaystyle=\delta_{\ell,0}\,\frac{m_{\tilde{q}}\,C_{F}}{8n^{3}}\alpha(m_{\tilde{q}})^{2}\left(\alpha_{b}^{\text{eff}}\right)^{3}\,.$
(75)

Gauge invariance fixes the gauge representations of $\tilde{q}$ to that of the
bottom quark in our model, hence the electromagnetic charge is $Q=-1/3$. The
addition of electromagnetic interactions leads to transitions between bound
states in dipole approximation but to no additional relevant decay or BSF
channels. The transition matrix elements are computed according to App. A.2
and enter the de-excitation rates as

$\displaystyle\Gamma_{\text{trans,\,SW}}^{n^{\prime}\ell^{\prime}\to
n\ell}=\frac{4Q^{2}_{\tilde{q}}\,\alpha_{\text{em}}}{3}(2\ell+1)(\Delta
E_{nn^{\prime}})^{3}\,|\bra{\psi_{n^{\prime}\ell^{\prime}}^{[\bf
1]}}\mathbf{r}\ket{\psi_{n\ell}^{[\bf 1]}}|^{2}\,,$ (76)

where the fine structure constant is $\alpha_{\text{em}}=1/128.9$, and

$\Delta
E_{nn^{\prime}}=\frac{m_{\tilde{q}}}{4}\left|\frac{\alpha_{b}^{\text{eff}}(n)}{n^{2}}-\frac{\alpha_{b}^{\text{eff}}(n^{\prime})}{n^{\prime
2}}\right|\,.$ (77)

##  Appendix C Relic abundance for dark QED

Figure 10:  Coupling strength $\alpha$ versus the dark matter mass $m$ shown
for various approximations, imposing the right relic abundance constraint.
Excited states with $n\leq 100,\ell\leq n-1$ and all dark electric dipole
transitions among them (about a million) shown as the blue solid line, are a
major correction to the original computation of capture into the ground state
only on top of the Sommerfeld effect (blue dashed), _cf._ Ref.
vonHarling:2014kha .

Here, we explore the parameter space of dark QED consistent with the relic
density measurement. We assume that the dark QED sector is in kinetic
equilibrium with the SM heat bath. Solving Eq. (11) for a large number of
points in the two-dimensional parameter space of the model, we compute the
coupling strength $\alpha$ as a function of the dark matter mass $m$ that
provides $\Omega_{\chi}h^{2}\simeq 0.12$ Planck:2018vyg . Figure 10 displays
the respective results under various approximations. While the known cases of
Sommerfeld enhanced annihilation and capture into the ground state further
improve the tree-level result, it is visible from the blue solid line that the
inclusion of all bound states with $n\leq 100$ and $\ell\leq n-1$ (about 5000
in total), and all possible electric dipole transitions among them (about
$10^{6}$ in total) results in corrections that are worth to us to report.

To arrive at this result, we have included spin-singlet and spin-triplet decay
of the s-wave bound states in the effective cross section, such that all
curves other than the blue solid line confirm the earlier result in Ref.
vonHarling:2014kha . We have observed convergence regarding the number of
included bound states and transitions. This completes the discussion within
the electric dipole operator picture, in particular, for ultra-soft processes
with one dark photon.

## References

  * (1) A. Mitridate, M. Redi, J. Smirnov and A. Strumia, _Cosmological Implications of Dark Matter Bound States_ , _JCAP_ 1705 (2017) 006, [1702.01141].
  * (2) S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo et al., _Closing the window on WIMP Dark Matter_ , _Eur. Phys. J. C_ 82 (2022) 31, [2107.09688].
  * (3) S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo et al., _The last complex WIMPs standing_ , _Eur. Phys. J. C_ 82 (2022) 992, [2205.04486].
  * (4) J. Ellis, F. Luo and K. A. Olive, _Gluino Coannihilation Revisited_ , _JHEP_ 09 (2015) 127, [1503.07142].
  * (5) J. Ellis, J. L. Evans, F. Luo and K. A. Olive, _Scenarios for Gluino Coannihilation_ , _JHEP_ 02 (2016) 071, [1510.03498].
  * (6) S. P. Liew and F. Luo, _Effects of QCD bound states on dark matter relic abundance_ , _JHEP_ 02 (2017) 091, [1611.08133].
  * (7) S. Kim and M. Laine, _Rapid thermal co-annihilation through bound states in QCD_ , _JHEP_ 07 (2016) 143, [1602.08105].
  * (8) J. Harz and K. Petraki, _Higgs Enhancement for the Dark Matter Relic Density_ , _Phys. Rev. D_ 97 (2018) 075041, [1711.03552].
  * (9) J. Harz and K. Petraki, _Radiative bound-state formation in unbroken perturbative non-Abelian theories and implications for dark matter_ , _JHEP_ 07 (2018) 096, [1805.01200].
  * (10) S. Biondini and M. Laine, _Thermal dark matter co-annihilating with a strongly interacting scalar_ , _JHEP_ 04 (2018) 072, [1801.05821].
  * (11) S. Biondini and S. Vogl, _Coloured coannihilations: Dark matter phenomenology meets non-relativistic EFTs_ , _JHEP_ 02 (2019) 016, [1811.02581].
  * (12) H. Fukuda, F. Luo and S. Shirai, _How Heavy can Neutralino Dark Matter be?_ , _JHEP_ 04 (2019) 107, [1812.02066].
  * (13) S. Biondini and S. Vogl, _Scalar dark matter coannihilating with a coloured fermion_ , _JHEP_ 11 (2019) 147, [1907.05766].
  * (14) M. Garny and J. Heisig, _Bound-state effects on dark matter coannihilation: Pushing the boundaries of conversion-driven freeze-out_ , _Phys. Rev. D_ 105 (2022) 055004, [2112.01499].
  * (15) M. Becker, E. Copello, J. Harz, K. A. Mohan and D. Sengupta, _Impact of Sommerfeld effect and bound state formation in simplified t-channel dark matter models_ , _JHEP_ 08 (2022) 145, [2203.04326].
  * (16) R. Oncala and K. Petraki, _Dark matter bound state formation via emission of a charged scalar_ , _JHEP_ 02 (2020) 036, [1911.02605].
  * (17) R. Oncala and K. Petraki, _Bound states of WIMP dark matter in Higgs-portal models. Part II. Thermal decoupling_ , _JHEP_ 08 (2021) 069, [2101.08667].
  * (18) B. von Harling and K. Petraki, _Bound-state formation for thermal relic dark matter and unitarity_ , _JCAP_ 12 (2014) 033, [1407.7874].
  * (19) A. Kamada, M. Yamada and T. T. Yanagida, _Unification for darkly charged dark matter_ , _Phys. Rev. D_ 102 (2020) 015012, [1908.00207].
  * (20) T. Binder, B. Blobel, J. Harz and K. Mukaida, _Dark matter bound-state formation at higher order: a non-equilibrium quantum field theory approach_ , _JHEP_ 09 (2020) 086, [2002.07145].
  * (21) S. Biondini, N. Brambilla, G. Qerimi and A. Vairo, _Effective field theories for dark matter pairs in the early universe: cross sections and widths_ , _JHEP_ 07 (2023) 006, [2304.00113].
  * (22) V. De Luca, A. Mitridate, M. Redi, J. Smirnov and A. Strumia, _Colored Dark Matter_ , _Phys. Rev. D_ 97 (2018) 115024, [1801.01135].
  * (23) P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer and J. Smirnov, _Thermal squeezeout of dark matter_ , _Phys. Rev. D_ 104 (2021) 095013, [2103.09827].
  * (24) T. Binder, K. Mukaida, B. Scheihing-Hitschfeld and X. Yao, _Non-Abelian electric field correlator at NLO for dark matter relic abundance and quarkonium transport_ , _JHEP_ 01 (2022) 137, [2107.03945].
  * (25) P. Asadi, E. D. Kramer, E. Kuflik, T. R. Slatyer and J. Smirnov, _Glueballs in a thermal squeezeout model_ , _JHEP_ 07 (2022) 006, [2203.15813].
  * (26) M. Becker, J. Harz, E. Copello, K. A. Mohan and D. Sengupta, _Implications of Nonperturbative Effects for Colored Dark Sectors_ , _LHEP_ 2023 (2023) 363.
  * (27) B. Assi and M. L. Wagman, _Baryons, multi-hadron systems, and composite dark matter in non-relativistic QCD_ , 2305.01685.
  * (28) P. Asadi, M. Baumgart, P. J. Fitzpatrick, E. Krupczak and T. R. Slatyer, _Capture and Decay of Electroweak WIMPonium_ , _JCAP_ 02 (2017) 005, [1610.07617].
  * (29) S. Biondini, P. Schicho and T. V. I. Tenkanen, _Strong electroweak phase transition in t-channel simplified dark matter models_ , _JCAP_ 10 (2022) 044, [2207.12207].
  * (30) M. B. Wise and Y. Zhang, _Stable Bound States of Asymmetric Dark Matter_ , _Phys. Rev. D_ 90 (2014) 055030, [1407.4121].
  * (31) K. Petraki, M. Postma and M. Wiechers, _Dark-matter bound states from Feynman diagrams_ , _JHEP_ 06 (2015) 128, [1505.00109].
  * (32) H. An, M. B. Wise and Y. Zhang, _Strong CMB Constraint On P-Wave Annihilating Dark Matter_ , _Phys. Lett. B_ 773 (2017) 121–124, [1606.02305].
  * (33) S. Biondini and V. Shtabovenko, _Bound-state formation, dissociation and decays of darkonium with potential non-relativistic Yukawa theory for scalar and pseudoscalar mediators_ , _JHEP_ 03 (2022) 172, [2112.10145].
  * (34) S. Biondini and V. Shtabovenko, _Non-relativistic and potential non-relativistic effective field theories for scalar mediators_ , _JHEP_ 08 (2021) 114, [2106.06472].
  * (35) M. Garny, J. Heisig, B. Lülf and S. Vogl, _Coannihilation without chemical equilibrium_ , _Phys. Rev._ D96 (2017) 103521, [1705.09292].
  * (36) T. Binder, A. Filimonova, K. Petraki and G. White, _Saha equilibrium for metastable bound states and dark matter freeze-out_ , _Phys. Lett. B_ 833 (2022) 137323, [2112.00042].
  * (37) L. Covi, J. E. Kim and L. Roszkowski, _Axinos as cold dark matter_ , _Phys. Rev. Lett._ 82 (1999) 4180–4183, [hep-ph/9905212].
  * (38) J. L. Feng, A. Rajaraman and F. Takayama, _SuperWIMP dark matter signals from the early universe_ , _Phys. Rev._ D68 (2003) 063504, [hep-ph/0306024].
  * (39) M. Garny and J. Heisig, _Interplay of super-WIMP and freeze-in production of dark matter_ , _Phys. Rev. D_ 98 (2018) 095031, [1809.10135].
  * (40) Q. Decant, J. Heisig, D. C. Hooper and L. Lopez-Honorez, _Lyman- $\alpha$ constraints on freeze-in and superWIMPs_, _JCAP_ 03 (2022) 041, [2111.09321].
  * (41) J. Bollig and S. Vogl, _Impact of bound states on non-thermal dark matter production_ , _JCAP_ 10 (2022) 031, [2112.01491].
  * (42) A. Pineda and J. Soto, _Effective field theory for ultrasoft momenta in NRQCD and NRQED_ , _Nucl. Phys. Proc. Suppl._ 64 (1998) 428–432, [hep-ph/9707481].
  * (43) M. Beneke, _New results on heavy quarks near threshold_ , in _3rd Workshop on Continuous Advances in QCD (QCD 98)_ , pp. 293–309, 6, 1998.  hep-ph/9806429.
  * (44) M. Beneke, _Perturbative heavy quark - anti-quark systems_ , _PoS_ hf8 (1999) 009, [hep-ph/9911490].
  * (45) N. Brambilla, A. Pineda, J. Soto and A. Vairo, _Potential NRQCD: An Effective theory for heavy quarkonium_ , _Nucl. Phys. B_ 566 (2000) 275, [hep-ph/9907240].
  * (46) S. Biondini, T. A. Chowdhury and S. Saad, _Non-perturbative effects for dark sectors with QCD portals_ , _JHEP_ 08 (2023) 002, [2306.09428].
  * (47) H. A. Bethe and E. E. Salpeter, _Quantum Mechanics of One- and Two-Electron Atoms_.  1957, 10.1007/978-3-662-12869-5.
  * (48) K. Petraki, M. Postma and J. de Vries, _Radiative bound-state-formation cross-sections for dark matter interacting via a Yukawa potential_ , _JHEP_ 04 (2017) 077, [1611.01394].
  * (49) H. A. Kramers, _Xciii. on the theory of x-ray absorption and of the continuous x-ray spectrum_ , _Philosophical Magazine Series 1_ 46 (1923) 836–871.
  * (50) V. M. Katkov and V. M. Strakhovenko, _Radiative recombination and its application in experiments on electron cooling_ , _Soviet Journal of Experimental and Theoretical Physics_ 48 (Oct., 1978) 639.
  * (51) K. M. Belotsky, E. A. Esipova and A. A. Kirillov, _On the classical description of the recombination of dark matter particles with a Coulomb-like interaction_ , _Phys. Lett. B_ 761 (2016) 81–86, [1506.03094].
  * (52) H. An, M. B. Wise and Y. Zhang, _Effects of Bound States on Dark Matter Annihilation_ , _Phys. Rev. D_ 93 (2016) 115020, [1604.01776].
  * (53) J. Smirnov and J. F. Beacom, _TeV-Scale Thermal WIMPs: Unitarity and its Consequences_ , _Phys. Rev. D_ 100 (2019) 043029, [1904.11503].
  * (54) I. Baldes and K. Petraki, _Asymmetric thermal-relic dark matter: Sommerfeld-enhanced freeze-out, annihilation signals and unitarity bounds_ , _JCAP_ 09 (2017) 028, [1703.00478].
  * (55) S. Bottaro and D. Redigolo, _The dark matter unitarity bound at NLO_ , 2305.01680.
  * (56) K. Griest and M. Kamionkowski, _Unitarity Limits on the Mass and Radius of Dark Matter Particles_ , _Phys. Rev. Lett._ 64 (1990) 615.
  * (57) R. J. Scherrer and M. S. Turner, _On the Relic, Cosmic Abundance of Stable Weakly Interacting Massive Particles_ , _Phys. Rev. D_ 33 (1986) 1585.
  * (58) J. Hisano, S. Matsumoto and M. M. Nojiri, _Explosive dark matter annihilation_ , _Phys. Rev. Lett._ 92 (2004) 031303, [hep-ph/0307216].
  * (59) J. Hisano, S. Matsumoto, M. M. Nojiri and O. Saito, _Non-perturbative effect on dark matter annihilation and gamma ray signature from galactic center_ , _Phys. Rev._ D71 (2005) 063528, [hep-ph/0412403].
  * (60) B. Konstantin, E. Ekaterina, D. Kalashnikov and A. Letunov, _Problems of the correspondence principle for the recombination cross section in dark plasma_ , 2011.04685.
  * (61) M. Garny, A. Ibarra and S. Vogl, _Signatures of Majorana dark matter with t-channel mediators_ , _Int. J. Mod. Phys._ D24 (2015) 1530019, [1503.01500].
  * (62) C. Arina, B. Fuks, J. Heisig, M. Krämer, L. Mantani and L. Panizzi, _A comprehensive exploration of t-channel simplified models of dark matter_ , 2307.10367.
  * (63) K. Griest and D. Seckel, _Three exceptions in the calculation of relic abundances_ , _Phys. Rev._ D43 (1991) 3191–3203.
  * (64) J. Edsjö and P. Gondolo, _Neutralino relic density including coannihilations_ , _Phys. Rev._ D56 (1997) 1879–1894, [hep-ph/9704361].
  * (65) R. T. D’Agnolo, D. Pappadopulo and J. T. Ruderman, _Fourth Exception in the Calculation of Relic Abundances_ , _Phys. Rev. Lett._ 119 (2017) 061102, [1705.08450].
  * (66) J. McDonald, _Thermally generated gauge singlet scalars as selfinteracting dark matter_ , _Phys. Rev. Lett._ 88 (2002) 091304, [hep-ph/0106249].
  * (67) T. Asaka, K. Ishiwata and T. Moroi, _Right-handed sneutrino as cold dark matter_ , _Phys. Rev._ D73 (2006) 051301, [hep-ph/0512118].
  * (68) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080, [0911.1120].
  * (69) A. Ibarra, A. Pierce, N. R. Shah and S. Vogl, _Anatomy of Coannihilation with a Scalar Top Partner_ , _Phys. Rev._ D91 (2015) 095018, [1501.03164].
  * (70) F. Herren and M. Steinhauser, _Version 3 of RunDec and CRunDec_ , _Comput. Phys. Commun._ 224 (2018) 333–345, [1703.03751].
  * (71) M. Cacciari, L. Del Debbio, J. R. Espinosa, A. D. Polosa and M. Testa, _A note on the fate of the Landau–Yang theorem in non-Abelian gauge theories_ , _Phys. Lett. B_ 753 (2016) 476–481, [1509.07853].
  * (72) Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * (73) J. Heeck, J. Heisig and A. Thapa, _Testing Dirac leptogenesis with the cosmic microwave background and proton decay_ , _Phys. Rev. D_ 108 (4, 2023) 035014, [2304.09893].
  * (74) W. Gordon, _Zur Berechnung der Matrizen beim Wasserstoffatom_ , _Annalen der Physik_ 394 (1929) 1031–1056.
  * (75) A. Sommerfeld, _Über die Beugung und Bremsung der Elektronen_ , _Annalen der Physik_ 403 (1931) 257–330.
  * (76) A. D. Sakharov, _Interaction of an Electron and Positron in Pair Production_ , _Zh. Eksp. Teor. Fiz._ 18 (1948) 631–635.
  * (77) J. Pirenne, _The proper field and the interaction of dirac particles. i-iii_ , _Arch. Sci. phys. nat_ 28 (1946) 233–272.
  * (78) J. A. Wheeler, _Polyelectrons_ , _Annals of the New York Academy of Sciences_ 48 (1946) 219–238.
  * (79) A. Ore and J. L. Powell, _Three-Photon Annihilation of an Electron-Positron Pair_ , _Physical Review_ 75 (Jun, 1949) 1696–1699.