1 Introduction

April 2018

{centering}

Thermal dark matter co-annihilating
with a strongly interacting scalar

S. Biondini and M. Laine

AEC, Institute for Theoretical Physics, University of Bern,
Sidlerstrasse 5, CH-3012 Bern, Switzerland

Abstract

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 scatterings, has been included. Making use of a non-relativistic effective theory framework and solving a plasma-modified Schrödinger equation, we address the effect of soft 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  TeV.

## 1 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.

## 2 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(2) but carries non-trivial QCD and hypercharge quantum numbers.111An SU(2) doublet would lead to similar results but a somewhat more complicated analysis [7]. The Majorana fermion is chosen as the dark matter particle, given that its low-energy scattering cross section is naturally suppressed, being -wave 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

 L = (2.1) − λ3η†ηH†H−yη†¯χa\tiny\rm{R}q−y∗¯qa\tiny\rm{L}χη.

The notation is reserved for the self-coupling of the Higgs doublet (). The chiral projectors , imply that only interacts with SU(2) singlet projections of quarks. We assume that the Yukawa coupling couples dominantly to one quark flavour only. The Yukawa coupling determining the mass of that flavour is denoted by , and the strong gauge coupling by . The free parameters of the model are then the two mass scales222More precisely, and refer to the renormalized parts of the masses appearing in the non-relativistic effective theory for and . The non-perturbative QCD contribution to is of the order which is smaller than the effects that we discuss below. and as well as the two “portal” couplings and that are assumed to be small at the scale .

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.

## 3 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)

 = g2sC\tiny\rm{F}T212Mη, (3.1) = 0, (3.2)

where . The real part is analogous to that for a heavy fermion [15]. The imaginary part vanishes because there is no phase space for the 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, , and consider the case . Then it can be verified that is modified by a correction of order , and by a correction of order , where is the Bose distribution. 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 (). In a plasma, colour-electric fields get Debye screened. We denote the Debye mass by . Parametrically, , where . 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 since IR sensitivity has now been regulated, we get (here and a blob stands for a HTL-resummed propagator)

 = g2sC\tiny\rm{F}T212Mη+g2sC\tiny\rm{F}2∫p1p2+m2\tiny\rm% {D} (3.3) = g2sC\tiny\rm{F}T212Mη−g2sC\tiny\rm{F}m\tiny\rm{D}8π, = −g2sC\tiny\rm{F}2∫pπTm2\tiny\rm{D}p(p2+m2\tiny\rm{D})2 (3.4) = −g2sC\tiny\rm{F}T8π.

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 , which is generally the case. The imaginary part in eq. (3.4), i.e. the interaction rate, reflects fast colour and phase-changing 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 . The HTL-resummed computation of the thermal mass correction (“static potential”) and interaction rate as a function of was carried out in refs. [21, 22, 23]. At leading non-trivial order the result can be expressed as

 G(r,t) t→+∞∼ G(r,0)exp{−i[V(r)−iΓ(r)]t}, (3.5) \;\parbox[c]{120.0pt}{\begin{picture}(120.0,25.% 0)(0.0,0.0)\SetWidth{1.0}\SetScale{1.0} \SetScale{0.7} \SetWidth{3.0} \DashLine(0,14)(120,14){3}\DashLine(0,-8)(120,-8){3}\SetWidth{1.0} \PhotonArc(60,15)(15,0,180){1.3} {6.283 15 mul 360 div 0 180 sub 0 180 sub mul sqrt mul Ldensity mul}\Photon(10% 0,-7)(100,13){1.3} {100 100 sub 100 100 sub mul -7 13 sub -7 13 sub mul add sqrt Ldensity mul}% \COval(60,30)(4,4)(0){Black}{Black}\COval(100,3)(3,3)(0){Black}{Black}% \LongArrow(7,-6)(7,12)\LongArrow(7,12)(7,-6)\Text(0,1.5)[c]{r} \end{picture}% }\;V(r) = −g2sC\tiny\rm{F}4π[m% \tiny\rm{D}+exp(−m\tiny\rm{D}r)r], (3.6) Γ(r) = g2sC\tiny\rm{F}T2π∫∞0dzz(z2+1)2[1−sin(zm\tiny\rm{D}r)zm\tiny\rm{D}r]. (3.7)

