IPMU19-0148
DESY 19-181

Rapid bound-state formation of Dark Matter in the Early Universe
Tobias Binder,1, ∗ Kyohei Mukaida,2, † and Kalliopi Petraki3, 4, ‡1

Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Deutsches Elektronen-Synchrotron (DESY), Notkestraße 85, Hamburg, D-22607 Germany
3
Sorbonne Université, CNRS, Laboratoire de Physique Théorique et Hautes Energies, LPTHE, F-75252 Paris, France
4
Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands(Dated: October 25, 2019)arXiv:1910.11288v1 [hep-ph] 24 Oct 20192

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 boundstates. The formation and subsequent decay of the latter further deplete the relic density during thefreeze-out process on top of the Sommerfeld enhancement, allowing for larger DM masses. Whileso far the bound-state formation was described via the emission of an on-shell mediator (W ± , 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 avirtually 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 isthat it consists of Weakly-Interacting-Massive-Particles(WIMPs) [1–3], which can explain the observed DM
abundance [4] in a natural way through the thermal
production mechanism. While strong upper bounds onthe 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 muchless 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 forpredicting the relic abundance accurately. In a seminal
work [5], it has been pointed out that already around athree 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 theannihilating DM particles to compensate for the effect(see, e.g., [9–17] and [18–24] for formal aspects).
Further quantum mechanical effects caused by attractive long-range interactions are bound-state solutions in∗ tobias.binder@ipmu.jp† kyohei.mukaida@desy.de‡ kpetraki@lpthe.jussieu.fr

the two-particle spectrum of WIMPs [25–27]. The formation of these bound states and their subsequent decay
into SM particles additionally depletes the relic density,
allowing for even heavier masses [28]. For every SM forcecarrier, possible DM scenarios have been found in theliterature where the inclusion of the SE and bound-state
formation (BSF) is relevant. Famous examples in Supersymmetric extensions of the SM are cases where DM coannihilates with closely mass degenerate partners, which
are either electroweakly or color charged [29–38]. In addition to the photon, gluon, Z or W ± boson induced
bound states in mentioned cases, also the Higgs particle [39, 40] can attractively contribute to confine DM intoa meta-stable bound state, even for otherwise repulsive
color octet states [41]. Furthermore, these quantum mechanical phenomena certainly play a role in bottom-up
motivated [42] scenarios, e.g., Self-Interacting DM [43]with light mediators [44–53], where long-range interactions are introduced by assumption.
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-dominatedepoch, an interesting question is by how much scattering
events with primordial plasma constituents in the radiation dominated epoch could additionally stimulate theconversion 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 theabsorption of an on-shell gluon is not the dominant process for temperature larger than the binding energy (see,e.g., [54, 55]). Instead it is the dissociation via light quarkand gluon scattering which dominates in that tempera-

2
Vµ

b(b̄)

χ
χ̄

χ
···

···

B

χ̄

···

···

B

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, thisautomatically 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 insightmight 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 isof 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 sucha case, it is known that BSF has only marginal effects
due to the kinematical block. However, inside a relativistic plasma background, bath-particle scattering canstimulate the BSF process by inducing a mediator excitation virtually. The virtual exchange implies that BSF
via bath-particle scattering has no kinematical block andmight 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 caseswith Z, W ± , H or other exotic interactions. The detailsof the model and the computation based on a generalized
BSF cross section [56], capturing higher-order processes
in a proper thermal field theory framework, are sharedin section II. The numerical evaluation of the thermallyaveraged 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:
L ⊃ −g χ̄γ µ χVµ − g b̄γ µ bVµ ,(1)where non-relativistic DM (χ) with mass mχ and theultra-relativistic primordial plasma particles (b) are DiracFermions. The mass of the Abelian vector mediator isassumed to fulfill mV . αmχ , where α ≡ g 2 /(4π) is thefine 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 orderprocesses in a proper thermal field theoretical framework,
is given by [56] (see also App. A):
Z 3
X µd p −+BSF
ν?
σnlm vrel =Tk,nlm (P )Tk,nlm(P ). (2)
D
(P
)
µν
(2π)3
spins

P = (P 0 = ∆E, p) is a four momentum, whose energy component is fixed by the positive quantity ∆E ≡Ek −Enlm . ∆E is the total energy emitted in the inelastic
conversion, i.e., relative kinetic energy Ek = k2 /(2µ) =2µvrel/2 of DM with reduced mass µ, plus the absolutevalue of the negative binding energy Enlm of the boundstate with quantum numbers nlm. The central Eq. (2)−+consist of a two-point correlation function Dµνof themassive vector field V , which is contracted with nonrelativistic (NR) scattering-bound state transition matrix elements at the Born-level. The matrix elements aredefined in dipole approximation (dip) as:NRµTk,nlm(P ) ≡ (gχ gχ̄ 4m2χ 2M )−1/2 Mµk,nlm dip ,(3)

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:Z
4/δ Mµk,nlm = d4 x eiP x hBnlm |g χ̄(x)γ µ χ(x)|Sk i.
(4)4/δ = (2π)4 δ 4 is the four-momentum conserving delta
function. Spin indices are implicit, see App. C for ourconvention of the states. Note that for non-Abelian theories, ref. [35] provides already many expressions of thistransition element, ready to be used with Eq. (2).
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 [57]:

µ 
ν
X µ
4πh(ζ, ξ) p · k̂p · k̂ν?
, (5)Tk,100 Tk,100= 3µ ∆E 2 ∆E k̂
∆E k̂
spins
with ∆E = Ek + |E100 | for the ground state, unit relativemomentum vector k̂, and we defined ζ ≡ α/vrel andξ ≡ αmχ /(2mV ). 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µ Tk,100=0
as a consequence of the global gauge symmetry. Dimensionless h(ζ, ξ) contains the dipole approximation of theoverlap integral, which can be computed only numericallyfor a massive mediator [57], see App. C for convention.
The two-point correlation function in coordinate space−+Dµν(x, y) ≡ hVµ (x)Vν (y)i, where h...i = Tr[e−Henv /T ...],

3
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 KuboMartin-Schwinger relation (see, e.g., [58]):−+
ρDµν(∆E, p) = [1 + fVeq (∆E)] Dµν(∆E, p).

(6)

The equilibrium phase-space distribution fVeq obeys BoseEinstein statistics. It is convenient to compute the specρ
tral function from the retarded correlator via Dµν= R
2= iDµν , since the latter obeys also in thermal field
theory the Dyson-Schwinger equation, given in momentum space byRR,0R,0 αβ R,0
Dµν
= Dµν
+ Dµα
ΠR Dβν + ... .

(7)

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 we analyze a particular self-energy in the following, where only the fermionic bath-particles appear,
it should be clear that also non-abelian interactions likegluon scattering (triple vertex) can be captured within
this formalism. For studying scalar mediators like the
Higgs Boson, one can just 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
R,0(P ) = −igµν [P 2 − m2V + isgn(P 0 )]−1 . Taking
by Dµν
2=[i...] of the latter to compute the spectral function,and inserting Eq. (6) together with Eq. (5) into the crosssection, the standard result for capture into the groundstate via the on-shell emission of a massive vector mediator (mBSF) is reproduced after performing elementary
integrals (cf. Eq.3.7a in [57]):
1/2

mBSF
σ100
vrel =

4h(ζ, ξ)∆E
sps (3 − sps ), (8)[1 + fVeq (∆E)]
3
3µ2

where sps ≡ 1 − m2V /∆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):ZR/ 1 γν P/ 2 ]×Πµν(∆E, p) = g 2 dΠ1 dΠ2 Tr[γµ P3i2/δ (p + p1 − p2 )eqeq[f (|p1 |)−fb (|p2 |)]+... . (9)∆E + |p1 | − |p2 | + i bThe on-shell integration over the ultra-relativistic mo3imenta is dΠi = (2π)d 3p2|p, and Pi = (|pi |, pi ). The seci|
ond line contains Fermi-Dirac distributions. The termshown 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 foundto be sub-dominant compared to the scattering case. Bytaking 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 crosssection via bath-particle scattering (bBSF) for a massive
vector mediator:ZdΩk bBSF
4h(ζ, ξ)∆E
bBSF
× R100.
(10)(σvrel ) =4π 1003µ3
The angular average arises in the thermal average of thecross section. We find it here for a computational reason meaningful to already perform. The first term isidentical to the one also appearing in the case of mBSF,
bBSF
see Eq. (8). R100
encompasses the details of the bathparticle scattering and is given in dimensionless coordinates y ≡ |p1 |/T by:


3Z ∞

α
T
1
1
bBSF
R100
=2×
dy y
1 − ∆E/T +y
2π ∆E
e +1
e
+1
0
("
# 

2

∆E
m2
4y(∆E/T + y)×y2 ++ y + 2V ln 1 +TT
m2V /T 2)
4y(∆E/T + y) m2V /T 2 + (∆E/T + 2y)2−.
(11)
m2V /T 2 + 4y(∆E/T + y)
The amount of integrals were reduced analytically down
to only one single remaining, without taking further approximations. The factor 2 originates from summing overparticle and anti-particle scattering, already containedin 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 forT  |E100 |. 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 automaticallytaken into account. In the Boltzmann framework they areabsent and arise from taking instead the imaginary partof the double mediator propagator in Eq. (7). For ourmodel, we find that their contribution for mV  |E100 |is negligible. However, notice that the additional interference terms regulate Eq. (11) in the limit mV → 0, bycanceling the forward scattering divergence of the bathparticles, as mathematically proven in ref. [56]. In thelatter 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, where the dominant one
can be model dependent.

4
Hence, forbidden channels can also not be accessed via
higher order processes contained in the mediator spectral
function.
The implications of a strongly enhanced bBSF rate on
the DM number density evolution can be qualitatively
understood by inspection of the Boltzmann equation,
given for one existing bound state by [28, 29]:


eq 
ṅχ + 3Hnχ = − hσ an vrel i + W eff nχ nχ̄ − neqχ nχ̄ , (12)

102

10

1

W eff (T ) ≡

10- 1

10

102

103

10

102

103

Figure 2. Comparison of thermally averaged cross sectionsfor SE annihilation (blue), ground-state capture via on-shell
mediator emission (yellow) and via bath-particle scattering(red) shown. All hσvrel i are normalized to (σvrel )0 ≡ πα2 /m2χ .

III.

NUMERICAL RESULTS ANDIMPLICATIONS

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 comparedin fig. 2. Two different cases are shown where in bothmV  |E100 |. In the left panel, only one bound stateclose to threshold exists, aiming to parametrically resemble Wino-neutralino DM scenarios. In the right panel,
the coupling has similar values as the strong couplingof the SM. The ground state is nearly Coulombic, aiming to provide insight into, e.g., co-annihilation withcolor-charged particles (neutralino-squark) and multiTeV WIMPs, typically residing in higher SU (2)L representations. Overall one can immediately recognize thatthe 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 mV & ∆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 theultra-relativistic bath particles nb ∝ T 3 . Second, while
mBSF is suppressed by the smallness of the dissipatedenergy (σ mBSF vrel ) ∝ |p| = ∆E [57], bBSF is enhanced
by it: (σ bBSF vrel ) ∝ 1/∆E 2 has a stronger dependence on
the inverse relative velocity. One reason for this is thatthe momentum of the off-shell V µ is not fixed by ∆E,
as is in mBSF. These features yield the relative scaling
RbBSF ∝ 1/∆E 3 . An important similarity, however, isthat both processes fulfill the same angular momentum
selection rules, seen from the factorization in Eq. (2).

dec BSF
Γ100
hσ100 vrel i.
BSF v ineq neq /neqdecΓ100 + hσ100
rel χ χ̄100

(13)

The effective depletion cross section W eff consists of thedecbound-state decay rate Γ100
, total thermally averaged
BSF
BSF cross section hσ100 vrel i and the dissociation rate
eq
eq
dis
BSF
Γ100
≡ hσ100
vrel ineqχ nχ̄ /n100 , which is the inverse process and via detailed balance directly related to total BSF
cross section.
A strongly enhanced bBSF cross section implies thatdec
dis
the ratio Γ100
/Γ100
is almost zero. In the limit, the system is in a phase of Saha ionization equilibrium (io. eq.),
where the effective depletion cross section W eff is i) independent of the total BSF cross section and dissociation
rate [58], and ii) its value is maximal for a fixed temperature and a given set of parameters. The differenceof W eff to its maximum value is for a fixed temperatureespecially large if the dissociation rate is comparable orless than the decay rate as, e.g., in the left panel of fig. 2
when considering mBSF only. We mark the critical transition point for bBSF by the green lines. To the right,
the decay rate exceeds the dissociation rate, DM starts todecouple from io. eq., and the comoving number density
gradually approaches a constant value.
The limitation of the cross section Eq. (2) is marked by
the gray line, where T is comparable to the typical momentum inside the bound state. At higher T , first of all
the dipole approximation breaks down. However, since
the system is robustly in io. eq., Eq. (12) only depends
on the SE annihilation and the decay rate [58], hence theactual size of BSF cross section is irrelevant. The impactof bath-particle scattering on the SE annihilation andbound-state decay rate in such a high temperature limit
is already well studied [33, 36, 37, 40, 58–62]. I.e., Eq. (7)can be resummed in the Hard-Thermal-Loop approximation [63], giving an effective in-medium potential for the
Schrödinger-like equation, including real part corrections(e.g., Debye screening), as well as an imaginary thermal
width [64, 65]. A most conservative estimate [59] showsthat the latter leads to an entire melting of the boundstate (see, e.g., fig.1 in [66] for data) for T above thegray line, where ultra-efficient bath-particle scatteringstrongly mixes the bound and the scattering states. The
formal key-point of this work is that the generalized BSF
cross section Eq. (2) [56] has overlap with the validity region of the formalism used in [33, 36, 37, 40, 58–62], where
the latter is strictly speaking limited to io. eq. [58, 62].
Combined, we achieved a complete procedure to accurately calculate thermal relics in the Early Universe for5
any temperature regime, ranging from the melting ofbound states at high T down to far below their decoupling from ionization equilibrium.
IV.

