Impact of leptonic unitarity and dark matter direct detection experiments on the NMSSM with inverse seesaw mechanism

TL;DR

This paper investigates how leptonic unitarity violations and dark matter detection experiments constrain the NMSSM with inverse seesaw, revealing strict bounds on neutrino Yukawa couplings and the potential impact of future experiments.

Contribution

It provides a comprehensive analysis of experimental constraints on the NMSSM with inverse seesaw, highlighting the complementarity and strictness of these bounds.

Findings

Leptonic unitarity and XENON-1T data limit neutrino Yukawa couplings.

Constraints are complementary and sometimes very strict.

Future LZ experiment could significantly impact the model.

Abstract

In the Next-to-Minimal Supersymmetric Standard Model with the inverse seesaw mechanism to generate neutrino masses, the lightest sneutrino may act as a feasible dark matter candidate in vast parameter space. In this case, the smallness of the leptonic unitarity violation and the recent XENON-1T experiment can limit the dark matter physics. In particular, they set upper bounds of the neutrino Yukawa couplings $\lambda_\nu$ and $Y_\nu$. We study such effects by encoding the constraints in a likelihood function and carrying out elaborated scans over the parameter space of the theory with the Nested Sampling algorithm. We show that these constraints are complementary to each other in limiting the theory, and in some cases, they are very strict. We also study the impact of the future LZ experiment on the theory.

Tables

Table 1. Table 1 : Specific configuration of the Higgs and neutralino sectors for the scenarios discussed in the text, and their prediction on the properties of the Higgs bosons and neutralinos such as the mass spectrum and the couplings of the Higgs bosons with different particles, C ¯ i ​ j ​ k subscript ¯ 𝐶 𝑖 𝑗 𝑘 \bar{C}_{ijk} , which are normalized to their corresponding SM predictions. Parameters with mass dimensions are in the unit of GeV. Other fixed parameters that are not listed in the table include m q ~ = 2000 ​ GeV subscript 𝑚 ~ 𝑞 2000 GeV m_{\tilde{q}}=2000~{}{\rm GeV} for flavor universal soft-breaking masses of squarks, M 1 = M 2 = 2000 ​ GeV subscript 𝑀 1 subscript 𝑀 2 2000 GeV M_{1}=M_{2}=2000~{}{\rm GeV} and M 3 = 5000 ​ GeV subscript 𝑀 3 5000 GeV M_{3}=5000~{}{\rm GeV} for gaugino masses, A i = 0 subscript 𝐴 𝑖 0 A_{i}=0 for all soft-breaking trilinear coefficients except for A λ subscript 𝐴 𝜆 A_{\lambda} and A t subscript 𝐴 𝑡 A_{t} , and [ Y ν ] 11 , 22 = 0.01 subscript delimited-[] subscript 𝑌 𝜈 11 22 0.01 [Y_{\nu}]_{11,22}=0.01 , [ λ ν ] 11 , 22 = 0.3 subscript delimited-[] subscript 𝜆 𝜈 11 22 0.3 [\lambda_{\nu}]_{11,22}=0.3 , and [ m ν ~ ] 11 , 22 = [ m x ~ ] 11 , 22 = 2000 ​ GeV subscript delimited-[] subscript 𝑚 ~ 𝜈 11 22 subscript delimited-[] subscript 𝑚 ~ 𝑥 11 22 2000 GeV [m_{\tilde{\nu}}]_{11,22}=[m_{\tilde{x}}]_{11,22}=2000~{}{\rm GeV} for the parameters of the first two generations of the sneutrinos. All these parameters are defined at the scale Q = 1000 ​ GeV 𝑄 1000 GeV Q=1000~{}{\rm GeV} . Besides, the Higgs masses and U i ​ j subscript 𝑈 𝑖 𝑗 U_{ij} are obtained by setting [ Y ν ] 33 = [ λ ν ] 33 = 0 subscript delimited-[] subscript 𝑌 𝜈 33 subscript delimited-[] subscript 𝜆 𝜈 33 0 [Y_{\nu}]_{33}=[\lambda_{\nu}]_{33}=0 , and sneutrino loop effects may slightly change them when varying the parameters in the sneutrino DM sector.

Light $h_s$ scenario with $A_{\lambda }=2000\mathrm{GeV}$						Massive $h_s$ scenario with $A_{\lambda }=2000\mathrm{GeV}$
$\tan \beta$	12.38	$\lambda$	0.24	$\kappa$	0.23	$\tan \beta$	28.46	$\lambda$	0.19	$\kappa$	0.60
$A_t$	2433	$A_{\kappa }$	-680.4	$\mu$	195.2	$A_t$	2363	$A_{\kappa }$	-120.4	$\mu$	328.0
$m_{h_s}$	120.4	$m_h$	125.1	$m_H$	2332	$m_h$	125.1	$m_{h_s}$	2042	$m_H$	5381
$m_{A_s}$	608.5	$m_{A_H}$	2331	$m_{H^{\pm }}$	2332	$m_{A_s}$	592.6	$m_{A_H}$	5381	$m_{H^{\pm }}$	5379
$m_{{\tilde{\chi }}_1^0}$	186.0	$m_{{\tilde{\chi }}_2^0}$	-206.1	$m_{{\tilde{\chi }}_1^{\pm }}$	197.9	$m_{{\tilde{\chi }}_1^0}$	318.8	$m_{{\tilde{\chi }}_2^0}$	-341.3	$m_{{\tilde{\chi }}_1^{\pm }}$	333.8
$U_{11}$	-0.018	$U_{12}$	-0.261	$U_{13}$	0.965	$U_{11}$	0.00001	$U_{12}$	-0.999	$U_{13}$	0.003
$U_{21}$	-0.005	$U_{22}$	0.965	$U_{23}$	0.260	$U_{21}$	0.006	$U_{22}$	-0.003	$U_{23}$	-0.999
$U_{31}$	0.999	$U_{32}$	0.000	$U_{33}$	0.019	$U_{31}$	-0.999	$U_{32}$	-0.000	$U_{33}$	-0.006
${\overline{C}}_{h_s{V}^{\ast }V}$	-0.261	${\overline{C}}_{h_s\overline{u}u}$	-0.263	${\overline{C}}_{h_s\overline{d}d}$	-0.038	${\overline{C}}_{h_s{V}^{\ast }V}$	-0.003	${\overline{C}}_{h_s\overline{u}u}$	-0.003	${\overline{C}}_{h_s\overline{d}d}$	-0.162
${\overline{C}}_{h{V}^{\ast }V}$	0.965	${\overline{C}}_{h\overline{u}u}$	0.965	${\overline{C}}_{h\overline{d}d}$	1.03	${\overline{C}}_{h{V}^{\ast }V}$	-0.999	${\overline{C}}_{h\overline{u}u}$	-0.999	${\overline{C}}_{h\overline{d}d}$	-1.00
${\overline{C}}_{H{V}^{\ast }V}$	0.000	${\overline{C}}_{H\overline{u}u}$	0.081	${\overline{C}}_{H\overline{d}d}$	-12.4	${\overline{C}}_{H{V}^{\ast }V}$	-0.000	${\overline{C}}_{H\overline{u}u}$	-0.035	${\overline{C}}_{H\overline{d}d}$	28.5
${\overline{C}}_{h_s{h}_s{h}_s}$	4.014	${\overline{C}}_{h_s{h}_{s}h}$	1.193	${\overline{C}}_{h_{s}hh}$	-0.460	${\overline{C}}_{h_s{h}_s{h}_s}$	-163.4	${\overline{C}}_{h_s{h}_{s}h}$	0.448	${\overline{C}}_{h_{s}hh}$	-1.591
$N_{13}$	-0.696	$N_{14}$	0.681	$N_{15}$	-0.147	$N_{13}$	0.675	$N_{14}$	-0.641	$N_{15}$	0.012
$N_{23}$	-0.699	$N_{24}$	-0.709	$N_{25}$	-0.055	$N_{23}$	-0.702	$N_{24}$	-0.709	$N_{25}$	-0.010

Table 2. Table 2 : Detailed information of the points in the light h s subscript ℎ 𝑠 h_{s} scenario (left side of the table) and the massive h s subscript ℎ 𝑠 h_{s} scenario (right side of the table) with the setting B μ X = 0 subscript 𝐵 subscript 𝜇 𝑋 0 B_{\mu_{X}}=0 . The number before each annihilation process represents the fraction of its contribution to the total DM annihilation cross section at the freeze-out temperature. Parameters in mass dimension are in unit of GeV, and the DM-nucleon scattering cross section are in unit of cm 2 superscript cm 2 {\rm cm^{2}} .

Light $h_s$ scenario					Massive $h_s$ scenario
	$P_1$		$P_2$			$P_3$		$P_4$
$Y_{\nu }$	0.060		0.035		$Y_{\nu }$	0.178		0.227
${\lambda }_{\nu }$	0.384		0.411		${\lambda }_{\nu }$	0.473		0.364
$A_{Y_{\nu }}$	731.6		-722.1		$A_{Y_{\nu }}$	297.8		225.2
$A_{{\lambda }_{\nu }}$	463.3		396.3		$A_{{\lambda }_{\nu }}$	-288.8		-295.0
$m_{\tilde{\nu }}$	495.7		348.4		$m_{\tilde{\nu }}$	191.8		375.9
$m_{\tilde{x}}$	217.6		333.1		$m_{\tilde{x}}$	269.6		442.7
$m_{\tilde{l}}$	857.2		816.5		$m_{\tilde{l}}$	875.3		568.9
$m_{{\tilde{\nu }}_1}$	132.6		179.4		$m_{{\tilde{\nu }}_1}$	313.0		299.6
$V_{11}$	-0.002		-0.007		$V_{11}$	-0.014		-0.031
$V_{12}$	0.526		-0.697		$V_{12}$	-0.717		-0.726
$V_{13}$	-0.851		0.717		$V_{13}$	0.697		0.687
$\Omega h^2$	0.108		0.120		$\Omega h^2$	0.119		0.126
$C_{{\tilde{\nu }}_1{\tilde{\nu }}_1^{\ast }{h}_s}$	34.80		30.95		$C_{{\tilde{\nu }}_1{\tilde{\nu }}_1^{\ast }{h}_s}$	-53.54		-134.9
$C_{{\tilde{\nu }}_1{\tilde{\nu }}_1^{\ast }h}$	8.097		7.481		$C_{{\tilde{\nu }}_1{\tilde{\nu }}_1^{\ast }h}$	3.694		4.927
${\sigma }_{{\tilde{\nu }}_1-p}^{SI}$	$3.64\times {10}^{-48}$		$1.90\times {10}^{-48}$		${\sigma }_{{\tilde{\nu }}_1-p}^{SI}$	$2.28\times {10}^{-48}$		$4.91\times {10}^{-47}$
${\sigma }_{{\tilde{\nu }}_1-n}^{SI}$	$4.34\times {10}^{-48}$		$8.00\times {10}^{-48}$		${\sigma }_{{\tilde{\nu }}_1-n}^{SI}$	$3.94\times {10}^{-46}$		$7.64\times {10}^{-45}$
			$16\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^0\to W^{+}{W}^{-}$		$11\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to d\overline{u}$	$42\%$	${\tilde{\nu }}_1{\tilde{\chi }}_1^0\to W^{+}{\tau }^{-}$
			$13\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to d\overline{u}$		$11\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to s\overline{c}$	$21\%$	${\tilde{\nu }}_1{\tilde{\chi }}_1^0\to Z{\nu }_{\tau }$
annihilation	$89\%$	${\tilde{\nu }}_1{\tilde{\nu }}_1^{\ast }\to h_s{h}_s$	$13\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to s\overline{c}$	annihilation	$11\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^0\to W^{+}{W}^{-}$	$19\%$	${\tilde{\nu }}_1{\tilde{\chi }}_1^0\to h{\nu }_{\tau }$
processes	$10\%$	${\tilde{\nu }}_1{\tilde{\nu }}_1^{\ast }\to h_{s}h$	$11\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^0\to ZZ$	processes	$8.9\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to b\overline{t}$	$1.6\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to d\overline{u}$
			$8.9\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to b\overline{t}$		$7.4\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^0\to ZZ$	$1.6\%$	${\tilde{\chi }}_1^0{\tilde{\chi }}_1^{-}\to s\overline{c}$
			…	. . .		…	. . .	…	. . .

