Dark quarkonium formation in the early universe

The relic abundance of heavy stable particles charged under a confining gauge group can be depleted by a second stage of annihilations near the deconfinement temperature. This proceeds via the formation of quarkonia-like states, in which the heavy pair subsequently annihilates. The size of the quarkonium formation cross section was the subject of some debate. We estimate this cross section in a simple toy model. The dominant process can be viewed as a rearrangement of the heavy and light quarks, leading to a geometric cross section of hadronic size. In contrast, processes in which only the heavy constituents are involved lead to mass-suppressed cross sections. These results apply to any scenario with bound states of sizes much larger than their inverse mass, such as U(1) models with charged particles of different masses, and can be used to construct ultra-heavy dark-matter models with masses above the na\"ive unitarity bound. They are also relevant for the cosmology of any stable colored relic.


Introduction
... and in the darkness bind them.

J. R. R. Tolkien
Stable "colored" particles, charged under QCD or a hidden confining gauge group, have been proposed as dark matter (DM) candidates , and are predicted in various extensions of the Standard Model [23][24][25]. Even in the simplest models, the cosmological history of colored relics is intriguing, and their present-day abundances have been the subject of some debate. The relic abundance of a heavy colored particle X is sensitive to the two inherent scales in the problem: its mass m X , and the confinement scale Λ D . If m X Λ D , the freeze-out of X proceeds via standard perturbative annihilations at temperatures T ∼ m X /30. However, at temperatures T ∼ Λ D , long after the perturbative X-X annihilations have shut off, the X relic abundance may be further reduced by interactions of hadronized Xs, whose size is set by 1/Λ D .
The annihilation process at T Λ D was described in a semiclassical approximation in Ref. [24]. At T ∼ Λ D , most of the Xs are in color-singlet heavy-light hadrons, which we label by H X . An X-hadron H X and anX-hadronH X experience a residual strong interaction whose effective range is ∼ 1/Λ D , much larger than the Compton wavelength ∼ 1/m X , but much smaller than the mean distance between X-hadrons at T ∼ Λ D . X-X annihilation then proceeds via the formation of X-X "quarkonia", which subsequently deexcite to the ground state. In the ground state, the X-X distance is of order the Compton wavelength, and the pair annihilates into light mesons, (dark) photons, or glueballs. In the following, we use parentheses to denote quarkonium-like states, and refer to the (XX) states as quarkonia.
The cross section for quarkonia formation was argued in Ref. [24] to be purely geometric. This certainly holds for the scattering cross section of two H X hadrons; however, it is less clear that it holds for the quarkonium production cross section, which requires a significant modification of the trajectories of the heavy particles. One semiclassical argument for a mass-suppressed cross section was described in Ref. [25]. To form a bound state, the X andX must lose energy and angular momentum. Classically, one can estimate the cross section by modeling the energy loss as Larmor radiation. This is proportional to the acceleration-squared, which scales as Λ 4 D /m 2 X . At T ∼ Λ D , the Xs are very slow, with speed v ∼ Λ D /m X . Thus it takes a long time for the hadrons to cross a distance 1/Λ D ; however, the total amount of radiation still scales as m −3/2 X and is suppressed by the large X mass.
Our goal in this paper is to quantify the cross section for dark quarkonia formation. This is, of course, a strong-coupling problem, so we will employ two simple toy models in which the calculation is tractable. As we will see, the results can be readily interpreted to infer the behavior of the cross section in the case of interest.
We consider a dark SU(N ) with two Dirac fermions X and q in the fundamental representation. We denote the SU(N ) confinement scale by Λ D , although much of our discussion applies to real QCD as well. X is heavy, with mass m X Λ D , while q is light, with m q Λ D . We denote the color-singlet heavy-light mesons by H X ≡ Xq and H X ≡Xq. The X andX, as well as their hadrons, are stable by virtue of a flavor symmetry.
We examine two prototypical contributions to quarkonia production in H X -H X collisions. The first is a radiation process, in which the "brown muck" is merely a spectator. To isolate the contribution of the heavy Xs, we invoke a dark U(1), under which X is charged while q is neutral. The heavy fermions X andX emit radiation in order to bind, and the relevant process is H X +H X → (XX) + ϕ [radiation by the Xs]. (1.1) Here ϕ is the dark photon. Since the photon is emitted by the heavy X, the cross section for this process can be calculated using non-relativistic QCD (NRQCD) with a simple potential modeling the SU(N ) interaction. We use the Cornell potential, with a cutoff at a distance of order 1/Λ D to simulate the screening by the brown muck. The resulting cross section is not geometric, but rather m X -suppressed, in accordance with the simple semiclassical estimate above.
In the second process, the brown muck plays a key role in the interaction, leading to a geometric cross section for quarkonium formation. This happens, for example, when the radiation is emitted by the brown muck itself, which, as a result, exerts a force on the heavy X. While we cannot reliably calculate the cross section in this case, we can nonetheless capture the brown muck dynamics by considering the limit m q > Λ D , in which quarkonium formation can be thought of as a rearrangement of the heavy and light quarks, The cross section for this process can be calculated in analogy with hydrogen-antihydrogen rearrangement into protonium and positronium. As we will see, for m q > Λ D , only the Coulombic states contribute. The result is a geometric cross section, which scales as the square of the Bohr radius a q = 1/(ᾱ D m q ). Thus, quarkonium production is effective at low temperatures not because of confinement per se, but because of the large hierarchy between a q (the size of H X ) and 1/m X (the Compton wavelength). We expect this result to persist as m q is dialed back below Λ D : the quarkonium cross section will continue to scale as the size of H X , which, in this case, is 1/Λ 2 D . As we will see, the geometric cross section arises from summing the contributions of many large (i.e., ∼ a q -sized) X-X bound states, for which the process is exothermic. These states cannot be dissociated, and will de-excite to the ground state, in which the X andX annihilate. We will not discuss the cosmology of a specific model in detail, but merely sketch the essentials, following Ref. [24]. Prior to the formation of the H X andH X mesons, the X particles annihilate and freeze out in the early universe with the standard relic density Following the second stage of annihilations, the H X relic abundance is given by (1.4) Some fraction of Xs remain in hadrons containing multiple Xs, such as baryons. The various final abundances are model dependent and we do not explore them in detail here. Still, the late re-annihilations give a new mechanism for generating the relic abundance of dark matter, which is now a function of the two scales m X and Λ D . This opens up many interesting directions to explore. We discuss some of the implications for cosmology in Sec. 5. In particular, the models can lead to a long era of matter domination between m X and Λ D . Note that the Xs can hadronize with light quarks q, and the potential between them is screened at large distances. The cosmology of these models is thus somewhat different from quirky models. The presence of light quarks is important for yet another reason: even if there are no photons in the theory, energy loss can proceed via the emission of light pions. In contrast, in models with a pure SU(N ) at low energies, the lightest particles are glueballs, whose mass is ∼ 7Λ D .
The formalism and the results in this paper can be applied more broadly. For example, it is applicable to any confined heavy relic-be it all the dark matter or a component thereof, such as gluinos in split supersymmetry, messengers in gauge-mediated supersymmetry breaking, and so on. This paper is structured as follows. The toy model is described in Sec. 2. The rearrangement calculation and its results are presented in Sec. 3. Section 4 focuses on the radiation process from the Xs in which the brown muck is a spectator. In Sec. 5 we consider the dynamics of the X-X bound states generated by these processes and further implications for cosmology. In the Appendix, we collect some useful results on the properties of the Cornell and linear potentials and their wavefunctions, and discuss the details of the derivation of the cross section used in Sec. 4.

