A very heavy sneutrino as viable thermal dark matter candidate in U(1)′ extensions of the MSSM

We study the Standard Model singlet (“right-handed”) sneutrino ν˜R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\tilde{\nu}}_R $$\end{document} dark matter in a class of U(1)′ extensions of the MSSM that originate from the breaking of the E6 gauge group. These models, which are referred to as UMSSM, contain three right-handed neutrino superfields plus an extra gauge boson Z′ and an additional SM singlet Higgs with mass ≃ MZ′, together with their superpartners. In the UMSSM the right sneutrino is charged under the extra U(1)′ gauge symmetry; it can therefore annihilate via gauge interactions. In particular, for Mν˜R≃MZ′/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {M}_{{\tilde{\nu}}_R}\simeq {M}_{Z\prime }/2 $$\end{document} the sneutrinos can annihilate by the exchange of (nearly) on-shell gauge or Higgs bosons. We focus on this region of parameter space. For some charge assignment we find viable thermal ν˜R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\tilde{\nu}}_R $$\end{document} dark matter for mass up to ∼ 43 TeV. This is the highest mass of a good thermal dark matter candidate in standard cosmology that has so far been found in an explicit model. Our result can also be applied to other models of spin−0 dark matter candidates annihilating through the resonant exchange of a scalar particle. These models cannot be tested at the LHC, nor in present or near-future direct detection experiments, but could lead to visible indirect detection signals in future Cherenkov telescopes.


Introduction
Supersymmetry (SUSY) is one of the best motivated theories to describe new physics beyond the Standard Model (SM) at the TeV scale. It introduces a space-time symmetry that relates bosons and fermions which can be used to cancel the quadratic divergences that appears in the radiative corrections of the masses of scalar bosons, providing thus a natural solution to the hierarchy problem of the SM. It allows the gauge couplings to unify at a certain grand unified scale in the vicinity of the Planck scale [1][2][3] and this can be seen as a clear hint that SUSY is the next step towards a grand unified theory (GUT) [4].
One of the most interesting features of low energy supersymmetric models with conserved R parity is that the lightest supersymmetric particle (LSP) is absolutely stable and behaves as a realistic weakly interacting massive particle (WIMP) dark matter (DM) candidate [5,6]. Since WIMPs have non-negligible interactions with SM particles, they can be searched for in a variety of ways. Direct WIMP search experiments look for the recoil of a nucleus after elastic WIMP scattering. These experiments have now begun to probe quite deeply into the parameter space of many WIMP models [7,8]. The limits from these experiments are strongest for WIMP masses around 30 to 50 GeV. For lighter WIMPs the JHEP04(2019)167 recoil energy of the struck nucleus might be below the experimental threshold, whereas the sensitivity to heavier WIMPs suffers because their flux decreases inversely to the mass.
It is therefore interesting to ask how heavy a WIMP can be. As long as no positive WIMP signal has been found, an upper bound on the WIMP mass can only be obtained within a specific production mechanism, i.e. within a specific cosmological model. In particular, nonthermal production from the decay of an even heavier, long-lived particle can reproduce the correct relic density for any WIMP mass, if the mass, lifetime and decay properties of the long-lived particle are chosen appropriately [9]. Here we stick to standard cosmology, where the WIMP is produced thermally from the hot gas of SM particles. The crucial observation is that the resulting relic density is inversely proportional to the annihilation cross section of the WIMP [10]. It has been known for nearly thirty years that the unitarity limit on the WIMP annihilation cross section leads to an upper bound on its mass [11]. Using the modern determination of the DM density [12], the result of [11] translates into the upper bound m χ ≤ 120 TeV . (1.2) While any elementary WIMP χ has to obey this bound, it is not very satisfying. Not only is the numerical value of the bound well above the range that can be probed even by planned colliders; a particle that interacts so strongly that the annihilation cross section saturates the unitarity limit can hardly be said to qualify as a WIMP. In order to put this into perspective, let us have a look at the upper bound on the WIMP mass in specific models.
An SU(2) non-singlet WIMP can annihilate into SU(2) gauge bosons with full SU(2) gauge strength. For a spin−1/2 fermion and using tree-level expressions for the cross section, this will reproduce the desired relic density (1.1) for m χ 1.1 TeV for a doublet (e.g., a higgsino-like neutralino in the MSSM [13]); about 2.5 TeV for a triplet (e.g., a wino-like neutralino in the MSSM [13]); and 4.4 TeV for a quintuplet [14]. Including large one-loop ("Sommerfeld") corrections increases the desired value of the quintuplet mass to about 9.6 TeV [15].
One way to increase the effective WIMP annihilation cross section is to allow for coannihilation with strongly interacting particles [16]. Co-annihilation happens if the WIMP is close in mass to another particle χ , and reactions of the kind χ + f ↔ χ + f , where f, f are SM particles, are not suppressed. In this case χχ and χ χ annihilation reactions effectively contribute to the χ annihilation cross section. If χ transforms non-trivially under SU(3) C , the χ χ annihilation cross section can be much larger than that for χχ initial states. On the other hand, χ then effectively also counts as Dark Matter, increasing the effective number of internal degrees of freedom of χ. For example, in the context of the MSSM, co-annihilation with a stop squark [17] can allow even SU(2) singlet (bino-like) DM up to about 3.3 TeV [18], or even up to ∼ 6 TeV if the mass splitting is so small that the lowest stoponium bound state has a mass below twice that of the bino [19]. Coannihilation with the gluino [20] can put this bound up to ∼ 8 TeV [21]. Very recently it JHEP04(2019)167 has been pointed out that nonperturbative co-annihilation effects after the QCD transition might allow neutralino masses as large as 100 TeV if the mass splitting is below the hadronic scale [22]; the exact value of the bound depends on non-perturbative physics which is not well under control.
The WIMP annihilation cross section can also be greatly increased if the WIMP mass is close to half the mass of a potential s-channel resonance R. Naively this can allow the cross section to (nearly) saturate the unitarity limit, if one is right on resonance. In fact the situation is not so simple [16], since the annihilation cross section has to be thermally averaged: because WIMPs still have sizable kinetic energy around the decoupling temperature, this average smears out the resonance. In the MSSM the potentially relevant resonances for heavy WIMPs are the heavy neutral Higgs bosons; in particular, neutralino annihilation through exchange of the CP-odd Higgs A can occur from an S-wave initial state [23]. However, the neutralino coupling to Higgs bosons is suppressed by gauginohiggsino mixing; it will thus only be close to full strength if the higgsino and gaugino mass parameters are both close to M A /2.
In this paper we therefore focus on models where the MSSM is extended by an extra U(1) group, yielding the "UMSSM". This not only gives rise to a new gauge boson Z , but also to an additional Higgs field s whose vacuum expectation value (VEV) breaks the additional U(1) . This can provide a natural solution to the µ problem of the MSSM [24] where the µ term is generated dynamically by the VEV of s [25]. Although this solution is similar to the one provided by the next-to-minimal supersymmetric standard model (NMSSM) [26], the UMSSM is free of the cosmological domain wall problem because the U(1) symmetry forbids the appearance of domain walls which are created by the Z 3 discrete symmetry of the NMSSM [27].
For concreteness we work in the E 6 inspired version of the UMSSM. Models of this kind were first studied more than 30 years ago in the wake of the first "superstring revolution" [28,29]. This framework allows to study a wide range of U(1) groups, since E 6 contains two U(1) factors beyond the SM gauge group.
We also add three SM singlet right-handed neutrino superfieldsN C i to the spectrum. Their fermionic members are needed to cancel anomalies related to the U(1) . Moreover, the scalar membersν R,i of these superfields make good WIMP candidates [30][31][32][33][34][35][36]. This is in contrast with the left-handed sneutrinos of the MSSM, which have been ruled out as DM candidates by direct WIMP searches because their scattering cross sections on nuclei are too large [37]. Right-handed sneutrinos have small scattering cross sections on nuclei. Moreover, being scalar SU(2) singlets, a right-handed sneutrino only has two degrees of freedom; in contrast, a higgsino-like neutralino, which also has unsuppressed couplings to the Z boson in many cases, effectively has eight (an SU(2) doublet of Dirac fermions, once co-annihilation has been included).
While the new Higgs superfieldŜ is a singlet under the SM gauge group, it is charged under U(1) . This forbids anŜ 3 term in the superpotential. Hence the quartic scalar interaction of this field is determined uniquely by its U(1) charge. As a result, the mass of the physical, CP-even Higgs boson h 3 is automatically very close to that of the Z boson, in the relevant limit M Z M Z . Hence for Mν R,1 M Z /2 the annihilation cross section

