Production of dark-matter bound states in the early universe by three-body recombination

The small-scale structure problems of the universe can be solved by self-interacting dark matter that becomes strongly interacting at low energy. A particularly predictive model for the self-interactions is resonant short-range interactions with an S-wave scattering length that is much larger than the range. The velocity dependence of the cross section in such a model provides an excellent fit to self-interaction cross sections inferred from dark-matter halos of galaxies and clusters of galaxies if the dark-matter mass is about 19 GeV and the scattering length is about 17 fm. Such a model makes definite predictions for the few-body physics of weakly bound clusters of the dark-matter particles. The formation of the two-body bound cluster is a bottleneck for the formation of larger bound clusters. We calculate the production of two-body bound clusters by three-body recombination in the early universe under the assumption that the dark matter particles are identical bosons, which is the most favorable case. If the dark-matter mass is 19 GeV and the scattering length is 17 fm, the fraction of dark matter in the form of two-body bound clusters can increase by as much as 4 orders of magnitude when the dark-matter temperature falls below the binding energy, but its present value remains less than 10−6. The present fraction can be increased to as large as 10−3 by relaxing the constraints from small-scale structure and decreasing the mass of the dark matter particle.


Introduction
The simplest paradigm for dark matter is that it consists of weakly interacting elementary particles. However visible matter consists not only of elementary particles, such as electrons, but also of composite particles, such as nuclei. Nuclei are clusters of protons and neutrons bound by residual forces from QCD. Protons and neutrons consist of quarks bound by the color force of QCD. Dark matter could also consist of composite particles. In particular, it could consist of "dark nucleons" and bound clusters of dark nucleons ("dark nuclei"). The dark nucleons could be elementary or composite, but they have an integer-valued conserved charge that we call "dark baryon number". Various models for dark-matter bound states have been discussed in the literature, including a near-threshold S-wave resonance [1][2][3], the exchange of a light mediating boson between elementary fermions, QCD-like structure in the dark sector, and other mechanisms . Light nuclei up to 7 Li are produced in the early universe by big-bang nucleosynthesis [35,36]. The relevant reactions are all 2body collisions of nuclei. Some of the 2-body reactions that produce a nucleus with larger baryon number than either of the colliding nuclei are rearrangement reactions, such as JHEP11(2018)084 d + d → 3 He + n. However the most important such reactions are radiative fusion reactions, such as p + d → 3 He + γ, in which the two incoming nuclei coalesce while radiating a photon to conserve energy and momentum. The effects of 3-body collisions are negligible in big-bang nucleosynthesis. Three-body collisions do play a role in stellar nucleosynthesis despite the relatively low density. In particular, the Hoyle reaction α + α + α → 12 C + γ provides a pathway around the bottleneck caused by the relatively large binding energy of the 4 He nucleus α.
If dark matter consists of dark nucleons that can form bound clusters, these dark nuclei can be produced in the early universe by "dark nucleosynthesis". Studies have shown that a sequence of dark nuclei with increasing dark baryon number can indeed be produced in the early universe [9,11,14,37]. The relevant few-body mechanisms were assumed to be 2-body radiative fusion reactions, in which two incoming dark nuclei coalesce while radiating a much lighter particle to conserve energy and momentum. If there is no such light particle, dark nuclei with larger dark baryon numbers must instead be built up through rearrangement collisions. If dark nuclei with dark baryon number 2 ("dark deuterons") have already been formed, dark nuclei with larger dark baryon numbers can be produced by rearrangement collisions of two dark nuclei, in which dark nucleons are transferred between the two colliding nuclei. However the production of the dark deuterons is a bottleneck that can only be overcome by collisions of 3 or more dark nucleons. The simplest such reaction is the 3-body recombination of three dark nucleons into a dark deuteron and a recoiling dark nucleon. Whether a significant population of dark deuterons can be produced in the early universe can only be determined by detailed calculations using specific models for the few-body physics.
One class of models for few-body physics that are extremely predictive are those with short-range interactions and an S-wave resonance very close to the scattering threshold for a pair of particles [38]. In these models, the elastic scattering cross section for a pair of particles has dramatic energy dependence that is completely determined by the particle mass m and the S-wave scattering length a, which is much larger that the range r 0 of the interactions. When the center-of-mass collision energy E decreases below the energy scale 1/mr 2 0 set by the range, the elastic cross section increases as 1/E, nearly saturating the Swave unitarity bound. The cross section levels off when E decreases below the energy scale 1/ma 2 , approaching a large value proportional to a 2 as E approaches 0. If a is positive, the S-wave resonance is a stable bound cluster. This weakly bound cluster is universal, in the sense that it has properties determined by a, including a small binding energy 1/ma 2 and large geometric size of order a. The universality in the two-particle sector for particles with a large scattering length extends to the 3-particle and higher sectors, although it can be more intricate [38]. It strongly constrains the universal bound clusters, whose binding energies are smaller than the energy scale set by the range. Whether there are universal bound clusters with 3 or more particles depends on the symmetries of the particles. The simplest case in which there are such universal bound clusters is identical spin-0 bosons. Universality also provides strong constraints on reaction rates in the 3-particle and higher sectors [38]. These reaction rates can display dramatic resonant enhancements at low energy. For example, the 3-body recombination rate can increase as 1/E 2 when the center-

JHEP11(2018)084
of-mass collision energy E decreases below the energy scale 1/mr 2 0 , and it can approach a large value proportional to a 4 in the low-energy limit.
Large low-energy cross sections for dark matter particles are motivated by discrepancies between observations of the small-scale structure of the universe and simulations based on collisionless cold dark matter [39][40][41][42][43][44]. Observations of dwarf galaxies are inconsistent with the cusp of dark matter at the center of a galaxy that is predicted by dark-matteronly simulations. Dark-matter-only simulations also imply that dwarf galaxies bound to the Milky Way should be denser than those that have been observed. Although other explanations for these problems have been proposed, they can all be solved by self-interacting dark matter that is strongly interacting at low energy [45,46]. Short-range interactions with a large scattering length provide a particularly predictive model of self-interactions that become strong at low energy [1][2][3].
If dark nucleons have short-range interactions with a large scattering length a, they have universal low-energy properties determined by a [1]. We denote the dark nucleon by d and a bound cluster of n dark nucleons by d n . If a is negative, there are no weakly bound clusters of two dark nucleons. If a is positive, there is one universal weakly bound cluster: the dark deuteron d 2 . If a pair of dark nucleons has annihilation channels, the scattering length a is complex with a small negative imaginary part. In addition to the elastic cross section and the binding energy of the dark deuteron, the annihilation rate of a pair of dark nucleons and the decay rate of the dark deuteron are also universal in the sense that they are determined by the complex scattering length [1]. In a direct detection experiment, the dark deuteron can scatter elastically from a target nucleus, or it can be broken up by the collision [2,3,47]. The low-energy cross sections for both processes are determined by a up to a multiplicative factor. Their dependence on the collision energy and the recoil angle provides interesting signatures for this simplest dark nucleus. The simplest reaction that can form a bound cluster is 3-body recombination: d + d + d → d 2 + d. In an expanding and cooling thermal system, such as the early universe, the decreasing number density will tend to suppress 3-body recombination while the decreasing temperature will tend to enhance it. Once dark deuterons are produced, the competing breakup reaction d 2 +d → d+d+d will destroy them. The net effect on the population of d 2 can only be determined by explicit calculations.
In this paper, we study 3-body recombination into dark deuterons during the Hubble expansion in the early universe under the assumption that the dark matter consists of dark nucleons that are identical bosons with a large positive scattering length, which is the most favorable case for the formation of universal bound clusters. We determine the mass m χ and the scattering length a of the dark nucleon that would be required to solve the small-scale structure problems of the universe. For these values of m χ and a, the fraction of dark matter in the form of dark deuterons can increase by orders of magnitude when the dark-matter temperature decreases to below the binding energy of the dark deuteron. However, we find that a significant population of dark deuterons cannot be produced during the Hubble expansion. Since the production of the dark deuteron is a bottleneck, larger dark nuclei will also not be formed. A much larger population of dark deuterons can be produced if the constraints from small-scale structure are relaxed and the mass of the dark matter particle is decreased.

