Hidden Higgs portal vector dark matter for the Galactic center gamma-ray excess from the two-step cascade annihilation, and muon g − 2

We have built a lepton-specific next-to-minimal two-Higgs-doublet-portal vector dark matter model. The vector dark matter in the hidden sector does not directly couple to the visible sector, but instead annihilates into the hidden Higgs bosons which decay through a small coupling into the CP-odd Higgs bosons. In this model, the Galactic center gamma-ray excess is mainly due to the 2-step cascade annihilation with τ’s in the final state. The obtained mass of the CP-odd Higgs A in the Galactic center excess fit can explain the muon g − 2 anomaly at the 2σ level without violating the stringent constraints from the lepton universality and τ decays. We show three different freeze-out types of the dark matter relic, called (i) the conventional WIMP dark matter, (ii) the unconventional WIMP dark matter and (iii) the cannibally co-decaying dark matter, depending on the magnitudes of the mixing angles between the hidden Higgs and visible two-Higgs doublets. The dark matter in the hidden sector is secluded from detections in the direct searches or colliders, while the dark matter annihilation signals are not suppressed in a general hidden sector dark matter model. We discuss the constraints from observations of the dwarf spheroidal galaxies and the Fermi-LAT projected sensitivity.


JHEP08(2018)099
The evidence for the dark matter (DM) in the Universe has been well established from various astronomical observations and cosmological measurements. The DM, which cannot be accounted for within the standard model (SM) scenario, indicates the existence of new physics. The attractive candidates for the DM are the so-called weakly interacting massive particles (WIMPs) which, having a weak scale mass and annihilating into SM particles via weak scale couplings, can provide a correct thermal relic abundance, following the Boltzmann suppression before freeze-out. Several collaborations have reported an excess of GeV gamma-rays near the region of the Galactic center (GC) [1][2][3][4][5][6][7][8][9][10][11], where the excess spectrum can be fitted using DM annihilation models. Although the excess result might be explained by the millisecond pulsars or some other astrophysical sources [12][13][14][15][16][17][18], the DM annihilation is a viable scenario from the particle physics point of view [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. However, the DM models capable of explaining the GC gamma-ray excess are increasingly constrained by direct detection experiments and measurements at colliders. Some ideas have been proposed that can avoid overproduced signals in the latter two experiments; for instance, the DM annihilates into bb through a pseudoscalar mediator exchange [19][20][21]35]. The present work is motivated by the idea called "secluded WIMP dark matter" in which the DM first annihilates to a pair of short-lived hidden mediators which subsequently decay into SM particles through very small couplings [29][30][31][32][33][34][35][36][37][38][39][40], so that it can easily evade the stringent constraints from current measurements, but provide observable gamma-ray signals in the indirect measurements.
In this paper, we will use the GC gamma-ray excess spectrum obtained by Calore, Cholis and Weniger (CCW) who analyzed the Fermi data in a consistent treatment of systematic uncertainties that came from the modeling of the Galactic diffuse emission background [8]. It is interesting to note that if the dark matter annihilates directly only into the τ + τ − , the spectral fit to the GC gamma-ray GeV excess gives the best-fit result with dark matter mass ∼ 9.5 GeV, and low-velocity annihilation cross section σv 0.37 × 10 −26 cm 3 s −1 , nevertheless corresponding to a lower p-value ∼ 0.05 for the goodnessof-fit test [8]. The direct DM annihilation to τ + τ − produces a gamma-ray spectrum which peaks sharply at a little higher energies. If the annihilation processes present some extra intermediate steps, the final state τ 's generated from the cascade decays are boosted, and therefore the resultant gamma-ray spectrum becomes broader and has a better fit to the GC excess observation [41,42]. We are interested in the two-step cascade annihilation process (see figure 1 for reference). The reason is that not only a much larger p-value ∼ 0.22 can be obtained, but also the fitted DM mass and annihilation cross section are enlarged by a factor of ∼ 4, compared with the direct annihilation to τ + τ − . The resulting annihilation JHEP08(2018)099 JHEP08(2018)099 sectors. In section 4, we first outline the approach of determining the gamma-ray spectrum from a 2-step cascade annihilation to the final state τ 's, and then describe the analysis and results concerning gamma-ray observations, compared with the relic abundance in the conventional WIMP dark matter scenario. Section 5 contains the analysis for the mixing angles between the hidden scalar and two neutron CP-even bosons in the two Higgs doublets. The discussions and conclusions are given in section 6.
2 Lepton-specific next-to-minimal 2HDM portal vector dark matter

The model
We consider a model with two Higgs doublets, Φ 1 , Φ 2 , and a complex scalar dark Higgs field Φ S . The CP-conserving potential for the Higgs sector is described by where Φ S is a singlet under the SM gauge fields. As usual, we have imposed a discrete Z 2 symmetry to the Higgs potential, such that Φ 1 → Φ 1 , Φ 2 → −Φ 2 , and Φ S → Φ S , under which the tree-level flavor changing neutral currents (FCNCs) are absent. The Z 2 symmetry is softly broken by the term containing m 2 12 . On the other hand, we have considered that Φ S is charged in the dark U dm (1) gauge group, while other Higgs fields and SM particles have no such quantum number. The U dm (1) group contains an abelian gauge boson, X µ . After spontaneous symmetry breaking, the vacuum expectation value (VEV) of Φ S generates a mass for X µ , and a discrete Z 2 symmetry: X µ → −X µ , Φ S → Φ * S , is still maintained, such that X µ is stable and can serve as a (vector) dark matter candidate.
The relevant kinetic terms in the dark sector are given by where X µν = ∂ µ X ν − ∂ ν X µ , and the covariant derivative is defined as with Q Φ S the U dm (1) charge of Φ S . After spontaneous symmetry breaking, we have where the imaginary part of Φ S is absorbed by the vector gauge boson (dark matter) due to the Z 2 symmetry: X µ → −X µ , and the vector gauge boson obtains a mass, m X = g X Q Φ S v S (see also refs. [23-27, 49, 50] for related discussions). In this paper, we will simply take JHEP08(2018)099 Q Φ S = 1; in other words, Q Φ S and g X are lumped together. The interacting terms of the dark sector is given by where the hidden scalar field h 3 will mix with the neutral scalars in the two Higgs doublets through the interaction given by eq. (2.1). After spontaneous symmetry breaking, the version of the Higgs sector becomes the next-to-minimal two-Higgs-doublet model (N2HDM). We decompose the Higgs doublet fields as where the vacuum expectation values (VEVs) of the doublets approximately satisfy v 2 1 + v 2 2 ≈ v 2 = (246 GeV) 2 , with v being the SM VEV. The scalar fields in Φ i and Φ S can be expressed in terms of mass eigenstates of physical Higgs states and Goldstone bosons as where (H, h) are the (heavy, light) Higgs CP-even scalars in the two Higgs doublets in the limit of δ, θ → 0, A the CP-odd scalar, H ± the two charged Higgs bosons, and (G ± , G 0 ) the Goldstone bosons corresponding to the longitudinal components of W ± and Z, respectively. Here β is the mixing angle of the charged bosons, and α is the mixing angle of neutral CPeven bosons (h 1 , h 2 ) in the limit of δ, θ → 0, where the former is defined as tan β = v 2 /v 1 . The theoretical requirements for the perturbativity, vacuum stability, and tree-level perturbative unitarity are given in appendix A. From the results for the square of masses of H ± and A, square of mass matrix of H, h and S, and the minimum conditions of the Higgs potential at the VEV, the quartic couplings λ i , with i ≡ 1, . . . , 8, can be rewritten in terms We show the relations in appendix B.