Description of the Toy Model
The minimal particle content in our models consists of two Dirac fermions, (q,q) and (X,X), in the fundamental representation of a dark SU(N ). In Sec. 4, we will assume that X andX are also charged under a U(1) gauge symmetry. To describe the X-X interaction, we turn to models of quarkonium [26,27]. The Cornell potential interpolates between the Coulombic QCD potential at small distances and the confining linear potential with string tension Λ D at large distances: where R is the distance between X andX, C = (C 1 + C 2 − C 12 )/2, and C i (C 12 ) are the quadratic Casimirs of the constituents (bound state). Since X is a fundamental of SU(N ) and we require a color-singlet bound state, C = C 1 = C 2 = (N 2 − 1)/(2N ). The deep bound states of the system are then Coulombic, while the shallow states are controlled by the linear potential. At large distances, the attractive potential is screened by the brown muck surrounding X andX. In QCD, for example, this distance is roughly the inverse of the string tension ≈ 400 MeV (see, e.g., Refs. [28,29]). In order to capture this screening, the potential is cut off at a distance R c , 1 whereᾱ D = Cα D , and V 0 is a constant. 2 The cutoff behavior will naturally emerge in the rearrangement calculation of Sec. 3, where we work in the calculable limit m q Λ D . For m q Λ D , the attractive potential is cut off at distances of order the Bohr radius of the heavy-light meson, 3 1 While this option is not pursued here, it may be interesting to use a temperature-dependent cutoff to qualitatively capture the screening effects of the quark-gluon plasma. These cause large X-X bound states to "dissolve" at finite temperatures [29][30][31]. 2 The choice of the constant V0 is of course a matter of convenience, and we will in fact choose different constants in the rearrangement and radiation calculations. 3 In this expression,ᾱD should be evaluated at the energy scale of the inverse Bohr radius.
Thus, for m q sufficiently large, the problem reduces to a purely Coulombic potential, and we can calculate the cross section in analogy with hydrogen-antihydrogen rearrangement into protonium and positronium. In Sec. 4 we will calculate quarkonia formation via radiation by the Xs. Here the linear part of the potential is important, and the cutoff R c is introduced by hand. As we will describe, the choice of R c will be motivated by a comparison with the masses of B and D mesons in the Standard Model.
We note that, in QCD, the string tension in the confined phase and the dimensional transmutation scale from the running of the QCD gauge coupling are approximately the same [29,32,33], whereas the deconfinement temperature is about a factor of two lower [34]. We will not be concerned with the lightest glueball state because its mass is a factor of about seven larger than the string tension [35].

