Bound-state formation for thermal relic dark matter and unitarity

We show that the relic abundance of thermal dark matter annihilating via a long-range interaction, is significantly affected by the formation and decay of dark matter bound states in the early universe, if the dark matter mass is above a few TeV. We determine the coupling required to obtain the observed dark matter density, taking into account both the direct 2-to-2 annihilations and the formation of bound states, and provide an analytical fit. We argue that the unitarity limit on the inelastic cross-section is realized only if dark matter annihilates via a long-range interaction, and we determine the upper bound on the mass of thermal-relic dark matter to be about 197 (139) TeV for (non)-self-conjugate dark matter.


I. INTRODUCTION
The annihilation cross-section of dark matter (DM) is largely what delimits the possibilities for its production and determines the expectations for its phenomenology. In the standard DM paradigm, DM is assumed to have been in chemical equilibrium with the thermal bath in the early universe, due to rapid annihilation and paircreation processes. As the universe expanded and cooled, these processes became inefficient, and the comoving DM density froze-out. In this scenario, the observed DM abundance determines the DM annihilation cross-section, provided that the DM particles are massive enough to have become non-relativistic at freeze-out. Weakly interacting massive particles naturally possess annihilation cross-sections in the vicinity of the required value [1][2][3][4][5][6][7][8][9], and constitute the standard candidate for thermal-relic DM; however, thermal-relic DM may also reside in a hidden sector [10].
Here, we focus on heavy DM, with mass m TeV. We discuss possibilities for the physics underlying the DM annihilation in this mass regime, and the couplings required to obtain the observed DM abundance from thermal freeze-out. Heavy DM is of particular interest in view of the upcoming 14 TeV run of the LHC, as well as future high-energy experiments, such as a 100 TeV collider [11]. It is already being probed by direct and indirect detection experiments [12][13][14][15][16]; in fact, DM with mass m 500 GeV has been invoked to explain the highenergy positron excess observed by PAMELA, Fermi and AMS [15][16][17]. The precise knowledge of the DM couplings required for efficient annihilation in the early universe is essential in interpreting the experimental data.
The DM annihilation processes may be either due to short-range interactions mediated by heavy particles, as for example the weak interactions of the Standard Model (SM), or due to long-range interactions mediated by light species. The latter possibility becomes of interest typically when DM is hypothesized to reside in a hidden sector [10]; however, even the weak interactions of the SM can manifest as long-range if DM is heavier than a few TeV [18,19]. For s-wave annihilation via a short-range interaction, the annihilation crosssection times relative velocity required to obtain the observed DM density from the freeze-out of thermal particles is σ ann v rel c 4.4 × 10 −26 cm 3 s −1 , assuming nonself-conjugate DM with mass above 10 GeV [20]. If DM couples to a light force mediator, the long-range interaction between two incoming DM particles distorts their wavepackets; this is the well-known Sommerfeld effect [21]. As a result, σ ann v rel is enhanced at low velocities. This Sommerfeld enhancement (SE) depends on the coupling strength of the DM to the light force mediator (c.f. Eq. (2)). The efficient annihilation of heavy DM in the early universe requires a large coupling, which renders SE significant during freeze-out. Indeed, the SE of the 2-to-2 annihilation processes affects the abundance of thermal-relic DM if the DM mass is m 800 GeV [22][23][24][25]; as a result, the coupling required to reproduce the observed DM density is lower than that estimated from σ ann v rel c in the absence of SE.
Attractive long-range interactions imply also the existence of bound states. Particle-antiparticle bound states, as well as bound states of self-conjugate identical particles, decay promptly into the force carrier particles; their formation thus contributes to the annihilation rate of their constituent species. In the early universe, the formation of DM bound states -a process which is also enhanced at low velocities by the Sommerfeld effectcan catalyze the DM annihilation, and affect the relation between the DM relic abundance and the DM couplings.
In this paper, we investigate this effect. For concreteness, we consider DM consisting of Dirac fermions X of mass m, charged under an unbroken dark U(1) gauge force. We show that the formation of dark positroniumlike states in the early universe and their subsequent decay into dark photons γ, affect the DM relic abundance for DM masses m few TeV. We calculate the dark finestructure constant α which yields the observed DM density, taking into account the formation of bound states, and their ionization and decay in the early universe.
The importance of considering long-range interactions is underscored by unitarity. It has long been shown that unitarity and the thermodynamics of the early universe set an upper bound on the mass of thermal-relic DM [26]. Indeed, unitarity implies an upper limit on the inelastic DM self-interaction cross-section which decreases with increasing DM mass; for s-wave annihilation, in the nonrelativistic regime, this is [26] where m is the DM mass and v rel is the relative velocity of the incoming DM particles. The inelastic cross-section includes, of course, all the processes which may result in the annihilation of DM particles. Notably, the dependence of (σ inel v rel ) max on v rel implies that the maximum value of the inelastic cross-section can be realized if the DM particles interact via a light or massless force carrier, giving rise to the Sommerfeld enhancement at low velocities exhibited by Eq. (1). In the following, we determine the unitarity bound on the mass of thermal DM, by employing Boltzmann equations to account properly for the late-time annihilations occurring as a result of the SE. Close to the unitarity limit, bound-state formation (BSF) is the most efficient annihilation channel. The effect of BSF on the DM relic abundance has been previously considered in Ref. [27], which found it to be negligible for DM masses m 10 TeV. While our results are reasonably consistent with this claim, we find some discrepancies in the way the cross-section for radiative formation of positronium-like states was adapted in [27] from [28]; this led to underestimating the effect. Here we show that BSF affects the relic abundance of DM with mass above a few TeV. We determine the unitarity bound on the mass of thermal-relic DM, and discuss BSF in relation to unitarity.

