Rapid bound-state formation of Dark Matter in the Early Universe

The thermal decoupling description of dark matter (DM) and co-annihilating partners is reconsidered. If DM is realized at around the TeV-mass region or above, even the heaviest electroweak force carriers could act as long-range forces, leading to the existence of meta-stable DM bound states. The formation and subsequent decay of the latter further deplete the relic density during the freeze-out process on top of the Sommerfeld enhancement, allowing for larger DM masses. While so far the bound-state formation was described via the emission of an on-shell mediator ($W^{\pm}$, $Z$, $H$, $g$, photon or exotic), we point out that this particular process does not have to be the dominant scattering-bound state conversion channel in general. If the mediator is coupled in a direct way to any relativistic species present in the Early Universe, the bound-state formation can efficiently occur through particle scattering, where a mediator is exchanged virtually. To demonstrate that such a virtually stimulated conversion process can dominate the on-shell emission even for all temperatures, we analyze a simplified model where DM is coupled to only one relativistic species in the primordial plasma through an electroweak-scale mediator. We find that the bound-state formation cross section via particle scattering can exceed the on-shell emission by up to several orders of magnitude.


I. INTRODUCTION
One of the leading dark matter hypothesis is that it consists of Weakly-Interacting-Massive-Particles (WIMPs) [1][2][3], which can explain the observed DM abundance [4] in a natural way through the thermal production mechanism. While strong upper bounds on the coupling strength of the WIMP to Standard Model (SM) particles, derived from direct, indirect and collider searches, rule out many electroweak mass-scale realizations in the thermal production scenario, the TeV mass region and above still remains an attractive and much less constrained possibility.
The heavier WIMPs are, the more important it is to include quantum mechanical effects induced by long-range interactions in the thermal decoupling description for predicting the relic abundance accurately. In a seminal work [5], it has been pointed out that already around a three TeV mass, it is possible that even the heaviest massive gauge bosons of the SM effectively act as attractive long-range forces between a WIMP pair, leading to an enhanced annihilation cross section during the freeze-out process. This so-called Sommerfeld enhancement (SE) [6] or Sakharov enhancement [7,8] lowers the predicted thermal relic abundance compared to a tree-level computation, which in turn allows for larger masses of the annihilating DM particles to compensate for the effect (see, e.g., [9][10][11][12][13][14][15][16][17] and [18][19][20][21][22][23][24] for formal aspects).
So far, the DM bound-state formation process in the Early Universe was described via the emission of an onshell mediator. While this resembles the situation of SM neutral hydrogen recombination in the matter-dominated epoch, an interesting question is by how much scattering events with primordial plasma constituents in the radiation dominated epoch could additionally stimulate the conversion process through virtual mediator exchanges.
In the case of annihilating heavy quarkonia, produced in a quark-gluon plasma (QGP) in heavy ion collision at, e.g., the Large-Hadron-Collider, it is well known that the dissociation of the heavy quarkonia bound states via the absorption of an on-shell gluon is not the dominant process for temperature larger than the binding energy (see, e.g., [58,59]). Instead it is the dissociation via light quark and gluon scattering which dominates in that tempera- Figure 1. Feynman diagrams for bound-state formation via on-shell mediator emission (left) and via bath-particle scattering (right) shown. Both processes contain also the diagram where the mediator is attached toχ.
ture regime. By the argument of detailed balance, this automatically implies that the dominant heavy quarkonia BSF channel must be via light parton scattering as well. Since the system of heavy quarks inside a QGP is similar to DM in the Early Universe (where the latter contains many more relativistic species), this insight might have already profound implications for DM models where light or massless mediators in co-annihilation scenarios are involved.
In this work, we investigate a mediator whose mass is of the order of the Z, W ± , or Higgs boson mass. The latter cases might be of particular interest, since the kinetic energy of WIMPs during the freeze-out process could be not enough to emit an on-shell massive boson at the electroweak scale in order to form a bound state. For such a case, it is known that BSF has only marginal effects due to the kinematical block. However, inside a relativistic plasma background, bath-particle scattering can stimulate the BSF process by inducing a mediator excitation virtually. The virtual exchange implies that BSF via bath-particle scattering has no kinematical block and might entirely dominate over the on-shell mediator emission.
As a proof of concept, this article investigates for the first time this possibility in a simplified DM model, which is aimed to parametrically resemble more realistic cases with Z, W ± , H or other exotic interactions. The details of the model and the computation based on a generalized BSF cross section [60], capturing higher-order processes in a proper thermal field theory framework, are shared in section II. The numerical evaluation of the thermally averaged quantities and the implications of strongly enhanced bound-state formation rates are discussed in section III. Section IV concludes this work.