CONCLUSIONThe 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 contributionof 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, wehave demonstrated that BSF via bath-particle scattering
can be highly active. In addition, it was emphasized thatthe latter can lead to ionization equilibrium, where theeffective depletion cross section is maximized for a given
temperature. The implications of our results are fromthe broader literature perspective (based on the on-shellemission), that bound-state effects in WIMP models canbecome more pronounced during the thermal freeze-out
process. Consequently, we conclude that DM could beheavier than previously expected, eventually informing
indirect and collider search strategies.
ACKNOWLEDGEMENTS

T.B. was supported by World Premier International
Research Center Initiative (WPI), MEXT, Japan. K.P.
was supported by the ANR ACHN 2015 grant (“TheIntricateDark” project), and by the NWO Vidi grant “Selfinteracting asymmetric dark matter”. We are thankful
to Burkhard Blobel, Peter Cox, Julia Harz and Shigeki
Matsumoto for very stimulating discussion, as well as toSatoshi Shirai for checking the manuscript.
Appendix A: Bound-state formation cross section

This appendix complements ref. [56], 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 orderin 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 abound state with quantum numbers nlm via the emission of an on-shell boson φ (here general), is given by:
Z
1
d3 p
d3 q
BSF
(σnlm vrel ) =4gχ gχ̄ m2χ (2π)3 2Ep (2π)3 2M
hi
4
× 1 + fφeq (Ep ) /δ (Kχ + Kχ̄ − P − Q)|M|2k,nlm . (A1)