Equations

$$W$$

$$L^{soft}$$

$${\mathcal{M}}_{11}^2$$

$${\mathcal{M}}_{12}^2$$

$${\mathcal{M}}_{13}^2$$

$${\mathcal{M}}_{22}^2$$

$${\mathcal{M}}_{23}^2$$

$${\mathcal{M}}_{33}^2$$

$$h_i=\sum _{j=1}^3{U}_{ij}{S}_j,$$

$${\bar{C}}_{h_i{V}^{\ast }V}=U_{i2},{\bar{C}}_{h_i\bar{u}u}=U_{i1}\cot \beta +U_{i2},{\bar{C}}_{h_i\bar{d}d}=-U_{i1}\tan \beta +U_{i2}.$$

$${\mathcal{M}}_{P,11}^2$$

$${\mathcal{M}}_{P,22}^2$$

$${\mathcal{M}}_{P,12}^2$$

$${\cal M}=\left(\begin{array}[]{ccccc}M_{1}&0&-\frac{g^{\prime}v_{d}}{\sqrt{2}}&\frac{g^{\prime}v_{u}}{\sqrt{2}}&0\\ &M_{2}&\frac{gv_{d}}{\sqrt{2}}&-\frac{gv_{u}}{\sqrt{2}}&0\\ &&0&-\mu&-\lambda v_{u}\\ &&&0&-\lambda v_{d}\\ &&&&\frac{2\kappa}{\lambda}\mu\end{array}\right),$$

$${\tilde{\chi }}_i^0=N_{i1}{\psi }_1^0+N_{i2}{\psi }_2^0+N_{i3}{\psi }_3^0+N_{i4}{\psi }_4^0+N_{i5}{\psi }_5^0.$$

$$\displaystyle M_{\mathrm{ISS}}=\left(\begin{array}[]{c c c}0&M_{D}^{T}&0\\ M_{D}&0&M_{R}\\ 0&M_{R}^{T}&\mu_{X}\end{array}\right)\,,$$

$$\displaystyle\left(U_{\nu}^{\dagger}\right)_{9\times 9}=\left(\begin{array}[]{cc}\hat{U}_{3\times 3}&X_{3\times 6}\\ Y_{6\times 3}&Z_{6\times 6}\end{array}\right).$$

$$M_{\nu }=\left[M_D^T{M}_R^{T^{-1}}\right]{\mu }_X\left[(M_R^{-1})M_D\right]+\mathcal{O}({\mu }_X^2)\equiv F{\mu }_X{F}^T+\mathcal{O}({\mu }_X^2).$$

$$U_{\mathrm{PMNS}}^T{M}_{\nu }{U}_{\mathrm{PMNS}}=\mathrm{diag}(m_{{\nu }_1},m_{{\nu }_2},m_{{\nu }_3}).$$

$$\hat{U}\simeq \left(1-\frac{1}{2}F{F}^{†}\right)U_{\mathrm{PMNS}}\equiv (1-\eta )U_{\mathrm{PMNS}},$$

$$\sqrt{2|\eta {|}_{ee}}<0.050,\sqrt{2|\eta {|}_{\mu \mu }}<0.021,\sqrt{2|\eta {|}_{\tau \tau }}<0.075,$$

$$\sqrt{2|\eta {|}_{e\mu }}<0.026,\sqrt{2|\eta {|}_{e\tau }}<0.052,\sqrt{2|\eta {|}_{\mu \tau }}<\mathrm{0.035.}$$

$${\mu }_X=M_R^T{m}_D^{T^{-1}}{\hat{U}}_{3\times 3}^{\ast }\mathrm{Diag}(m_{{\nu }_1},m_{{\nu }_2},m_{{\nu }_3}){\hat{U}}_{3\times 3}^{†}{m_D}^{-1}{M}_R,$$

$$\left|\frac{[{\lambda }_{\nu }{]}_{11}\mu }{[Y_{\nu }{]}_{11}\lambda v_u}\right|>14.1,\left|\frac{[{\lambda }_{\nu }{]}_{22}\mu }{[Y_{\nu }{]}_{22}\lambda v_u}\right|>33.7,\left|\frac{[{\lambda }_{\nu }{]}_{33}\mu }{[Y_{\nu }{]}_{33}\lambda v_u}\right|>\mathrm{9.4.}$$

$$\delta {\mu }_X\sim \frac{1}{16{\pi }^2}\frac{{\lambda }_{\nu }{M}^2{B}_{{\mu }_X}{M}^2{\lambda }_{\nu }}{M_{\mathrm{SUSY}}^5},$$

$$\delta {\mu }_X(\mathrm{keV})\sim 63.3\times {\left(\frac{{\lambda }_{\nu }}{0.1}\right)}^2{\left(\frac{{\mathrm{M}}^2}{{\mathrm{M}}_{\mathrm{SUSY}}^2}\right)}^2\left(\frac{{\mathrm{B}}_{{\mu }_{\mathrm{X}}}/Ge{\mathrm{V}}^2}{{\mathrm{M}}_{\mathrm{SUSY}}/GeV}\right),$$

$${\tilde{\nu }}_L=\frac{1}{\sqrt{2}}\left({\varphi }_1+i{\sigma }_1\right),{\tilde{\nu }}_R^{\ast }=\frac{1}{\sqrt{2}}\left({\varphi }_2+i{\sigma }_2\right),\tilde{x}=\frac{1}{\sqrt{2}}\left({\varphi }_3+i{\sigma }_3\right),$$

$$\displaystyle m^{2}_{\tilde{\nu}}=\left(\begin{array}[]{ccc}m_{11}&m_{12}&m_{13}\\ m_{12}^{\ast}&m_{22}&m_{23}\\ m_{13}^{\ast}&m_{23}^{\ast}&m_{33}\end{array}\right),$$

$$m_{11}$$

$$m_{12}$$

$$m_{13}$$

$$m_{22}$$

$$m_{23}$$

$$m_{33}$$

$$C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1{h}_i}=C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1\mathrm{Re}[H_d^0]}{\tilde{U}}_{i1}+C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1\mathrm{Re}[H_u^0]}{\tilde{U}}_{i2}+C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1\mathrm{Re}[S]}{\tilde{U}}_{i3},$$

$$C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1\mathrm{Re}[H_d^0]}$$

$$C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1\mathrm{Re}[H_u^0]}$$

$$C_{{\tilde{\nu }}_1^{\ast }{\tilde{\nu }}_1\mathrm{Re}[S]}$$

Full text

1 Introduction

As the most popular ultraviolet-complete Beyond Standard Model, the Minimal Supersymmetric Standard Model (MSSM) with R-parity conservation predicts two kinds of electric neutral, possibly stable and weakly interactive massive particles, namely, sneutrino and neutralino, which may act as dark matter (DM) candidates [1, 2]. In the 1990s, it was proven that the left-handed sneutrino as the lightest supersymmetric particle (LSP) predicted a much smaller relic abundance than the measured value as well as an unacceptably tremendous DM-nucleon scattering rate due to its interaction with the Z boson [3, 4]. This fact made the lightest neutralino (usually with the bino field as its dominant component) the only reasonable DM candidate, and consequently, it was studied intensively since then. However, with the rapid progress in DM direct detection (DD) experiments in recent years, the candidate became more and more tightly limited by the experiments [5, 6, 7, 8] assuming that it was fully responsible for the measured relic density and the higgsino mass $\mu$ was less than 300 GeV, which was favored to predict the Z boson mass naturally [9]. These conclusions apply to the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [10, 11, 12], where the sneutrinos are purely left-handed, and the neutralino DM candidate may be either bino- or singlino-dominated [13]. In this context, we revived the idea of the sneutrino DM in a series of works [14, 8, 15, 16]. In particular, motivated by the phenomenology of the neutrino oscillations, we extended the NMSSM with the inverse seesaw mechanism by introducing two types of gauge singlet chiral superfields $\hat{\nu}_{R}$ and $\hat{X}$ for each generation matter, which have lepton numbers -1 and 1, respectively, and their fermion components corresponded to the massive neutrinos in literatures [14]. Subsequently, we studied in detail whether the $\tilde{\nu}_{R}$ (the scalar component of $\hat{\nu}_{R}$) or $\tilde{x}$ (the scalar component of $\hat{X}$) dominated sneutrino could act as a feasible DM candidate [14]. We were interested in the inverse seesaw mechanism because it was a TeV scale physics to account for the neutrino oscillations and maybe experimentally testable soon. We showed by both analytic formulas and numerical calculations that the resulting theory (abbreviated as ISS-NMSSM hereafter) was one of the most economic framework to generate the neutrino mass and, meanwhile, to reconcile the DM DD experiments naturally [14, 8]. We add that, besides us, a lot of authors have showed interest in the sneutrino DM in recent years [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36], but none of them considered the same theoretical framework as ours.

In the NMSSM, the introduction of the singlet field $\hat{S}$ can solve the $\mu$ problem of the MSSM [13], enhance the theoretical prediction of the SM-like Higgs boson mass [37, 38, 39], as well as enrich the phenomenology of the NMSSM significantly (see, for example, Ref. [40, 41, 42, 43, 44, 45]). In the ISS-NMSSM, the $\hat{S}$ field also plays extraordinary roles in generating the massive neutrino mass by the Yukawa interaction $\lambda_{\nu}\hat{S}\hat{\nu}_{R}\hat{X}$ and making the sneutrino DM compatible with various measurements, especially the DM DD experiments [14]. There are at least two aspects in manifesting the latter role. One is that the newly introduced heavy neutrino superfields are singlet under the gauge group of the SM model. Thus, they can interact directly with $\hat{S}$ by the Yukawa couplings [14]. In this case, the sneutrino DM candidate $\tilde{\nu}_{1}$, the singlet dominated scalars $h_{s}$ and $A_{s}$, and the massive neutrinos $\nu_{h}$ compose a roughly secluded DM sector where the annihilations $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to A_{s}A_{s},h_{s}h_{s},\nu_{h}\bar{\nu}_{h}$ can produce the measured relic density (In the ISS-NMSSM, these annihilations proceed by quartic scalar interactions, $s$-channel exchange of $h_{s}$ and $t/u$-channel exchange of the sneutrinos or the singlino-dominated neutralino). Since this sector communicates with the SM sector by the small singlet-doublet Higgs mixing (dubbed by Higgs-portal in literatures [46]) and/or by the massive neutrinos (neutrino-portal [48, 47, 49, 50]), the scattering of the DM with nucleons is naturally suppressed, which coincides with current DM DD results. The other aspect is that the singlet-dominated Higgs scalars can mediate the transition between $\tilde{\nu}_{1}$ pair and the higgsino pair, and consequently, these particles were in thermal equilibrium in early Universe before their freeze-out from the thermal bath. If their mass splitting is less than about $10\%$, the number density of the higgsinos can track that of $\tilde{\nu}_{1}$ during the freeze-out [51] (in literatures such a phenomenon was called co-annihilation [52]). Since, in this case, the couplings of $\tilde{\nu}_{1}$ with SM particles is usually very weak, the scattering is again naturally suppressed. We emphasize that, in either case, the suppression of the scattering prefers a small higgsino mass that appears in the coupling of $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}$ state with Higgs bosons, and hence, there is no tension between the DM DD experiments and the naturalness for the mass of the Z boson [14].

In the ISS-NMSSM, the rates of the DM annihilation and the DM-nucleon scattering depend on the coupling strength of $\tilde{\nu}_{1}$ interacting with Higgs fields, i.e., the Yukawa couplings $\lambda_{\nu}$ and $Y_{\nu}$ (the coefficient for $\hat{\nu}_{L}\cdot\hat{H}_{u}\hat{\nu}_{R}$ interaction) and their corresponding soft-breaking trilinear parameters $A_{\lambda_{\nu}}$ and $A_{Y_{\nu}}$. They also depend on the Higgs mass spectrum and the mixing between the Higgs fields that are ultimately determined by the parameters in the Higgs sector [14]. As such, the DM physics is quite complicated and is difficult to understand intuitively. This fact inspired us to study the theory from different aspects, e.g., from the features of the DM-nucleon scattering [14, 8], and its capability to explain the muon anomalous magnetic momentum [53] or other anomalies at the LHC [16]. In this work, we noted that a large $\lambda_{\nu}$ or $Y_{\nu}$ can enhance the DM-nucleon scattering rate significantly, so the recent XENON-1T experiment should limit them [54]. We also noted that the upper bound on the unitarity violation in neutrino sector sets a specific correlation between the couplings $\lambda_{\nu}$ and $Y_{\nu}$ [55], which can limit the parameter space of the ISS-NMSSM. Since these issues were not studied before, we decided to survey the impact of the leptonic unitarity and current and future DM DD experiments on the sneutrino DM sector. We will show that they are complementary to each other in limiting the theory, and in some cases, the constraints are rather tight. It is evident that such a study helps improve the understanding of the theory, and may be treated as a preliminary work before more comprehensive studies in the future.

