Thermal dark matter co-annihilating with a strongly interacting scalar

Recently many investigations have considered Majorana dark matter co-annihilating with bound states formed by a strongly interacting scalar field. However only the gluon radiation contribution to bound state formation and dissociation, which at high temperatures is subleading to soft 2->2 scatterings, has been included. Making use of a non-relativistic effective theory framework and solving a plasma-modified Schrodinger equation, we address the effect of soft 2->2 scatterings as well as the thermal dissociation of bound states. We argue that the mass splitting between the Majorana and scalar field has in general both a lower and an upper bound, and that the dark matter mass scale can be pushed at least up to 5...6 TeV.


Introduction
Negative results from direct and indirect detection experiments and collider searches pose a challenge for many minimal dark matter models. This has led to the construction of less minimal models. In the latter, the cross sections probed experimentally are not directly related to the cross section affecting freeze-out dynamics in the early universe. Therefore experimental bounds might be respected while at the same time maintaining the correct cosmological abundance.
If the dark matter particles are massive and interact strongly enough with the Standard Model to have been in equilibrium with it at some time in the early universe, the basic feature that is needed for the above task is a strongly temperature-dependent annihilation cross section. At low temperatures, the cross section needs to be very small, to satisfy the non-observation bounds from indirect detection. In the early universe, the cross section needs to be large enough to keep dark matter in chemical equilibrium for a long while, reducing its number density and thereby evading overclosure of the universe.
An example of a possible scenario along these lines is to postulate a model in which the dark sector consists of two particle species. The lighter one is the true dark matter, long-lived and interacting very weakly. In contrast, the heavier one could interact strongly and act as an efficient dilution channel for the overall abundance at high temperatures (cf., e.g., ref. [1]).
If the heavy species interacts strongly, as in QCD, this scenario could lead to rather rich phenomenology. Strongly interacting particles form generally bound states. Because of the associated binding energy, their thermal abundance is larger than that for the same particles in a scattering state. Bound states may annihilate efficiently because the two particles are close to each other. Though often alluded to previously, a more intensive study of bound-state effects on the freeze-out dynamics has only started a few years ago (cf., e.g., refs. [2,3]).
Recently, we have participated in developing a non-perturbative formalism for addressing the thermal annihilation of non-relativistic particles [4,5]. The formalism was already applied to a first full model, which did not include strongly interacting particles but nevertheless displayed weakly bound states [6]. The purpose of the current paper is to apply the formalism to a strongly interacting model that has been much discussed in recent literature.
Our plan is as follows. Having introduced the model in sec. 2, we review some salient features concerning its thermal behaviour in sec. 3. The main technical ingredients of our analysis are specified in sec. 4: the operators responsible for the hard annihilation event; the spectral functions describing the soft initial-state effects that influence the annihilation; as well as generalized "Sommerfeld factors" which capture the effect both of bound and scattering states on the thermal annihilation cross sections. The cosmological evolution equations are integrated in sec. 5, whereas conclusions and an outlook are offered in sec. 6.

Model
The model considered consists of the Standard Model extended by a gauge singlet Majorana fermion (χ) as well as a scalar field (η) which is singlet in SU L (2) but carries non-trivial QCD and hypercharge quantum numbers. 1 The Majorana fermion is chosen as the dark matter particle, given that its low-energy scattering cross section is naturally suppressed, being pwave at tree level [8]. In the MSSM language, the Majorana fermion could be a bino-like neutralino and the scalar a right-handed stop or sbottom. However, for generality we do not fix couplings to their MSSM values. The hypercharge coupling of the scalar is generally omitted, as its effects are subleading compared with QCD effects.
The Lagrangian for this extension of the Standard Model can be expressed as (2.1) The notation λ 1 is reserved for the self-coupling of the Higgs doublet (H). The chiral projectors a L = (½ − γ 5 )/2, a R = (½ + γ 5 )/2 imply that χ only interacts with SU L (2) singlet projections of quarks. We assume that the Yukawa coupling y couples dominantly to one quark flavour only. The Yukawa coupling determining the mass of that flavour is denoted by h, and the strong gauge coupling by g s . The free parameters of the model are then the two mass scales 2 M and ∆M ≡ M η − M as well as the two "portal" couplings λ 3 and |y| 2 that are assumed to be small at the MS scaleμ ∼ 2M .
In the MSSM context, the importance of co-annihilations in such a model was stressed long ago [1]. Sommerfeld enhancements from QCD interactions were included in refs. [9,10], however without the consideration of bound states. Similar theoretical ingredients were applied to the generalized model in ref. [11]. Direct, indirect and collider constraints on the generalized model were reviewed in ref. [7]. More recently, bound-state effects have been approximately included in this model [12][13][14], as a single additional degree of freedom in a set of Boltzmann equations, a treatment which we aim to improve upon in the following.

