Dark matter constraints and the neutralino sector of the scNMSSM

The neutralino sector of the semi-constrained next-to-minimal supersymmetric standard model is explored under recent experimental constraints, with special attention to Dark Matter limits. The effects of Dark Matter relic density upper and lower bounds, and recent direct detection constraints on spin-independent and dependent cross-sections, are thoroughly analyzed. Particularly, we show which regions of the parameter space are ruled-out due to the different Dark Matter constraints and the corresponding model-specific parameters: $\lambda, \kappa, A_{\lambda}$ and $A_{\kappa}$. Additionally, we analyze all annihilation and co-annihilation processes (with heavier neutralinos and charginos) that contribute to the dark matter relic density. The mass components of the dark matter candidate, the lightest neutralino $\tilde{\chi}_1^0$, are studied, and the decays of heavy neutralinos and charginos, especially $\tilde{\chi}_2^0$ and $\tilde{\chi}_1^+$, into the lightest neutralino are probed. We impose semi-universal boundary conditions at the GUT scale, and require a moderate range of $\tan{\beta} \lesssim 10$. We find that the allowed parameter space is associated with a heavy mass spectrum in general and that the lightest neutralino is mostly Higgsino with a mass range that resides mostly between 1000 GeV and 1500 GeV.


Introduction
The nature of Dark Matter (DM) continues to be a mystery, even though there is an abundance of evidence for its existence [1][2][3][4][5][6]. A well-motivated candidate for DM is the neutralino LSP in SUSY models [7]. However, such models are being challenged by recent null results from experiments, including DM searches [8][9][10][11][12][13]. The latter can be categorized into collider, direct detection, and indirect detection searches. Collider searchers rely on detecting missing energy events that could arise from decays into DM particles, indirect detection depends on observations of astronomical events, such as gamma ray lines, and direct detection searches rely on observing possible elastic scattering events between DM particles and atomic nuclei (a comprehensive review is given in [14]). One of the most considered low-scale SUSY models is the NMSSM (reviewed in [15,16]), which opens new possibilities given that its neutralino sector is richer than that in the MSSM. A fact that is attributed to the inclusion of a SM-Singlet superfield, which is absent in the MSSM. The inclusion of such a singlet is motivated by explaining the origin of the µ-term in the MSSM and increasing the predicted SM Higgs mass at tree-level.
The issue of DM in the NMSSM has been and still is the center of attention of researchers in the field [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. The case of low-mass neutralino LSP has been investigated with attention to the interesting cases large singlino-higgsino mixture, or singlino-like DM [23,[33][34][35][36]. Additionally, the implications of direct and indirect detection for the NMSSM and its extensions were detailed in [22,[37][38][39][40][41][42][43]. Most recently, the authors of [44] studied the NMSSM and DM, where low-scale input parameters were used, and the attention is targeted towards future searches. However, most of these works do not include the most recent constraints from DM searches, especially those coming from direct detection. Moreover, from a phenomenology point of view, it is interesting to understand in-details how current limits affect the NMSSM with the well-motivated GUT boundary conditions, in terms of the affected regions of the parameters, and more physically, which channels contribute to DM relic density, the branching ratio of chargino/neutralino to the LSP.
Considering all this, we explore, exploiting state-of-art tools incorporating latest DM constraints, the parameter space of the NMSSM under semi-constrained boundary conditions at the GUT scale, given that a major motivation for SUSY models is to provide gauge coupling unification. Particularly, we analyze how relic density abundance constraints, as well as recent direct detection limits affect the model, without particular focus on specific regions except for relatively moderate tan β 10, which motivated by enhancing the tree-level prediction of the SM-like Higgs mass in the model. We present the portions of the parameter space that pass all known constraints (apart from g-2 of the muon), and show those regions which fail the aforementioned DM constraints, providing a deeper insight into the model that usually overlooked in recent literature.
The paper is organized as follows. In Section 2, a brief overview of the superpotential of the NMSSM and the neutralino mass matrix is provided, and this is where the necessary parameters of the model are introduced. Next, Section 3 specifies the boundary conditions of model, and presents our methodology of probing the parameter space of the model, and the relevant non-DM and DM constraints taken into account. In Section 4, the results are presented and divided into subsections corresponding to the various aspects of the neutralino section that are considered as mentioned earlier. Section 5 is where we discuss the results, comparing them to recent literature, and conclude the paper.

The NMSSM and its Neutralino Sector
The NMSSM is characterized by the following superpotential, where W¡ µ MSSM encodes the Yukawa couplings as in the MSSM. The second term replaces the µ-term in the MSSM so that it is dynamically generated when the singlet superfieldŜ acquires a vacuum expectation value (VEV). λ is the coupling betweenŜ and the up and down Higgs superfieldsĤ u ,Ĥ d . The third term represents the cubic self-coupling ofŜ. This κ term is included to maintain a Z 3 invariance in the model.
In contrast to the MSSM, the soft-SUSY breaking lagrangian of the NMSSM contains extra terms associated with trilinears terms A λ and A κ , as well as a soft-SUSY breaking mass term for the singlet. The Higgs sector of the model receives extra contributions from the singlet, and the number of Higgs bosons is increased to seven particles (two charged, 3 CP-even, and 2 CP-odd). Furthermore, the neutralino sector is also affected by the presence of the new NMSSM-specific terms. As opposed to the 4 × 4 mass matrix in the MSSM, the neutralino mass matrix in the NMSSM is 5 × 5, which reads, where M 1,2 represent gaugino masses, while µ = λs is the effective mu-term. g 1,2 denote gauge couplings, v u,d , s indicate VEVs of the up-and down-Higgs fields, and singlet field. The rest of parameters were defined previously. Upon diagonalization, the mass eigenstates are ordered from small to large asχ 0 1 , . . . ,χ 0 5 , where the first is the lightest neutralino, which is usually the LSP of the model and the dark matter candidate.
Finally, the NMSSM can be specified at the GUT scale by a number of input parameters, which are: a soft-SUSY breaking scalar mass parameter (m 0 ), a soft-SUSY breaking gaugino mass parameter (m 1/2 , a soft-SUSY breaking trilinear parameter A 0 , the ratio between the VEVs of up-and down-Higgs fields tan β = v u /v d , an effective Higgs/Higgsino mass parameter µ eff = λs, a singlet-doublet coupling parameter λ, and a singlet coupling parameter κ. It is desirable to understand how NMSSM-specific parameters behave and affect the model's predictions, and therefore, it is often the case that one treats the two trilinear soft-SUSY breaking parameters A λ and A κ differently from A 0 .

The parameter space
NMSSMTools (v5.5.2) was used to generate the mass spectrum of the model. NMSSMTools utilizes MicrOmega 2 , which in turn, calculates the relevant DM observables of the model. All theoretical and experimental constraints where taken into account apart from the (g − 2) µ anomaly. These constraints include limits on the decays of both the K and B mesons, on the SM Higgs mass, couplings, and decays (including exotic decays), on supersymmetric particles. Theoretical constraints include reaching a global minimum, no negative physical mass values, successful EWSB in the Higgs sector.
Moreover, the package was modified to our purpose, which is to analyze each DM constraint separately (after imposing non-DM constraints). This enable us to understand and report on how DM limits affects the considered model. However, and as far as DM constraints are concerned, in this paper, we have restricted our analysis to constraints on DM relic density and direct detection.
With this in mind, we explore the parameter space focusing on specific regions where tan β 10, but allow λ to vary as shown in Table 1. We specify the inputs of the model at the GUT scale in what is commonly called the scNMSSM as explained in the previous Section. The ranges of the input parameters are shown in Table 1. Table 1 It is important to note that the ranges for the trilinear parameters are varied separately, and the same is true for λ and κ. The latter parameters are specified in the low-scale, and m 2 hi (i = u, d) is calculated by NMSSMTools at the GUT-scale given the specified input parameters (for more details about how NMSSMTools computes the spectrum we refer the reader to the relevant literature [45][46][47][48][49]).
Finally, it is worth mentioning that the latest version of NMSSMTools utilizes MicrOMEGAsv5.0, and includes a number of recent and important constraints on Higgs physics [50][51][52], as well as DM limits on spin independent crosssections from DarkSide-50 [53] and CRESST-II [54], and limits on spin dependent cross-sections from XENON1T [55] and PICO60 [56]. In Fig. 1, it can be seen in the top-left panel that the ruled-out points are associated with no DM candidate in the model. Basically, this reflects results where the LSP is a charged particle (hence cannot be a DM candidate), or that the input parameters did not lead to successful solutions causing the algorithm of MicrOMEGAs to stop. This issue becomes more relevant as both A λ and A κ increase, especially for A λ > −2000 GeV. On the other hand, the top-right panel shows the ruled-out points which predict too large values of the relic density Ωh 2 0.131, while the middle-left panel shows red points with Ωh 2 0.107. Comparing the two plots indicates that in the scanned parameter space, most the problematic points (passing all non-DM constraints still) are associated with too small values for DM relic density 3 . As for cross-sections, the middle-left panel presents models predicting a DM direct detection rate that is too large. That is a larger than experimentally allowed spin-independent cross-section σ SI , which resembles the density of ruled-out points due to Ωh 2 being too small in the previously described plot. This shows that both of these limits play a major effect in constraining the scNMSSM, and it warns against common practice of imposing the upper limit on Ωh 2 only, since even if one accepts a point with too small relic density, such a point might be ruled-out by direct detection constraints. Next, the ruled-out points due to spin-dependent cross-sections, σ n SD and σ p SD being larger than the experimental limits are shown in the bottom-left panel and bottom-right panel, respectively. Clearly, such constraints have smaller impact, but can eliminate parts of the model where both A λ and A κ are positive. In all the previous sub-figures, it is clear that DM constraints heavily affect regions where both trilinears are positive and large. Lastly, and focusing on the allowed points (shown in green), the Figure shows that when A λ < 0, only negative A κ can lead to successful predictions. However, as A λ becomes positive, A κ can take positive values.

Mappiong DM constraints into the models' parameters
To facilitate comparison and acquire a deeper insight, the same results are presented in the λ-κ plane, as is shown in Fig. 2. Some of the data passing all non-DM restrictions still failed to produce a successful DM candidate (top-left panel), this took place particularly in the region 0.4 < λ < 0.7 and 0.3 < κ < 0.6. Other points failed to comply with the upper limit on Ωh 2 (top-right panel), which occurs in all of the parameter space, but as clearly shown, regions where λ 0.55 with κ 0.2 suffer from this problem as no allowed points can be found there (given our choices of inputs). The middle-left panel shows the case where Ωh 2 is too small. Comparing this result with that of too large σ SI , it is clear that these two problems have the main impact as stated before. This impact is reflected in the λ-κ plane for large values of λ accompanied with small values of κ. In the same region, we can see in both the bottom-left, and -right panels that the limits of spin-dependent cross-sections affect the same "bad" region.
Finally, Fig. 3 shows the results in the common (m 0 -m 1/2 ) plane, which in turn provides further understanding of the how DM constraints impact the parameter space of the model, and enables comparison with other low-scale SUSY models such as the MSSM. As can be seen in the top-left panel, the issue of not obtaining a successful DM candidate is concentrated in regions where m 0 1200 GeV. On the other hand, most of the ruled-out points due to predicting large Ωh 2 are associated with m 1/2 < 3000 GeV (top-right panel). Moreover, the two main problems causing the ruling-out of the red points, that is small Ωh 2 and large σ SI , are common in the whole parameter space (see middle panels), whereas the limits on the spin-dependent cross-sections are mainly present in the region where m 1/2 < 3000 GeV. However, these are somewhat less relevant than the previous two constraints which occur in the same area. These results also show that the good points (in green) are concentrated in regions where m 1/2 > 3000 GeV. Such a large value for m 1/2 , which in turn is needed to satisfy other non-DM limits, indicates a heavy mass spectrum, which is certainly the case in the studied parameter space.

Properties of the mass of the neutralino DM
Now that we have an understanding of how recent DM constraints affects the parameters of the model. It is crucial to analyze the DM particle in terms of its mass content. In the allowed parameter space, the LSP represents a valid DM candidate, and has a mass range that is mostly concentrated between 1000 GeV and 1200 GeV. Lighter and heavier values still exist in our results, but such cases are sparse. 1 ≤ 1149 GeV. The contribution corresponding to the lower mass is about 0.98, while that associated with the higher mass is roughly 0.97. The maximum Singlino contribution is found to be approximately equal to 0.99 at mχ0 1 ≈ 1106 GeV. In contrast, the minimum contribution is found to be roughly 0.92 at mχ0 1 ≈ 912 GeV. Meanwhile, the region with a lower Singlino contribution occurs between 691 ≤ mχ0 1 ≤ 1526 GeV, where the Singlino component is found to be about 0.0026, and 0.022 for the lower and higher mass, respectively. The maximum contribution in this region is about 0.27 at mχ0 1 ≈ 1200 GeV. However, most of the data points lie between 1025 ≤ mχ0 1 ≤ 1290 GeV and cover a range of Singlino contributions between ∼ 0.0002 − 0.06, which is a negligible contribution. Hence, the scanned parameter space, it is very difficult to find a Singlino-like DM candidate.    The Bino components are given in the bottom-right panel of Fig. 4, where it is clear that the data points split into two regions with higher and lower contribution to the LSP mass. The region with large Bino contribution occurs at 451 ≤ mχ0 1 ≤ 1294 GeV, where the value of the component is 0.99 (for the small end of the mass) and 0.98 (for the large end of the mass). The maximum Bino contribution is almost unity and occurs at mχ0 1 = 923, 1015 GeV. However, as can be seen in the Figure, most of the allowed points lie outside of this high Bino contribution region. Indeed, the region with low Bino contribution occurs at 313 ≤ mχ0 1 ≤ 1526 GeV, where the component takes the values 0.005, and 0.053, respectively. The maximum Bino contribution in this "low" region is about 0.2, which lies at mχ0 1 = 1431 GeV. The area where the allowed data is concentrated is the mass range 1032 ≤ mχ0 1 ≤ 1346 GeV, where the Bino contribution assumes values roughly between 0.03 − 0.12.
Finally, the results above show that in the model under consideration, the LSP is mostly Higgsino. The bino and singlino components are small and somewhat comparable in most of the parameter space, expect very rarely where both can be close to unity (in such rare cases, the Higgsino component is very small). The wino component is always negligible.

Contributions to dark matter relic density
There are several annihilation and co-annihilation processes that contribute to the dark matter relic density. Here, we present the results for those processes that can contribute more than 1%. In the relevant plots, all of the points with zero y-axis correspond to data points that have a contribution of less than 1%. The exact values were not included in the output routine since these are very small (and hence there were set to zero in the plotting command), but those values are indeed taken into account in MicrOMEGAs and NMSSMTools for the calculation of the thermally averaged cross-section and the corresponding relic density. In addition, the co-annihilation modesχ 0 1χ + 1 → ud, cs are illustrated in the right panel of the second row, and the left panel of the third row, respectively. Here, the data points spread between a lower value of mχ0 1 = 913 GeV, corresponding to the contribution 1.1%, and a higher value of mχ0 1 = 1528 GeV, corresponding to the contribution 5%. The maximum contribution is 9.4% which corresponds to mχ0 1 = 1019 GeV, the minimum contribution is 1.1% corresponding to mχ0 1 = 913 GeV. Most of the allowed points are located in the range mχ0 1 = 1029 − 1330 GeV, with a contribution range of 9 − 7%.
The co-annihilation modeχ 0 1χ + 1 → tb is given in the right panel of the third row. The allowed points spread between a lower mχ0 1 = 580 GeV, corresponding to the contributions 1.45%, and a higher value mχ0 1 = 1526 GeV, corresponding to the contribution 3.46%. The maximum contribution is 8.36% which corresponds to mχ0 1 = 914 GeV, while the minimum contribution is 1.01% corresponding to mχ0 1 = 1195 GeV. Most of the data points are located at mχ0 1 = 1023 − 1400 GeV with a contribution range 7.31 − 3.41%. Furthermore, the bottom-left panel shows the co-annihilation channelχ 0 1χ + 1 → H 2 W + . In this case, the data points spread between mχ0 1 = 914, GeV corresponding to the contribution 3.62%, and mχ0 1 = 1514 GeV, corresponding to the contribution 1.66%. The maximum contribution is 6.63% which corresponds to mχ0 1 = 1114 GeV. The majority of data points are located at mχ0 1 = 1069 − 1324 GeV, with a contribution range 4.59 − 0.998%. The last co-annihilation mode considered here isχ 0 1χ + 1 → A 1 W + , which is shown in the bottom-right panel. The data points spread between mχ0 1 = 952 GeV, corresponding to the contribution 6.35%, and mχ0 1 = 1526 GeV, corresponding to the contribution 3.1%. The maximum contribution is 7.1% associated with mχ0 1 = 1317 GeV. Most data points are located at mχ0 1 = 1043 − 1338 GeV, with a contribution range 6.82 − 0.97%.
4.4χ 0 2 andχ ± 1 decays to the LSP There are several decay channels producing an LSP (along with other particles), however, the decays of the chargino and the second-lightest neutralino are of special importance. Table 3, in the Appendix, lists all available branching ratios for a point in the allowed parameter space withχ 0 1 ≈ 1100 GeV. Fig. 6 shows that forχ 0 2 decays, the radiative decayχ 0 2 →χ 0 1 γ can be significant in regions of the parameter space, where the branching ratio can reach up to ∼ 0.8 for a nearly degenerate mass of ∼ 1200 GeV for bothχ 0 2 andχ 0 1 . On the other hand, an important type of decay channel forχ + 1 turns out to beχ + 1 →χ 0 1 ud/cs, where the branching ratio lies in the range 0.35 ≤ BR(χ + 1 →χ 0 1 ud) ≤ 0.38 for almost all of the allowed parameter space. Other decay modes of both considered sparticles into the LSP can be important as well (see Table 3 in the Appendix). For instance, the decayχ 0 2 →χ 0 1 Z can be important but in a very rare cases where the branching ratio is equal to 1 (where difference in mass between both neutralinos is around 100 GeV), which does not represent the region where the majority of allowed points lie. Additionally, decays ofχ 0 2 to a pair of SM qq and LSP can be significant, while Figure 5: Contributions of annihilation and co-annihilation processes to dark matter relic density.
weakly mediated decays ofχ + 1 to an LSP and a pair of ff (where such a pair can consist of a lepton/antilepton and an antineutrino/neutrino) can also be important with branching ratio above 0.3. Moreover, it is found that for a mass range of 1154 ≤ mχ0 2 ≤ 1517 GeV, the decay channelχ 0 2 →χ 0 1 H 1 has a branching ratio that is close to unity 0.94 ≤ BR(χ 0 2 →χ 0 1 H 1 ) ≤ 1 for a small number of successful points (where the mass difference between the two neutralinos is 300 GeV minimum). The results of the branching ratios forχ + 1 →χ 0 1 W + indicate that the branching ratio equal to 1 mostly in the range mχ+ 1 = 1151 − 1513 GeV for a limited number of successful points (where the mass difference between the chargino and the neutralino is 300 GeV minimum). The lower value BR(χ + 1 →χ 0 1 W + ) = 0.033 corresponding to mχ+ 1 = 849 GeV.

Discussion and conclusions
The results in the previous Section show that the NMSSM with semi-constrained GUT boundary conditions can explain current experimental observations including DM. However, the mass spectrum of the new particles, including the DM candidate (the LSP) is quite heavy. As far as the input parameters are concerned, it is observed that large values of m 0 and m 1/2 are required. As such inputs are increased, the whole SUSY spectrum increases, including the LSP dark matter. There are several reasons for this situation, the most crucial one is satisfying the bounds from direct searchers for SUSY particles. One can still look for regions where the parameters of the model conspire to give a lighter spectrum (e.g. regions with specific relations between the gluino mass mg or the stop mass mt 1 and mχ0 1 ), but we do not consider these fine-tuned cases in this paper.
Additionally, the effects of DM constraints on other NMSSM-specific parameters were considered. These parameters include λ and κ, where it is found that they can assume values between 0.1 and below 0.7, however, the majority of allowed points reside in regions where λ 0.6 and κ 0.3. Furthermore, the other NMSSM-specific parameters A λ and A κ , and for most of the allowed parameter space, are found to lies mostly in regions where A κ is negative.
In meanwhile, it is observed that the mass content of the DM particle (the neutralino LSP) comes mostly from the Higgsino component, where in most of the parameter space it that component assumes the value 0.7, for a range of neutralino mass between 1000 GeV and 1500 GeV. The singlino component is found to be quite small (between 0 and 0.1) in the majority of the parameter space. However, this does not entirely rule-out the option of finding regions where the singlino component could be large. Indeed, it is noticed that very rare points in the scanned parameter space could achieve such large values, but such a case requires a dedicated and separate study that is beyond the scope of this paper. On the other hand, the Bino component is found to be somewhat small with values reaching up to 0.15 in most of the allowed parameter space. Again, very few points are found to have a large Bino components, and investigating such regions is also an interesting question for future research. The last component, the Wino, is found to be small throughout the allowed parameter space.
Additionally, it is observed that in most of the allowed parameter space, the annihilation of the LSP, as a mechanism for generating the required DM relic density, is not sufficient in the model. To achieve the observed value, co-annihilation must be taken into account. Indeed, it is found that co-annihilation processes with the charginoχ + 1 constitute a major contributor to the dark matter relic density (given the specific allowed parameter space). Specifically, co-annihilation processes leading to ud/cs/tt, H 2 W + and A 1 W + are found to provide most of the contributions.
As for decays of heavier neutralinos and charginos to the LSP, we see that the decays ofχ 0 2 andχ + 1 are crucial given that they can be close to the mass ofχ 0 1 . Indeed, it is observed that the radiative decay ofχ 0 2 can be very significant accounting for about 80% of its decays, especially in regions where it has a nearly degenerate mass withχ 0 1 . which happens to be the region where most allowed points are found and the common mass between the two is about 1200 GeV. Other channels can be significant, albeit very rarely as we discussed in the results Section.