The Rearrangement Process
At temperatures below Λ D , the heavy Xs are mostly found in H X (Xq) andH X (Xq) mesons. These mesons can further deplete through H X -H X scattering into (XX) quarkonia plus light hadrons. For m q < Λ D , the calculation of the cross section for this process requires the full machinery of perturbative NRQCD [36] and is extremely difficult; we will limit ourselves to the case m q Λ D . This puts us firmly in the non-relativistic limit, in which quarkonium production can be thought of as rearrangement of the four partons, For m X m q , the wavefunctions of the system can be calculated in the Born-Oppenheimer approximation, as in hydrogen-antihydrogen rearrangement into protonium and positronium [37][38][39]. We will closely follow this calculation, applying it to the near-threshold energies of interest.
If the semiclassical arguments in Ref. [24] are correct, the cross section is expected to be geometric when the temperature is comparable to the binding energy of H X , with no m X suppression. We verify this in the following calculation.

Setup
As discussed above, for m q sufficiently larger than Λ D , only the Coulombic (XX) states contribute. We will later comment on the validity of this approach as m q is taken below Λ D .
The full interacting Hamiltonian of our system is written as the sum where 3) Figure 1: Coordinate system used in the calculation of the rearrangement process.
Here R is the vector fromX to X and r q , rq are the positions of q,q relative to the X-X center-of-mass (CM), respectively, as shown in Fig. 1. The potentials V qX , Vq X , V qX , VqX , V qq , and V XX are the usual Coulomb potentials (with the relevant sign for same/opposite color quarks taken into account in Eq. (3.3)): Since we assume that X-X are in a color-singlet configuration, this factor is the same for the six potentials. The calculation of the rearrangement cross section involves a subtlety well known to nuclear physicists: the asymptotic in and out states are not eigenstates of the same free Hamiltonian, but rather eigenstates of two different interacting Hamiltonians. This is different from conventional non-relativistic scattering where lim t→±∞ H tot = H free . In our case, the infinite past Hamiltonian is while the infinite future Hamiltonian is The scattering cross section is then calculated in the multi-channel formalism. By solving the Lippmann-Schwinger equation for multi-channel scattering (see, e.g., Ref. [40]), we get the simple formula for the cross section: where k f and k i are the momenta of the final and initial states in the CM frame (see below) and the transition matrix element is where Ψ f , Ψ i are the final-and initial-state wavefunctions and H tr = H tot − H out , as can be seen from Eqs. (3.3) and (3.6). Note that in this representation of the cross section, the outgoing states Ψ f are eigenstates of H out , while the incoming states Ψ i are eigenstates of the full Hamiltonian H tot . Below, we discuss these states in more detail.

