Scope of strongly self-interacting thermal WIMPs in a minimal $U(1)_D$ extension and its future prospects

In this work we have considered a minimal extension of Standard Model by a local $U(1)$ gauge group in order to accommodate a stable (fermionic) Dark Matter (DM) candidate. We have focussed on parameter regions where DM possesses adequate self interaction, owing to the presence of a light scalar mediator (the dark Higgs), alleviating some of the tensions in the small-scale structures. We have studied the scenario in the light of a variety of data, mostly from dark matter direct searches, collider searches and flavour physics experiments, with an attempt to constrain the interactions of the standard model (SM) particles with the ones in the Dark Sector (DS). Assuming a small gauge kinetic mixing parameter, we find that for rather heavy DM %$\gtrsim \mathcal{O}(1-10)\,\, {\rm GeV}$%, the most stringent bound on the mixing angle of the Dark Higgs with the SM Higgs boson comes from dark matter direct detection experiments, while for lighter DM, LHC constraints become more relevant. Note that, due to the presence of very light mediators the usual realisation of direct detection constraints in terms of momentum independent cross sections had to be reevaluated for our scenario. In addition, we find that the smallness of the relevant portal couplings, as dictated by data, critically suppress the viability of DM production by the standard"freeze-out"mechanism in such simplified scenarios. In particular, the viable DM masses are $\lesssim \mathcal{O}(2)$ GeV $i.e.$ in the regions where direct detection limits tend to become weak. For heavier DM with large self-interactions, we hence conclude that non-thermal production mechanisms are favoured. Lastly, future collider reach of such a simplified scenario has also been studied in detail.


Introduction
The evidence for existence of non-relativistic non-luminous Dark Matter (DM) has been overwhelming. Several astrophysical and cosmological observations have indicated the presence of DM [1]. Within the paradigm of the standard model of cosmology (ΛCDM) recent measurements of the Cosmic Microwave Background Radiation (CMBR) estimated that about 26% of the energy budget of our Universe consists of DM [2,3] 1 . CMBR, together with Big Bang Nucleosynthesis (BBN) requires DM to be non-baryonic. Astro-physical objects, especially primordial black holes (PBHs) have been considered as DM candidates. While such PBHs, for certain windows of mass 2 can account for the entire energy budget of DM, it has been shown that adequate production of PBHs after (single-field slow-roll) inflation can be difficult [5,6]. However, the standard model (SM) of particle physics, does not incorporate any suitable DM candidate. Various extensions of the SM has been considered in the literature [7,8]. A generic aspect in DM models concerns about the stability of DM. Constraints from structure formation, indirect searches [9][10][11] and CMB [12] require DM to be very long-lived (life-time τ 10 26 s depending on the decay modes). The simplest models introduce a global symmetry to prevent the decay of DM. However, it has been argued that such global symmetries may be broken due to gravitational effects [13][14][15]. This, in turn, can induce Planck-scale suppressed effective terms contributing to the decay of the DM [16]. In this article, therefore, we extend the SM with a simple continuous local symmetry U (1). We further assume that the DM is charged under this symmetry, and therefore, is stable. Earlier works that also explored simplest extensions of Standard Model with similar phenomenological signatures can be found in [17][18][19][20][21][22][23][24] Recent results from the Large Hadron Collider (LHC), so far, present no convincing evidence for new physics, see e.g. [25][26][27] for implications on DM models. Further, no evidence for particle DM has emerged from direct [28][29][30] (and indirect [31,32]) searches of DM. The Standard Model (SM) of particle physics remains a good low energy effective description. Stringent constraints on most well-motivated new physics scenarios present two distinct possibilities: any new physics may be beyond the energy reach of the LHC, and perhaps can only be probed indirectly through their contribution to the higher dimensional effective operators; or new physics, if exists at (or below) the electroweak scale, may be very weakly coupled to the SM. In the subsequent discussion, we will assume the latter and consider the Dark Sector (DS) particles to be accessible at LHC.
There is an additional motivation for this consideration. Note that in spite of the enormous success of ΛCDM in the cosmological scales, there have been concerns, in particular when it comes to the small scale structures (see e.g. [33][34][35] for recent reviews). The most notable ones include Core-vs-cusp [36,37] and Too-big-to-fail [38][39][40][41] 3 . It has been argued that self-interacting DM can also simultaneously resolve these issues [46,47], see also [48,49] for a review. Although this requires a very large selfinteraction cross-section, it remains consistent with the present constraints on DM self-interaction [50][51][52][53] 4 5 6 . Further, it has been argued that a velocity dependent self-interaction cross-section is favoured [60]. In the presence of light mediators, it is possible to generate large (velocity dependent) self-interaction among DM particles, thanks to the Sommerfeld enhanced self-scattering cross-section [49,[61][62][63][64]. In the present work, we will ensure that the U (1) extended DS consists of at least one light mediator which facilitate adequate self-interaction among the DM particles.Imposing these criteria, we will investigate the implications on this simple DS. Further, combining constraints from collider searches, flavour physics, beam dump experiments and direct detection of DM, we will try to put restrictions on its interaction with the SM sector.
This article is organised as follows : In Sec. 2 we describe the minimal U (1) D model and its particle content in detail. In Sec. 3 the parameter space relevant for the study of large self interaction of dark matter is motivated quantitatively. Next, concentrating on the light H 1 region as motivated from self-interactions, impact of several experiments (like beam dump, flavour physics and collider) on our parameter space of interest is described studied in Sec. 4 and Sec. 5. Sec. 6 deals with bounds from dark matter phenomenology and discusses its implications. In sec.7 we have presented a detailed study about the future prospects of our scenario from a collider perspective. Finally we summarize and conclude in Sec. 8. 3 Missing satellite problem [42,43] has also been extensively discussed. However, the discovery of faint dwarfs has lead to dissolution of this issue [44,45]. 4 See also [48] and references there for a compilation of constraints from bullet-cluster [54][55][56], ellipticity of the halos and substructure mergers. 5 An observed separation between the stars of a galaxy and the DM halo, while the galaxy falls into the core of a galaxy cluster Abel 3827, has been recently observed [57]. When interpreted as due to DM self-interaction, this leads to σ m 1.5(3) cm 2 gm −1 for contact (long-range) interactions [58].
However, the data has been argued to be also consistent with standard collisionless DM [59]. 6 While baryonic feedback can possibly address some of these issues, it has been argued to face difficulties in the context of field galaxies. For a review see [48].

Description of the Model
As described in the introduction, we consider a simple extension of SM where the stability of DM has been attributed to a local symmetry. However, any unbroken local symmetry would imply the existence of a massless gauge boson. Cosmological observations together with the success of BBN constrain the presence of such a boson through the tight limits of extra relativistic dof. This demands the U (1) D to be broken, and we employ Higgs mechanism to achieve this 7 .
We have assumed minimal particle content for a Dark Sector with a local U (1) D gauge symmetry. Two Weyl fermions needed to be introduced so that the gauge anomaly associated with U (1) D group cancels, and the theory remains anomaly free. We will further elaborate on the motivation behind this particular choice. The Dark Sector is assumed to have the following particle content: In Table 1, ξ 1 and ξ 2 repre-  where "DS" and "portal" denote Dark Sector and mediator respectively. The two component Weyl spinors ξ 1 and ξ 2 can be expressed as a four component fermion χ = (ξ 1 ,ξ 2 ) T , whereξ j = −iσ 2 ξ * j . With this notation, we have, where γ µ in the Weyl representation has been assumed. Further, in the lagrangian above D µ = ∂ µ − ig D q D Z µ D and χ c = iγ 2 γ 0χ T = (ξ 2 ,ξ 1 ) T . Note that a U (1) D charge of 1(−1) for χ (χ) together with a charge of ∓2 for ϕ ensures the invariance of the Yukawa term under the guage group U (1) D .
The Dark Sector interacts with the SM sector via the following gauge invariant terms.

3)
Although none of the particles in the low energy spectrum are charged under both U (1) Y and U (1) D gauge symmetry, we still keep the gauge invariant phenomenologically viable kinetic mixing term [69] in the lagrangian, which can possibly be generated due to high-scale physics. We parametrize the mixing term with parameter g . The scalar field ϕ can couple with the SM Higgs via the usual portal interaction term. We will elaborate on the constraints on both mixing terms in a subsequent discussion.

