A Sub-GeV Low Mass Hidden Dark Sector of $SU(2)_H \times U(1)_X$

We present a detailed study of the non-abelian vector dark matter candidate $W^\prime$ with a MeV-GeV low mass range, accompanied by a dark photon $A^\prime$ and a dark $Z^\prime$ of similar masses, in the context of a gauged two-Higgs-doublet model with the hidden gauge group that has the same structure as the Standard Model electroweak gauge group. The stability of dark matter is protected by an accidental discrete $Z_2$ symmetry ($h$-parity) which was usually imposed ad hoc by hand. We examine the model by taking into account various experimental constraints including dark photon searches at NA48, NA64, E141, $\nu$-CAL, BaBar and LHCb experiments, electroweak precision data from LEP, relic density from Planck satellite, direct (indirect) detection of dark matter from CRESST-III, DarkSide-50, XENON1T (Fermi-LAT), and collider physics from the LHC. The theoretical requirements of bounded from below of the scalar potential and tree level perturbative unitarity of the scalar sector are also imposed. The viable parameter space of the model consistent with all the constraints is exhibited. While a dark $Z^\prime$ can be the dominant contribution in the relic density due to resonant annihilation of dark matter, a dark photon is crucial to dark matter direct detection. We also demonstrate that the parameter space can be further probed by various sub-GeV direct dark matter experimental searches at CDEX, NEWS-G and SuperCDMS in the near future.


I. INTRODUCTION
The stability of heavy dark matter (DM) is usually implemented in many particle physics models beyond the standard model (SM) by imposing a discrete Z 2 symmetry in the Lagrangian. These models include the simplest scalar phantom model by adding just a singlet Z 2 -odd scalar field (real or complex) to the SM [1], the popular inert two-Higgs-doublet model (I2HDM) [2,3], the minimal supergravity standard model with R-parity (MSSM) [4][5][6], little Higgs model with T -parity (LHM) [7][8][9] etc. In I2HDM, the DM candidate can be either the CP-even or -odd scalar residing in the second Z 2 -odd Higgs doublet. Many detailed analysis of DM phenomenology in the scalar phantom models and I2HDM can be found in the literature in [10][11][12] and [13][14][15][16][17][18] respectively. In MSSM, the lightest supersymmetric particle (LSP) for the DM candidate can be the spin 0 sneutrino (the superpartner of neutrino) or the lightest spin 1/2 neutralino [19] (in general a linear combination of two gauginos and two Higgsinos [20]). We note that in some low energy supergravity models, the LSP can be the spin 3/2 gravitino, the superpartner of graviton. For a review of supersymmetric dark matter, see for example [21]. In LHM, the spin 1 T -odd partner of the photon can be the DM candidate whose collider implication was studied in [22]. There are also well-motivated non-abelian dark matter models based on additional gauge group like SU (2) [23][24][25][26][27][28], in which the extra spin 1 gauge boson W can be a DM candidate. Moreover, instead of specifying an underlying dark matter model, one can also use the effective dark matter theory approach [29][30][31] to discuss various dark matter phenomenologies [32][33][34][35][36][37].
Recently a gauged two-Higgs-doublet model (G2HDM) based on an extended electroweak gauge group SU (2) L ×U (1) Y ×SU (2) H ×U (1) X , was proposed [38], in which a hidden discrete Z 2 symmetry (h-parity) [39] arises naturally as an accidental symmetry rather than imposed by hand. This discrete symmetry ensures the stability of the DM candidate in G2HDM, which can be either a complex scalar (which in general is a linear combination of various fields in the model) or a heavy neutrino ν H or an extra gauge boson W (p,m) , all of which have odd h-parity. We note that, unlike the left-right symmetric model [40,41], the W (p,m) in G2HDM do not carry electric charge.
The novel idea of G2HDM, as compared with many variants of general 2HDM [42], is that the two Higgs doublets H 1 and H 2 of SU (2) L are grouped into a fundamental representation of a new gauge group SU (2) H . Consistency checks of the model were scrutinized for the scalar and gauge sectors in [43] and [44] respectively. In [39], a detailed phenomenological analysis of the scalar DM candidate in G2HDM was carried out. In general, the scalar DM candidate is a complex field made up of a linear combination of the components from the inert Higgs doublet H 2 and the SU (2) H doublet Φ H and triplet ∆ H . By performing a detailed parameter scan in the model it was demonstrated [39] that only the triplet-like DM is favored when all the constraints from the relic density, direct and indirect searches for the DM are taken into account. As discussed in [38], the triplet ∆ H plays the primary role as a trigger for spontaneous symmetry breaking of the model down to U (1) EM . We note also that this triplet can have topological implications from the hidden sector. Since ∆ H is an adjoint representation of SU (2) H , there exists magnetic monopole [45,46] and dyon [47] solutions in the hidden sector which may play the role of topological stable DM [48]. However, the triplet is not required to generate realistic mass spectra for all the particles in G2HDM.
Therefore if one omits this triplet scalar, the scalar DM candidate is no longer favorable in the parameter space of G2HDM according to the analysis in [39]. However, as mentioned above, there are two other alternative DM candidates in the model. In this paper, we will show that the non-abelian gauge boson W (p,m) associated with SU (2) H can be a viable DM as well. In particular we will focus on the low mass DM scenario in the MeV-GeV range which has attracted a lot of attention in recent years. This paper is organized as follows. In Sec. II, we will review some salient features of the simplified G2HDM without introducing the Higgs triplet field ∆ H of the extra SU (2) H . The theoretical constraints on the Higgs potential, electroweak precision data for the Z-boson mass shift from Large Electron-Positron Collider (LEP) data, dark photon constraints from various low energy experiments and the 125 GeV Higgs data constraints from the Large Hadron Collider (LHC) are discussed in Sec. III. We then turn to the experimental constraints of dark matter physics in Sec. IV. The cosmological relic density from Planck satellite [49], underground direct detection constraints from CRESST III [50], DarkSide-50 [51] and XENON1T [52], astrophysical gamma-ray indirect detection constraints from Fermi-LAT [53,54] and mono-jet constraints from LHC [55][56][57] for the sub-GeV dark matter W (p,m) in G2HDM are studied in Secs. IV A, IV B, IV C and IV D respectively. Our numerical results are presented in Sec. V. We conclude in Sec. VI. In Appendix A, we discuss the mixing effects in the gauge fixings of the model in a general renormalizable R ξ gauge. This work can be regarded as an expanded detailed version of the compact and partial results presented in [58].

