Scope of 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 flavor 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 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 effective operator approach to the direct detection is inapplicable here and these constraints have been re-evaluated 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 ≲O2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \lesssim \mathcal{O}(2) $$\end{document} 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 favored. Lastly, future collider reach of such a simplified scenario has also been studied in detail.


JHEP05(2019)177
also [50,51] for a review. Although this requires a very large self-interaction cross-section, it remains consistent with the present constraints on DM self-interaction [52][53][54][55]. 4,5,6 Further, it has been argued that a velocity dependent self-interaction cross-section is favoured [62]. In the presence of light mediators, it is possible to generate large (velocity dependent) selfinteraction among DM particles, thanks to the Sommerfeld enhanced self scattering crosssection [51,[63][64][65][66]. 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, flavor 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 section 2 we describe the minimal U(1) D model and its particle content in detail. In section 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, flavor physics and collider) on our parameter space of interest is described studied in section 4 and section 5. Section 6 deals with bounds from dark matter phenomenology and discusses its implications. In section 7 we have presented a detailed study about the future prospects of our scenario from a collider perspective. Finally we summarize and conclude in section 8.
2 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 represents left-chiral Weyl 4 See also [50] and references there for a compilation of constraints from bullet-cluster [56][57][58], 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 [59]. When interpreted as due to DM self-interaction, this leads to σ m 1.5(3) cm 2 gm −1 for contact (long-range) interactions [60]. However, the data has been argued to be also consistent with standard collisionless DM [61]. 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 [50]. 7 An Abelian gauge boson can also get massive via Stuckelberg mechanism [67][68][69][70], which essentially assumes the neutral scalar component to be very heavy and therefore decoupled from the spectrum. However, we will not consider that possibility here.
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 gauge group U(1) D . The Dark Sector interacts with the SM sector via the following gauge invariant terms.
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 [71] in the lagrangian, which can possibly be generated due to highscale 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 JHEP05(2019)177 G SM × U(1) D . 8 Around the minimum of the scalar potential, these scalar fields can be described as, (2.6) and Here, ϕ R (ϕ I ) and h (σ) are the CP-even (CP-odd) parts of φ and H respectively and σ + is a complex scalar. In the unitary gauge, where the low energy degrees of freedom (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,

13)
H 2 = sin θ mix ϕ R + cos θ mix h. (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 electroweak 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, The kinetic term is canonical in the basis V µ = ζV µ , with 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,

JHEP05(2019)177
In this context, we have used the approach described in refs [72,73]. We have used SARAH and SPheno [74,75] 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 electroweak precision data in general, constraints g [76,77]. 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 electroweak symmetry breaking (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 above expression, 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.

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 self-interacting 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 eigen- value M − f vϕ is positive. The corresponding mass eigenstate is • 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.
Following these two points of view, we consider 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 [66]. 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.

JHEP05(2019)177
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 [66,78,79]: . Analytical results can also be obtained in the resonance region in between by approximating the Yukawa potential by Hulthen potential [80].
σ 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 [81,82]. From the X-ray and the lensing observations of Bullet cluster we further have σ self /M DM 1 cm 2 gm −1 [58] 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 figure 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 figure 2 (left). This was however done for a fixed α D (≡ f 2 /4π) = 0.001. The right hand panel of figure 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 [83,84]. This is because, BBN JHEP05(2019)177  . 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 [56][57][58] 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.

JHEP05(2019)177
places strong upper bounds on presence of relativistic particles when 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 self-interaction 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 g . 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 . The scan range is inspired by the region of parameter space shown in figure 2 (b) which corresponds to those sectors which facilitate a dark matter with moderate to strong self-interactions. It can be observed from figure 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 (eq. 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)

Constraints
The U(1) D parameter space discussed so far gets constrained by a multitude of experimental search results.  shown in figure 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 studied 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 [86]. In the context of our analysis, the normalized H 1 V V coupling (ζ) is directly proportional to sin θ mix , and the upper limits derived in [86] [86], 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 [87], where κ is the ratio of production cross-section of H 1 in the Higgs strahlung process to its SM value. The yellow colored points represents the allowed region of parameter space. study by OPAL collaboration [87], 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 Higgs strahlung process to that of SM Higgs production in the Higgs strahlung 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 figure 4. We find that the upper limits from [87] 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 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 10 It may be noted that decay channels originating from the light scalar Higgs (MH 1 2 GeV) require a careful treatment owing to the uncertainties associated with partial decay widths of H1 into hadronic channels. Within 2mπ MH 1 2 GeV, the uncertainties associated with theoretical calculations remain significant [91] (in the second printing (1990) by Perseus Books in the collection Frontiers in physics, no. 80, a number of errors and omissions are corrected and the references at the end of each chapter are updated. A paperback reprint of the 1990 edition has been published in 2000), 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 H1 → γγ, e + e − , µ + µ − from [92]. 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), and the relevant Higgs boson decay modes:

JHEP05(2019)177
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 Model /σ H 2 SM in eq. (5.1) can be approximated 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.

JHEP05(2019)177
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 most recent Higgs signal strength constraints, tabulated in table II and table III of [94]. 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 figure 5. We would like to note that the parameter space points corresponding to figure 5 (left) have been generated from the parameter space specified in eq. (4.3), except for f (f 0.1) and M (M 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 figure 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, , and is therefore, proportional to λ mix v H . It can be observed from figure 5 (left) that the current constraints from Higgs signal strength measurements allow the SM-like Higgs (H 2 ) to have 15% of non-SM branching fraction.
In the low g limit, ratio of Br(    [95]. This search probed the mass range of 1 GeV ≤ M A ≤ 60 GeV, and furnishes the strongest bounds for lower values of M A or M H 1 (in the range 1 to 8 GeV). We show these limits as a solid black colored line in figure 6, where the allowed parameter space points have been plotted in the σ H 2 /σ H SM × Br(H 2 → H 1 H 1 ) plane. It can be observed that the current LHC limit excludes a fraction of our parameter space points in the low H 1 mass region, M H 1 2.5 GeV, however, it would require an improvement of roughly ∼ 3 orders of magnitude in order to become sensitive enough for the M H 1 4 GeV region. The ∼ 36 f b −1 dataset was also used by ATLAS to probe the 2b2µ final state produced via H SM → SS → 2b2µ [96] and 4b final state originating via H SM → SS → 4b [97], 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 [98]. 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τ [ [105,106], 2τ 2b [107].
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 [95]. 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 [95] 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 [95]) JHEP05(2019)177 on our parameter space in figure 7. 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 figure 7 are excluded by the current direct search constraints while the blue colored points are still allowed. These search limits exclude sin θ mix upto ∼ 0.006 for M Z D ∼ 1 − 2 GeV. The impact becomes less severe at M Z D 4 GeV, where sin θ mix is excluded upto ∼ 0.01

Constraints from B-factories and beam dump experiments
The LHCb collaboration [108] 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 [109], 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 [110], MAMI [111], NA-64 [112], E141 [113], E774 [114], KEK [115], HADES [116], and Mini-BooNE [117]. Some of the upcoming experiments like DarkLight [118] and APEX [119], are expected to improve upon the sensitivity to g by ∼ 2 orders of magnitude.
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 [120] 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 [126] by reinterpreting the analysis scheme of [127]. In the context of our analysis, we use the corresponding exclusion contour shown in figure 1 of [126] and show the excluded parameter space points in purple color in figure 8. The CHARM collaboration [121] 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 translated and presented as exclusion contours in the M H 1 − sin θ mix plane in [126], which we directly use in our analysis, and the excluded parameter space points have been shown in yellow color in figure 8. 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 [122] has also been extracted from figure 1 of [126], and excludes the parameter space region corresponding to M H 1 5 GeV and sin 2 θ mix 10 −5 .

Constraints from cosmology and astro-physics
Further, constraints on light mediators can come from Big Bang Nucleosynthesis (BBN) [84] and Cosmic Microwave Background Radiation (CMBR) [2,128,129]. The effect on BBN comes from mainly due to the contribution of the light mediator ( 1 MeV) to the energy density at the onset of BBN, as well as from the late decay of the mediator during the BBN. For our context, we restrict the mediator mass 10 MeV. When the decay into µ, π JHEP05(2019)177 are kinematically forbidden, such a mediator decays mostly into e + e − , and the constraints have been discussed in refs. [92,[130][131][132][133]. The effect on CMB (on the MeV mediators relevant for our context), is mainly due to the enhanced DM annihilation rate during the recombination epoch [128]. However, note that, in our case DM annihilating into a pair of light Higgs is a p-wave process, and thus, would have additional velocity suppression during the recombination epoch. Also, an abundance of the light mediator can lead to late kinetic decoupling of the DM, and thus can lead to a reduction in the small scale power, thus attracting constraints from Lyman-α observations [134][135][136][137]. The presence of light (late decaying) mediators can also enhance the rate of supernovae cooling and such particles may be constrained from the same [138].

Dark matter aspects
Having constrained a substantial part of the parameter space from flavor 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 since 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 as 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 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 DM-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 [140], 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.

JHEP05(2019)177
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 to the right hand side of eq. (6.3) [141]. 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 [142] to incorporate the above. We have then used this code to compute a Poisson likelihood with the predicted number of events for a specific direct detection experiment with a particular DM benchmark, and compare that with the likelihood for background-only hypothesis to constrain the model parameters f sin θ mix . This is shown in figure 9 (left) along with the self-interaction crosssection 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 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 figure 10. In case of low mass DM (2 GeV M DM 4 GeV), CDMSLite [143] puts the dominant constraint. 13 For M DM 5 GeV onwards Xenon-1T [139] puts forward the most stringent limits. 14 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 [145], XENON-10 [146], XENON-100 [147], SuperCDMS [148] and DarkSide-50 [149] have measured dark matter electron scattering 13 For even lower masses, CRESST-III [144] has better sensitivity, however, we have not explored the region MDM 1 GeV.
14 Note that Lux [31] and Panda-II [32] also lead to comparable constraint as Xenon-1T [139]. cross-sections and have put constraints on the latter. As mentioned in [150], 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: σ e f 2 y 2 e sin 2 θ mix 8 ln

JHEP05(2019)177
For a 100 MeV dark matter, cross-section with electrons turn out to be roughly ∼ 5.6 × 10 −50 cm 2 for f ∼ O(0.1) and sin θ mix ∼ 10 −5 . These values of f and sin θ mix also gives rise to sizeable self-interaction as discussed earlier.
Secondly, DM direct detection experiments are not sensitive once the cross-section becomes comparable to the neutrino nucleon coherent scattering cross-section [151]. This imposes a lower limit on constraining f sin θ mix . For example, with a 50 GeV dark matter and O(1) GeV H 1 mass, we found earlier that f sin θ mix 10 −6 . With future direct detection probes, for this particular choice of parameters, this constraint can be further lowered by roughly one order of magnitude before it hits the neutrino floor. With lighter mediator mass, however, f sin θ mix can be lowered further (by a few orders) before it merges with the neutrino cross-section (see figure 10).

(Thermal) relic density of dark matter
From the discussions in the previous sections, we find that in order to satisfy constraints from direct detection experiments and allow sizeable self-interactions at the same time, the required value of sin θ mix is actually very small ( 10 −5 − 10 −7 , for a suitable f , where 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 SM sector and eventually can come to equilibrium with the photon bath. At high temperature, well above EWSB, the most important interaction leading to effective thermalization of the DM and the SM sectors is governed by the coupling λ mix (H 2 H 2 ↔ H 1 H 1 type of interactions). Note that, sin θ mix depends on vevs via eq. (2.10), and thus would be negligible well above EWSB. 15 The cross-section of H 2 H 2 ↔ H 1 H 1 is given by σ ∼ λ 2 mix 32π s , where we have neglected the thermal masses of the scalars and treated these as relativistic species. The two systems can thermalize if the rate of interaction can exceed that of the Hubble expansion at a (sufficiently) high temperature. 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 γ ( 160 GeV) [153,154], Γ/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 (figure 11). Hence for values smaller than this, the two sectors will fail to thermalize and usual techniques of calculation of relic density by freeze-out mechanism will lead to incorrect results.
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 10 −3 f 1. This in turn constrains the scalar mixing angle (depending on the dark 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 figure 12, 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. 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 figure 13. From this, we can easily conclude that light dark matter is our best bet if one wishes 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:

JHEP05(2019)177
Thermal Benchmark  satisfying thermal relic density (see, for example, table 2). In such cases, for points with correct relic density, we have f ∼ M DM (in TeV)/1 TeV (assuming only t and u channels are the dominant ones, which is the case here because of small values of λ φ ). However, for general dark matter masses, correct relic density and large selfinteractions cannot always be satisfied by adjusting this single parameter f , and may lead to over abundance.
• χ 1 χ 1 ↔ Z D Z D : if it is over abundant, in order to obtain the correct relic abundance, we in principle can resort to a new channel employing the extra dark gauge boson Z D . This type of channels are also a generic feature of minimal extensions of SM by a new local gauge group. The cross-section is now governed by both f and g D (σ ∼ f 2 g 2 D or g 4 D , depending on which channel has the dominant contribution). 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 selfinteraction between the dark matter particles. The mass of the dark matter particle is also constrained now such that M DM > M Z D .
In our case, from figure 2, we however find that for M DM 10 GeV and M H 1 10 MeV, those points with sizeable self-interactions (and consequently with f ∼ O(0.1)), the χ 1 χ 1 ↔ H 1 H 1 cross-section is itself too large rendering the dark matter under abundant. Thus, having this extra channel depending on g D will not be helpful in practice to fix the relic density.
• χ 1 χ 2 ↔ A B (co-annihilation): thermal relic density in principle can also be dominated by co-annihilation cross-sections of χ 1 (DM) and χ 2 annihilating to A and B.
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 . On the other hand, self-interaction demands for large values of the Yukawa coupling f . Thus co-annihilation fails to be the dominant contributor to the relic density of dark matter in our case.

JHEP05(2019)177
To summarize, thermal dark matters which are strongly self-interacting and also respect constraints from direct detection experiments, as well as have observable collider signatures, should be light, i.e., O(5) GeV. To get a better feel of the numbers let us consider a dark matter of order O(500) GeV. From figure 9 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 crosssection but not too light so as to violate bounds from BBN ( O(10) MeV). Also, for sizeable self-interaction at these dark matter and mediator masses, we found that 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 (thermalizability) 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 sizeable self-interactions although only with small dark matter masses. 16 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 section 5.

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 [157]. 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 [157]. 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 figure 14. The blue colored points in figure 14 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 section 5.4, we show the projected reach of two relevant future experiments, namely, SHiP [158] and LZ, in figure 15. 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 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 section 7.3.4.

Future prospects at HL-LHC
We begin this subsection by a discussion of the future prospects of the light scalar mediator (H 1 ). We also perform a search for Z D in the H 2 → Z D Z D → 4µ channel and derive upper limits on σ ggH 2 /σ ggH SM × Br(H 2 → Z D Z D ) for the case of HL-LHC. We also present a search strategy to probe Z D in the 4µ + / E T final state, produced via cascade decay of χ + χ + pair produced from a H 2 . Finally, we conclude this subsection by a discussion of the prospects of a late decaying Z D at the HL-LHC and HE-LHC.   [97]. Apart from these analyses, CMS has also looked for direct production of a light pseudoscalar Higgs boson in the dimuon decay channel using LHC 7 TeV data [159] and the limits are not very strong yet. For lower values of M A or M H 1 (in the range 1 to 8 GeV) much stronger bounds have been obtained by ATLAS in 4µ channel [95] (shown as black dashed line in figure 6). Rest of the channels like 2µ2b [96,101], 2µ2τ [100-102], 2τ 2b [99, 107], etc. are sensitive for higher mass range (typically 15 GeV < M A < 62.5 GeV).

JHEP05(2019)177
The possibility to probe a light pseudoscalar particle from Higgs decays at the 14 TeV LHC has been studied in [160] 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 [161] (see figure 7 of [161]). The projected upper limit on the branching ratio of H 2 → H 1 H 1 is roughly 5% [161], and has been represented in figure 16 as a dashed red line. 17 The color code of figure 6 has been followed in figure 16.
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.  [95]. The dashed red line corresponds to the future projection for 14 TeV, 300 f b −1 LHC derived in [161] using the results from [160], assuming σ H2 /σ HSM = 1.

JHEP05(2019)177
In [162], 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 [163]. 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 [162,163] are not sensitive for the parameter space of our interest (M H 2 < 10 GeV).

H
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 integrated luminosity of 3000 f b −1 . This mass range of Z D is particularly motivated by the preferred choice of a small v D which is favourable for thermalization (as discussed in section 6.2), which eventually results in a small Z D mass. Below M Z D 2 GeV, prospects of observability in collider searches becomes very weak. 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 [164], with M H 2 fixed at 125 GeV. The subsequent decay process, followed by showering and hadronization is performed by Pythia-6 [165]. The 4l background has been generated by MadGraph5 at the partonic level, while Pythia-6 has been used for showering.

JHEP05(2019)177
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 JHEP05(2019)177  figure 17, 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 [95], 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 figure 18. 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 figure 18. 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 [95]. 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 (M OS = 2.846 − 3.9861 GeV) or Υ resonance (M Υ(2s,3s,4s) = 8.7603 − 11.1052 GeV), are vetoed. The selection cuts are summarized in table. 3. The tt and bbll backgrounds were also evaluated, and no events survived the selection cuts.

JHEP05(2019)177
Fraction of events surviving selection cuts from  4. Upper limits are obtained on the ratio of production cross-section of SM-like H 2 through ggF mode (σ ggH 2 ) and the ggF production cross-sections in SM (σ ggH SM ) 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 figure 19. The branching ratio for Z D → µ + µ − has been taken from [92]. ATLAS has also performed a search for Z D in the H 2 → Z D Z D → 4µ final state using the LHC Run-II data collected at an integrated luminosity of ∼ 36.1 f b −1 [95], and has derived upper bounds on σ ggH 2 /σ ggH SM × Br(H 2 → Z D Z D ) at 95% CL. We have represented these upper bounds as dashed red lines in figure 19. It can be observed that the upper limits JHEP05(2019)177 for HL-LHC presented in figure 19 are roughly ∼ 10 − 15 times stronger than the current bounds. These projected limits for HL-LHC would be capable to probe sin θ mix up to sin 0.002 (0.007) in the M Z D ∼ 2 GeV (M Z D 4 GeV) region.
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 [166][167][168] 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 figure 19 reveal that BP1 could be marginally probed at HL-LHC. 18 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 [169,170]. 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 case, χ − constitutes a source of / E T . 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 section 7.3.2, 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 [164], while Pythia-6 [165] 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 section 7.3.2 is applied here as well, and an event is required to have exactly four isolated muons in the final state. We JHEP05(2019)177 T represents the highest p T muon and p 4 T represents the lowest p T muon. When compared to the case of section 7.3.2, 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 (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 figure 21, following the color code of figure 20. ∆R max (∆R min ) represents the maximum (minimum) value of ∆R between all possible muon pairs. It can be observed from figure 21 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 section 7.3.2. An event is required to have exactly two SFOS pairs, with the invariant mass of the leading and sub-leading SFOS pairs required to be within the range 0. 88      In the recent times, non-traditional search methodologies for beyond the Standard Models (BSM) physics have started garnering popularity. Search for the late decaying long-lived particles (LLP) which feature a characteristic secondary vertex has emerged as one such avenue. In general, a particle is said to be long-lived if its proper decay length exceeds 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 figure 2 of [171]).
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 the radii 554 mm − 1082 mm. The ECAL segment JHEP05(2019)177 Figure 22. 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.
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 figure 22 (following the ATLAS detector geometry), JHEP05(2019)177 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 figure 22 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 figure 22 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 figure 23. Here, we show the total number of gg →  to the similar quantities in figure 22. A higher boost in the case of 27 TeV collision results in an upward shift of color pattern in figure 23 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 boson. 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, 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 restricted 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

JHEP05(2019)177
if there exists a substantial parameter space for a thermal dark matter with large selfinteractions. The latter is desirable since it solves various 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 entire 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 flavor 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 sizeable 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 the 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 sizeable self-interactions). Our findings suggest that such points on the parameter space satisfying both the limits are very rare and generally prefers 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 freezeout [172] and freeze-in [173][174][175][176][177] to achieve the correct relic density, although the latter requires substantially small portal couplings and would be difficult to probe. λ 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 find that it is difficult to probe the MeV scale dark H 1 in our model via future collider experiments. However, the dark gauge boson Z D 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 into 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 , 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.