The Yukawa sectors
The type-X Yukawa interactions are imposed a Z 2 symmetry only to the right-handed quarks, u R → −u R and d R → −d R . Thus, the Yukawa Lagrangian, describing the interactions of the Higgs doublets to the SM fermions, is given by Table 1. The tree level Yukawa couplings of the neutral type-X N2HDM Higgs bosons, keeping terms linear in sin θ and sin δ, with respect to that of the SM Higgs.
whereΦ 2 = iσ 2 Φ * 2 , and y i is the 3 × 3 Yukawa matrix. In terms of the mass eigenstates of the scalar bosons, the Yukawa interaction terms can be rewritten by Keeping small terms linear in sin θ and sin δ, the Yukawa couplings of the neutral Higgs bosons in the type-X N2HDM, normalized with respect to the SM Higgs, are given in table 1, where For the type-X Yukawa interactions, the normalized lepton Yukawa coupling to the SM Higgs is given by Considering the LHC data but without muon g − 2 constraint, the allowed parameters, consistent with the alignment limit of s β−α → 1, lie in two different regions [43]; in one region, the g h couplings (→ 1) have values near the SM ones, while in the other region which is called the wrong-sign region, the g h → −1 has opposite sign to the SM Higgs couplings to V V , (normalized) g hV V,hZZ = s β−α → 1, and to the quark pair, g hf f = c α /s β = s β−α + c β−α /t β → 1 for a large tan β satisfying 2t β c β−α ∼ 2. Only the wrong-sign region is favored by the muon g − 2 measurement [43] (see also the following discussion in this work). On the other hand, for the type-X Yukawa interactions, the couplings, g A and g H ∝ tan β, are enhanced by a large tan β, while g Aqq and g Hqq ∝ 1/ tan β are suppressed, where q ≡ quark. Therefore, as shown in figure 1, the two-step cascade DM annihilation process via the on-shell pseudoscalar boson into the SM particles are dominated by τ 's in the final states for a large tan β. Note that the DM annihilation into tau's cannot be through the heavier on-shell neutral Higgs, which is kinematically forbidden, because, as shown in this work, its mass m H ∼ 300 GeV is much larger than the DM mass. Note also that, in contrast with the type-X case, for the type-II Yukawa interactions, because both the JHEP08(2018)099 Figure 1. The main DM annihilation process of two-step cascades, X µ X ν → SS → 4A's → 8τ 's, relevant to the GC gamma-ray excess. The red shaded region contains interactions shown in figure 2. down-type quark and lepton couplings of the heavier neutral Higgs boson are enhanced by tan β, that model will be severely constrained by the extra Higgs search at the LHC and by the flavor physics [51].
In the present work, we study that the vector DM (X) first annihilates into the unstable hidden Higgs bosons (S), as shown in figure 2, and then the S dominantly decays into the pseudoscalar pair. In the following section, we will give the triple and quartic Higgs couplings, which are relevant to the XX → SS and S → AA processes. Moreover, these couplings are also relevant to the Boltzmann equations, which will be discussed in section 5.2.

The triple and quartic Higgs couplings
We consider that the GC gamma-ray excess originates from the two-step cascade DM annihilation, so that the mixing angles, δ and θ, are small, and only terms linear in sin δ and sin θ are kept in the effective couplings. However, because tan β needs to be larger (∼ 35) in the present case, we also keep the terms, which are quadratic in these two angles and involve tan β. The constraints on δ and θ will be discussed in section 5.
The Lagrangian, containing triple and quartic Higgs couplings, are relevant to the present study. Using the result given in appendix B, we can express these couplings in where, neglecting all terms suppressed by s 2 δ , s 2 θ , and s δ s θ except that enhanced by t β , the couplings are given by Here and in the following, we adopt the abbreviations: s σ ≡ sin σ, c σ ≡ cos σ, and t σ ≡ tan σ.