II. THE SIMPLIFIED G2HDM MODEL
In this section, we discuss a simplified version of the G2HDM first proposed in [38]. In particular, we will remove the triplet scalar ∆ H in the original model because it is not absolutely required for a realistic particle spectra and the number of free parameters in the scalar potential can be reduced significantly. We note that the Yukawa couplings are not affected by this simplification since the triplet does not couple to the fermions in the model.

A. Particle Content
The gauge group of the simplified G2HDM is the same as in [38], In Table I, we summarize the matter content and their quantum number assignments in G2HDM. At the minimum risk of confusion, we will continue refer this model as G2HDM to avoid cluttering throughout the paper with the adjective word "simplified". What we are really concerned about is the electroweak part of G, so the color group SU (3) C is not relevant in what follows.

B. Higgs Potential and Spontaneous Symmetry Breaking
The most general Higgs potential which is invariant under where (α, β, γ, δ) and (i, j) refer to the SU (2) H and SU (2) L indices respectively, all of which run from 1 to 2, and H αi = H * αi . To facilitate spontaneous symmetry breaking (SSB) and obtain the particle mass spectra of the model, we shift the fields based on our conventional wisdom where v and v Φ are the vacuum expectation values (VEV) of H 1 and Φ H fields respectively.
H 2 is the inert doublet in G2HDM and hence does not have VEV. Naively we would think that the Goldstone bosons G + , G p H , G 0 and G 0 H will be absorbed by the longitudinal components of W + , W p , W 3 and W 3 respectively. In Appendix A we will show that the last three Goldstone fields have mixing effects with other fields in the scalar sector that makes the situation more interesting but a little bit more complicated.
Substituting the scalar field decomposition of Eq. (2) into the scalar potential Eq. (1) and then minimize the potential, one can obtain the solutions of VEVs as follows Equivalently, we can use the minimization conditions to trade µ 2 H and µ 2 Φ with v and v Φ as

C. Scalar Mass Spectrum
In the S = {h, φ 2 } basis the mass matrix is given as One can use an orthogonal transformation O S , which can be parametrized as where with h 1 being identified as the 125 GeV SM-like Higgs boson and h 2 as a heavier scalar boson. The mass squared eigenvalues of Eq. (7) are given by In the basis of S = {G p H , H 0 * 2 }, we obtain the mass matrix: Similarly, this mass matrix can be diagonalized by an orthogonal matrix, where which gives We note that, in Eq. (15), the zero eigenvalue corresponds to the Nambu-Goldstone boson mass eigenstateG p H , while the other eigenvalue is the mass of a new dark scalar boson.
The charged Higgs boson mass is given as The Goldstone bosons G 0 and G 0 H are massless. The above scalar mass spectrum is derived in the so-called 't Hooft-Landau gauge. We will discuss further the mixing effects of the Goldstone bosons with other scalar fields in a general renormalizable R ξ gauge in Appendix A.
We note that h 1,2 , G 0 and G 0 H are even under h-parity, whileG p H , D and H ± are odd [39].