JHEP11(2018)084
In section 2, we summarize the universal 2-body physics of particles with a large scattering length. We also determine the mass and the scattering length of the dark nucleon that would be required to solve small-scale structure problems of the universe. In section 3, we summarize the universal 3-body physics of identical bosons with large scattering length. In section 4, we present results for the rate constants for many-body systems of identical bosons with large scattering length in thermal equilibrium. In section 5, we consider the formation of dark deuterons by 3-body recombination during the Hubble expansion of the early universe. We calculate the fraction of dark matter in the form of dark deuterons as a function of the red shift. Our results are summarized and discussed in section 6.

Universal two-body physics with large scattering length
In this section, we summarize the universal two-body physics of particles with short-range self-interactions and a large scattering length. We determine the mass and the large scattering length of a dark nucleon that would be required to solve the small-scale structure problems of the universe.

Two-body physics
Atomic physics has provided a strong impetus for developing the universal few-body physics of particles with large scattering lengths [38]. There are naturally occurring atoms with large scattering lengths, such as the 4 He atom. There are other atoms whose scattering lengths can be controlled and made arbitrarily large by using Feshbach resonances [48]. In this subsection, we use the concise language of atomic physics for the particles with large scattering lengths and their bound clusters. The particle d is referred to as an atom, and the two-body bound cluster d 2 is called a dimer. We make factors of Planck's constant explicit.
We denote the mass of the atom d by m. The atom has short-range self-interactions with range r 0 and an S-wave scattering length a that is much larger than r 0 . The range and the scattering length provide a high energy scale E 0 = 2 /mr 2 0 and a low energy scale E 2 = 2 /ma 2 . At energies well below E 0 , the two-body physics is universal in the sense that it is completely determined by a. It depends on the nature of the particles and on the details of their short-range interactions only through a. The universal behavior becomes exact in the zero-range limit r 0 → 0. In this limit, all higher partial-wave interactions go to 0, so two-body scattering is purely S-wave.
The universal region for the scattering of two atoms is when the collision energy E, which is the kinetic energy in the center-of-mass frame, is well below E 0 . The universal elastic scattering cross section for identical bosons is If the two colliding atoms are distinguishable particles, such as the two spin states of a spin-1 2 fermion, the numerator is replaced by 4π. The cross section has dramatic energy dependence. When the collision energy E decreases below E 0 , the elastic cross section JHEP11(2018)084 increases in accordance with eq. (2.1). In the scaling region E 2 E E 0 , the cross section nearly saturates the S-wave unitarity bound 8π 2 /mE. As the energy decreases below E 2 , the cross section levels off and approaches its maximum value 8πa 2 as E → 0. In the limit a → ±∞, the scaling behavior 8π 2 /mE extends down to arbitrarily low energy. Since this cross section saturates the S-wave unitarity bound, the limit a → ±∞ is called the unitary limit.
A universal bound state is one that has properties determined by a. Its binding energy per pair of particles must be less than E 0 . Whether or not there is a universal dimer d 2 depends on the sign of a. If a < 0, there is no universal dimer. If a > 0, there is a single universal dimer. The universal binding energy of d 2 in the zero-range limit is (2.2) A beautiful example in atomic physics of a boson with a large scattering length is the 4 He atom. Its scattering length is about 200 a 0 , where a 0 is the Bohr radius. The scattering length is larger than the effective range by about a factor of 15, so the cross section increases at low energies by more than two orders of magnitude. The 4 He dimer is a universal two-body bound state with the tiny binding energy E 2 = 1.4 × 10 −7 eV. The 4 He dimer was first observed in 1993 using electron impact ionization [49]. The universal low-energy behavior of particles with a large scattering length is illustrated even more dramatically by experiments with ultracold trapped atoms. The scattering length a of the atoms can be controlled and made arbitrarily large by tuning the magnetic field to a Feshbach resonance [48]. Thus the binding energy of the universal dimer can be controlled and made arbitrarily small.
If the atoms have inelastic scattering channels, the scattering length a is complex with a negative imaginary part. If all the inelastic scattering channels have energy release large compared to E 2 , the inclusive inelastic cross section is also universal and determined by a. The universal inelastic scattering cross section for identical bosons is We have assumed the imaginary part of 1/a is tiny compared to the real part of 1/a, in which case the imaginary part can be ignored except in the numerator where it appears as a multiplicative factor. lnelastic atom-atom scattering channels are also decay channels for the dimer. The universal expression for the decay rate is The imaginary part of a should be ignored in the denominator. The energy Γ 2 is twice the imaginary part of the complex binding energy given by eq. (2.2) with complex a. Note that the imaginary part of 1/a cancels in the ratio of the inelastic cross section in eq. Self-interaction reaction rate v σ elastic for dark matter particles as a function of the mean velocity v . The data points are results from Kaplinghat, Tulin, and Yu for dwarf galaxies (red), low-surface-brightness galaxies (blue), and galaxy clusters (green) [53]. The curves are the best fit for a dark-photon model with α = 1/137 [53] (dashed) and the best fit to eq. (2.5) (solid). The diagonal lines are for energy-independent cross sections.