We organize this work as follows. In section 2, we briefly introduce the theory of the ISS-NMSSM. In section 3, we describe the strategy to study the constraints, present numerical results and reveal the underlying physics. Finally, we draw our conclusions in section 4.

2 NMSSM with inverse seesaw mechanism

Since the ISS-NMSSM has been introduced in detail in [14, 8], we only recapitulate its key features in this section.

2.1 Basics of the ISS-NMSSM

The renormalizable superpotential and the soft breaking terms of the ISS-NMSSM take following form [14]

[TABLE]

where $W_{\rm MSSM}$ and $L_{\rm MSSM}^{soft}$ represent the corresponding terms of the MSSM without the $\mu$-term. The terms in the first brackets on the right side of each equation make up the Lagrangian of the NMSSM that involves the Higgs coupling coefficients $\lambda$ and $\kappa$ and their soft-breaking parameters $A_{\lambda}$ and $A_{\kappa}$. The terms in the second brackets are needed to implement the supersymmetric inverse seesaw mechanism. Coefficients such as the neutrino mass term $\mu_{X}$, the Yukawa couplings $\lambda_{\nu}$ and $Y_{\nu}$, and the soft-breaking parameters $A_{\lambda_{\nu}}$, $A_{Y_{\nu}}$, $B_{\mu_{X}}$, $m_{\tilde{\nu}}$, and $m_{\tilde{x}}$ are all $3\times 3$ matrices in flavor space. Besides, among the parameters in the superpotential, only the matrix $\mu_{X}$ is dimensional, and it parameterizes the effect of lepton number violation (LNV). Since this matric arises from the integration of massive particles in the high-energy ultraviolet theory with LNV interactions (see, for example, [56, 58, 57]), its magnitude should be small. Based on a similar perspective, the matrix $B_{\mu_{X}}$ is also theoretically favored to be suppressed.

It is the same as the NMSSM that the ISS-NMSSM predicts three CP-even Higgs bosons, two CP-odd Higgs bosons, a pair of charged Higgs bosons, and five neutralinos. Throughout this work, we take $\lambda$, $\kappa$, $\tan\beta\equiv v_{u}/v_{d}$, $A_{\lambda}$, $A_{\kappa}$, and $\mu\equiv\lambda v_{s}/\sqrt{2}$ as inputs of the Higgs sector, where $v_{u}\equiv\sqrt{2}\langle H_{u}\rangle$, $v_{d}\equiv\sqrt{2}\langle H_{d}\rangle$, and $v_{s}\equiv\sqrt{2}\langle S\rangle$ denote the vacuum expectation values (vev) of the fields $H_{u}$, $H_{d}$, and $S$, respectively. The elements of the CP-even Higgs fields’ squared mass matrix in the bases ($S_{1}\equiv\cos\beta{\rm Re}[H_{u}^{0}]-\sin\beta{\rm Re}[H_{d}^{0}]$, $S_{2}\equiv\sin\beta{\rm Re}[H_{u}^{0}]+\cos\beta{\rm Re}[H_{d}^{0}]$, $S_{3}\equiv{\rm Re}[S]$) are given by [13, 39]

[TABLE]

where $S_{1}$ denotes the heavy doublet Higgs field with a vanishing vev, $S_{2}$ represents the SM Higgs field with its vev $v\equiv 246~{}{\rm GeV}$, ${\cal M}_{22}$ is the mass of $S_{2}$ at tree level without considering its mixing with the other bases, and ${\cal M}_{23}$ characterizes the mixing of $S_{2}$ with the singlet field $S_{3}$.

The squared mass matrix in Eq. (2.1) can be diagonalized by a unitary matrix $U$, and its eigenstates $h_{i}$ with $i=1,2,3$ are obtained by

[TABLE]

where $h_{i}$ are labelled in an ascending mass order, i.e. $m_{h_{1}}<m_{h_{2}}<m_{h_{3}}$. Then the couplings of $h_{i}$ to vector bosons $W$ and $Z$ and fermions $u$ and $d$ quarks, which are normalized to their SM predictions, take the following form [59]

[TABLE]

Obviously, $U_{i2}\simeq 1$ if the components of the particle $h_{i}$ are far dominated by $S_{2}$, and consequently, $\bar{C}_{h_{i}V^{\ast}V}\simeq\bar{C}_{h_{i}\bar{u}u}\simeq\bar{C}_{h_{i}\bar{d}d}\simeq 1$. We call $h_{i}$ as the SM-like Higgs boson.

Similarly, the elements of the CP-odd Higgs fields’ squared mass matrix are [13]

[TABLE]

in the bases ($A\equiv\cos\beta{\rm Im}[H_{u}^{0}]+\sin\beta{\rm Im}[H_{d}^{0}]$, ${\rm Im}[S]$). As a result, the two CP-odd mass eigenstates $A_{1}$ and $A_{2}$ are the mixtures of $A$ and ${\rm Im}[S]$. The charged Higgs are given by $H^{\pm}=\cos\beta H_{u}^{\pm}+\sin\beta H_{d}^{\pm}$, and their masses are $m_{H^{\pm}}=2\mu(\lambda A_{\lambda}+\kappa\mu)/(\lambda\sin 2\beta)+v^{2}(g^{2}/2-\lambda^{2})$.

Concerning the neutralinos, they are the mixtures of the bino field $\tilde{B}^{0}$, the wino field $\tilde{W}^{0}$, the Higgsinos fields $\tilde{H}_{d}^{0}$ and $\tilde{H}_{u}^{0}$, and the singlino field $\tilde{S}^{0}$. In the bases $\psi^{0}=(-i\tilde{B}^{0},-i\tilde{W}^{0},\tilde{H}_{d}^{0},\tilde{H}_{u}^{0},\tilde{S}^{0})$, their mass matrix is given by [13]

[TABLE]

where $M_{1}$ and $M_{2}$ are soft breaking masses of the gauginos. It can be diagonalized by a rotation matrix $N$ so that the mass eigenstates are

[TABLE]

It is evident that $N_{i3}$ and $N_{i4}$ characterize the $\tilde{H}_{d}^{0}$ and $\tilde{H}_{u}^{0}$ components in $\tilde{\chi}_{i}^{0}$, respectively, and $N_{i5}$ denotes the singlino component.

In this work, we use the following features in the Higgs and neutralino sectors:

• •

A CP-even state corresponds to the SM-like Higgs boson discovered at the LHC. This state is favored to be ${\rm Re}[H_{u}^{0}]$-dominated by the LHC data when $\tan\beta\gg 1$, and its mass may be significantly affected by the interaction $\lambda\,\hat{s}\,\hat{H}_{u}\cdot\,\hat{H}_{d}$, the doublet-singlet Higgs mixing as well as the radiative correction from top/stop loops [37, 38, 39]. In the following, we denote this state as $h$.

• •

In most cases, the heavy doublet-dominated CP-even state is mainly composed of the field ${\rm Re}[H_{d}^{0}]$. It roughly degenerates in mass with the doublet-dominated CP-odd state and also with the charged states. The LHC search for extra Higgs bosons and the $B$-physics measurements requires these states to be heavier than about $500~{}{\rm GeV}$ [60]. We represent them by $H$, $A_{H}$, and $H^{\pm}$.

• •

Concerning the singlet-dominated states, they may be very light without conflicting with any collider constraint. As we introduced before, these states may appear as the final state of the sneutrino pair annihilation or mediate the annihilation, and thus, they can play a vital role in the sneutrino DM physics. In this work, we label these states by $h_{s}$ and $A_{s}$.

• •

The lightest neutralino $\tilde{\chi}_{1}^{0}$ is Higgsino dominated if $|\mu|<|M_{1}|,|M_{2}|$ and $|2\kappa/\lambda|<1$. In this case, $|N_{13}|\simeq|N_{14}|\simeq\sqrt{2}/2$.

We add that, to study the property of the sneutrino DM, we consider two benchmark scenarios where all the input parameters for the Higgs and neutralino sectors are fixed. The details of the scenarios are presented in Table 1. For the first scenario, $h_{s}$ and the SM-like Higgs boson $h$ correspond to the lightest and the next-to-lightest CP-even Higgs bosons $h_{1}$ and $h_{2}$. The $S_{2}$ component of $h_{s}$ is measured by the rotation element $U_{12}$, which is determined by the elements ${\cal{M}}_{23}^{2}$ and ${\cal{M}}_{33}^{2}$ in Eq.(2.1). We dub this scenario light $h_{s}$ scenario. By contrast, we call the second scenario as the massive $h_{s}$ scenario. It predicts $h=h_{1}$, $h_{s}=h_{2}$, and $U_{22}$ to characterize the $S_{2}$ component in $h_{s}$. Besides, we note that triple Higgs interactions may play an essential role in the sneutrino DM annihilation. So in addition to the couplings $\bar{C}_{h_{i}V^{\ast}V}$, $\bar{C}_{h_{i}\bar{u}u}$ and $\bar{C}_{h_{i}\bar{d}d}$, we also list in Table 1 the coupling strengths for $h_{s}h_{s}h_{s}$, $h_{s}h_{s}h$ and $h_{s}hh$ interactions, which are normalized to the triple Higgs coupling in the SM and denoted by $\bar{C}_{h_{s}h_{s}h_{s}}$, $\bar{C}_{h_{s}h_{s}h}$ and $\bar{C}_{h_{s}hh}$, respectively. These strengths are obtained by the formulas in [13]. They are characterized by $|\bar{C}_{h_{s}h_{s}h_{s}}|\gg|\bar{C}_{h_{s}h_{s}h}|,|\bar{C}_{h_{s}hh}|$, which is evident by the superpotential and the soft breaking terms of the ISS-NMSSM.

2.2 Leptonic Unitarity

In the interaction bases $(\nu_{L},\nu_{R}^{\ast},x)$, the neutrino mass matrix is given by [14]

[TABLE]

where both the Dirac mass $M_{D}=Y_{\nu}v_{u}/\sqrt{2}$ and the Majorana mass $M_{R}=\lambda_{\nu}v_{s}/\sqrt{2}$ are $3\times 3$ matrix in the flavor space. One can diagonalize this mass matrix by a $9\times 9$ unitary matrix $U_{\nu}$ to obtain three light neutrinos $\nu_{i}$ $(i=1,2,3)$ and six massive neutrinos $\nu_{h}$ as mass eigenstates, i.e., $U^{\ast}_{\nu}M_{\mathrm{ISS}}U^{\dagger}_{\nu}={\rm diag}(m_{\nu_{i}},m_{\nu_{h}})$, and decompose $U_{\nu}$ into the following blocks:

[TABLE]

The sub-matrix $\hat{U}_{3\times 3}$ encodes the neutrino oscillation information and it is determined by the neutrino experimental results.

Alternatively, one can get the analytic expression of the light active neutrinos’ mass matrix from Eq. (2.10) in the limit $\|\mu_{X}\|\ll\|M_{D}\|\ll\|M_{R}\|$, where $\|M\|$ is defined by $\|M\|\equiv\sqrt{{\rm Tr}(M^{\dagger}M)}$ for an arbitrary matrix $M$. The result is

[TABLE]

where $F=M_{D}^{T}M_{R}^{T^{-1}}$, and its elements’ magnitude is of the order $\|M_{D}\|/\|M_{R}\|$. This $3\times 3$ matrix can be diagonalized by the unitary Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, i.e.,

[TABLE]

Due to the mixings among the states $(\nu_{L},\nu_{R}^{\ast},x)$, the matrix $\hat{U}$ in Eq.(2.13) does not coincide with $U_{\rm PMNS}$. Instead, they are related by [61]