D. Gauge Sector
After SSB, the W ± gauge boson of SU (2) L remains the same as in SM with its mass given by m W = gv/2. The SU (2) H gauge boson W (p,m) receives mass from H 1 and Φ 2 given by Note that W (p,m) are electrically neutral and thus do not mix with the SM W ± . In addition, obtains immediately a zero eigenvalue identified as the SM photon and a 3 × 3 sub-matrix in the basis of Z SM , W 3 , X given by [44], where g, g , g H and g X are the gauge couplings of SU (2) L , U (1) Y , SU (2) H and U (1) X respectively, and m Z SM = v g 2 + g 2 /2 is the SM Z boson mass expression. The mass matrix in Eq. (19) can be diagonalized by an orthogonal rotation matrix O so that 1 In this analysis, we arrange the neutral gauge boson masses as m A < m Z < m Z m Z SM with Z identified as the physical Z boson with mass 91.1876 ± 0.0021 GeV [59]-the first heavy neutral vector gauge boson discovered in 1983 at the Super Proton Synchrotron at CERN. Two more massive neutral vector gauge bosons are predicted in G2HDM.
By means of a few assumptions motivated by our expectations, it is possible to draw a few conclusions from Eq. (19). First, the new gauge couplings g H and g X are expected to be much smaller than the SM g and g to avoid large effects on the very precise measurements of the Z properties. Second, the scale of v Φ is expected to be larger than v given that it characterizes the scale of new physics and is directly related to the masses of beyond the SM (BSM) states. By neglecting any term composed by a product of any three or more of g H , g X and v 2 /v 2 Φ it is possible to put Eq. (19) into a block diagonal matrix where only W 3 and X mix resulting in the approximation Moreover, the 2×2 squared mass matrix of W 3 and X can be easily diagonalized and somewhat simple approximations can be found. Since we want the hierarchy m A < m Z < In our study, the couplings of the extra gauge bosons to the SM charged leptons and to quarks u and d will be important for dark photon constraints and direct detection of DM.
The vectorial and axial parts of their couplings are given by where i = 1, 2, 3 and Z(1) ≡ Z, Z(2) ≡ Z , Z(3) ≡ A . Using the relations between mixing matrix elements and mixing angles in Eqs. (2.6) to (2.9) of Ref. [44], we can find the following relation 2 which can be used to simplify the above vectorial and axial couplings to obtain where r i = m Z(i) /m Z SM . In this form, it is obvious that the axial couplings magnitude is expected to be smaller than the ratio of squared masses r 2 i , which, e.g., for m Z(i) = 1 GeV is already close to 10 −4 . In contrast, the vectorial couplings have an r i -independent part whose size is controlled by the mixing matrix element O 1i that also affects axial couplings.
Therefore, for sufficiently light Z and A , the contributions from the axial couplings are expected to be subleading.
The neutral current interactions induced by W (p,m) are given by where with V H u , V H d , V H ν and V H e being the new unitary mass rotation matrices for the fermions. The neutral current interactions induced by Z, Z , A bosons can be found in Ref. [44].

E. Free Parameters
From the scalar potential we recognize 9 free parameters including couplings and VEVs.
Of those 9 parameters, µ 2 H and µ 2 Φ can be related to other parameters using the minimization conditions according to Eqs. (5) and (6), leaving only 7 free parameters. Furthermore, Eq. (10) can be used to relate the parameters λ H , λ Φ and λ HΦ to the physical squared masses m 2 h 1 , m 2 h 2 and the mixing angle θ 1 obtaining the following relations Using the mass of D, Eq. (16), we can express λ HΦ in the form Using this last expression and the mass of the charged Higgs of Eq. (17), we can write λ H as Finally, we can use Eq.
Using the expressions in this subsection allow us to trade six model parameters with five physical squared masses and one mixing angle, The remaining free parameters of the model are the heavy fermion masses m f H , the Stueckelberg mass M X , and the gauge couplings g X and g H . Considering that the mass of the Higgs, m h 1 , has already been measured [59], we are left with a total of 8 free parameters plus the masses for 12 heavy hidden fermions. The effects of heavy hidden fermions in complex scalar dark matter phenomenology in G2HDM was analyzed in [60].

III. CONSTRAINTS
In this section, we examine the model using various constraints including the theoretical constraints of the scalar potential, electroweak precision data, dark photon physics and Higgs measurements at LHC. DM constraints are presented separately in the next section.

A. Theoretical Constraints
The theoretical constraints on the original G2HDM were studied in Ref. [43]. Here, we follow closely on the steps of that work but remove the scalar triplet and related parameters from the original model.
Vacuum Stability: To make sure the scalar potential is bounded from below, the sum of all quartic terms in the scalar potential needs to be positive. In the same way as Ref. [43] we use copositivity conditions given by the following constraints where λ H (η) ≡ λ H + ηλ H and λ HΦ (ξ) ≡ λ HΦ + ξλ HΦ . The conditions of Eq. (42) have to be met for any value of ξ and η in the ranges 0 ≤ ξ ≤ 1 and −1 ≤ η ≤ 0.
Perturbative Unitarity: It is required that the parameter space remains within the perturbative limits. To this end, here we compute the 2→2 scalar scattering amplitudes induced by the quartic couplings. The 2→2 processes induced by vertices from scalar cubic couplings and gauge interactions are suppressed by large momentum exchange in their propagators [43].
Perturbative unitarity requires