Parametric forms of thermal masses and interaction rates
The coloured scalars are responsible for most of the annihilations during thermal freeze-out. We start by reviewing the thermal mass corrections and interaction rates that they experience.
The important point is that, because of Bose enhancement, the gluonic contributions are infrared (IR) sensitive, and need to be properly resummed for a correct result.
As a first step, consider a naive (i.e. unresummed) computation of the self-energy of the coloured scalar. Evaluating the (retarded) self-energy at the on-shell point yields (the line " " stands for η and the wiggly line for a gluon) The real part is analogous to that for a heavy fermion [15]. The imaginary part vanishes because there is no phase space for the 1 ↔ 2 process.
However, at high temperatures these naive results are misleading. Perhaps the simplest way to see this is to replace the scalar in the loop by a particle with a different mass, M η + ∆M , and consider the case ∆M ≪ πT ≪ M η . Then it can be verified that Re Π R /M η is modified by a correction of order ∼ g 2 s C F ∆M , and Im Π R /M η by a correction of order In other words, the result in eq. (3.2) seems to change qualitatively because Bose enhancement of the soft contribution compensates against the phase-space suppression.
The correct treatment of the Bose-enhanced IR contribution requires resummation. The heavy scalars are almost static, and interact mostly with colour-electric fields (A a 0 ). In a plasma, colour-electric fields get Debye screened. We denote the Debye mass by m D . Parametrically, m D ∼ g s T , where g s ≡ √ 4πα s . The proper inclusion of Debye screening in a gauge theory requires Hard Thermal Loop (HTL) resummation [16][17][18][19]. Recomputing the 1-loop self-energy with HTL propagators, and setting ∆M → 0 since IR sensitivity has now been regulated, we get (here p ≡ |p| and a blob stands for a HTL-resummed propagator) The new contribution in eq. (3.3), originating from the Debye-screened Coulomb self-energy, is known as the Salpeter correction (cf. ref. [20] for a review). It dominates over the other mass correction if T < ∼ g s M η , which is generally the case. The imaginary part in eq. (3.4), i.e. the interaction rate, reflects fast colour and phase-changing 2 → 2 scatterings off light medium particles. It was first derived for the case of a heavy quark [16]. We finally replace the coloured scalar by a pair of heavy scalars, separated by a distance r. The HTL-resummed computation of the thermal mass correction ("static potential") and interaction rate as a function of r was carried out in refs. [21][22][23]. At leading non-trivial order the result can be expressed as As a crosscheck, for r → ∞ twice the results of eqs.
This can be compared with the 1 ↔ 2 gluon radiation contribution, ∼ g 2 s C F (∆E) 3 r 2 n B (∆E) [23], where ∆E is the energy difference between the singlet and octet potentials. At high temperatures, when m D , πT ≫ ∆E, the 2 → 2 contribution dominates over the 1 ↔ 2 one.
In order to determine the spectral function of the scalar pair, characterizing the states that appear in the scalar-antiscalar sector of the Fock space, V (r) and Γ(r) can be inserted into a time-dependent Schrödinger equation satisfied by the appropriate Green's function [24]. More details are given in sec. 4. We have checked numerically that, in accordance with theoretical expectations [25], the states originating from this solution respect the qualitative pattern seen above for Γ(r), namely that at high temperatures the width from 2 → 2 reactions dominates over the gluo-dissociation contribution.
We close this section by considering another essential ingredient of the framework, namely the rate at which Majorana dark matter particles convert into the coloured scalars. Once more, this rate is dominated by 2 → 2 scatterings, and obtaining the correct result requires HTL resummation. Setting for simplicity the external momentum to zero, we find (the thick line is the Majorana fermion and the arrowed line the quark flavour with which it interacts, treated for simplicity as massless in vacuum which is a good approximation if m vac < ∼ πT ) where the last line applies under the assumption ∆M ≪ m q , πT ≪ √ T M . The thermal quark mass m q , originating from the phase space integral of the light plasma particles off which the 2 → 2 scattering takes place, is The rate in eqs. (3.8) and (3.9) is faster than the Hubble rate in a broad temperature range, e.g. down to M/T > ∼ 3000 for y = 0.3 and ∆M/M < ∼ 0.01. It does fall out of equilibrium when T ≪ ∆M , however transitions to virtual bound-state constituents may continue and form presumably the relevant concern. Non-equilibrium effects have been discussed in ref. [26].