II. CROSS-SECTIONS AND RATES
The direct annihilation of DM and the formation of DM bound states both contribute to the inelastic scattering of DM. Once they form, bound states may either decay into dark photons, or get ionized by the ambient radiation. The DM relic abundance depends on the balance of these processes. Below we list the pertinent crosssections and rates.
Bound states can form via emission of a dark photon. The cross-section times relative velocity can be conveniently cast in the form, At ζ 1, S BSF (ζ) 2 10 πζ/(3 exp 4). Clearly, at large ζ, BSF becomes more rapid than the direct DM annihilation into dark photons, with S BSF (ζ)/S ann (ζ) 3.1. Moreover, capture to higher-level states yields an overall enhancement factor for BSF of at most π 2 /6 1.6 [27]. However, because the capture into the ground state dominates, and the excited states are longer-lived, in the following we consider only capture to the ground state.
Note that in the non-relativistic regime, neglecting the spin-orbit coupling, the BSF cross-section is independent of the spin configuration of the incoming particles, which remains conserved in the process. σ BSF is the crosssection for any such process. The spin-averaged crosssections for the formation of para-and ortho-states are σ BSF,↑↓ = σ BSF /4 and σ BSF,↑↑ = 3σ BSF /4 respectively.

C. Thermal average
To estimate the effect of annihilations and BSF on the DM abundance, we need to average the respective rates over the momentum distribution of DM in the early universe. It will be convenient to define the time variables where T is the temperature of the dark plasma, and T X is the temperature of the DM particles. We discuss their relation in the next section. Assuming a Maxwellian velocity distribution for the DM particles, the thermally averaged cross-sections for annihilation is σ ann v rel = σ 0Sann (z X ), where [24] S ann (z X ) = For BSF, we shall include, for completeness, the Bose enhancement due to the final-state dark photon, which is emitted with energy ω ∆+µv 2 rel /2. The Bose enhancement remains important for T ω, i.e. typically until after freeze-out, though during this time the ionization of bound states is still rapid and impedes efficient DM annihilation via BSF (see below). The BSF rate is pro- At z, z X 1,S BSF 2 11 √ πz X /(3 exp 4) andS ann 4 √ πz X .S ann andS BSF are shown in Fig. 1.
The ionization cross-section of the bound states, σ ion , is related to σ BSF by the Milne relation, σ ion /σ BSF = µ 2 v 2 rel /(2ω 2 ) [31], where the factor of 2 counts the photon polarizations and ω ∆ + µv 2 rel /2 is the photon energy. The thermally averaged ionization rate is where III. TIMELINE DM remains in chemical equilibrium with the dark photons due to annihilation, BSF and the inverse processes. As usual, we define the freeze-out of these processes as the time when the DM abundance differs from the equilibrium value by a factor of order 1. We estimate x = x f at freeze-out by adapting the standard result [6] to incorporate the SE of σ ann and σ BSF , where z f = α 2 x f /4, g X = 2 are the degrees of freedom of the DM species (the antiparticles are counted separately), and g * is the number of relativistic degrees of freedom.
Typically, x f ∼ 25. Because of the SE, DM annihilations (direct and via BSF) continue to be significant after freeze-out. To account for this, we integrate the Boltzmann equations for z > z f , as described below.
After chemical decoupling, the DM particles remain in kinetic equilibrium with the dark photons via Thomson scattering, typically until rather late. The kinetic decoupling of the DM from the dark radiation occurs at z = z kd ∼ 10 2 (α/0.02) 3 (TeV/m) 1/2 [24]. At z z kd , z X = z. After kinetic decoupling, provided that the effect of any residual interactions is negligible, z X = z 2 /z kd [6].
Here, we assume that the dark photons are at the same temperature as the plasma of SM particles. This is a viable possibility during the freeze-out of DM with mass m 100 GeV. Indeed, DM then freezes-out at a temperature T f = m/x f ∼ m/25 4 GeV, before the QCD phase transition. Provided that the dark radiation has decoupled from the SM at that time, the subsequent decoupling of the QCD degrees of freedom reheats the SM plasma, and leads to the dark radiation having a lower temperature than ordinary photons; this ensures that the BBN and CMB constraints on the total relativistic energy of the universe are satisfied. It is, of course, straightforward to generalize our calculation to the case that the dark radiation bath and the SM plasma are at different temperatures during DM freeze-out.
Using the cross-sections and rates of the previous section, in table I we outline the order in which various effects become significant. Importantly, the efficiency of BSF in annihilating the thermal population of DM depends not only on σ BSF , but also on the balance between the decay and ionization of the bound states which are formed.