The Higgs Sector
The scalar potential is given by, The stability conditions can be read off as follows : Since G SM × U (1) D is broken spontaneously, we require µ 2 H < 0 and µ 2 ϕ < 0. Let v h and v ϕ denote the vevs of the scalar fields responsible for the spontaneous breaking of G SM × U (1) D . 8 Around the minimum of the scalar potential, these scalar fields can be described as, Here, φ R (h) and φ I (σ) are the CP-even (CP-odd) parts of φ and H respectively and σ + is a complex scalar. In the unitary gauge, where the low energy d.o.f constitutes only the physical fields, the form of the corresponding expressions may simply be obtained by setting ϕ I , σ + and σ to 0 in equations (2.6) and (2.7).
The squared mass matrix for the scalar fields, in the basis {ϕ R , h}, is given by, The orthonormal mixing matrix N s , is then given by, (2.10) The mass eigenvalues can simply be obtained as the elements of the diagonal matrix (M D s ) 2 = N T s M 2 s N s and the corresponding eigenvalues give the physical masses of the CP-even neutral scalar particles, The mass eigenstates are given by, 12) or more explicitly, (2.14) We will denote the mass eigenstate corresponding to the lightest eigenvalue M H 1 as H 1 and that to the heavier one (M H 2 ) as H 2 . In this convention cos θ mix simply denotes the ϕ R content in the lightest mass eigenstate. As we will consider a scenario with a light H 1 , H 2 would be the SM-like Higgs boson. therefore, cos θ mix denotes the SM Higgs content in H 2 as well as ϕ R content in H 1 .

The Gauge Boson Sector
The gauge group in the present context is : G SM × U (1) D , where G SM denotes the SM gauge group. In this section we focus on the electro-weak gauge bosons and implication of the gauge kinetic mixing with the U (1) D gauge boson. The relevant lagrangian for the neutral gauge boson sector is given by, where, The above equations are written in the (non-canonically normalized) basis ,V µ = (Ŵ 3 µB µẐDµ ) T . Further, where we have kept only terms upto O( 2 g ). The invariance of the gauge covariant derivatives imply that the gauge couplings G in the modified basis are given by, , andĜ = Diagonal(g 2 , g 1 , g D ). In this basis, the mass matrix takes the following form, In this context, we have used the approach described in refs [70,71]. We have used SARAH and SPheno [72,73] to numerically diagonalize the mass matrix and thus obtain the respective mass eigenstates and the corresponding mixing matrix. The mass eigenstates include a zero mass state, corresponding to the unbroken electro-magnetic gauge group U (1), and two massive states representing Z and Z D bosons respectively. Needless to mention, the mass of the Z boson, and the electro-weak precision data in general, constraints g [74,75]. For the present work, we will mainly consider small enough g only ensuring prompt decay of Z D at high energy colliders.

The Dark Sector
After EWSB, ϕ R assumes a vev (v ϕ ) generating a Majorana mass term in the dark fermionic sector. The mass term for the fermions, then, can be expressed as : where This mass matrix can be diagonalized by an orthonormal matrix The mass of the physical states can be obtained from the eigenvalues of M DM and are given by M ± = |M ± f v ϕ |. The corresponding physical states (mass eigenstates) are The state with the lower mass i.e. M − is our dark matter candidate and from now on its mass will be denoted by M DM for definiteness. 9 The lagrangian L DS , in the mass basis, is given by : In the expression above H 1 and H 2 denote the ϕ R -like and h-like states (only true for small mixing), respectively, in the mass eigenbasis as discussed before. 9 Note that we have used the absolute value of the smallest eigenvalue, since the sign of the fermion mass can simply be rotated away by a chiral rotation. For example, with M > f v ϕ and M > 0, the smallest eigenvalue M − f v ϕ is positive. The corresponding mass eigenstate is 1 3 Self-interaction of dark matter: Allowed regions A study of self interaction in our dark matter scenario can give rise to novel features and lead to interesting results. As was discussed in the introduction, strongly selfinteracting dark matter can help to alleviate many small scale structure problems. In our model there are two Majorana fermions, an extra gauge boson and an extra Higgs-boson. Hence, it is generic to any such U (1) gauge theories to have vertices like • When M H 1 is much lighter than the dark matter mass, but M Z D is heavier.
Sommerfeld enhancement takes place by exchanging the light Higgs boson only by the usual ladder diagrams.
• When Z D is lighter than the new Higgs, then enhancement proceeds by exchanging only this light extra gauge boson.
• When both Z D and the new Higgs are of comparable masses then enhancement should in principle occur due to both the mediators.
The case when both the new gauge boson and the extra Higgs are very light may be interesting from a theoretical point of view. But, by making both of these particles light (lighter than the DM) simultaneously, we will lose out on the robust collider signatures that we also intend to explore in this work. So we do not probe those parameter regions where both of the mediators are light. Next, let us discuss about the second case. Since the new scalar is not very light, this implies that λ mix v D is not small (from the expression of the mass of the new scalar). But on the other hand, we demand that M Z D ∼ g D v D should be much smaller than that of M H 1 . To satisfy both of these conditions we have to take resort to very small dark gauge boson couplings. But the enhancement factor also depends on g D (through χ + χ − Z D vertices). Hence in this case, the enhancement that we were expecting by introducing light gauge boson mediator is nullified by the presence of small gauge couplings. So, practically speaking, Sommerfeld enhancement with only light Z D is not a good a choice either.
From these two points of view, it hence seems to take the first case as our preferred choice for studying strong self interactions of dark matters.
In calculating the self-interacting cross section we closely followed the analytical expressions presented in [64]. In the Born limit (α D M DM /M H 1 1), the cross section is given by : where, α D = f 2 /4π and v is the virial velocity of galaxies. Outside the Born regime (α D M DM /M H 1 1) non-perturbative effects become important. Analytical results can be obtained in the classical limit (M DM v/M H 1 1). For an attractive potential [64,76,77] : . Analytical results can also be obtained in the resonance region in between by approximating the Yukawa potential by Hulthen potential [78].
σ Hulthén where the l = 0 phase shift is given in terms of the Γ-function by and κ ≈ 1.6 is a dimensionless number. The bound on self-interacting dark matter is given in terms of σ self /M DM . This should be (0.1 − 10) cm 2 gm −1 when v 0 ∼ 10 km/sec to alleviate the Core-vs-cusp problem of the dwarf spheroidal galaxies [79,80]. From the X-ray and the lensing observations of Bullet cluster we further have σ self /M DM 1 cm 2 gm −1 [56] when v 0 = 1000 km/sec.
The self-interaction cross section in units of the dark matter mass (for a fixed α D and M DM ) is plotted with respect to the mediator mass in Fig. 1. Since we are concentrating on the low mediator mass regime in this study, the mass of the light Higgs was varied up to 10 GeV throughout this work. If we now perform a scan by varying both the dark matter mass (0.1 GeV-1 TeV) and the mediator mass (0.1 MeV-10 GeV), and the allowed points (i.e. points with 0.1 cm 2 gm −1 σ self /M DM 10 cm 2 gm −1 ) are plotted, then we arrive at Fig. 2 (left). This was however done for a fixed α D (≡ f 2 /4π) = 0.001. The right hand panel of Fig. 2 shows the more general case with varying α D . The values of α D are shown in the colour bar. We see immediately that for light mediators (with mass reaching up to 10 GeV) and for a wide range of dark matter masses, the allowed range of α D is quite large. Note that, self interaction can be achieved even for very small values of α D (∼ 10 −7 ) when the mediator is extremely light (∼ 100 keV). Here, we have taken care of the bounds on the self-interaction cross-section arising from the Bullet Cluster.
So, to summarise, in this work we will be concentrating on self interaction of dark matter mediated by light scalar mediators alone. But, this mediator (H 1 ) cannot be extremely light in order to respect the constraints arising from BBN [81,82]. This is because, BBN places strong upper bounds on presence of relativistic particles when Figure 2: Left: Region for allowed self-interaction with varying dark matter and mediator mass for a fixed α D . v 0 , the virialized velocity is fixed to 10 km/sec. The region within the black contour denotes the part of the parameter space ruled from Bullet cluster [54][55][56] constraint (see text for more details). Right: Region of parameter space where moderate to strong self-interaction of dark matter is allowed to alleviate Core-vs-cusp, Too-big-to-fail problems, with varying α D (0.1 cm 2 gm −1 σ self /M DM 10 cm 2 gm −1 ). The variation of α D is shown in the color bar. The allowed points shown in this figure satisfies the constraint arising from Bullet Cluster. temperature of the universe is ∼ 1 MeV. Hence, for all practical purposes, M H 1 10 MeV would be admissible. Therefore, the mass range of the scalar H 1 for suitable self-interaction is expected to be : 10 MeV M H 1 10 GeV. We will concentrate on this range of M H 1 throughout the rest of the work. The upper bound on selfinteraction cross section from Bullet Cluster was also taken into consideration during our calculation.