The thermal average for all cross sections in our work
follows the standard convention h(...)i = (...) if (...) is
velocity independent. Final state integration is over theemitted boson four momenta P = (P 0 , p), and over thedark matter bound state four momentum Q = (Q0 , 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 isn/δ = (2π)n δ n . The squared amplitude |M|2k,nlm for thebound-state formation process in Born approximation issummed 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 firstrewritten into four-momentum integration as
Z
h
i
d3 p
eq
1
+
f
(E
)
=
p
φ
(2π)3 2Ep
Z 4 h
i

d P
eq
0
1+f
(P
)
θ(P 0 )sgn(P 0 )/δ(P 2 − m2φ ) , (A2)
φ
(2π)4
Z
d3 q
=
(2π)3 2M


Z
d4 Q
1
0
2
/
δ
(Q
−
M
−
|Q|
/(2M
))
.
(A3)
(2π)4 2M
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:
1
BSF
(σnlm
vrel ) = 2
N
Z 4 h
i

d P
eq
0
×
θ(P 0 )sgn(P 0 )/δ(P 2 − m2φ )
1
+
f
(P
)
φ
4
(2π)


K2
(K − P)2× /δ P 0 − ∆E −+|M|2k,nlm . (A4)4mχ4mχHere, we defined N = (gχ gχ̄ 4m2χ 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 whereK ≡ kχ + kχ̄ and k ≡ (kχ − kχ̄ )/2 with |k| = µvrel .
In the zero momemtum-recoil approximation, the delta
function can be estimated as: 1(K − P)2
K2+∼ δ P 0 − ∆E (A5)
δ P 0 − ∆E −4mχ4mχ1 In the zero-momentum-recoil approximation it is assumed thatK 
p 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, |E100 |]) 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 |E100 |. This captures the entire region of
the freeze-out we are interested in, since DM decouples from thebound states when temperature is already less than the groundstate energy, i.e., T ∼ |E100 |  α2 |E100 |. For electroweak couplings we have checked the zero-recoil approximation in the hightemperature regime mχ /T = 10 and noticed that the difference