B. Electroweak Constraints
A comprehensive study on the electroweak precision constraints in the original model has been preformed in Ref. [44]. We found that the Z mass shift is the most stringent among all the electroweak precision constraints. In particular, for the parameter space of interest in our analysis, i.e., g H ∼ g X 1 and (m 2 where g Z = g/ cos θ W with θ W being the Weinberg angle. Following the methodology of Ref. [61], we can obtain the experimental uncertainty of the Z mass as where t W = tan θ W and ∆r is the radiative correction. Using the PDG values [59] of m W ± δm W = 80.387 ± 0.016 GeV, ∆r ± δ∆r = 0.03652 ∓ 0.00021 ± 0.00007 and sin θ W = 0.22343, and by requiring |∆m Z | < |δm Z |, one obtains an upper bound on g H and g X

C. Dark Photon
The light boson A can be treated as a dark photon and, therefore, dark photon constraints have to be applied. In particular, due to the vertex A ¯ , where represents any charged lepton, it is expected that for a sufficiently large coupling it should be possible to observe the A resonance in the invariant mass distribution of e + e − and µ + µ − . Dark photon experiments constrain the size of the coupling via a parameter ε . In the decay width Γ (A → ) the parameter ε appears as [62] Γ where µ = 2m /m A < 1 since this decay channel only opens for m A > 2m . In the G2HDM, the parameter ε at tree level is given by where v A and a A are given in Sec. 2 of Ref. [44]. From Ref. [44] we know that v A and a A are the same for all the charged leptons and, thus, the only distinction in ε between different flavors of leptons comes from µ . As mentioned in Sec. II D, for light enough A the axial coupling will be negligible and ε is expected to be nearly independent of µ , as is usually the case in models with dark photon. It is important to mention that, since Z is also expected to be light, the dark photon experimental limits can also be applied as above with A → Z in Eq. (50) and {v A , a A } → {v Z , a Z } in Eq. (51). However, since A is lighter by definition it is expected to be more strongly constrained.
The existing limits on ε for a dark photon mass m A > 1 MeV are displayed on the top pane of Fig. 10 in Ref. [62].
D. Higgs Collider Data

Higgs Boson Mass
As aforementioned, h 1 is identified as the observed Higgs boson at the LHC. In this analysis, we take the mass of Higgs boson as m h 1 = 125.10 ± 0.14 GeV [59].

Higgs Decays into Diphoton
The decay rate for h 1 → γγ is given by Here, and in what follows, we define τ i = 4 m 2 i /m 2 h 1 where i indicates which particle is running inside the loop. N c is the color factor, 1 for leptons and 3 for quarks. The coupling g h 1 H + H − corresponds to the h 1 H + H − vertex and is given by The well-known loop function f (x) is The signal strength parameter for the Higgs boson produced from the gluon-gluon fusion (ggH) can be obtained as where the superscript SM refers to the SM Higgs boson h. The decay width of h 1 into two gluons is given by [38] where with q SM and q H refer to the SM quarks and the new colored fermions. The latest measurement of this signal strength is given by ATLAS as 0.96±0.14 [71].
The decay width of Higgs boson to SM fermions is given by The signal strength for ggH production is then given by Note that this expression is independent of the fermion flavor and thus we compare against the best measured signal strength given by the decay into a pair of τ + τ − : µ τ τ ggH = 1.05 +0. 53 −0.47 [72].

Invisible Higgs Decay
If m h 1 > 2 m W , the Higgs boson can decay invisibly into a pair of W (p,m) . The decay In our parameter choice, we will assume 2 m ν H > m h 1 so that h 1 does not decay into a pair of ν H . The branching ratio of invisible Higgs decay is then given by Recently, the ATLAS collaboration reported the most stringent constraint on the invisible decays of the Higgs produced via vector boson fusion. Assuming that the Higgs boson production cross section is comparable to the SM, the ATLAS collaboration set the limit BR(h 1 → inv) < 0.13 at 95% C.L. [73].

A. Relic Density
The DM scenario presented here works similarly to the very well known WIMP DM.
The DM candidate W (p,m) begins in thermal equilibrium with the particle species in the hot primordial soup in the early universe before starting to freeze-out due to the expansion of the for the main annihilation processes as given in Fig. 1 can be computed straightforwardly and give rise to the following total cross section for each final fermion pair where The annihilation mediated by the Z is suppressed by a factor of O 2 21 required to be small mostly by measurements on the decay width of the Z and the decay branching fractions that limit the Z → W p W m process. As mentioned above, the channel mediated by the Z is also suppressed by combinations of gauge couplings that make most of the size of v Z f and a Z f . However, these suppressions are not as strong as the suppression in other channels and when we include the effects from Z resonance it is possible to bring the relic density to its expected value of Ωh 2 ∼ 0.1.
In our study, we will consider the measured value of Ωh 2 = 0.120 ± 0.001 as given by the Planck collaboration [49].

B. Direct Detection
Due to the small coupling between the DM candidate, W (p,m) , with the SM-like states where h j mediated processes are suppressed by a factor of g 2 H m 2 q /v 2 in the cross section with an additional m −4 h j suppression from the propagator since these interactions happen via t-channel. Therefore, we are only left with the processes mediated by Z, Z and A in the t-channel. Usually, for direct detection processes the momentum exchange is considered to be very small and therefore t is expected to be small as well. This will result in amplitudes suppressed by the inverse squared of the mass of the mediator meaning that the light states, Z and A , will be less suppressed. In the approximation where the momentum exchange is smaller than the mass of the mediator, we can write the interaction between DM and light quark q as a contact interaction given by Furthermore, for the axial part of the interaction with the quark, in the small momentum exchange limit, only the space components of γ ν remain but these components are suppressed by the W (p,m) momentum due to the derivatives ∂ ν W (p,m) in Eq. (70) [74,75]. This, together with axial couplings that are comparably much smaller than the vectorial ones results in an spin dependent cross section that is expected to be several orders of magnitude smaller than the spin independent one.
From Eq. (70), it is clear that the A mediated process is expected to dominate the cross The case where both mediators participate equally is expected to happen only through fine tuning of masses and mixings. Therefore, we expect the cross section with the nucleons to be mostly mediated by either A or Z . The elastic cross section between W (p,m) and a nucleon, N , is given by where i = 2 or 3 depending on the dominant mediator according to the discussion above, is the reduced DM-isotope nucleus mass and f p and f n are effective couplings of the DM with protons and neutrons, respectively. The atomic number is Z atom and the isotope dependent variables η k and A k are the abundance and mass number of the k th target isotope, respectively. Direct detection experiments usually report the number in Eq. (71) assuming isospin conservation, i.e., f p = f n . In that case, it is straightforward to see that the ratio of the sums over isotopes reduces to 1 and σ SI W N = σ SI W p . However, in our case the couplings between quarks, u and d, and the gauge bosons, Z and A , are all different due to their distinct SM charges leading to isospin violation (ISV), i.e., f p = f n . Following Refs. [76,77], we can rescale the reported experimental limit, σ limit → σ limit × σ SI W p /σ SI W N to account for ISV effects and use it to limit σ SI W p as given by Eq. (72). This rescaling depends on the mass of DM, the atomic numbers and the ratio f n /f p , and, therefore, will be different for different points in the parameter space.
To constraint the W (p,m) -proton cross section we will use the most recent upper limits set by the experiments CRESST III [50], DarkSide-50 [51] and XENON1T [52].

C. Indirect Detection
Due to DM annihilation before freeze out happening through the resonance of an otherwise suppressed channel, the annihilation of DM in the present-after the shift in energy from the early to the current Universe-loses the resonance resulting in a very low annihilation cross section. We have checked that the value of the total annihilation cross section in G2HDM at the present time is of order 10 −32 cm 3 ·s −1 or below, much lower than the canonical limits set for various channels by Fermi-LAT data [53,54]. The occurrence of energetic jets with large missing transverse momentum has been searched by ATLAS [55,56] and CMS [57] collaborations. However, the observed results are overall in agreement with the SM predictions and only exclusion limits have been reported.
In our model, the process pp → W p W m j can give rise to mono-jet events at the LHC.
To analyze the mono-jet signal at the LHC, we choose two benchmark points (BPs) which are shown in Table II. These two BPs satisfy all theoretical, Higgs data and DM constraints.
We evaluate the signal process cross section using MadGraph 5 [78] with precuts for jets p j T > 30 GeV and |η j | < 2.8, and for the missing transverse momentum p miss T > 100 GeV. It turns out the production cross sections with the precuts are about 3.0 fb and 3.8 fb for BP 1 and 2 respectively, and dominated by Z and Z mediated diagrams. We generate 10 4 events for the pp → W p W m j process and recast ATLAS mono-jet search [56] using MadAnalysis 5 [79].
The most sensitive signal region is found to be in the window E miss T ∈ (700, 800) GeV (the signal region EM7 in Ref. [56]). The 95% C.L. exclusion limits on the production cross section are 400 fb and 680 fb for BP 1 and 2 respectively, which are much larger than the signal expected from the model. Therefore, the LHC with luminosity of 139 fb −1 is not sensitive enough to search for mono-jet events from this model. However, the model can be probed by mono-jet searches at future hadron colliders such as the High-Luminosity Large Hadron Collider (HL-LHC) [80], the High-Energy Large Hadron Collider (HE-LHC) [81] and the Future Circular hadron-hadron Collider (FCC-hh) [82].

A. Methodology
The masses of the gauge bosons, m Z , m Z and m A , as well as the constraints from Sec. III are calculated through our own fortran codes, except for the Higgs invisible decay which can be calculated together with the DM constraints. The DM constraints of Sec. IV, in particular relic density, direct detection and indirect detection are calculated using micrOMEGAs [83] and a set of model files generated by FeynRules [84]. For the invisible decay branching ratio of the Higgs, we take advantage of the use of CalcHEP [85] within micrOMEGAs to calculate the decay width along with the rest of the DM constraints just mentioned.
All the points outside the theoretical constraints of Sec. III A are simply rejected. By the same token, the dark photon constraints are used to reject any parameter combination of ε e,µ -m A or ε e,µ -m Z located inside the currently excluded regions. The rest of the constraints in Sec. III are summed into a total χ 2 that also includes relic density and direct detection cross section. In the case of direct detection experiments, where a limit is reported at a 95% C.L. with null-signal assumption, we use a χ 2 DD of the form where the 4.61 factors allows χ 2 DD = 4.61 when we are exactly at the 95% C.L. of this two-dimensional limit. 3 In mass ranges where more than one limit exists we take the one with the largest χ 2 DD . Note that, due to ISV, the largest χ 2 for direct detection may not correspond to the experiment with the smallest cross section. Since direct detection limits are reported assuming f p = f n in Eq. (71), it is possible for ISV (f p = f n ) to produce some amount of cancellation or enhancement of the limits depending on the atoms used in the detector. Calculating the cross section in the way described in Sec. IV B allows us to account for ISV and the atoms used in different experiments.
In the case of Higgs invisible decay branching fraction, where a limit is reported with a 95% C.L., the appropriate χ 2 inv is given by where, similarly to direct detection, the 2.71 factor allows for χ 2 inv = 2.71 when our result is exactly at the reported 95% C.L. in the one-dimensional case.

Parameter [units]
Range  To sample the parameter space we use the affine invariant Markov Chain Monte Carlo (MCMC) ensemble sampler emcee [86] which presents advantages such as fast calculation of parameter distributions in several dimensions. The initial prior and ranges of each parameter are contained in Table III. In particular, the parameters m W , M X , g H and g X are scanned in base-10 logarithmic scale. This is mostly because we expect these parameters to be small but different from zero and that their effects depend heavily on their orders of magnitude.
For the masses of the heavy fermions, we expect their contributions to be heavily suppressed by the requirement that m f H = O(1 TeV). Therefore we consider all of them degenerated with a mass of m f H = 3 TeV putting them safely above any current search for heavy fermionic states. The rest of the parameters are scanned with a uniform prior in linear scale.
To guarantee that our final distributions are independent of the initial points we perform several small runs collecting O(10 4 ) points for each run using O(100) walkers. The initial points for the walkers are always allowed by theoretical and dark photon constraints but otherwise random inside the prior. After checking that the final distributions are consistent between different runs, we perform a large scan with 300 walkers collecting 160,000 points after burn-in and thinning.

B. Numerical Results
We present the numerical results for visualization in Figs. 2-5. To follow the discussion below more smoothly, we suggest our readers to read the captions of these figures first and then view and compare them in parallel.
The most notable feature of both panes of Fig. 2 is the band-shaped allowed region. In the case of (m W , g H ) plane shown in the left pane, the band is caused by the relation between relic density and cross section, Ωh 2 ∝ 1/ σv . Considering that we have σ ∝ g 2 H m 2 W /s 2 from Eq. (67), assuming s ∼ 4 m 2 W we have that g 2 H m 2 W /s 2 ∼ g 2 H /(16 m 2 W ) resulting in Ωh 2 ∝ m 2 W /g 2 H . This means that to keep a constant relic density, m W and g H have to keep a linear relationship as displayed in the left pane of Fig. 2. Deviations from this band result in the relic density going either above or below the value measured by the Planck satellite.
In the case of the right pane of Fig. 2, the band can be explained by the possibility of having a resonant annihilation of W (p,m) mediated by the Z . First note that Eq. (23) directly relates m A and M X and, as can be seen in the right pane of Fig. 3, m A is required to be mostly below 0.2 GeV due to the LHCb results, thus limiting also the size of M X .
Then, from Eq. (22) we know that the term with m 2 W factor dominates over the term with M 2 X . Finally, resonant annihilation is achieved for m Z ≈ 2 m W meaning 1 + 4g 2 X /g 2 H ≈ 4 or g 2 X /g 2 H ≈ 3/4, resulting in the band seen in the right pane of Fig. 2. Again, large deviations from this band result in too much or not enough annihilation to achieve the correct relic density. Due to this g X -g H correlation, the two-dimensional allowed regions in Fig. 4 for all the distributions involving g X and g H have similar shapes. Here it is important to mention that exact resonance, m Z = 2 m W , would result in too much annihilation and, therefore, in Ωh 2 well below Planck's measurement.
Besides their band-shaped tendency, both panes in Fig. 2 are bounded in their topright and bottom-left corners by the DM direct detection and dark photon constraints, respectively. We know that the direct detection cross section grows with g 2 H as seen in Eq. (72), therefore, it is expected to see it setting an upper bound on g H . Furthermore, as can be seen in the left pane of Fig. 3, direct detection experiments practically create a wall that limits the size of m W from above. The effects of this limit are reflected in the upper bound of m W in the left pane of Fig. 2. In the case of the region disfavored by dark photon searches, this is mostly due to the ν-CAL I experiment limiting ε from below as seen in the right pane (olive green shaded zone) of Fig. 3. The ε coupling limit is passed to g H through the vectorial and axial couplings v A f and a A f that depend on it, resulting on the lower limit on g H that can be seen in both panes of Fig. 2. In the left pane of Fig. 3 we show the allowed region projected on the (m W , σ SI W p ) plane, where σ SI W p is the cross section for spin-independent scattering on a proton. The dark (light) blue shaded zone represents the 1σ (2σ) allowed region. The current DM direct detection measurements from CRESST III (green) [50], DarkSide-50 (orange) [51] and XENON1T (brown) [52] constrain the DM mass to remain below ∼ 2 GeV. A small part of the 2σ allowed region lies below the neutrino floor (light orange), where the coherent neutrino-nucleus scattering would dominate over any DM signal. Additionally, we show that experiments in the near future such as NEWS-G [87], SuperCDMS [88] and CDEX [89] can further probe our allowed parameter space, in particular for m W > ∼ 0.3 GeV with NEWS-G and down to σ SI W p ∼ 10 −44 cm 2 with SuperCDMS and CDEX. The right pane in Fig. 3 shows the allowed region projected on the (m A , ε) plane with the coupling ε ≡ ε . Various experimental limits from dark photon searches are displayed in color shaded zones including LHCb (green) [63], BaBar (pink) [64], NA48 (purple) [65], NA64 (light brown) [66], E141 (magenta) [67] and ν-CAL I (light green) [68,69]. The dilepton searches at the LHCb, BaBar and NA48 put upper limits of ε 10 −3 for m A > ∼ 0.03 GeV, especially LHCb which sets a strong limit on ε at 0.2 GeV < m A < 0.5 GeV causing a concave region in the 2σ allowed region at this mass range. We note that this concave region due to LHCb corresponds to the concave region at (m W , σ SI W p ) ∼ (1 GeV, 10 −42 cm 2 ) in the left pane of the same figure. The LHCb long lived dark photon search constraints [63] are also shown by the two isolated green shaded islands around equals 2 × 10 −5 . On the other hand, the beam dump experiments NA64, E141 and ν-CAL I close the available space for smaller ε and lighter m A setting lower bounds of m A > 0.02 GeV and ε > ∼ 2 × 10 −5 . The lower limit on ε for m A > 0.05 GeV is due to the DM relic density measured by the We present two-dimensional projections of the allowed region for our BSM parameters in Fig. 4. The dark, light and lighter blue zones indicate the 1σ, 2σ and 3σ allowed regions respectively. One important thing to note is that the m h 2 , m D and m H ± masses show an apparent upper bound at 1σ and 2σ. This upper bound actually depends on the maximum value chosen for the prior of these three parameters and has no physical meaning. In the case of m h 2 , the apparent limit is due to the reduced θ 1 for large m h 2 seen in the (m h 2 , θ 1 ) subfigure while for m D and m H ± it is due to their near degeneracy shown in the (m D , m H ± ) subfigure.
We have checked that changing the maximum scanned value for these three parameters does not change the rest of the distributions, except for θ 1 where, understandably, larger m h 2 sharpens the peak at θ 1 = 0 where h 1 = h SM exactly. The near degeneracy for m D and m H ± is due to their comparably small mass squared difference m 2 D − m 2 H ± = v 2 (λ HΦ + λ H )/2. Given that (λ HΦ + λ H )/2 is bounded by unitarity constraints, the mass squared difference is expected to remain O(v 2 ) or less, meaning that as m D and m H ± grow away from v their proportional difference rapidly grows smaller. Their near degeneracy is also noticeable in all their two-dimensional distributions since the distributions become nearly identical as m D and m H ± grow larger.
Another interesting feature is that the charged Higgs mass distributions reveal the presence of a lower limit around m H ± ≈ 400 GeV where the contours show that the distribution falls rapidly. This is due to the constraint from the Higgs decays into diphoton as shown in Eq. (52). Moreover, due to the relation between M X and m A that can be inferred from Eq. (23), the distributions in the (M X , g X ) and (M X , g H ) subfigures (second column from the right in (Fig. 4)) are close to the distribution shown for the dark photon constraint in the right pane of Fig. 3.
We present the marginalized one-dimensional distributions for the most relevant parameters in Fig After the above long discussions of our numerical analysis, perhaps a high level summary is useful.
• We have focused our numerical analysis on the parameter space of the model that can lead to a sub-GeV W (p,m) DM with a mass range of MeV−GeV.
• The viable domains of the parameter space are summarized in Fig. 4, where the 1-3σ contours allowed by all existing experimental constraints are shown for any combination of two of the 8 free parameters while the rest of the parameters are marginalized.
All the masses of the new heavy fermions in the model, required by anomaly cancellation, are set at 3 TeV. For the two cases of (m W , g H ) and (g H , g X ), more detailed information on the boundaries of the contours due to the different constraints imposed are exhibited in Fig. 2.
• It is both clear and exciting to see from the allowed 1σ and 2σ contours on the plane of m W versus σ SI W p displayed in the left pane of Fig. 3 that future experiments like NEWS-G and SuperCDMS can put more stringent constraints on the model. Indeed, about 1/3 of the current viable parameter space in the (m W , σ SI W p ) plane would be facing challenge. On the other hand, the region where the projected CDEX sensitivity can reach is already disfavored in G2HDM. Note that a small portion of the 2σ allowed parameter space in the (m W , σ SI W p ) plane is overlapping with the neutrino floor.
• For the allowed contours of dark photon coupling ε and mass m A exhibited in the right pane of Fig. 3, future upgraded NA64 experiment with 5 × 10 12 electrons-on-target, proof-of-principle experiment AWAKE run 2 with 10 16 electrons-on-target, and next generation B-factory experiment Belle II would probe more than 1/2 of the current viable domain in the (m A , ε) plane.
• The correlation between the DM and dark photon physics and their constraints in

VI. CONCLUSION
In this work, we studied a simplified version of the G2HDM [38]. In previous works, it was demonstrated that the original G2HDM successfully explains dark matter while keeping a parameter space consistent with theoretical and experimental expectations [39,43,44].
The simplifications considered in this study reduced the size of the scalar potential and the parameter space by removing the scalar SU (2)  We started with the usual theoretical checks on the scalar potential ensuring that the minimum is stable and that the couplings remain unitary at tree level. Given that the LHC has been closing in on the detailed properties of the 125 GeV Higgs, we check that our scalar sector provides a particle, h 1 , that matches the mass and decays that have been measured.
The same can be said in the case of the gauge sector and the Z with its properties already very well measured at LEP. In the case of the light Z and A we constrained their interactions with SM leptons by checking against the regions excluded by dark photon searches in LHCb, BaBar, NA48, NA64, E141 and ν-CAL I. Finally, we required our DM candidate, W (p,m) , to have the correct relic density measured by Planck satellite with a spin-independent direct detection cross section below the limits found by CRESST-III, DarkSide-50 and XENON1T.
In our numerical analysis, we uncovered some interesting features of the parameter space, such as a correlation between the couplings g H and g X with superweak size O(10 −5 − 10 −3 ) and that most of our results lie inside the gap between various dark photon explorations.
The latter becomes more important when we consider that this gap is projected to be further explored in the future with the upgrades to NA64 and AWAKE and the B-factory Belle II.
This would reduce the allowed parameter space nearly by half, and therefore remove a good portion of our allowed parameter space in the lighter A region. Moreover, future DM direct detection experiments like NEWS-G and SuperCDMS may reduce the parameter space by exploring the regions with heavier W (p,m) and larger cross section.
To summarize, we found that the simplified G2HDM developed in this work provides a viable vector DM candidate with mass down to O(10 −2 ) GeV. All the predictions in the scalar and gauge sectors are in good agreement with current observations. Importantly, both new vector states, Z and A , play key roles for DM observables. Besides the possibility of detecting the W (p,m) in DM direct detection experiments, the dark photon, A , is predicted to be well positioned for future observations that may reach m A ∼ 0.1 GeV. This work demonstrates that the G2HDM is not only a successful and competitive dark matter formulation but can also serve as a starting point with diverse exploration possibilities.
Here S and M X are the Stueckelberg field and mass respectively. The middle term gives X a mass M X , while the last term indicates the mixing of S with the longitudinal component of X µ . Note that S is massless.
Recall that the first two terms in Eq. (A8) are SM-like and they can be combined as Thus the total mixing terms including both Eq. (A8) and the last term in Eq. (A9) are Note that, as one would expect, the photon field A µ doesn't enter in Eq. (A11) which is coming entirely from SSB and Stueckelberg mechanism. The photon field remains massless, it has no associated Goldstone boson and its general gauge fixing term is simply given by − (∂ µ A µ ) 2 /2ξ γ as in the SM case. So does the massless gluon field, whose gauge fixing term is − (∂ µ A aµ ) 2 /2ξ g where a = 1, . . . , 8 is the adjoint index of the color group SU (3) C .
As mentioned in the text, the neutral vector gauge bosons Z SM , W 3 and X are in general mixed together according to  where O is an orthogonal matrix. In our numerical scan of the parameter space in this work, we have Z 1 = Z Z SM which is very close to the SM Z-boson, and Z 2 = Z and Z 3 = A , both of which are lighter than the Z. We will use Z i in this appendix instead of Z, Z and A .
In terms of the physical fields Z i , Eq. (A11) becomes Here we have defined three physical Goldstone fields absorbed by the longitudinal components of the three physical Z i as where the coefficients C ij can be read off from the first two lines in Eq. (A13), namely To cancel the mixing term in the last expression in Eq. (A13), we need the following gauge fixing term, the procedure of Faddeev-Popov quantization can be proceeded as usual. Of course physical observables if computed correctly should be independent of all the otherwise arbitrary gauge fixing parameters ξs discussed in this Appendix! In an ideal world, one would compute things using arbitrary ξs and show the dependence of ξs is completely dropped out at the end for any physical observable. In practice, things are hardly get done that way.