As a crosscheck, for twice the results of eqs. (3.3) and (3.4) are reproduced.

The interaction rate in eq. (3.7) can again be traced back to scatterings. At short distances, up to logarithms, . This can be compared with the gluon radiation contribution,  [23], where is the energy difference between the singlet and octet potentials. At high temperatures, when , the contribution dominates over the 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, and 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 , namely that at high temperatures the width from 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 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 )

 = −|y|2Nc8M∫pπm2qn\tiny\rm{F{}}(ΔM+p22M)p(p2+m2q) (3.8) ≈ −|y|2Ncm2q64πMln(1.76388MTm2q), (3.9)

where the last line applies under the assumption . The thermal quark mass , originating from the phase space integral of the light plasma particles off which the scattering takes place, is

 m2q=2g2sC\tiny\rm{F}∫qn\tiny% \rm{B{}}(q)+n\tiny\rm{F{}}(q)q=g2sT2C% \tiny\rm{F}4. (3.10)

The rate in eqs. (3.8) and (3.9) is faster than the Hubble rate in a broad temperature range, e.g. down to for and . It does fall out of equilibrium when , however transitions to virtual bound-state constituents may continue and form presumably the relevant concern. Non-equilibrium effects have been discussed in ref. [26].

## 4 Quantitative framework for estimating the annihilation rate

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

### 4.1 Non-relativistic fields

The basic premise of our framework is to make use of the non-relativistic approximation, assuming that , , , where is the dark matter mass and 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

 (4.1)

The non-relativistic fields and transform in the fundamental representation of SU(), with colour indices denoted by . The Majorana spinor is simplest to handle by choosing the standard representation for the Dirac matrices, i.e. . Then

 χ=(ψe−iMt−iσ2ψ∗eiMt),¯χ=(ψ†eiMt,−ψTiσ2e−iMt), (4.2)

where the Grassmannian spinor has two spin components, labelled by . Only the left-chiral projection of participates in interactions according to eq. (2.1).

In the following, we generally set whenever possible. The influence of (and its thermal modification) is discussed in sec. 4.3.

### 4.2 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 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, . Therefore a possible spin-dependent operator can be reduced to a spin-independent one: .

At leading order in an expansion in , the absorptive operators read

 L\scriptsize abs = i{c1ψ†pψ†qψqψp+c2(ψ†pϕ†αψpϕα+ψ†pφ†αψpφα) (4.3) + c3ϕ†αφ†αφβϕβ+c4ϕ†αφ†βφγϕδTaαβTaγδ+c5(ϕ†αϕ†βϕβϕα+φ†αφ†βφβφα)}.

Here are Hermitean generators of SU(). In the partial wave language, the operators in eq. (4.3) correspond to -wave annihilations, whereas -wave annihilations would lead to operators of . At leading order in couplings, the coefficients read

 c1=0,c2=|y|2(|h|2+g2sC\tiny\rm{F})128πM2, c3=132πM2(λ23+g4sC\tiny\rm{F}Nc),c4=g4s(N2c−4)64πM2Nc,c5=|y|4128πM2. (4.4)

A non-zero value of may be generated at higher orders. To minimize the magnitude of higher-order effects, the couplings should be evaluated at the renormalization scale . We note that gets contributions from the “Majorana-like” processes and in fig. 1, but not from the “Dirac-like” amplitude .

### 4.3 Number density, effective cross section, evolution equation