6and the BSF cross section can be brought into the form:
Z

d4 P
1
BSF
/δ P 0 − ∆E
(σnlm
vrel ) = 2
4
N
(2π)
h
i

eq
0
× 1 + fφ (P ) sgn(P 0 )/δ(P 2 − m2φ ) |M|2k,nlm . (A6)Here, the theta function θ(P 0 ) was dropped, since thedelta function in the first line already ensures positivityof the emitted energy. Assuming a massive vector field,
one can decompose the amplitude as
XX
Mµk,nlm Mν?
|M|2k,nlm =µ (P, σ)?ν (P, σ)k,nlmσ

P µP ν
= − gµν −
m2φ

spins

!
X

Appendix B: Retarded self-energy expression

The retarded self-energy is defined in terms of greaterand lesser self-energies as
 −+

R
+−
Πµν
(x − y) = θ(x0 − y 0 ) Πµν(x − y)−Πµν(x − y) ,
(B1)
+−
Πµν
(x − y) = g 2 Tr[γµ B +− (x − y)γν B −+ (y − x)], (B2)−+Πµν(x − y) = g 2 Tr[γµ B −+ (x − y)γν B +− (y − x)], (B3)

where the two-point functions of the fermionic bathparticles areMµk,nlm Mν?
k,nlm , (A7)spins