[TABLE]

where $\eta=\frac{1}{2}FF^{\dagger}$ is a measure of the non-unitarity for the matrix $\hat{U}$. A recent global fit of the theory to low energy experimental data reveals that [61]

[TABLE]

We call these inequalities as the leptonic unitarity constraint.

Eq. (2.14) indicates that the tininess of the active neutrino masses in the inverse seesaw mechanism is due to the smallness of the lepton-number violating matrix $\mu_{X}$ and the suppression factor $\|M_{D}\|^{2}/\|M_{R}\|^{2}$. For $\|\mu_{X}\|\lesssim{\cal O}({\rm keV})$ and $\|M_{R}\|\sim{\cal O}$(TeV), the magnitude of the Dirac Yukawa coupling $Y_{\nu}$ may reach order one in predicting $m_{\nu_{i}}\sim 0.1~{}{\rm eV}$. However, a large $Y_{\nu}$ may conflict with the unitarity constraint once the Majorana mass $M_{R}$ is specified. So the constraint must be taken into account in phenomenological study.

In the following, we will discuss the application of the unitarity constraint in DM physics. After noticing that the neutrino oscillation phenomenon can be explained by choosing an appropriate $\mu_{X}$ [62, 55], we assume flavor diagonal $Y_{\nu}$ and $\lambda_{\nu}$ to simplify the DM physics (see discussion below). We then determine $\mu_{X}$ by the formula [62, 55]

[TABLE]

where $m_{\nu_{i}}$ and $\hat{U}\simeq U_{\rm PMNS}$ take the values extracted from relevant neutrino experiments. With the assumption, the neutrino oscillation is solely attributed to the non-diagonality of $\mu_{X}$, and the unitarity constraint in Eq. (2.17) becomes

[TABLE]

These inequalities reveal that the ratio $[\lambda_{\nu}]_{33}/[Y_{\nu}]_{33}$ may be significantly smaller than $[\lambda_{\nu}]_{11}/[Y_{\nu}]_{11}$ and $[\lambda_{\nu}]_{22}/[Y_{\nu}]_{22}$ for fixed $\lambda$, $\mu$, and $v_{u}$, or equally speaking, $[Y_{\nu}]_{33}$ may be much larger than $[Y_{\nu}]_{11}$ and $[Y_{\nu}]_{22}$ when $\lambda_{\nu}$ is proportional to identity matrix.

Concerning the LNV coefficients $\mu_{X}$ and $B_{\mu_{X}}$, one should note two points. One is that $B_{\mu_{X}}$ can induce an effective $\mu_{X}$ through sneutrino-singlino loops to significantly affect the active neutrino masses by Eq. (2.14). We estimate the correction by the mass insertion method, which was widely used in B physics study. We find

[TABLE]

where $M^{2}$ parameterizes the mixing of the field $\tilde{\nu}_{R}^{\ast}$ with the field $\tilde{x}$, and $M_{\rm SUSY}$ represents the sparticles’ mass scale. Under the premise that $Y_{\nu}$, $\lambda_{\nu}$, and $B_{\mu_{X}}$ are flavor diagonal, $M^{2}$ can be roughly flavor diagonal, too (see the discussion in the next section). So one can study the correction in one generation case. The result is

[TABLE]

which indicates that $B_{\mu_{X}}/{\rm GeV^{2}}$ may be comparable with $M_{\rm SUSY}/{\rm GeV}$ in getting $\delta\mu_{X}\sim 1~{}{\rm keV}$ for $\lambda_{\nu}=0.1$ and $M^{2}\lesssim 0.1\times M_{\rm SUSY}^{2}$. Alternatively, if $\lambda_{\nu}=0.5$ and $M^{2}=0.1\times M_{\rm SUSY}^{2}$, the approximation requires $B_{\mu_{X}}/{\rm GeV^{2}}<0.1\times M_{\rm SUSY}/{\rm GeV}$ to get $\delta\mu_{X}\sim 1~{}{\rm keV}$. These estimations provide an upper bound on $B_{\mu_{X}}$’s magnitude. In our study, we limit $B_{\mu_{X}}\leq 100~{}{\rm GeV^{2}}$ for simplicity. The other point is that the LNV coefficients may induce sizable neutrinoless double beta decay since the inverse seesaw scale may be around several hundred GeV and the Yukawa couplings $Y_{\mu}$ and $\lambda_{\nu}$ may be moderately large. As indicated in [63], because $\mu_{X}$ is related to the active neutrino mass, the decay rate is below current experiment sensitivity when the massive neutrinos are heavier than $1~{}{\rm GeV}$. So there is no need to consider the constraint in our study.

2.3 Properties of Sneutrino Dark Matter

If the sneutrino fields are decomposed into CP-even and CP-odd parts

[TABLE]

the squared mass of the CP-even fields is given by

[TABLE]

in the bases $(\phi_{1},\phi_{2},\phi_{3})$, where

[TABLE]

These formulas indicate the following facts:

• •

The squared mass is a $9\times 9$ matrix in three-generation $(\phi_{1},\phi_{2},\phi_{3})$ bases. It involves a series of $3\times 3$ matrices in the flavor space, such as $Y_{\nu}$, $\lambda_{\nu}$, $A_{Y_{\nu}}$, $A_{\lambda_{\nu}}$, $\mu_{X}$, $B_{\mu_{X}}$, $m_{\tilde{l}}$, $m_{\tilde{\nu}}$, and $m_{\tilde{x}}$. Among these matrices, only $\mu_{X}$ must be flavor non-diagonal to account for the neutrino oscillations, but since its magnitude is less than 10 keV [62], it can be neglected. Thus, if there is no flavor mixings for the other matrices, the squared mass is flavor diagonal, and one can adopt one-generation $(\phi_{1},\phi_{2},\phi_{3})$ bases in studying the mass. In this work, we only consider the third generation sneutrinos as DM candidates. This is motivated by that both the unitarity bound and the LHC constraint in sparticle search are weakest for the third generation [14]. When we mention the sneutrino parameters hereafter, we are actually referring to their 33 elements. Under the assumption, the squared mass is diagonalized by a $3\times 3$ unitary matrix $V$, which parameterizes the chiral mixings between the fields $\phi_{1}$, $\phi_{2}$ and $\phi_{3}$. Consequently, the sneutrino mass eigenstates are given by $\tilde{\nu}_{R,i}=V_{ij}\phi_{j}$ with $i,j=1$, 2, and 3. we add that $Y_{\nu}$ and $\lambda_{\nu}$ are real and positive numbers after properly rotating the phase of the fields $\hat{\nu}_{R}$ and $\hat{X}$.

• •

The mixing of $\phi_{1}$ with the other fields is determined by $Y_{\nu}$ and $A_{\nu}$. As $Y_{\nu}$ approaches zero, $|m_{12}|$ and $|m_{13}|$ diminish monotonically, and so is $|V_{11}|$ which represents the $\tilde{\nu}_{L}$ component in the lightest sneutrino state $\tilde{\nu}_{R,1}$. In the extreme case $Y_{\nu}=0$, all these quantities vanish and $\tilde{\nu}_{R,1}$ is merely the mixture of $\phi_{2}$ and $\phi_{3}$. Furthermore, if $\lambda_{\nu}/\lambda$ is moderately large, the first term in $m_{22}$ and $m_{33}$ may be far dominant over the other contributions so that $m_{22}\simeq m_{33}$. This results in maximal mixing between $\phi_{2}$ and $\phi_{3}$ and $\tilde{\nu}_{R,1}\simeq 1/\sqrt{2}[\phi_{2}-{\rm sgn}(m_{23})\phi_{3}]$ [14]. This is a case encountered frequently in our study.

Similarly, one may adopt the one-generation $(\sigma_{1},\sigma_{2},\sigma_{3})$ bases to study the CP-odd sneutrino’s mass, which is the same as Eq. (2.25) except for the substitution $B_{\mu_{X}}\to-B_{\mu_{X}}$. The mass eigenstates are then given by $\tilde{\nu}_{I,i}=V_{ij}^{\prime}\sigma_{j}$ ($i,j=1,2,3$) where $V^{\prime}$ denotes the rotation of the CP-odd fields. Given that $B_{\mu_{X}}$ represents the degree of the LNV and is theoretically small, we are particularly interested in the following two cases:

• •

The extreme case of $B_{\mu_{X}}=0$ where any CP-even sneutrino state is accompanied by a mass-degenerate CP-odd state. In this case, any sneutrino mass eigenstate corresponds to a complex field and it has an anti-particle [64]. Concerning the sneutrino DM $\tilde{\nu}_{1}$, we have $\tilde{\nu}_{R,1}\equiv{\rm Re}[\tilde{\nu}_{1}]$, $\tilde{\nu}_{I,1}\equiv{\rm Im}[\tilde{\nu}_{1}]$, $V_{ij}=V^{\prime}_{ij}$, and $\tilde{\nu}_{1}$ and its anti-particle $\tilde{\nu}_{1}^{\ast}$ contribute equally to the relic density. This case is actually a two-component DM theory. It is notable that the $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}Z$ coupling is proportional to $|V_{11}|^{2}$ and it contributes to the scattering of $\tilde{\nu}_{1}$ with nucleons. This effect is important when $|V_{11}|\sim 0.01$ (discussed below). It is also notable that the $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}A_{i}$ coupling vanishes since it is induced only by the LNV effect.

• •

A more general case satisfying $|B_{\mu_{X}}|\lesssim 100~{}{\rm GeV^{2}}$. It has four distinctive features. First, since $m_{\tilde{\nu}_{R,i}}>m_{\tilde{\nu}_{I,i}}$ when $B_{\mu_{X}}>0$, the DM candidate $\tilde{\nu}_{1}$ is identified as the $\tilde{\nu}_{I,1}$ state with a definite CP number -1. The opposite conclusion applies to $B_{\mu_{X}}<0$ case. Second, any CP-even state is slightly different from its CP-odd partner in mass, e.g., $m_{\tilde{\nu}_{R,1}}-m_{\tilde{\nu}_{I,1}}\simeq 0.2~{}{\rm GeV}$ when $B_{\mu_{X}}=100~{}{\rm GeV^{2}}$ and $m_{\tilde{\nu}_{R,1}}=100~{}{\rm GeV}$, and so are the rotations $V$ and $V^{\prime}$. These sneutrino states compose a pseudo-complex particle  [67, 58, 68]. Third, given the approximate mass degeneracy, $\tilde{\nu}_{R,1}$ and $\tilde{\nu}_{I,1}$ always co-annihilated in early universe to affect the DM density. We will discuss this issue later. Finally, $Z$ boson does not mediate the DM-nucleon scattering any more since it couples only to a pair of sneutrino states with opposite CP numbers. It also contributes little to the DM annihilation because the $\tilde{\nu}_{R,1}\tilde{\nu}_{I,1}Z$ coupling is suppressed by a factor $V_{11}^{\ast}V^{\prime}_{11}\simeq|V_{11}|^{2}$.

We fix $B_{\mu_{X}}=0$ or $B_{\mu_{X}}=-100~{}{\rm GeV^{2}}$ in this work. In either case, the $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}h_{i}$ ($h_{i}=h_{s},h,H$) coupling coefficient is given by

[TABLE]

where $\tilde{U}$ diagonalizes the CP-even Higgs fields’ squared mass in $({\rm Re}[H_{d}^{0}],{\rm Re}[H_{u}^{0}],{\rm Re}[S])$ bases [38, 39], and $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}s}$ on the right side denotes the sneutrino coupling to the scalar field $s$. For the one-generation sneutrino case, $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}s}$ is given by

[TABLE]

These formulas indicate that the parameters $Y_{\nu}$, $\lambda_{\nu}$, $A_{Y_{\nu}}$, and $A_{\lambda_{\nu}}$ affect not only the sneutrino interactions but also their mass spectrum and mixing. In particular, a large $\lambda_{\nu}$ or $Y_{\nu}$ can enhance the coupling significantly. Instead, the soft-breaking masses $m_{\tilde{\nu}}^{2}$ and $m_{\tilde{x}}^{2}$ affect only the latter property. For typical values of the parameters in Eq. (2.27), e.g., $\tan\beta\gg 1$, $|V_{11}|<0.1$, $Y_{\nu},\kappa,\lambda,\lambda_{\nu}\sim{\cal{O}}(0.1)$ and $\lambda_{\nu}v_{s},\lambda v_{s},A_{Y_{\nu}},A_{\lambda_{\nu}}\sim{\cal{O}}({\rm 100~{}GeV})$, $C_{\tilde{\nu}_{1}\tilde{\nu}_{1}s}$ is approximated by