Dark matter parameters
The small-scale structure problems of the universe can be solved by self-interacting dark matter that becomes strongly interacting at low energies [46,[50][51][52]. In ref. [53], Kaplinghat, Tulin, and Yu determined self-interaction reaction rates v σ elastic for dark matter particles from astrophysical data on dwarf galaxies, low-surface-brightness galaxies, and galaxy clusters [54][55][56][57]. Their data points are shown as a function of the mean relative velocity v of the dark atoms in figure 1. In the galaxies, v ranges from about 20 km/s to about 200 km/s. The values of v σ elastic for the galaxies only are roughly compatible with an energy-independent cross section with σ elastic /m = 2 cm 2 /g. In the galaxy clusters, v is about 2000 km/s. The values of v σ elastic for the clusters only are compatible with an energy-independent cross section with σ elastic /m = 0.1 cm 2 /g. To fit the results for both the galaxies and the clusters requires a cross section that increases dramatically with decreasing velocity. The results for v σ elastic versus v can be fit by a dark-photon model with three parameters: the dark matter mass m χ , the dark photon mass µ, and the coupling constant α for a Yukawa potential. Kaplinghat

JHEP11(2018)084
The results for the self-interaction reaction rates in ref. [53] can be fit equally well by a short-range interaction model with a large scattering length. We assume dark nucleons are identical spin-0 bosons with a large real and positive scattering length. The parameters required to describe the universal two-body physics of dark nucleons are their mass m χ and the scattering length a. The elastic cross section is given in eq. (2.1). The reaction rate as a function of the relative velocity v is (2.5) Our fit to the data points for v σ elastic versus v shown in figure 1 gives The curve for the best fit with m χ = 19 GeV and a = ±17 fm is shown in figure 1. The binding energy E 2 = 1/(m χ a 2 ) of the dark deuteron is predicted to be 7.1 keV. This is also the value of the collision energy m χ v 2 /4 where v σ elastic (v) has a maximum as a function of v. The relative velocity v at the maximum is about 300 km/s. The reaction rate in eq. (2.5) must remain accurate for v beyond the values of v for galaxy clusters. At some larger velocity scale v 0 = 2/m χ r 0 set by the range r 0 of self-interactions, the reaction rate in eq. (2.5) may cross over to that for an energy-independent cross section, which is a diagonal line in figure 1. Assuming the crossover does not occur until v is at least 3 times larger than v for galaxy clusters, the energy-independent cross section must satisfy σ elastic /m χ < 0.01 cm 2 /g. The range r 0 must be less than 0.5 fm, and the energy scale E 0 = 1/m χ r 2 0 set by the range must be greater than 200 MeV. One can obtain a very similar curve in figure 1 with spin-1 fermions, for which the factor 8π in eq. (2.5) is replaced by 4π. The best-fit parameters for the mass and scattering length are 15 GeV and ±22 fm. This mass is the same as that obtained in the dark photon model of ref. [53]. The same scattering length could be obtained in that model by tuning either the dark photon mass or the Yukawa coupling constant. The mapping from the parameters of a model of dark fermions with gauge bosons to the scattering length is extensively discussed in refs. [58][59][60]. The mapping from the parameters of a more fundamental dark matter model with bosons to the scattering length could be as nontrivial as the mapping from the parameters of QCD to the large neutron scattering length.
An upper bound on the elastic cross section for dark matter particles has been obtained from the Bullet Cluster, which is the result of a collision of two galaxies with a relative velocity estimated to be O(1000) km/s [61][62][63][64][65][66][67]. The apparent absence of significant scattering from the two dark matter halos implies an upper bound on the elastic cross section for the dark matter particles divided by their mass. If the dark matter particles have an energyindependent cross section, the upper bound on σ elastic /m χ is roughly 1 cm 2 /g [68][69][70]. The curves in figure 1 at v = 1000 km/s are compatible with this bound.
Another constraint on the elastic cross section for dark matter particles can be obtained by demanding that self-scattering removes the cusp in the dark matter distribution JHEP11(2018)084 at the center of dwarf galaxies that is predicted by the ΛCDM model. If the dark matter particles have an energy-independent cross section, this condition provides an estimate of σ elastic /m χ that is roughly 1 cm 2 /g [46]. A typical mean velocity of dark matter particles in a dwarf galaxy is 10 km/s. The curves in figure 1 at v = 10 km/s are compatible with this estimate.
Note that in our fit to eq. (2.5) we have only considered the data compiled by Kaplinghat et al. in ref. [53]. We are aware that there are analyses of astrophysical systems that exhibit agreement with collisionless dark matter (see e.g. refs. [68,71]). There are caveats to these analyses, since the galaxy-dark matter offsets are predicted by strongly interacting dark matter are small [72] and since the concentration parameter of dwarf galaxies may be higher than assumed. 1 Whether or not self-interacting dark matter is required in astrophysical systems requires more research. We will use the fit parameters in eq. (2.6) to illustrate the near-threshold S-wave resonance model, but we will also consider values of the parameters that do not solve the small-scale structure problems.

Universal three-body physics of identical bosons
In this section, we summarize the universal three-body physics of identical bosons with a large scattering length, which is surprisingly intricate [38]. The 3-body physics depends strongly on the scattering length a. It also depends log-periodically on a 3-body parameter κ * that can be determined from the binding energy of a universal bound 3-body cluster.
In this section, we use the concise language of atomic physics for the particles and their bound clusters. The particle d is referred to as an atom. A two-body bound cluster d 2 , a three-body bound cluster d 3 , and a four-body bound cluster d 4 are called a dimer, a trimer, and a tetramer, respectively. We make factors of Planck's constant explicit.

Trimer spectrum
The remarkable nature of trimers composed of identical bosons with a large scattering length was first realized by Vitaly Efimov. In 1970, Efimov pointed out that in the unitary limit where a is infinitely large, there is a sequence of infinitely many trimers whose binding energies have an accumulation point at the 3-atom scattering threshold [73]. The ratio of the binding energies of two successive trimers is the square of a universal number λ 0 = 22.694. The order of magnitude of the binding energy of the most deeply bound Efimov trimer is the energy scale E 0 = 2 /mr 2 0 set by the range. The discrete spectrum of Efimov trimers in the unitary limit a = ±∞ implies that few-body physics in the zero-range limit must depend on a 3-body parameter. The Efimov trimers can be labeled by an integer n. A convenient choice for the 3-body parameter is the binding wave number κ * in the unitary limit a = ±∞ of some arbitrarily chosen Efimov trimer labelled by n * . In the unitary limit, the binding energies of the other Efimov trimers differ by integer powers of λ 2 0 ≈ 515:

JHEP11(2018)084
If the binding wave number of a different Efimov trimer was chosen as the 3-body parameter, the value of κ * would differ by a multiplicative factor that is an integer power of λ 0 . Since κ * can only be defined modulo multiplicative factors of λ 0 , few-body physics can only depend log-periodically on κ * . In particular, 3-body reaction rates must be functions of a and κ * that are invariant under replacing κ * by λ 0 κ * . The binding energies of Efimov trimers are smooth functions of the inverse scattering length 1/a [74]. If the scattering length is not infinitely large, there are only a finite number of Efimov trimers. As 1/a decreases through negative values, Efimov trimers disappear through the 3-atom scattering threshold at critical values of a that differ by multiplicative factors of λ 0 . As 1/a increases through positive values, Efimov trimers disappear through the atom-dimer scattering threshold at critical values of a that differ by multiplicative factors of λ 0 . The Efimov trimer whose binding momentum in the unitary limit is κ * disappears through the 3-atom scattering threshold E = 0 at the negative scattering length a − = −1.508 κ −1 * [75], and it disappears through the atom-dimer scattering threshold E = −E 2 at the positive scattering length a * ≈ 0.07076 κ −1 * [38]. Given any large scattering length a, the energy of one Efimov trimer can be used to determine κ * and the binding energies of the other Efimov trimers can then be predicted. The number of Efimov trimers is not predicted, because the binding energy of the deepest Efimov trimer is determined by the range r 0 .
In atomic physics, the two 4 He trimers are beautiful examples of Efimov trimers. The binding energy of the more deeply bound 4 He trimer is about 1.1 × 10 −5 eV. It was first observed using diffraction from a transmission grating [76]. The binding energy of the more weakly bound 4 He trimer relative to the 3-atom threshold is about 2.3 × 10 −7 eV, which is about a factor of 2 larger than that of the 4 He dimer. It has been observed only recently using Coulomb explosion imaging [77]. The first Efimov trimer observed in cold atom physics was a 133 Cs trimer observed in 2008 as a resonance in the atom loss rate from 3-body recombination [78].

Dimer-atom scattering
The dimer-atom scattering processes are elastic scattering (d 2 + d → d 2 + d) and dimerbreakup scattering (d 2 + d → d + d + d). The collision energy, which is the total kinetic energy of the atom and dimer in the center-of-momentum frame, is where k is the relative momentum of the atom and dimer. The partial wave expansion for the elastic scattering amplitude is The phase shifts δ J (k) are dimensionless functions of k. The scattering is purely elastic for energies between the atom-dimer threshold E = 0 and the dimer-breakup threshold E = E 2 . The phase shifts are therefore real for E < E 2 and complex for E > E 2 . The

JHEP11(2018)084
cross sections for elastic scattering and for breakup scattering can be expressed in terms of the phase shifts: In the universal regime where the energy E is much smaller than the energy scale E 0 set by the range, the only relevant interaction parameters are the scattering length a and the 3-body parameter κ * . The S-wave phase shift δ 0 (k) is a dimensionless function of ka and aκ * that depends only log-periodically on aκ * . It can be expressed in the form [38] where s 0 = π/ log λ 0 ≈ 1.00624 is a universal constant and a + is an alternative 3-body parameter that differs from κ −1 * by a multiplicative factor: a + = 0.3165 κ −1 * [75]. The dimensionless functions s 11 (x), s 12 (x), and s 22 (x) are complex-valued functions of the scaling variable x = ka. In terms of this variable, the dimer-breakup threshold 3, the scaling functions s 11 (x), s 12 (x), and s 22 (x) are entries of a 2 × 2 unitary matrix. For x > 2/ √ 3, they are entries of a 2 × 2 submatrix of a 3 × 3 unitary matrix. In ref. [79], they were calculated numerically over the range of x from 10 −1 to 10 +1 . The phase shifts δ J (k) for the higher partial waves are also dimensionless functions of x = ka. They are real for x < 2/ √ 3 and complex for x > 2/ √ 3. In ref. [79], the phase shifts for J from 1 to 6 were calculated numerically over the range of x from 10 −1 to 10 +1 . The results of ref. [79] allow the elastic cross section and the breakup cross section to be calculated for collision energies up to 100 E 2 , where E 2 is the dimer binding energy in eq. (2.2).
The breakup cross section is shown in figure 2 as bands whose envelope corresponds to minimizing and maximizing the cross section with respect to a + . The upper band is the total cross section, and the lower band is the contribution from J = 0. The curves inside the lower band are for 8 values of a + /a between 1 and λ 0 that are equally spaced on a log scale. At E = 100 E 2 , the sum of the higher partial waves is larger than the maximum J = 0 contribution by more than an order of magnitude. The behavior of the individual partial waves at large E is consistent with decreasing as 1/E 2 , but the sum of the partial waves is consistent with decreasing as 1/E. This is the scaling behavior of the cross section at high energy that is expected from dimensional analysis, given that the only energy scale from interactions is E 2 . Both the upper and lower bands are extrapolated beyond 100 E 2 by fitting to power-law scalings between 50 E 2 and 100 E 2 . The breakup cross section in the high-energy limit can be approximated as where the coefficient is c 2 ≈ 35. where a * = 0.0708 κ −1 * . The cross section at E = 0 depends log-periodically on a/a * , and its value ranges from 0 to ∞. It diverges at a * = a and at other values of a * that differ from a by an integer power of λ 0 = 22.7, because there is an Efimov trimer at the atom-dimer threshold. The cross section at E = 0 vanishes at a * = 2.63 a and at other values of a * that differ from 2.63 a by an integer power of λ 0 , because there is destructive interference between two scattering pathways. The S-wave contribution to the elastic cross section at the dimer-breakup threshold is known analytically. It can be approximated with an accuracy of better than 1% by [38] σ (J=0) where a + ≈ 0.3165 κ −1 * . This contribution vanishes at a + = a and at other values of a + that differ from a by an integer power of λ 0 , because there is perfect destructive interference between two reaction pathways. The S-wave contribution to the breakup cross section at an energy E just above the dimer-breakup threshold E 2 is [79] (3.9)

JHEP11(2018)084
The coefficient C 3 (a/a + ) in the prefactor depends log-periodically on a. A completely analytic expression for this coefficient has been derived by Macek, Ovchinnikov, and Gasaneo [80]. It can be expressed as where s 0 = π/ log λ 0 = 1.00624 and a + ≈ 0.3165 κ −1 * [75]. It vanishes at a + = a and at other values of a + that differ from a by an integer power of λ 0 = 22.7. At these values of a + , there is perfect destructive interference between two recombination pathways. The coefficient in (3.10) can be approximated with an error of less than 1% by the simpler expression C 3 (a/a + ) ≈ 67.1 sin 2 [s 0 log(a/a + )]. (3.11)

Three-body recombination
Three-body recombination is a reaction in which the collision of three atoms results in the formation of a dimer: The reaction rate depends on the wave vectors k 1 , k 2 , and k 3 of the three colliding atoms, but not on the total wave vector k 1 + k 2 + k 3 . It can be expressed as a function of the Jacobi wave vectors defined by k 12 = k 1 − k 2 and k 3,12 = k 3 − 1 2 (k 1 + k 2 ). The collision energy E is the kinetic energy in the center-of-momentum frame: . (3.12) The recombination rate can be expressed as a function of E and 5 dimensionless hyperangles consisting of the spherical angles of the two Jacobi vectors and arctan( √ 3 k 12 /2k 3,12 ). The hyperangular average of the recombination rate is a function of E only. It can be expressed in terms of the breakup cross section at the kinetic energy E 2 + E [79]: In the universal regime where the collision energy is much smaller than the energy scale 2 /mr 2 0 set by the range, the only relevant interaction parameters are the scattering length a and the 3-body parameter κ * . A completely analytic expression for the three-body recombination rate at zero collision energy has been derived by Macek, Ovchinnikov, and Gasaneo [80]. It can be expressed as (3.14) The coefficient C 3 (a/a + ) depends log-periodically on a/a + , and it can be approximated by the expression in eq. (3.11). The recombination rate at E = 0 vanishes at a + = a and at other values of a + that differ from a by an integer power of λ 0 . It has its maximum value 67.1 a 4 /m at a + = 4.76 a and at other values of a + that differ from 4.76 a by an integer power of λ 0 .

Four-body physics and beyond
In 2004, Platter, Hammer, and Meissner predicted the existence of universal 4-body bound clusters composed of identical bosons with large scattering length [81,82]. Their binding energies were mapped out as functions of a by von Stecher, D'Incao, and Greene [83]. In an experiment with 133 Cs atoms in 2009, the dramatic increase of the 4-body recombination rate at low temperature near a specific value of a was used to discover the first such universal tetramer [84]. There is theoretical evidence for universal bound clusters of 5, 6, and even more identical bosons with a large scattering length [85].

Rate coefficients at thermal equilibrium
In this section, we give expressions for the rate coefficients for few-body reactions for identical bosons with large scattering lengths in thermal equilibrium. We use the concise language of atomic physics for the bosons and their bound clusters. We consider a gas of atoms with number density n 1 and dimers with number density n 2 in kinetic equilibrium at temperature T but not necessarily in chemical equilibrium. For simplicity, we assume the gas is sufficiently dilute that the Bose-Einstein momentum distributions of the atoms and dimers can be approximated by Maxwell-Boltzmann distributions.

Inelastic atom-atom scattering
We assume all the inelastic atom-atom scattering channels have energy release large compared to E 2 , and that the energetic particles produced by the reaction have scattering cross sections with an atom that are small compared to the elastic atom-atom cross section. The particles produced by an inelastic reaction can therefore be ignored, and the only effect of the reaction is to decrease the number of atoms by 2. In a homogeneous system, the rate at which the number density n 1 of atoms decreases from inelastic atom-atom scattering is proportional to n 2 1 : d dt n 1 = −2K 1 (T ) n 2 1 . (4.1) The rate coefficient K 1 (T ) depends on the temperature and can be expressed as a weighted integral over the inelastic cross section: Upon inserting the universal approximation to the inelastic cross section for identical bosons in eq. (2.3), we obtain an analytic result: where the dimensionless function g(t) Figure 3 shows the rate coefficient K 1 (T ) and its limiting behaviors: 16π Im[a]/m in the low-T limit and 32π(E 2 /kT ) Im[a]/m in the high-T limit.

Dimer breakup
The dimer-breakup reaction dd 2 → ddd decreases the number of dimers by 1 and increases the number of atoms by 2. We assume the final-state atoms are thermalized by the elastic atom-atom scattering. In a homogeneous system, the rate at which the number density n 2 of dimers decreases from dimer-breakup scattering is proportional to n 1 n 2 : The rate coefficient K 2 (T ) depends on the temperature and can be expressed as a Boltzmann average of the dimer-breakup cross section: The universal approximation to the dimer-breakup cross section is given in eq. (3.4). The universal results for the atom-dimer phase shifts δ J (k) in ref. [79] are obtained up to about 100 E 2 . The breakup cross section is extended above 100 E 2 as shown in figure 2 by fitting to the power-law behavior in eq. (3.6). This allows K 2 (T ) to be calculated for all temperatures kT up to the scale E 0 set by the range. The results for the rate coefficient are shown in figure 4 as bands whose envelope corresponds to minimizing and maximizing the rate with respect to a + . The upper band is the total rate coefficient, and the lower band is the contribution from J = 0. The curves inside the lower band are for 8 values of a + /a between 1 and λ 0 that are equally spaced on a log scale.
We can obtain a simple analytic approximation for K 2 (T ) in the low-temperature limit kT E 2 . In this limit, the breakup cross section in eq.  contribution. The limiting behavior of the rate coefficient is where λ T = (2π 2 /mkT ) 1/2 is the thermal wavelength and C 3 (a/a + ) is the coefficient of a 4 /m in the 3-body recombination rate at zero collision energy in eq. (3.14). This coefficient can be accurately approximated by eq. (3.11). Note that the dimer-breakup rate coefficient in eq. (4.6) is exponentially suppressed by the Boltzmann factor. We can also obtain a simple analytic approximation for K 2 (T ) in the scaling region, where kT is much larger than E 2 and much smaller than the energy scale E 0 = 2 /mr 2 0 set by the range. In the scaling region E 2 kT E 0 , the breakup cross section in eq. (3.9) is dominated by the higher partial-wave contributions. Figure 4 shows that at kT = 100 E 2 , the sum of the higher partial waves is already more than an order of magnitude larger than the maximum J = 0 contribution. The dependence on the S-wave scattering length a can therefore be neglected. Since the interactions provide no other length scales smaller than the range, the dependence of the rate coefficient on T can be determined up to a numerical coefficient by dimensional analysis: where c 2 ≈ 35 is the same coefficient as in eq. (3.6). The extrapolation in T provided by the scaling behavior in eq. (4.7) is shown as a dashed line in figure 4.

Three-body recombination
The three-body recombination reaction ddd → dd 2 increases the number of dimers by 1 and decreases the number of atoms by 2. We assume the final-state atom and the finalstate dimer are thermalized by elastic atom-atom scattering and by elastic atom-dimer scattering, respectively. In a homogeneous system, the rate at which the number density n 2 of dimers increases from 3-body recombination is proportional to n 3 1 : The rate coefficient K 3 (T ) depends on the temperature and can be expressed as a Boltzmann average of the three-body recombination rate: where R is the hyperangular average of the 3-body recombination rate, which is a function of the collision energy E only. The factor of 1/3! compensates for the overcounting of 3body states of the 3 identical bosons in the Boltzmann average. The rate coefficient can be expressed as a weighted integral over the dimer-breakup cross section: (4.10) We can use eq. (4.5) to express K 3 (T ) in terms of K 2 (T ): where λ T = (2π 2 /mkT ) 1/2 is the thermal wavelength. The universal approximation to the dimer-breakup cross section is given in eq. (3.4). The universal results for the atomdimer phase shifts δ J (k) in ref. [79] are obtained up to about 100 E 2 and the breakup cross section is extended above 100 E 2 in figure 2 by fitting to the power-law behavior. This allows the recombination rate coefficient K 3 (T ) to be calculated for all temperatures kT up to scale E 0 set by the range. The results are shown in figure 5 as bands whose envelopes corresponding to minimizing and maximizing the rates with respect to a + . The upper band is the total rate coefficient, and the lower band is the contribution from J = 0. The curves inside the lower band are for 8 values of a + /a between 1 and λ 0 that are equally spaced on a log scale.
We can obtain simple analytic approximations for the 3-body recombination rate coefficient by using the relation between K 2 (T ) and K 3 (T ) in eq. (4.11) and the analytic approximations for K 2 (T ) in eqs. (4.6) and (4.7). In the low-temperature limit kT E 2 , the rate coefficient approaches a constant that depends log-periodically on a/a + : where C 3 (a/a + ) is the coefficient in eq. (3.10), which can be accurately approximated by eq. (3.11). In the scaling region E 2 kT E 0 , the rate coefficient scales as the power of temperature required by dimensional analysis: where c 2 ≈ 35 is the same coefficient as in eq. (3.6). The extrapolation in T provided by the scaling behavior in eq. (4.13) is shown as a dashed line in figure 5.
In experiments with ultracold trapped atoms, the atoms form an extremely dilute gas in the sense that the typical interatom spacing is much larger than the range of the interactions between atoms: n 1 1/3 r 0 1, where n 1 is the density-weighted average of the number density. Three-body recombination can be important in these experiments, because the dimer and atom in the final state often have enough kinetic energy to escape from the trapping potential. In that case, every recombination event results in the loss of three atoms. In an experiment with 133 Cs atoms in 2005, the dramatic increase of the 3-body recombination rate at low temperature when the scattering length was tuned to near the negative value a − = −1.5 κ −1 * was used to discover an Efimov trimer [78].

Early universe
In this section, we study the production of dark deuterons through three-body recombination of dark nucleons during the Hubble expansion of the early universe under the assumption that the dark nucleons are identical bosons with a large positive scattering

JHEP11(2018)084
length. We calculate the fraction of dark matter in the form of dark deuterons as a function of the redshift.

Rate equations
After the decoupling of dark matter from ordinary matter, the densities of dark nucleons and larger dark nuclei evolve in thermal equilibrium until they are captured by the gravitational potential wells of galaxies. The time evolution is due to the Hubble expansion and to reactions among the dark nuclei. Assuming that the larger dark nuclei are weakly bound, the density and temperature at decoupling are large enough that any larger dark nucleus that is formed is immediately broken up by a collision with a dark nucleon. Thus we can take as an initial condition that the dark matter consists entirely of dark nucleons at the decoupling time.
Given an initial state consisting only of dark nucleons, larger dark nuclei can be formed by N -body recombination reactions in which N dark nuclei collide and some of them form bound states. At sufficiently low dark nucleon number density n 1 , the N -body recombination rate is proportional to n N 1 . Thus if a dark deuteron d 2 exists, the most favorable reaction is 3-body recombination (d+d+d ↔ d 2 +d). Once dark deuterons have been produced, larger dark nuclei can be formed by rearrangement collisions, such as d n + d 2 → d n+1 + d. The formation of dark deuterons is a bottleneck that must be overcome by 3-body recombination in order to form the larger dark nuclei. We wish to determine whether this bottleneck can be overcome in the early universe when the dark matter is still in thermal equilibrium. To answer this, we can ignore dark nuclei d n with n ≥ 3 and consider only the time evolution for dark nucleons and dark deuterons. The only reactions we need to take into account are 3-body recombination and the dark deuteron breakup reaction (d 2 + d ↔ d + d + d). We wish to determine whether a significant population of dark deuterons can be generated in the early universe.
We denote the number densities of the dark nucleon and the dark deuteron by n 1 (t) and n 2 (t). We assume the number densities of dark nuclei with larger dark baryon number are negligible, so the total dark baryon number density is n dark (t) = n 1 (t) + 2n 2 (t). (5.1) The time evolution equations for n 1 (t) and n 2 (t) obtained from the Boltzmann equation are where H is the Hubble function, K 3 (T ), K 2 (T ), and K 1 (T ) are temperature-dependent event rate coefficients, and Γ 2 is the dark deuteron decay rate. The Hubble function H(t) depends on time, being determined by the scale factor a(t) of the universe: H = d ln(a)/dt. The rate coefficients are functions of the temperature T (t) of the dark matter, which also depends on time.

JHEP11(2018)084
We neglect the effects of the annihilation of dark nucleons into ordinary matter. We therefore set K 1 (T ) = 0 and Γ 2 = 0 in the rate equations in eqs. (5.2). If there were such an annihilation process, it would decrease n 1 through annihilation collisions of two dark nucleons and it would decrease n 2 through the annihilation of the two constituents of the dark deuteron. The rates for both processes are determined by the same parameter Im[1/a], which appears as a multiplicative parameter in both K 1 (T ) and Γ 2 . When Γ 2 is much larger than 3H(t), the number density n 2 of dark deuterons decreases exponentially. Any dimers that have been produced by 3-body recombination would decay quickly on a cosmological time scale. The net effect is that n 2 would remain essentially 0, and the decrease in n 1 would be given by the 3-body recombination term in eq. (5.2a) only. Since we ignore the annihilation of dark nucleons, our results for the number density of dark deuterons can be interpreted as upper bounds.
If we ignore the annihilation terms in the evolution equations in eqs. (5.2), we get a simpler equation for the total dark baryon number density: .

(5.4)
The time evolution of the total dark baryon number density does not depend on the dark matter interactions; it is just diluted by the Hubble expansion. It is convenient to use the redshift z as an alternative time variable. The redshift is related to the scale factor a by 1 + z(t) = a(0)/a(t). The solution for the total dark baryon number density in eq. (5.4) can be expressed as where ρ cdm = 2.23 × 10 −30 g/cm 3 is the present average mass density of dark matter in the universe [86] and m χ is the mass of the dark nucleon. The Hubble function in terms of redshift is given by where the Hubble constant is H 0 = 67.8 km s −1 Mpc −1 and the fractions of the critical density of the Universe for CMB photons (Ω γ ), matter (Ω m ), and dark energy (Ω Λ ) are 5.38 × 10 −5 , 0.308, and 0.692, respectively [86].
Since the dark nucleons are nonrelativistic after the decoupling, their temperature T (z) is proportional to the square of their average momentum [87]. On the other hand, the temperature T γ (z) of the photons is proportional to their average momentum. The Hubble expansion changes the momentum of a particle by a factor of 1 + z. Thus the two temperatures are different functions of the redshift:

JHEP11(2018)084
where T (0) is the present temperature of dark matter that has not been captured by gravitational potential wells and T cmb = 2.73 K is the present temperature of the photons. At decoupling, the dark matter and ordinary matter are in thermal equilibrium: T (z dc ) = T γ (z dc ), where z dc is the redshift at decoupling. We are not displaying the dependence on the Standard Model degrees of freedom in these and the following expressions for the temperature for simplicity. The variation due to the three-body parameter is larger than the effects from the decreasing number of relativistic degrees of freedom. The dark matter temperature is therefore The decoupling redshift can be expressed as z dc ≈ T (z dc )/T cmb . If the thermal decoupling of dark matter and ordinary matter occurs not long after their chemical decoupling, the decoupling temperature is given approximately by the dark-matter mass multiplied by a constant: kT (z dc ) ≈ m χ /20 [88]. The resulting estimate for the decoupling redshift is The mass fraction of the dark matter in the form of dark deuterons is If we ignore the annihilation terms in the evolution equations in eqs. (5.2), the dark deuteron fraction satisfies the differential equation We have used dt = −[(1 + z)H] −1 dz to replace the time derivative by a redshift derivative. Given H(z), n dark (z), and T (z) in eqs. (5.6), (5.5), and (5.8), our problem reduces to solving this single differential equation for f 2 (z) subject to the initial condition f 2 (z dc ) = 0, where z dc is given in eq. (5.9).

Approximation in scaling and threshold regions
The evolution equation for the dark deuteron fraction with redshift in eq. (5.11) involves the rate coefficients K 2 (T ) and K 3 (T ). If dark nucleons are identical bosons with a large positive scattering length and if kT is much smaller than the energy scale E 0 set by the range, the rate coefficients are given in eqs. (4.5) and (4.10). The rate coefficients have simple behavior in the low-temperature limit kT E 2 , where E 2 = 1/(m χ a 2 ) is the dark deuteron binding energy, and in the scaling region E 2 kT E 0 , where E 0 = 1/(m χ r 2 0 ) is the energy scale set by the range. We can use those results to determine the qualitative behavior of the dark deuteron fraction in those regions.
In the scaling region E 2 kT E 0 , the rate coefficients K 2 (T ) and K 3 (T ) have the limiting behaviors given in eqs. (4.7) and (4.13). They scale as λ T and λ 4 T , respectively, where λ T is the thermal wavelength, which is proportional to (1 + z) −1 :

JHEP11(2018)084
Since n dark is proportional to (1 + z) 3 , the products K 2 n dark and K 3 n 2 dark are both proportional to (1 + z) 2 . Thus there can be an equilibrium value of f 2 for which the two terms on the right side of eq. (5.11) cancel. The ratio of K 3 (T ) and K 2 (T ) is given in eq. (4.11). The equilibrium fraction satisfies If the equilibrium value of f 2 is much less than 1, it can be approximated by the right side of eq. (5.13). Upon inserting the estimate for the decoupling redshift in eq. (5.9), the right side reduces to (ρ cdm /m χ )/(kT cmb ) 3 multiplied by a numerical constant. Since ρ cdm /(kT cmb ) 3 = 0.74 eV, the equilibrium value of f 2 is tiny as long as m χ is orders of magnitude larger than 1 eV. In the low-temperature region kT E 2 , the rate coefficients K 2 (T ) and K 3 (T ) have the limiting behaviors given in eqs. (4.6) and (4.12). They are proportional to λ −6 T e −λ 2 T /a 2 and λ −6 T , respectively. The dark deuteron breakup is exponentially suppressed by the Boltzmann factor, so the breakup term in the rate equation can be dropped. If the value of f 2 is much less than 1, we need to keep only the leading terms in f 2 in the recombination term in eq. (5.11). The rate equation then simplifies to When z 10 4 , the Hubble function in eq. (5.6) can be approximated as H(z) ≈ H 0 Ω 1/2 γ z 2 . The solution of eq. (5.14) is then where f 2 (0) is the present dark-deuteron fraction for dark matter that has not been captured by gravitational potential wells and C 3 (a/a + ) can be approximated by eq. (3.11).
The approach to f 2 (0) is predicted to be z 4 multiplied by a coefficient whose dependence on a is C 3 (a/a + ) a 4 . The value of f 2 (0) should be determined by a boundary condition from the region of larger z where kT (z) is comparable to E 2 . We are unable to determine f 2 (0) analytically, but it should depend log-periodically on the three-body parameter a + .
Our numerical results for f (0) are consistent with an expression linear in C 3 . The evolution equation for the dark deuteron fraction f 2 (z) in eq. (5.11) with the rate coefficients K 2 (T ) and K 3 (T ) in eqs. (4.5) and (4.10) applies all the way back to the decoupling redshift provided the decoupling temperature is smaller than the scale E 0 set by the range: kT (z dc ) < E 0 . If E 0 is smaller than kT (z dc ), the simple expressions for K 2 (T ) and K 3 (T ) in eqs. (4.5) and (4.10) are not applicable until kT decreases to below E 0 . However once kT enters the scaling region E 2 kT E 0 , f 2 will be driven quickly to the equilibrium value given by eq. (5.13). Thus the present value of f 2 is completely determined by the scattering length a and the 3-body parameter a + provided only that the decoupling temperature is much larger than the scale E 2 = 1/(m χ a 2 ). This condition kT (z dc ) E 2 can be expressed approximately as m χ a √ 20.

Numerical results
Assuming the dark nucleons are identical bosons with a large scattering length, the fewbody parameters are the dark nucleon mass m χ , the scattering length a, and the three-body parameter a + . For the mass and the scattering length, we use values that can solve smallscale structure problems of the universe. The values that give the best fit to the data points for v σ elastic versus v in figure 1 are m χ = 19 GeV and a = 17 fm. Since the 3-body parameter a + is only defined modulo multiplicative factors of λ 0 ≈ 22.69, the complete range of possibilities is covered by varying a + from a to 22.69 a.
To determine the dark-deuteron mass fraction f 2 (z) as a function of the redshift z, we solve the differential equation in eq. (5.11) subject to the initial condition f 2 (z dc ) = 0. Given the mass m χ = 19 GeV, the decoupling redshift in eq. (5.9) is z dc ≈ 4 × 10 12 . We want to determine whether a significant fraction f 2 is ever generated during the subsequent time evolution.
In figure 6, we show the dark deuteron fraction f 2 as a function of a red-shift variable z dc /z on a log scale. This variable increases from 1 at the decoupling time to infinity at the present time. The band in figure 6 corresponds to minimizing and maximizing f 2 with respect to a + . The individual curves are for eight values of a + that are equally spaced on a log scale between a and 22.69 a. As z decreases from z dc , the dark deuteron fraction f 2 (z) increases very quickly to a plateau of about 4 × 10 −11 from thermal equilibrium between 3body recombination and dark deuteron breakup. That equilibrium value is consistent with the estimate for the scaling region in eq. (5.13). We could therefore just as well take the initial value of f 2 at the decoupling red shift to be the equilibrium value given by eq. (5.13). When the dark matter temperature T decreases below the dark deuteron binding energy 1/(m χ a 2 ) = 7.1 keV, there is a dramatic increase in f 2 by 3 or 4 orders of magnitude. This feature is expected from the exponential suppression of the breakup process in eq. (4.6) and from the z 4 dependence at late times that is predicted by eq. (5.15). The dark-deuteron fraction plateaus at a value f 2 (0) that depends log-periodically on the 3-body parameter a + and can be approximated by f 2 (0) = (4.6 × 10 −8 ) + (6.7 × 10 −9 ) C 3 (a/a + ). (5.16) The maximum value of f 2 (0) from varying a + is larger than the minimum value by a factor of 11. The fraction f 2 (0) has its minimum when the value of a + /a is just a few percent lower than 1 (or equivalently λ 0 = 22.69), which is the value for which there is total destructive interference in the 3-body recombination rate at zero temperature. It has its maximum when the value of a + /a is just a few percent lower than λ The results for f 2 (z) shown in figure 6 are for parameters m χ = 19 GeV and a = 17 fm that solve small-scale structure problems of the universe, as illustrated in figure 1. However there are also mechanisms involving baryonic physics that can solve or at least ameliorate the small-scale structure problems. We first relax the constraint on m χ and a by ignoring the results for v σ elastic versus v from clusters of galaxies. The results in figure 1 from galaxies only are roughly compatible with a cross section that at low velocities approaches σ elastic /m = 2 cm 2 /g. This requires the constraint 8πa 2 /m χ = 2 cm 2 /g. Given m χ , the scattering length a is determined. The results for the dark deuteron fraction f 2 (0) at late times as a function of m χ are shown in figure 7. If m χ is too small, the decoupling temper-JHEP11(2018)084 ature kT (z dc ) ≈ m χ /20 cannot be much larger than the binding energy E 2 = 1/(m χ a 2 ). In this case, the temperatures after decoupling do not include a scaling region in which T E 2 , so f 2 is not determined by a and a + only, but has additional sensitivity to the range r 0 . This additional sensitivity to r 0 is avoided if m χ a √ 20, which implies m χ 0.4 GeV. This requires m χ to be well above the shaded region in figure 7. In the unshaded region, the fraction f 2 (0) scales roughly as m −2.5 χ . At m χ = 1 GeV, the range of f 2 (0) from varying a + is from 6 × 10 −5 to 2 × 10 −3 . Thus a dark deuteron fraction as large as 10 −3 is possible if the dark-matter elastic cross section at low velocities is 2 cm 2 /g.
If the small-scale structure problems of galaxies are ameliorated by mechanisms involving baryonic physics, the constraint on m χ and a becomes the inequality 8πa 2 /m χ < 2 cm 2 /g. At a given value of m χ , this allows the scattering length to be decreased. This can only decrease the dark deuteron fraction f 2 (0) at late times.
The results presented above assume the kinetic decoupling temperature T kdc is very close to the chemical decoupling temperature T dc [89]. We now consider the case when T kdc is significantly smaller than T dc . The value of T kdc depends on the interactions between dark matter and ordinary matter. We treat it as an unknown parameter and simply describe how f 2 (0) scales with T kdc /T dc . During thermal equilibrium, the dark matter temperature is the same as the photon temperature in eq. (5.7b), while after kinetic decoupling the temperature is quadratic in z as in eq. (5.7a). The dark matter temperature can be written as where z kdc ≈ T kdc /T cmb − 1 is the redshift at the kinetic decoupling temperature. By inserting eq. (5.17) into eq. (5.11) instead of eq. (5.8), we obtain the fraction with kinetic decoupling taken into account. Figure 8 shows the fraction f 2 (z = 0) at late times as a function of T kdc /T dc for the masses 1 GeV (upper band) and 19 GeV (lower band). As shown in figure 8, as T kdc decreases, the fraction decreases, scaling roughly as (T kdc /T dc ) 1.9 .

Discussion and conclusion
The predictions of ΛCDM cosmology face a number of challenges at small scales. These small-scale structure problems can be resolved either by effects of baryons on structure formation or by novel dark matter dynamics. Self-interacting dark matter is a paradigm that can solve the small-scale structure problems in ΛCDM cosmology while remaining consistent with other cosmological data. Perhaps the most predictive model of self-interacting dark matter involves a near-threshold S-wave resonance that produces a large scattering length [1]. This near-threshold resonance is a bound state if the scattering length is positive. The signatures for this two-body bound state (darkonium or dark deuteron) in direct detection and directional detection experiments have been studied [2,3]. In this paper, we have studied the production rate of the dark deuteron in the early universe.
We first compared the predictions of the near-threshold S-wave resonance model with cross sections for self-interacting dark matter in different astrophysical objects determined by Kaplinghat, Tulin and Yu [53]. They showed that a dark-photon model with three adjustable parameters can reproduce the velocity dependence of the self-interaction cross section, which spans two orders of magnitude in velocity [53]. We find that the nearthreshold S-wave resonance model provides an equally good fit to these astrophysical data (see figure 1) with only two adjustable parameters: the mass m χ of the dark-matter particle and the scattering length a. The best-fit values are m χ = 19 GeV and a = ±17 fm.
We have assumed the dark nucleons are identical bosons with a large positive scattering length. The smallest universal bound cluster is the dark deuteron d 2 . The simplest reaction that can form this bound cluster is 3-body recombination into the dark deuteron: d+d+d → d 2 + d. The three-body recombination rate is a function of the mass m χ , the scattering length a, and a three-body parameter a + , with the dependence on a + being log-periodic with discrete scaling factor 22.7. If the temperature at decoupling is much larger than the binding energy of the dark deuteron, the present fraction f 2 (0) of dark matter in the form of dark deuterons is completely determined by these three parameters. For m χ = 19 GeV and a = 17 fm, the fraction f 2 (z) at early red shifts has an equilibrium value of about 4 × 10 −11 . When the dark-matter temperature decreases to below the binding energy of the dark deuteron, which occurs at a red shift z ≈ 10 10 , f 2 (z) increases by orders of magnitude to a value between 4 × 10 −8 and 5 × 10 −7 that depends on a + .
The present fraction f 2 (0) of dark matter in the form of dark deuterons can be increased by relaxing the constraint on m χ and a from solving small-scale structure problems and decreasing m χ . However the decoupling temperature must be much larger than the binding energy of the dark deuteron for f 2 (0) to be insensitive to the range of the interactions.

JHEP11(2018)084
Given this constraint, f 2 (0) cannot be larger than about 10 −3 . If the system remains in thermal equilibrium longer after chemical decoupling, the fraction f 2 (0) decreases, scaling approximately as the 2nd power of the ratio of the temperatures for kinetic and chemical decoupling. We conclude that a significant population of dark deuterons cannot be produced in the early universe by 3-body recombination of dark matter particles with a large scattering length. Since the production of dark deuterons is a bottleneck for the formation of larger bound clusters, we conclude that the formation of bound clusters in the early universe would require additional microphysics. An example is a light mediator that allows radiative fusion reactions.
If the large scattering length a is negative, the smallest universal bound clusters are Efimov clusters d 3 ("dark tritons"). The simplest reaction that can form bound clusters is 4-body recombination into a dark triton: d + d + d + d → d 3 + d. The rate for 4-body recombination is suppressed compared to the rate for 3-body recombination by an additional factor of the number density of dark matter particles. Since a significant population of dark deuterons cannot be produced in the early universe by 3-body recombination, a significant population of dark tritons cannot be produced by 4-body recombination either. Since the production of dark tritons is a bottleneck for the production of larger dark nuclei, a significant number of dark nuclei will not be formed in the early universe if the dark nucleons are identical bosons with a large negative scattering length.
Identical bosons are not the only types of particles for which there is dramatic enhancement of the 3-body recombination rate at low temperature when the scattering length is large. The degree to which 3-body recombination is enhanced depends on the symmetries and mass ratios of the particles with large scattering lengths. Three-body recombination requires the three particles to come within a distance of order the de Broglie wavelength of the final-state particles, which is of order 1/a if the collision energy is small. For identical bosons, the 3-body recombination rate K 3 (T ) in the low-temperature limit is proportional to a 4 . If the dark matter consists of the two spin states of a spin-1 2 fermion, K 3 (T → 0) is suppressed by (r 0 /a) 2 , where r 0 is the range, because the Pauli exclusion principle suppresses the contribution from the region where the separations of the three fermions are all of order a. If the dark matter consists of the four spin states of two spin-1 2 fermions, there is no such suppression and K 3 (T → 0) is proportional to a 4 .
We can also show that a significant fraction of dark deuterons cannot form once the dark matter particles fall inside the gravitational potential well of a galaxy. It is easy to put an upper bound on the rate of increase in the dark deuteron fraction in the Milky Way from 3-body recombination. The maximum possible rate of increase in n 2 is (67.1 a 4 /m χ )n 3 1 . The dark matter mass density in the solar system, which is about 8 kpc from the center of the Milky Way, is m χ n 1 = 0.3 GeV/cm 3 . If feedback between strongly interacting dark matter and baryons is taken into account, the radius of the dark matter core of the Milky Way may be about 0.3 kpc [90]. The dark matter mass density in the core of the Milky Way may be about m χ n 1 = 8 GeV/cm 3 . For m χ = 19 GeV and a = 17 fm, the maximum rate of increase of f 2 is about 10 −51 /s. The age of the Milky Way is about 10 Gyr ≈ 3 × 10 17 s, so we see that the dark deuteron fraction remains negligible. If we relax the constraints on m χ and a from solving small-scale structure problems but keep the binding energy of the JHEP11(2018)084 dark deuteron small compared to the decoupling temperature, the rate of increase of f 2 can be made larger at most by about an order of magnitude. Dwarf galaxies can have higher dark matter densities than the Milky Way, but the rate of increase in the dark deuteron fraction is small in those systems too.
Although 3-body recombination of dark matter particles is unable to build up a large fraction of dark deuterons in the early universe, it may still have a significant effect on dark matter annihilation. If a pair of dark matter particles has an annihilation scattering channel, the constituents of a dark deuteron will eventually annihilate once the dark deuteron is formed. Three-body recombination therefore provides an additional annihilation channel. If the dark matter particles have a large scattering length a, the annihilation scattering cross section and the dark deuteron decay rate are both determined by m χ and a up to a multiplicative constant that cancels in their ratio. The resonant enhancement of annihilation scattering can induce a second period of dark matter annihilation after the thermal freezeout [91]. The effects of reannihilation have been explored in a dark photon model [91]. A near-threshold S-wave resonance model provides a more predictive framework in which the effects of reannihilation through 3-body recombination can also be easily explored.