where  is the polarization vector. Note that in ourmodel, 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 elementin Eq. (A6), one may notice the lesser two-point correlation function in thermal equilibrium Eq. (6) in the free
limit:hi−+,0Dµν(P 0 , p) = 1 + fφeq (P 0 ) ×"!
#
P µP ν
0 /22−gµν +
sgn(P )δ (P − mφ )(A8)
m2φ
After integrating Eq. (A6) over P 0 , and replacing the freelesser function by interacting one, Eq. (2) is recovered.
The justification of this replacement from first principles
can be found in ref. [56]. Here, we would like to argue
from the Kadanoff-Baym equation that this “upgrade” isexpected. As derived in ref. [67], the DM collision termin the forward direction can be determined through the
DM self-energy component Σ −+ . In the 1-PI formalism,−+. This is another inthe Σ −+ contains interacting Dµν
dication 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 schemeof 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). In non-abelian cases, caremust be taken with the definition of two-particles states,
see QCD works on pNRQFT. This is because a gluon for
example can be emitted also from the exchanged gluonof the DM pair.

is about only one per-cent compared to the exact result (butstill assuming dipole approximation). Also note that ref. [56] assumes the zero-momentum-recoil approximation to be valid fromthe beginning, since the derivation is based on pNRQFT.

−+
Bij(x − y) ≡ hbi (x)b̄j (y)i,

(B4)

+−
Bij
(x − y) ≡ −hb̄j (y)bi (x)i,(B5)

and their free Fourier transforms are in equilibrium
2/
B +− (K) = −K(2π)δ(K) −θ(−K 0 ) + fbeq (|K 0 |) , (B6)2/B −+ (K) = −K(2π)δ(K) −θ(+K 0 ) + fbeq (|K 0 |) . (B7)The Fourier transform of the retarded self-energy is thenRΠµν(P ) = g 2