The incoming and outgoing wavefunctions
We wish to express our incoming and outgoing wavefunctions in the factorized form In the final state this factorization is exact, since the outgoing X-onium and q-onium are asymptotically non-interacting. In other words: with The final state therefore trivially factorizes as the product of a plane wave for the outgoing (qq) and the Coulomb bound states ψ XX nlm (R) (an eigenstate of H XX ) and ψ qq 100 (r q , rq) (an eigenstate of H qq ). For concreteness, we assume that the final-state q-onium is in its ground state and the (XX) is static.
In contrast with the outgoing state, the incoming state is an eigenstate of the full Hamiltonian H tot , so we naïvely do not expect it to factorize as in Eq. (3.9). However, we can use the Born-Oppenheimer approximation to express it in this factorized form, Since we are in the limit m X m q , this is a very good approximation: at any given X-X distance R, q andq will quickly adjust their configuration, and their wavefunction ψ qq i can therefore be calculated by integrating out X andX and treating them as sources for the light quarks. This gives the energy and wavefunction of the light quarks for a fixed separation R between the heavy Xs as solutions to the eigenvalue problem Substituting this back into the full Schrödinger equation and neglecting derivatives of ψ qq i with respect to the X-X coordinates, one obtains the equation for the X-X wavefunction, (3.14) The effective potential for the X-X system is then The incoming effective potential V in in units ofᾱ 2 D m q for the X-X system in the Born-Oppenheimer approximation (blue solid), as a function of X-X separation in units of the Bohr radius a q = 1/(ᾱ D m q ). Also shown is the Coulomb potential for the (XX) quarkonium (red dashed).
V BO (R) should interpolate between twice the binding energy of H X , 2E b ≡ᾱ 2 D m q , at large R, and the q-onium binding energy, −ᾱ 2 D m q /4, at small R. Unlike in molecules, for which the Coulombic repulsion of the nuclei must be overcome, here the two heavy particles attract each other, so we do not expect a significant potential barrier. These naïve expectations are borne out in the calculation of Ref. [41]. Since V BO (R) does not depend on the initial energy of the system or on the mass m X , we can extract V BO (R) from Ref. [41]. We plot V in in Fig. 2: as expected, the effects of the light quarks captured in V BO (R) set in for R of order a q . Their main effect is to screen the X-X interaction at large R; in practice, this happens for R ∼ 2a q .
Since V in (R) approaches a constant at large R, the X-X wavefunction at large distances (R ≥ 4a q ) is the standard free-particle solution, (3.17) The wavefunctions for R ≤ 4a q are found numerically, while their normalization is fixed by matching to Eq. (3.17) at R = 4a q . Some examples for the incoming and outgoing wavefunctions are given in Fig. 3.

The matrix element for rearrangement
Using the factorized incoming and outgoing wavefunctions, the transition matrix element Eq. (3.8) can be written in position space as where We will assume that the angular part of T (R) is constant. This is justified when the (qq) is in the ground state, and in the short-distance approximation for the plane wave of the (qq) relative to the (XX). The second condition is broken when k 2 f becomes large enough, where we expect an O(1) correction. We neglect this correction in this work, since we are mostly interested in the parametric behavior of this process.
It is easy to see that T (R) is appreciable only for R of order a q . For R a q , this can be seen by substituting H tr = V BO (R) − H qq in Eq. (3.19). ψ qq f is an eigenfunction of H qq , so one is left with the overlap integral of ψ qq f and the asymptotic ψ qq i at large R. These solve the same Hamiltonian with different eigenvalues: the former is a bound state and the latter is a continuum state. For small R, the initial-and final-state wavefunctions in the integrand do overlap, but H tr tends to zero.
The full calculation of T (R) is complicated. In fact, the Born-Oppenheimer approximation breaks down for R 0.74a q , where the q andq are no longer bound to their respective X (see, e.g., Ref. [39]). Still, we can use this approximation to get a rough estimate of the cross section. In particular, T (R) is independent of the mass m X in the Born-Oppenheimer approximation, so we can extract T (R) from Ref. [39]. 4 The result in the relevant range (R ∼ a q ) can be parametrized as , since