II. MODEL AND CROSS SECTIONS
We consider a rather simple three-parameter model: where non-relativistic DM (χ) with mass m χ and the ultra-relativistic primordial plasma particles (b) are Dirac Fermions. The mass of the Abelian vector mediator is assumed to fulfill m V αm χ , where α ≡ g 2 /(4π) is the fine structure constant, such that bound-state solutions exist in the DM two-particle spectrum.
The generalized bound-state formation cross section, which includes the on-shell emission and the bathparticle scattering in fig. 1, as well as other higher order processes in a proper thermal field theoretical framework, is given by [60] (see also App. A): is the four momentum of the vector boson V , whose energy component is fixed by the positive quantity ∆E ≡ E k − E nlm . ∆E is the total energy emitted in the inelastic conversion, i.e., relative kinetic energy E k = k 2 /(2µ) = µv 2 rel /2 of DM with reduced mass µ, plus the absolute value of the negative binding energy E nlm of the bound state with quantum numbers nlm. The central Eq. (2) consist of a two-point correlation function D −+ µν of the vector boson, which is contracted with non-relativistic (NR) scattering-bound state transition matrix elements at the Born-level. The integral over p picks out the various physical poles appearing in D −+ µν . The matrix elements are defined in dipole approximation (dip) as: where M is the bound-state mass and initial DM spin d.o.f. are g χ gχ = 4. M µ is defined as the Fourier transform of the Scattering-Bound state transition matrix element of the current operator, here given by: / δ 4 = (2π) 4 δ 4 is the four-momentum conserving delta function. Spin indices are implicit, see App. C for our convention of the states. For the rest of this work, the dominant direct capture into the ground state is considered. Contracting the spinors in Eq. (4), taking the non-relativistic and dipole limit, we obtain [61]: with ∆E = E k + |E 100 | for the ground state, unit relativemomentum vectork, and we defined ζ ≡ α/v rel and ξ ≡ αm χ /(2m V ). Note that once the Born transition matrix element Eq. (4) is determined [e.g., Eq. (5)] it can be reused in Eq.
(2) at any order of the perturbative expansion of the two-point correlation function. As a check, T fulfills the current conservation P µ T µ k,100 = 0 as a consequence of the global symmetry. Dimensionless h(ζ, ξ) contains the dipole approximation of the overlap integral, which can be computed only numerically for a massive mediator [61], see App. C for convention.
The two-point correlation function in coordinate space encodes all interactions with the primordial plasma environment and hence the information of the bath-particles. It is related to the spectral function D ρ by the Kubo-Martin-Schwinger relation (see, e.g., [62]): The equilibrium phase-space distribution f eq V obeys Bose-Einstein statistics. It is convenient to compute the spectral function from the retarded correlator via D ρ µν = 2 iD R µν , since the latter obeys also in thermal field theory the Dyson-Schwinger equation, given in momentum space by In the following, we show that BSF via on-shell mediator emission is reproduced from the free term, while BSF via bath-particle scattering is contained in the interaction term with the retarded self-energy Π R . For both processes, Eq. (5) can be reused to calculate the BSF cross section Eq. (2), thanks to the factorization. Although a self-energy with fermions is analyzed in this work, the formalism can capture non-abelian interactions, e.g. gluon scattering (triple vertex). For non-Abelian theories, ref. [35] provides many expressions of the transition element Eq. (3), ready to be used with Eq. (2). For scalar mediators, one can drop "µν" everywhere and adjust Eq. (4) to Yukawa interactions.

A. Capture via on-shell mediator emission
To the lowest order in perturbation theory the retarded propagator of the massive vector mediator is given by ..] of the latter to compute the spectral function, and inserting Eq. (6) together with Eq. (5) into the cross section, the standard result for capture into the ground state via the on-shell emission of a massive vector mediator (mBSF) is reproduced after performing elementary integrals (cf. Eq.3.7a in [61]): where s ps ≡ 1 − m 2 V /∆E 2 implies the kinematical suppression.