Z/ 1 γν P/ 2 ]×dΠ1 dΠ2 Tr[γµ P

(B8)

3

i/δ (p + p1 − p2 )[f eq (|p1 |)−fbeq (|p2 |)]+
0
P + |p1 | − |p2 | + i b3i/δ (p + p2 − p1 )[f eq (|p2 |)−fbeq (|p1 |)]+
P 0 + |p2 | − |p1 | + i b
3i/δ (p − p1 − p2 )[1 − fbeq (|p1 |)−fbeq (|p2 |)]+
0P − |p1 | − |p2 | + i3i/δ (p + p1 + p2 )
eq
eq
[−1
+
f
(|p
|)+f
(|p
|)]
.
1
2
b
b
P 0 + |p1 | + |p2 | + iThe first and second term inside the brackets corresponds
to particle and anti-particle scattering, when taking theimaginary 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, sincein our case P 0 = ∆E is positive, leading to the fact thatthe cut of the propagator is zero. The following identity
was used for computing Eq. (11):[1 + fVeq (∆E)] [fbeq (|p1 |) − fbeq (|p1 | + ∆E)] =fbeq (|p1 |) [1 − fbeq (∆E + |p1 |)] ,(B9)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.

7
Appendix C: States convention and overlap integral

In terms of particle and anti-particle creation operators, the states in Eq. (4) follow the convention:√ Z d3 p0s0 †r0 †ψnlm (p0 )aQ/2+p|Bnlm i = 2M
0 bQ/2−p0 |0i, (C1)3(2π)
Z 3 00q
d p
r†|Sk i = 4m2χψk (p00 )as†
K/2+p00 bK/2−p00 |0i. (C2)(2π)3Q 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 thelong-range interactions. They have energy Eigenvalues
2
Ek = k2 /(2µ) = µvrel/2 and negative binding energy
Enlm , and are normalized to plane waves in the case of
scattering states, while for the bound states the spatialintegral over the absolute squared wave function is normalized to unity. The scattering state is a coherent sum
over all partial waves.
In Eq. (5) we defined h(ζ, ξ) in terms of the overlap
integral:pZh(ζ, ξ) k̂ i?
d3 r ψ100(r)ri ψk (r) ,≡κ1/2 µ ∆E(C3)[1] B. W. Lee and S. Weinberg, Phys.Rev.Lett. 39, 165(1977).
[2] J. R. Ellis, J. S. Hagelin, D. V. Nanopoulos, K. A. Olive,and M. Srednicki, Particle physics and cosmology: Dark
matter, Nucl. Phys. B238, 453 (1984), [,223(1983)].
[3] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, Eur.
Phys. J. C78, 203 (2018), arXiv:1703.07364 [hep-ph].
[4] P. A. R. Ade et al. (Planck), Astron. Astrophys. 594,
A13 (2016), arXiv:1502.01589 [astro-ph.CO].
[5] J. Hisano, S. Matsumoto, M. Nagai, O. Saito, and
M. Senami, Phys.Lett. B646, 34 (2007), arXiv:hepph/0610249 [hep-ph].
[6] A. Sommerfeld, Annalen der Physik 403, 257.
[7] A. D. Sakharov, Zh. Eksp. Teor. Fiz. 18, 631 (1948),[Usp. Fiz. Nauk161,no.5,29(1991)].
[8] K. Belotsky, D. Fargion, M. Khlopov, and R. V. Konoplich, Phys. Atom. Nucl. 71, 147 (2008), arXiv:hepph/0411093 [hep-ph].
[9] A. Freitas, Phys. Lett. B652, 280 (2007),
arXiv:0705.4027 [hep-ph].
[10] M. Cirelli, A. Strumia, and M. Tamburini, Nucl.Phys.
B787, 152 (2007), arXiv:0706.4071 [hep-ph].
[11] C. F. Berger, L. Covi, S. Kraml, and F. Palorini, JCAP
0810, 005 (2008), arXiv:0807.0211 [hep-ph].

where κ ≡ µα is the Bohr-momentum. We have extracted a specific combination of dimensional quantities
to make h(ζ, ξ) dimensionless. For convention, the analytic expression for h(ζ, ξ) in the Coulomb limit ξ → ∞is given by [57]:ζ6h(ζ, ∞) = 26 πS(ζ)
e−4ζacotζ ,(1 + ζ 2 )3

2πζ
S(ζ) =.
1 − e−2πζ

(C4)
(C5)

Appendix D: Sommerfeld enhanced annihilation anddecay rate

The total Sommerfeld enhanced s-wave annihilation
cross section, as well as the total s-wave decay rate ofbound states for the spin-singlet Sl=0 (χχ̄)[Bl=0 (χχ̄)] →V V and spin triplet Sl=0 (χχ̄)[Bl=0 (χχ̄)] → bb̄ reactions,
are given by:
2πα2|ψk,l=0 (r = 0)|2 ,
4µ2
2πα2
dec
Γn00
=
|ψn00 (r = 0)|2 .
4µ2(σ an vrel ) =

(D1)
(D2)

The cross section and decay width factorize into the shortdistance term and the reduced two-body wave function
2
evaluated at the origin. Note that ψn00
carries mass dimension three from normalization.

[12] A. Hryczuk, R. Iengo, and P. Ullio, JHEP 1103, 069(2011), arXiv:1010.2172 [hep-ph].
[13] A.
Hryczuk,
Phys.Lett.
B699,
271
(2011),
arXiv:1102.4295 [hep-ph].
[14] K. Harigaya, K. Kaneta, and S. Matsumoto, Phys. Rev.
D89, 115021 (2014), arXiv:1403.0715 [hep-ph].
[15] S. El Hedri, A. Kaminska, and M. de Vries, Eur. Phys.
J. C77, 622 (2017), arXiv:1612.02825 [hep-ph].
[16] M. Beneke, A. Bharucha, A. Hryczuk, S. Recksiegel, andP. Ruiz-Femenia, JHEP 01, 002 (2017), arXiv:1611.00804[hep-ph].
[17] M. Beneke, A. Bharucha, F. Dighera, C. Hellmann,
A. Hryczuk, S. Recksiegel, and P. Ruiz-Femenia, JHEP
03, 119 (2016), arXiv:1601.04718 [hep-ph].
[18] T. R. Slatyer, JCAP 1002, 028 (2010), arXiv:0910.5713[hep-ph].
[19] M. Beneke, C. Hellmann, and P. Ruiz-Femenia, JHEP
05, 115 (2015), arXiv:1411.6924 [hep-ph].
[20] K. Blum, R. Sato, and T. R. Slatyer, JCAP 1606, 021(2016), arXiv:1603.01383 [hep-ph].
[21] E. Braaten, E. Johnson, and H. Zhang, JHEP 11, 108(2017), arXiv:1706.02253 [hep-ph].
[22] E. Braaten, E. Johnson, and H. Zhang, JHEP 05, 062(2018), arXiv:1712.07142 [hep-ph].
[23] E. Braaten, E. Johnson, and H. Zhang, JHEP 02, 150(2018), arXiv:1708.07155 [hep-ph].

8[24] Y.-L. Tang and G.-L. Zhou, Phys. Rev. D99, 036016(2019), arXiv:1806.10124 [hep-ph].
[25] M. Pospelov and A. Ritz, Phys.Lett. B671, 391 (2009),
arXiv:0810.1502 [hep-ph].
[26] J. D. March-Russell and S. M. West, Phys.Lett. B676,
133 (2009), arXiv:0812.0559 [astro-ph].
[27] W. Shepherd, T. M. Tait, and G. Zaharijas, Phys.Rev.
D79, 055022 (2009), arXiv:0901.2125 [hep-ph].
[28] B. von Harling and K. Petraki, JCAP 12, 033 (2014),
arXiv:1407.7874 [hep-ph].
[29] J. Ellis, J. L. Evans, F. Luo, and K. A. Olive, JHEP 02,
071 (2016), arXiv:1510.03498 [hep-ph].
[30] S. P. Liew and F. Luo, JHEP 02, 091 (2017),
arXiv:1611.08133 [hep-ph].
[31] P. Asadi, M. Baumgart, P. J. Fitzpatrick, E. Krupczak,
and T. R. Slatyer, JCAP 1702, 005 (2017),
arXiv:1610.07617 [hep-ph].
[32] E. Johnson, E. Braaten,and H. Zhang (2016)arXiv:1611.06212 [hep-ph].
[33] S. Kim and M. Laine, JHEP 07, 143 (2016),
arXiv:1602.08105 [hep-ph].
[34] A. Mitridate, M. Redi, J. Smirnov, and A. Strumia,
JCAP 1705, 006 (2017), arXiv:1702.01141 [hep-ph].
[35] J. Harz and K. Petraki, JHEP 07, 096 (2018),
arXiv:1805.01200 [hep-ph].
[36] S. Biondini and M. Laine, JHEP 04, 072 (2018),
arXiv:1801.05821 [hep-ph].
[37] S. Biondini and S. Vogl, JHEP 02, 016 (2019),
arXiv:1811.02581 [hep-ph].
[38] H. Fukuda, F. Luo, and S. Shirai, JHEP 04, 107 (2019),
arXiv:1812.02066 [hep-ph].
[39] J. Harz and K. Petraki, Phys. Rev. D97, 075041 (2018),
arXiv:1711.03552 [hep-ph].
[40] S. Biondini, JHEP 06, 104 (2018), arXiv:1805.00353[hep-ph].
[41] J. Harz and K. Petraki, JHEP 04, 130 (2019),
arXiv:1901.10030 [hep-ph].
[42] A. Kamada, M. Kaplinghat, A. B. Pace, and H.-B. Yu,(2016), arXiv:1611.02716 [astro-ph.GA].
[43] M. R. Buckley and P. J. Fox, Phys.Rev. D81, 083522(2010), arXiv:0911.3898 [hep-ph].
[44] L. G. van den Aarssen, T. Bringmann, and C. Pfrommer,
Phys.Rev.Lett. 109, 231301 (2012), arXiv:1205.5809[astro-ph.CO].
[45] S. Tulin, H.-B. Yu, and K. M. Zurek, Phys.Rev. D87,
115007 (2013), arXiv:1302.3898 [hep-ph].
[46] T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, andP. Walia, (2016), arXiv:1612.00845 [hep-ph].
[47] H. An, M. B. Wise, and Y. Zhang, Phys. Rev. D93,
115020 (2016), arXiv:1604.01776 [hep-ph].
[48] I. Baldes and K. Petraki, JCAP 1709, 028 (2017),
arXiv:1703.00478 [hep-ph].
[49] I. Baldes, M. Cirelli, P. Panci, K. Petraki, F. Sala, and
M. Taoso, SciPost Phys. 4, 041 (2018), arXiv:1712.07489[hep-ph].
[50] A. Kamada, K. Kaneta, K. Yanagi, and H.-B. Yu, JHEP
06, 117 (2018), arXiv:1805.00651 [hep-ph].
[51] A. Kamada, M. Yamada, and T. T. Yanagida, JHEP
03, 021 (2019), arXiv:1811.02567 [hep-ph].
[52] A. Kamada, M. Yamada, and T. T. Yanagida, (2019),
arXiv:1908.00207 [hep-ph].
[53] P. Ko, T. Matsui,
and Y.-L. Tang,(2019),
arXiv:1910.04311 [hep-ph].

[54] N. Brambilla, M. A. Escobedo, J. Ghiglieri, and A. Vairo,
JHEP 05, 130 (2013), arXiv:1303.6097 [hep-ph].
[55] J. Hong and S. H. Lee, Phys. Rev. C99, 034905 (2019),
arXiv:1811.07607 [nucl-th].
[56] T. Binder et al., in preparation.
[57] K. Petraki, M. Postma, and J. de Vries, JHEP 04, 077(2017), arXiv:1611.01394 [hep-ph].
[58] T. Binder, L. Covi, and K. Mukaida, Phys. Rev. D98,
115023 (2018), arXiv:1808.06472 [hep-ph].
[59] S. Kim and M. Laine, JCAP 1701, 013 (2017),
arXiv:1609.00474 [hep-ph].
[60] S. Biondini and M. Laine, JHEP 08, 047 (2017),
arXiv:1706.01894 [hep-ph].
[61] S. Kim and M. Laine, (2019), arXiv:1904.07882 [hep-ph].
[62] S. Biondini, S. Kim,
and M. Laine,(2019),
arXiv:1908.07541 [hep-ph].
[63] E. Braaten and R. D. Pisarski, Nucl. Phys. B337, 569(1990).
[64] M. Laine, O. Philipsen, P. Romatschke, and M. Tassler,
JHEP 03, 054 (2007), arXiv:hep-ph/0611300 [hep-ph].
[65] N. Brambilla, J. Ghiglieri, A. Vairo, and P. Petreczky,
Phys. Rev. D78, 014017 (2008), arXiv:0804.0993 [hepph].
[66] A. M. Sirunyan et al. (CMS), Phys. Lett. B790, 270(2019), arXiv:1805.09215 [hep-ex].
[67] M. Beneke, F. Dighera, and A. Hryczuk, JHEP 10, 045(2014), [Erratum: JHEP07,106(2016)], arXiv:1409.3049[hep-ph].