Rearrangement results
We calculate the cross section of Eq. (3.7) for different masses m X and incoming kinetic energies E i , keeping a q fixed. In the approximation we use (see Sec. 3.3), the angular part of the overlap integral is translated into a selection rule l XX = l i ≡ l. The breakdown of the cross section into partial waves-or (XX) angular momenta-is given in Fig. 4 for high incoming kinetic energy E i = 0.6E b andᾱ D = 1/137. We see that it vanishes above some maximal l, which corresponds to the classical angular momentum It is also interesting to examine the branching fraction σ nl /σ l for some initial partial wave to form an (XX) quarkonium of definite binding energy. In Fig. 5 (left panel), we show this branching fraction for l = 0, 14, as a function of the final kinetic energy for two values of the initial kinetic energy E i /E b = 0.6, 0.06.
We see that quarkonium formation by rearrangement is an exothermic process: the kinetic energy in the final state does not vanish even when the initial momentum is taken   to zero. Therefore the inverse process shuts off at low energies E b . Furthermore, only quarkonia with binding energies around E b are produced at T ∼ E b . The cross sections drop to zero for large final-state energies corresponding to (XX) binding energies above E b . Thus, deep Coulombic bound states with binding energies E XX b ∼ᾱ 2 D m X are not formed in the rearrangement process. Correspondingly, the bound states produced are large, with size ∼ a q . This behavior is clearly exhibited in Fig. 5 (right panel), where we show the cross section as a function of the mean quarkonium radius R XX .
The total cross section for an initial energy of order the H X binding energy is plotted in Fig. 6 (blue line) as a function of m X . Indeed, the cross section is geometric, σ ∼ a 2 q , and is independent of m X to a very good approximation. It is interesting to compare the partial-wave contributions with the unitarity bound, We therefore also plot in this figure several individual partial-wave contributions normalized by 2l + 1 (green lines), compared to 4π/k 2 i (red dashed line). Clearly, for initial-state energies close to the binding energy of H X , the cross section for each partial wave lies close to the unitarity bound. Summing over all the partial waves up to l max , Thus the total cross section is geometric, and scales with the area of the H X bound state. In the low-temperature limit, we have checked that only s-wave processes are non-vanishing and the cross section scales as as demonstrated in Refs. [37,39]. The above results apply to pure U(1) models. In the context of an SU(N ), we have explicitly seen that for m q above Λ D , the light quarks truncate the X-X attraction for R a q via V BO , long before the linear potential sets in. This justifies neglecting the linear potential in the rearrangement calculation.
We can now turn to the limit of interest, m q below Λ D . As we have seen above, the cross section scales with the size of the H X bound state, a q 1/m X , thanks to the large number of partial waves contributing. This behavior is not special to the purely Coulombic case. In fact, the Coulombic contribution gives a conservative estimate of the cross section generated by the Cornell potential. Thus we expect a geometric cross section for m q below Λ D , with the Bohr radius a q replaced by 1/Λ D . The X-X bound states produced via rearrangement at T ∼ E b are of size ∼ a q , much larger than the Compton wavelength of X. However, since the process is exothermic, these X-X bound states cannot be dissociated at T E b . In the case at hand, since there are light pions (or mesons) in the theory, nothing impedes the relaxation of these states to the ground state, in which the X-X pair annihilates.