[TABLE]

It is estimated that $|C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[H_{d}^{0}]}|,|C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[H_{u}^{0}]}|\lesssim 10~{}{\rm GeV}$ and $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[S]}\lesssim 100~{}{\rm GeV}$ in most cases, which reflects that $|C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[S]}|$ may be much larger than the other two couplings. The basic reason is that $\tilde{\nu}_{1}$ is a singlet-dominated scalar, so it can couple directly to the field $S$ and the mass dimension of $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[S]}$ is induced by $v_{s}$ or $A_{\lambda_{\nu}}$. By contrast, the other couplings emerge only after the electroweak symmetry breaking when $V_{11}=0$, and their mass dimension originates from $v_{u}$.

2.4 Relic Density of Sneutrino Dark Matter

In the $B_{\mu_{X}}=0$ case, both $\tilde{\nu}_{1}$ and $\tilde{\nu}_{1}^{\ast}$ act as the DM candidate. Their annihilation includes those initiated by $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}$, $\tilde{\nu}_{1}\tilde{\nu}_{1}$, and $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}^{\ast}$ state, and the co-annihilation of $\tilde{\nu}_{1}$ and $\tilde{\nu}_{1}^{\ast}$ with the other sparticles. Considering the numerousness of the annihilation channels and the complexity of this issue, we will only discuss the channels frequently met in our study (see footnote 2 of this work for more details), which are [14, 65]:

• (1)

$\tilde{\nu}_{1}\tilde{H},\tilde{\nu}_{1}^{\ast}\tilde{H}\rightarrow XY$ and $\tilde{H}\tilde{H}^{\prime}\rightarrow X^{\prime}Y^{\prime}$, where $\tilde{H}$ and $\tilde{H}^{\prime}$ denote Higgsino-dominated neutralinos or charginos, and $X^{(\prime)}$ and $Y^{(\prime)}$ represent any possible SM particles, the massive neutrinos or the Higgs bosons if the kinematics are accessible. More specifically, the channels $\tilde{\nu}_{1}\tilde{H}\rightarrow Wl,Z\nu,h\nu$ ($l$ and $\nu$ denote any possible lepton and neutrino, respectively) proceed by the $s$-channel exchange of neutrinos, and the $t/u$ channel exchange of sleptons or sneutrinos. The processes $\tilde{H}\tilde{H}^{\prime}\rightarrow f\bar{f}^{\prime},VV^{\prime},hV$ ($f$ and $f^{\prime}$ denote quarks or leptons, and $V$ and $V^{\prime}$ represent SM vector bosons) proceed by the $s$-channel exchange of vector bosons or Higgs bosons, and the $t/u$ channel exchange of sfermions, neutralinos or charginos. This annihilation mechanism is called co-annihilation [51, 52].

• (2)

$\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\rightarrow ss^{\ast}$ ($s$ denotes a light Higgs boson), which proceeds through any relevant quartic scalar coupling, the $s$-channel exchange of CP-even Higgs bosons, and the $t/u$-channel exchange of sneutrinos.

• (3)

$\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\rightarrow\nu_{h}\bar{\nu}_{h}$ via the $s$-channel exchange of CP-even Higgs bosons or the $t/u$-channel exchange of neutralinos, where $\nu_{h}$ denotes a massive neutrino.

• (4)

$\tilde{\nu}_{1}\tilde{\nu}_{1}\rightarrow\nu_{h}\nu_{h}$ and $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}^{\ast}\rightarrow\bar{\nu}_{h}\bar{\nu}_{h}$, which mainly proceed through the $t/u$-channel exchange of a singlino-dominated neutralino due to its majorana nature.

• (5)

$\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\rightarrow VV^{\ast}$, $Vs$, $f\bar{f}$, which proceeds mainly by the $s$-channel exchange of CP-even Higgs bosons. They are important if one of the bosons is at resonance.

Under specific parameter configurations, these channels can be responsible for the DM density precisely measured by the Planck experiment [66]. In this aspect, we have the following observations (see footnote 2 for more explanations):

• •

In most cases, the DMs annihilated mainly through the co-annihilation to get the measured density. This mechanism works only when the mass splitting between $\tilde{H}$ and $\tilde{\nu}_{1}$ is less than about $10\%$, and a specific channel’s contribution to the density depends not only on its cross-section but also on the mass splitting. To illustrate this point, we assume that the DM annihilations comprise those initiated by $\tilde{\nu}_{1}\tilde{\nu}_{1}$, $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}$, $\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}^{\ast}$, $\tilde{\nu}_{1}\tilde{\chi}_{1}^{0}$, $\tilde{\nu}_{1}^{\ast}\tilde{\chi}_{1}^{0}$, and $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ states, and denote the cross-sections of these channels by $\sigma_{AB}$ with $A,B=\tilde{\nu}_{1},\tilde{\nu}_{1}^{\ast},\tilde{\chi}_{1}^{0}$. The effective annihilation rate at temperature $T$ is then given by [52]

[TABLE]

where $\Delta\equiv(m_{\tilde{\chi}_{1}^{0}}-m_{\tilde{\nu}_{1}^{0}})/m_{\tilde{\nu}_{1}^{0}}$ parameterizes the mass splitting and $x\equiv m_{\tilde{\nu}_{1}^{0}}/T$. This formula indicates that the $\tilde{\nu}_{1}\tilde{\chi}_{1}^{0}$ and $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ channel’s contributions are suppressed by factors $e^{-x\Delta}$ and $e^{-2x\Delta}$, respectively. So they become less and less critical as $m_{\tilde{\nu}_{1}^{0}}$ deviates from $m_{\tilde{\chi}_{1}^{0}}$. Besides, the formulae of the density in [52] indicate that the density depends on the sneutrino parameters only through $m_{\tilde{\nu}_{1}}$ and $\sigma_{eff}$. In the extreme case of $\sigma_{\tilde{\nu}_{1}\tilde{\nu}_{1}}\simeq\sigma_{\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}}\simeq\sigma_{\tilde{\nu}_{1}\tilde{\chi}^{0}_{1}}\simeq 0$ realized when $\lambda_{\nu}$ and $Y_{\nu}$ are sufficiently small, the $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ annihilation is solely responsible for the measured density through tuning the value of $m_{\tilde{\nu}_{1}}$. This situation was intensively studied in [8]. We will present such examples in Section III.

• •

Barring the co-annihilation, $\tilde{\nu}_{1}\tilde{\nu}_{1}\rightarrow ss^{\ast}$ is usually the most crucial channel in affecting the density if the kinematics are accessible. In particular, the process $\tilde{\nu}_{1}\tilde{\nu}_{1}\rightarrow h_{s}h_{s}$ can be solely responsible for the measured density if the Yukawa coupling $\lambda_{\nu}$ is moderately large. We exemplify this point by considering the light $h_{s}$ scenario in Table 1. From the Higgs boson and sneutrino mass spectrum and the $\tilde{\nu}_{1}$’s couplings to $h_{s}$, one can learn that the annihilation proceeds mainly by the $s$-channel exchange of $h_{s}$, $t/u$-channel exchange of $\tilde{\nu}_{1}$, and $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}h_{s}h_{s}$ quartic scalar coupling. As a result, the cross-section of the annihilation near the freeze-out temperature is approximated by [14, 65]

[TABLE]

where

[TABLE]

The measured density then requires $a+3b/25\simeq 4.6\times 10^{-26}{\rm cm}^{2}$ because we are considering a two-component DM theory [69, 70]. This requirement limits the $\tilde{\nu}_{1}$’s couplings to $h_{s}$ or for the fixed parameters in Table 1, ultimately the Yukawa coupling $\lambda_{\nu}$ since the cross-section is very sensitive to $\lambda_{\nu}$. We estimate that $\lambda_{\nu}\sim 0.4$ for $m_{\tilde{\nu}_{1}}=130~{}{\rm GeV}$ can account for the measured density.

• •

The process $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to\nu_{h}\bar{\nu}_{h}$ could be responsible for the density when $m_{\tilde{\nu}_{1}}>\nu_{h}$, $m_{\tilde{\nu}_{1}}<h_{s}$, and the co-annihilation mechanism did not work. This process proceeded mainly by the $s$-channel exchange of $h_{s}$, and consequently, the cross-section at the freeze-out temperature $T_{f}$ takes the following form:

[TABLE]

which implies that the density limits non-trivially $\lambda_{\nu}$, $m_{\tilde{\nu}_{1}}$ and $m_{h_{s}}$.

• •

About the other channels, they usually played a minor role in determining the density. So we leave the discussion of them in our future works.

Concerning the $B_{\mu_{X}}\neq 0$ case, either $\tilde{\nu}_{R,1}$ or $\tilde{\nu}_{I,1}$ acts as the DM candidate. Since the mass splitting between the DM $\tilde{\nu}_{1}$ and $\tilde{\nu}_{1}^{\prime}$ (the partner of $\tilde{\nu}_{1}$ with a different CP number) is small, $\tilde{\nu}_{1}$ always co-annihilated with $\tilde{\nu}_{1}^{\prime}$ to get the measured density. The relevant annihilation included $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\prime}$ and $\tilde{\nu}_{1}^{\prime}\tilde{\nu}_{1}^{\prime}$ initiated processes, and they proceeded in a way similar to the previous discussion. We confirmed that the density is insensitive to $B_{\mu_{X}}$ for $|B_{\mu_{X}}|\leq 100~{}{\rm GeV^{2}}$, which can be inferred from Eq. (2.29). We also verified that the cross-section of the DM annihilation today is insensitive to $B_{\mu_{X}}$.

2.5 DM-nucleon Scattering

In the $B_{\mu_{X}}\neq 0$ case, the scattering of $\tilde{\nu}_{1}$ with nucleon $N$ ($N=p,n$) proceeds by the $t/u$-channel exchange of the CP-even Higgs bosons. Consequently, the spin independent (SI) cross-section is given by [14]

[TABLE]

where $\mu_{\rm red}=m_{N}/(1+m_{N}^{2}/m_{\tilde{\nu}_{1}}^{2})$ ($N=p,n$) represents the reduced mass of the nucleon with $m_{\tilde{\nu}_{1}}$, $F^{N}_{u}=f^{N}_{u}+\frac{4}{27}f^{N}_{G}$ and $F^{N}_{d}=f^{N}_{d}+f^{N}_{s}+\frac{2}{27}f^{N}_{G}$ are nucleon form factors with $f^{N}_{q}=m_{N}^{-1}\left<N|m_{q}q\bar{q}|N\right>$ and $f^{N}_{G}=1-\sum_{q}f^{N}_{q}$ for $q=u,d,s$ [2]. With the default setting of the package micrOMEGAs [71, 72, 73] for nucleon sigma terms, i.e., $\sigma_{\pi N}=34~{}{\rm MeV}$ and $\sigma_{0}=42~{}{\rm MeV}$ [74]111It is notable that $\sigma_{0}$ was replaced by the strangeness-nucleon sigma term, $\sigma_{s}\equiv m_{s}/(m_{u}+m_{d})\times(\sigma_{\pi N}-\sigma_{0})\simeq 12.4\times(\sigma_{\pi N}-\sigma_{0})$, in recent calculation of the nucleon form factor [72]. Compared with the previous calculation, this treatment changes significantly the strange quark content in nucleon $N$, $f_{s}^{N}$, but it change little $F_{u}^{N}$ and $F_{d}^{N}$. , one can conclude $F_{u}^{p}\simeq 0.15$ and $F_{d}^{p}\simeq 0.14$ for protons. Instead, if $\sigma_{\pi N}=59~{}{\rm MeV}$ [75, 76, 77] and $\sigma_{0}=57~{}{\rm MeV}$ [78] are adopted, the form factors become $F_{u}^{p}\simeq 0.16$ and $F_{d}^{p}\simeq 0.13$. These results reflect that different choices of $\sigma_{\pi N}$ and $\sigma_{0}$ can induce uncertainties of ${\cal{O}}(10\%)$ in $F_{u}^{p}$ and $F_{d}^{p}$, and it does not drastically change the cross-section. Besides, the default setting also predicts $F_{u}^{n}\simeq 0.15$ and $F_{d}^{n}\simeq 0.14$ for neutrons, which implies the relation $\sigma^{\rm SI}_{\tilde{\nu}_{1}-p}\simeq\sigma^{\rm SI}_{\tilde{\nu}_{1}-n}$ for the Higgs-mediated contribution.