IV. BOLTZMANN EQUATIONS
Let Y X ≡ n X /s, Y ↑↓ ≡ n ↑↓ /s and Y ↑↑ ≡ n ↑↑ /s, where n X , n ↑↓ and n ↑↑ are the number densities of the unbound X particles, the para-and the ortho-bound states respectively. s = (2π 2 /45)g * S T 3 is the entropy density of the universe, with g * S being the entropic relativistic degrees of freedom. The abundances of the unbound and bound DM particles are governed by the coupled Boltzmann equations where Y eq = 0.14 (g X /g * S ) x 3/2 e −x is the equilibrium number density of the X particles normalized to s, c 1 ≡ π/45 M Pl ∆ σ 0 (g * S / √ g * ), and c 2 ≡ 45/(4π 3 g * )(α 5 µ M Pl /∆ 2 ). We take g * = g * S = 108.75 to account for the SM plus the two dark-photon degrees of freedom, and assume that g * and g * S remain constant.
We numerically integrate Eqs. (10)- (12), starting from z = z i = 5(α 2 /4), until z = z s = 100 z kd , using the initial conditions (14) Increasing z i up to z f and varying z s within reasonable limits has a negligible effect. The thermal equilibrium values for Y ↑↓ (z i ) and Y ↑↑ (z i ) are warranted because BSF gets into equilibrium before freeze-out for the couplings of interest; nevertheless, choosing instead vanishing initial values does not appreciably change the result. We checked that the BSF and ionization rates appearing in Eqs. (10) - (12) cancel each other when Y X , Y ↑↓ and Y ↑↑ are equal to their equilibrium values. The fractional DM relic density is Ω X = 2mY X (z s )s 0 /ρ c , where the factor of 2 accounts for the sum of X andX. ρ c 4.9 × 10 −6 GeV cm −3 and s 0 2795 cm −3 [34] are the critical energy density and the entropy density of the universe today.
We numerically evaluate Y X (z s ), and determine α as a function of m, such that the observed DM density, Ω X 0.26 [34], is reproduced. We present our results in Fig. 2, where we compare them with those obtained by neglecting BSF, and by neglecting both BSF and the SE of the 2-to-2 annihilations. As can be seen, the effect of BSF is significant for m few TeV. The relic DM density can also be estimated without employing the coupled differential equations (10) -(12), yet incorporating the effect of BSF. BSF contributes effectively to the DM annihilation, provided that the bound states which are formed decay into dark photons faster than the ambient radiation can reionize them into their constituents. We may thus define an effective thermally averaged SE factor (c.f. the timeline in table I) , 0.28 z and c α f ion (z) We can then estimate the DM relic abundance from the evolution equation for z z f . (We ignore the pair-creation processes, which become unimportant soon after freeze-out.) Equation (16) can be analytically integrated to give Using Eq. (17), we find α as a function of m such that the observed DM density is reproduced. The results are in agreement with those obtained from solving the coupled Boltzmann equations, to better than 2% accuracy.