The Radiation Process: Spectator Brown Muck
In the rearrangement process described above, the brown muck plays a central role. It is instructive to contrast this with a process in which the brown muck is merely a spectator. As we will see, in this case, the cross section scales with m X .
To isolate the dynamics of the heavy quarks, we take X to be charged under a dark U(1), while the light quarks are neutral. The relevant bound states are large, of size ∼ 1/Λ D , and are described by the linear part of the Cornell potential. To calculate (XX) quarkonium production for T Λ D , we can therefore neglect the Coulombic part of the potential. This is also consistent with previous studies showing that, for pure U(1) models with no light charged particles, (XX) bound-state formation gives only mild modifications of the X relic abundance [42,43]. For bound-state formation to deplete the X abundance by orders of magnitude, a new scale is required. In the radiation process considered here, the new scale is Λ D .
The radiative quarkonium production process is then where ϕ is a (dark) photon that couples only to X; however, our results below apply more generally to other light particles which can be emitted by X. Unlike in the previous section, here the light quarks q are relativistic. We will follow the field-theoretic formalism for computing bound-state formation cross sections with long-range interactions detailed in Ref. [44]. Alternatively, these results can be obtained using the standard non-relativistic QM approach for calculating transition amplitudes, treating the photon as a classical field [45,46]. The first step is to derive the spectrum and two-particle wavefunctions that describe the bound and scattering states of the X-X system. While the light quarks do not actively participate in the radiation process, they screen the heavy Xs at large distances. This is captured by the cutoff R c in the potential of Eq. (2.2), and leads to a continuum of H X -H X states with energies above the open H X -H X -production threshold. Roughly speaking, the hadron mass is given by the sum of the heavy constituent masses, with each light quark or gluon contributing about Λ D to the mass. More precisely, where κ Λ D is an O(1) constant [47]. The spectrum of bottom and charm mesons in QCD suggests κ Λ D Λ D ∼ 600 MeV, with Λ D ∼ 400 MeV, so in the following, we set κ Λ D = 1.5.
To estimate the cutoff R c , we use the fact that the maximal H X binding energy, E max b , coincides with the onset of the continuum, Since the maximal bound state energy of the linear potential is we set the cutoff to In summary, the potential we consider is with R c given by Eq. (4.5) and V 0 of Eq. (2.2) chosen as zero for convenience. Defining which are the characteristic splittings in energy and size between successive states, the radial part of the wavefunction, χ ln , solves Using the semiclassical approximation, we can estimate the maximal angular momentum l max of the bound states and the energy min of the lowest bound state with a given l. The lowest energy bound state for each l classically corresponds to a minimum of the effective potential; the position x l min and the energy V l eff (x l min ) of the minimum must satisfy x l min < x c and V l eff (x l min ) < 0, which result in In Appendix A.1, we collect some results for the effective potential and radial wavefunctions for various choices of the parameters.

Radiation results
The cross section for H X +H X → (XX) lmn + ϕ in the CM frame of the initial state is given by where P ϕ is the three-momentum of the radiated light state and, assuming that ϕ is massless, where E k is the kinetic energy of the initial state. We calculate this cross section in Appendix A. It is useful to write the cross section in terms of dimensionless quantities as where e X is the U(1) charge of X, k = E k /E 0 , J k,ln = ( k − ln ) 3 (l + 1) |I k,l+1→ln | 2 + l |I k,l−1→ln | 2 , (4.14) and I is the radial wavefunction overlap integral We plot J k,ln for several values of l in Fig. 7, fixing m X /Λ D = 125 and the initial kinetic energy E k /Λ D = 0.05. As expected, the large, shallow bound states give the largest contributions.
The total thermally-averaged cross section v rel σ (see Appendix A.3) is shown as a function of the temperature in Fig. 8 (left panel) for several choices of m X . We also show σ for the same parameters (right panel). The cross section is clearly dependent on m X , and decreases as X becomes heavier. In fact, for T 0.1Λ D , the scaling is well described by v rel σ ∝ m −2 X and σ ∝ (m 3 X Λ D ) −1/2 , which agrees with the semiclassical estimate 10 -2 0.05 0.1 0.5 1 5 10 10 -5 10 -4 10 -3 in Sec. 1. Thus, for the high mass region of interest, m X Λ D , this contribution is negligible compared to processes mediated by the brown muck, such as the rearrangement process. We expect this qualitative behavior to persist regardless of the spin of the radiated particles.