Quantitative framework for estimating the annihilation rate
We now present a framework for computing (co-)annihilation rates in the model of sec. 2.

Non-relativistic fields
The basic premise of our framework is to make use of the non-relativistic approximation, assuming that πT , m top , ∆M ≪ M , where M is the dark matter mass and ∆M = M η − M is the mass splitting within the dark sector. This simplification opens up the avenue to a non-relativistic effective field theory investigation of soft initial-state effects. In the non-relativistic limit, the interaction picture field operator of the coloured scalar is expressed as The non-relativistic fields φ and ϕ † transform in the fundamental representation of SU(N c ), with colour indices denoted by α, β, γ, δ, ... . The Majorana spinor χ is simplest to handle by choosing the standard representation for the Dirac matrices, i.e. γ 0 = diag(½, −½). Then where the Grassmannian spinor ψ has two spin components, labelled by p, q, r, s, ... . Only the left-chiral projection of χ participates in interactions according to eq. (2.1).
In the following, we generally set M η → M whenever possible. The influence of ∆M = 0 (and its thermal modification) is discussed in sec. 4.3.

Imaginary parts of 4-particle operators
The first step is to determine annihilation cross sections for all possible processes with dark matter initial states. The leading order Feynman diagrams are shown in fig. 1. According to the optical theorem, the amplitudes squared |M| 2 can be expressed as an imaginary (or "absorptive") contribution to an effective Lagrangian [27]. An important simplification in the Majorana case follows from the identity satisfied by Pauli matrices, σ k pq σ k rs = 2δ ps δ qr − δ pq δ rs . Therefore a possible spin-dependent operator can be reduced to a spin-independent one: ψ † p ψ † r ψ s ψ q σ k pq σ k rs = −3ψ † p ψ † q ψ q ψ p . At leading order in an expansion in 1/M 2 , the absorptive operators read Here T a are Hermitean generators of SU(N c ). In the partial wave language, the operators in eq. (4.4) A non-zero value of c 1 may be generated at higher orders. To minimize the magnitude of higher-order effects, the couplings should be evaluated at the MS renormalization scalē µ ∼ 2M . We note that c 5 gets contributions from the "Majorana-like" processes M i and M j in fig. 1, but not from the "Dirac-like" amplitude M h .