To clarify the features of the cross-section, we assume $m_{H^{\pm}}\gtrsim 1{\rm TeV}$ and integrate out the heavy doublet Higgs field. As a result, the CP-even Higgs sector at the electroweak scale contains only the SM Higgs field $S_{2}={\rm\sin\beta}{\rm Re}[H_{u}^{0}]+{\rm\cos\beta}{\rm Re}[H_{d}^{0}]$ and the singlet field ${\rm Re}[S]$. We then calculate the scattering amplitude by the mass insertion method to get the following result:

[TABLE]

where $\delta=m_{h}^{2}/m_{h_{s}}^{2}-1$, and $\theta$ is the mixing angle of the $S_{2}$ field with ${\rm Re}[S]$ to form mass eigenstates. This formula reveals that if the terms in the second brackets are on the order of 0.1, which can be achieved if $\lambda_{\nu}$ and/or $Y_{\nu}$ in Eq. (2.28) are sufficiently large, the cross-section may reach the sensitivity of the recent XENON-1T experiment [54]. We will discuss this issue later.

Concerning the $B_{\mu_{X}}=0$ case, where the DM corresponds to a complex field, the $Z$-boson also mediates the elastic scattering of the DM with nucleons. Since the total SI cross-section in this case is obtained by averaging over $\tilde{\nu}_{1}N$ and $\tilde{\nu}_{1}^{\ast}N$ scatterings and the interferences between the $Z$- and the Higgs-exchange diagrams for the two scatterings have opposite signs [79], the SI cross-section is given by [4]

[TABLE]

where $\sigma_{N}^{h}$ is the same as before and the $Z$-mediated contributions are

[TABLE]

with $G_{F}$ and $\theta_{W}$ denoting the Fermi constant and the weak angle, respectively. Since $\sigma_{n}^{Z}$ is larger than $\sigma_{p}^{Z}$ by a factor around 100, $\sigma_{n}^{\rm SI}$ may differ significantly from $\sigma_{p}^{\rm SI}$. Correspondingly, one may define the effective cross-section for the coherent scattering of the DMs with xenon nucleus as $\sigma_{\rm eff}^{\rm SI}=(\sigma^{\rm SI}_{\tilde{\nu}_{1}-Xe}+\sigma^{\rm SI}_{\tilde{\nu}_{1}^{\ast}-Xe})/(2A^{2})$, where $A$ denotes the mass number of the xenon nucleus, and calculate it by

[TABLE]

where the three coefficients on the right side are obtained by averaging the abundance of different xenon isotopes in nature. It is evident that the effective cross-section is identical to $\sigma^{\rm SI}_{p}$ if $\sigma_{p}^{\rm SI}=\sigma_{n}^{\rm SI}$, and it is related directly with the bound of the XENON-1T experiment [54].

Before concluding the introduction of the sneutrino DM, we add that its spin dependent cross-section is always zero, and its SI cross-section is usually much smaller than that of the neutralino DM in the MSSM and NMSSM, which was discussed in detail in Ref.[14, 8]. As a result, the extension is readily consistent with the XENON-1T experiment except for large $\lambda_{\nu}$ and/or $Y_{\nu}$ case studied in this work.

3 Constraints on sneutrino DM sector

In this section, we clarify the impact of the leptonic unitarity and current and future DM DD experiments on the sneutrino DM sector under the premise that the theory predicts the right density and the photon spectrum from the DM annihilation in dwarf galaxies is compatible with the Fermi-LAT observation. Since the singlet-dominated Higgs boson, $h_{s}$, plays a vital role in the density and the DM-nucleon scattering, we study the DM physics in both the light and the massive $h_{s}$ scenarios in Table 1. We emphasize that fixing the parameters in the Higgs and neutralino sectors can simplify greatly the analysis of the impact and make the underlying physics clear. We also emphasize that the two scenarios were obtained by scanning intensively the parameters in the Higgs and DM sectors222With the parameter scan strategy reported in [15] for the Type-I NMSSM, we explored the parameter space of the ISS-NMSSM which takes $\tan\beta$, $\lambda$, $\kappa$, $A_{t}$, $A_{\kappa}$, $\mu$, $\lambda_{\nu}$, $Y_{\nu}$, $A_{\lambda_{\nu}}$, $A_{Y_{\nu}}$, $m_{\tilde{\nu}}$, and $m_{\tilde{x}}$ as inputs. For either the light or the massive $h_{s}$ scenario, we have studied more than fifty million samples. The settings in Table 1 were chosen from the samples that best fit the experimental data. We will present the analysis of these samples elsewhere.. They agree well with the latest Higgs data of the LHC if the exotic decays $h\to\nu_{h}\bar{\nu}_{h},\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}$ are kinematically forbidden. This was confirmed by the packages HiggsSignal-2.4.0 [80] and HiggsBounds-5.7.0 [81].

3.1 Research strategy

The procedure of our study is as follows. We constructed a likelihood function of the DM physics to guide sophisticated scans over the sneutrino parameters for either scenario. With the samples obtained in the scans, we plotted the profile likelihood map in different two-dimensional planes to illustrate its features and underlying physics. We express the likelihood function as

[TABLE]

where $\mathcal{L}_{\Omega_{\tilde{\nu}_{1}}}$, $\mathcal{L}_{\rm DD}$, $\mathcal{L}_{\rm ID}$, and $\mathcal{L}_{Unitary}$ describe the relic density, the current XENON-1T experiment [54] or the future LZ experiment [82], the Fermi-LAT observation of dwarf galaxies, and the unitarity constraint, respectively. They are given by

• •

$\mathcal{L}_{\Omega_{\tilde{\nu}_{1}}}$ is Gaussian distributed, i.e.,

[TABLE]

where $\Omega_{th}$ denotes the theoretical prediction of the density $\Omega_{\tilde{\nu}_{1}}h^{2}$, $\Omega_{\rm obs}=0.120$ represents its experimental central value [66], and $\sigma=0.1\times\Omega_{\rm obs}$ is the total (including both theoretical and experimental) uncertainty of the density.

• •

$\mathcal{L}_{\rm DD}$ takes a Gaussian distributed form with a mean value of zero [83]:

[TABLE]

In this formula, $\sigma^{\rm SI}_{\rm eff}$ is defined in Eq. (2.35) and its error bar $\delta_{\sigma}$ is evaluated by $\delta_{\sigma}=\sqrt{UL_{\sigma}^{2}/1.64^{2}+(0.2\sigma^{\rm SI}_{\rm eff})^{2}}$, where $UL_{\sigma}$ denotes experimental upper limits on the scattering cross-section at $90\%$ C.L. and $0.2\sigma^{\rm SI}_{\rm eff}$ parameterizes the theoretical uncertainty of $\sigma^{\rm SI}_{\rm eff}$.

• •

$\mathcal{L}_{\rm ID}$ is calculated by the likelihood function proposed in [84, 85] with the data of the Fermi-LAT collaboration taken from [86, 87].

• •

The likelihood function of the unitarity constraint in Eq. (2.18) is as follows:

[TABLE]

where $r\equiv\lambda_{\nu}\mu/(Y_{\nu}\lambda v_{u})$ and $0.2r$ parameterizes total uncertainties.

In addition, we abandoned samples that open up the decays $h\to\nu_{h}\bar{\nu}_{h},\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}$333In fitting the ISS-NMSSM to experimental data, the total likelihood function is calculated by $\mathcal{L}_{\rm tot}=\mathcal{L}_{\rm Higgs}\times\mathcal{L}_{\rm DM}$, where $\mathcal{L}_{\rm Higgs}$ represents the Higgs physics function. Given $\chi^{2}\equiv-2\ln\mathcal{L}$, one can infer that $\chi^{2}_{\rm tot}=\chi^{2}_{\rm Higgs}+\chi^{2}_{\rm DM}$, and the $2\sigma$ confidence interval defined below Eq. (3.5) satisfies

$\displaystyle\chi^{2}_{\rm tot}-\chi^{2}_{\rm tot,min}\equiv\delta\chi^{2}_{\rm Higgs}+\delta\chi^{2}_{\rm DM}\leq 6.18,$

where $\delta\chi^{2}_{\rm Higgs}\equiv\chi^{2}_{\rm Higgs}-\chi^{2}_{\rm Higgs,min}$ and $\delta\chi^{2}_{\rm DM}\equiv\chi^{2}_{\rm DM}-\chi^{2}_{\rm DM,min}$ . For the Higgs parameter settings in Table 1, $\delta\chi^{2}_{\rm Higgs}$ vanishes because the settings correspond to the scenarios’ best points when the decays are forbidden. It increases if the decays are open and their branching ratios gradually increase. Since we hope to determine the interval only by the DM physics in the following study, we kinematically shut down the decays. . In practice, this was completed by setting the likelihood value to be $e^{-100}$ if any of the decays were kinematically accessible.

To make the conclusions in this study complete, we adopt the MultiNest algorithm [88, 89] to implement the scans. We take the prior probability density function (PDF) of the input parameters uniformly distributed and set nlive parameter of the algorithm to be 10000. This parameter represents the number of active or live points used to determine the iso-likelihood contour in each scan’s iteration [88, 89]. The larger it is, the more elaborated the scan becomes. The output of the scans includes the Bayesian evidence defined by

[TABLE]

where $P(\Theta|M)$ represents the prior PDF of the inputs $\Theta=(\Theta_{1},\Theta_{2},\cdots)$ in a model $M$, and $P(D|O(M,\Theta))\equiv\mathcal{L}(\Theta)$ denotes the likelihood function involving theoretical predictions of observables $O$ and their experimental measurements $D$. Computationally, the evidence is an averaged likelihood that depends on the priors of the theory’s input. In comparing different scenarios of the theory, the larger $Z$ is, the more readily the corresponding scenario is consistent with the data.

The output of the scan also includes the profile likelihood (PL) defined in frequentist statistics as the most significant likelihood value [90, 15]. For example, two-dimensional (2D) PL is defined by

[TABLE]

where the maximization is obtained by varying the parameters other than $\Theta_{A}$ and $\Theta_{B}$. The PL reflects the preference of the theory on the parameter $(\Theta_{A},\Theta_{B})$, or in other words, the capability of the parameter to account for experimental data. Sequentially, one can introduce the concept of confidence interval (CI) to classify the parameter region by how well the points in it fit the data. For example, the $1\sigma$ and $2\sigma$ CIs for the 2D PL are defined by satisfying $\chi^{2}-\chi^{2}_{\rm min}\leq 2.3$ and $\chi^{2}-\chi^{2}_{\rm min}\leq 6.18$, respectively, where $\chi^{2}\equiv-2\ln\mathcal{L}(\Theta_{A},\Theta_{B})$ and $\chi^{2}_{\rm min}$ is the minimal value of $\chi^{2}$ for the samples obtained in the scan.

In this work, we utilized the package SARAH-4.11.0 [91, 92, 93] to build the model file of the ISS-NMSSM, the SPheno-4.0.3 [94] code to generate its particle spectrum, and the package MicrOMEGAs 4.3.4 [95, 71, 73] to calculate the DM observables.

3.2 Results for the light $h_{s}$ scenario

Given the information of the light $h_{s}$ scenario in Table 1, the Higgs-mediated SI cross-section in Eq. (2.32) is approximated by

[TABLE]

Since $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[S]}$ may be larger than $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[H_{u}^{0}]}$ and $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[H_{d}^{0}]}$ by two orders, the first term in the brackets can be comparable with the other contributions. To clarify the impact of the unitarity and DM DD experiments on the theory, we performed four independent scans over the following parameter space:

[TABLE]