Implications for Cosmology
We have found that, at temperatures below the confinement scale Λ D , (XX) bound states are formed with a geometric cross section with no m X suppression. These bound states are of size ∼ 1/Λ D , but since the process is exothermic, they cannot be dissociated, and eventually de-excite to the ground state, in which the X-X pair annihilates. The rate for this de-excitation process depends mainly on the light degrees of freedom. For the large (XX) bound states produced, the level splittings are of order (Λ D /m X ) 1/3 Λ D , so we need massless photons or light pions in order to have allowed transitions (as in the models we considered here). In the case where the DM is charged under real QCD, the rate should be sizable, so we expect any (XX) bound states to quickly decay to light particles. 5 As a result, the abundance of H X hadrons is depleted by this second stage of annihilations, down to [24] which can be much smaller than the relic density from perturbative X-X annihilations earlier in the thermal history. There are, however, various different hadronic states (in addition to H X ) in which X can survive [1]. These include other types of single-X hadrons such as baryonic Xq N −1 hadrons, double-X states such as baryonic XXq N −2 , and up to purely heavy baryons X N . Thus, in general, our simple toy model can produce multi-component DM with different masses and relic abundances.
We leave a detailed investigation of the parameter space of the models for future study, but note a few qualitative features. 6 The cross section for producing double-X states should be comparable to the quarkonium cross section we calculated, albeit smaller by O(1) factors because of the smaller binding energies in this case. However, for N > 2, the double-X hadrons contain some light quark(s), so their size is 1/Λ D . They can therefore interact quite efficiently with hadrons containingX to form bound states containing both X and X, where X-X pairs can quickly annihilate. Meanwhile, they can also interact with H X hadrons and form triple-X hadrons, and this chain may go on until pure-X (or pure-X) baryons are formed. All these processes should have geometric cross sections. The pure-X baryons eventually de-excite to the ground state, the size of which is much smaller than 1/Λ D . Their interactions with other hadrons therefore shut off. As a result, the Xs inside the baryons remain as stable relics, while the other Xs may effectively be annihilated. A systematic analysis necessitates solving the Boltzmann equation with these dynamics, which is beyond the scope of this work. However, it is safe to assume that a sizable fraction of the original Xs can be preserved in a variety of stable hadronic relics.
For DM charged under real QCD, heavy-light hadrons, with size 1/Λ D , are subject to stringent constraints, but they are efficiently depleted at T Λ D as we have shown. In contrast, the relic abundance of X 3 baryons may be just somewhat smaller than the original X relic abundance. Thus the fundamental Xs we considered here may give a different realization of colored DM [21], but whether the models are viable requires a more detailed analysis. Note that the scenario considered in Ref. [21], namely a heavy Dirac adjoint X, is non-generic, in that two Xs can form a stable singlet with no additional light quarks.
In summary, the simple models considered here typically give rise to several components of DM composites of different masses. For certain choices of m X and Λ D , these can exhibit self-interactions, transitions between different excited states, and, depending on the coupling to the Standard Model, modified direct and indirect detection cross sections.
In some variants of these models, the DM abundance can be significantly depleted at Λ D , leading to a long era of matter domination between m X and Λ D .

Conclusions
In this paper we have considered the cosmological dynamics of bound states that are much larger than their inverse mass, taking as an example Xq mesons in a confining theory where X is much heavier than q and the confinement scale. We calculated the cross section for quarkonium production from heavy-light meson scattering. The cross section is geometric, and scales with the area of the incident heavy-light mesons. The relic density of heavy-light X-mesons is therefore efficiently diluted with rates much higher than the s-wave unitarity bound due to the participation of many partial waves in the process. We also find that the process is mainly mediated by the effective interaction of the light quarks. In contrast, processes in which only the heavy constituents participate have mass-suppressed rates. It is amusing to note that if the lifetime of B-mesons were longer, such processes could be experimentally measured.
where the dimensionless coordinate x is defined below Eq. (4.8). The scattering-state wavefunctions can be expanded as For k alongẑ this simplifies to The radial wavefunctions χ ln (x) and χ kl (x) solve Note that this dimensionless problem has only two parameters: x c and the effective Coulomb strength a D ≡ (m X /Λ D ) 2/3ᾱ D . The radial wavefunctions are zero at the origin, χ ln (kl) (0) = 0, and satisfy the normalization conditions so that ψ lm,n |ψ l m ,n = δ ll δ mm δ nn and φ k |φ k = (2π) 3 δ (3) (k − k ). Figure 9 shows V eff (x) in Eq. (A.5) for x c = 15 with a D = 7 (left) and a D = 0 (right). Because of the cutoff, the angular momentum quantum number l of bound states has an upper bound l max given in Eq. (4.10). This is confirmed in the numerical results. 7 The bound-state energy levels are shown in Fig. 10 for a cutoff Cornell (left) and linear (right) potentials. In Ref. [48], the semiclassical approximation is used to obtain the energy 7 In fact, as one can see in Fig. 9, the effective potential for very large l may have a minimum with V eff (x0) ≥ 0, and it produces wavefunctions with > 0 which are mostly confined to the region x < xc.
or in our notation, since We reproduce this result in Fig. 10b. For a Cornell potential (Fig. 10a), deep states are governed by the Coulomb force and therefore obey the well-known Coulombic energy levels where the principal quantum number N is given by n+l. For shallower states, the Coulomb force is negligible and the linear potential governs the spectrum. 8 We show some examples of bound-state and scattering-state wavefunctions that solve Eq. (A.4) with the cutoff linear potential (Eq. (A.5) with a D = 0) in the top and bottom rows of Fig. 11, respectively. The bound state with (l, n) = (5, 10) is found to be the shallowest bound state, but it should be emphasized that this is an accident due to the cutoff being just above the energy of this state. In general, shallower states with smaller l have wavefunctions that tend to penetrate beyond x > x c . For the scattering states, wavefunctions of states with larger and smaller l have penetrate further into the region x < x c .