Number density, effective cross section, evolution equation
Within Boltzmann equations the overall dark matter abundance evolves as [28][29][30] n = − σ eff v n 2 − n 2 eq , (4.5) whereṅ is the covariant time derivative in an expanding background. To go beyond the quasiparticle approximation underlying the Boltzmann approach, the effective cross section can be re-interpreted as a chemical equilibration rate, Γ chem , and then defined on the non-perturbative level within linear response theory [31]. Furthermore, within the nonrelativistic effective theory, Γ chem can be related to the thermal expectation value of L abs from eq. (4.3) [4]. These relations can be expressed as Im L abs . (4.6) In our model the number density amounts to The mass difference ∆M T gets a vacuum contribution, ∆M = M η − M , and a thermal correction from eq. (3.3) as well as from a similar tadpole involving λ 3 , Note that the negative Salpeter correction may cancel against the positive terms. At leading order the Debye mass parameter amounts to 9) where N f is the number of quark flavours (cf. ref. [32] for higher orders). The effective values of g s and N f are changed with the temperature, as reviewed in appendix A. For future reference we define a "tree-level" effective cross section, σ eff v (0) , by evaluating the thermal expectation value Im L abs at leading order and then making use of eqs. (4.6) and (4.7). Wick contracting the indices in eq. (4.3) leads to (4.10)

Plasma-modified Schrödinger equation and generalized Sommerfeld factors
Going beyond leading order, we evaluate Im L abs as elaborated upon in ref. [5], expressing it as a Laplace transform of a spectral function characterizing the dynamics of the dark matter particles before their annihilation. Denoting by E ′ the energy of the relative motion and by k the momentum of the center-of-mass motion, this implies where α 2 M ≪ Λ ≪ M restricts the average to the non-relativistic regime. 3 The spectral functions are obtained as imaginary parts of Green's functions, 4 Here V i contains a negative imaginary part, and N i is a normalization factor giving the number of contractions related to the operator that c i multiplies in eq. (4.3): If the potentials V i (r) were r-independent and with an infinitesimal imaginary part, i.e. V i (r) = Re V i (∞)−i0 + , they would only induce mass shifts. In this case the spectral functions can be determined analytically, This form can be used for defining generalized Sommerfeld factors: Then eq. (4.11) combined with eqs. (4.6) and (4.7) leads to a generalization of eq. (4.10), If a potential V i (r) leads to a bound state, whose width is much smaller than the binding energy, the corresponding generalized Sommerfeld factor can be computed in analytic form. In this case eq. (4.12) can be solved in a spectral representation, resulting in 3 Some elaboration about the need to introduce such a cutoff can be found in ref. [5]. In practice, we choose Λ ≃ 2α 2 M , and have verified that making it e.g. 2-3 times larger plays no role on our numerical resolution. where ψ j are the bound state wave functions. Inserting into eq. (4.16), the contribution of the jth bound state toS i reads This becomes (exponentially) large when T ≪ α 2 s M , however chemical equilibrium is lost in the dark sector at low T , which imposes an effective cutoff on the growth (cf. secs. 5 and 6).

Thermal potentials
In order to write down the potentials V i (r) appearing in eq. (4.12), let us define The integrand in eq. (4.20) corresponds to the static limit of the time-ordered HTL-resummed temporal gluon propagator. Then we find The structure V 3 (r) equals the combination V (r) − iΓ(r) shown in eqs. (3.5)-(3.7), whereas C F Re[v(0)] yields the Salpeter part of ∆M T in eq. (4.8). The potential V 3 (r) corresponds to a singlet potential, V 4 (r) to an octet potential, and V 5 (r) to a particle-particle potential, relevant because of the presence of a particle-particle annihilation channel generated by Majorana exchange (cf. the discussion around the end of sec. 4.2). We note in passing that at T < 160 GeV, when the Higgs mechanism is operative, additional potentials can be generated, particularly through the Higgs portal coupling λ 3 in eq. (2.1) (cf. e.g. ref. [33]). However the coefficients of these potentials are suppressed by ∼ λ 2 3 v 2 /M 2 , where v is the Higgs expectation value. Given that we consider M ≥ 2 TeV, we expect their contributions to be negligible compared with QCD effects and have not included them. We also note that an r-dependence can be generated for V 2 (r) through quark exchange, however this is suppressed by ∼ |y| 2 σ · ∇/M . For a practical use of eq. (4.21), numerical values are needed for the parameters g 2 s and m 2 D . We relegate a discussion of this point into appendix A. Let us however note that we Here ω ≡ 2M + E ′ . At low temperatures a dense spectrum of bound states can be observed, which gradually "melts away" as the temperature increases. Right: The generalized Sommerfeld factors, eq. (4.16), corresponding to the annihilation of the coloured scalars via different channels.
restrict to temperatures T > ∼ 1 GeV, so that the real part of the potential contains no trace of a string tension [34]. Furthermore, in accordance with the low-temperature gluon-radiation contribution specified below eq. (3.7) and with general arguments presented in ref. [6], the imaginary part of the potential is multiplied by the Boltzmann factor e −|E ′ |/T for E ′ < 0.