JHEP04(2019)167
of the lightest right-handed sneutrinoν R,1 is enhanced by two resonances. Out of those, the exchange of h 3 is more important since it can be accessed from an S-wave initial state. For a complex scalar, Z exchange is accessible only from a P -wave initial state, which suppresses the thermally averaged cross section. Notice that the h 3νR,iν * R,i coupling contains terms that are proportional to the VEV of s, which sets the scale of the Z mass; for Mν R,1 M Z /2 this dimensionful coupling therefore does not lead to a suppression of the cross section. Finally, the couplings of h 3 to the doublet Higgs bosons can be tuned by varying a trilinear soft breaking term. This gives another handle to maximize the thermally averagedν R,1 annihilation cross section in the resonance region.
The remainder of this paper is organized as follows. In section 2 we describe the theoretical framework of the UMSSM and discuss its particle content, with a particular emphasis on the gauge, Higgs, sneutrino and neutralino sectors. We describe the calculation of the relic density, and explain our procedure to minimize it in section 3. In section 4 we first present the results of our numerical analysis for two specific U(1) models derived from E 6 and then, after following the same procedure in other U(1) models, we show the distribution of the DM upper limit of the RH sneutrino mass in the whole UMSSM; prospects of probing such scenarios experimentally are also discussed. Finally, section 5 summarizes and concludes the paper.

Model description
We focus on Abelian extensions of the MSSM with gauge group SU(3) C ×SU(2) L ×U(1) Y × U(1) , which can result from the breaking of the E 6 gauge symmetry [28,29]. In other words, it can be seen as the low energy limit of a -possibly string-inspired -E 6 grand unified gauge theory. E 6 contains SO(10) × U(1) ψ and, since SO(10) can be decomposed into SU(5) × U(1) χ where SU(5) contains the entire gauge group of the SM, one can break Here we assume that only one extra U(1) factor survives at the relevant energy scale, which in general is a linear combination of U(1) ψ and U(1) χ , parameterized by a mixing angle θ E 6 [29] . The U(1) charges of all the fields contained in the model are then given by where Q ψ and Q χ are the charges associated to the gauge groups U(1) ψ and U(1) χ , respectively. In addition to the new vector superfieldB and the MSSM superfields, the UMSSM contains one electroweak singlet supermultipletŜ ≡ (s,s), with a scalar field s that breaks the U(1) gauge symmetry, and three RH neutrino supermultipletsN c i ≡ (ν c R , ν c R ) i . In table 1 we give the U(1) charge of all relevant matter and Higgs fields in the UMSSM for certain values of the mixing angle θ E 6 . It should be noted that U(1) ψ and U(1) χ are both anomaly-free over complete (fermionic) representations of E 6 . Since U(1) χ is a subgroup JHEP04(2019)167   of SO (10), which is also anomaly-free over complete representations of SO (10), and the SM fermions plus the right-handed neutrino complete the 16-dimensional representation of SO(10), U(1) χ is anomaly-free within the fermion content we show in the table. However, U(1) ψ will be anomaly-free only after we include the "exotic" fermions that are contained in the 27-dimensional representation of E 6 , but are not contained in the 16 of SO(10).
Here we assume that these exotic superfields are too heavy to affect the calculation of thẽ ν R,1 relic density. We will see that this assumption is not essential for our result. Figure 1 shows these charges as functions of the mixing angle θ E 6 . We identify by vertical lines values of θ E 6 that generate the well-known U(1) groups denoted by U(1) ψ , U(1) N , U(1) I , U(1) S , U(1) χ and U(1) η . The black curve in figure 1 shows that for θ E 6 = arctan glets, and do not couple to any potential s-channel resonance. Similarly, for θ E 6 = 0, i.e. U(1) = U(1) χ , the charge ofŜ vanishes; in that case s cannot be used to break the gauge symmetry, i.e. the field content we have chosen is not sufficient to achieve the complete breaking of the (extended) electroweak gauge symmetry down to U(1) QED . All other values of θ E 6 are acceptable for us.
The superpotential of the UMSSM contains, besides the MSSM superpotential without µ term, a term that couples the extra singlet superfield to the two doublet Higgs superfields; this term is always allowed, since it is part of the gauge invariant 27 3 of E 6 . The superpotential also contains Yukawa couplings for the neutrinos. We thus have: where · stands for the antisymmetric SU(2) invariant product of two doublets. The neutrino Yukawa coupling Y ν is a 3×3 matrix in generation space and λ is a dimensionless coupling. Note that for θ E 6 = 0 the U(1) symmetry forbids both bilinearN C iN C j and trilinear SN C iN C j terms in the superpotential. In this model the neutrinos therefore obtain pure Dirac masses, which means that the Yukawa couplings Y ν,ij must be of order 10 −11 or less; in our numerical analysis we therefore set Y ν = 0.
The electroweak and the U(1) gauge symmetries are spontaneously broken when, in the minimum of the scalar potential, the real parts of the doublet and singlet Higgs fields acquire non-zero vacuum expectation values. These fields are expanded as We define tan β = vu v d and v = v 2 d + v 2 u exactly as in the MSSM; this describes the breaking of the SU(2) L × U(1) Y symmetry, and makes subleading contributions to the breaking of U(1) . The latter is mostly accomplished by the VEV of s. The coupling λ in eq. (2.3) then generates an effective µ−term: As well known, supersymmetry needs to be broken. We parameterize this by soft breaking terms [38]:

JHEP04(2019)167
Here we have used the notation of SPheno [39,40]. The soft scalar masses and the soft trilinear parameters of the sfermions are again 3 × 3 matrices in generation space. In the UMSSM, the Bµ term of the MSSM is induced by the T λ term after the breaking of the U(1) gauge symmetry.
In the following subsections we discuss those parts of the spectrum in a bit more detail that are important for our calculation. These are the sfermions, in particular sneutrinos; the massive gauge bosons; the Higgs bosons; and the neutralinos. The lightest right-handed sneutrino is assumed to be the LSP, which annihilates chiefly through the exchange of massive gauge and Higgs bosons in the s-channel. Requiring the lightest neutralino to be sufficiently heavier than the lightest right-handed sneutrino gives important constraints on the parameter space. The mass matrices in these subsections have been obtained with the help of the computer code SARAH [41][42][43]; many of these results can also be found in refs. [34][35][36].

Sfermions
In the UMSSM, the U(1) gauge symmetry induces some new D-term contributions to the masses of all sfermions with nonvanishing U(1) charges. These modify the diagonal entries of the MSSM sfermion mass matrices: where g is the U(1) gauge coupling and F ∈ {Q, L, The first two terms on the right-hand side (r.h.s.) of eq. (2.7) are therefore essentially negligible. However, due to the contribution ∝ v 2 s these D-terms can dominate the sfermion masses. Moreover, depending on the value of θ E 6 these terms can be positive or negative. For example, figure 1 shows that for arctan √ 15 < θ E 6 < π 2 all the sfermion masses receive positive corrections, the corrections to the RH sneutrino masses being the smallest ones. In contrast, for 0 < θ E 6 < arctan √ 15 the D-term contribution to the RH sneutrino masses is negative. For θ E 6 < 0 the RH sneutrino masses again receive positive corrections from this D-term.
The tree-level sneutrino mass matrix written in the basis (ν L ,ν R ) is The 3 × 3 sub-matrices along the diagonal are given by: where g 1 and g 2 are the U(1) Y and SU(2) L gauge couplings, respectively. As noted earlier, the neutrino Yukawa couplings have to be very small. We therefore set Y ν = T ν = 0, so

JHEP04(2019)167
that the 6×6 matrix (2.8) decomposes into two 3×3 matrices. 1 Since all interactions of thẽ ν R fields are due to U(1) gauge interactions which are the same for all generations, we can without loss of generality assume that the matrix m 2 N C of soft breaking masses is diagonal. The physical masses of the RH sneutrinos are then simply given by m 2 Our LSP candidate is the lightest of the threeν R states, which we callν R,1 .

Gauge bosons
The UMSSM contains three neutral gauge bosons, from the SU(2) L , U(1) Y and U(1) , respectively. As in the SM and MSSM, after symmetry breaking one linear combination of the neutral SU(2) L and U(1) Y gauge bosons remains massless; this is the photon. The orthogonal state Z 0 mixes with the U(1) gauge boson Z 0 via a 2 × 2 mass matrix: Recall that g 2 , g 1 and g are the gauge couplings associated to SU(2) L , U(1) Y and U(1) , respectively. The eigenstates Z and Z of this mass matrix can be written as: 14) The mixing angle α ZZ is given by The masses of the physical states are Note that the off-diagonal entry ∆ Z in eq. (2.10) is of order v 2 . We are interested in Z masses in excess of 10 TeV, which implies v 2 , which is automatically below current limits [7] if M Z ≥ 10 TeV. Moreover, mass mixing increases the mass of the physical Z boson only by a term of order M 4 Z /M 3 Z , which is less JHEP04(2019)167 than 0.1 MeV for M Z ≥ 10 TeV. To excellent approximation we can therefore identify the physical Z mass with M Z 0 given in eq. (2.13), with the last term ∝ v 2 s being the by far dominant one.
Recall from the discussion of the previous subsection that the mass of the right-handed sneutrinos can get a large positive contribution from the U(1) D-term for some range of θ E 6 . In fact, from eqs. (2.7) and (2.13) together with the charges listed in table 1 we find that this D-term contribution exceeds ( For this range of θ E 6 one therefore needs a negative squared soft breaking contribution m 2Ñ Note that we neglect kinetic Z −Z mixing [36,44]. In the present context this is a loop effect caused by the mass splitting of members of the 27 of E 6 . This induces small changes of the couplings of the physical Z boson, which have little effect on our result; besides, this loop effect should be treated on the same footing as other one-loop corrections.

The Higgs sector
The Higgs sector of the UMSSM contains two complex SU(2) L doublets H u,d and the complex singlet s. Four degrees of freedom get "eaten" by the longitudinal components of W ± , Z and Z . This leaves three neutral CP-even Higgs bosons h i , i ∈ {1, 2, 3}, one CP-odd Higgs boson A and two charged Higgs bosons H ± as physical states. After solving the minimization conditions of the scalar potential for the soft breaking masses of the Higgs fields, the symmetric 3 × 3 mass matrix for the neutral CP-even states in the basis (φ d , φ u , φ s ) has the following tree-level elements: In general the eigenstates and eigenvalues of this mass matrix have to be obtained numerically. We denote the mass eigenstates by h 1 , h 2 , h 3 , ordered in mass.
The tree level mass of the single physical neutral CP-odd state is

JHEP04(2019)167
In our sign convention, tan β and v s are positive in the minimum of the potential; eq. (2. 19) then implies that T λ must also be positive. As in the MSSM M 2 A differs from the squared mass of the physical charged Higgs boson only by terms of order v 2 : (2.20) Both A and H ± are constructed from the components of H u and H d , without any admixture of s.
Recall that we are interested in the limit v s v. Eqs. (2.18) show that the entries mixing the SU(2) L singlet with the doublets are of order v s v or T λ v, whereas the diagonal mass of the singlet is of order v 2 s . The mixing between singlet and doublet states is therefore small. Moreover, we will work in the MSSM-like decoupling limit M 2 A M 2 Z , which ensures that the lightest neutral CP-even Higgs boson has couplings close to those of the SM Higgs. Its mass can then approximately be written as [36,45]: (2.21) The first term on the r.h.s. is as in the MSSM. The second term is an F -term contribution that also appears in the NMSSM, while the third term is due to the U(1) D− term. These terms are positive. The second line is due to mixing between singlet and doublet states; note that this mixing always reduces the mass of the lighter eigenstate, but increases the mass of the heaviest state h 3 . As well known, the mass of h 1 also receives sizable loop corrections, in particular from the top-stop sector [46,47]; we will briefly discuss them below when we describe our numerical procedures. As noted above, in the limit v s v the mixing between singlet and doublet states can to first approximation be neglected. Here we chose the heaviest state to be (mostly) singlet. From the last eq. (2.18) and eq. (2.13) we derive the important result Here we have assumed |T λ | ≤ v s because for larger values of |T λ | the mass of the heavy doublet Higgs can exceed the mass of the singlet state. As we will see, in the region of parameter space that minimizes theν R,1 relic density we need M A < Mν R,1 . Eq. (2.22) leads to an h 3 − Z mass splitting of order M 2 Z /M Z , which is below 1 GeV for M Z > 10 TeV. Loop corrections induce significantly larger mass splittings, with M Z > M h 3 ; however, the splitting still amounts to less than 1% in the relevant region of parameter space, which is well below the typical kinetic energy of WIMPs in the epoch around their decoupling from the thermal bath. We thus arrive at the important result that Mν R,1 M Z automatically implies Mν R,1 M h 3 in our set-up, so thatν R,1 annihilation is enhanced by two nearby resonances.

Neutralinos
The neutralino sector is formed by the fermionic components of the neutral vector and Higgs supermultiplets. So, in addition to the neutralino sector of the MSSM, the UMSSM has another gaugino state associated with the U(1) gauge symmetry and a singlino state that comes from the extra scalar supermultipletŜ. The neutralino mass matrix written in the basis λB,W 0 ,H 0 d ,H 0 u ,S, λB is: This matrix is diagonalized by a unitary 6 × 6 matrix N which gives the mass eigenstates (in order of increasing mass)χ 0 1 ,χ 0 2 ,χ 0 3 ,χ 0 4 ,χ 0 5 ,χ 0 6 as a linear combinations of the current eigenstates. We have ignored a possible (gauge invariant) mixedBB mass term [48,49].
Note that the singlet higgsino (singlino for short)S and the U(1) gauginoB mix strongly, through an entry of order v s . On the other hand, these two new states mix with the MSSM only through entries of order v. Therefore the eigenvalues of the lower-right 2×2 submatrix in eq. (2.23) are to good approximation also eigenvalues of the entire neutralino mass matrix. Note that the smaller of these two eigenvalues decreases with increasing M 4 . Requiring this eigenvalue to be larger than Mν R,1 M Z /2 therefore implies Moreover, the smallest mass of the MSSM-like states should also be larger than M Z /2, which implies We have used eqs. (2.5) and (2.13) in the derivation of the last inequality. We finally note that the chargino sector of the UMSSM is identical to that of the MSSM, with µ → µ eff .
3 Minimizing the relic abundance of the right-handed sneutrino As described in the Introduction, we want to find the upper bound on the mass of the lightest RH sneutrinoν R,1 from the requirement that it makes a good thermal WIMP in standard cosmology. As well known [10], under the stated assumptions the WIMP relic density is essentially inversely proportional to the thermal average of its annihilation cross section into lighter particles; these can be SM particles or Higgs bosons of the extended sector. The upper bound on Mν R,1 will therefore be saturated for combinations of parameters that maximize the thermally averagedν R,1ν * R,1 annihilation cross section. All relevant couplings of the RH sneutrinos are proportional to the U(1) gauge coupling g . In particular, two RH sneutrinos can annihilate into two neutrinos through exchange of a U(1) gaugino. This, and similar reactions where one or both particles in the initial and final state are replaced by antiparticles, are typical electroweak 2 → 2 reactions without enhancement factors. They will therefore not allow RH sneutrino masses in the multi-TeV range.

JHEP04(2019)167
In contrast,ν R,1ν * R,1 annihilation through Z and scalar h 3 exchange can be resonantly enhanced if Mν R,1 M Z /2; recall that M h 3 M Z is automatic in our set-up, if h 3 is mostly an SM singlet, as we assume. Note that the Z exchange can only contribute if the sneutrinos are in a P -wave. This suppresses the thermal average of the cross section by a factor ≥ 7. For comparable couplings, h 3 exchange, which is depicted in figure 2, is therefore more important.
In the h 3 resonance region the annihilation cross section scales like Since the h 3νR,1ν * R,1 coupling, denoted by C 1 in figure 2, originates from the U(1) D-term, it is proportional to the product of S and N C charges: These charges are determined uniquely once the angle θ E 6 has been fixed. The denominator in eq. (3.1) results from dimensional arguments, using the fact that there is essentially only one relevant mass scale once the resonance condition has been imposed. 2 Note that near the resonance the annihilation cross section is effectively only O(α ), not O(α 2 ), where α = g 2 /(4π). Theν R,1ν * R,1 annihilation cross section is then larger than typical co-annihilation cross sections, if the latter are not resonantly enhanced. Even co-annihilation with a superparticle that can also annihilate resonantly (e.g., a higgsinolike neutralino) will not increase the effective annihilation cross section, but will increase the effective number of degrees of freedom per dark matter particle g χ . As a result, we find that co-annihilation reduces the upper bound on Mν R,1 . For example, if all three RH JHEP04(2019)167 sneutrinos have the same mass, the upper bound on this mass decreases by a factor of √ 3, since the annihilation cross section has to be increased by a factor of 3 in order to compensate the increase of g χ . We therefore require that the lightest neutralino is at least 20% heavier thanν R,1 .
As noted earlier, the initial-state coupling C 1 in figure 2 is essentially fixed by θ E 6 . The upper bound on Mν R,1 for given θ E 6 can therefore be found by optimizing the final state couplings. We find that the relic density is minimized if the effective final state coupling C 2 , defined more precisely below, is of the same order as C 1 . This can be understood as follows. For much larger values of C 2 the width of h 3 increases, which reduces the cross section. On the other hand, since the peak of the thermally averaged cross section is reached for Mν R,1 slightly below M h 3 /2 [16], h 3 →ν R,1ν * R,1 decays are allowed, and dominate the total h 3 width if C 2 C 1 ; in this case increasing C 2 will clearly increase the cross section, i.e. reduce the relic density.
The only sizable couplings of the singlet-like Higgs state h 3 to particles with even Rparity (i.e., to particles possibly lighter than the LSPν R,1 ) are to members of the Higgs doublets. h 3 couples to H u and H d through the U(1) D-term, with contributions ∝ g 2 Q S Q Hu,H d v s ; through F -terms associated to the coupling λ, with contributions ∝ λ 2 v s ; and through a trilinear soft breaking term, with contributions ∝ T λ . In the decoupling limit M 2 A M 2 Z the relevant couplings are given by: (2) L is effectively unbroken. The couplings of h 3 to two members of the heavy doublet containing the physical states H ± , h 2 and A therefore are all the same, see eq. (3.3), as are the couplings to the light doublet containing h 1 and the would-be Goldstone modes G 0 and G ± , see eq. (3.4); finally, eq. (3.5) describes the common coupling to one member of the heavy doublet and one member of the light doublet. Of course, the would-be Goldstone modes are not physical particles; however, again since M h 3 v the production of physical longitudinal gauge bosons can to very good approximation be described as production of the corresponding Goldstone states. This is the celebrated equivalence theorem [50][51][52][53]. 3

JHEP04(2019)167
We find numerically that theν R,1 relic density is minimized when h 3 decays into two members of the heavy Higgs doublet are allowed. From eqs. (2.19) and (2.20) we see that this requires This implies that the singlet-like state is indeed the heaviest physical Higgs boson.
We can now define an effective final-state coupling C 2 for the diagram shown in figure 2: Here we have included the kinematic factors into the effective coupling, using the same mass Since the contribution from h 3 exchange is accessible from an S-wave initial state, it peaks for DM mass very close to M h 3 /2 where one needs quite small velocity to get exactly to the pole s = M 2 h 3 ; at such a small velocity, the Z exchange contribution, which can only be accessed from a P -wave initial state, is quite suppressed. As a consequence, near the peak of the thermally averaged total cross section the h 3 exchange processes always contributes more than 90% to the total, whereas the Z exchange contribution shrinks as we approach the peak. The latter reaches its maximum at a larger difference between M Z and 2Mν R,1 , but its contribution exceeds 10% of the total only if 2Mν R,1 is at least 3% below M Z , or else above the resonance. Note also that the annihilation into pairs of SM fermions via Z exchange is completely determined by θ E 6 . In principle we could contemplate annihilation into exotic fermions, members of 27 of E 6 that are required for anomaly cancellation, as noted in section 2.1. However, the contribution from the SM fermions already sums to an effective final state coupling which is considerably larger than the initial state coupling; this helps to explain why the Z contribution is always subdominant. Adding additional final states therefore reduces the Z exchange contribution to theν R,1 annihilation cross section even further. This justifies our assumption that the exotic fermions are too heavy to affect the calculation of theν R,1 relic density.
Finally, all other processes of the model contribute at most 1% to the thermally averaged total cross section in the resonance region. This shows that the parameters that describe the rest of the spectrum are irrelevant to our calculation, as long asν R,1 is the LSP and sufficiently separated in mass from the other superparticles to avoid co-annihilation. These parameters were therefore kept fixed in the numerical results presented below.

JHEP04(2019)167 4 Numerical results
We are now ready to present numerical results. We will first describe our procedure. Then we discuss two choices for θ E 6 , i.e. for the U(1) charges, before generalizing to the entire range of possible values of this mixing angle.

Procedure
We have used the Mathematica package SARAH [41][42][43] to generate routines for the precise numerical calculation of the spectrum with SPheno [39,40]. This code calculates by default the pole masses of all supersymmetric particles and their corresponding mixing matrices at the full one-loop level in the DR scheme. SPheno also includes in its calculation all important two-loop corrections to the masses of neutral Higgs bosons [54][55][56]. The dark matter relic density and the dark matter nucleon scattering cross section relevant for direct detection experiments are computed with micrOMEGAS-4.2.5 [57]. The mass spectrum generated by SPheno is passed to micrOMEGAS-4.2.5 through the SLHA+ functionality [58] of CalcHep [59,60]. The numerical scans were performed by combining the different codes using the Mathematica tool SSP [61] for which SARAH already writes an input template.
SARAH can generate two different types of templates that can be used as input files for SPheno. One is the high scale input, where the gauge couplings and the soft SUSY breaking parameters are unified at a certain GUT scale and their renormalization group (RG) evolution between the electroweak, SUSY breaking and GUT scale is included. The other one is the low scale input where the gauge couplings, VEVs, superpotential and soft SUSY breaking parameters of the model are all free input parameters that are given at a specific renormalization scale near the sparticle masses, in which case no RG running to the GUT scale is needed. In this template the SM gauge couplings are given at the electroweak scale and evolve to the SUSY scale through their RGEs. The dark matter phenomenology of a model in the WIMP context is usually well studied at low energies; moreover, acceptable low energy phenomenology for both the U(1) ψ and the U(1) η model in the limit where the singlet Higgs decouples works much better with nonuniversal boundary conditions [62]. Finally, a bound that is valid for general low-scale values of the relevant parameters will also hold (but can perhaps not be saturated) in constrained scenarios.
In our work we therefore define the relevant free parameters of the UMSSM directly at the SUSY mass scale, which is defined as the geometric mean of the two stop masses. We created new model files for different versions of the UMSSM to be used in SARAH and SPheno where all the U(1) charges are written in terms of the U(1) mixing angle θ E 6 using eq. (2.1).
Our goal is to find the upper bound on the mass of the lightest RH sneutrino, and therefore on M Z M h 3 . We argued in section 3 that co-annihilation would weaken the bound. We therefore have to make sure that all other superparticles are sufficiently heavy so that they do not play a role in the calculation in the relic density. The precise values of their masses are then irrelevant to us. We therefore fix the soft mass parameters of the gauginos and sfermions to certain values well above Mν R,1 ; recall from eq. (2.24) that this implies an upper bound on the mass M 4 of the U(1) gaugino. As noted in section 2 we set JHEP04(2019)167 Y ν = 0, since the small values of the neutrino masses force them to be negligible for the calculation of the relic density. We also set most of the scalar trilinear couplings to zero, except the top trilinear coupling T t which we use together with tan β and M 3 to keep the SM Higgs mass in the range 125 ±3 GeV, where the uncertainty is dominated by the theory error [63]. Since we are interested in superparticle masses in excess of 10 TeV, the correct value of M h 1 can be obtained with a relatively small value of tan β, which we also fix.
As already noted in the previous section, all relevant interactions ofν R,1 scale (either linearly or quadratically) with the U(1) gauge coupling g . Since our set-up is inspired by gauge unification, we set this coupling equal to the U(1) Y coupling in GUT normalization, i.e.
Note also that the charges in table 1 are normalized such that Y 2 , where the sum runs over a complete 27-dimensional representation of E 6 [29]. We will later comment on how the upper bound on Mν R,1 changes when g is varied.
Recalling that we work in a basis where the matrix m 2Ñ C is diagonal, with m 2Ñ C ,11 being its smallest element, the remaining relevant free parameters are thus: All these parameters are related to the extended sector that the UMSSM has in addition to the MSSM. Since the mixing angle θ E 6 defines the U(1) gauge group, we want to determine the upper bound on the mass of the lightest RH sneutrino as a function of θ E 6 .
We will see below that this will also allow to derive the absolute upper bound, valid for all versions of the UMSSM.
From the discussion of the previous section we know that the first two of the parameters listed in (4.2) are strongly correlated by the requirement that Mν R,1 is close to M Z /2. More precisely, the minimal relic density is found if the RH sneutrino mass is very roughly one h 3 decay width below the nominal pole position, the exact distance depending on the couplings C 1 and C 2 ; this shift from the pole position is due to the finite kinetic energy of the sneutrinos at temperatures around the decoupling temperature [16].
The parameters λ and T λ have to satisfy some bounds. First, requiring the mass of the SU(2) L higgsinos to be at least 20% larger than M Z /2 leads to the lower bound where we have used eqs. (2.5) and (2.13). Moreover, T λ has to satisfy the upper bound (3.6), so that pairs of the heavy SU(2) L doublet Higgs bosons can be produced inν R,1 annihilation with Mν R,1 M Z /2. Having fixed tan β and T λ , the effective final state coupling C 2 defined in eq. (3.7) depends only on λ, which is constrained by eq. (4.3); fortunately this still leaves us enough freedom to vary C 2 over a sufficient range.

JHEP04(2019)167
The bound on the lightest RH sneutrino mass for a given value of θ E 6 can then be obtained as follows. We start by choosing some value of M h 3 M Z in the tens of TeV range. Note that this fixes the coupling C 1 , since we have already fixed g and θ E 6 and hence the charge Q N C . We then minimize the relic density for that value of M h 3 by varying the soft-breaking contribution to the sneutrino mass and λ; as noted in section 3, the minimum is reached when the physical RH sneutrino mass is just slightly below M Z /2, and C 2 is close to the initial state coupling C 1 of eq. (3.2). If the resulting relic density (Ωh 2 ) 1 is very close to the measured value of eq. (1.1), we have found the upper bound on M Z and hence on Mν R,1 . Otherwise, we change the value of M h 3 by the factor 0.12/(Ωh 2 ) 1 , and repeat the procedure. Since the minimal relic density to good approximation scales like TeV; recall that the physical squared sfermion masses also receive D-term contributions, which amount to M 2 Z /8 in this model. In this model the two Higgs doublets have the same U(1) charge, and the product Q Hu Q S is negative. As a result, the λ 2 and the g 2 terms in the diagonal couplings given in eqs. (3.3) and (3.4) tend to cancel, while the contribution ∝ g 2 to the off-diagonal couplings given in eq. (3.5) vanishes. The contributions involving these off-diagonal couplings are therefore subdominant. The largest contribution usually comes from final states involving two heavy SU(2) L doublet Higgs bosons, but the contributions from two light states (including the longitudinal modes of the gauge bosons) are not much smaller. Moreover, due to this cancellation we need relatively large values of λ; the numerical results shown below have been obtained by varying it in the range from 0.32 to 0.46. Figure 3a depicts the relic abundance of the RH sneutrino as a function of Mν R,1 for different values of the mass of the singlet Higgs boson. All the curves show a pronounced minimum when Mν R,1 is very close to but below M h 3 /2. The blue and the green curves are for v s = 59 TeV and thus have the same coupling C 1 and (approximately) the same mass of the singlet Higgs, but the blue curve has a smaller value of C 2 . This reduces the width of h 3 as well as the annihilation cross section away from the resonance, and therefore leads to a narrower minimum.
In figure 3b we show the dependence of the relic density on the ratio of couplings C 2 /C 1 for fixed mediator masses close to the resonance. This confirms our expectations from the previous section: if C 2 is significantly larger than C 1 , the relic density increases with C 2 because the increase of the mediator decay width over-compensates the increased coupling strength in the total annihilation cross section. If relic density because it increases the normalization of the annihilation cross section. Note that the relic density curve is fairly flat over some range of C 2 /C 1 . Moreover, the optimal choice of C 2 /C 1 also depends somewhat on how far Mν R,1 is below M h 3 /2. Altogether, for given M h 3 there is an extended 1-dimensional domain in the (Mν R,1 , C 2 /C 1 ) plane over which the relic density is quite close to its absolute minimum. This simplifies our task of minimization. Note also that we calculate the annihilation cross section only at tree-level; a change of the predicted relic density that is smaller than a couple of percent is therefore not really physically significant.
The parameters of the blue curve in figure 3b in fact are very close to those that maximize Mν R,1 within the U(1) ψ model, under the assumption thatν R,1 was in thermal equilibrium in standard cosmology. M max ν R,1 11.5 TeV corresponds to an upper bound on M h 3 and M Z of about 23.0 TeV. This is clearly beyond the reach of the LHC, and might even stretch the capabilities of proposed 100 TeV pp colliders.
Recall that all left-handed SM (anti)fermions have the same U(1) ψ charge. As a result, in the absence of Z − Z mixing the Z ff couplings are purely axial vector couplings, for all SM fermions f . Z exchange can therefore only contribute to spin-dependent WIMPnucleon scattering in this model. Since our WIMP candidate doesn't have any spin, Z exchange does not contribute at all. Once Z −Z mixing is included, Z exchange contributes a term of order Mν R,1 M N sin α ZZ /M 2 Z ∝ Mν R,1 M N /M 2 Z to the matrix element forν R,1 N scattering, while the mixing-induced Z exchange contribution is suppressed by another factor M 2 Z /M 2 Z ; here M N is the mass of the nucleon. There is also a small contribution from the light SM-like Higgs boson h 1 , which is very roughly of order M 2 N /M 2 h 1 . As a result the scattering cross section on nucleons is very small, below 10 −13 pb for the scenario that maximizes Mν R,1 . For the given large WIMP mass, this is not only several orders of magnitude below the current bound, but also well below the background from coherent neutrino scattering ("neutrino floor").

The U(1) η model
We now consider a value of θ E 6 with a larger U(1) charge of the right-handed neutrino superfields. This increases the coupling C 1 for given M Z , and thus theν R,1 annihilation cross section for given masses, which in turn will lead to a weaker upper limit on Mν R,1 from the requirement that theν R,1 relic density not be too large.
In our analysis we therefore choose the SUSY breaking scale to be 50 TeV and we fix tan β = 2.2, and m 2 Q = 1.28 × 10 9 GeV 2 · 1, m 2 U C = 1.45 × 10 9 GeV 2 · 1, m 2 D C = 3.0 × 10 9 GeV 2 · 1, m 2 L = 3.0 × 10 9 GeV 2 · 1, m 2 E C = 1.28 × 10 9 GeV 2 · 1, m 2Ñ In this model the two Higgs doublets have different U(1) charges; hence there is a sizable gauge contribution to the off-diagonal couplings of eq. (3.5). The Higgs doublet charges again have the opposite sign as the charge of S, leading to cancellations between the λ 2 and g 2 terms in the diagonal couplings (3.3) and (3.4). This cancellation is particularly strong for the coupling to two light states, so that for the interesting range of λ the most important final states involve two heavy SU(2) L doublets, although final states with one light and one heavy boson are also significant. Partly because of this, and partly because the coefficients of the g 2 terms are smaller than in the U(1) ψ model, smaller values of the coupling λ are required; the numerical results below have been obtained with λ ∈ [0.260, 0.352].
In figure 4 we again show the dependence of the relic density on the mass of the lightest RH sneutrino (left) and on the ratio of couplings C 2 /C 1 (right). The qualitative behavior is similar to that in the U(1) ψ model depicted in figure  Since Q Q = Q U C = Q D C in this model, there is no vector coupling of the Z to up quarks, but such a coupling does exist for down quarks. Hence now the Z exchange contribution to the matrix element for elastic scattering ofν R,1 on nucleons is comparable to that of Z exchange once Z −Z mixing has been included, and the h 1 exchange contribution has roughly the same size as in the U(1) ψ model. The totalν R,1 N scattering cross sections are again below 10 −13 pb, for parameters near the upper bound on Mν R,1 . Since our WIMP candidate is now even heavier than in the U(1) ψ model, this is even more below the current constraints as well as below the neutrino floor.

The general UMSSM
In this subsection we investigate in more detail how the upper bound on Mν R,1 depends on θ E 6 . To this extent we have applied the procedure outlined in subsection 4.1, and applied to two specific U(1) models in subsections 4.2 and 4.3, to several additional U(1) models, each with a different value of θ E 6 .
The results are shown in figure 5, where we plot the upper bound on the mass of the lightest RH sneutrino as a function of the absolute value of the product g Q N C . In order of increasing |Q N C |, the six red points correspond to the following choices of θ E 6 : Note that the first point has a vanishing U(1) charge for the N C superfields, i.e. the resonance enhancement of the annihilation cross section does not work in this case. We

JHEP04(2019)167
checked that the cross section for elasticν R,1 N scattering is well below the experimental bound for all other points.
Evidently the upper bound on Mν R,1 scales essentially linearly with Q N C ; recall that g has been fixed to 5/3g 1 here. This linear dependence can be understood as follows. The h 3νR,1ν * R,1 coupling can be written as g Q N C M Z 2g Q N C Mν R,1 . Moreover, we saw above that the maximal sneutrino mass is allowed if the effective final-state coupling C 2 is similar to C 1 ; it is therefore also proportional to Q N C . Therefore at the point where the bound is saturated, the h 3 decay width scales like |C 1 | 2 M h 3 ∝ g 2 (Q N C ) 2 Mν R,1 , where we have again used that near the resonance all relevant masses are proportional to Mν R,1 . Note finally that for a narrow resonance -such as h 3 , for the relevant parameter choices -the thermal average over the annihilation cross section scales like 1/(M h 3 Γ h 3 ) [16]. Altogether we thus have The linear relation between the upper bound on Mν R,1 and Q N C then follows from the fact that the thermally averaged annihilation cross section essentially fixes the relic density.
Note that here Q N C always comes with a factor g ; indeed, for a U(1) gauge interaction only the product of gauge coupling and charge is well defined. The linear dependence of the bound on Mν R,1 on Q N C for fixed g depicted in figure 5 can therefore also be interpreted as linear dependence of the bound on the product g Q N C . A fit to the points in figure 5 gives: M max ν R,1 = (0.071 + 113.477g |Q N C |) TeV . Finally, we recall from eq. (2.17) that for θ E 6 between − arctan √ 15 and 0 one needs a negative squared soft breaking mass in order to have Mν R,1 M Z /2. Since theN C superfields appear in the superpotential (2.3) only multiplied with the tiny couplings Y ν , this superpotential will not allow to generate negative squared soft breaking masses for sneutrinos via renormalization group running starting from positive values at some high scale. If we insist on positive squared soft breaking mass for allν R fields the upper bound on |Q N C | is reduced to 5/8 0.79, in which case the bound on Mν R,1 is reduced to about 42 TeV. We note, however, that theN C superfields can have sizable couplings to some of the exotic color triplets that reside in the 27-dimensional representation [28]. Recalling that at least some of these exotic fermions are usually required for anomaly cancellation it should not be too difficult to construct a UV complete model that allows negative squared soft breaking terms for (some)ν R at the SUSY mass scale.

Prospects for detection
Clearly spectra near the upper bound presented in the previous subsection are not accessible to searches at the LHC, nor even to a proposed 100 TeV pp collider.
As already noted for the U(1) η and U(1) ψ models theν R,1 nucleon scattering cross section is very small. The very large Z mass suppresses the Z exchange contribution; as we saw in section 2.3 it also suppresses Z −Z mixing, so that the Z exchange contribution also scales like M −2 Z . The contribution due to the exchange of the singlet-like Higgs boson (h 3 in our analysis) is suppressed by the very large value of M h 3 as well as the tiny h 3 qq couplings, which solely result from mixing between singlet and doublet Higgs bosons. Finally, the contribution from the exchange of the doublet Higgs bosons, in particular of the 125 GeV state h 1 , is suppressed by the small size of the h 1νR,1ν * R,1 coupling, which is of order g v Mν R,1 , as well as the rather small h 1 qq couplings, which are much smaller than gauge couplings. As a result, theν R,1 nucleon scattering cross section, and hence the signal rate in direct WIMP detection experiments, is well below the neutrino-induced background; recall that this "neutrino floor" increases ∝ Mν R,1 since the WIMP flux, and hence the event rate for a given cross section, scales ∝ 1/Mν R,1 .
The best chance to test these scenarios therefore comes from indirect detection. Naively one expects the cross section for annihilation from an S-wave initial state to be essentially independent of temperature, in which case the correct thermal relic density implies σv 2.4 · 10 −26 cm 3 /s 0.8 pb · c [64,65]. However, as pointed out in [66,67] this can change significantly in the resonance region; here the thermally averaged annihilation cross section can be significantly higher in today's universe than at the time of WIMP decoupling.

JHEP04(2019)167
This is illustrated in figure 6 for the parameter choice that saturates the upper bound on Mν R,1 in the U(1) η model. Here we show the thermally averagedν R,1ν * R,1 annihilation cross section times relative velocity as function of the scaled inverse temperature x = Mν R,1 /T . 4 We see that for a quite extended range of temperatures around the decoupling temperature, σv grows almost linearly with x. This is because Mν R,1 is only slightly below the nominally resonant value M h 3 /2; by reducing the temperature the fraction of the velocity distribution that falls within approximately one h 3 decay width of the pole therefore at first increases.
Today's relic density is essentially inversely proportional to the "annihilation integral", defined as [16] J (4.6) An annihilation cross section that grows significantly for x > x F therefore has to be compensated by a smaller value of σv (x F ) in order to keep the relic density constant. As a result, in our scenarios the annihilation cross section at decoupling is actually significantly smaller than for typical S-wave annihilation.
Because for parameters that saturate the upper bound on Mν R,1 the right-handed sneutrino mass is somewhat below M h 3 /2, for very large x, i.e. very small temperature, the thermally averaged annihilation cross section starts to decrease again. However, for the parameters of figure 6 it asymptotes to a value that is still about three times larger than the "canonical" thermal WIMP annihilating from an S-wave initial state. As shown in refs. [66,67] this enhancement factor strongly depends on 2Mν R,1 − M h 3 ; it can be even larger for slightly smaller sneutrino masses that are even closer to M h 3 /2.
The WIMP annihilation rate in today's universe scales like the square of the WIMP number density. This means that the flux of annihilation products scales like 1/M 2 ν R,1 ; for parameters (nearly) saturating our upper bound on the sneutrino mass it is thus too small to be detectable by space-based observatories like FermiLAT [68], simply because of their small size. Recall also that our sneutrinos annihilate into (longitudinal) gauge or Higgs bosons, and thus mostly into multi-hadron final states. This leads to a continuous photon spectrum which, for parameters near the upper bound on the sneutrino mass, extends well into the TeV region. Photons of this energy can be detected by Cherenkov telescopes on the ground, via their air showers. Note also that the astrophysical cosmic ray background drops even faster than E −2 with increasing energy E of the cosmic rays; the signal to background ratio therefore actually improves with increasing WIMP mass. Indeed, simulations show that at least for a favorable distribution of dark matter particles near the center of our galaxy, the continuum photon flux of multi-TeV WIMPs annihilating with the canonical thermal cross section should be detectable by the Cherenkov Telescope Array [69].

JHEP04(2019)167 5 Summary and conclusions
In this paper we analyzed the UMSSM, i.e. extensions of the minimal supersymmetrized Standard Model that contain an additional U(1) gauge group as well as additional righthanded (RH) neutrino superfields which are singlets under the SM gauge group but carry U(1) charge. We assume that U(1) is a subgroup of E 6 , which has been suggested as an (effective) GUT group, e.g. in the context of early superstring phenomenology. In this case the lightest RH sneutrinoν R,1 can be a good dark matter candidate.
We found that even within minimal cosmology, and fixing the U(1) gauge strength to be equal to that of the hypercharge interaction of the (MS)SM (in GUT normalization),ν R,1 masses of tens of TeV are possible. For given U(1) charges the bound on Mν R,1 is saturated ifν R,1 can annihilate resonantly through the exchange of both the new Z gauge boson and of the new Higgs boson h 3 associated with the spontaneous breaking of U(1) ; note that M Z M h 3 automatically in this model. Scalar h 3 exchange is more important since Z exchange can only occur from a P -wave initial state. The h 3νR,1ν * R,1 coupling is fixed by the U(1) charge Q N C of the right-handed neutrinos, but the h 3 couplings to the relevant final states can be tuned independently, allowing a further maximization of the annihilation cross section. In our analysis we used SU(2) L doublet Higgs bosons as well as longitudinal W and Z bosons as final states. While the light SU(2) doublet Higgs states, including the longitudinal W and Z modes, are always accessible, we could have replaced the heavy Higgs doublet in the final state by some exotic fermions which in most cases are required to cancel anomalies. The only requirement is that the effective final state coupling of h 3 should be tunable to values close to its coupling toν R,1 . Since the Z exchange contribution is basically fixed by θ E 6 , and non-resonant contributions are negligible for Mν R,1 ∼ M Z /2, most of the many free parameters of this model, which describe the sfermion and gaugino sectors, are essentially irrelevant to us. The only requirement is that these superparticles are sufficiently heavy to avoid co-annihilation, which would increase the relic density in our case.
We found that the final upper bound on Mν R,1 is essentially proportional to the product g |Q N C |, where g is the U(1) gauge coupling. Within the context of theories unifiable into E 6 this leads to an absolute upper bound on Mν R,1 of about 43.8 TeV. In other words, in this fairly well motivated set-up we can find a thermal WIMP candidate with mass less than a factor of three below the bound derived from unitarity [11]. This is to be contrasted with an upper bound on the mass of a neutralino WIMP in the MSSM of about 8 TeV for unsuppressed co-annihilation with gluinos [21]. In a rather more exotic model featuring a WIMP residing in the quintuplet representation of SU(2) a WIMP mass of up to 9.6 TeV is allowed [15].
Of course, this mechanism requires some amount of finetuning: the mass of the WIMP needs to be just below half the mass of the s-channel mediator. We find that typically the predicted WIMP relic density increases by a factor of 2 when the WIMP mass is reduced by between 1 and 3% from its optimal value. In contrast, the recent proposal to allow thermal WIMP masses near 100 TeV via non-perturbative co-annihilation requires finetuning to less than 1 part in 10 5 [22].

JHEP04(2019)167
We also note that our very heavy WIMP candidates have very small scattering cross sections on nuclei, at least two orders of magnitude below the neutrino floor. This shows that both collider searches and direct WIMP searches are still quite far away from decisively probing this reasonably well motivated WIMP candidate. On the other hand, we argued that indirect signals for WIMP annihilation might be detectable by future Cherenkov telescopes. Our analysis thus motivates extending the search for a continuous spectrum of photons from WIMP annihilation into the multi-TeV range.
While the result (4.5) has been derived within UMSSM models that can emerge as the low-energy limit of E 6 Grand Unification, it should hold much more generally. To that end g |Q N C | should be replaced by g χχφ /m φ , where χ is a complex scalar WIMP annihilating through the near resonant exchange of the real scalar φ, g χχφ being the (dimensionful) χχ * φ coupling. In order to saturate our bound the couplings of φ to the relevant final states should be tunable such that the effective final state coupling, which we called C 2 in section 3, should be comparable to the initial-state coupling g χχφ . In this case the algorithm we used to find the upper limit on Mν R,1 , see subsection 4.1, can directly be applied to finding the upper bound on M χ . We finally note that M χ can be increased by another factor of √ 2 if χ is a real scalar.