V. UNITARITY AND CRITICAL COUPLING
For α 0.54, the sum σ ann +σ BSF becomes equal to the unitarity bound on the inelastic cross-section in Eq. (1). For this α, the observed DM abundance is reproduced if m equals This is the unitarity bound on the mass of thermal relic DM consisting of Dirac fermions. From this, we deduce the corresponding bound on Majorana DM, m M uni = √ 2 m D uni 197 TeV. The bounds for complex and real scalar DM are the same as for Dirac and Majorana DM respectively. Of course, if there is significant entropy release in the universe after DM annihilations become inefficient, or if the dark radiation is at a lower temperature than the SM plasma during freeze-out, the bounds are relaxed accordingly.
The preceding analysis takes into account the belated DM freeze-out and the DM annihilations occurring after that point, due to the SE of σ ann and σ BSF at low velocities. The SE results in m uni being larger than what would otherwise be expected. For comparison, using Ω DM h 2 0.12 [34], the analysis of Ref. [26] would give m D uni 340 TeV · (Ω DM h 2 /2) 1/2 83 TeV for annihilation without SE, where the factor 1/ √ 2 translates the bound from Majorana to Dirac DM. Reference [26] also considered the case of annihilation with SE (although they deemed it improbable). Their estimate for this case would be m D uni 550 TeV · (Ω DM h 2 /2) 1/2 135 TeV, which is consistent with Eq. (18). (Indeed, close to the unitarity limit, the effect of the ionization of the bound states is negligible.) We may estimate the α = α uni for which σ inel = σ inel, max . Assuming that σ inel = σ ann , α uni 0.86. For σ inel = σ ann + σ BSF , α uni 0.54. α uni provides an estimate of the range of validity of the approximation used in evaluating the inelastic cross-section, albeit not necessarily the most stringent. According to Gribov [35,36], gauge theories have a critical coupling above which the Coulomb interaction between fermions becomes strong enough to cause a rearrangement of the perturbative vacuum. In QED, this is α crit = π(1− 2/3) 0.58 [35][36][37][38]. It is intriguing that when taking into account both direct annihilation and BSF, α uni turns out to be very close to (but lower than) α crit . On the other hand, when considering 2-to-2 annihilation only, α uni appears to be larger than α crit . If the vacuum structure of the theory changes at α crit , then this indicates a priori that the direct annihilation cannot account for the total inelastic cross-section at such large couplings; indeed, BSF dominates.

VI. CONSTRAINTS
The model and the parametric regime investigated here, as well as other similar scenarios, are viable with respect to observational constraints. In particular, the coupling of DM to a light or massless particle mediates DM self-interactions. The most stringent bounds on these interactions arise from the observed ellipticity of Milky-Way-sized haloes. Although significant uncertainties exist, recent simulations and observations estimate the upper bound on the momentum-transfer scattering crosssection per mass of DM to be σ mt /m 2 barn/GeV [39][40][41][42]. In the present model, where θ min > 0 encodes the effect of screening due to the Debye length in neutral plasma and/or due to a finite mediator mass. Taking ln[csc(θ min /2)] ∼ 10, and v rel ∼ 250 km/s 8 × 10 −4 for a Milky-Way-sized halo, we estimate which satisfies the existing constraints in the mass range of interest. (Note that because of the SE of σ ann and of the effect of BSF on the DM freeze-out, the observed DM abundance is obtained for σ 0 < σ ann v rel c , for m TeV.) Additional constraints may arise if the dark sector couples to SM particles. Exploring these constraints is beyond the scope of the present work.

VII. CONCLUSION
We demonstrated that the formation of bound states in the early universe significantly enhances the annihilation rate of thermal DM with mass above a few TeV, and affects its relic abundance. We argued that DM can be as heavy as unitarity permits only if it annihilates via a long-range interaction; we determined the unitarity bound to be about 139 TeV for non-self conjugate DM, and showed that BSF is the dominant annihilation channel in this regime.
Here we focused on DM consisting of Dirac fermions annihilating into massless dark photons. We determined the coupling strength which yields the observed DM abundance; our results are presented in Fig. 2. We found that the required coupling α is up to about 40% smaller than estimated when BSF is ignored, with the largest discrepancy arising close to the unitarity limit. Our analysis by powers of α/(4π) which is always small for the range of α considered here. The same is true if we include direct annihilation to a larger number of photons, which is additionally suppressed by phase space. may be extended to other types of DM interacting via a light but massive scalar or vector boson. This is particularly compelling for heavy DM coupled to the weak interactions of the SM, which can be probed at the LHC and future high-energy colliders. It is also important for hidden-sector DM, which may yield observable highenergy astrophysical signals. Indeed, the accurate interpretation of the experimental results necessitates a precise knowledge of the couplings which yield the observed DM abundance.
In fact, the significance of these couplings is even broader. Thermalized species which annihilate more efficiently can account for the entirety of DM, provided that they carry a particle-antiparticle asymmetry [43][44][45][46][47]. Asymmetric DM provides a dynamical explanation for the similarity of the dark and the ordinary matter abundances [43]; furthermore, it is a very good host of self-interacting DM [48], which is favored by observations of the galactic structure [49]. On the other hand, DM which annihilates less efficiently may have been pro-duced only non-thermally, otherwise it would overclose the universe; candidates in this category are sterile neutrinos [50][51][52][53][54][55][56], and axions [57][58][59][60][61]. The above possibilities arise, of course, within DM theories of very different structure, i.e. very different beyond-SM physics. The precise value of the couplings which produce the observed DM abundance in the symmetric thermal-relic scenario sets the border in the parameter space between structurally different DM theories and is an important quantity for DM physics even beyond that specific paradigm.