Numerical evaluations
Having determined the spectral functions from eqs. (4.12) and (4.13) and the generalized Sommerfeld factors from eq. (4.16) or (4.19), the effective cross section is obtained from eq. (4.17). Subsequently eq. (4.5) can be integrated for the dark matter abundance. As usual we define a yield parameter as Y ≡ n/s, where s is the entropy density, and change variables from time to z ≡ M/T , whereby eq. (4.5) becomes Here m Pl is the Planck mass, e is the energy density, and c is the heat capacity, for which we use values from ref. [36] (cf. also ref. [37]). The final value Y (z final ) yields the energy fraction Ω dm h 2 = Y (z final ) M /[3.645 × 10 −9 GeV]. We integrate eq. (5.1) up to z final = 10 3 . At around these temperatures, depending on the value of ∆M/M , the processes of interest have either ceased to be active, or are falling out of chemical equilibrium, because their rates are suppressed by e −∆M/T ≪ 1. Therefore they cannot be reliably addressed within the current framework.
In fig. 2(left) we show the spectral function ρ 3 corresponding to the attractive channel, displaying a dense spectrum of bound states at low temperatures. The corresponding generalized Sommerfeld factor, obtained from eq. (4.16), is shown in fig. 2(right). An exponential increase is observed at low temperatures, as indicated by eq. (4.19). The repulsive channels also show a modest increase at very low temperatures, due to the fact that the spectral function extends below the threshold at finite temperature [6]. Examples of results obtained by integrating eq. (5.1) are shown in fig. 3. In particular, it can be observed how a very efficient annihilation sets in at low temperatures, if ∆M is small so that bound states of coloured scalars are lighter than scattering states of Majorana fermions. Finally, fig. 4 shows slices of the parameter space leading to the correct dark matter abundance.
In the plots the Yukawa couplings have been set to the stop-like values y = 0.3, h = 1.0. However these couplings only have a modest effect if chosen otherwise, because they do not affect the coefficient c 3 appearing the attractive channel, cf. eq. (4.4). As an example, setting h = 0.0 increases the abundance typically by ∼ 5%, cf. fig. 3. The most important role is played by the coupling λ 3 . For c 3 this coupling has been evaluated at the scaleμ = 2M , whereas for collider phenomenology its value at a scaleμ ∼ m H would be more relevant. The latter can be obtained from eq. (A.7), and is some tens of percent smaller than λ 3 (2M ). We stress that, as shown by eq. (A.7), Yukawa couplings always generate a non-zero value for λ 3 through renormalization group running.

