Multi-scalar signature of self-interacting dark matter in the NMSSM and beyond

This work studies the self-interacting dark matter (SIDM) scenario in the general NMSSM and beyond, where the dark matter is a Majorana fermion and the force mediator is a scalar boson. An improved analytical expression for the dark matter (DM) self-interacting cross section which takes into account the Born level effects is proposed. Due to the large couplings and light mediator in SIDM scenario, the DM/mediator will go through multiple branchings if they are produced with high energy. Based on the Monte Carlo simulation of the showers in the DM sector, we obtain the multiplicities and the spectra of the DM/mediator from the Higgsino production and decay at the LHC for our benchmark points.


Introduction
The weakly interacting massive particle being the dark matter (DM) candidate can successfully explain the large-scale structure of the universe from the galactic scale to the cosmological scale. However, there have been tensions between the N-body simulations of collisionless cold DM and astrophysical observations on small scale structure of the universe, dubbed small scale structure problem [1]. The issues can be resolved if the DM has strong self-interactions [2,3]. On the other hand, the self-interacting DM (SIDM) scenario is stringently constrained in high-velocity system such as galaxy clusters [4][5][6]. A viable SIDM scenario requires the DM self-scattering cross section that is suppressed at higher velocity and increases toward smaller velocity. Such feature can be implemented if the DM self-interaction is mediated by a light (O(1) MeV) scalar or vector particle [6][7][8][9][10][11][12][13][14][15][16][17].
The strongly interacting DM with scalar mediator can be naturally realized in the nextto-minimal supersymmetric standard model (NMSSM) [18,19]. In the small λ limit, the singlet sector (includes a CP-even scalar, a CP-odd pseudoscalar and a singlino) is nearly decoupled from the SM sector, such that they can be light while evade all experimental searches. The Peccei-Quinn limit of the NMSSM with light singlet sector has been studied in Ref. [20]. However, the SIDM scenario requires large κ coupling. In the Z 3 invariant NMSSM, the singlino mass is proportional to κ/λ, which is sizable in the limit of λ κ. The general NMSSM is required to address the SIDM scenario [21][22][23][24][25]. Meanwhile, the correct DM relic density can be set via thermal freeze-out with DM annihilating into scalar mediators.
Given light mediator and DM, together with relatively large couplings in the SIDM scenario, the production of the mediator/DM at the electroweak scale or higher will be followed by copious emissions of the mediator in the collinear direction due to the logarithmic enhancement. This is analogous to the QED/QCD parton showering. The phenomenology of the DM showering in models with extra dark gauge group has been studied in Ref. [26][27][28][29][30], where the force mediators are massless dark photon/gluons (similar to the hidden valleys scenario [31]). In the NMSSM, the mediator is the CP-even massive scalar and the DM is the massive Majorana fermion. The properties (divergence behaviors) of the splitting function are quite different from the massless vector mediator ones. More detailed studies of the DM showering with massive states have been performed in Ref. [32,33], inspired by the studies for the electroweak showing [34,35] In this paper, we focus on the SIDM scenario in the general NMSSM. An improved analytical estimation for the DM self-interacting cross section which takes into account the Born level effects is proposed. We will demonstrate that the analytical expression matches well with the numerical solution, in terms of calculations on the scanned points in the NMSSM. Two benchmark points (in the NMSSM) that address the small scale structure problem with correct relic density of SIDM and satisfy other phenomenological constraints will be provided for the first time. They are featured by light Higgsinos and large probability of φ → φφ splitting. The multiplicity and the spectrum of the scalar mediator in the Higgsino production (and decay) process will be studied in detail, based on Monte Carlo simulation of the showers in the DM sector. However, κ can not be large in the NMSSM in order to implement correct DM relic density, which suppresses the DM splitting. Two benchmark points beyond the NMSSM with κ = 2.5 and address the small scale structure issue are proposed. The DM splitting (χ → φχ) becomes also significant in this case. A generic issue of SIDM scenario with simple DM sector is the tension between Big Bang Nucleosynthesis (BBN) limits and dark matter direct detection constraints [22,[36][37][38][39]. Possible solutions require extended particle content and couplings. And the collider signature of SIDM scenario will be highly depend on the forms of extensions.
The rest of the paper is organized as follows. In Sec. 2, we will discuss how to calculate the non-relativistic DM self-interaction cross section with analytic method and numerical method. In particular, an improved analytical estimation for the DM self-interacting cross section is proposed. In Sec. 3, focusing on the general NMSSM model, we study the parameter space that addresses the small scale structure problem and is consistent with other phenomenological constraints. The showering of the DM and the mediator is discussed in Sec. 4. Finally, in Sec. 5, we discuss the possible LHC signatures of the SIDM in the NMSSM and beyond, based on four selected benchmark points.

Self-interacting dark matter -scalar mediator
The scattering between DMs (χ) with a scalar mediator (φ) in the non-relativistic limit is controlled by an attractive Yukawa potential where α χ = κ 2 /(2π) and √ 2κ is the coupling of φχχ. The scattering amplitude is where δ l is the phase shift for a partial wave l. It can be obtained by solving the Schrödinger equation for given potential V (r), and k = m χ v/2 with v being the relative velocity between the DMs in the scattering.
To describe the scattering between distinguishable particles, the transfer cross section σ T and the viscosity cross section σ V are usually used [40], which are defined as The relation between the differential scattering cross section dσ/dΩ and the phase shift δ l is given by (2 + 1)e iδ P (cos θ) sin δ 2 .

(2.4)
So it can be calculated that When describing the scattering between identical particles, only the viscosity cross section σ V is useful [41]. For Majorana fermion DM scattering, the spatial wave function should be symmetric (antisymmetric) when the total spin of the DM pair is 0 (1). Thus, σ V is replaced by σ VS and σ VA in the symmetric case and antisymmetric case, respectively, which are (2.8) Note that a symmetry factor 1/2 is inserted in integrals to avoid the double-counting in scattering of two identical particles. In the following analysis, we assume the DMs participating the scattering are unpolarized and refer to σ V as the one averaging over the all spins: To solve the Schrödinger equation, it is useful to define the new variables: In general, there is no analytic solution to the Schrödinger equation with the Yukawa potential on the a − b plane. However, in the Born regime where b 1, computed perturbatively in α χ , the scattering amplitude can be expressed as at the leading order, from which we can calculate the leading order DM scattering cross sections as where t = ab. Besides, in the Born regime, we can also use [42] e iδ l sin to estimate the δ l at the leading order, where j l is the spherical Bessel function.
In the quantum regime where t 1, the s-wave scattering is dominant, i.e. |δ 0 | |δ l | for l > 0. By taking the Hulthén approximation, δ 0 for the attractive Yukawa potential is given by [13] δ Hulthén where c ≈ b/1.6. So we can use [13,41]  as approximate expressions in the quantum regime.
Recent study in Ref. [41] provides the analytic approximations of σ T , σ VS , and σ VA for both attractive and repulsive Yukawa potentials in the semi-classical regime where t 1. The cross sections in this regime are strongly depend on β = 1/(2a 2 b). For the attractive Yukawa potential, the expressions for cross sections are summarized as where n = 1 2 n = 3 2 for σ Clas VS σ Clas VA and K i (i = 0, 1, 2) stands for the modified Bessel functions of the second kind.
The σ T and the spin averaged σ V on the whole a−b plane can be obtained by combining the expressions above. They are  The transfer cross section σ T (×k 2 /(4π)), the symmetric viscosity cross section σ VS (×k 2 /(4π)) and the asymmetric viscosity cross section σ VA (×k 2 /(4π)) on the a-b plane, calculated by Eq. 2.25 (upper panels) and calculated by numerically solving the Schrödinger equation (low panels).
Meanwhile, following the method proposed in Ref. [13], the Schrödinger equation with Yukawa potential can be solved numerically. In the lower panels of Fig. 1, we show the σ T k 2 /(4π), σ VS k 2 /(4π), and σ VA k 2 /(4π) obtained with the numerical calculation. Comparing the analytical and numerical results, we can observe that the analytical estimations match reasonably well with those from numerical calculations. A more specific comparison based on the scanned parameter points in the NMSSM will be given later in this work. As for the asymmetric viscosity cross section, the analytical calculation in the quantum regime raises zero, which is in consistent with the numerical calculation as most points in the quantum regime give σ VA k 2 /(4π) 10 −8 and it is much smaller than σ T k 2 /(4π) and σ VS k 2 /(4π).

The NMSSM and SIDM
In this section, we will realize the SIDM scenario in the NMSSM, considering the phenomenological constraints. In particular, the relic density and DM direct detection constraints will be discussed in detail.

The NMSSM -masses and couplings of the singlet sector
The NMSSM is a well-motivated extension of the MSSM by a gauge singlet chiral superfield S. Its most general superpotential is: where W Yukawa describes the Yukawa couplings of quark and lepton superfields. We choose µ = 0 in this work, following the convention in Ref. [19] and NMSSMtools [50,51]. The corresponding soft SUSY breaking terms are After the electroweak symmetry breaking, the scalar fields H u , H d and S obtain vacuum expectation values v u , v d and s, respectively. The elements of the CP-even scalar mass matrix square M 2 S in the basis (H d , H u , S) can be written as follows (only those relevant to the singlet are shown): with µ eff = λs and B eff = A λ + κs. We denote the mass eigenstates of the mixing of the CP-even scalars from H u , H d , and S as And the mass matrix of neutralino in the basis (H d ,H u ,S) is read as The mass eigenstates of the neutralinos are denoted as In the limit where the mixing between the singlet scalar (S = s+S h +iSa √ 2 ) and doublet Higgs fields (H u , H d ), as well as the mixing between singlino (S) and Higgsinos (H u ,H d ) are small, we have φ = S h and χ =S with masses given by And the couplings in the singlet sector can be evaluated as [19]:

SIDM in the NMSSM
There are three parameters that control the DM self-scattering, m χ , m φ and the coupling V φχχ . In the NMSSM with nearly decoupled singlet sector, we have Since the singlino-like DM is a Majorana fermion, the spin averaged viscosity cross section σ V is used to describe the non-relativistic scattering between self-interacting DMs and account for the small scale structure of the universe as have been discussed in Sec. 2. The pseudoscalar mediator can not induce large DM selfinteraction [38], as there will be further contributions to the Yukawa potential scaling as m 2 φ /m 2 χ × e −m φ r /r n with n ≥ 2. /m χ (blue). The green lines correspond to the parameters that give the correct relic density. The red, purple and orange contours indicate the splitting probabilities of χ → φχ for initial χ energies (Q max ) 1 TeV, 100 GeV and 22 GeV, respectively. More details for the splitting will be discussed in Sec. 4. Three different masses of the scalar mediator are considered.
According to the simulations in Ref. [13,43], σ V /m χ ∼ 1 − 10 cm 2 /g on dwarf scales (the characteristic velocity is 10 km/s) is needed to solve the core-vs-cusp and too-bigto-fail problems, while the constraints on cluster scales (the characteristic velocity is 1000 km/s) require σ V /m χ ∼ 0.1 − 1 cm 2 /g. In Fig. 2, we show the regions on the κ-m χ plane, where the self-scattering cross section of SIDM is consistent with the simulations on dwarf and cluster scales. The mass of the scalar mediator m φ is chosen to be 0.01 GeV, 0.1 GeV, and 1 GeV, respectively. In the upper left panel where m φ = 0.01 GeV, the blue points which can address the small scale structure problems are in the classical regime on cluster scales (v ∼ 1000 km/s) as they satisfy t > 1 when m χ > 6 GeV. And they are in the quantum regime on dwarf scales (v ∼ 10 km/s) because t < 1 is fulfilled when m χ < 600 GeV. Further calculations show that the blue region shrinks toward smaller m χ when we increase m φ due to the constraints on cluster scale and disappears when m φ is large than 0.07 GeV. As heavy DM scenarios are stringently constrained by DM direct detection experiments, we will focus on DM mass 20 GeV in parameter scan. In the upper right panel where m φ = 0.01 GeV, the blue points are in the quantum regime on both dwarf and cluster scales. In the lower left panel where m φ = 0.1 GeV, the blue points come from the quantum regime (on dwarf and cluster scales). The multi-band structure of the blue region corresponds to the quantum regime in the a − b plane in Fig. 1, which contains many peaks and valleys. It is shown that there is an upper limit on m χ , which can be understood by that t = vm χ /(2m φ ) < 1 (m χ < 2m φ /v) is satisfied in the quantum regime. In the lower right panel where m φ = 1 GeV, the blue points also come from the quantum regime (on dwarf and cluster scales) and thus also have an upper limit on m χ . The blue belts in the m φ = 1 GeV case are narrower than those in the m φ = 0.1 GeV case because larger value of σ V is required to explain the small scale structure for heavier m φ , corresponding to narrower region in the a − b plane, as shown in Fig. 1. Due to the similar reason, the coupling κ needs to be at least ∼ 0.7 in this case.

DM relic density and direct detection constraints
In the NMSSM with almost decoupled singlet sector, the correct DM relic density can be achieved through the χχ → φφ annihilations in s-channel, t-channel, and u-channel. For t-channel and u-channel annihilation, only the coupling V φχχ is relevant. As for s-channel annihilation, V φφφ comes into play. Note that with negligible V φφφ , the density of SIDM tends to be under-abundant, because a large κ and a light φ (m φ < m χ ) is required for solving the small structure problems. Using the triple scalar interaction in NMSSM we can get the total DM annihilation cross section as The relic density of DM can be estimated as [44] where g * ≈ 10 and x f ≈ 10. It is noted that d gets its minimum value (the relic density is maximal) as Because the SIDM is under-abundant in most cases, this limit is useful to test whether the NMSSM with freeze-out mechanism can explain the current DM relic density. If the Ω χ h 2 with A κ ∼ (11/3)m χ is still less than ∼ 0.1, either the model needs to be extended or other DM production mechanisms are required. In Fig. 2, the green lines correspond to the relations of κ and m χ such that the DM relic density calculated by Eq. 3.12 (taking A κ = (11/3)m χ ) is equal to 0.12 [45]. When m φ is greater than 0.1 GeV, the green lines will not cross the blue regions and will further departure away as m φ increases, as the small scale structure problem requires much larger κ than the correct relic density does. However, the green line crosses the blue regions when m φ is 0.01 GeV and smaller, which means in this case, the small scale structure and the DM relic density can be addressed simultaneously. By fixing the value of m χ to make Ωh 2 = 0.12 according to Eq. 3.12, Fig. 3 shows the region where the small scale structure problem can be addressed on the m φ − κ plane. There is an upper limit for m φ (∼ 0.09 GeV). The selected points above (below) the dashed line (t = 1, v = 1000 km/s) belong to the classical (quantum) regime on the cluster scales (v ∼ 1000 km/s), and all the selected points are in the quantum regime on the dwarf scales. Moreover, many DM underground direct detection experiments have put stringent limit on SUSY DM models. For the singlino-like DM in the NMSSM, the nucleon-DM scattering cross section is proportional to the mixings between the singlet scalar and the Higgs boson (the one observed at the LHC, which is H u -like): And the magnitude of the spin independent proton-DM scattering cross section can be estimated by the following equation [46]: GeV. Since the SIDM scenario requests a relatively large κ and light scalar mediator φ, in order to suppress the nucleon-DM scattering cross section below ∼ 10 −4 pb, an extremely small λ is required, i.e. λ 10 −7 , unless the A λ is tuned appropriately to implement an exact cancellation. As the SIDM mass is typically O(10 −1 − 10 1 ) GeV , the most stringent limits come from XENON1T [47], CRESST [48] and DarkSide50 [49] experiments.

Parameter scan
The parameters in the NMSSM most relevant to the scalar sector and neutralino sector are We perform random scan of these parameters in the ranges: Note that the λ is scanned logarithmically (i.e. λ = 10 r and r ∈ [−8, −5] is a uniform random number), in order to focus on the small λ region. To implement the SIDM scenario in the NMSSM, some of the parameters are scanned in much smaller regions:
The scanning is performed with the package NMSSMtools [50,51], imposing the following phenomenological constraints as preselections: • All LEP and Tevatron constraints that are implemented in the NMSSMtools.
• The lightest CP-even scalar h 1 and the lightest neutralino being singlet-like, i.e. with singlet component of the mixing matrix greater than 0.9.
To implement a realistic SIDM scenario in the NMSSM, the parameter λ as small as ∼ 10 −7 is required, such that the DM sector is isolated from the SM sector and the DMnucleon cross section is suppressed. However, with λ 10 −7 and µ eff 100 GeV (required by the chargino search at the LEP), the VEV of the singlet field s = µ eff /λ ∼ 10 9 GeV. Moreover, the SIDM scenario favors light DM and large κ. According to Eq. 3.19 and Eq. 3.20, µ ∼ 10 8 GeV and ξ S ∼ 10 25 GeV 3 (given κ ∼ O(0.1)). The light CP-even scalar (small M 2 S,33 ) requires almost exact cancellation between numbers of order 10 16 . Although the parameter ξ S is set to implement this relation, it is numerically unstable, because we are using double precision floats, the machine precision of which are around 10 −16 . So in practice, we scan the  The results of the parameter scan are demonstrated in Fig. 4, where the grey boxes correspond to the points that pass the preselections. Moreover, The small scale structure constraint (σ v=10 km/s V /m χ ∈ [1, 10] cm 2 /g and σ v=1000 km/s V /m χ ∈ [0.1, 1] cm 2 /g), DM relic density constraint (Ωh 2 ∈ [0.09, 0.11]) and DM direct detection constraints (XENON1T, CRESST and DarkSide50) are applied in order. Points that pass those constraints are shown with different colors. Due to the scanning range that we choose for λ, the typical spin independent cross section between DM and proton σ SI (pχ) is ∼ 10 −7 -10 −3 pb. So the DM direct detection constraints are stringent. An even smaller λ is still possible since it is not relevant for implementing the SIDM scenario. However, we will show later that smaller λ renders long-lived singlet scalar, which will not provide any detectable signatures at colliders, except missing transverse momentum. The relic density constraint excludes all of the points with dark matter mass m χ 1 GeV, because Ωh 2 ∝ m 4 χ according to Eq. 3.12. Taking into account the direct detection constraints, the survival points have DM mass ∼ 1-5 GeV. It is interesting to observe in the middle panel of Fig. 4 that the µ eff (i.e. the Higgsino mass) is bounded from above, especially when the DM mass is heavier than 2 GeV. In order to address the small scale structure problem, the κ/m φ needs to be large. However, relic density and DM direct detection constraints favor small value of κ. They put an upper limit on m φ , i.e. m φ 0.02 GeV.  The spin averaged viscosity cross sections for the points presented in Fig. 4 are obtained by numerically solving the Schrödinger equation with Yukawa potential [13]. In Sec. 2, we have also presented analytical expression for the same variables. For comparison, in the left panel of Fig. 5, we show the distributions of the relative difference between the analytically and numerically calculated σ V /m χ . There are more points with σ V /m χ < σ Comb V /m χ than the opposite case. The relative deviation between two methods is ∼ O(10)%. As discussed in Sec. 2, the σ V can vary across several orders of magnitudes for different a and b, such amount of deviation is acceptable for the purpose of estimation. The velocity dependence of the DM self-interacting cross section is illustrated in the right panel of Fig. 5. It can be found that the direct detection constraints are in contradiction with the velocity dependent feature. Parameter points that pass all DM direct detection constraints can only have the ratio between the σ V /m χ on dwarf and cluster scales at most around 3.

Splitting functions for singlet scalar and singlino
The SIDM scenario requires light scalar mediator and relatively large DM self-coupling. Such parameter setup may lead to DM/mediator showering: the particle can split multiple times during its propagation.
Considering a time-like branching process A → B + C with the off-shell particle A being in the final state of a preceding hard process, the opening angle between A and B (C) is denoted by θ b (θ c ). In the high energy limit, the case that the daughters B and C are nearly collinear to the parent particle A dominates the branching rather than the non-collinear case, as the scattering amplitude is proportional to We parameterize the collinear time-like branching as follows. θ = θ b + θ c is the open angle and satisfies θ 1. The particle mass, energy, and momentum are m i , E i , and Because P a ≈ P b + P c in the collinear branching, we use z = Pc Pa andz = 1 − z representing the momentum fractions of A taken up by C and B, respectively. It is not hard to verify that Q 2 ≈ P b P c θ 2 = P 2 a zzθ 2 in the collinear branching. For the hard process with a final state particle A being the parent particle of a collinear time-like branching, the differential cross section can be expressed as where dP A→B+C is the differential splitting function for the A → B + C. Using the parameters defined above, the splitting function can be expressed as where N = 2 when B and C are identical particles and N = 1 when B and C are different. The |M split | 2 is the spin-averaged matrix-element square for the A → B + C branching process, which can be computed from the amputated A → B + C Feynman diagram with on-shell polarization vectors.
Process In the singlet sector of the NMSSM, branching processes of the singlino/singlet scalar include χ → φ + χ, φ → χ + χ, and φ → φ + φ. The |M split | 2 of these processes which are depend on the helicities of the fermions are summarized in Tab. 1, giving The evolution of the final-state radiation (FSR) is dominated by the splitting functions. For the possible time-like branching of a parent particle A, the famous Sudakov form factor describes A's probability of evolving from Q max to Q 0 without branching, where the allowed z range (z min (Q), z min (Q)) at Q depends on kinematics and is given by We study the evolution of the FSR by a numerical Monte Carlo method with Markov chain based on the Sudakov factors of χ and φ. The evolution is operated by running from a high virtuality scale Q max , chosen to be the CM-frame energy of the hard partonic process, down to a low scale Q min with small Q steps. If the branching of A → B + C takes place at some Q, the evolution will be carried on with both the daughters B and C. In Fig. 2, we plot the contours of the probability of the χ → φ + χ branching taking place, which is computed as P χ→φχ = 1 − ∆ χ (Q max ; m χ + m φ ), with m φ =0.01, 0.1, 1.0 GeV. The branching probability increases with κ and can be above 0.7 when κ ∼ 2.5 and m φ = 0.01 GeV. Lower initial energy Q max leads to smaller branching probability, especially when Q max approaches m χ + m φ , which is the threshold for the branching. When m φ increases, the branching probability drops quickly because the factor of 1/ Q 2 − m 2 χ 2 in the splitting function dP χ→φ+χ (z,Q) dz d ln Q 2 decreases dramatically with m φ when Q approaches Q min = m χ +m φ .

Production of the singlet scalar at the LHC
The couplings between the MSSM sector and singlet sector are suppressed by the tiny λ. However, the particles in the singlet sector can still be produced through the decay of SUSY particles at the LHC. In particular, as shown in Fig. 4, the effective µ parameter is find to be 1 TeV for m χ 10 GeV (note that we scan the µ eff in the range [300,2000] GeV). Thus, relatively light Higgsinos are predicted. Assuming R-parity conservation, the light Higgsinos can be pair produced (in terms ofH ±H ∓ ,H ± χ 2,3 , χ 2,3 χ 2,3 ) at the LHC through s-channel gauge boson exchange with the SM gauge couplings.
At the tree level, the mass splitting between two neutral Higgsinos is ∼ m 2 Z /M 1,2 , and the splitting between charged Higgsino and the lighter neutral Higgsino is approximately half of that. Moreover, radiative corrections further induce ∼ O(100) MeV mass difference between charged and neutral states [52]. As a result, the heavier neutral Higgsino and the charged Higgsino will dominantly decay through χ 3 → Z * (→ f f )χ 2 andH ± → W * (→ f f )χ 2 , which typically happen in time scale of O(10 −15 ) second (prompt decay inside the detector).
The lightest Higgsino goes through two body decays: χ 2 → Zχ, χ 2 → H 2 χ and χ 2 → φχ with branching ratio proportional to λ 2 , λ 2 and κ 2 λ 2 , respectively. Here the H 2 is the SM-like Higgs boson. With λ ∼ O(10 −7 ), the typical lifetime of the lightest Higgsino is O(10 −10 ) second, so that it decays inside the detector. A relatively large κ that is used to address the small scale structure problem, also leads to relatively large branching ratio of χ 2 → φχ, which provide a unique opportunity to probe the singlet sector of the SIDM scenario in the NMSSM at colliders.
In Fig. 6, we present the decay branching ratios of the neutral Higgsinos for the scanned points. As before, the small scale structure constraints, the DM relic density constraint and the DM direct detection constraints are applied in order. For κ 0.1, the lightest Higgsino will dominantly decay into Z/H 2 +χ. As the ATLAS collaboration has conducted searches for Higgsinos decaying into Z boson and Higgs boson, some of the points can be excluded already. In upper panels, the constraints from ATLAS searches for four or more charged leptons [53], at least three b-tagged jets [54] and two photon + missing transverse momentum [55] are projected on. We can find the ATLAS search for four or more charged leptons excludes the points with Higgsino mass 400 GeV, while the constraints on χ 2 → H 2 χ channel are much milder. Points that evade the LHC constraints and all others are shown by green dots in lower plots. The Higgsino masses for allowed points are ∼ 400 − 1000 GeV and the corresponding production cross sections [56,57] are ∼ 88.7 − 0.969 fb (summed cross sections ofH ±H ∓ ,H ± χ 2,3 , χ 2,3 χ 2,3 productions at the 13 TeV LHC). The branching ratios of χ 2 → φχ for allowed points are ∼ O(10 −3 ). So, at the high-luminosity LHC, the signal rate of the Higgsino production with subsequent decay χ 2 → φχ in the SIDM scenario of NMSSM is sizable. In the lower right panel, we show the the summed branching ratios of χ 3 → φχ and χ 3 → φχ 2 for the heavier neutral Higgsino (χ 3 ). As have been discussed above, the χ 3 is dominated by the three body decay with off-shell gauge boson. Its branching ratio to the singlet-like scalar is highly suppressed by the λ coupling.

Benchmark points
For illustration purpose, we provide two benchmark points for the SIDM scenario of NMSSM in Tab 2. With light mediator and relatively large κ, the spin averaged viscosity cross sections for both points are 1.0 cm 2 /g on dwarf scales, thus address the small scale structure problem. The production rates of the Higgsinos pair at the 13 TeV LHC with one of the Higgsino decaying through χ 2 → φχ for benchmark point (i) and (ii)    The most remarkable collider signature of those two benchmark points is the showers of the singlino and the singlet scalar, which are produced by the Higgsino decays, leading to multiple scalars/DMs in the final state. The κ is bounded from above due to the DM relic density constraints. So the DM splitting, as well as the scalar splitting into DM pairs is suppressed. According to our simulation, this kind of splittings happens less than twice per thousand events for our benchmark points. On the other hand, the probability of φ → φφ splitting can be sizable, since the coupling V φφφ can be large (as given in Eq. 3.9).
In the left panel of Fig. 7, we present the multiplicity of φ and χ from the process with subsequent showers for our two benchmark points. Note the lightest Higgsino χ 2 can be produced either from its direct production or from the decay of heavier Higgsinos. The final state φ particles have two origins: from the χ (denoted by χ → φ) and the φ (denoted by φ → φ) splittings in the process 5.1. For both benchmark points, the final state φ particles are mostly from the φ splitting. As the coupling |V φφφ | are 1.40 GeV and 45.9 GeV for benchmark point (i) and (ii), respectively. The benchmark point (ii) has much higher φ multiplicity than the benchmark point (i) due to its larger |V φφφ |. Although it is difficult to obtain a precise analytic expression for the φ multiplicity distribution, the rough estimation for the order of the average multiplicity of φ can be obtained bȳ wheren φ;φ→φ stands for the multiplicity of φ coming from the branching of φ. When taking Q max = 500 GeV, the Eq. 5.2 givesn φ ∼ 28 and 2642 for the benchmark point (i) and (ii), respectively. In the right panel of Fig. 7, the spectra of the final state χ and φ for two benchmark points are shown. The p T spectra of χ are harder than that of φ, since the φ goes through multiple splittings while χ does not split. Due to the same reason, the p T (φ) spectrum of benchmark point (ii) is much softer than that of benchmark point (i).

Beyond the NMSSM
The κ coupling can be as large as O(1) when addressing the small scale structure problem. However, in the NMSSM, in order to implement the correct relic density with freeze-out mechanism, we find κ 0.2 is required (as shown in the right panel of Fig. 4). Moreover, applying all phenomenological constraints leads to sizable V φφφ coupling. Thus the φ splitting is copious while the χ splitting is highly suppressed in the NMSSM.
In this subsection, we try to study the splitting of singlet/singlino in a simplified model framework where the small scale structure problem can be addressed while leaving out the DM relic density and direct detection constraints (this may be realized in a UV-complete model, where the DM is not produced from the thermal freeze-out mechanism). As a result, the κ can be large and the coupling of V φφφ is a free parameter. For illustration, we consider two benchmark points: • Benchmark (iii): m φ = 0.1 GeV, m χ = 5.65 GeV, κ = 2.5, A κ = 11/3m χ . Since the DM self-interacting viscosity cross section is irrelevant to the A κ parameter, both points have σ v=10 km/s V /m χ = 6.16 cm 2 /g and σ v=1000 km/s V /m χ ∼ 0.938 cm 2 /g, addressing the small scale structure problem. Given κ ∼ 2.5, the branching ratio of the χ 2 → φχ channel is dominant, since its partial width is proportional to κ 2 . So the production rate for the φ and χ is much enhanced. Moreover, with sizable κ, the χ → φχ and φ → χχ splittings are no longer negligible. In Fig. 8, we present the multiplicities and p T spectra for the singlet scalar φ and the singlino χ. Similar as before, we are considering the process in Eq. 5.1 for their production, and the mass of the lightest Higgsino is set to 350 GeV. Because the |V φφφ | coupling is 133 GeV for benchmark point (iii) and vanishes for benchmark point (iv), we can expect that the φ multiplicity from φ → φ of benchmark point (iii) is high. There is also a certain amount of φ coming from the χ splitting, due to the large κ. So the φ multiplicity of benchmark point (iv) can also reach ∼ 5 in spite of vanishing V φφφ . Two benchmark points have similar distributions of χ multiplicity, and the χ particles is less abundant comparing to the φ particles. As for the p T spectrum, the φ of the benchmark point (iii) is much softer than φ of the benchmark point (iv) and χ in both points, mainly because of the high multiplicity of φ such that the energy of the original φ is shared among its daughters.

Decay of the singlet scalar
The light scalar meditator in SIDM scenarios usually has mass less than twice of the pion mass, i.e. m φ 270 MeV. In this case, it mainly decays into muons and electrons. The leading order decay width is given by where θ hφ is the mixing angle in the scalar sector (given by Eq. 3.15 in NMSSM), and With θ hφ ∝ λ ∼ 10 −7 , the typical lifetime of the singlet scalar (with mass ∼ 10 MeV) can reach O(10 4 ) seconds. As a result, the Higgsino decay and the subsequent showering inside the detector can only produce invisible final states. Moreover, the existence of a light long-lived scalar may spoil the success of the Big Bang Nucleosynthesis (BBN). In fact, the tension between the BBN (which requires lifetime of the mediator τ φ 1 second) and direct detection constraints (which favor small mediator-SM Higgs mixing) is commonly exist in SIDM scenarios with scalar mediator. There have been many solutions proposed to alleviate the tension. One class of solutions assumes different DM thermal history [15,25,58]. For example, in Ref. [15], they assume the scalar mediator to be stable and annihilate into an extra scalar, constituting a subdominant DM. In the second class of solutions, the DM direct detection rate is suppressed by introducing CP-violating scalar sector [38] or inelastic DM-nucleon scattering [37], so that a large singlet scalar-SM Higgs mixing is allowed, i.e. short lifetime of the scalar. And the simplest solution is to introduce new particles and interactions [22,36,39] such that the scalar mediator can decay almost promptly through new channels. Different solutions will lead to dramatically different phenomena. For the solutions in later two classes, the scalar mediator can decay into visible final states inside the detector. Each χ or φ from the Higgsino decay will lead to a jet-like object with the identities of its constituents depending on the decay modes of the scalar mediator. Taking the φ → e + e − as an example (which is the dominant decay channel for our benchmark points), the corresponding signatures of χ and φ from Higgsino decay for each benchmark point are listed in Tab. 3. Due to the small κ of benchmark point (i) and (ii), the χ does not split and simply behaves as missing transverse energy ( / E T ) at the detector. After the shower, the φ for benchmark point (i) will produce around O(10) φ particles in the final state (through φ → φφ splitting). Both φ and χ for benchmark point (iv) will produce a few φ particles in the final state (through χ → φχ splitting). This leads to lepton-jet   [27,[59][60][61] for each of the initial φ and χ. The ATLAS collaboration has searched for both prompt and displaced lepton jets [62,63], aiming to the mediator mass 0.1 GeV. Moreover, the mediator multiplicity in each lepton jet is not larger than four. So those existing searches are not optimal for our benchmark points. For benchmark point (ii) and (iii), the φ → φφ splitting is so copious such that the transverse momenta of final state leptons are below the selection threshold in lepton jet searches (which requires energy of electron to be greater than 10 GeV). As a result, the signatures of φ and χ become Soft Unclustered Energy Patterns (SUEP). There have been some studies on searching for SUEP [64][65][66], focusing on hadronic final states. The detailed study for the collider search of our benchmark points will be addressed in future works.

Conclusion
We study the SIDM scenario in the general NMSSM and beyond, where the DM is a Majorana fermion and the force mediator is a scalar boson. Due to the relatively large couplings and light scalar mediator in this scenario, the DM/mediator will go through multiple branchings if they are produced with high energy, leading to the signature of multiple scalars in the final state.
The DM self-interaction cross section is calculated in both analytical and numerical ways. In particular, an improved analytical estimation for the DM self-interacting cross section which takes into account the Born level effects is proposed. Based on the scanned points in the NMSSM, we demonstrate that the analytical expression matches well with the numerical solution. In most case, the relative difference of DM self-interacting cross section between two methods is ∼ O(10)%.
Two benchmark points in the general NMSSM that address the small scale structure problem with correct DM relic density and satisfy other phenomenological constraints are explicitly given for the first time (one of the benchmark points is challenged by the DM direct searches). They are featured by relatively small µ parameter, light singlet-like scalar mediator and relatively large triple scalar interaction. In order to have enough DM relic density, κ 0.2 is not allowed in the NMSSM if the DM is lighter than 20 GeV. Otherwise, the annihilation of χχ → φφ will dilute DM density efficiently. The branching ratio of the lighter neutral Higgsino decaying into the singlet scalar and singlino is proportional to κ 2 , which is around 10 −3 . The production rates of the Higgsinos pair at the 13 TeV LHC with at least one of the Higgsino decaying into the singlet scalar and singlino are 0.108 fb and 8.95 fb for benchmark point (i) and (ii), respectively.
A Monte Carlo simulation of the DM/mediator showers are implemented by building the multiple branchings as a Markov process based on the Sudakov form factors. For the two benchmark points in the NMSSM, κ ∼ 0.1, which means the splitting of χ → φχ is suppressed. On the other hand, the splitting probability of φ → φφ is high. The number of φ (with highest probability) from the Higgsino production and decay at the LHC can reach ∼ 20 for benchmark point (i) and is even more for benchmark point (ii). Two benchmark points beyond the NMSSM with κ = 2.5 and address the small scale structure issue are proposed as well, to illustrate the cases with sizable χ → φχ splitting. When φ → φφ is turned off (corresponding to benchmark point (iv)), number of φ (with highest probability) in the final state is around 3 ∼ 4. The p T spectra of DM/mediator are also shown for benchmark points.
Finally, we comment on the tension between the Big Bang Nucleosynthesis (BBN) limits and dark matter direct detection constraints for the SIDM scenario with light longlived mediator. Extensions to the simple singlet-dominant DM sector are required to alleviate the tension. Depending on the form of possible extension, the energetic φ and χ from Higgsino decay can induce remarkable signatures at the LHC.