B. Capture via bath-particle scattering
The retarded self-energy for ultra-relativistic fermionic particles is given by (see App. B): The on-shell integration over the ultra-relativistic momenta is dΠ i = d 3 pi (2π) 3 2|pi| , and P i = (|p i |, p i ). The second line contains Fermi-Dirac distributions. The term shown corresponds to bath-particle scattering. The omitted term contains an off-shell decay process of the vector mediator into a bath-particle pair, which we have found to be sub-dominant compared to the scattering case. By taking the imaginary part of the propagator in the second line to compute the spectral function, and together with the transition elements in Eq. (5), we have for the first time computed the DM bound-state formation cross section via bath-particle scattering (bBSF) for a massive vector mediator: The angular average arises in the thermal average of the cross section. We find it here for a computational reason meaningful to already perform. The first term is identical to the one also appearing in the case of mBSF, see Eq. (8). R bBSF 100 encompasses the details of the bathparticle scattering and is given in dimensionless coordinates y ≡ |p 1 |/T by: The amount of integrals were reduced analytically down to only one single remaining, without taking further approximations. The factor 2 originates from summing over particle and anti-particle scattering, already contained in Eq. (9). We would like to stress the fact that bBSF has no kinematical block, and has an additional temperature and velocity-enhancement factor (T /∆E) 3 for T |E 100 |. It was explicitly checked that Eq. (10) is identical to the expression one would obtain in the Boltzmann formalism. The benefit of starting from a proper thermal field theory definition Eq. (2) is that the contribution of additional interference terms is automatically taken into account. In the Boltzmann framework they are absent and arise from taking instead the imaginary part of the double mediator propagator in Eq. (7). For our model, we find that their contribution for m V |E 100 | is negligible. However, notice that the additional interference terms regulate Eq. (11) in the limit m V → 0, by canceling the forward scattering divergence of the bathparticles, as mathematically proven in ref. [60]. In the latter work, it is also shown that the UV-divergent vacuum parts in Eq. (9) do not contribute significantly after standard renormalization. This shows the power of Eq. (2). It contains all possible processes up to a certain order in the coupling expansion.