Survey of parameter space
The latest results from the SM Higgs coupling measurements still allow the SM Higgs to have non-standard couplings. In the context of the U (1) D model considered in this work, such non-standard decay modes of the SM-like Higgs boson could arise in three possible ways, viz., decay into a pair of light Higgs bosons (H 2 → H 1 H 1 ), decay into a pair of dark gauge bosons (H 2 → Z D Z D ) and decay into the Majorana fermions (H 2 → χ + χ + , χ − χ − ). The H 2 → H 1 H 1 decay width is crucially controlled by sin θ mix , while the H 2 → Z D Z D decay process also has a direct dependence on 2 . Throughout this analysis, we restrict the mass of Z D within M Z D < 60 GeV. Thus, the H 2 → H 1 H 1 and H 2 → Z D Z D decay modes always remain kinematically feasible.
We perform a random scan over the 9-dimensional parameter space to capture the underlying physics of the U (1) D model. In this respect, we scan over the following parameters : λ H , λ φ , λ mix , g D , g , v D , M , and f .
With the main focus on exploring the self-interaction motivated parameter region, we perform a random scan over the following range of input parameters. The scan range is inspired by the region of parameter space shown in Fig. 2 (b) which corresponds to those sectors which facilitate a dark matter with moderate to strong self-interactions. It can be observed from Fig. 2 that a self-interacting DM up to mass 1 TeV could be accommodated with M H 1 < 10 GeV. The mass of H 1 (Eqn. 2.11) is determined by λ φ , λ mix and v D , and the range of these parameters is chosen in such a way that a significant fraction of scanned points populate the M H 1 10 GeV region. Correspondingly, M has also been varied from 1 GeV − 1000 GeV to obtain the DM mass in that range. Since we are interested in studying the collider prospects of Z D , whose mass is given by g D × v D , the choice of g D is governed primarily by the requirement of M Z D within 2 GeV − 60 GeV. −0.10 < λ H < −0.16, − 10 −8 < λ φ < −0.7, 10 −12 < λ mix < 0.3 (4.1) 10 −6 < g D < 0.6, 10 −4 < g < 10 −10 (4.2) 1 GeV < v D < 1000 GeV, 10 −6 < f < 10 −1 , 1 GeV < M < 500 GeV (4.3)

Constraints
The U (1) D parameter space discussed so far gets constrained by a multitude of experimental search results.  Fig. 3 in the λ mix − M H 1 plane. The color palette represents the value of sin θ mix , and exhibits its direct proportionality with λ mix . The effect of constraints on sin θ mix from light H 1 searches at LEP has been analyzed in this section. We also evaluate the implications from Higgs signal strength measurements through a global χ 2 analysis, by combining, both, 8 and 13 TeV LHC results. The status of our parameter space in light of direct Z D searches at the LHC is studies as well. The U (1) D parameter space under study also receives strong constraints from measurements at B-factories and various beam dump experiments. In the remainder of this section, we will analyse in detail, the implications on our parameter space from each of these constraints.
Upper bounds on H 1 V V (V = W ± , Z) coupling normalized to that for SM, at 95% C.L., taking into account, both, LEPI and LEPII data has been derived in [84]. In  [84], where ζ is the normalized H 1 V V coupling with respect to SM. The green colored points are excluded by the upper limits on κ derived by OPAL [85], where κ is the ratio of production cross-section of H 1 in the Higgs strahlung process to its SM value. the context of our analysis, the normalized H 1 V V coupling (ζ) is directly proportional to sin θ mix , and the upper limits derived in [84] exclude the blue colored points in Fig. 4. It can be observed from Fig. 4 that the upper limits on H 1 V V coupling excludes sin θ mix 0.5, for all values of M H 1 obtained in our parameter space scan. The study by OPAL collboration [85], using the LEPI and LEPII dataset, has derived upper bounds on κ at 95% C.L., for M H 1 = 10 −6 − 100 GeV, where κ is the ratio of production cross-section of the new light scalar in the Higgsstrahlung process to that of SM Higgs production in the Higgsstrahlung process, under the assumption that the mass of the new light scalar is equal to the SM Higgs. Within the U (1) D model considered here, κ is proportional to sin 2 θ mix , and the corresponding upper bounds, upon being implemented on our parameter space, excludes the green colored region in Fig. 4. We find that the upper limits from [85] exerts a relatively stronger constraint on the parameter space and excludes sin θ mix 0.2, over the entire range of M H 1 10 obtained in our parameter 10 It may be noted that decay channels originating from the light scalar Higgs (M H1 2 GeV) space scan.

Constraints from LHC Higgs signal strength measurements
The ATLAS and CMS collaborations have performed numerous measurements of the coupling of the 125 GeV Higgs boson using LHC Run-I and Run-II datasets. Results from these measurements are usually presented through signal strength observables, which considers the most important Higgs production modes at LHC : gluon fusion mode (ggF ), vector boson fusion (V BF ), associated production with vector bosons (V H 2 , V = W ± , Z) and associated production with a top-antitop pair (tth). The relevant Higgs boson decay modes are H 2 → bb, τ + τ − , W W , ZZ, γγ. The signal strength variable is defined as Here, σ i H 2 corresponds to the Higgs production cross-section in the i th mode (i = ggF , V BF , V H 2 and ttH 2 ), and Br j H 2 corresponds to the branching fraction of the Higgs in the j th decay mode (j = bb, τ + τ − , W W , ZZ, γγ). (σ i H 2 ) SM and (Br i H 2 ) SM corresponds to the SM counterparts. In the current analysis, the heavier scalar Higgs boson, H 2 , is required to be consistent with the SM 125 GeV Higgs boson, thereby requiring it to be dominantly doublet-like. However, mixing between the doublet and singlet Higgs fields renders a small singlet admixture in H 2 as well. The coupling of H 2 with the SM particles remain similar to that of the case of SM Higgs boson, except for an additional suppression (by cos θ mix ). Consequently, the ttH 2 vertex in the ggF production mode of Higgs, and the V V H 2 vertex in the V BF , V H 2 and ttH 2 production modes of Higgs, gets an additional factor of cos θ mix , resulting in the Higgs production cross-section acquiring a cos 2 θ mix suppression. As a result, the ratio of σ i H 2 M odel /σ H 2 SM in Eqn. 5.1 becomes equal to cos 2 θ mix . The branching fraction of H 2 to SM final states also gets affected by the presence of new non-SM decay modes. As specified in the previous section, H 2 has three possible non-SM decay modes, and the relative interplay of input parameters determine the partial decay width in each of the SM channel.
In this study, the signal strength constraints have been imposed upon the scanned parameter set through a global χ 2 analysis 11 performed by taking into account the require a careful treatment owing to the uncertainties associated with partial decay widths of H 1 into hadronic channels. Within 2m π M H1 2 GeV, the uncertainties associated with theoretical calculations remain significant [89], and the corresponding partial decay widths are computed using low energy effective theories of QCD. In our analysis, we have considered the branching ratios of H 1 → γγ, e + e − , µ + µ − from [90]. 11 This analysis set up has been validated in [91,92]. most recent Higgs signal strength constraints, tabulated in Table II and Table III of [92]. The value of χ 2 is computed as where, x i corresponds to the best-fit value of the observable derived through experimental measurements,x i corresponds to the value of the observable computed for the current model, and ∆x i refers to the error associated with the experimental measurement. In the context of this study, x i represents the best-fit value of the signal strength observables. The value of χ 2 was computed for all scanned parameter space points by combining 28 signal strength observables from LHC Run-I data and 18 observables from LHC Run-II data (∼ 15 f b −1 and ∼ 36 f b −1 ) 12 , and the lowest value of χ 2 was determined (represented by χ 2 min ). Allowing 2σ uncertainty, we choose parameter space points which lie within χ 2 min + 6.18. The implications of the global χ 2 analysis are shown in Fig. 5. We would like to note that the parameter space points corresponding to Fig. 5 (left) have been generated from the parameter space specified in Eqn.4.3, except for f and M , whose values were fixed in order to generate M DM = 500 GeV. Thus, H 2 → H 1 H 1 and H 2 → Z D Z D are the only non-SM decay modes for H 2 . We show the correlation between these two non-SM decay modes in Fig. 5 (left), where the grey points (color palette points) have been excluded (allowed) by the global χ 2 analysis. The H 2 Z D Z D coupling emerges from the D µ φ † D µ φ term, when one of the singlet Higgs field receives a vacuum expectation value. The covariant derivative contains a term ∝ g D Z D , and the φ field yields a term proportional to sin θ mix H 2 , resulting in the H 2 Z D Z D coupling to become proportional to g 2 D × sin θ mix . Another contribution arises from the SM term D µ H † D µ H through Z − Z D mixing. However, this term is proportional to 2 g and since we have restricted ourselves to small values of g , contributions from this term can be safely ignored. The H 2 H 1 H 1 coupling manifests from the quadratic Higgs mixing term in the scalar potential,    [93]. This search probed the mass range of 1 GeV < M A < 60 GeV. The same dataset was also used by ATLAS to probe the 2b2µ final state produced via H SM → SS → 2b2µ [94] and 4b final state originating via H SM → SS → 4b [95],  [93]. All parameter space points in this figure satisfy the Higgs mass constraints. The blue colored points are the ones which are still allowed by the