where $m_{\tilde{l}}$ denoted the common soft breaking mass of three-generation sleptons and its lower bound was motivated by the non-observation of slepton signals at the LHC Run-II. For the first scan, we fixed $B_{\mu_{X}}=-100~{}{\rm GeV}^{2}$ and used the XENON-1T’s bound on the SI cross-section to calculate the $\mathcal{L}_{\rm DD}$. The second scan was same as the first one except that we adopted the sensitivity of the LZ experiment. The last two scans differed from the previous ones only in that we set $B_{\mu_{X}}=0$. As explained before, the setting induces an additional $Z$-mediated contribution to the DM-nucleon scattering so that the constraints of the DD experiments are strengthened.

With the samples obtained in the scans, we show different 2D PL maps in Figures 1 to 5. Fig. 1 and 2 plot the CIs on $\lambda_{\nu}-m_{\tilde{\nu}_{1}}$ and $\sigma^{\rm SI}_{\tilde{\nu}_{1}-p}-m_{\tilde{\nu}_{1}}$ planes. They show the following features:

• •

$m_{\tilde{\nu}_{1}}$ is concentrated on the range from $120$ to $181~{}{\rm GeV}$. Specifically, $m_{\tilde{\nu}_{1}}$ is close to $m_{\tilde{\chi}_{1}^{0}}$ for $172~{}{\rm GeV}\lesssim m_{\tilde{\nu}_{1}}\lesssim 181~{}{\rm GeV}$, and the DM achieves the correct density mainly through the $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ annihilation (see discussions about the co-annihilation in Section 2.4 and details of the points in subsequent Table 2). In this case, the density is insensitive to the parameter $\lambda_{\nu}$. Thus, $\lambda_{\nu}$ varies within a broad range from 0.15 to 0.6 in Fig. 1, where the lower limit forbids the decay $h\to\nu_{h}\bar{\nu}_{h}$ kinematically, and the upper bound comes from the DM DD experiments (discussed below). For the other mass range, the DM obtains the correct density mainly through the annihilations $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to h_{s}h_{s},h_{s}h,hh$. This requires $\lambda_{\nu}\gtrsim 0.26$, which can be understood from the discussion of Eq. (2.30).

• •

Fig. 2 indicates that the SI cross-section of the DM-nucleon scattering may be as low as $10^{-49}~{}{\rm cm^{2}}$ over the entire mass range. It reflects that the theory has multiple mechanisms to suppress the scattering, which becomes evident by the approximation in Eq. (3.6) and was recently emphasized in [8].

• •

Although $\lambda_{\nu}>0.6$ is allowed by the setting in Eq. (3.7), it is upper bounded by 0.56, 0.56, 0.50, and 0.45 for the $2\sigma$ CIs of the four cases, respectively. Careful comparisons of the left and right panels revealed that it was due to the DD experiments’ constraint on the co-annihilation region. Besides, we studied the Bayesian evidences $Z_{i}$ ($i=1,2,3,4$) of the four cases and found $\ln Z_{1}=-55.6$, $\delta_{12}\equiv\ln Z_{1}-\ln Z_{2}=1.0$, $\delta_{13}\equiv\ln Z_{1}-\ln Z_{3}=0.54$, and $\delta_{34}\equiv\ln Z_{3}-\ln Z_{4}=1.43$. These results reveal at least two facts. On the one side, the Jeffreys’ scale $\delta_{13}$ [96, 97] reflects that current XENON-1T experiment has no significant preference of the $B_{\mu_{X}}\neq 0$ case to the $B_{\mu_{X}}=0$ case [98]. On the other side, $\delta_{12}$ and $\delta_{34}$ show that the Bayesian evidence (or equivalently the averaged $\mathcal{L}_{\rm DM}$) is reduced by a factor of more than $40\%$. It implies that a sizable portion of the parameter space will become disfavored once the future LZ experiment improves the XENON-1T’s sensitivity by 50 times. This feature is also reflected in Fig. 1 and Fig. 2 by the sizable shrink of the $1\sigma$ CIs.

In order to better understand Fig. 1, we describe how we obtained it. From Eq. (3.5), the 2D PL $\mathcal{L}(\lambda_{\nu},m_{\tilde{\nu}_{1}})$ is given by

[TABLE]

In plotting the figure, we implemented the maximization over the parameters $Y_{\nu}$, $A_{\lambda_{\nu}}$, $A_{Y_{\nu}}$, $m_{\tilde{\nu}}$, $m_{\tilde{x}}$, and $m_{\tilde{l}}$ by three steps. First, we split the $\lambda_{\nu}-m_{\tilde{\nu}_{1}}$ plane into $80\times 80$ equal boxes, i.e., we divided each dimension of the plane by 80 regular bins. Second, we fit the samples obtained in the scan into each box. Consequently, samples in each box correspond to roughly equal $\lambda_{\nu}$ and $m_{\tilde{\nu}_{1}}$, even though the other parameters may differ significantly. Finally, we select the maximum likelihood value of the samples in each box as the PL value. These procedures imply that the CIs are not necessarily contiguous, instead they usually distributed in isolated islands [90, 15]. Besides, we emphasize that $\chi^{2}_{\rm min}\simeq 0$ for the best point in the scans. This is because the DM experiments are independent and consistent with each other, and the ISS-NMSSM can explain them well.

Next, we study 2D PL on $Y_{\nu}-\lambda_{\nu}$ plane. The results are shown in Fig. 3 where the red dashed line denotes the correlation $\lambda_{\nu}\mu/(Y_{\nu}\lambda v_{u})=9.4$ or equivalently $\lambda_{\nu}=2.9Y_{\nu}$ from the unitarity constraint. This figure shows that $Y_{\nu}$ is maximized at 0.17 when $\lambda_{\nu}\simeq 0.52$ and it is upper bounded only by the unitarity. The reason is that the unitarity requires $\lambda_{\nu}\gtrsim 2.9Y_{\nu}$, so the SI cross-section is much more sensitive to $\lambda_{\nu}$ than to $Y_{\nu}$. Consequently, the DD experiments set the upper bound of $\lambda_{\nu}$ and by contrast, the unitarity limits $Y_{\nu}$.

We also plot 2D PLs on $Y_{\nu}-m_{\tilde{l}}$ and $V_{11}-m_{\tilde{l}}$ planes in Fig. 4 and Fig. 5, respectively. Fig. 4 indicates that the $2\sigma$ CI in each panel occupies a roughly rectangular area on the $Y_{\nu}-m_{\tilde{l}}$ plane. This result reflects that $\mathcal{L}_{\rm DM}$ is insensitive to parameter $m_{\tilde{l}}$. It can be understood from the following two aspects. One is that $\mathcal{L}_{\rm DM}$ relies on $m_{\tilde{l}}$ mainly through $V_{11}$ by the $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}h_{i}$ coupling in Eq. (2.27). The other is that $m_{\tilde{l}}$ and $V_{11}$ are weakly correlated, which can be inferred by the expression of $m_{12}$ and $m_{13}$ in Eq. (2.26) and is shown numerically in Fig. 5 and Fig. 10. Specifically, for the $B_{\mu_{X}}\neq 0$ case, both the annihilation and the scattering are insensitive to $V_{11}$ since its magnitude is small, and so is $\mathcal{L}_{\rm DM}$. This property determines that the allowed range of $Y_{\nu}$ is roughly independent of $m_{\tilde{l}}$, and thus explains the rectangular shape. For the $B_{\mu_{X}}=0$ case, although the effective cross-section in Eq. (2.35) is sensitive to $V_{11}$ by the formula in Eq. (2.34), the XENON-1T experiment has required $|V_{11}|\lesssim 0.02$ and this upper bound is very insensitive to $m_{\tilde{l}}$. In this case, one may replace $m_{\tilde{l}}$ by $V_{11}$ as a theoretical input so that $\mathcal{L}_{\rm DM}$ does not depend on $m_{\tilde{l}}$ any more. This feature again leads to the conclusion that the allowed range of $Y_{\nu}$ is roughly independent of $m_{\tilde{l}}$. We add that the tight experimental constraint on the mixing $V_{11}$ for the $B_{\mu_{X}}=0$ case was also discussed in [99]. We also add that one may fix $m_{\tilde{l}}$ in performing global fit of the ISS-NMSSM to experimental data due to the insensitivity of $\mathcal{L}_{\rm DM}$ to $m_{\tilde{l}}$. Such a treatment affects little the generality of the fit results.

In Table 2, we present the details of two points to illustrate the scenario’s features further. For the point P1, the DMs annihilated mainly by $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to h_{s}h_{s},h_{s}h$ to get the density. The process $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to hh$ is unimportant because $|C_{h_{s}hh}|$ is significantly smaller than $|C_{h_{s}h_{s}h_{s}}|$ and $|C_{h_{s}h_{s}h}|$, and also because the phrase space of the final state is relatively small. By contrast, the DMs got their right relic density mainly by the Higgsino pair annihilation for the point P2, and the mass splitting is $\Delta\equiv m_{\tilde{\chi}_{1}^{0}}-m_{\tilde{\nu}_{1}}\simeq 7~{}{\rm GeV}$. We confirmed that, due to the specific parameter setting of P2, there is cancellation between different contributions to the process $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to h_{s}h_{s}$, and consequently, its effect is negligibly small. Besides, both the points predict $Y_{\nu}\sim 0.01$. As a result, $V_{11}$’s magnitude is only a few thousandths, and the DM-neutron scattering rate is not much larger than the DM-proton scattering rate. We verified that, once we set $B_{\mu_{X}}=-100~{}{\rm GeV}$, the two rates became roughly equal.

3.3 Results for the massive $h_{s}$ scenario

In the massive $h_{s}$ scenario, the Higgs-mediated SI cross-section is given by

[TABLE]

when one takes the parameters in Table 1. In large $\lambda_{\nu}$ and $Y_{\nu}$ case, e.g., $\lambda_{\nu}\gtrsim 0.4$ and $Y_{\nu}\gtrsim 0.4$, the typical sizes of $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[S]}$ and $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[H_{u}^{0}]}$ are $100~{}{\rm GeV}$ and $10~{}{\rm GeV}$, respectively. Thus, the first term in the brackets is no longer more critical than the other terms, and the $\sigma^{\rm SI}_{\tilde{\nu}_{1}-N}$ for $m_{\tilde{\nu}_{1}}\simeq 300~{}{\rm GeV}$ may reach $10^{-46}~{}{\rm cm^{2}}$ only in optimal cases. Consequently, the XENON-1T experiment scarcely limit the $B_{\mu_{X}}\neq 0$ case. This situation is significantly different from the light $h_{s}$ scenario.

Similar to the analysis of the light $h_{s}$ scenario, we performed four independent scans over the parameter region in Eq. (3.7), and projected the PL onto different planes. The results are presented from Fig. 6 to Fig. 10 in a way similar to those for the light $h_{s}$ scenario. These figures indicate the following facts:

• •

Since the unitarity for the parameters in Table 1 requires only $\lambda_{\nu}\gtrsim 1.3~{}Y_{\nu}$, $Y_{\nu}$ may be comparable with $\lambda_{\nu}$ in size. As a result, the SI cross-section is sensitive to both $\lambda_{\nu}$ and $Y_{\nu}$, which is different from the light $h_{s}$ scenario.

• •

Since the $B_{\mu_{X}}\neq 0$ case is hardly limited by the XENON-1T experiment, both $\lambda_{\nu}$ and $Y_{\nu}$ may be larger than 0.4, which is shown in the upper left panel of Fig. 8. However, with the experimental sensitivity improved or the $Z$-mediated contribution considered in the $B_{\mu_{X}}=0$ case, the DM DD experiments become powerful enough to limit $\lambda_{\nu}$ and $Y_{\nu}$. In this case, $Y_{\nu}\gtrsim 0.4$ may contradict the experiments, which is indicated by the other panels of Fig. 8.

• •

$\tilde{\nu}_{1}$ obtained the correct density through the co-annihilation with $\tilde{\chi}_{1}^{0}$, which is reflected by the range of $m_{\tilde{\nu}_{1}}$ in Fig. 6. We will take the points P3 and P4 in Table 2 as examples to show more details of the annihilation later.