The decay widths for the S and A bosons
The partial decay widths for S and A are relevant to the studies of the relic density and indirect detection searches. The partial decay widths of the S boson are 25) where N f c ≡ 3 (1) for quarks (leptons), and f S ( The partial decay widths of the CP-odd boson A are given by with the couplings g Aτ τ = g Aµµ = tan β, |g Aqq | = 1/ tan β in the consideration of the type-X Yukawa interactions. In the present case, because we take into account the large tan β(∼ 35) and m A = 15 ∼ 20 GeV (see later discussions), we therefore have Br(A → τ τ ) 1, Br(A → µµ) Br(A → τ τ ) × (m µ /m τ ) 2 0.0035, and neglect A →qq and A → gg due to the 1/ tan 2 β suppression in the decay rates.

Use of the parameters under the experimental and theoretical constraints
For the CP-conserving LN2HDM, we adopt the observed Higgs resonance as one of the CP-even scalars: h with mass m h = 125.09 GeV [52], v 2 ≡ ( √ 2G F ) −1 (246 GeV) 2 . Compared with the SM, the interactions contain 11 more independent parameters. We take the following remaining parameters as inputs: In this parametrization, the tree-level couplings of the Higgs bosons to SM particles are functions of tan β, β −α, θ, and δ. In the following, we will experimentally and theoretically constrain the parameters relevant to the two-Higgs doublets and Yukawa sectors.

Experimental considerations
In the present paper, we consider that the GC gamma excess is mostly due to the two-step cascade annihilation of the dark matter into the final state τ 's, for which the main process is schematically shown in figure 1, where the shaded region denotes the interactions, depicted in figure 2, and is relevant to the DM annihilation cross section. As shown in the following GC gamma-ray excess study that m A ∼ 10−20 GeV < m h /2, we therefore need to consider the constraint on Br(h → AA), which is proportional to the square of λ hAA . The magnitude of λ hAA can be constrained from the measurement of Br(h → AA → 4τ ) Br(h → AA), of which the current upper bound [53,54] is about 0.2-0.4 for 8 ≤ m A ≤ 30 GeV, resulting in |λ hAA | < 1.95 × 10 −2 . Following ref. [44], we will take λ hAA = 0, i.e., Using the results in eqs. (3.3) and (3.4), the normalized Yukawa coupling of the SM Higgs to the lepton pair can be expressed as In this case, the alignment limit, s β−α → 1, reproduces the wrong-sign SM coupling g h → −1. Taking the combination of the ATLAS and CMS h → τ τ measurements, the signal strength reads µ τ τ ≡ (σ h · BR) obs τ τ /(σ · BR) SM τ τ = 1.11 +0.24 −0.22 [55], which is defined as the observed product of the SM-like Higgs production cross section and the decay branching ratio h → τ τ , normalized to the corresponding SM value. The corresponding requirement for m A ≤ 20 GeV is |g hτ τ | < 1.26 at 2σ confidence level (C.L.), such that we have M 245 GeV.
The masses of Higgs bosons can be constrained by the electroweak precision measurements. Such new physics effects, which contribute the gauge vacuum polarization at the one-loop level, can be described by three oblique parameters, S, U and T . We adopt the definition of these parameters, originally introduced by Peskin and Takeuchi [56,57]. Taking the limit s β−α → 1, |m H ± − m H | m H and m A m H , and keeping terms linear in sin θ and sin δ, we obtain three oblique parameters from that given in ref. [58], where a general multi-Higgs-doublet model was studied. The results are collected in appendix C. In the limit that we take, the formulas are consistent with those in the two-Higgs doublet model, i.e. the correction due to the hidden Higgs boson S is negligible, and the results

JHEP08(2018)099
approximately read where the T parameter is especially sensitive to the mass splitting, m H ± − m H . For the values m H ≈ 300 GeV and |m H ± − m H | ∼ O(10) GeV, the theoretical prediction is consistent with that from the data fit which gives [73] S = 0.05 ± 0.10, T = 0.08 ± 0.12, U = 0.02 ± 0.10. (3.7)

Theoretical considerations
For this model, we need to have m H ∼ m H ± ∼ M m A . To satisfy the perturbative bound, we impose the absolute values of all the quartic couplings to be less than 4π. We can easily make the estimate on the mass bound for the heavy Higgs as follows. From eq. (B.5), GeV. For small mixing angles, θ and δ, which are relevant to the present work, the tree-level perturbative unitarity, as the case of the type-X 2HDM, gives m H ± 700 GeV [44]. On the other hand, the vacuum stability and perturbativity could be broken when we consider this model at higher scale, for which, again, in the limit of small θ and δ, the related bound is the same as the type-X 2HDM, and given by m H ± (400) 310 GeV for the cutoff scale Λ (10) 100 TeV [22,44].
Neglecting the terms with power higher than that linear in s θ and s δ , and taking the limit s β−α → 1 and t β 1, we have, from eq. (B.1), that where the second line is obtained by using the relations given in eqs. The gamma-ray spectrum originating from the two-step cascade dark matter annihilations: determining m X , m S , and m A The differential gamma-ray flux, arising from the two-step cascade annihilations of the vector DM, can be expressed by where the J-factor is the integral of the DM density squared along the line of sight (l.o.s.) and over the solid angle ∆Ω that covers the region of interest (ROI), and σv f and (dN f γ /dE) X are the low-velocity averaged annihilation cross section and the gamma-ray spectrum produced per annihilation with final state f , respectively. For illustration, the dominant process is depicted in figure 1, where the final states are τ 's, which mainly arise from the process, σv τ σv XX→SS ×Br(S → AA)×Br(A → τ τ ), with Br(S → AA) 1 and Br(A → τ τ ) 1. Following the method given in ref. [41], we can perform two-step Lorentz boosts to transform the gamma-ray spectrum given in the A boson rest frame, (dN τ γ /dE) A , to the XX center of mass (CM) frame; we first boost the spectrum to the S rest frame and then to the CM frame of the XX pair. For (dN τ γ /dE) A , we will use the PPPC4DMID result [59,60], which was generated by using PYTHIA 8.1 [61]. Thus, with E, E 1 , and E 0 being the photon energies in the XX CM frame, S rest frame, and A rest frame, respectively.

The Galactic center gamma-ray excess
We use the GC gamma-ray excess spectrum obtained by CCW [8], who have studied Fermi-LAT data covering the energy range 300 MeV-500 GeV in the inner Galaxy, where the ROI extended to a 40 • × 40 • square region around the GC with the inner latitude less than 2 • JHEP08(2018)099 masked out. The systematic uncertainties of the background have been taken into account by CCW through a large number of Galactic diffuse emission models. To see whether the present vector DM model can meet the observation, we perform a goodness-of-fit test, by calculating the χ 2 test statistic, where 24 energy bins are adopted in the range 300 MeV−−500 GeV, dΦ γ /dE i and (dΦ γ /dE i ) obs are the model-predicted and observed flux in the ith bin, respectively. Here the covariance Σ ij contains the uncorrelated statistical error, and correlated uncertainties, of which the latter is composed of the empirical model systematics and residual systematics. Although CCW performed the analysis using the older Fermi dataset, however, it was shown in ref. [62] that the results have very little changes between Fermi Pass 7 and newer Pass 8 data. This appreciable difference at low energies might be due to the modeling for the point sources in various datasets [11,62]. On the other hand, it is interesting to note that the central values of the low energy spectrum given by [11] seem to be smaller than that obtained by CCW. If so, the best-fit DM mass will become larger compared with the present result. Two physical parameters, σv and m X , can thus be obtained from the fit. The value of σv is sensitive to the form of the Galactic DM density distribution, for which we use a generalized Navarro-Frenk-White (gNFW) halo profile [63,64], where the scale radius r s = 20 kpc, r is the distance to the GC, −γ is the inner log slope of the halo density near the GC, and ρ is the local DM density at r = 8.5 kpc, which is the radial distance of the Sun from the GC. We will take γ = 1.2 and ρ = 0.357 GeV as the canonical values. However, the uncertainties about the local dark matter density and the halo distribution near the GC remains large. The resulting annihilation cross section in the fit due to the variation of γ ∈ [1.1, 1.3] and ρ ∈ [0.2, 0.6] GeV will be discussed later.

The constraint from dwarf spheroidal observations
In the present analysis, we will use the combined gamma-ray data of 28 confirmed and 17 candidate dwarf spheroidal galaxies (dSphs), recently reported by the Fermi-LAT and DES Collaborations [65,66]. Compared with the earlier Fermi-Lat analysis [67], where some point-like sources were modeled as extended ones, a consistent analysis across 45 targets were presented in ref. [65], and a limit weaker by a factor of ∼ 1.5 were obtained in the low DM mass region ( 70 GeV). Because there is no gamma-ray signal detected so far from this kind of objects, a bound on the DM annihilation can thus be set. We perform a combined likelihood analysis of 45 confirmed and candidate dSphs with 6 years of Fermi-LAT Pass 8 data in the energy range from 500 MeV to 500 GeV. The

JHEP08(2018)099
log-likelihood test statistic (TS) is given by with N dSph = 45 and the profile likelihood of an individual target k, where N bin = 24 are the numbers of bins, L ki is the i-th binned likelihood of the target k [66], and the J-factor likelihood for a target k is modeled by a normal distribution [68], Here, J k is the expected J-factor of a target k, while the nominal value J o,k together with its error σ k is the spectroscopically determined value when possible, or the predicted one from the distance scaling relationship with an uncertainty of 0.6 dex, otherwise [65]. For a given m X , σv andJ k are the maximum likelihood estimators (MLEs), which maximize k=N dSph k=1 ln L k . When σv is fixed to a given value,Ĵ k are the conditional MLEs of the nuisance parameters. We can obtain the 95% confidence level (C.L.) limit on low-velocity annihilation cross section σv from the null measurement by increasing its value from σv until TS = 2.71.

Results
In figure 1, we have depicted the two-step cascade DM annihilation process into the final state τ 's, which is the dominant mechanism to explain the GC gamma-ray excess in the present study. The DM annihilation diagrams, relevant to the indirect search and also to the relic abundance, are shown in figure 2; the resulting cross sections and related discussions are collected in appendix D.
The DM two-step cascade annihilation process, increasing the final gamma-ray multiplicity and therefore resulting in a broader gamma-ray spectrum, provides a better fit to the GC data compared with that obtained from the DM annihilation directly into the tau pair. For illustration, in figure 3, we show the GC gamma-ray energy spectrum [8] compared our the best-fit model prediction. The corresponding GC fitted regions, featuring by the p-values, are shown in figure 4 on the plane of m X and low-velocity annihilation cross section for four cases of (m S , m A ) = (0.5m X , 0.2m X ), (0.7m X , 0.3m X ), (0.9m X , 0.36m X ), and (0.95m X , 0.45m X ), where ρ = 0.357 GeV/cm 3 and γ = 1.2 have been adopted, so that the best GC fit is consistent with the WIMP relic abundance. The gauge coupling constant g X in this model can thus be determined, and given as a function of m X in figure 5. Note that, because the decay width of S is much less than m S and, on the other hand, the annihilation XX → SS via an s-channel h or H exchange is highly suppressed, these two effects can be negligible. The detailed discussion will be given in section 5.1. We find that this model can provide a good fit to the GC gamma-ray excess spectrum for the regions JHEP08(2018)099 10 -5 Figure 3. Comparison of the spectrum of the Galactic Center gamma-ray excess [8] vs. the best-fit results for DM annihilating through a two-step cascade to the final state τ 's, where, for the former, the statistical and systematical errors are shown by error bars and orange rectangles, respectively, while for the latter, the solid (blue), dashed (brown), dotted (magenta), and dotdashed (green) curves are for the cases of (m S , with Compared with the DM annihilation directly into the tau pair, the two-step cascades increase the final gamma-ray multiplicity by a factor of ∼ 4, such that if the dark matter mass is still the same, the resulting energy of the GeV photon peak will be reduced by a factor ∼ 4. Therefore, to fit the observed GeV gamma-ray excess, we need to increase m X , i.e. the initial energy, by a factor ∼ 4 in magnitude. On other hand, having the resulting changes for the final photon multiplicity and m X , we thus know from eq. (4.1) that the annihilation cross section also needs to be enlarged by a factor of ∼ 4 to fit the gamma-ray spectrum.
Figures 4 and 5 show the current 95% C.L. upper bound and projected limit from the gamma-ray observations of dSphs. The projected limit approximately rescales with the square root of the data size and the square root of the number of targets [69]. Following the estimate given in ref. [70], we conservatively assume that the 15-yr gamma-ray emission data can be successfully collected from the observation of 60 dSphs. Thus, the projected sensitivity on the σv will be further improved by a factor of ∼ 1.83, and, as shown in figures 4 and 5, the present model is very likely to be probed in the near future.
In figures 4 and 5, we also show the region which is allowed by the correct DM relic abundance in the conventionally thermal WIMP scenario, for which we have rescaled the thermally averaged annihilation cross section at freeze-out temperature to its corresponding value defined at the low-velocity limit, where the effective number of degrees of freedom  (DoF) g * 85.5 corresponding to T 1.9 GeV (and x ≡ m X /T 22) is adopted [71,72]. Note that for a too small coupling constant λ SAA , the S particle cannot maintain its chemical equilibrium with the thermal bath, such that the dark sector particles to be out of equilibrium with the bath when T m X(S) . For this case, we will show in section 5.2 that the allowed DM annihilation cross section to have the relic abundance could be (much) larger than that in the conventional WIMP scenario.
Three remarks for the gamma-ray fit are in order here. First, for the two-step cascade process, the kinematical condition, m X ≥ m S ≥ 2m A ≥ 4m τ , needs to be satisfied. Second, considering variation of the local dark matter density from 0.357 to 0.2 (or 0.6) GeV/cm 3 and the halo slope γ from 1.2 to 1.1 (or 1.3), the low-velocity annihilation cross section and g X in the GC gamma-ray fit would be further raised (or lowered) by factors of 4.10 and ∼ 4.10 1/4 (or 0.27 and ∼ 0.27 1/4 ), respectively. Third, the uncertainties of the observed J values of dSphs are subject to determination of DM mass profile which is assumed to be spherically symmetric and to have negligible binary motions [68].  Figure 5. Same as figure 4, but in the (m X , g X ) plane, where the blue shaded region delineated with the blue line provides a good fit to the GC gamma-ray data with p-value ≥ 0.05.
5 Determination of mixing angles, θ and δ

Constraints from invisible Higgs decays and two-step cascade annihilation
Taking into account the GC gamma-ray excess which results from the two-step cascade annihilations of the vector dark matter into the final state τ 's in this model, we have found that the mass of dark matter lies within 25 ∼ 50 GeV, along with m S m X and m A m S /2, as shown in figures 4 and 5. Therefore, because m h > 2m S , the decay h → SS is allowed and is followed by S → AA and subsequently A → τ + τ − . This is relevant for search for the exotic Higgs decay with 8 τ 's in the final state. The exotic Higgs decays with 4 τ 's or other modes in the final state was discussed in ref. [54]. Here and in the following sections, we will take parameters, as a benchmark in the discussions, and we have λ hAA 0, Br(h → AA) 0 due to the adopted values of tan β and β − α. In figure 6, we show the contour plot for Br(h → SS) on the (θ, δ) plane, in comparison with a 95% C.L. limit: Br(h → beyond SM) < 34%, which was fitted with the Higgs produced via SM couplings [73].
A more stringent constraint can be obtained by requiring that the two-step cascade annihilation to the final state τ 's is dominant over the one-step cascade process described JHEP08(2018)099 by the s-channel XX → AA followed by A → τ + τ − . For this requirement, we will therefore restrict σv XX→AA / σv XX→SS be to less than 0.05. Because the s-channel contributes about 15% to the total XX → SS cross section, we thus need to have |λ SAA /λ SSS | 0.48 and |λ SAA | 0.022 (see also figure 6), where λ SSS = −3g X m 2 S /(m X v) −0.046. Imposing the constraints required by Br(h → beyond SM) < 34% and |λ SAA | 0.022, we then obtain |θ| 0.001, |δ| 0.088. (5.2) There may exist points with λ SAA 0 on the (θ, δ) plane for δ = 0. If we avoid this tiny region around that points, where the other decay modes are also highly suppressed due to the very small values of θ and δ, we always have Br(S → AA) 1. In figure 6, the contour plot for the S → AA decay width is shown on the (θ, δ) plane. In this (two-step cascade annihilation dominant) case, we have Γ S /m S 6.7 × 10 −4 for |θ| 0.001, where Γ S is the total width of the S boson; the value of the width can thus be negligible in the calculation. On the other hand, compared with XX → SS via an s-channel S exchange, as shown in figure 2(b), if the mediator is replaced by h (or H), the resulting cross section is suppressed not only by the propagator of the heavier boson but also by the couplings squared: (s δ λ hSS /λ SSS ) 2 (or (s θ λ HSS /λ SSS ) 2 ), where s δ (or s θ ) comes from the X-X-h (or X-X-H) vertex, and the coupling ratios λ hSS /λ SSS and λ HSS /λ SSS are shown in figure 6. Therefore, the annihilation XX → SS via a heavier mediator, h or H, is negligible in the calculation; the conclusion is also valid for the 0-step cascade annihilation via an s-channel h (or H) exchange since these cross sections are also suppressed by a factor of s 2 δ (or s 2 θ ) resulted from the X-X-h (or X-X-H) vertex.
In figure 6, we show the contour results of λ SAA and λ SAA /λ SSS on the (θ, δ) plane. As indicated in the relevant (θ, δ) region, the λ SAA is much more sensitive to the variation of θ, compared with its dependence on δ. In the following analysis, the δ is simply set to be zero, and the dependence of the DM relic density on the effective coupling can be related to the variation of θ. If taking δ = 0, the constraint from the two-step cascade annihilation gives |θ| 0.00043. Our conclusion can be easily extended to the case with δ = 0.

JHEP08(2018)099
and n S , for X and S, respectively, are described by the coupled Boltzmann equations, where n eq i is the equilibrium number density for the particle "i", and σv and Γ are respectively the thermally averaged cross section and decay rate, corresponding to the specific process denoted in the subscript. For the present vector DM+type-X N2HDM, the interaction between the CP-odd Higgs boson A and lepton τ , via A ↔ τ + τ − and AA ↔ τ + τ − , is strong enough to maintain the chemical and thermal equilibrium between A and the SM particles in the early Universe until the DM is completely freeze-out, i.e., n A = n eq A ; in other words, during the relevant epoch, like other SM particles, the A boson is in thermal equilibrium with the reservoir. On the other hand, the interaction strength between S and A depends on the effective coupling λ SAA , which is a function of θ and δ (see also figure 6).
To solve eqs. (5.3) and (5.4), we define the normalized yields, where Y X ≡ n X /s and Y S ≡ n S /s are respectively the dark matter and mediator number densities normalized by the total entropy density, g is the effectively total number of relativistic DoF, and x ≡ m X /T is the variable that will be used instead of time t. Here g eff and h eff are the effective DoF for the energy density and entropy density, respectively [71]. Using the new defined quantities, we can rewrite these two Boltzmann equations into the following forms where M pl = (8πG) −1/2 = 2.44 × 10 18 GeV is the reduced Planck mass, and the values of y i in equilibrium are given by

JHEP08(2018)099
with the subscript index "i" ≡ X or S, and g X = 3 and g S = 1 being the internal degrees of freedom of the X (dark matter) and S (mediator) particles, respectively. Using the relations, which are inferred from the Boltzmann equation [74], (y eq X ) 2 σv XX→SS = (y eq S ) 2 σv SS→XX , (5.10) where Γ S→AA is the S → AA decays width given in the rest frame of S, and K i is the modified Bessel function of second kind, we can further recast eq. (5.8) into the form, Because both the annihilation processes, SS → AA and XX → SS, occur through the swave, for simplicity, in the following discussion we approximate σv SS→AA σv (0) SS→AA and σv XX→SS σv (0) XX→SS using their leading values in v → 0 limit. Note that in the Boltzmann equation, the contribution due to the SS ↔ AA interaction is usually much smaller than the process S ↔ AA; in other words, in eq. (5.12), the second term of the right hand side (r.h.s. ) is negligible, especially when the co-decay occurs with |λ SAA | 1.03 × 10 −8 , corresponding to |θ| 2 × 10 −10 if δ = 0.
After DM freeze out, so that y 2 X (y eq X · y S /y eq S ) 2 , the Boltzmann equation given in eq. (5.7) can reduce to Integrating this approximate equation from the freeze-out epoch x f (= m X /T f ) until very late times x ∞ ( x f ), one can obtain where y ∞ y f . In the following study, along with δ = 0, we use the parameters given in eq. the corresponding mixing angle θ satisfies 7 × 10 −10 < |θ| < 0.00043 if taking δ = 0. As an example shown in figure 7 (a), using (θ, δ) = (0.00043, 0), which corresponds to λ SAA −0.022, we obtain x f = 22. For the values of |λ SAA | 0.022, we have found that σv XX→AA / σv XX→SS 0.05, so that the GC gamma-ray excess originating from the two-step cascade annihilation to the final state τ 's is still dominant over the one-step cascade process described by XX → AA followed by A → τ + τ − .
(ii) The unconventional WIMP dark matter. If |θ| is less than 7 × 10 −10 but still larger than 2 × 10 −10 , the nonrelativistic dark sector particles start to decouple from the thermal bath as for x ∼ 1. However, for this case, the dark sector can reach again thermal equilibrium with the reservoir before the DM freeze out, so that the dark matter is still WIMP-like, and has the same freeze-out temperature and thermally averaged annihilation cross section as the WIMP case. As an example, using (θ, δ) = (2 × 10 −10 , 0), which corresponds to λ SAA −1.03 × 10 −8 , we show the result in figure 7 (b). See also figure 9.
(iii) The co-decaying dark matter. This scenario, for which the corresponding mixing angle is |θ| < 2×10 −10 if taking δ = 0, is characterized by y S y eq S and y 2 X (y eq X ·y S /y eq S ) 2 at x = x f . For this type of the scenario, featuring by a small |λ SAA | 1.03 × 10 −8 , when the dark sector particles become to be nonrelativistic (x ≈ 1), they start to decouple from the thermal reservoir, and their total yields, y X + y S , (which is related to the total number density of the dark sector in the co-moving frame) tend to remain constant. Then after a time interval, ∆t Γ ≈ 2H −1 (during the radiation dominated epoch) ∼ Γ −1 S , for a degenerate case m X ≈ m S , the dark matter X as well as the mediator S undergoes an exponential decay until freeze-out, of which at the temperature x f , the Hubble expansion rate H becomes larger than the XX → SS annihilation rate. This process is described by the "co-decaying dark matter" mechanism which was proposed by Dror, Kuflik, and Ng [46] for the degenerate case. In the following, we will further discuss a generic case, including the degenerate and non-degenerate ones. We will show that if ∆t Γ is large enough, the exponential suppression for y X can occur much earlier than y S due to a significantly suppressed up-scattering rate for the process SS → XX, such that the DM freeze-out time is earlier than the time that S undergoes an exponential decay. See also figures 7 (c)-(d).
We then further estimate the value of ∆t Γ for the co-decaying DM. There are two different cases: (i) m X = m S , and (ii) m X > m S . Considering the case of m X = m S , we sum eqs. (5.7) and (5.12), and neglect the second term of r.h.s. of the latter equation, Taking the initial conditions: y X (1) = y eq X (1), y S (1) = y eq S (1), and y eq X (1)/y eq S (1) = g X /g S = 3, and approximating y X + y S 4y S , we get the solution to eq. (5.15), for which the lifetime of S is ∆t Γ 4 Γ −1 S→AA . We find that the result shown in eq. (5.17) or eq. (5.18) can be a good approximation for the case with x Γ 20, where the hidden Higgs S can have a sufficient time to satisfy the approximation y S − y eq S ≈ y S which has been taken in eq. (5.15).
For the case of m X > m S , the value of x Γ depends not only on θ but also on the mass difference of X and S. If the difference of m X and m S is sizable enough, then the down-scattering rate, XX → SS, could be significantly larger than the up-scattering rate, SS → XX. Under this condition and after a sufficient time with x Γ 20, one could have JHEP08(2018)099 y eq X /y eq S 1 and y X y S , so that the second and third terms on the right hand side of eq. (5.12) are negligible, and this equation can thus be approximated as From this equation, we get where the effective initial S yield, y in S (1), can be approximated as y X (1) + y S (1) 4y S (1), because the down scattering is larger than the upper one, and the mostly initial X could scatter into S before the S boson undergoes an exponential decay. Therefore, if m X − m S is sizable enough, we have x Γ 1 + 4.8/C and y S (x Γ )/y S (1) 4e −Γ S→AA ∆t Γ , for which the lifetime of S is ∆t Γ 2.4 Γ −1 S→AA . In general, for a co-decay process with x Γ 20, we get 1 + 4.8/C x Γ 1 + 8/C. Taking the same parameters which have been used, but adopting m S as a free parameter, we find that x Γ = 1 + 4.8/C is a good approximation provided that 2m A < m S 35 GeV. Figures 7 (c) and (d) show the values of x Γ to be about 51 and 117, respectively, where the latter satisfies x Γ ≈ x f . For a process with a much larger x Γ as shown in figures 7 (c) and (d), the inequality y S > y X , i.e. n S > n X , becomes much more noticeable and, due to a suppressed up-scattering rate, the exponential suppression for n X occurs much earlier than n S . It was stressed in refs. [47,48] that when the hidden sector is out of the thermal equilibrium with the bath, the hidden particles may undergo so-called cannibalism before DM freezeout. Here we include the cannibal interactions, XXS ↔ SS, SSS ↔ XX, XXX ↔ XS, SSS ↔ SS and XSS ↔ XS, which are number changing annihialtions, in the Boltzmann equations, dn X dt + 3Hn X = · · · − σv 2 XXS→SS n 2 X n S − (n eq X ) 2 n 2 S (n eq S )

JHEP08(2018)099
where the dots are all terms of the right hand side of eqs. (5.3) and (5.4), respectively, and σv 2 is the thermally averaged annihilation cross section for the 3 → 2 cannibal process. One can refer to ref. [75] for a general from of 3 → 2 scattering rates. In each cannibal term of the above equations, the factor (including the relative sign) in front of the cross section equals to ∆n i /N !, where 1/N ! is for avoiding the double counting for the initial number density in the reaction with N being the number of the identical particles of the initial states for the relevant cross section, and ∆n i is the number change for the hidden particle with i ≡ X for eq. (5.3) or ≡ S for (5.4); for instance, for the process SSS → XX, we have N = 3, but ∆n X = 2, ∆n S = −3, resulting in the factor 1/3 and −1/2 shown in the corresponding terms in eqs. (5.21) and (5.22), respectively.
To calculate 3 → 2 thermally averaged annihilation cross sections for the nonrelativistic dark sector particles with x 1, we take the low-velocity approximation, σv 2 σv 2 , by neglecting the correction of O(T /m X ): where the value of σv XX→SS used in the numerical analysis satisfies eqs. (5.30) and (5.31).
As shown in figure 7 (b), (c) and (d), for the case of |θ| < 7 × 10 −10 , when x 1, the dark sector particles become nonrelativistic and are kinetically decoupled from the thermal bath due to small S ↔ AA and SS ↔ AA interaction rates, so that the comoving number density of the total hidden sector particles remain constant before the time that the hidden Higgs bosons undergo the exponential decay. However, in the present case, the cannibal annihilations cannot be neglected. We take eq. (5.21) as an example to give a qualitative analysis on cannibalization as follows. The left hand side is of order n X H, while the right hand side due to XXS → SS is of order n 2 X n S σv 2 XXS→SS . Therefore, if the XXS → SS reaction rate is much larger than the expansion rate, the only way to maintain the equality of eq. (5.21) is to have n X = n eq X and n S = n eq S , about which one can obtain the same conclusion using either SSS → XX, XXX → XS, SSS → SS, XSS → XS, or eq. (5.22) in doing a similar analysis.
Using the same parameters as in section 5.2.1, and further including the 3 ↔ 2 interactions in the Boltzmann equations, in figure 8 we show the results of the normalized yields. Compared with figure 7 (c) and (d), at x Γ , which has been obtained in eq. (5.20), the corresponding value of y S (and also the number density of S) reduces 2 orders of magnitude due to the cannibal effect. Figure 8 In concluding this section, we would like to discuss the relations and constraints between the parameters and observables. The XX → SS annihilation cross section and x f can be related to each other via the DM relic abundance. The dimensionless density parameter of the present-day DM relic abundance, determined to be Ω DM = (0.1198 ± 0.0026)/h 2 from the observations [73,76], is given by [71] where Y ∞ X , related to y ∞ X via eq. (5.5), is the post-freeze-out value of Y X , ρ c = 3H 2 0 /(8πG) is the present critical density, h 0.678 is the present-day Hubble constant H 0 in units of 100 km s −1 Mpc −1 , s 0 = 2891 cm −3 is the present-day entropy, and In figure 9, we show x f as a function of θ, where the relations between x f and the XX → SS annihilation cross section σv , and between θ and Γ(S → AA) are depicted. For simplicity, we have approximated the thermally averaged DM annihilation cross section by using its leading s-wave value in v → 0 limit. Note that the vertical axis for the annihilation cross section has been rescaled by use of the temperature-dependent g * given in ref. [72]. Moreover, in this lepton-specific (type-X) N2HDM portal vector dark matter model, the total width of S satisfies Γ S Γ(S → AA) with taking δ = 0.
In figure 9, if θ is less than 0.00043, denoted by the vertical dotted (blue) line, the GC gamma-ray annihilation is dominated by the 2-step cascade DM annihilation. For The range of θ between the two vertical dashed (red) lines denotes that the DM X as well as the mediator S can be again in thermal equilibrium with the thermal bath before freeze out. The horizontal dot-dashed (purple) line shows the maximum value of the annihilation cross section that could still account for the GC gamma-ray excess. The 95% C.L. upper limit and project sensitivity for observation of dSphs are denoted as horizontal dashed (green) and long-dashed (brown) lines, respectively. For a θ smaller than that denoted by the vertical solid (magenta) line, we have x Γ > x f . the exponential suppression for y X can occur much earlier than that for y S . Moreover, the cannibal annihilations also reduce the value of x f . Therefore, as shown in figure 8(d), we can have x Γ x f . In figure 9, for θ < 7.5 × 10 −12 denoted by vertical solid (magenta) line, we have x Γ > x f . It should be noted that the approximation for the temperature dependence of DoF given in ref. [72] breaks down during the QCD phase transition which may occur in the temperature range of 150 − 400 MeV, i.e. corresponding to x ∼ 100 − 267 for m X = 40 GeV. For simplicity, the freeze-out results shown in figures 7 (d) and 9 are assumed to occur before the QCD phase transition.
We then further discuss the big bang nucleosynthesis (BBN) and cosmic microwave background (CMB) constraints. The requirement of avoiding the n/p ratio and 4 He abundance to deviate from the standard BBN [77] is to have Γ −1 S 1 sec, which imposes a quite relaxed bound of |λ SAA | 2.7 × 10 −13 , |θ| 5.3 × 10 −15 with taking δ = 0. As for the Planck result for the CMB [78] which provides a probe of the DM annihilation at the time of recombination, t CMB ∼ 380,000 yrs, the resulting bound is highly insensitively dependent on the number of cascade steps (see figure 11 of ref. [42]), and gives σv (7 − 19) × 10 −26 cm 3 /s. The current CMB constraint is modestly weaker than that given by the observations of dSphs. Note that by varying the local dark matter density from 0.357 GeV/cm 3 to 0.2 GeV/cm 3 and the halo slope γ from 1.2 to 1.1, the annihilation cross section allowed by the GC data fit, shown in figure 4, would be further extended by a factor of ∼ 4.10 upward. In figure 9, the horizontal dot-dashed (purple) line depicts the 95% C.L. upper limit on the annihilation cross section ( σv 1.15 × 10 −25 cm 3 /s), corresponding to x f 116, to account for the GC gamma-ray excess due to variations of ρ and γ. In figure 9, we also show the current constraint from the observations of the dSphs, and the projected sensitivity for observations of 60 dSphs with 15-year data collection. The dSph projected limit on σv is likely to be improved by a factor of ∼ 1.83.

Discussions and conclusions
Before making conclusion, we would like to discuss the parameter space in favor of the measurement of the muon anomalous magnetic moment a µ ≡ (g − 2)/2 [73,79], which can especially constrain the parameters, m A and tan β. The discrepancy between the experiments and SM prediction is more than 3σ confidence level [80,81]: In this secluded DM model, the mixing angles between the hidden Higgs boson and the visible two-Higgs doublets are very small, such that the dominant contributions to a µ are from the visible sectors. It was shown in ref. [82] that the two-loop Barr-Zee diagrams can give sizable contributions to a µ . For our model, in contrast with the one-loop result, the two-loop diagram containing an A propagator with m A being O(10) GeV gives positive contributions to a µ and its magnitude is even larger than that of one-loop. We collect the relevant a µ calculations in the lepton-specific N2HDM in appendix E.1.
In figure 10, we show the g − 2 allowed regions at 1σ C.L. and 2σ C.L. on the m A -tan β plane vs. 95% C.L. and 99% C.L. upper limits from the lepton universality and τ decays, where we impose λ hAA = 0 to have Br(h → AA) = 0, and then show the results for two different values of √ M 2 = m H = m H ± = 300 and 400 GeV. A larger m H can suppress the JHEP08(2018)099 negative two-loop correction, whereas the output seems to be insensitive to variation of its value. We also take m S = 35 GeV as input. However, because the result is insensitive to the values of the θ and δ due to their smallness, the hidden Higgs boson is secluded from the SM related phenomenology, just as the present case; thus the g − 2 result that we have shown here is basically consistent with that in the type-X 2HDM. The τ decays and lepton universality, which was stressed in ref. [44], can provide a stringent bound on parameter space favored by g − 2. The relevant formulas and data are summarized in appendix E.2. As shown in figure 10, the result, capable of accounting for the muon g − 2 anomaly at the 2σ level but without violating the constraint from the lepton universality and τ decays at 95% C.L., favors the parameter space in which 8 m A 37 GeV and 30 tan β 55 for m H ∼ 300-400 GeV.
It is interesting to note that the B s → µ + µ − data disfavor m A 10 GeV for m H = m H ± = 300 GeV, but the constraint is less restrictive for larger m H and m H ± [44]. Moreover, a precise measurement can also provide a stringent constraint on m A and M 2 (see eq. (3.5) for reference).
We give our conclusion as follows. We have built a hidden sector dark matter model with dark discrete Z 2 symmetry which is the remnant of the spontaneously broken dark U(1) dm gauge group. After symmetry breaking, the hidden sector contains vector dark matter, along with a real hidden Higgs mediating the DM interactions to the visible sector through the mixing among the neutral Higgs bosons in the model, while the visible sector is the lepton-specific 2HDM.
In this model, the GC gamma-ray excess is mainly due to the 2-step cascade annihilation describing that the vector dark matter particles annihilate to the pairs of hidden scalars, S, which subsequently undergo the decays S → AA and then A → τ − τ + . We find the parameter space with m X ∼ 25 − 50 GeV, m A ∼ 3.6 − 25 GeV, and m X m S 2m A can provide a good fit to the GC gamma-ray excess spectrum. The obtained mass of the CP-odd Higgs A in the GC excess fit can explain the muon g − 2 anomaly at the 2σ level without violating the stringent constraints from the lepton universality and τ decays.
Provided that the GC gamma-ray excess is generated by the 2-step cascade annihilation, the interplay of the hidden sector and visible sector mainly arises from the interaction S ↔ AA which is dominant over SS ↔ AA, where the interaction coupling depends on the mixing angles between the hidden Higgs boson and visible ones. We have shown that three different freeze-out types for the DM relic abundance, depending on such mixing angles, may occur, and correspond to the magnitudes of coupling constant λ SAA to be (i) in between 3.60 × 10 −8 and 0.022, (ii) in between 1.03 × 10 −8 and 3.60 × 10 −8 , and (iii) less than 1.03 × 10 −8 .
For the type (i), the hidden Higgs (mediator) can be in chemical and thermal equilibrium with the bath, and the DM particles, consistent the thermal WIMP scenario, are Boltzmann suppressed until freeze out. For the type (ii), when the temperature of the thermal bath decreases with T m X , m S , the nonrelativistic particles in the dark sector decouple from the bath. Nevertheless, the dark sector can be again in thermal equilibrium with the bath before freeze out, so that the dark matter particles behave like WIMPs finally. See figures 8(b) and 9 for example. For the type (iii) which is a typical case of JHEP08(2018)099 the cannibally co-decaying DM, the dark sector decouples from the thermal background at the temperature below their masses, but undergoes a cannibal phase first. After that, the total comoving number density of the dark sector will not exponential suppressed until the time that the mediator decaying to the CP-odd Higgs boson pair occurs. For this type, the exponential suppression for the comoving dark matter number density could occur much earlier than that for the mediator due to a significantly kinematical suppression on the up-scattering (SS → XX) rate. See section 5.2 for details.
The vector dark matter in the hidden sector does not directly couple to the visible sector, 1 but instead annihilates into the"short-lived" hidden Higgs bosons which decay through a small coupling into the CP-odd Higgs bosons. Considering the BBN bound, the lifetime of the short-lived hidden Higgs is required to be less than 1 sec, which imposes a quite relaxed constraint on the coupling |λ SAA | 2.7 × 10 −13 . Therefore, due to the very small coupling constant, |λ SAA | < 0.022, the DM in the hidden sector is secluded from detections in the direct searches or colliders even though the resultant DM annihilation cross section in the case of type (iii) could be much larger than that in the thermal WIMP scenario, as shown in figure 9.
However, the DM annihilation signals are not suppressed in a general hidden sector model. We have shown the constraints from the observations of dSphs and from the 15-year Fermi-LAT projected sensitivity of 60 dSphs, where the projected limit on σv might be improved by a factor of ∼ 1.83. The observations of dSphs thus provide a promising way to test this hidden dark matter model in the near future.
with s being the invariant mass of the DM pair, v lab being the dark matter relative velocity in the rest frame of one of the incoming particles, and the effective coupling g SSS ≡ vλ SSS −3g X m 2 S /m X . The thermally averaged annihilation cross section σv Møl is equivalent to σv lab [71], and can be approximated as [35] σv Møl 2x 3/2 √ π provided that x(≡ m X /T ) 1, where = (s − 4m 2 X )/(4m 2 X ). For the indirect searches, assuming that the dark matter particles are locally in thermal equilibrium, we have x −1 = v 2 p /(2c 2 ), with v p the most probable speed of the DM distribution. The value of v p /c is about O(10 −5 ) in the dwarf spheroidal satellite galaxies [88][89][90][91][92] and O(10 −3 ) in our Galactic center [93,94]. In the low-velocity limit, we can have σv Møl ∼ = (σv lab ) LV by taking s → 4m 2 X , D.2 The annihilation processes, SS → AA and AA → τ + τ − In the derivation of the coupled Boltzmann equations, given in eqs. (5.3) and (5.4), one needs to take into account the chemical equilibrium among the relevant particles. The chemical equilibrium between the CP-even scalar S and CP-odd scalar A particles depends on the magnitudes of the decay width for S → AA and annihilation cross section for JHEP08(2018)099 SS → AA. The former is described in eqs. (2.23)-(2.26), and the latter is given by with g SAA ≡ vλ SAA and g SSS ≡ vλ SSS . The chemical equilibrium between A and τ depends on the magnitudes of the A → τ + τ − decay width and AA → τ + τ − annihilation cross section, for which the former is described in eqs. (2.27) and (2.28), and the latter is given by Sτ τ s (s − m 2 S ) 2 + Γ 2 S m 2 S , (D. 8) whereḡ Aτ τ ≡ (m τ /v)g Aτ τ andḡ Sτ τ ≡ (m τ /v)g Sτ τ with g Aτ τ = t β and g Sτ τ = (−c α s θ + s α s δ )/c β , as shown in table 1.
E The muon g − 2 and constraints from τ decays and lepton universality E.1 The muon g − 2 The contributions to a µ in the present type-X N2HDM are approximately by the following one-and two-loop results. The one-loop contributions are given by [95] ∆a ( with g H ± ≡ tan β, g j ≡ g jµµ for S, H and A being the normalized Yukawa couplings given in table 1, and f S (r) = f h (r) = f H (r) − ln r − 7/6 , f A (r) ln r + 11/6 , f H ± (r) −1/6 .

JHEP08(2018)099
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.