Constraints from LHC direct searches
where S is a spin-zero boson, in the mass range 18 GeV < M S < 62 GeV. The CMS collaboration has also searched the 4µ final state produced via H SM → SS → 4µ, and derived upper limits in the mass range 0.25 GeV < M S < 8.5 GeV. This search used the 13 TeV dataset collected at L ∼ 35.9 f b −1 [96]. The same dataset has also been used by CMS to search for exotic decays of the Higgs boson to a pair of light A in the 2b2τ [97] and 2µ2τ [98] final state, focussing on the mass range 15 GeV < M A < 62 GeV. ATLAS and CMS have also performed similar searches using the LHC Run-I dataset for a multitude of final states : 4τ [99], 2µ2b [99], 2µ2τ [99,100], 4µ [101,102], 4τ [103,104], 2τ 2b [105].
The ATLAS collaboration has searched for the dark vector boson (Z D ) in the mass range of 1 GeV M Z D 60 GeV, where Z D can be produced as pair (Z D Z D ) or in association with SM Z boson and eventually the Z D Z D /ZZ D pair decays into 4l(l = e, µ) final state [93]. Upper limits were obtained on the quantity σ/σ SM × Br(H 2 → ZZ D ) or σ/σ SM ×Br(H 2 → Z D Z D ), where σ and σ SM are the production cross-sections of the SM-like Higgs boson in the NP and SM scenarios. It is to be noted that the limits obtained in [93] assume Br(Z D → e + e − ) ∼ 50% and Br(Z D → µ + µ − ) ∼ 50%, resulting in 4l = 4e(25%), 2e2µ(50%), 4µ(50%), and therefore, a correct scaling is required while evaluating the implication of these constraints on the NP model under consideration. The decay processes : H 2 → ZZ D and H 2 → Z D Z D , depend on the kinetic mixing factor ( g ), and are independent of the mixing between the Higgs doublet from the SM and the singlet Higgs from the dark sector. As a result, these decay channels serve an excellent probe of g . However, in the case of H 2 → ZZ D → 4l search channel, the SM H SM → ZZ * → 4l process offers a strong irreducible background, and eventually dilutes the resolution between the signal and background, rendering these search limits insensitive to the low g region which is relevant to our parameter space. On the other hand, the search channel, H 2 → Z D Z D → 4l stands on an advantageous ground owing to the possibility of application of SM Z-boson veto, which significantly improves the signal sensitivity as compared to the earlier case. We show the implications from the current Z D limits from H 2 → Z D Z D searches (from [93]) on our parameter space in Fig. 6. The vertical axis corresponds to the σ/σ SM × Br(H 2 → Z D Z D ) and the horizontal axis corresponds to M Z D . The yellow colored points in Fig. 6 are excluded by the current direct search constraints while the blue colored points are still allowed.

Constraints from B-factories and beam dump experiments
The LHCb collaboration [106] has derived upper limits on the kinetic mixing factor at 90% C.L., covering 214 MeV M Z D 70 GeV, using the LHC data set collected at an integrated luminosity of 1.6 f b −1 at √ s = 13 TeV. For M Z D 10 GeV, this search offers the strongest upper limits on g , among all other contemporary dark photon experiments, and excludes g 10 −3 for M Z D ∼ 10.6 GeV. The limit becomes slightly weaker towards higher M Z D and excludes g 2 × 10 −2 at M Z D ∼ 70 GeV. At M Z D below 10 GeV, the most stringent constraints are offered by BaBar [107], which exclude g 10 −3 in the mass range 0.25 GeV < M Z D < 10 GeV. BaBar performed the search in the e − e + → Z D γ channel, while assuming that the Z D predominantly decays invisibly.
In the current analysis, our scanned set of parameter space points have M Z D in between 2 GeV and 60 GeV, while g has been scanned up to 10 −4 . Thus, the current constraints on g derived from dark vector boson searches do not affect our parameter space. In addition to the Z D searches by LHCb and BaBar, there are numerous beam dump experiments which have also probed a light vector gauge boson. However, the searches by the beam dump experiments mostly concentrate in the mass range of the order O(M eV ). Some of these experiments are KLOE [108], MAMI [109], NA-64 [110], E141 [111], E774 [112], KEK [113], HADES [114], and MiniBooNE [115]. Some of the upcoming experiments like DarkLight [116] and APEX [117], are expected to improve upon the sensitivity to g by ∼ 2 orders of magnitude.  [118][119][120] and flavor physics experiments [121][122][123].
Light scalar boson searches in B-factory and beam dump experiments also yield exclusion contours on sin θ mix as a function of M H 1 . The E949 experiment [118] probed the kaon decay process, K + → π + νν, in the pion momentum range 140 MeV < p π < 199 MeV. These search limits have been translated to the M H 1 − sin 2 θ mix plane in [124] by re-interpreting the analysis scheme of [125]. In the context of our analysis, we use the corresponding exclusion contour shown in Fig. 1 of [124] and show the excluded parameter space points in purple color in Fig. 7. The CHARM collaboration [119] performed a search for axion like particles using a 400 GeV proton beam from CERN-SPS dumped onto a copper target. The corresponding limits have also been been translated and presented as exclusion contours in the M H 1 − sin θ mix plane in [124], which we directly use in our analysis, and the excluded parameter space points have been shown in yellow color in Fig. 7. We would like to note that the sensitivities of CHARM and E949 experiments overlap in the M H 1 250 MeV region with E949 exerting more stringent constraints below M H 1 40 MeV. Results from the search for weakly interacting massive particles by the SuperCDMS collaboration [120] has also been extracted from Fig. 1 of [124], and excludes the parameter space region corresponding to M H 1 5 GeV and sin 2 θ mix 10 −5 . These excluded points have been shown in light brown color in Fig. 7. The B-factories exert the strongest constraints on sin 2 θ mix in the intermediate light Higgs mass region, 400 MeV M H 1 5 GeV. The search for H 1 performed by LHCb collaboration [121] in the decay process : B 0 → K + π − H 1 , with the H 1 eventually decaying into a di-muon pair, excludes sin 2 θ mix 10 −7 in the mass range of 214 MeV M H 1 4 GeV. The corresponding exclusion contour has been taken from Fig. 1 of [124] and the excluded parameter space points have been shown in red color in Fig. 7. The measurements in B → KH 1 channel by BELLE [122,123] and LHCb have also been translated into limits on sin 2 θ mix in [124]. The parameter space points excluded by those are shown in sky blue color in Fig. 7. It can be observed from

Dark Matter Aspects
Having constrained a substantial part of the parameter space from flavour physics, beam dump and collider experiments, we now turn to dark matter phenomenology, which, as we shall see, will help us to put more stringent limits on our model thereby enhancing its predictability.

