Bound-state effects for dark matter with Higgs-like mediators

In this paper we study the impact of a scalar exchange on the dark matter relic abundance by solving a plasma-modified Schroedinger equation. A simplified model is considered where a Majorana dark matter fermion is embedded in a U(1)$'$ extension of the Standard Model and couples with a dark Higgs via a Yukawa interaction. We find that the dark-Higgs exchange can increase the overclosure bounds significantly. For the largest (smallest) value of the Yukawa coupling examined in this work, the dark matter mass is lifted from 5 TeV (0.55 TeV) to 27 TeV (0.70 TeV).


Introduction
Dark matter is one of the main open problems in the realm of cosmology and particle physics. If dark matter is assumed to be a particle rather than an astrophysical object, the hypothesis of a weakly interacting massive particle (WIMP) has been certainly the most studied. This choice does not fix a unique candidate though, on the contrary a plethora of possible dark matter particles are available [1,2]. The quest for a successful candidate poses interesting connections between the machinery of quantum field theory, needed to calculate dark matter annihilation and scattering rates, and the many constraints imposed from the astrophysical and Earth-based experimental measurements. This has resulted in highly constrained scenarios: the viable parameter space of a given model is often in tension with that needed to reproduce the observed dark matter relic abundance via the so called freeze-out mechanism (see e.g. ref. [3] for a comprehensive status on WIMPs). Here, the key ingredient is the annihilation cross section of dark matter pairs that enters a Boltzmann equation and eventually determines the freeze-out abundance [4][5][6]. The latter has to match with the accurate measurement of the dark-matter energy density Ω DM h 2 = 0.1186 ± 0.0020 [7].
Recently, simplified models have been suggested for the interpretation of beyond the Standard Model searches at colliders, direct and indirect detection experiments [8][9][10]. In this framework, rather than considering a fully fledged theory, bounds and constraints are set on a simple model that captures the most relevant physics. Reinterpreting the experimental results in terms of simplified models, strong lower bounds are currently being set by recent analyses at the Large Hadron Collider (LHC) [11] and the XENON1T experiment [12,13] that look for the footprint of a new massive particle. Within simplified models, one is able to classify in a systematic way the nature of new degrees of freedom that may play the role of a dark matter particle, together with accompanying particles of the new physics model. Indeed in many cases, the so-called mediators act as portals between the dark and visible sector (it is also possible to have more than one mediator), preserve unitarity and gauge invariance, and enrich the phenomenology.
When moving to such realistic particle models, some processes may occur that call for revisiting the standard relic abundance calculation, i.e. the derivation of the annihilation cross section in the early universe. For example, potential-like interactions are induced by a sufficiently light vector or scalar mediator (lighter than the dark matter mass) together with the possibility of bound-state formation. For a mediator mass comparable with the dark matter mass, coannihilations can play an important role and the mediator can itself experience soft interactions if coupled with light Standard Model degrees of freedom. Thermal masses and thermal interaction rates may also be important, the latter can lead to boundstate formation/dissociation in a thermal bath. The inclusion of some of these effects has led to substantial revision of the overclosure bound for a given dark matter model, namely the largest value of the particle mass compatible with the observed dark matter energy density.
In particular, the electroweak gauge boson exchange and gluon exchange can be important and the corresponding Sommerfeld enhancement has been included in the annihilation cross section in many studies, e.g. [14][15][16][17]. The inclusion of bound-state effects in the annihilation process through a Boltzmann equation is rather non-trivial and different approaches have been put forward lately [18][19][20][21][22][23][24][25].
A non-perturbative formalism for addressing the thermal annihilation of non-relativistic particles has been developed quite recently [21,26]. In this context, the thermally averaged annihilation cross section is obtained in terms of a chemical equilibration rate [27], the latter extracted from correlators evaluated in equilibrium and independent of the assumptions typical of a Boltzmann description. The key ingredient is the imaginary part of a two-point Green's function, namely a spectral function. The advantage of using such an approach is twofold: (i) the spectral function can be determined by solving a thermallymodified Schrödinger equation with static potentials that comprise several in-medium effects like virtual and real scatterings; (ii) the appearance of bound states is naturally described in this framework and the need of complicated bound-state production and dissociation rates is avoided. This formalism has been applied to the Inert Doublet Model and to a simplified model comprising a Majorana fermion coannihilating with a strongly interacting scalar, where weakly and strongly bound states appear respectively [24,25].
Potential-like interactions arise naturally when considering a fermion or a scalar dark matter coupled to gauge bosons (due to the trilinear vertex in the covariant derivative). However, it is also possible to have a scalar exchange between dark matter pairs, such as the Standard Model Higgs boson or the corresponding Higgs boson of the new physics model. In the latter case, we refer to it as dark Higgs throughout the paper. The effect of the Higgs boson exchange has been studied for the Inert Doublet Model with a focus on dark-matter annihilations leading to gamma ray signals [28], together with an estimate of the impact on cross sections in the early universe. Similar analyses have been carried out for scalar and fermionic dark matter with a Higgs portal [17,[29][30][31]. In all cases, the Sommerfeld effect has been studied that affects the dark matter pair wave function at zero temperature. In this work, we aim to apply the aforementioned finite-temperature formalism [21,26] to assess the formation of bound state induced by a scalar exchange besides the Sommerfeld enhancement. We shall work in the framework of simplified models. The bulk of the analysis is carried out for a model with a spontaneously broken U(1) gauge symmetry that contains a Majorana dark matter fermion, a dark gauge boson and a dark Higgs [32][33][34][35]. In addition, we elaborate on an another model of recent interest, namely a Majorana dark matter coannihilating with a coloured scalar charged under QCD and interacting with the SM Higgs boson [36].
The plan of the paper is the following. In section 2.1, we discuss the simplified model that we focus on, i.e. a U(1) extension of the SM. In section 2.2 the thermally averaged annihilation cross section is presented within an effective field theory approach. Then we derive the non-relativistic Lagrangian in section 2.3, the thermal potentials are given in section 2.4, whereas the plasma-modified Schrödinger equation is discussed in section 2.4 together with numerical outputs for the overclsoure bound. We consider other simplified models where a Higgs exchange can appear in section 3. Finally some conclusions and discussion are offered in section 4.

Majorana fermion dark matter and U(1) gauge symmetry
We want to study dark matter models where a scalar field can be exchanged between the dark matter particles. As a well-motivated and interesting example, we pick the simplified model recently described in refs. [34,35] that realizes perturbativity and gauge invariance at the same time.

Model description and light-mediators regime
The model contains a dark Higgs and a dark gauge boson in addition to a Majorana fermion dark matter (the latter is assumed to be the actual dark matter particle that contributes to the present universe energy density). The dark Higgs provides the mass of both the dark matter fermion and the dark gauge boson via the spontaneous breaking of the U(1) symmetry. Portal couplings induce an interaction between the dark and the SM sector (scalar mixing and gauge boson mixing). The Lagrangian of the model reads [34,35] where χ is a Majorana fermion field, V µ is the dark gauge boson, S is the dark Higgs field, H is the Standard Model Higgs doublet, Then f is a generic Standard Model fermion that couples via a vector current with the U(1) boson and q f is the corresponding charge. The fermion dark matter couples to the dark gauge boson with an axial-vector current (the vector current vanishes for a Majorana fermion and this choice helps in suppressing direct detection cross section with respect to the Dirac case). The covariant derivative acting on the dark Higgs field reads D µ = ∂ µ + ie q s V µ . In order to write the gauge invariant mass term for the Majorana dark matter in the first line of eq. (2.1), we have to require q s = −2q χ [34]. Then we define g χ ≡ e q χ and therefore D µ = ∂ µ − 2ig χ V µ . In the following we neglect the portal couplings λ hs and κ.
An important observation is that the couplings between the dark matter and the dark Higgs and the dark matter and the dark gauge bosons are not independent [34,35]. Indeed, after the U(1) symmetry breaking, S = (w + s + iϕ)/ √ 2, the two masses read in the T = 0 and they are related to each other as According to the global analysis given in [35], the model is rather unconstrained by experiments in the region where M χ > m V , m s . This is also the situation where one expects the dark-Higgs and dark-vector exchange to have some impact. Moreover, from eq. (2.3), one can see that requiring M χ m V implies y χ g χ . This suggests that the coupling between the dark Higgs and the dark matter is larger than the one between the dark matter and the dark gauge boson. This is a hint to motivate the inspection of dark-Higgs exchange diagram, see figure 1. Furthermore, we also ask the dark matter to be heavier than the scalar mass. We can use the relation in the T = 0 limit 1 and then pick the appropriate values for the couplings to fix the desired ratio M χ /m s 1. Let us stress that, in this particular model, the dark matter mass is provided by the spontaneous breaking of the U(1) gauge symmetry. Therefore, only the broken phase is relevant to us in order to study the freeze-out mechanism: the dark matter has to acquire a finite mass M χ , attain thermal equilibrium and enter a non-relativistic regime when its mass drops below the plasma temperature. Eventually it decouples around T ∼ M χ /25...M χ /20 like in the standard WIMP scenario. However, we notice that in the case λ hs = 0 the dark and Standard Model Higgs expectation values are coupled and their evolution with temperatures may not be trivial (see appendix in ref. [34] for more details on the scalar mixing).

Dark-matter annihilations in a thermal bath
Our aim is to describe accurately dark matter pair annihilations and include systematically near-threshold effects in a finite temperature environment, most importantly bound-state formation. Soft exchange processes are mediated by the dark Higgs and the gauge boson.
First, let us summarize the framework of the freeze-out of a heavy thermal relic that puts us in a deep non-relativistic regime. The dark matter particles are kept in chemical equilibrium through interactions with the thermal bath until T M χ ≡ M and gradually freeze out at temperatures T ∼ M/25. Annihilations continue even during later stages where the dark matter particles are still in kinetic equilibrium. In this situation most of the energy of a dark matter particle is given by its mass and, for non-relativistic species, the typical momentum is |p| = √ M T = M T /M . One usually identifies an average velocity v ≡ T /M , which is smaller than unity in the regime of interest. Therefore, the degrees of freedom during freezeout annihilations are non-relativistic Majorana fermions, for which M T , light Standard Model and dark particles (the dark Higgs and gauge boson).
In order to make manifest the non-relativistic nature of the dark matter, one may write down a non-relativistic Lagrangian from the start. Moreover, non-relativistic particle annihilations can be described by four-particle operators O i , arranged as an expansion in 1/M 2 . The prototype for such effective field theory (EFT) is the well-known non-relativistic QCD (NRQCD) [37]. The small parameter of the effective field theory is the average velocity v 1 of the heavy particles, here the dark-matter Majorana fermions. In the EFT language, hard energy/momentum modes of order M are integrated out from the fundamental theory (2.1). We write the low-energy Lagrangian explicitly in the next section 2.3. The major benefit of the EFT formulation is to separate two classes of processes: those occurring at the hard scale M , and those typical of the soft scales, either thermal or non-relativistic. Indeed, given the large energy release in the annihilation process, the typical distance scales are much smaller than those introduced by the thermal plasma, i.e. ∆x ∼ 1/M 1/T . Many scales remain still dynamical in the so-obtained low-energy theory: the thermal scales πT and gT , and the non-relativistic scales M α and M α 2 . 2 At smaller energy scales, the heavy pairs can be sensitive to medium effects and a quantum statistical interpretation of pair annihilations is desirable. Since dark matter particles are slowly moving, repeated soft interactions can occur that are mediated by the dark Higgs and dark gauge boson. These interactions, that can modify the wave function of the annihilating dark matter pair, happen in a thermal bath. Hence, correlators should be evaluated within finite temperature field theory. It comes as the main strength of the approach exploited here [21,[24][25][26] to recast the partition function of the annihilating pair as the thermal expectation value of the four-particle operators. This way one can dynamically account for the whole two-particle spectrum, both scattering and bound states properly weighted by the corresponding Boltzmann factor, and include near-threshold soft effects for which T ∼ M α 2 . Bound states have an effect of order unity for such temperatures that is reflected in the Boltzmann factor of the annihilating pair (cfr. eq. (2.16)). For a more detailed and comprehensive discussion see refs. [21,24].
In summary, we shall compute the thermally averaged annihilation cross section as σv = T . First, we have to derive the matching coefficients c i of non-relativistic four-particle operators O i that create and annihilate dark matter pairs. In a second stage, we shall compute the thermal average of the very same four-particle operators O i that amounts to solve a thermally modified Schrödinger equation for the dark matter pair with the thermal potentials of the mediators (see section 2.4 and 2.5). Finally, the extraction of the corresponding spectral function comprises the information on the annihilating states in the statistical ensemble, i.e. scattering states and bound states.

Non-relativistic Lagrangian
In this section we outline the vertices between the heavy Majorana dark matter and the light degrees of freedom, namely the dark gauge boson V µ , the Goldstone boson ϕ and the dark Higgs s in the low-energy theory. This is the field theory that comprises energy modes with typical energies smaller than the dark matter mass. In addition we also write the fourparticle operators describing the heavy Majorana fermion pair annihilations. We write the non-relativistic Majorana fermion as follows (we choose the standard parametrization of the Dirac matrices 3 ) where the Grassmanian spinor ψ has two components. Starting from the interaction Lagrangian after the U(1) symmetry breaking, the terms which have no fast oscillations read The superscript stands for non-relativistic (NR) and σ i are the Pauli matrices. The Majorana fermion does not show any interaction with the temporal component of the gauge boson V 0 , 3 We take the following assignment: at variance with what happens in the case of heavy Dirac fermions interacting with gauge bosons via vector like currents, such as in the well-known Heavy Quark Effective Theory (HQEFT) [38], NRQCD [37] and potential NRQCD [39,40]. In our case, only the spatial components of the gauge field interact with the non-relativistic spinor. An EFT approach for NR Majorana fermions has been previously introduced in ref. [41]. Then we write down the absorpative Lagrangian that comprises the four-particle operators of the effective theory. Dealing with a Majorana fermion, there is only one operator at order 1/M 2 which describes the dark matter annihilation [25] and we find the following matching coefficient 4 According to the optical theorem, the imaginary part of the one-loop diagrams with fourparticles external legs is equivalent to the matrix element squared of the annihilation processes of the type χχ → a b, where a and b are generic light degrees of freedom the heavy particles can annihilate into. Matching the four-point Green's function of the fundamental theory onto that of the low-energy theory fixes the coefficient given in eq. (2.8). This procedure is well established in the realm of non-relativistic effective field theories for QCD [37]. Since we are working with vanishing portal couplings, the possible final states are combinations of the real scalar, Goldstone boson and gauge boson referred to as dark terminators [34,35]. The annihilation cross section in the free case reads simply σv (0) = 2c 1 . For general orientation on the dark matter masses that provide the correct relic density, we anticipate some benchmark values to be M ≈ 0.5, 2, 5 TeV for y χ = 0.5, 1, 1.5 respectively and for g χ = y χ /10.

Scalar and vector induced potentials
The dark Higgs and the dark gauge boson can be exchanged between the dark matter pairs. If these particles are sufficiently lighter than the dark matter mass, they can induce sizeable effects on the scattering states, namely the Sommerfeld effect, and below threshold effects, i.e. a bound state spectrum. Moreover, thermal effects can enter such dynamics and we include them in two respects. First, we use the scalar and gauge boson propagator in the Hard Thermal Loop (HTL) approximation [42][43][44][45]. In general the so-obtained propagators contain both a thermal mass and a finite thermal width that account for virtual and real scatterings with light degrees of freedom in the thermal plasma. Second, dark matter pairs interact in a statistical background and, therefore, their dynamics is properly described by correlators evaluated in a finite temperature field theory. These can be expressed in terms of a spectral function at T = 0 that exhibits a smoothing between the bound state spectrum and the scattering states.
Since the Majorana dark matter fermion couples to the spatial components of the vector boson (see eq. (2.7)), the relevant self-energy is Π ij . It is well known that in the static limit Π ij vanishes, namely there is no thermal mass nor imaginary part at one-loop order for the spatial gauge fields. However, the gauge field has a "thermal mass" through the temperature dependence of the dark-Higgs expectation value. The temperature dependent dark-Higgs expectation value reads (see appendix B for details) from which we define m 2 V = 4g 2 χ w 2 T . When the temperature is such that w T ≤ 0 the U(1) symmetry is restored and the mass of the dark gauge boson vanishes accordingly.
As far as the dark-Higgs propagator at finite temperature is concerned, we notice that no imaginary part arises in the HTL static limit. Only a finite thermal mass appears that is related to the expectation value already written in eq. (2.10). The dark-Higgs propagator reads, in the static limit and in the imaginary-time formalism lim ω→0 i s s T (ω, k) = 1 k 2 + 2λ s w 2 where k ≡ |k| and m 2 s (T ) ≡ 2λ s w 2 T . We recall that by requiring small mediator masses, i.e. M χ m V , m s , implies the condition y χ g χ , λ s . Hence, the interaction between the dark fermion and the dark scalar is parametrically more relevant than that involving the dark fermion and the gauge boson. Therefore, we focus on the interactions induced by fermion-scalar vertex in (2.7) and we consider the corresponding diagrams in figure 1. The corresponding thermal propagator is given in eq. (2.11).
With the definition of an auxiliary potential function we write the dark-Higgs potential obtained from the diagrams in figure 1 as where we defined α y = y 2 χ /(4π). We notice that the r-dependent part is attractive and the r-independent part provides an overall negative correction to the dark matter pair self-energy, given that m s (T ) < m s . This can be traced back to a mass correction for a single dark matter particle. Moreover, the r-independent part are linearly divergent, therefore the corresponding vacuum counterterms are defined such that lim r→∞ V 1 (r) = 0 at T = 0 [24].
In the next section we study the modification to the annihilation rate induced by the potential written in eq. (2.13).

Plasma-modified Schrödinger equation and overclosure bound
In order to compute the annihilation rate for a dark matter pair as part of a thermal bath, we use the formalism developed in refs. [21,26] and already applied for two dark matter models in refs. [24,25]. At the core of the method is the extraction of a spectral function from the imaginary part of Green's functions where the thermal potential is the one given in eq. (2.13) and N i refers to the number of contractions of the four-particle operator. In this case there is only one operator with N 1 = 2, see eq. (2.8). In the potential induced by the dark Higgs there is no imaginary part within the approximation adopted in this work. However, we allow for a small imaginary part in the potential, i.e. V 1 − iΓ, in order to extract the spectral function and we set it to Γ ≈ (10 −6 -10 −5 )M . 5 In this model we have to study a single spectral function corresponding to the annihilating Majorana fermion pair. Since there is no thermal width due to the Landau damping, the shape of the spectral function is rather insensitive to the value of the temperature. The dependence on the temperature enters the thermal dark-Higgs mass and the couplings. The thermal mass m s (T ) differs from the in-vacuum mass by up to 10% depending the temperature and the model parameters. As far as y χ is concerned, one has to evaluate it in a broad range of energy scales, namely µ ≈ πT, e −γ E /r in the thermal potential (2.13) and at µ = 2M in the matching coefficient of the hard annihilation in eq. (2.9). However, the main feature of the Laplace transform for the spectral function preserves the importance of bound states, if there are any (see eq. (2.16) below). Whereas for T larger than the binding energy the main contribution to the annihilation rate is given by the above threshold region, i.e. Sommerfeld factors with appropriate thermal masses accounted for, the bound state region dominates for 5 In practice the value of Γ is is chosen to obtain numerical stability while keeping it as small as possible in order not to introduce fictitious effects. That said, it is possible to consider the decay width of the dark matter pair in the bound state. This choice has been made in the literature, see e.g. [28]. However, it does not differ much from our choice since Γ T =0 ≈ α 5 y M/2.
where α 2 y M Λ M is a cutoff restricting the average to the non-relativistic regime [21,26]. In figure 2 we show the spectral function close to threshold for three different choices of the ratio between the dark matter fermion and the dark-Higgs masses M/m s . The Yukawa coupling is chosen to be y χ = 1.5 (it corresponds to α y ≈ 0.18, pretty close to the largest value considered in ref. [17], but smaller than the maximum value considered in ref. [35], i.e. y χ = 2). A running Yukawa coupling has been included and it plays a role in a better estimation of the generalized Sommerfeld factors. Indeed, energy scales smaller than the hard annihilation scale are relevant in the Schrödinger equation, e.g. α y M and πT . The Yukawa coupling y χ decrease with the energy (see appendix A for details) at variance with what happens in QCD and for the gluon exchange. We look at a temperature around the freeze-out region, namely M/T = 25. A bound state appears and is more prominent for smaller mediator masses, respectively M/m s = 20 and M/m s = 50 for the dashed-red and dot-dashed brown lines. The corresponding Sommerfeld factors, as defined in (2.16), are shown in the right panel of figure 2. Now we can proceed to the determination of the freeze-out abundance. Within a Boltzmann equation the dark matter abundance evolves as [4,5] (we label n χ ≡ n) n = − σv (n 2 − n 2 eq ) , where the generalized Sommerfeld factorS 1 is extracted from the corresponding spectral function as in (2.16) and c 1 is from eq. (2.9). Then we define the usual yield parameter Y ≡ n/s, where s is the entropy density, and change variables from time to z ≡ M/T . Therefore eq. (2.17) becomes where m Pl is the Planck mass, e is the energy density, and c is the heat capacity, for which we use values from ref. [46]. In figure 3 we show the overclosure bounds obtained with free cross sections and those accounting for the dark-Higgs exchange. On the left plot we set y χ (2M ) = 1.0 whereas y χ (2M ) = 1.5 in the right plot. In the latter case, the dark matter mass that reproduces Ω DM h 2 = 0.1186 ± 0.0020 is lifted from M = 5. ref. [24], up to larger effects observed in the case of strong interactions [25].

Other simplified models with Higgs-like exchange
In this section we address a different simplified model that comprises a trilinear vertex between a Higgs field and a dark matter pair. The model we have in mind comprises a Majorana dark matter particle coannihilating with a coloured scalar, the latter charged under QCD (see [36] for a review of the model). Besides the interactions with gluons and the corresponding potentials, additional effects induced by the Higgs exchange can appear. We are not going to derive the overclosure bounds as systematically as in the previous case, however we make contact with some of the results derived in section 2 when possible.
We divide the discussion by following two different implementations of the interaction between the coloured scalar and the Higgs boson. First, we stick to the model we studied in [25]. In this case the interaction reads where η is the coloured scalar, H the Standard Model Higgs doublet and λ 3 a scalar coupling. This Lagrangian leads to an interaction between the coloured scalar and the Standard Model Higgs that is suppressed by v/M η after the electroweak symmetry breaking. We want to assess whether relevant contributions to the generalized Sommerfeld factors can arise from this particular realization.
Second, we start with the model as written in ref. [17], where h is taken to be a real scalar (possibly the Higgs boson) and there is no 1/M η suppression after one expands η in the non-relativistic modes. In this second option, the coupling between the real scalar and the coloured scalar is taken to be proportional to M η from the beginning (motivated by some SUSY arguments [47,48]).

Case 1
The impact of the gluon exchange for this model has been extensively studied [20,22,23,25,[49][50][51]. The effect is particularly relevant when the mass splitting between the Majorana dark matter and the coloured scalar is small (M η ≡ M + ∆M with ∆M M ), so that the dark matter abundance is actually controlled by that of the η particles, the latter experiencing strong interactions. Then Sommerfeld effects, decohering scatterings, bound state formation/dissociation have been included in the derivation of the freeze-out abundance. However, at temperatures T 160 GeV, a trilinear coupling between the coloured scalar and the Higgs boson is established.
After the electroweak symmetry breaking the trilinear vertex is given by the following non-relativistic Lagrangian which is obtained from (3.1) when expanding η = (φe −iM t + ϕ † e iM t )/ √ 2M in terms of the non-relativistic fields φ and ϕ, where φ(ϕ) annihilates a particle (antiparticle). Then h is the real scalar field corresponding to the Standard Model Higgs boson after the symmetry breaking. This vertex is suppressed by a factor v/M with respect to the gluon induced one. The vertex is controlled by the temperature dependent Higgs expectation value, namely [21] v .
At temperatures larger than T c ≈ 160 GeV the Higgs mechanism melts away and the trilinear vertex inducing the Higgs exchanges does as well. The Higgs thermal mass squared reads m 2 h (T ) = 2λv 2 T . The potential induced by the Higgs exchange can be calculated from the four-particle operators [25] At variance with the r-dependent potentials induced by the gluon, that are different for each color representation, we obtain the same scalar contribution for all the operators where we define the auxiliary thermal potential as . (3.7) The exchange diagram gives an attractive potential for all the annihilation channels, on the contrary to what happens for the gluon exchange that induces a repulsive potential in the octet and η-η operators (second, third and fourth operator in the second line of (3.5)). We define an effective coupling Exploiting the renormalization group equations (RGEs) derived in ref. [25], we explore some possibilities according to different combinations for the model couplings (these are λ 2 , λ 3 and y [36]). Moreover, we use v T as given in eq. (3.4). From figure 5 (left panel), one may see that even in the case λ 3 = π, which is a rather large value, we obtain at most α eff ≈ 0.01. Based on our previous study involving such weak-interaction values for the coupling strength [24], we conclude that the effect of the Higgs exchange cannot compete with the gluon exchange in this realization of the model. Indeed, the corresponding generalized Sommerfeld factors induces an increase of the overclosure bound of about 10%, whereas QCD strong interactions give an increase of about 200% for the same mass splitting of the co-annihilating species (see figure 6 in ref. [24] and figure 3 in ref. [25] for ∆M = 5 × 10 −3 ).

Case 2
The main reason for the smallness of the effective coupling in eq. (3.8) is the ratio v/M originating from non-relativistic Lagrangian (3.3). The simplified model considered in ref. [17] is such that this suppression is absent and the coupling between the coloured scalar and the Higgs boson is taken to be in the range α h ∈ [0.02, 0.2], where α h = g 2 h /(16π) [17,52]. The largest and smallest value correspond almost to what we have considered in Section 2, namely y χ = 1.5 and y χ = 0.5 that give α y ≈ 0.18 and α y ≈ 0.02. At this point the analysis carried out in section 2 can help in estimating the effect of the Higgs enhancement.
Let us start writing the cross section where both the gluon and Higgs exchange are included. We neglect the Majorana fermion p-wave suppressed operator and corresponding contribution to the cross section, and we then obtain where we split the thermally averaged Sommerfeld factors as (3.10) Here we writeS h in order to signal that the Higgs-induced potentials have the same form as those given in eq. (3.6), however the auxiliary potential reads in this case . (3.11) Once more we notice that the Higgs-induced generalized Sommerfeld factor is the same for all the operators and larger than unity. The thermal mass splitting entering eq. (3.9) has been derived in ref [25]. 6 In figure 5 we compare the thermally averaged Sommerfeld factors induced by the gluon exchange (we take the results from ref. [25]) with those coming from the scalar exchangē S h , for a dark matter mass of M = 3 TeV. Fixing the scalar mass to the Higgs boson mass m h = 125 GeV, we obtain the ratio M/m h = 24. We can borrow the results from the previous model, where we studied the generalized Sommerfeld factors for different dark-matter dark-Higgs mass ratio and for different Yukawa couplings. One can see that the attractive gluon exchange is already more important than the Higgs exchange for α h = 0.08 (that corresponds to y χ ≈ 1), whereas the two processes provide a rather similar generalized Sommerfeld factor 6 At the level of this study we do not need the mass splitting, however for completeness let us mention that ∆MT = ∆M +g 2 s CF T 2 /(12M )−αsCF mD(T )/2+α h (m h (T )−m h )/2, where αs = g 2 s /(4π) and α h = g 2 h /(16π). The Higgs contribution is different from ref. [25] due to the different Lagrangian.
for α h = 0.18. We note in passing that the Higgs exchange dominates over the gluon exchange in the octet channel (dashed-blue line in figure 5) as already noted in ref. [17].
A comment is in order. The mediator mass that leads to a Yukawa potential has a different origin in the two cases. On one hand, the gluon mass is purely thermal, i.e. m D ≈ g s T , and of order 10 2 GeV for temperatures around the freeze-out for a dark matter mass at a TeV range. On the other hand, the Higgs mass is m h (T ) = √ 2λv T with v T from eq. (3.4). Here the thermal contributions play a little role around the freeze-out temperature, making the in-vacuum mass the relevant mass scale of the exchanged particle.

Conclusions and discussion
In this paper we have studied the impact of a scalar exchange on the dark matter relic density. In order to quantify such an effect, we considered a simplified model with a Majorana dark matter fermion charged under a new U(1) gauge group. The dark sector is made of a dark gauge boson and a dark Higgs boson in addition to the Majorana fermion. The dark vector and scalar are the model mediators, and can possibly interact with Standard Model particles. We restrict our study to the case of vanishing portal couplings (κ = λ hs = 0 in eq. (2.1)) and we assume light mediators, M m s , m V . The latter assumption implies that the coupling between the dark matter and the dark scalar is larger than that between the dark matter and the gauge boson.
Profiting from an effecting field theory framework, we derived the non-relativistic Lagrangian that describes the interaction between the heavy dark matter fermion and the light degrees of freedom. The impact of the dark-Higgs exchange is taken into account by solving a thermally modified Schrödinger equation and extracting the generalized Sommerfeld factors from a spectral function. Then the Boltzmann equation is solved with the corresponding annihilation cross section. We scan over the Yukawa coupling y χ ∈ [0.5, 1.5] and for the dark-Higgs mass m s ∈ [M/50, M/10]. Going to lighter scalar masses, a larger impact of the scalar exchange is observed and bound states appear for sufficiently large values of the Yukawa couplings. Our results complement previous works where the scalar exchange has been considered and we add a possible treatment of bound-state effects. As already observed in the case of the gauge boson exchange (weak gauge bosons and gluons), we find that the dark matter mass reproducing the observed relic abundance is shifted to larger values with respect to the tree level one. However, this enhancement depends crucially on the parameters of the model at hand (see figures 3 and 4), namely the coupling y χ and the scalar mass m s . The generalized Sommerfeld factors obtained in this model are effective down to low temperatures because there is no suppression given by mass splittings with any coannihilating specie, i.e. e −∆/T . These results are collected in section 2.
In addition, we compared the generalized Sommerfeld factors coming from the interactions with gluons and a Standard Model Higgs boson for a different simplified model in section 3.
Here, the impact of the scalar exchange depends on how the interaction between the coloured scalar and the real scalar is implemented (see section 3.1 and 3.2). We find that the Higgs exchange can induce an effect as large as the gluon exchange and lead to bound-states formation.
Finally, let us remark that the scalar exchange can affect the overcloure bounds significantly and should be then included in the relic density calculation. For the simplified model with a Majorana dark matter and a dark Higgs, the dark matter mass is lifted from (0.55, 2.2, 5.1) TeV to (0.70, 4.1, 27.0) TeV respectively for three benchmark values y χ = (0.5, 1, 1.5) considered in this work. The parameter space that reproduces the observed dark matter abundance is rather modified and the overclosure bound is pushed to larger masses. In light of these results, it seems worth exploring the impact on direct and indirect searches for the very same model in order to better assess the reach of present and upcoming experiments (such as XENON1T [12], DARWIN [53] and CTA [54]) in the medium/high mass range.

Acknowledgements
This work was supported by the Swiss National Science Foundation (SNF) under grant 200020-168988. The author thanks Mikko Laine for valuable discussions and comments on the manuscript.

A Renormalization Group Equations
In this section we present the results for the running couplings relevant for the U(1) model. In particular, we need the coupling y χ for a broad range of energies: from E ∼ 2M , the typical energy scale of the hard annihilations, to E ∼ πT, e −γ E /r, where the very same coupling is evaluated in the thermal potential. RGEs for a similar model have been derived in ref. [55].
We use dimensional regularization in the MS with D = 4 − 2ε and compute the diagrams in the Feynman gauge. The RGEs read at one loop for the relevant couplings µ dg 2 In the running for the g χ , that is fixed by the wave-function renormalization of the vector boson, we show the different contributions explicitly (dark matter fermion, dark scalar and Standard Model fermion). In our case we have N F = 1 and N s = 1 and we set q f = 0 for all the N f Standard Model fermions. In the numerical evaluation we simply impose g χ = y χ /10, in order to satisfy the relation y χ g χ .

B Dark-Higgs self-energy
The dark scalar self-energy at finite temperature has been used in the body of the paper. The thermal self-energies for a dark gauge boson and dark scalar in a model very similar to the one we studied here can be found in ref. [26]. As far as the dark-Higgs self-energy is concerned, the diagrams are shown in figure 6 for the Feynman gauge (ghosts and Goldstone bosons are included in the diagrams). Our result agrees with that in ref. [26], upon the change e → 2e ≡ 2g χ . The self-energy in the imaginary time formalism reads in the HTL limit Π s = 8g 2 χ (D − 1) + 8λ s Then the finite temperature dark-Higgs expectation value is that gives eq. (2.10) and where m s is the T = 0 scalar mass.