We confirmed that $\delta_{12}\equiv\ln Z_{1}-\ln Z_{2}=0.2$, $\delta_{13}\equiv\ln Z_{1}-\ln Z_{3}=0.62$ and $\delta_{34}\equiv\ln Z_{3}-\ln Z_{4}=0.56$ in the massive $h_{s}$ scenario. Similar to the analysis of the light $h_{s}$ scenario, the smallness of $\delta_{12}$ and $\delta_{34}$ reflects that the LZ experiment can not improve the constraint of the XENON-1T experiment on the scenario significantly, and the smallness of $\delta_{13}$ reflects that the XENON-1T experiment does not show significant preference of the $B_{\mu_{X}}\neq 0$ case to the $B_{\mu_{X}}=0$ case.

• •

Concerning the other features of the massive $h_{s}$ scenario, such as the suppression of the SI cross-section and the correlation of $m_{\tilde{l}}$ with $Y_{\nu}$ and $V_{11}$, they are similar to those of the light $h_{s}$ scenario. We do not discuss them anymore.

Next, let us study two representative points, P3 and P4, of the massive $h_{s}$ scenario in Table 2. For the former point, it is the annihilation of the Higgsino pair that is responsible for the measured density, and the corresponding mass splitting is about $5~{}{\rm GeV}$. By contrast, the $\tilde{\nu}_{1}^{0}\tilde{H}$ annihilation mainly accounts for the latter point density, and the mass splitting reaches about $19~{}{\rm GeV}$. The difference is caused by the fact that P4 takes a relatively large $Y_{\nu}$ and a lighter $m_{\tilde{l}}$, making the $\tilde{\nu}_{1}^{0}\tilde{H}$ annihilation more critical. Besides, it is notable that both the points predict $Y_{\nu}\gtrsim 0.18$ to induce a sizable $\tilde{\nu}_{L}$ component in $\tilde{\nu}_{1}$, e.g., $|V_{11}|>0.01$. Consequently, $Z$ boson can mediate a large DM-neutron scattering so that $\sigma_{\tilde{\chi}_{1}-n}^{\rm SI}\gg\sigma_{\tilde{\chi}_{1}-p}^{\rm SI}$. Such a significant difference disappears if one sets $B_{\mu_{X}}\neq 0$.

In summary, both $\lambda_{\nu}$ and $Y_{\nu}$ are more constrained in the light $h_{s}$ scenario than in the massive $h_{s}$ scenario. The unitarity always plays a vital role in limiting $Y_{\nu}$ except for the case shown in the last panel of Fig. 8, where the LZ experiment may be more critical in limiting $Y_{\nu}$. We emphasize that the tight DD constraint on the $B_{\mu_{X}}\neq 0$ case of the light $h_{s}$ scenario arises from that $h_{s}$ is light and it contains sizable doublet components. In this case, the coupling $C_{\tilde{\nu}_{1}^{\ast}\tilde{\nu}_{1}{\rm Re}[S]}$ contributes significantly to the scattering rate.

Before we end this section, we emphasize that the parameter points discussed in this work are consistent with the LHC results in searching for sparticles. Specifically, for the parameters in Table 1, it is evident that the LHC fails to detect gluinos and squarks because these particles are too massive. Concerning the Higgsino-dominated particles, they may be detectable at the 8 TeV and 13 TeV LHC since their production rates reach $100~{}{\rm fb}$. We scrutinized the property of the points in $B_{\mu_{X}}=0$ case and found that they all predict

[TABLE]

due to the Yukawa interaction $Y_{\nu}\,\hat{l}\cdot\hat{H}_{u}\,\hat{\nu}_{R}$ in the superpotential. In this case, the most promising way to explore the two scenarios at the LHC is to search the Di-$\tau$ plus missing momentum signal through the process $pp\to\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\to(\tau^{\pm}E_{\rm T}^{\rm Miss})(\tau^{\mp}E_{\rm T}^{\rm Miss})$ [14, 15]. So far, the ATLAS collaboration has finished three independent analyses of the signal based on $20.3~{}{\rm fb}$ data at the 8 TeV LHC [100], $36.1~{}{\rm fb}$ data at the 13 TeV LHC [101], and $139~{}{\rm fb}$ data at the 13 TeV LHC [102], respectively. We repeated these analyses by elaborated Monte Carlo simulations, like what we did for the first two analyses in [14, 15]. We found that the tightest constraint on the two scenarios comes from the last analysis, and its efficiency in detecting the signal decreases gradually as the gap between $m_{\tilde{\nu}_{1}}$ and $m_{\tilde{\chi}_{1}^{\pm}}$ becomes narrow. As far as the light and massive $h_{s}$ scenarios are concerned, the analysis can not exclude at $95\%$ confidence level the points satisfying $m_{\tilde{\nu}_{1}}\gtrsim 100~{}{\rm GeV}$ and $m_{\tilde{\nu}_{1}}\gtrsim 200~{}{\rm GeV}$, respectively. So we conclude that the LHC analyses do not affect the results presented in this work.

4 Conclusion

Motivated by the increasingly tight limitation of the DM DD experiments on the traditional neutralino DM in the natural MSSM and NMSSM, we extended the NMSSM by the inverse seesaw mechanism to generate the neutrino mass in our previous studies [14, 16, 8], and studied the feasibility that the lightest sneutrino acts as a DM candidate. A remarkable conclusion for the theory is that experimental constraints from both the collider and DM search experiments are relaxed significantly. Consequently, large parameter space of the NMSSM that has been experimentally excluded resurrects as physical points in the extended theory. In particular, the higgsino mass may be around 100 GeV to predict Z-boson mass naturally. This feature makes the extension attractive and worthy of a careful study.

We realized that sizable neutrino Yukawa couplings $\lambda_{\nu}$ and $Y_{\nu}$ contributed significantly to the DM-nucleon scattering rate. Thus, the recent XENON-1T experiment could limit them. We also realized that the unitarity in the neutrino sector set a specific correlation between the couplings $\lambda_{\nu}$ and $Y_{\nu}$, which in return limited the parameter space of the ISS-NMSSM. Since these issues were not studied before, we investigated the impact of the leptonic unitarity and current and future DM DD experiments on the sneutrino DM sector in this work. Specially, we considered the light and massive $h_{s}$ scenarios after noticing that the singlet dominated Higgs plays a vital role in both the DM annihilation and the DM-nucleon scattering. For each scenario, we studied the $B_{\mu_{X}}\neq 0$ and $B_{\mu_{X}}=0$ case separately. Their difference comes from that Z boson can mediate the DM-nucleon scattering for the $B_{\mu_{X}}=0$ case, and thus, the experimental constraints on it are much tighter.

In this work, we encoded the experimental constraints in a likelihood function and performed sophisticated scans over the vast parameter space of the model by the Nested Sampling method. The results of our study are summarized as follows:

• •

The XENON-1T experiment set an upper bound on the couplings $\lambda_{\nu}$ and $Y_{\nu}$, and the future LZ experiment will improve the bound significantly. The limitation is powerful when $h_{s}$ is light and contains sizable doublet components.

• •

As an useful complement to the DM DD experiments, the unitarity always plays a vital role in limiting $Y_{\nu}$. It becomes more and more powerful when $v_{s}$ approaches $v$ from top to bottom.

• •

The parameter space favored by the DM experiments shows a weak dependence on the left-handed slepton soft mass $m_{\tilde{l}}$. This property implies that one may fix $m_{\tilde{l}}$ in surveying the phenomenology of the ISS-NMSSM by scanning intensively its parameters and considering various experimental constraints. This treatment does not affect the comprehensiveness of the results.

• •

The DM experiments tightly limit the left-handed sneutrino component in the sneutrino DM, e.g., if one considers the XENON-1T experiment’s results, $|V_{11}|\lesssim 0.15$ for the $B_{\mu_{X}}\neq 0$ case and $|V_{11}|\lesssim 0.02$ for the $B_{\mu_{X}}=0$ case; these upper bounds become 0.10 and 0.01, respectively, if one adopts the LZ experiment’s sensitivity.

Finally, we briefly discuss the phenomenology of the ISS-NMSSM. The sparticles’s signal in this theory may be distinct from those in traditional supersymmetric theories, and so is the strategy to look for them at the LHC. This feature can be understood as follows: since the sneutrino DM carries a lepton number, and in most cases has feeble interactions with particles other than the singlet-dominated Higgs boson and the massive neutrinos, the sparticle’s decay chain is usually long, and its final state contains at least one $\tau$ or $\nu_{\tau}$. In addition, the decay branching ratio depends not only on particle mass spectrum but also on new Higgs couplings, such as $Y_{\nu}$ and $\lambda_{\nu}$. As a result, sparticle’s phenomenology is quite complicated [8, 16]. Depending on the mechanism by which the DM obtained the correct density, one usually encounters the following two situations:

• •

The DM co-annihilated with the Higgsino-dominated particles. This situation requires the mass splitting $\Delta\equiv m_{\tilde{\chi}_{1}^{0}}-m_{\tilde{\nu}_{1}}$ to be less than about $10~{}{\rm GeV}$. Consequently, the Higgsino-dominated particles usually appear as missing momentum at the LHC due to the roughly degenerate mass spectrum. As pointed out in [8], this situation’s phenomenology may mimic that of the NMSSM with the Higgsino-dominated $\tilde{\chi}_{1}^{0}$ as a DM candidate.

• •

The singlet-dominated particles $\tilde{\nu}_{1}$, $h_{s}$, $A_{s}$, and $\nu_{h}$ compose a secluded DM sector where the DM was mainly annihilated by any of the channels $\tilde{\nu}_{1}\tilde{\nu}_{1}^{\ast}\to A_{s}A_{s},h_{s}h_{s},\nu_{h}\bar{\nu}_{h}$. It communicates with the SM sector by the Higgs-portal or the neutrino-portal. As we introduced before, this situation constrains the Yukawa coupling $\lambda_{\nu}$ tightly in getting the measured density, but it has no limitation on the splitting between $m_{\tilde{\nu}_{1}}$ and the Higgsino mass. As mentioned before, the signals of the sparticles in this situation are complicated. However, systematic researches on this subject are still absent.

We suggest experimentalists to look for the $2\tau$ plus missing momentum signal of the process $pp\to\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\to(\tau^{\pm}E_{\rm T}^{\rm Miss})(\tau^{\mp}E_{\rm T}^{\rm Miss})$ in testing the theory. Unlike the colored sparticles that may be very massive, light Higgsinos are favored by natural electroweak symmetry breaking. As a result, they are expected to be richly produced at the LHC. For the secluded DM case, ATLAS analyses have excluded some parameter space discussed at the end of the last section. With the advent of the LHC’s high luminosity phase, more parameter space will be explored. For example, we once compared the ATLAS analyses of the signal at the 13 TeV LHC with $36.1~{}{\rm fb^{-1}}$ and $139~{}{\rm fb^{-1}}$ data [101, 102]. We found the excluded region on $m_{\tilde{\nu}_{1}}-m_{\tilde{\chi}_{1}^{\pm}}$ plane expanded from $m_{\tilde{\chi}_{1}^{0}}\lesssim 45~{}{\rm GeV}$ to $m_{\tilde{\chi}_{1}^{0}}\lesssim 110~{}{\rm GeV}$ for $m_{\tilde{\chi}_{1}^{\pm}}=200~{}{\rm GeV}$, and from $m_{\tilde{\chi}_{1}^{0}}\lesssim 120~{}{\rm GeV}$ to $m_{\tilde{\chi}_{1}^{0}}\lesssim 200~{}{\rm GeV}$ for $m_{\tilde{\chi}_{1}^{\pm}}=300~{}{\rm GeV}$. Concerning the co-annihilation case, it is hard for the LHC to detect the signal due to the compressed spectrum, but the future International Linear Collider may be capable of doing such a job (see, for example, the study in [103] for the compressed spectrum case). We emphasize that, different from the prediction of the MSSM, $m_{\tilde{\chi}_{1}^{\pm}}$ may be significantly larger than $m_{\tilde{\chi}_{1}^{0}}$ in the ISS-NMSSM due to the mixing of $\tilde{H}_{u,d}$ with $\tilde{S}$ in Eq. (2.5). As a result, the splitting between $m_{\tilde{\chi}_{1}^{\pm}}$ and $m_{\tilde{\nu}_{1}}$ can reach $20~{}{\rm GeV}$ (see the points in Table 2), and it becomes even larger as the parameter $\lambda$ increases. This feature is beneficial for the signal’s detection.

Acknowledgement

This work is supported by the National Natural Science Foundation of China (NNSFC) under grant No. 11575053 and 12075076.