III. NUMERICAL RESULTS AND IMPLICATIONS
The thermal average of the BSF cross sections via onshell mediator emission Eq. (8) and via bath-particle scattering Eq. (10), as well as of the Sommerfeld enhanced annihilation cross section (App. D) are compared in fig. 2. Two different cases are shown where in both m V |E 100 |. In the left panel, only one bound state close to threshold exists, aiming to parametrically resemble Wino-neutralino DM scenarios. In the right panel, the coupling has similar values as the strong coupling of the SM. The ground state is nearly Coulombic, aiming to provide insight into, e.g., co-annihilation with color-charged particles (neutralino-squark) and multi-TeV WIMPs, typically residing in higher SU (2) L representations. Overall one can immediately recognize that the results show the proof of concept, i.e., bBSF can entirely dominate over mBSF. In sharp contrast to mBSF, bBSF is kinematically accessible even if m V ∆E. The reason why bBSF is strongly enhanced compared to mBSF, can be qualitatively explained as in the following. First, notice from Eq. (10) that bBSF is enhanced at high temperatures by the large number density of the ultra-relativistic bath particles n b ∝ T 3 . Second, while mBSF is suppressed by the smallness of the dissipated energy (σ mBSF v rel ) ∝ |p| = ∆E [61], bBSF is enhanced by it: (σ bBSF v rel ) ∝ 1/∆E 2 has a stronger dependence on the inverse relative velocity. One reason for this is that the momentum of the off-shell V µ is not fixed by ∆E, as is in mBSF. These features yield the relative scaling R bBSF ∝ 1/∆E 3 . An important similarity, however, is that both processes fulfill the same angular momentum selection rules, seen from the factorization in Eq. (2).
Hence, forbidden channels can also not be accessed via higher order processes contained in the mediator spectral function.
The implications of a large bBSF rate on the evolution of the DM density can be seen from the Boltzmann equation, including the lowest bound states as [28,29]: .
The depletion of the DM number density depends on the SE annihilation cross section and the effective cross section W . The latter quantity stores all the information of the bound states and consists of the total thermally averaged BSF cross section, weighted by branching fractions containing the decay rate of spin-singlet and triplet bound states into two V and a pair of bath particles, respectively (see App. D). The factors 1/4 and 3/4 account for the spin multiplicities. The dissociation rate is related via detailed balance to the BSF cross section as Γ dis 100 = σ BSF 100 v rel (n eq χ ) 2 /n eq 100 , where n eq 100 is the ground state equilibrium number density with four spin d.o.f.
The effective cross section has two asymptotic regimes: where Γ dec 100 = (Γ dec 100,S + 3Γ dec 100,T )/4. A large BSF cross section implies that Γ dec 100,i Γ dis 100 at early times. In this first asymptotic regime, the system is in a phase of Saha ionization equilibrium (io. eq.), where W is i) independent of the BSF cross section [62], and ii) maximal for a given temperature and bound-state decay rate. The end of the initial io. eq. phase is marked by the green lines in Fig. 2, where Γ dec 100,i = Γ dis 100 for the spin-singlet and triplet states in our model. Importantly, this transition depends on the value of total BSF cross section contained in Γ dis 100 . In the second asymptotic regime at later times, the number density depletion depends on the actual BSF cross section, since the bound states immediately decay without being dissociated back into the scattering states.
From this discussion one sees that an additional BSF channel can significantly increase the total depletion. First, a large bBSF contribution implies that io. eq. is sustained until much lower temperatures than if mBSF was the only capture process. Secondly, bBSF can significantly enhance or entirely dominate the total BSF rate after the end of io. eq., even if briefly so, as seen in Fig. 2.
The dipole approximation for bBSF breaks down at T larger than the typical momentum inside the bound state, marked by the gray line in Fig. 2. However, since the system is robustly in io. eq., the actual magnitude of the BSF cross section is irrelevant. It is possible though, that the thermal bath at high T affects the SE annihilation cross section and the bound state decay rate, which are the only quantities that determine the collision term in Eq. (12) in io. eq. The impact of the primordial plasma environment on the SE annihilation and bound-state decay rate is already well studied within a different thermal field theory approach [33,36,37,39,41,[62][63][64][65][66]. In particular, Eq. (7) can be resummed in the Hard-Thermal-Loop approximation [67], giving an effective in-medium potential for a Schrödinger-like equation, including real part corrections (e.g., Debye screening), as well as an imaginary thermal width [68,69]. A most conservative estimate [63] shows that the latter leads to an entire melting of the bound states (see, e.g., fig. 1 in [70] for data) for T to the left of the gray line, where ultra-efficient bath-particle scattering strongly mixes the bound and scattering states.
The formal key-point of this work is that the generalized BSF cross section Eq. (2) [60] has overlap with the validity region of the formalism used in [33,36,37,39,41,[62][63][64][65][66], where the latter is strictly speaking limited to io. eq. [62,66]. Combined, we have achieved a complete procedure to accurately calculate thermal relics in the Early Universe for any temperature regime, ranging from the melting of bound states at high T down to far below their decoupling from ionization equilibrium.

IV. CONCLUSION
The Early Universe consists of more than about 100 relativistic d.o.f., which can in principle all contribute to the DM bound-state formation process via particle scattering. In this work, we have analyzed the contribution of only one single species of relativistic fermions, coupled to DM through a mediator at the electroweak mass scale. While in the standard on-shell emission computation, the BSF rate can be highly suppressed for such a case, we have demonstrated that BSF via bath-particle scattering can be highly active. In addition, it was emphasized that the latter can lead to ionization equilibrium, where the effective depletion cross section is maximized for a given temperature. The implications of our results are from the broader literature perspective (based on the on-shell emission), that bound-state effects in WIMP models can become more pronounced during the thermal freeze-out process. Consequently, we conclude that DM could be heavier than previously expected, eventually informing indirect and collider search strategies.
This appendix complements ref. [60], where the generalized BSF cross section Eq. (2) is derived from pNRQFT in the framework of non-equilibrium quantum field theory. Here, we start from the more commonly known Boltzmann framework and derive Eq. (2) at the free order in the two-point correlator, where thermal field theory and the Boltzmann approach can be brought together.
In the standard Boltzmann framework, the cross section for describing the conversion of a χχ pair into a bound state with quantum numbers nlm via the emission of an on-shell boson φ (here general), is given by: The thermal average for all cross sections in our work follows the standard convention (...) = (...) if (...) is velocity independent. Final state integration is over the emitted boson four momenta P = (P 0 , p), and over the dark matter bound state four momentum Q = (Q 0 , q).
The four momenta of the initial DM two-body state, as well as of the final bound state with mass M are described non-relativistically. The slashed convention is / δ n = (2π) n δ n . The squared amplitude |M| 2 k,nlm for the bound-state formation process in Born approximation is summed over initial and final state internal degrees of freedom.
To allow for general dispersion relations of the mediator at the end, the on-shell momentum integration is first rewritten into four-momentum integration as One can recognize that the mediator and bound-state dispersion relation are here on the mass shell. Within this notation the BSF cross section can be written as: Here, we defined N = (g χ gχ4m 2 χ 2M ) 1/2 , see Eq.
(3), and used the four-momentum conserving delta function in Eq. (A1) to perform the Q integration. Furthermore, we switched to the DM CM-momenta coordinates where K ≡ k χ + kχ and k ≡ (k χ − kχ)/2 with |k| = µv rel . In the zero momemtum-recoil approximation, the delta function can be estimated as: 1 and the BSF cross section can be brought into the form: Here, the theta function θ(P 0 ) was dropped, since the delta function in the first line already ensures positivity of the emitted energy. Assuming a massive vector field, one can decompose the amplitude as where is the polarization vector. Note that in our model, M µ k,nlm here is identical to what has been defined already in Eq. (4). By factoring out the polarization tensor dependence from the squared matrix element in Eq. (A6), one may notice the lesser two-point correlation function in thermal equilibrium Eq. (6) in the free limit: After integrating Eq. (A6) over P 0 , and replacing the free lesser function by interacting one, Eq. (2) is recovered. 1 In the zero-momentum-recoil approximation it is assumed that K p, which is justified since typical DM CM-momentum K (∼ mχT ) is much larger than the momentum of the emitted particle p (∼ ∆E ∼ Max[T, |E 100 |]) in the bound-state formation process. Already around the freeze-out temperature, where DM enters the non-relativistic regime this holds, and for temperature larger than ∼ α 2 |E 100 |. This captures the entire region of the freeze-out we are interested in, since DM decouples from the bound states when temperature is already less than the ground state energy, i.e., T ∼ |E 100 | α 2 |E 100 |. For electroweak couplings we have checked the zero-recoil approximation in the high temperature regime mχ/T = 10 and noticed that the difference is about only one per-cent compared to the exact result (but still assuming dipole approximation). Also note that ref. [60] assumes the zero-momentum-recoil approximation to be valid from the beginning, since the derivation is based on pNRQFT.
The justification of this replacement from first principles can be found in ref. [60]. Here, we would like to argue from the Kadanoff-Baym equation that this "upgrade" is expected. As derived in ref. [71], the DM collision term in the forward direction can be determined through the DM self-energy component Σ −+ . In the 1-PI formalism, the Σ −+ contains interacting D −+ µν . This is another indication that the upgrade from vacuum correlations in the Boltzmann formalism to the interacting correlation functions in the thermal field theory approach is indeed justified. Finally we would like to remark that this scheme of derivation presented in this Appendix also holds for scalar mediators, where at the end, one has to just drop the Greek indices in Eq. (2).

Appendix B: Retarded self-energy expression
The retarded self-energy is defined in terms of greater and lesser self-energies as where the two-point functions of the fermionic bathparticles are and their free Fourier transforms are in equilibrium The Fourier transform of the retarded self-energy is then The first and second term inside the brackets corresponds to particle and anti-particle scattering, when taking the imaginary part of the propagator. Third and fourth line corresponds to off-shell decay into a bath-pair, as well as the reverse process. The latter does not contribute, since in our case P 0 = ∆E is positive, leading to the fact that the cut of the propagator is zero. The following identity was used for computing Eq. (11): which makes the particle scattering more apparent. The Bose-enhancement factor comes from Eq. (6). Similar identities can also be derived for the off-shell decay.
Appendix C: States convention and overlap integral In terms of particle and anti-particle creation operators, the states in Eq. (4) follow the convention: Q is the momentum of the bound and K is the CMmomentum of the scattering state. Due to the spinsum in Eq. (2), we consider for simplicity an unpolarized BSF cross section which is however not necessary to do. The mass norms from the states cancel out in Eq. (3). The reduced spin-independent two-body wave functions encode the plane-wave distortion due to the long-range interactions. They have energy Eigenvalues E k = k 2 /(2µ) = µv 2 rel /2 and negative binding energy E nlm , and are normalized to plane waves in the case of scattering states, while for the bound states the spatial integral over the absolute squared wave function is normalized to unity. The scattering state is a coherent sum over all partial waves.