Within Boltzmann equations the overall dark matter abundance evolves as [28, 29, 30]

 ˙n=−⟨σ\scriptsize effv⟩(n2−n2\scriptsize 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, , and then defined on the non-perturbative level within linear response theory [31]. Furthermore, within the non-relativistic effective theory, can be related to the thermal expectation value of from eq. (4.3[4]. These relations can be expressed as

 ⟨σ\scriptsize effv⟩=Γ% \scriptsize chem2n\scriptsize eq=4n2% \scriptsize eq⟨ImL\scriptsize abs% ⟩. (4.6)

In our model the number density amounts to

 n\scriptsize eq=∫pe−Ep/T{2+2Nce−ΔMT/T},Ep≡M+p22M. (4.7)

The mass difference gets a vacuum contribution, , and a thermal correction from eq. (3.3) as well as from a similar tadpole involving ,

 ΔMT≡ΔM+(g2sC\tiny\rm{F}+λ3)T212M−g2sC\tiny\rm{F}m\tiny\rm{D}8π. (4.8)

Note that the negative Salpeter correction may cancel against the positive terms. At leading order the Debye mass parameter amounts to

 m\tiny\rm{D}=gsT√Nc3+Nf6, (4.9)

where is the number of quark flavours (cf. ref. [32] for higher orders). The effective values of and are changed with the temperature, as reviewed in appendix A.

For future reference we define a “tree-level” effective cross section, , by evaluating the thermal expectation value at leading order and then making use of eqs. (4.6) and (4.7). Wick contracting the indices in eq. (4.3) leads to

 ⟨σ\scriptsize effv⟩(0)=2c1+4c2Nce−ΔMT/T+[c3+c4C\tiny\rm{F}+2c5(Nc+1)]Nce−2ΔMT/T(1+Nce−ΔMT/T)2. (4.10)

### 4.4 Plasma-modified Schrödinger equation and generalized Sommerfeld factors

Going beyond leading order, we evaluate 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  the energy of the relative motion and by  the momentum of the center-of-mass motion, this implies

 ⟨ImL\scriptsize abs% ⟩ ≈ ∫ke−2MT−k24MT∫∞−ΛdE′πe−E′/T∑iciρi(E′) (4.11) = (MTπ)3/2e−2M/T∫∞−ΛdE′πe−E′/T∑iciρi(E′),

where restricts the average to the non-relativistic regime.333Some elaboration about the need to introduce such a cutoff can be found in ref. [5]. In practice, we choose , and have verified that making it e.g. 2-3 times larger plays no role on our numerical resolution. The spectral functions are obtained as imaginary parts of Green’s functions,444At higher orders in the non-relativistic expansion, kinetic terms and potentials suppressed by powers of could be added. In addition, operators suppressed by should be added in eq. (4.3).

 [−∇2rM+Vi(r)−E′]Gi(E′;r,r′) = Niδ(3)(r−r′), (4.12) limr,r′→0ImGi(E′;r,r′) = ρi(E′). (4.13)

Here contains a negative imaginary part, and is a normalization factor giving the number of contractions related to the operator that multiplies in eq. (4.3):

 N1≡2,N2≡4Nc,N3≡Nc,N4≡NcC\tiny\rm{F},N5≡2Nc(Nc+1). (4.14)

If the potentials were -independent and with an infinitesimal imaginary part, i.e. , they would only induce mass shifts. In this case the spectral functions can be determined analytically,

 ρ(0)i(E′)=NiM324πθ(E′−ReVi(∞))√E′−ReVi(∞). (4.15)

This form can be used for defining generalized Sommerfeld factors:

 ¯Si≡∫∞−ΛdE′πe−E′/Tρi(E′)∫∞−ΛdE′πe−E′/Tρ(0)i(E′)=(4πMT)32∫∞−ΛdE′πe[ReVi(∞)−E′]/Tρi(E′)Ni. (4.16)

Then eq. (4.11) combined with eqs. (4.6) and (4.7) leads to a generalization of eq. (4.10),

 ⟨σ\scriptsize effv⟩=2c1+4c2Nce−ΔMT/T+[c3¯S3+c4¯S4C\tiny\rm{F}+2c5¯S5(Nc+1)]Nce−2ΔMT/T(1+Nce−ΔMT/T)2. (4.17)

If a potential 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

 ρi(E′)=πNi∑j|ψj(0)|2δ(E′−E′j)∫d3r|ψj(r)|2, (4.18)

where are the bound state wave functions. Inserting into eq. (4.16), the contribution of the th bound state to reads

 Δj¯Si=(4πMT)32|ψj(0)|2e[ReVi(∞)−E′j]/T∫d3r|ψj(r)|2. (4.19)

This becomes (exponentially) large when , however chemical equilibrium is lost in the dark sector at low , which imposes an effective cutoff on the growth (cf. secs. 5 and 6).

### 4.5 Thermal potentials

In order to write down the potentials appearing in eq. (4.12), let us define

 v\tiny\rm{ }(r) ≡ (4.20) = g2s2×⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩exp(−m\tiny\rm{D}r)4πr−iT2πm\tiny% \rm{D}r∫∞0dzsin(zm\tiny\rm{D}r)(z2+1)2,r>0−m\tiny\rm{D}4π−iT4π,r=0. (4.21)

The integrand in eq. (4.20) corresponds to the static limit of the time-ordered HTL-resummed temporal gluon propagator. Then we find

 V1(r) = 0,V2(r)=C\tiny\rm{F}v(0),V3(r)=2C\tiny\rm{F}[v(0)−v(r)], (4.22) V4(r) = 2C\tiny\rm{F}v(0)+v(r)Nc,V5(r)=2C\tiny\rm{F}v(0)+(Nc−1)v(r)Nc. (4.23)

The structure equals the combination shown in eqs. (3.5)–(3.7), whereas yields the Salpeter part of in eq. (4.8). The potential corresponds to a singlet potential, to an octet potential, and 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  GeV, when the Higgs mechanism is operative, additional potentials can be generated, particularly through the Higgs portal coupling in eq. (2.1) (cf. e.g. ref. [33]). However the coefficients of these potentials are suppressed by , where is the Higgs expectation value. Given that we consider  TeV, we expect their contributions to be negligible compared with QCD effects and have not included them. We also note that an -dependence can be generated for through quark exchange, however this is suppressed by .

For a practical use of eq. (4.21), numerical values are needed for the parameters and . We relegate a discussion of this point into appendix A. Let us however note that we restrict to temperatures  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 for .

## 5 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 , where is the entropy density, and change variables from time to , whereby eq. (4.5) becomes

 Y′(z)=−⟨σ\scriptsize effv⟩Mm% \scriptsize Pl×c(T)√24πe(T)×Y2(z)−Y2\scriptsize eq(z)z2∣∣ ∣ ∣∣T=M/z. (5.1)

Here is the Planck mass, is the energy density, and is the heat capacity, for which we use values from ref. [36] (cf. also ref. [37]). The final value yields the energy fraction .

We integrate eq. (5.1) up to . At around these temperatures, depending on the value of , the processes of interest have either ceased to be active, or are falling out of chemical equilibrium, because their rates are suppressed by . Therefore they cannot be reliably addressed within the current framework.

In fig. 2(left) we show the spectral function 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 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 . However these couplings only have a modest effect if chosen otherwise, because they do not affect the coefficient  appearing the attractive channel, cf. eq. (4.4). As an example, setting increases the abundance typically by %, cf. fig. 3. The most important role is played by the coupling . For this coupling has been evaluated at the scale , whereas for collider phenomenology its value at a scale would be more relevant. The latter can be obtained from eq. (A.7), and is some tens of percent smaller than . We stress that, as shown by eq. (A.7), Yukawa couplings always generate a non-zero value for through renormalization group running.

## 6 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 -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 [9, 10, 11, 7, 12, 13, 14] by incorporating the full spectrum of thermally broadened bound states as well as the effect of soft scatterings. In general such scatterings dominate interaction rates at small mass splittings, because they are not phase-space suppressed in the same way as 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  TeV, provided that the mass splitting is small, , and that the “Higgs portal” coupling between the coloured scalar and the Higgs doublet is substantial. We recall that in supersymmetric theories, is proportional to the quark Yukawa coupling squared, , 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, , is negative. If it overcompensates for the mass difference, so that , 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  could be substantially larger than in fig. 4, for instance  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, and . The roles that these play are rather different. The value of at the scale influences the coefficient which mediates the most efficient annihilations, cf. eq. (4.4). In contrast affects the rate of transitions between the Majorana fermions and coloured scalars, cf. eq. (3.9), as well as the running of , cf. eq. (A.7). As long as 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.

## Acknowledgements

This work was supported by the Swiss National Science Foundation (SNF) under grant 200020-168988.

## Appendix Appendix A Fixing of vacuum and thermal couplings

We start by listing the 1-loop renormalization group (RG) equations satisfied by the model of sec. 2. Apart from the couplings shown in eq. (2.1), the Higgs self-coupling , the Higgs mass parameter , the weak and strong gauge couplings , and the Yukawa coupling  associated with the quark flavour  appear. The hypercharge coupling is omitted for simplicity. The number of colours is denoted by , and , whereas , and refer to the numbers of fermion generations, strongly interacting scalar triplets, and weakly interacting scalar doublets, respectively.

Parametrizing the renormalization scale through

 t≡ln¯μ2, (A.1)

we find

 ∂tμ2\tiny\rm{H} = 1(4π)2{[6λ1−9g2w4+|h|2Nc]μ2\tiny\rm{H}+λ3NcM2η}, (A.2) ∂tM2η = 1(4π)2{[2λ2(Nc+1)−3g2sC\tiny\rm{F}+|y|2]M2η+2λ3μ2\tiny\rm{H}−2|y|2M2}, (A.3) ∂tM2 = 1(4π)2{|y|2NcM2}, (A.4) ∂tλ1 = 1(4π)2{[12λ1−9g2w2+2|h|2Nc]λ1+λ23Nc2+9g4w16−|h|4Nc}, (A.5) ∂tλ2 = 1(4π)2{[2λ2(Nc+4)−6g2sC\tiny\rm{F}+2|y|2]λ2 (A.6) +λ23+3(N3c+N2c−4Nc+2)g4s8N2c−|y|4}, ∂tλ3 = 1(4π)2{[6λ1+2λ2(Nc+1)+2λ3−9g2w4−3g2sC\tiny\rm{F}+|y|2+|h|2Nc]λ3 (A.7) −2|h|2|y|2}, ∂t|y|2 = |y|2(4π)2{|y|2(Nc+3)2+|h|2−3g2sC\tiny\rm{F}}, (A.8) ∂t|h|2 = |h|2(4π)2{|h|2(2Nc+3)2+|y|22−9g2w4−6g2sC\tiny\rm{F}}, (A.9) ∂tg2w = g4w(4π)2{n\tiny\rm{W}6+4n\tiny\rm{G}3−223}, (A.10) ∂tg2s = g4s(4π)2{n\tiny\rm{S}6+4n\tiny\rm{G}3−11Nc3}. (A.11)

We note in particular that a non-zero value is generated for by the running induced by Yukawa couplings, cf. eq. (A.7).

The only coupling that we need at a scale is the strong coupling. Since it has a large influence, we evaluate it at 2-loop level for (nowadays running is known up to 5-loop level [39, 40, 41]). Denoting by the number of flavours and setting for brevity, the 2-loop running is given by

 ∂tas = −{β0a2s+β1a3s+…}, (A.12) as ≡ g2s4π2,β0=14{11−2Nf3},β1=142{102−38Nf3}. (A.13)

The value of is changed when a quark mass threshold is crossed at , where continuity is imposed. The initial value is . For , 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 . Since parametrically only the scales 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]

 m2\tiny\rm{D}≡T2[g2sα¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny\rm{MS}% \scriptsize E4+g4s(4π)2α¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny\rm{MS}\scriptsize E6+O(g6s)],g2\tiny\rm{E}≡g2s+g4s(4π)2α¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯% \tiny\rm{MS}\scriptsize E7+O(g6s). (A.14)

For general masses, only and are available at present:

 α¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny\rm{MS}\scriptsize E4 = Nc3+4Nf∑i=1[F2(m2iT2,0)−m2i(4π)2T2F3(m2iT2,0)], (A.15) α¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny\rm{MS}\scriptsize E7 = 22Nc3[ln(¯μeγ\tiny\rm{E}4πT)+122]−23Nf∑i=1[θ(¯μ−mi)ln(¯μ2m2i)+F3(m2iT2,0)]. (A.16)

Here the functions read (; chemical potentials have been set to zero)

 F2(y,0) ≡ 14π2∫∞0