Revisiting (s)neutrino dark matter in natural SUSY

We study natural supersymmetric scenarios with light right-handed neutrino superfields, and consider the possibility of having either a neutrino or a sneutrino as a dark matter candidate. For the former, we evaluate the possibility of having SUSY corrections on the ν4 → ν`γ decay rate, such that the NuStar bounds are relaxed. We find that corrections are too small. For sneutrino dark matter, we consider thermal and non-thermal production, taking into account freeze-out, freeze-in and super-WIMP mechanisms. For the non-thermal case, we find that the ν̃R can reproduce the observed relic density by adjusting the R-sneutrino mass and Yukawa couplings. For the thermal case, we find the need to extend the model in order to enhance sneutrino annihilations, which we exemplify in a model with an extended gauge symmetry.


I. INTRODUCTION
In the past decades, supersymmetric (SUSY) models have been the most popular candidates for physics beyond the Standard Model (BSM). This popularity has been very well justified, given the capacity of the Minimal Supersymmetric Standard Model (MSSM) of solving the hierarchy problem, providing a dark matter (DM) candidate, and achieving gauge coupling unification, among other features.
Unfortunately, in the last years there has been a lack of statistically significant signals of new physics in both collider and dark matter experiments. These constraints affect all BSM models, and in the case of SUSY they imply that either the new particles are all much heavier than expected, or that the SUSY spectrum is much more complicated than what was initially expected. This situation motivates the relaxation of assumptions typically taken in past works, such as supergravity-inspired spectra, attempting at the same time to keep most of the attractive features of such models.
In particular, we are interested in preserving naturalness as best as possible. This has important consequences, as we know from fine-tuning arguments [1] that the µ parameter should be close to the electroweak breaking scale. This implies that the lightest neutralinos χ 0 1,2 and charginoχ ± 1 should be of higgsino type, which leads to them being nearly mass degenerate [2].
Another feature that we wish to keep is having a good DM candidate. On the MSSM, the DM candidates are the lightest sneutrinoν L and neutralinoχ 0 , which are so-called weakly interacting massive particles (WIMPs). Theν L has no difficulty in reproducing the correct relic density for GeV-scale masses, but in this case the spin-independent DM-nucleon crosssection lies above the current bounds, such as XENON1T [3]. On the other hand, if theχ 0 is of Higgsino type, one finds that for O(100 GeV)-scale masses the correct relic density cannot be reproduced [4]. Given this situation, if we insist on having a natural SUSY solution to the DM problem, it is necessary to expand the model.
Having theν R as an LSP is also of phenomenological interest in the context of colliders [17,19,[23][24][25][26][27][28]. In ref. [26] it was found that there exists a region of the parameter space that has not been probed so far by the LHC, characterized by having aν R LSP, with light sleptons and higgsinos, as well as heavy gauginos. This region was also characterized by having large Y ν couplings, which are interesting for two reasons. First, having large Yukawa couplings allows us to avoid collider bounds on long-lived sleptons. And second, this case could be probed in future experiments, such as those searching for charged lepton flavour violation, lepton number violation, or heavy neutral leptons at colliders. It is thus of our interest to probe further this region of the parameter space, and find out what conditions do the ν R andν R have to follow in order to reproduce the correct relic density. This is the main motivation behind this paper.
To this end, in Sec. II, we present the main features of a minimal model where only theν R superfields are added. We study the conditions where thermal equilibrium can be attained, and explore thermal freeze-out as well as super-WIMP and freeze-in mechanisms.
In Sec. III we add a U (1) B−L gauge symmetry, which is spontaneously broken by additional superfields. We find the new superfields in this model can act as mediators which help thẽ ν R attain thermal equilibrium.

A. Model definition
As mentioned previously, we add three sterile neutrino superfieldsν Rk (k = 4, 5, 6) to the MSSM particle content, and assume conserved R-parity. With this, the superpotential reads as The corresponding soft SUSY breaking terms are given by For only one family ofν L /ν R , the seesaw mechanism determines the size of the Yukawa couplings in terms of the neutrino masses, Y ν ∼ ( √ 2m M )/v u . If the heavy neutrino masses M are on the GeV scale, one expects very small Yukawa couplings. For example, assuming a light neutrino mass m 1 = 10 −3 eV and tan β = 6, then, for M = 1 GeV we have Y ν = 6×10 −9 .
Nevertheless, when having more than one generation ofν R , it is possible to enhance the Yukawas. In this case, it is useful to implement a Casa-Ibarra-like parametrization [29,30].
This describes neutrino mixing in terms of the active-light mixing matrix U PMNS , all six neutrino masses, and the orthogonal R matrix: where s ij and c ij are the sines and cosines of new complex angles, ρ ij + iγ ij . The γ ij lead to hyperbolic functions, which exponentially enhance the Yukawas.
If only one of the γ ij is non-zero, the Yukawas can be expressed in a simple way. For example, by taking only γ 56 different from zero, and normal mass ordering, we find: where a = e, µ, τ , and m 1 , m 2 , m 3 (M 4 , M 5 , M 6 ) are the light (heavy) neutrino masses.
Notice that only the couplings ofν R5 andν R6 are enhanced, withν R4 following the standard Seesaw expectation. By taking γ 56 ≈ 3, 5.3, 7.6, 9.9, the elements Y a5 and Y a6 are enhanced by a factor 10, 10 2 , 10 3 and 10 4 , respectively. Switching on the other angles does not change the generic feature that the couplings of two of the heavy neutrinos are enhanced with respect to the third one. This is a consequence of the fact that one needs an even number of right-handed neutrinos forming pseudo-Dirac neutrinos. In the following, we shall use this scenario, setting all other γ ij and ρ ij to zero.
It is well known that when enhancing the heavy neutrino couplings one can have large contributions to neutrinoless double beta decays from ν R5 and ν R6 . The non-observation of this process forces the latter to have almost degenerate masses, M 5 ≈ M 6 [31][32][33][34], that is, they form a pseudo-Dirac pair. This statement holds in the presence of R-parity conserving SUSY, as no tree-level SUSY contributions to this process exist.
Let us now turn to theν sector. As was discussed in [26], if we neglect Bν and Y ν M † R terms, and assume vanishing CP-violation in the sneutrino mass matrix, the real and imaginary parts of the sneutrino fields remain aligned. This means that we can work directly with complex (ν L ,ν R ), with the mass matrix having the following leading terms: Diagonalization of this mass matrix leads to the lightest mass eigenstateν 1 , which in our framework is the LSP. Theν 1 state will be a superposition ofν R andν L . Thus, for one generation, we can haveν 1 = cosθν R + sinθν L with L-R mixing angle: (2.9) Here we have defined T ν ≡ Y ν A ν . The mixing is strongly suppressed by Y ν and, given our assumptions on the size of µ, can only be enhanced by taking a very large A ν , or tuning the masses such that the denominator vanishes. In the following we shall not consider any of these possibilities, such that all L-R mixing effectively vanishes.
In this work, for definiteness, we take the neutrino and SUSY parameters as shown in Tab. I, with all CP phases equal to zero. Oscillation parameters can be found in [35,36].
M 4 will be allowed to vary between 1 keV and 1 MeV within our results. For the other heavy neutrino masses given in Tab. I, direct search bounds [37] restrict γ 56 8. Taking as example γ 56 = 3 (7), we find |(Y ν ) a5 | = |(Y ν ) a6 | ∼ 10 −7 (10 −5 ). On the sneutrino soft sector, only mL and mν R are non-zero. Given the strong flavour constraints coming from processes such as µ → eγ, we take the soft masses flavour diagonal. This implies that all mixing effects are negligible, allowing us to identify the mass eigenstates with the interaction eigenstates.
In particular,ν 1 will be one of the threeν R states.

B. Neutrino Dark matter
It is well known that a sterile neutrino with mass in the keV range can serve as potential DM candidate via oscillations with the left-handed neutrinos [38,39]. However, it is unstable, as it can decay into the lighter active neutrinos through the following processes: The first final state can be detected via the resulting photon, e.g. by cosmological observations. In particular, data from the NuStar collaboration puts severe constraints on the allowed parameter space, as has been shown recently in [40]. From their Fig. 5, considering only the usual W -l-loop mediated contribution, one gets an upper bound for the quantity of few times 10 −11 for M 4 = 7 keV and up to around 10 −14 for M 4 = 50 keV. We find that the sin 2 2θ predicted in our scenario is always well above the corresponding limit, as can be seen from the values given in Tab. II (see also [41] for a recent discussion on the Seesaw prediction).
We note that it has been assumed in [40] that the decay into a photon occurs with probability 1. However, when taking the three body decays into three active neutrinos via a virtual Z-boson into account, we find it actually has a branching ratio of about 0.1, which weakens the bound, e.g. somewhat larger values of sin 2 2θ are allowed. Nevertheless, the bound would need to be relaxed by at least four orders of magnitude for a tiny region around M 4 = 7 keV to be marginally allowed. This would require a careful re-analysis of the investigation carried out in [40] which is beyond the scope of this paper.
Still, one might wonder if the branching ratio into the photon could be modified via loops containing SUSY particles 2 , such as those shown in Fig. 1. We have varied the soft SUSY parameters in the ranges 450 , as well as 100 ≤ |µ| ≤ 500 GeV. We find at most a 10 per-cent variation in the branching ratio for scenarios with light sleptons and charginos, and small tan β. The main reason for this is due to the experimental bounds on the masses of the SUSY particles, which imply that the corresponding loops are suppressed with respect to the W -l contributions. Thus, we conclude that the presence of SUSY does not affect the NuStar bounds.
In light of this negative result, one has two options to avoid the bounds. The first option is to increase M 4 above the MeV scale, such that ν 4 decays quickly and leaves no signal. The second option is to decrease m 1 under 10 −9 eV, which suppresses the ν 4 Yukawa coupling (see Eq. (2.4)). This in turn decreases the mixing with the active neutrinos. Nevertheless, in both options ν 4 can no longer be a good dark matter candidate, in the first case it decays too soon, while on the second case it is unlikely to be produced via oscillations with the left-handed neutrinos [38,39]. We thus discard ν 4 as a candidate for dark matter in this model.

C. Sneutrino Dark matter
The calculation of theν 1 relic density depends on whether it is in thermal equilibrium with the primordeal plasma, or not. For this to happen, one requires that at some temperature T we have: where H(T ) is the Hubble constant, σv T is the thermal average of the cross-section times velocity, and n(T ) is the number density of theν 1 . Ultimately, a relevant factor in this condition is the size its couplings with other thermal particles, which are in turn determined by the size of the Y ν .
If theν 1 couplings are large enough such that it can be in thermal equilibrium at some temperature, the final value of Ωh 2 follows the freeze-out mechanism, see e.g. [43] and refs. therein. Here, the expansion of the Universe leads to a point where Eq. (2.13) does not hold, and theν 1 decouples from the plasma. The relic density depends on the number density at the freeze-out temperature, as well as on (co-)annihilations with other thermal particles.
If the couplings are too small, theν 1 relic density proceeds from decays or annihilation of other particles in the primordeal plasma. In the former case, the value of the relic density is mainly due to decays of the next to lightest SUSY particle (NLSP). For example, thermal neutralinos could decay after they freeze-out viaχ 0 →ν 1 ν. In this situation theν 1 is said to follow the super-WIMP mechanism [12,[44][45][46], and its relic density is proportional to the NLSP yield after freeze-out. In contrast, if the annihilations are more important than decays, the mechanism is that of freeze-in [47]. Here, theν 1 is created throughout the thermal history of the Universe by annihilations of all thermal particles, not only the NLSP.
In both of the scenarios described above, theν 1 is called a feebly interacting massive particle (FIMP), and one assumes that it was not generated by other means at earlier times, such as during an inflationary period.
The naive seesaw expectation is that a lightν 1 should have very small Y ν couplings, and should thus be non-thermal. An important result of this procedure is that the NLSP can be long lived, which can lead to displaced vertices, disappearing tracks or heavy metastable charged particle signals at colliders [19,23,48].
If the neutrino mass generation mechanism allows for larger Yukawas, theν R can be thermal [13]. From Eq. (2.13), we see that as long as σv T n(T )/H(T ) > 1 for some value of T , these particles will have been in thermal equilibrium at some point in the history of the Universe. Thus, as a first step, we evaluate the required size of Y ν for this condition to hold.
To this end, we used SARAH 4.14.0 [49][50][51][52][53] to implement the model of [26] in SPheno 4.0.3 [54,55], which calculates the mass spectrum and branching ratios. We used SSP 1.2.5 [56] to carry out the parameter variation. With the output, the variable σv T n(T )/H(T ) was calculated with Micromegas 5.0.9 [57]. Here one has to keep two details in mind. First of all, the program only provides σv ann T , the (co-)annihilation ofν R with themselves and other SUSY particles into SM final states. In addition, Micromegas provides the total n eq (T ) from all SUSY particles, instead of the exclusive one forν 1 . We do not consider this a serious problem, as n eq (T ) is dominated by the LSP. Still, in order to be certain of our results, in the following we shall consider σv ann T n eq (T )/H(T ) < 0.1 as non-thermal, and σv ann T n eq (T )/H(T ) > 10 as thermal. Our results are shown in Fig. 2, where we evaluate σv ann T n eq (T )/H(T ) as a function of T . Curves are shown for different SUSY masses, as well as different Yukawa couplings. We on the left and right panels, respectively. We fix the SUSY masses as in Tab. I, and probeν R soft masses of 100, 200 and 300 GeV (blue, green and orange, respectively). We also include MSSM slepton soft masses between 150 and 450 GeV, with the restriction mL = mẼ > mν R , such that the R-sneutrino is always the LSP.
In the Figure, we find that all curves of σv ann T n eq (T )/H(T ) increase very quickly with T up to a peak, with a subsequent drop as temperatures increase. This decrease is not surprising, since H ∼ T 2 , n eq ∼ T 3 and σv T is expected to fall as 1/T 2 for large tem- , where we find no temperature which satisfies the condition σv ann T n eq (T )/H(T ) > 1 and theν R is certainly a FIMP. In the following, we explore the status of DM in the setup of ref. [26]. As mentioned in the Introduction, we seek to enhance Yukawa couplings in order to avoid bounds on longlived particles, and to hope for experimental correlations between the ν R andν R sectors. To achieve this, we set γ 56 = 7, soν R5 andν R6 are thermal particles. For simplicity, we also fix mν R5 = mν R6 . In contrast, the remaining eigenstateν R4 is left with small Yukawa couplings, such that it is a non-thermalν R .
If theν R4 mass is larger than that ofν R5 orν R6 , then the latter is a thermal LSP, and the relic density (Ωh 2 ) th is generated via freeze-out. In this case we find that the value of (Ωh 2 ) th calculated by Micromegas is extremely large, ruling out the whole parameter space.
Alternatively, ifν R4 is the LSP, the relic density is obtained from a combination of super-WIMP (Ωh 2 ) dec and freeze-in (Ωh 2 ) in processes 3 . The super-WIMP contribution is obtained from: where m th N LSP is the mass of the lightest thermal SUSY partner. Thus, one can adjust this contribution to any desired value by choosing an appropriate mν 1 . In this case it is not necessary for the thermalν R to be NLSP, it is perfectly possible to have a thermalν L orh NLSP which later decays to theν 1 .
The super-WIMP contribution by itself cannot explain the observed relic density in two situations. First, since (Ωh 2 ) dec ≤ (Ωh 2 ) th , it is impossible to get the correct relic density if (Ωh 2 ) th is too small (this can be the case for bothν L andh NLSP). Second, it is also possible for (Ωh 2 ) th to be so large that the feebleν 1 mass must be lower than the ν 4 mass.
This situation would require a negative soft m 2 ν R4 , which we shall discard 4 . Nevertheless, only the latter case is a problem, since a small (Ωh 2 ) dec can be complemented by a larger (Ωh 2 ) in .
We show the results of this super-WIMP-only analysis in Fig. 3. We vary mL = mẼ and mν R ≡ mν R5,6 , which are the soft masses of thermal particles. The black contour lines shown the maximum allowed mass ofν 1 , such that the super-WIMP mechanism saturates the relic density, that is, (Ωh 2 ) dec = 0.12. We find that mν 1 spans from 1 keV to almost 300 GeV.
On the gray regions of the Figure, the thermal L-sleptons and higgsinos have large (co-)annihilation cross-sections, leading to a too small (Ωh 2 ) th . This means that the super-WIMP mechanism cannot explain by itself the observed relic density, so there is no upper bound on mν 1 apart from the requirement of being the LSP.
The brown and red regions are excluded. On the brown region, the feebleν 1 requires a mass lower than the minimum of M 4 = 1 keV, so (Ωh 2 ) dec is too big. This region becomes larger if one decides to use a heavier M 4 to avoid the NuStar constraint mentioned in Sec. II B. Moreover, the red region was ruled out by LHC searches in Ref. [26].
Due toν 1 being a FIMP, it is important to be aware of the thermal NLSP lifetime. If longer than one second, it can be subject to Big Bang Nucleosynthesis (BBN) constraints [19,[59][60][61]. These vary depending on the lifetime and decay modes on the NLSP, and are left for a future study. However, for completeness, in It is important to take into account that this information on the lifetime is strictly valid for the selected values of m 1 and M 4 , since the |(Y ν ) a4 | 2 coupling ofν R4 depend on them (see Eq. (2.4)). For smaller m 1 or M 4 , the lifetime will increase.
We now turn to the calculation of the freeze-in contribution, (Ωh 2 ) fi , which is again calculated with Micromegas. Throughout the evaluated parameter space, we find that if we take mν 1 equal to the maximum value allowed by the super-WIMP mechanism, the freeze-in contribution can greatly exceed the observed value.
One way of suppressing freeze-in is by reducing the lightest neutrino mass m 1 , which decreases (Y ν ) a4 . Due to the long calculation time, we show results for only one representative point in Fig. 4, with M 4 = 1 MeV, mL = 323 GeV and mν R = 302 GeV, whereν R5 is the NLSP. We see that reducing (Y ν ) a4 increases theν R5 lifetime, as shown in the red curve of the Figure. Given the large lifetime, assuming that BBN constraints are avoided, it is 5 For this choice of parameters, the brown region should be extended to the 1 MeV contour. still important to check that the time for recombination is not exceeded, else theν R5 would decay intoν 1 after the cosmic microwave background was emitted. In this case, the relic density observed by Planck [62] would correspond to (Ωh 2 ) th , ruling out the scenario. For our example, we find that (Ωh 2 ) fi becomes subdominant when m 1 10 −9 eV, and that the recombination constraint is satisfied. This solution is also attractive as no additional ν R DM is produced, since both M 4 is large and (Y ν ) a4 is tiny.
Alternatively, it is possible to have a smaller mν 1 , such that the super-WIMP contribution is negligible. In this case, one needs to adjust the lightest neutrino mass, such that (Ωh 2 ) fi reproduces the observations. Such a scenario is shown on the yellow curve of Fig. 4 In principle one could go to even lower mν 1 and M 4 values, such as in the keV range.
However, for these points one needs large (Y ν ) a4 values to get the correct Ωh 2 . This means that it is likely to have a ν 4 DM component, which could be in potential conflict with the NuStar bounds.
Thus, we can conclude that in this scenario, where (Y ν ) a4 is suppressed and (Y ν ) a5 , (Y ν ) a6 are enhanced, the lightest R-sneutrino can be a good non-thermal dark matter candidate, either by super-WIMP or freeze-in mechanisms, with BBN and the time for recombination as constraints.
As a final remark, we must comment that other works attempted to have a thermalν R by generating a large mixing with theν L [14,16,17,27]. The main idea was to allow the mass eigenstate to interact via gauge interactions, with a relative suppression of the overall coupling due to the mixing. This would allow the sneutrino to remain thermal, but at the same time avoid direct detection bounds. For this approach to work, it was required to have sinθ 10 −2 .
As we see from Eq. (2.9), the L-R mixing is suppressed by Y ν . It is clear that even if we increased γ ij to their maximum values allowed by experiment, an additional enhancement is required to bring the mixing to the required level. Earlier works have considered an extremely large A ν , such that Y ν A ν becomes of the order of the other soft terms [63], e.g. the slepton mass parameters. However, this creates the danger of a charge breaking minimum similar to the well-known problem of charge and color breaking minima within the MSSM, see e.g. [64] for an extensive discussion on tree-level constraints.
In our case, we need to assess if the addition of the newν R can lead scalar fields other than the Higgs to acquire a vacuum expectation value (vev ). In particular, we will evaluate the following vev pattern for one generation of sleptons and sneutrinos: 15) as this corresponds to a D-flat direction in the scalar potential. Assuming that all other fields do not acquire a vev and that all parameters are real, one gets for the potential at tree level In order to avoid this direction in the potential from being equal or lower than zero, a sufficient way is requiring a negative discriminant in the solutions for α of the equation V (α) = 0. We get the following bound 6 which implies it is A ν , and not T ν , who can be at most of the order of soft terms. Thus L-R mixing is not a safe procedure to solve the DM problem in this context. We note for completeness that this tree-level estimate can get significant loop corrections, as has been in shown in similar models [66,67]. However, the estimate gives the correct order of magnitude, which is sufficient in our case. We have explicitly checked that this estimate is correct within a factor two for a few points using the package Vevacious [68].
We have seen in the previous Section that in the minimal model it is not possible to have a thermalν R as an LSP, as it will yield a too large contribution to Ωh 2 . We therefore seek an extension where the R-sneutrino can be a thermal relic, potentially giving the correct relic density, without impacting too much the high-energy and precision phenomenology, in particular in view of the collider constraints obtained in [26]. Several possibilities have been considered in a similar context in the literature: an additional U (1) gauge factor [69][70][71][72], left-right symmetric models [73], the NMSSM [74,75] or via additional F-and D-terms as they occur in models of hybrid inflation [76].
Taking a U (1) B−L extension as an example, we will show how this can be achieved in the parameter space we are interested in here. For this, we first briefly summarize the main aspects of the model presented in [77]. Here, the MSSM particle content is extended by three new types of superfields. First, one has a B vector superfield associated to the U (1) B−L symmetry. Second, one adds two new Higgs-like chiral superfields,η andη, carrying B − L number ±2, whose scalar components break U (1) B−L and provide mass to the Z boson.
These chiral superfields will be called bileptons in the following. Note, that the symmetry breaking is such that R-parity is conserved 7 . Finally, the anomaly-cancellation requires the existence of three right-handed neutrino superfieldsν R [78]. These have masses around 10-100 GeV if the U (1) B−L breaking scale is of the order 1-100 TeV.
The superpotential is given by and one has the additional soft SUSY-breaking terms: Without loss of generality one can take B µ and B µ to be real, as in the MSSM. The extended gauge group breaks to SU (3) C ⊗ U (1) em as beside the Higgs fields also the bileptons receive vev s, denoted by v η and vη, see ref. [77] for details. We define tan β = vη vη in analogy to the ratio of the MSSM vev s. In the following we will neglect all effects concerning gauge kinetic mixing as they do not impact our findings, e.g. we will set M BB = 0.
As in the MSSM one can always get a SM-like light Higgs boson h, with mass of 125 GeV, assuming third generation squarks sufficiently heavy and having a sizable mixing in the stop sector. Beside h this model contains a second scalar h , which is light provided tan β is close to 1 [77]. As the limits on Z are in the multi-TeV range, this state is mainly an SM gauge-singlet, since the mixing scales like v u / v 2 η + v 2 η . However, since h can still have a non-negligible mixing with h , we use HiggsBounds [79] to check if its properties are compatible with the experimental results.
Neglecting mixing effects between the MSSM Higgs bosons and the bileptons, one gets for the mass of h at tree level [77]: Dark matter aspects of this model have been discussed in [72] where various candidates were evaluated, including sneutrinos. It has been shown that in case of a CP-even sneutrino it is possible to reproduce the correct relic density through a Higgs funnel mechanism, similar to other related models (see discussion in [18]). This mechanism is also applicable to CP-odd sneutrinos, being a consequence of the D-term and F -term couplings of the h to theν R .
One issue that we want to address here is whether the Higgs funnel mechanism still explains the observed Ωh 2 if we enforce a small or even negligible mass splitting between the scalar and pseudoscalar sneutrino states.
As it is well know, funnel regions imply particular relation between various masses, e.g.
here we shall require mνS 1 m h /2. Moreover, having a small splitting between the scalar and pseudoscalar implies A x sin β µ cos β as can be seen from Eqs. (3.4) and (3.5). We take as input the mass of the pseudoscalar sneutrino mνP 1 , tan β , m Z and fix the mass difference betweenν S 1 andν P 1 to 0.5 GeV. The values of the parameters, which we fix for these investigations, are given in Tab. III.
As an example we show Ωh 2 as a function of tan β in Fig. 5. In the Figure,  region, excluded due to tachyonic states, this channel is also closed and the dominant process isνν → ν 4 ν 4 . Since it is suppressed by a small Y x , the relic density increases considerably.
Thus, similar to the case of the MSSM, the Higgs funnel is only possible within a small strip in parameter space. In the case shown we find exactly two working values of tan β .
Interestingly enough, both are consistent with the experimental constraints, which we have checked using HiggsBounds. We note for completeness, that we find the soft parameter m 2 ν R < 0 when tan β < 1. However, when matching this model to the one of the previous Section, the D-term contribution from the Z yields a positive value. 9 In the points explaining the relic density, h has about one per-cent admixture with the SM-like Higgs boson.
The mixing among the Higgs bosons actually imply that the SM-like Higgs boson also yields a funnel region. This is shown in Fig. 6 where we vary the sneutrino mass parameters but keep the Higgs masses as at the points in Fig. 5 yielding the correct relic density. On the right (left) plot we display the case for tan β = 0.972 (tan β = 1.032). For these values, we actually have four (two) possible sneutrino masses for the smaller (larger) value of tan β which yield the correct relic density. These correspond to scenarios where either the gauge-singlet h or the SM-like h give rise to the Higgs funnel.
Last but not least we note for completeness that the Z leads to a thermal distribution of the right-handed neutrinos in the early universe. However, in the preferred keV mass range for explaining the correct relic density, it is actually a warm dark matter candidate, which is excluded by observations of the Lyman-α forest [80][81][82]. To avoid this, as before, one needs to set a relatively large M 4 , such that the ν 4 decay quickly. Thus, the right-handed neutrino cannot contribute to the DM in this model.

IV. CONCLUSIONS
In the last years there has been increasing interest in scenarios where a light right-handed neutrino could explain the observed relic dark matter density. This has also been motivated by the fact that no hint of a WIMP dark matter particle has been found so far. At the same time there has also been increasing interest in so-called natural SUSY scenarios where the higgsinos are the LSP. While in these scenarios one can avoid fine-tuning in the Higgs sector, one cannot explain the observed relic density as the higgsinos annihilate very effectively.
We have thus combined both ideas and considered in a first step a supersymmetric scenario where the MSSM is extended by right-handed neutrino superfields where the Yukawa couplings are enhanced by an inverse seesaw structure. Here we attempted to address several questions: (i) Is it possible for SUSY loops to modify the lifetime of the lightest right-handed neutrino ν 4 ? (ii) Is it possible for an R-sneutrinoν R to be a thermal relic in such scenarios?
(iii) If non-thermal, under which conditions would dark matter production occur?
We have found that the answer to the first two questions is negative in the more minimal model. (i) The LHC bounds on the SUSY particles are already so strong that the corresponding contributions to the decay ν 4 → γν i are suppressed, leading to a change of the partial widths of at most 10 per-cent with respect to the SM contribution. Consequently, since the ν 4 should have a mass of a few keV to be an acceptable DM candidate, we find this scenario to be heavily constrained by cosmological data. (ii) In scenarios whereν R are the LSP the Yukawa couplings can be at most of the order 10 −5 . This size is sufficient for theν R to thermalize, but the corresponding relic density is too large by many orders of magnitude. The only exception are regions where a sizable mixing with theν L is present, as has already been noticed in the literature. However, we have seen that these regions lead to charge breaking minima and, thus, are excluded.
We then find that the only possibility for R-sneutrino dark matter, in the minimal model, is (iii) having a non-thermal FIMP, produced via super-WIMP or freeze-in mechanisms. We find that the observed relic density can be reproduced through both mechanisms, with the super-WIMP depending on the ratio ofν 1 and NLSP masses. The freeze-in is in addition sensitive to (Y ν ) a4 . By modifying these parameters, one can make either super-WIMP or freeze-in dominant. We have also found that the NLSP lifetime can be significantly large, which could be in conflict with BBN.
Although theν R cannot be a thermal DM candidate in the MSSM+ν R framework, as discussed above, it might well be that this is only an effective model, e.g. an additional gauge group might be realized in the multi-TeV range. As an example we have considered the case of an additional U (1) B−L and extended previous studies by the mass hierarchy studied in the MSSM+ν R . The additional Z in the multi-TeV range implies that the ν R andν R get thermalized in the early Universe. In this class of models a second light Higgs boson is possible, which mixes somewhat with the SM-like Higgs boson. We have shown this allows for an explanation of the observed relic density via a Higgs funnel.