Prospects of Direct Detection
The direct detection experiments pose severe constraints on interactions of DM with nucleons. In the present context, we will assume H 1 to be very light, as required to enhance the self-interaction cross-section of the DM. Further, scenarios with both heavy DM and light DM (χ − ) will be discussed.
Note that, the contribution from only H 1 , H 2 mediated processes are important where v esc 544 km/s denotes the escape velocity of our galaxy, the incoming DM particle (χ − ) would not have enough kinetic energy to excite the heavier state, leading to kinematic suppression of the Z D mediated t-channel process. The differential event-rate at a detector, as a function of the nuclear recoil energy E R , is given by, where, n T is the number of target nuclei in the detector material, ρ χ − is the local density of DM halo ( 0.3GeVcm −3 ) and dσ(v, E R ) dE R is the scattering cross-section with a nucleus. Further, f E ( v) denotes the velocity distribution of the DM with respect to earth and can be related to the velocity distribution f ( v) of DM in the galactic halo where v E denotes the velocity of earth in the galactic rest frame.
We will assume f ( v) to be a Maxwell-Boltzmann distribution with velocity dispersion v 0 = 220 km/s and a cut-off set to v esc . The minimum velocity v min corresponding to a recoil energy E R to the target nucleus is given by, where m T denotes the mass of the target nucleus and µ = m T M DM (m T + M DM ) −1 denotes the reduced mass of the nucleus-DM. The interaction with a nucleus with atomic number A and charge Z, then, is given by, is the χ − nucleus cross-section at zero momentum transfer. Also, f p and f n denotes couplings with p (proton) and n (neutron) respectively. We have In the above expression λ q denotes the effective coupling of χ − with the quark q in the limit of small momentum transfer, and is given by f y q sin θ mix cos θ mix 1 M 2 where y q denotes the Yukawa coupling for quark q. f N q denotes the contribution of quark q to the mass m N of nucleon N . While the light quarks contribute to the nucleon masses directly, the heavy quark contributions to f N appears through the loop-induced interactions with gluons. These are given by, Note that, since H 1 and H 2 mediated t-channel processes contribute, the s quark content of the nucleon can be of significant importance. Following micrOMEGAs4.3.5 [127], we have used σ πN = 34 MeV and σ s = 42 MeV to determine the quark contents of the nucleon. Further, F (q) denotes the nuclear form factor corresponding to a momentum transfer q. However, since we are interested in the presence of a very light H 1 , as required for the Sommerfeld enhancement of the self-interaction cross-section, the mediator mass can be comparable to or even smaller than the typical momentum transfer of O(100) keV (for ∼ O(100) GeV DM). In such cases, use of σ T SI , as described above, overestimates the direct detection constraint on f sin θ mix . In order to account for the same, we have introduced a factor M 4 H 1 (q 2 + M 2 H 1 ) 2 and multiplied the same to the usual Helm form factor (F (k) in eq 6.3). For k max M H 1 , this additional term simply reduces to 1, while for k M H 1 this ensures the correct momentum dependence of the DM-nucleus interaction cross-section. Note that depending on the detector threshold, the minimum momentum transfer q min is different, while the maximum possible momentum transfer q max is set by v esc . Typically q min falls well below 1 MeV, and assuming H 1 thermalizes with the SM particles, q min M H 1 would contribute to the additional relativistic d.o.f during BBN and would be in tension with the relevant observations. We have modified the publicly available code DDCalc [128] to incorporate the above. and computed the constraints on f sin θ mix from direct detection experiments. This is shown in Fig. 8 (left) along with the self-interaction cross section in the colour bar.
The range of f used in the scan is from 10 −5 to 1. On the right hand side, we show the same plot but this time only the points with strong self interaction is plotted. The suitable range of f corresponding to this is shown in the adjacent colour bar. The direct detection constraints becomes more and more rigid as we increase the dark matter mass being most tightly constraining at M DM ∼ 50 GeV. Then the bound weakens gradually as we increase the mass. This fact is also reflected in the figure through the arrangement of the different exclusion lines. The line corresponding to the 50 GeV dark matter rules out the largest volume of the parameter space, whereas the 500 GeV dark matter corresponds to a looser constraint than 50 GeV (but tighter than 5 GeV) as expected.
The improved technique of calculation of direct detection constraints opens up a large part of the parameter space as opposed to the conventional calculation. A comparison between the two methods in shown in Fig. 9. In case of low mass DM  Before moving forward some important points are in order. Firstly, till now we were discussing about direct detection experiments which measures nuclear recoil when a dark matter particle hits them. Unfortunately, if dark matter mass goes down below (O(300) MeV) then nuclear recoils are not an effective method to detect dark matter since recoils energies become pretty low. The effective method is these low mass regions is to measure electron recoils. Experiments like SENSEI [131], XENON-10 [132], XENON-100 [133], SuperCDMS [134] and DarkSide-50 [135] have measured dark matter electron scattering cross sections and have put constraints on the latter. As mentioned in [136], the upper bound on dark matter electron cross section on a O(100) MeV dark matter is ∼ 10 −38 cm 2 (with form factor F DM set to unity). The dark matter electron cross section in our model goes as ∼ f 2 y 2 e sin 2 2θ mix , where y e is the electron Yukawa coupling and hence expected to be quite small. For the case of a light mediator (but M H 1 > 10 MeV to avoid BBN constraints), the dominant contribution to the dark matter electron cross section comes from the t-channel process. For M DM 1 GeV but greater that M H 1 and M e , the DM-electron cross section at small center of mass momentum is approximately given by : For a typical f that gives sizable self-interaction i.e. f ∼ O(0.1) and sin θ mix ∼ 10 −5 , a 100 MeV dark matter the dark matter has a cross section of ∼ 5.6 × 10 −50 cm 2 with electrons and hence much below the upper limit set by the DM-electron scattering experiments.
On the other hand the cross section for an very light dark matter (M DM ∼ 10 keV, and consequently much lighter than the mediator as well as the electron) is given by the following expression : 13 For even lower masses CRESST-III [130] has better sensitivity, however, we have not explored the region M DM 1 GeV.
But such low masses are beyond the reach of DM-electron scattering experiments and no limits on this cross section exists as far as the direct detection experiments are concerned. Secondly, we should be careful about the presence of the neutrino floor [137], going below which would render direct detection signals to be meaningless. Hence f sin θ mix can not be lowered indefinitely. For a 50 GeV dark matter and O(1) GeV H 1 mass, we found earlier that f sin θ mix 10 −6 . This f sin θ mix can be further lowered roughly by ∼ 24.5 times before it hits the neutrino floor for this dark matter and mediator mass.