A.2 Bound-state formation cross section in the dipole approximation
We calculate the matrix element for the bound-state formation by a vector-mediated interaction, H X +H X → (XX) lm,n + ϕ (A. 10) with the X-ϕ interaction given by We follow the approach and notation in Refs. [44,49]. From now on, we focus on the linear potential and set a D = 0 because, as we will see, the radiative bound-state formation process favors shallow bound states, for which the Coulomb force is negligible. The bound-state formation cross section is given by [44] (v rel σ) BSF | · M k→lm,n | 2 (A. 12) in the CM frame, where k, P ϕ , and are the relative momentum of the initial state, the momentum of the radiated light state ϕ, and its polarization vector, respectively, and v rel |k| m X /2 , E k is kinetic energy of the initial scattering state. The matrix element is whereφ k (p) andψ lm,n (p) are the momentum-space equivalents of Eqs. (A.1)-(A.2). It is convenient to expand the matrix element in partial waves (cf. Eq. (A.2)): After partial integration, the derivative acts on ψ * lm,n (r), and the integral over d 3 r is trivial using the delta function. Dotting with the polarization vector and performing another partial integration (the derivative of the cosine is zero as · P ϕ = 0) yields · M k,l →lm,n = 4ie X √ m X d 3 r cos P ϕ 2 · r · ψ * lm,n ∇φ k,l − φ k,l ∇ψ * lm,n , (A. 15) where ψ * lm,n and φ k,l are functions of r. The quantity in brackets above appears in the difference of the Schrödinger equations for the scattering-and bound-state wavefunctions, − 1 m X ∇ · (ψ * lm,n ∇φ k,l − φ k,l ∇ψ * lm,n ) = (E k − E nl )ψ * lm,n φ k,l . (A.16) Using the identity 9 d 3 r (∇ · F ) r = − d 3 r F (∇ · r) = − d 3 r F (A. 17) with F = ψ * lm,n ∇φ k,l − φ k,l ∇ψ * lm,n and substituting into Eq. (A.15), we obtain 10 · M k,l →lm,n = 4ie X √ m X d 3 r ( · r) m X (E k − E ln )ψ * lm,n φ k,l cos P ϕ 2 · r + P ϕ 2 · (ψ * lm,n ∇φ k,l − φ k,l ∇ψ * lm,n ) sin P ϕ 2 · r .
(A. 18) We are interested in temperatures T Λ D for which P ϕ m X , and thus · M k,l →lm,n 4ie X m 3 X (E k − E ln ) d 3 r ( · r)ψ * lm,n φ k,l cos Also, we can evaluate P ϕ · r/2 as Note that as the integrand contains a bound-state wavefunction, the integral has support only for r R c . Also, for T Λ D , the overlaps of the bound and scattering states are larger for shallower bound states. Combining these, we can approximate P ϕ · r/2 1 in order to employ the dipole approximation cos(P ϕ · r/2) → 1, which simplifies the cross section to with J k,ln defined in Eq. (4.14).

A.3 Thermally-averaged cross section
Here we briefly describe the procedure to calculate the thermally-averaged cross section, v rel σ (β), as a function of the inverse temperature β = T −1 . As in the previous discussion, we denote the momenta of the initial particles as k 1 , k 2 .
As we are interested in T m X , the kinetic distributions of H X andH X are given by the Maxwell-Boltzmann distribution, which is normalized so that d 3 p (2π) 3 f MB (p) = 1. With this distribution, the thermallyaveraged cross section is given by