Conclusions
We have investigated a simple extension of the Standard Model, cf. sec. 2, which has become popular as a prototypical fix to the increasingly stringent empirical constraints placed on "WIMP"-like frameworks. In this model dark matter consists of a Majorana fermion, which only has a p-wave annihilation cross section at tree level, helping to respect experimental non-observation constraints from indirect detection. The Majorana fermion has a Yukawa interaction with a QCD-charged scalar field (such as a right-handed stop or sbottom in the MSSM) and a Standard Model quark. For large masses and small mass splittings between the Majorana fermion and the scalar field, the best sensitivity for discovering the Majorana fermion appears to be direct detection by XENON1T [7], enhanced by resonant scattering off quarks through scalar exchange, even if interactions with top or bottom quarks are much less constrained than those with up or down quarks. Despite its simplicity, the model displays rich physics in the early universe. We have extended previous investigations [7,[9][10][11][12][13][14] by incorporating the full spectrum of thermally broadened bound states as well as the effect of soft 2 ↔ 2 scatterings. In general such scatterings dominate interaction rates at small mass splittings, because they are not phasespace suppressed in the same way as 1 ↔ 2 scatterings are, cf. sec. 3.
The reason that the model leads to a viable cosmology is that at high temperatures dark matter annihilates efficiently through the scalar channel, guaranteeing that its overall abundance remains low. The fast annihilations proceed particularly through bound states formed by the scalars, cf. fig. 2. As shown in fig. 4, the model can be phenomenologically viable for masses up to M ∼ 5...6 TeV, provided that the mass splitting is small, ∆M/M < 5 × 10 −3 , and that the "Higgs portal" coupling λ 3 between the coloured scalar and the Higgs doublet is substantial. We recall that in supersymmetric theories, λ 3 is proportional to the quark Yukawa coupling squared, λ 3 ∼ |h| 2 , and therefore indeed large if we identify the coloured scalar as a right-handed stop. Actually, similar arguments but a somewhat more complicated analysis are expected to apply to a left-handed stop as well (cf. e.g. ref. [38]).
We believe that the mass splitting should not be too small, however. The non-relativistic binding energy of the lightest bound state, E ′ 1 , is negative. If it overcompensates for the mass difference, so that 2∆M + E ′ 1 < 0, the lightest two-particle states in the dark sector are the bound states formed by the coloured scalars. However these states are short-lived. Therefore it seems possible that (almost) all dark matter converts into the scalars and gets subsequently annihilated, so that the model may not be viable as an explanation for the observed dark matter abundance. This domain has been excluded through the grey bands in fig. 4. If we close eyes to this concern and assume that chemical equilibrium is maintained, then the value of M could be substantially larger than in fig. 4, for instance M ∼ 8 TeV as shown in fig. 3, and even more if we integrate down to lower temperatures.
We end by remarking that the model contains two portal couplings, λ 3 and y. The roles that these play are rather different. The value of λ 3 at the scaleμ = 2M influences the coefficient c 3 which mediates the most efficient annihilations, cf. eq. (4.4). In contrast y affects the rate of transitions between the Majorana fermions and coloured scalars, cf. eq. (3.9), as well as the running of λ 3 , cf. eq. (A.7). As long as y is not miniscule, so that the rate in eq. (3.9) remains in equilibrium, it has in practice little influence on our main results in fig. 4.
The only coupling that we need at a scaleμ ≪ M is the strong coupling. Since it has a large influence, we evaluate it at 2-loop level forμ ≤ M (nowadays running is known up to 5-loop level [39][40][41]). Denoting by N f the number of flavours and setting N c = 3 for brevity, the 2-loop running is given by The value of N f = 3, ..., 6 is changed when a quark mass threshold is crossed atμ = m i , where continuity is imposed. The initial value is α s (m Z ) ≃ 0.118. Forμ > M , the contribution of the coloured scalar is added and we switch over to 1-loop running, i.e. eq. (A.11). When we evaluate the static potential, a wide range of distance scales appears. At short distances, inspired by refs. [42,43], we evaluate the 2-loop coupling at the scaleμ = e −γ E /r. Since parametrically only the scales αM ≪ M play a role in the Schrödinger equation, the running does not include the coloured scalar in this domain.
At large distances, we employ effective thermal couplings. In the absence of NLO computations for thermal quarkonium observables, we adopt effective couplings from another context, that of dimensionally reduced field theories [44,45]. There the Debye mass parameter and an "electrostatic" coupling are expressed as [46]  For general masses, only α MS E4 and α MS E7 are available at present: Here the functions read (n F (x) ≡ 1/(e x + 1); chemical potentials have been set to zero) Given that α MS E6 is not currently known for general masses, we estimate inserting here PDG values for the quark masses [47]. The scale parameter is set toμ = 2πT .