(Thermal) Relic density of dark matter
From our discussions in the previous sections, we find that to satisfy constraints from direct detection experiments and at the same time allowing for sizable self interactions, the required value of sin θ mix is actually very small ( 10 −5 − 10 −7 , for a suitable f , while the latter is fixed from considerations of self-interactions). With such small values of mixing angle, the standard procedure for calculation of relic density (assuming prior thermalization and a subsequent freeze-out) is placed under scrutiny and demands some in depth analysis before proceeding further. In this work, we have already assumed that the kinetic mixing between Z and Z D is very small (denoted by g ). Hence, the scalar mixing angle (θ mix ) is the only possibility through which the dark sector can communicate with the standard model sector and eventually can come to equilibrium with the photon bath. The two parameters in our model that can control this thermalization effectively are λ mix (H 2 H 2 ↔ H 1 H 1 type of interactions) and sin θ mix (ZZ ↔ Z D Z D type of interactions). Although sin θ mix depends on λ mix , but considering only relevant values of the latter is not sufficient alone to give us an idea about the mixing angle required for thermalization. This is because the processes that depend solely on sin θ mix (e.g. ZZ ↔ Z D Z D ) can also play an important role in equilibrating the dark and standard model sectors. However, before EWSB (T EWSB γ ∼ 153 GeV), the scalar mixing angle had no significant role as such. Then λ mix is the only relevant parameter that can bring the two systems to equilibrium. Before EWSB, at very high temperatures, the scalar masses can be neglected. The cross section of H 2 H 2 ← H 1 H 1 is given by σ ∼ λ 2 mix 32π s . The two systems can thermalize if the rate of annihilation can exceed that of the Hubble expansion at some (sufficiently high) temperatures. Mathematically, the rate of annihilation is given by Γ = n eq σv , where n eq is the equilibrium number density of the particle under consideration and v is the relative velocity of the annihilating particles. If at some temperatures T T EWSB γ , Γ/H ∼ 1, then we can safely conclude that the dark and the visible sectors did thermalize in the early universe. Here H is the Hubble expansion rate and is given by 1.67 √ g T 2 M Pl , where g is the degree of freedom and M Pl is the Planck mass. In our discussion we have taken T ∼ 10 5 GeV. We find, that for the two sectors to be in equilibrium at such temperatures λ mix 10 −5 (Fig. 10). After EWSB, H 2 H 2 ↔ H 1 H 1 annihilations can compete with the Hubble rate if λ mix 7 × 10 −6 . However not all interactions are dependant on λ mix exclusively. For example ZZ ↔ Z D Z D depends only on sin θ mix . So even if λ mix is small, the dark and the standard model sector can come to equilibrium through the aforementioned interaction, and this is in turn can constrain sin θ mix . For temperatures ∼ 100 GeV, we find for Γ(ZZ ↔ Z D Z D ) H, sin θ mix 7.5 × 10 −6 . Hence for mixing angles below this value, the two sectors will fail to thermalize through this channel. Thus in the case where the dark and SM sectors fail to thermalize in the early universe (i.e both λ mix and sin θ mix are small) usual techniques of calculation of relic density by freeze-out mechanism will lead to incorrect results. However, as discussed above, keeping one of them large can in principle lead to successful thermalization of the two sectors.
In the previous section, we derived constraints on f sin θ mix from direct detection experiments and plotted points which were simultaneously allowed by it and also can lead to sizeable self-interactions. The range of f found from such considerations was matter mass). But, we know, Hence, we can try to visualize the part of the parameter space (in M Z D − λ mix plane, say) that can give rise to a thermal relic. Two values of f and two different dark matter masses were chosen for demonstration purpose. It is evident from Fig. 11, that with increasing f the allowed region decreases for a given dark matter mass. On the whole, light dark matters are more favoured with respect to the heavier ones in the thermal scenario. The value of the dark gauge coupling g D was however arbitrarily fixed to 0.2 in these plots. If we lower g D , the allowed parameter space decreases as expected.
In principle, g D can be fixed from considerations of correct thermal relic (for e.g see Table. 2 Finally, analogous to the plot in the previous section, we show the range of variation of relic density alongside that of λ mix for the points allowed by direct detection experiments and at the same time giving rise to sizeable self-interaction in Fig. 12. From this, we can easily conclude that light dark matter is our best bet if wish to stick to the (usual) regions of thermal relic density.
The dominant contributing channel to the thermal relic density can be broadly classified into three classes : • χ 1 χ 1 ↔ H 1 H 1 : The cross section solely depends on f (σ ∼ f 4 ). This type of channel can in general be present in simple extensions of SM (for example in cases where we have only a scalar portal mediating the dark sector and the visible sector). We have seen that the requirement of strong self-interaction among the dark matter particles pushes f towards higher values (f ∼ O(0.1)). So, this may lead to under abundance in usual circumstances. But in our model since the dark matter χ − is a Majorana particle, this type of annihilations are CP-odd and hence p-wave suppressed consequently satisfying thermal relic density. The sole controlling parameter is the Yukawa coupling f which also allows for a sizable self interaction cross section. But for general dark matter masses, correct relic density and large self-interactions cannot always be satisfied by adjusting this single parameter f , and leads to over abundance. The cross section is now governed by both f and g D (σ ∼ f 2 g 2 D ). The dark gauge coupling can be tuned to adjust the present day relic density to the observed value. Since g D does not contribute to the self-interaction of dark matter, this parameter can be freely chosen to fix the relic, while f is fixed to a value that gives strong self-interaction between the dark matter particles. The mass of the dark matter particle is also constrained now such that M DM > M Z D .
• χ 1 χ 2 ↔ SM SM (co-annihilation) : Thermal relic density in principle can also be dominated by co-annihilation cross sections of χ − (DM) and χ + . The impact of co-annihilation is determined by the mass difference of the two dark states. It is strongest if the mass splitting is small. The mass splitting is given by ∆m χ = √ 2f v D . But, from our earlier discussions we came to the conclusion that to stick to the thermal scenarios, its best to adhere to low mass dark matter. On the other hand, self-interaction demands for large values of the Yukawa Figure 12: Plots showing the variation of of relic density along with the quartic coupling, λ mix in the colour bar for points allowed by direct detection constraints and also having sizeable self interactions.
coupling f . Thus co-annihilation fails to be the dominant contributor to the relic density of dark matter if the mass of the latter is taken to be small. This can also be understood clearly as follows : the principle term controlling the thermally averaged co-annihilation cross section is the Boltzmann factor given by ∼ e M DM . Here, T F O is the freeze out temperature of the dark matter. Hence higher the dark matter mass and lower the mass splitting, the more important are these co-annihilation channels in their contribution towards the dark matter relic density. From Fig. 8, we find that for massive dark matters the upper bound on f sin θ mix is quite stringent. But for sizable self interactions f cannot be very small. This forces sin θ mix towards smaller values (from direct detection constraints) and consequently demands a small v D (since λ mix 10 −5 for a thermal scenario). But very small v D will lead to very light Z D thereby making them unsuitable for probing via colliders. Hence no such points exists in the parameter space where co-annihilation is the dominant contributor to the relic density and at the same time, it also respects all the other requirements (like direct detection, sizable self interaction, thermal dark matter and not too light Z D ).
To summarize, thermal dark matters which are strongly self interacting and also Thermal Benchmark   Fig. 8 we find that f sin θ mix 6 × 10 −9 for M H 1 ∼ O(10) MeV. The mediator should be light to allow for large self-interaction cross section but not too light so as to violate bounds from BBN ( O(10) MeV). Also, we found that for sizeable self interaction at these dark matter and mediator masses f ∼ O(0.1). Taking v D to be as low as ∼ O(1) GeV, we find λ mix 7.7 × 10 −6 and that falls below our derived limit of thermal threshold. Hence for high mass dark matter, we are almost out of points in the parameter space that satisfies all of our aforementioned desired criteria. Light dark matters are hence more favourable in our scenario. Even for a O(5) GeV dark matter, direct detection experiments force f v D 0.25. For reasons previously mentioned, with f ∼ 0.1 we find v D 2.5. Along with this, requirement of λ mix 10 −5 (thermalisability) and not too small Z D mass restricts our allowed parameter space quite heavily.
Next, as a concrete example, we present a specific benchmark in Table. 2 which illustrates that for specific choice of parameter values one can indeed have a perfectly thermal scenario satisfying all the constraints and also having a sizable self-interactions although only with small dark matter masses. 15 In the next section we will discuss about the future detection prospects of our scenario in colliders.

Future searches
From the discussions in the perspective of dark matter phenomenology, we have thus been able to convince ourselves that for a thermal dark matter with large self interactions we need very light H 1 as the mediator ( ∼ O(1 GeV)). In the following sections we will investigate about the future of prospects and detectability of our model from a collider perspective. We probe the prospect of future collider experiments and estimate their effect on the parameter space which survives the current constraints as discussed in Sec. 5. We begin by a discussion about the future search for the light scalar mediator (H 1 ).
The possibility to probe a light pseudoscalar particle from Higgs decays at the 14 TeV LHC has been studied in [141] for 2µ2τ final state. This search covered the mass range : 2M τ < M A < 2M b . Using the upper limit obtained from this analysis, expected future reach for 14 TeV LHC with integrated luminosity L = 300 f b −1 have been translated in [142] (see Fig. 7 of [142]). The projected upper limit on H 2 → H 1 H 1 branching ratios is roughly 5% [142]. For our case we find that H 2 → H 1 H 1 branching ratios lie in the range 10 −3 to 10 −6 , which will be beyond the reach of HL-LHC in 2µ2τ final state.
In [143], the expected future reach at the 14 TeV LHC has been studied for H 2 → AA → 4b and H 2 → AA → 2b2τ final states, where the Higgs boson is produced in association with a W or Z boson. This analysis is sensitive for M A > 10 GeV and 4b final state is more promising than the 2b2τ channel. The potential of an exotic Higgs decay search for H 2 → XX → 2b2µ (X = H 1 /A) in the mass range of 15 to 60 GeV has been presented in [144]. It is found that Br(H 2 → 2X → 2b2µ) can be constrained at the few ×10 −5 level at the HL-LHC. Both these analyses [143,144] are not sensitive for the parameter space of our interest (M H 2 < 10 GeV).

Future prospects at ILC
Uncertainties in Higgs boson couplings to various SM final states from a combination of global fit of Higgs coupling measurements at the ILC have been presented in [145]. Combination of the results of Higgs coupling measurements from the 300 fb −1 run of LHC, and 500 GeV run of ILC, may lead to an error of 0.3% in HZZ coupling [145]. Correspondingly, we test the impact of this constraint on our parameter space by choosing those points which generate µ ggF ZZ * in the range of 0.997−1.003. Such parameter space points have been shown in orange color in Fig. 13. The blue colored points in Fig. 13 correspond to the entire set of parameter space scanned by us. It can be observed that the projected reach of ILC is sin θ mix 7 · 10 −3 which is roughly 1 order of magnitude stronger than the current LEP results.

Exclusions from future SHiP and LZ experiments
Continuing our discussion of Sec. 5.4, we show the projected reach of two relevant future experiments, namely, SHiP [146] and LZ, in Fig. 14. The projected sensitivity of SHiP is remarkably strong when compared to the current experiments. For M H 1 ∼ 1 GeV, the current limits exclude sin θ mix above ∼ 10 −3 , whereas, SHiP is projected to probe  (Fig. 7). The green colored points will be within the projected reach of SHiP experiment while the brown colored points could be probed by the proposed LZ.
until sin θ mix ∼ 10 −5 . LZ is expected to gain effectiveness in the O(GeV ) region, and is expected to improve upon the existing sensitivity by around an order of magnitude.
As seen from the findings of the preceding three subsections, a hunt for a light scalar mediator in future colliders seem to be challenging. Our model however also possesses an O(1) GeV Z D and dark matter. Hence, we can instead try searching for this dark photon in future colliders. It is important to note that the analysis performed in the next two subsections would comply only in the case of a promptly decaying Z D boson. The total decay width of Z D is proportional to 2 g . A value of 2 g 10 −12 results in a lifetime which is compatible with the regime of prompt decay. Smaller values of g ( 10 −6 ) results in decay lengths of the order of ∼ 10 −5 m or higher, resulting in a late decaying phenomena. The generic collider features of a late decaying Z D has been discussed in Sec. 7.6.

H 2 → Z D Z D → 4µ at HL-LHC
A detailed search analysis for Z D is presented (M Z D = 2 − 12 GeV), using the process pp → H 2 → Z D Z D → 4µ final state, in the context of 14 TeV run of LHC at an inte-grated luminosity of 3000 f b −1 . At the LHC, the Higgs boson is dominantly produced via the gluon fusion mode (ggF ). The ggF mode overshadows the other modes of Higgs production, such as vector boson fusion (V BF ), associated production with b quarks (bbH 2 ) and associated production with top quarks (ttH 2 ), and therefore, in the current analysis, we consider only the ggF mode of Higgs production.
The signal sample involves the process : gg → H 2 → Z D Z D → 4µ. Dominant backgrounds arise from electroweak 4l production, and the H 2 → ZZ * decay process. The gg → H 2 samples have been generated using MadGraph5 [147], with M H 2 fixed at 125 GeV. The subsequent decay process, followed by showering and hadronization is performed by Pythia-6 [148]. The 4l background has been generated by MadGraph5 at the partonic level, while Pythia-6 has been used for showering.
The final state muons are required to have transverse momentum, p T > 2.6 GeV and must lie within a pseudorapidity range of |η| < 4.0. Isolation condition requires the sum of transverse momentum of other tracks, excluding the leading four muons, within a cone of radius ∆R = 0.2 around the muon to be less than 30% of the p T of the muon. Furthermore, if the ∆R between the muon and the reconstructed jet is less than 0.4, then the muon is removed. Event selection requires exactly 4 muons to be present in the final state. Before moving on to discuss the choice of selection cuts, we discuss the kinematic distribution of the four muons in the final state. In this context, we choose three different values of M Z D = 2, 5, 8 GeV to generate the signal. The p T distribution of the four muons from the background and the benchmark signals are shown in Fig. 15, where, the 4l and H 2 → ZZ * backgrounds have been described by brown and orange colored lines, respectively. We would like to mention that the muons have been p T ordered, and therefore, p 1 T corresponds to the highest p T muon while p 4 T represents the lowest p T muon. The blue, green and black colored lines correspond to the p T distribution of the signal samples with M Z D = 2, 5, 8 GeV, respectively. The signal samples and the H 2 → ZZ * background yields a similar p T distribution, while the 4l background peaks at very low values of p T . Taking motivation from a similar search analysis by the ATLAS collaboration [93], and the p T distribution of the 4l background, we demand the highest p T muon to have p T > 20 GeV. The second and third leading p T muons are required to have p 2 T > 15 GeV and p 3 T > 10 GeV, respectively. The angular distributions of the final state muons also provide an additional control to filter out the signal events. In this context, we show the ∆R distribution between all some of the possible final state muon pairs in Fig. 16. Before specifying the selection cuts on the ∆R variables, it is important to discuss the criteria for formation of same flavor opposite sign (SFOS) muon pairs. The four muons in the final state are required to form two SFOS pairs. Each quadruplet of muon per event can result in two separate combinations of di-SFOS pairs. The di-SFOS pair which results in smaller difference of SFOS muon invariant masses is chosen to be the correct di-SFOS pair. The SFOS pair with invariant mass closest to the Z boson mass is referred to as the leading SFOS pair and its invariant mass will be represented by M 12 ., while M 34 represents the invariant mass of the sub-leading SFOS pair. Henceforth in this subsection, we will refer to the muons in the leading SFOS pairs as µ 1 and µ 2 , while muons in the other SFOS pair will be referred to as µ 3 and µ 4 . Now, we go back to the discussions on the ∆R distributions shown in Fig. 16. It can be observed that the backgrounds mostly peak around ∆R ∼ 1 − 2.5 region, while the signal events peak at two distinct regions, ∆R ∼ 0 and ∆R ∼ 3. The peak in the ∆R ∼ 0 arise from muon pairs which originate from the same Z D pair, and therefore, mostly correspond to the same SFOS pair, while the ∆R ∼ 3 peak manifests from the muons belonging to different SFOS pairs. Deriving motivation from this observation, we impose a lower cut on ∆R between muon pairs from separate SFOS pairs, ∆R 13 , ∆R 14 , ∆R 23 , ∆R 24 > 2.
In addition, the invariant mass of the SFOS pairs, M 12 and M 34 are required to lie within the range 0.88 − 20 GeV, and the ratio M 34 /M 12 is required to be greater than ∼ 0.85, following [93]. The invariant mass of the four isolated muons (M 4µ ) is also required to lie within M 4µ = 120 − 130 GeV. We would like to mention that events with any of the SFOS pairs having invariant mass in the range of J/Ψ resonance Selection cuts (a). Exactly 4 muons in final state.  Table 3: Selection cuts for the cut-based analysis in the 4µ final state, following [93].  Table. 3. The tt and bbll backgrounds were also evaluated, and no events survived the selection cuts.
We   Table. 4. Upper limits are obtained on the production cross-section of SM-like H 2 through ggF mode times the branching fraction H 2 → Z D Z D at 5σ and 2σ, assuming zero systematic uncertainty, for the case of 14 TeV LHC operating an at integrated luminosity of 3000 fb −1 , as shown in Fig. 17. The branching ratio for Z D → µ + µ − has been taken from [90]. The upper limits for HL-LHC presented in Fig. 17 is roughly ∼ 10 − 15 times stronger than the current bounds shown in Fig. 6. We conclude this subsection by evaluating the prospect of exclusion/discovery of the allowed benchmark point (BP1) of Table. 2 by HL-LHC using the search limits derived in this subsection. BP1 furnishes a value of Br(H 2 → Z D Z D ) ∼ 5 · 10 −6 for M Z D ∼ 2.4 GeV, while σ(gg → H 2 ) attains a value of cos 2 θ mix · σ(gg → H 2 ) SM , where, σ(gg → H 2 ) SM corresponds to the SM value of H 2 production cross-section in the ggF mode. At NNLO+NNLL level, σ(gg → H 2 ) SM = 39.56 +7.32% −8.38% fb [149][150][151] for the 14 TeV run of LHC. Thus, σ(gg → H 2 → Z D Z D ) attains a value of ∼ 0.2 for √ s = 14 TeV. A comparison with the corresponding upper limit derived in Fig. 17 reveal that BP1 could be marginally probed at HL-LHC 16 . These searches have the potential to yield strong exclusion limits in the context of future high energy/ high luminosity colliders. Directly produced Z D process, pp → Z D → µµ, could be another viable mode of probing Z D . However, the signal is marred by a huge continuum Drell-Yan background. Consequently, the search strategies would involve an efficient treatment of low p T muons. Some recent studies have focused on the case of such low p T muons through scouting techniques [152,153]. A more precise understanding of the background and derivation of further improved scouting techniques might help in alleviating pp → Z D → µµ as an important probe of Z D in the future runs of LHC.

H
The SM-like Higgs boson, H 2 , can also undergo decay into a pair of χ + χ + , which can result in a 4µ+E / T final state, through cascade decay via In this subsection, we present a search strategy for Z D in the H 2 → χ + χ + → Z D χ − Z D χ − channel, in context of a 14 TeV LHC with L = 3000 fb −1 . Similar to the case of Sec. 7.4, we only consider the ggF mode of Higgs production. gg → H 2 → χ + χ + → Z D χ − Z D χ − constitutes the signal, where, gg → H 2 samples have been generated with MadGraph [147], while Pythia-6 [148] has been used to perform the cascade decay and showering. The benchmark point shown in Table. 2, with M χ + = 5.8187 GeV and M χ − = 2.4 GeV, has been chosen to perform this search.
The muon isolation and selection criteria specified in Sec. 7.4 is applied here as well, and an event is required to have exactly four isolated muons in the final state. We show the p T distribution of the final state muons, corresponding to the relevant backgrounds (same color code as Fig. 15) and two benchmark signal events (M Z D = 2.0 GeV and M Z D = 2.8 GeV), in Fig. 18, where, p 1 T represents the highest p T muon and p 4 T represents the lowest p T muon. When compared to the case of Sec. 7.4, the p T distribution of the signal events peak at a relatively smaller value on account of a residual E / T . Consequently, we apply slightly weaker trigger cuts on the muon p T . The muon with the highest p T is required to have p T > 15 GeV, while the second Figure 18: Transverse momentum (p T ) distribution of the final state muons at the partonic level. The blue and green colored lines correspond to the signal event generated with different Z D masses, M Z D = 2, and 2.80 GeV, respectively. The orange and the brown colored line represents the p T distribution of the ZZ * → 4l and electroweak 4l background.
(third) leading muon is required to have p T > 10 GeV(5 GeV). Here again, the angular variables, ∆R, provides an additional control in improving upon the signal significance. In this context, we show the ∆R max and ∆R min distributions in Fig. 19, following the color code of Fig. 18. ∆R max (∆R min ) represents the maximum (minimum) value of ∆R between all possible muon pairs. It can be observed from Fig. 19 that the ∆R max distribution for the 4l background falls off earlier than the signal samples. Consequently, we impose ∆R max > 2.0. In addition, a lower limit on ∆R min > 0.5 is also imposed. The construction of SFOS pair is done following the strategy specified in Sec. 7.4. An event is  Table. 5. The HL-LHC upper limits on σ(gg → H 2 ) × Br(H 2 → χ + χ + → χ − Z D χ − Z D ) corresponding to a signal significance of 5σ and 2σ are shown in Table. 6. For the sake of comparison, the upper limits on σ(gg → H 2 → Z D Z D ), derived in Sec. 7.4, are also shown in the same table.
The case of H 2 → Z D Z D (discussed in Sec. 7.4) , furnishes stronger limits as compared to the current case, and the reason can be attributed to the possibility of invariant mass reconstruction of H 2 (M 4µ ) in the previous case. We would like to note that in the case of H 2 → Z D Z D , we had imposed a selection cut on M 4µ and had restricted it within 120 − 130 GeV, which was extremely efficient in filtering out the background.     cτ > 10 −4 m. Z D is one such viable LLP candidate within the framework of the U (1) D model considered in this analysis. The proper decay length of Z D has an inverse square dependence on the kinetic mixing factor g , and can become a potential LLP candidate for smaller values of g . For example, at M Z D ∼ 2 GeV, a value of g (10 −5 − 10 −6 ) will result in the Z D to become a LLP (see Fig. 2 of [154]). Before proceeding ahead with the details of the analysis, we lay down a brief description of the existing segmentation of the ATLAS detector. The ATLAS detector can be broadly classified into four different segments, viz., the tracker region, the electromagnetic calorimeter (ECAL), the hadronic calorimeter (HCAL) and the muon spectrometer. The tracker region can be further subdivided into three major segments. The first segment involves the Pixel detector, which contains three sub-layers at a radii of 33.25 mm, 50.5 mm, and 88.5 mm. The next segment within the tracker is the SCT which can be further classified into three segments with radii 299 mm, 371 mm, and 443 mm. The final segment in the tracker is the TRT which spans Figure 20: The vertical and horizontal axes represent the radii from the beam axis and proper decay length, respectively, of the LLP. The color palette corresponds to the percentage of events which will undergo decay in the corresponding segments, for the case of 14 TeV LHC with an integrated luminosity of 3000 fb −1 . The segmentation along the vertical axis corresponds to the existing geometry of the ATLAS detector. the radii 554 mm − 1082 mm. The ECAL segment lies roughly in between a radii of ∼ 1300 mm−2100 mm, while the HCAL scans across a radii of ∼ 2285 mm−3815 mm. The final segment of the ATLAS detector, the muon spectrometer, extends from a radii of ∼ 4100 mm all the way up to 10000 mm.
The LLP can undergo decay in different segments of the LHC, based on its kinematic distribution and its decay length. The decay length is given by l d = βγcτ , where, β is the boost of the particle and is defined as the ratio of its velocity (v) to the speed of light (c): β = v/c, whereas, cτ is the proper lifetime of the particle in its rest frame, and γ is the relativistic factor defined as : γ = 1 1 − β 2 . It can be inferred from the form of the decay length that identification of the secondary vertex within the detector (which corresponds to the point of decay of the LLP) and measurement of the boost factor, can also be used to estimate the proper lifetime of the LLP. Within a typical detector, if N 0 be the number of long-lived particles produced with a proper life time τ i and a mean life-time of τ , then the exponential distribution gives the total number of decay events, N = N 0 e −τ i /τ . Consequently, the fraction of long-lived particles undergoing decay in different segments of the detector will be a characteristic reflection of the proper decay length.
To visualize this effect, we used Pythia-6 to simulate the signal pp → H 2 → Z D Z D with Z D being a LLP. Event samples were generated by varying the decay length of Z D over the range from 0.5 mm to 10000 mm, and 20, 000 events were generated for specific choice of decay length. We show the proper decay length (in mm) along the horizontal axis and the radii (in cm) along the vertical axis in Fig. 20 Table 7: Percentage of events undergoing decay in different segments of the detector corresponding to the ggF and ZH 2 production modes of H 2 , for the case of 14 TeV LHC. The difference is the manifestation of different kinematics for the two cases.
detector geometry), with the color palettes being representative of the percentage of events decaying within the concerned segment. The black colored regions correspond to the void regions within the detector and also represents those segments which do not register an event decay. The event samples for Fig. 20 have been simulated for the case of 14 TeV LHC. It can be observed that particles with proper decay length up to ∼ 10 cm mostly decay before reaching the Pixel detector. Particles with large proper decay length can be observed to have a fairly uniform probability of decaying throughout the detectors. We would like to note that, in order to correctly simulate the effect of a smeared distribution, we generated events with proper decay lengths up to 20000 mm and then truncated the horizontal range in Fig. 20 at 10000 mm. Observation of the LLP at different segments of the detector can be instrumental in motivating the design of segment-specific search analyses. Reconstruction of the same type of particle at the tracker and the ECAL involves vastly different strategies, and different search analyses would be required at these two different segments. Under such circumstances, having a generic idea about the probability of decay within a certain segment of the detector can prove instrumental in performing more specific and focused searches. We perform a similar study for the case of a 27 TeV LHC machine as well, assuming the current geometry of ATLAS detector, and the corresponding results are represented in Fig. 21. Here, we show the total number of gg → H 2 → Br(H 2 → Z D Z D ) events, assuming σ(gg → H 2 ) = 140 pb, Br(H 2 → Z D Z D ) = 10 −5 , and an integrated luminosity of 15 ab −1 , in the color palette. The horizontal and vertical axes in Fig. 21 correspond to the similar quantities in Fig. 20. A higher boost in the case of 27 TeV collision results in an upward shift of color pattern in Fig. 21 as compared to the previous case.
The preceding discussion assumed the Higgs production in the ggF channel. The production of Higgs in other production channels will result in alteration of the segmentwise decay fraction, owing to different boosts in the transverse direction. In this respect, we briefly explore the case of ZH 2 production, and the LLP being produced from decay of H 2 , in the context of 14 TeV LHC. We present a comparison between the fraction of events undergoing decay within various segments of the detector, for the case of H 2 being produced via ggF mode and in ZH 2 mode, in Table. 7.

Summary and Conclusion
In this work, we have explored an U (1) D -gauge extension of Standard Model from a phenomenological perspective. Our model has two Majorana fermions to render the theory anomaly free. The lightest of them serves as a dark matter candidate. As generic to any gauge theory, our model has an extra neutral Z like boson (Z D ) and since we have also agreed to employ Higgs mechanism to make the particles massive, there also exists an extra neutral scalar (H 1 ) besides the usual Standard Model Higgs as well. Motivation of studying strong dark matter self interactions led us to restrict ourselves only to light scalar mediators (10 MeV M H 1 10 GeV). Here we have studied the specific case where Sommerfeld enhancement is mediated only through one light mediator (H 1 ). In general, both the Z D and H 1 could have been light and led to large self interactions, but very light Z D would have been devoid of any collider signatures.
For simplicity, throughout this work we have kept the kinetic mixing between Z and Z D to small values. However it plays a role in the study of collider signatures arising from the prompt decay of Z D . The primary goal of this work was to investigate if there exists a substantial parameter space for a thermal dark matter with large selfinteractions. The latter is desirable since it solves some small scale structure issues as already mentioned earlier. The range of M H 1 was hence motivated by the requirement of large self-interaction cross sections. We have systematically checked constraints on scalar mixing angle (as a function of M H 1 ) arising from LEP data, B-factories and beam dump experiments. We found that LEP data rules out sin θ mix 0.2 for the range of M H 1 considered in this work. Higgs signal strength measurements on the other hand constrained sin θ mix to 0.1 in our parameter space of interest. Beam dump and flavour physics experiments put a tighter bound of sin θ mix 10 −4 on almost the whole of our parameter space. LHC analyses mainly dealt with Standard Model Higgs decaying to 4l channels and provided us with a probe to the otherwise small gauge kinetic mixing ( g ).
From the point of view of dark matter phenomenology we however have the most stringent constraints on the scalar mixing angle (for a given Yukawa coupling f ). The Yukawa coupling f is the sole controlling parameter that gives rise to sizable self interactions while direct detection experiments tend to put limits on f sin θ mix . Due to the presence of light mediators (H 1 ), the standard calculation of DM-nucleon cross section had to be refined using momentum dependent propagators. On the other hand, we have shown that condition of thermalization of the dark and visible sectors sets a lower limit to λ mix . This competes with the upper bound on λ mix that arises from the direct detection experiments (for a given v D and a suitable f which will give rise to sizable self interactions). Our findings suggest that such points on the parameter space satisfying both the limits are very rare and automatically pushes us to low dark matter masses (∼ O(1)) GeV where constraints from direct detection is considerably weaker. Alternatively, we can also try to do away with the usual thermal dark matter scenario and probe into other exotic mechanisms like dark freeze-out [155] and freeze-in [156][157][158] to help us get the correct relic density. λ mix will have no lower limit in these cases and hence a large portion of the parameter space can be recovered.
As discussed earlier, we found that it is very difficult to probe the MeV scale dark H 1 in our model via future collider experiments. Our model however also predicts an extra dark gauge boson Z D which may be probed at future HL-LHC. We have presented the reach of Z D from H 2 decay via 4µ final state and obtained projected upper limits on σ H 2 × Br(H 2 → Z D Z D ) for HL-LHC. We have also looked for the scenario where H 2 decays into a pair of χ + (χ + → Z D χ − ) and finally results into 4µ+ MET signal and found that this channel gives weaker limit compared to H 2 → Z D Z D → 4µ channel. For very small values of g , it is observed that Z D becomes a long lived particle (LLP) and different search strategies are required for the LLP scenarios. This non-traditional prospects of LLP (Z D ) at 14 TeV and 27 TeV high luminosity LHC have also been discussed.
Thus, solely from a data driven perspective, in a general and minimal U (1)-gauge extension of Standard Model we were able to restrict the presence of a heavy thermal dark matter candidate and this motivates (quantitatively) future explorations along these uncharted avenues.