A potentially detectable gamma-ray line in the Fermi Galactic center excess — in light of one-step cascade annihilations of secluded (vector) dark matter via the Higgs portal

We show the presence of a potentially detectable gamma-ray line in the Fermi Galactic center excess in light of the secluded (vector) dark matter (DM) model in which the hidden scalar, nearly degenerate with DM in mass, mediates the interaction of the secluded DM with the Standard Model (SM) due to its mixing with the SM Higgs. We find that the parameter region mX ∈ [60, 132] GeV can provide a good fit to the Fermi Galactic center gamma-ray excess spectrum, appearing a prominent gamma-ray line with the energy ∈ [30, 66] GeV. The best fit gives mX ≃ mS≃ 86 GeV with a p-value = 0.42, so that the resultant gamma-ray line, arising from the decay of the scalar mediator into γγ, peaks at 43 GeV. We derive constraints on the annihilation cross section from the Fermi- LAT gamma-ray line search, gamma-ray observations of the Fermi-LAT dwarf spheroidal galaxies, and Planck cosmic microwave background measurement. For the secluded vector DM model, the parameter space constrained by the current XENON1T and future LUX- ZEPLIN is shown. Finally, for the mixing angle between the Higgs sectors, we discuss its lower bound, which is required by the big bang nucleosynthesis constraint and relevant to the hidden sector decoupling temperature.

1 Introduction 2 2 The model 4 3 Formulation of the gamma spectrum with prominent lines arising from one-step cascade DM annihilations 6 3.1 Basic formula of the differential gamma-ray flux originating from the DM annihilation 6 3.2 Formulations of the spectra arising from various channels of one-step cascade annihilations in the DM rest frame 7 3.2.1 Gamma-ray spectrum generated from S → SM SM 7 3.2.2 Gamma-ray spectrum generated from the three-body decay S → V V * → V f 1f2 7 3.2.3 Gamma-ray spectrum generated from S → γγ 10 3.2.4 Gamma-ray spectrum generated from the decay S → Zγ

Introduction
The existence of non-baryonic dark matter (DM) is evident from various cosmological observations and measurements [1,2]. Moreover, the majority of the matter density in our Universe is dominated by the DM. Currently, one of the favorable DM candidates is the so-called weakly interacting massive particles (WIMPs). For this scenario, dark matter, having mass of order GeV-TeV and interacting with the Standard Model (SM) particles at the electroweak scale, can give the correct relic abundance today. Meanwhile, the nonrelativistic WIMPs, following Boltzmann suppression, remains thermal equilibrium with the bath until freeze-out. However, the DM models built based on the WIMPs scenario are increasingly constrained due to the null results from the direct detection and collider experiments. Instead, a paradigm of DM was proposed to suggest that (WIMP) dark matter is secluded within one of the hidden sectors, and is very weakly coupled to the visible sector via a metastable mediator which is lighter than the DM [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. As such, the secluded DM may become undetectable or hard to detect in colliders and underground direct searches, but can still produce viable signals in the indirect experiments [14,16,17].
For the indirect DM searches, a number of studies have confirmed an excess of few-GeV gamma-rays from the region around the Galactic center (GC) and suggested that the excess emission could arise from the DM annihilation [18][19][20][21][22][23][24][25][26][27][28][29]. The signal origin of GC excess is not conclusive yet. Several interpretations, recently proposed from the astrophysical point of view, suggested that the excess can be better correlated with stellar over-density in the Galactic bulge and the nuclear stellar bulge [30,31], or described by point sources [32,33]. In this paper, we will focus on the secluded DM scenario for explaining the GC gamma-ray excess. In such a model, compared with the WIMP case of direct annihilations to the SM, a multi-step cascade DM annihilation can accommodate a higher DM mass, allow a larger cross section in the fit, and broadens the spectrum of the secondary particles.
Not only for a conventional WIMP model but also for a secluded DM model, gammaray lines are very likely to be expected at the loop level. Thus, the gamma-ray line signal directly/indirectly reveals the particle nature of the underlying theory of dark matter. Moreover, it provides a striking signature which could be clearly distinguished from astrophysical backgrounds. It is interesting to note that the direct DM annihilation to the SM Higgs pair, h h, gives a moderately good fit to the Fermi GC excess spectrum but with a p-value = 0.13 at best, as long as the produced h is approximately at rest [26] (cf. p-value = 0.17 obtained in ref. [34]). In this case, a detectable width of the gamma line with energy m h /2 62.6 GeV, about half of mass of DM, is very sensitive to the Lorentz-boost from the Higgs rest frame to the DM center of mass (com) frame. A distinguishable line signal also depends on the energy resolution of the detector.
Motivated by the above gamma line results [26,34], in this paper, we consider a secluded DM model, where the DM mainly annihilates into a pair of lighter scalar mediators, S. For this secluded DM model, a generic case of m DM m S can have a good fit to the GC spectrum. Instead, here we focus on the study of the GC gamma-ray spectrum with prominent lines which could be detectable. To have a clear gamma-line signal, we take -2 -

JHEP07(2020)148
into account the case that the scalar mediator is nearly degenerate with DM in mass. As will be shown in section 4, when the both masses of the DM and scalar mediator are about 86 GeV, the p-value of the best GC spectrum fit can be as large as 0. 42.
The mediator, a mixture of the hidden sector Higgs and the SM Higgs, has a mass lighter than the 125 GeV Higgs observed at the LHC, so that the resulting cascade decays can soften the gamma-ray spectrum to have a better fit to the peak at 1-3 GeV. For illustration, in figure 1, we show one-step cascade annihilations of the secluded dark matter into a pair of scalar mediators which subsequently decay to the SM final states. In indirect searches, the qualitative fit of the gamma-ray spectrum is relevant to the decay channels of the mediator as well as the DM mass, which determines the initially kinetic energy of the mediator and thus the boost factor for the spectrum, while the DM annihilation cross section plays as an overall factor in the fit. For determination of the DM relic abundance today, the thermodynamic evolution of the hidden sector before freeze-out depends on the strength of couplings between the mediator and SM [17]. If the couplings are small enough, the hidden sector can be kinetically decoupled from the bath before it becomes nonrelativistic. As such the freeze-out DM annihilation cross section required to give a correct relic density could be boosted above the conventional WIMP value [17,[35][36][37].
In the analysis, we will use the Fermi GC gamma-ray excess spectrum obtained by Calore, Cholis, and Weniger (CCW) [25]. 1 The result of CCW is based on the Fermi Pass 7 dataset, 2 of which the energy resolution is about 10% [39]. In the parameter plane of the DM annihilation cross section and DM mass, that is relevant to the spectral line(s) generated from the Higgs portal one-step cascade annihilation of DM, we will further show the current bounds imposed by the Fermi-LAT observations of dwarf spheroidal galaxies (dSphs) [40], by Fermi-LAT gamma-ray line search in the region around the GC [41], by the Planck cosmic microwave background (CMB) [2], and by direct detections [42,43]. The Fermi-LAT projected sensitivity with as much as 15 years of data [44] as well as the high energy resolution detectors from forthcoming experiments [45][46][47] is capable of exploring the considered parameter space. Thus, the Higgs portal scenario is very likely to be testable in the near future. See the details in section 4.
To be more specific, we will consider a simplest secluded vector dark matter model in which the vector DM interacts with the SM mainly through the scalar mediator, which is a hidden physical Higgs state resulting from an extremely small mixing angle between the dark sector scalar singlet and the SM Higgs. However, one should note that the determination of the gamma line is nothing to do with the fundamental property of DM, but is related to the Higgs portal.
The layout of this paper is as follows. In section 2, we present a renormalizable vector DM model in which the dark sector described by the U(1) X gauge symmetry contains 1 It should be note that the nature of the GC excess is under active debate. Besides the scenario that GC excess might arise from the DM annihilation, some newly proposed diffuse models could provide an even improved fit to the data by including various astrophysical phenomenologies, e.g., models correlating the excess with stellar over-density of the Galactic bulge [30,31], or with point sources [32,33]. 2 The extracted GC spectra do not have obvious difference among Fermi Pass 7 and Pass 8 datasets [38].
However, their results at low energies can have appreciable difference, depending on event selections of the point sources in various datasets.  a gauge vector boson and a complex scalar. Compared with the SM, four additional parameters, including the DM and mediator masses (m X , m S ), DM-mediator coupling constant (g dm ), and Higgs mixing angle (α), are introduced. In section 3, we outline the formulation of the gamma-ray spectrum with prominent lines, arising from a one-step cascade annihilation of DM to scalar mediators, which subsequently decay into SM particles through very small couplings, owing to the tiny Higgs mixing angle. In order to have a correct spectrum fit, for the mediator mass range m V m S 2m V , not only the usual two-body decay modes but also the three-body decay modes, S → V V * → V f 1f2 with V ≡ W, Z, need to be taken into account. Furthermore, we calculate the expected the gamma lines originating from S → γγ, γZ, where the continuum spectrum resulting from the Z decay is also considered in the S → γZ decay. In section 4, we present the main analysis. In section 5, we discuss the constraint on the mixing angle of the two scalar sectors from the thermodynamic point of view, and the scale-dependence of vacuum stability for the secluded vector DM model. We conclude in section 6.

The model
We consider the simplest abelian vector dark matter model, which is renormalizable. In this model, the vector dark matter, X, associated with a dark U(1) X gauge symmetry, interacts with the SM particles vis the Higgs portal, which originates from the mixture of the SM Higgs and the hidden complex scalar (Φ S ). In addition to the usual SM part, the relevant Lagrangian, involving the dark kinetic terms and scalar potentials, are described by Higgs doublet, and Φ S is the hidden complex scalar with a U X (1) charge assignment Q Φ S . In the following, we will simply use Q Φ S = 1. After spontaneous symmetry breaking, DM gets its mass, m X = g dm Q Φ S v S , and the CP-odd state, σ s , is absorbed to be the longitudinal component of X, where the Z 2 symmetry, X µ → −X µ and Φ S → Φ * S , is preserved, so that DM is stabilized. Under this Z 2 symmetry, all other fields are even. The scalar fields (φ h , φ s ) can be rewritten in terms of mass eigenstates of physical Higgses (h, S) as 4) and the mass term in the Lagrangian is given by 5) and the abbreviations, s α ≡ sin α and c α ≡ cos α, are used here and in the following.
In the analysis, we will use v H 246 GeV and m h = 125.18 GeV [48] as inputs, and take m X , m S , g dm and α as the four independent parameters, i.e., the remaining λ S , λ H , λ HS , and v S can be parametrized in terms of these parameters.
In figure 2, we show the branching ratios of the hidden scalar, S, in the range m S < 200 GeV, where a good fit of a photon spectrum showing gamma-line to the GC gammaexcess data can be obtained and will be discussed in the following analysis. For m S < 2m h , because all the decay widths of S are proportional to sin 2 α, its decay branching ratios are thus independent of the value of α. The relevant formulas for decay widths of the hidden scalar S are collected in appendix A.
We consider the small α region, where the DM annihilation is dominated by XX → SS, while XX → hh is negligible. More detailed discussions can be found in appendix B. Moreover, we consider only a sliver region of the masses, where (m X −m S ) m X , resulting in the produced S to be close to rest, can account for the gamma-ray line phenomenon, and the value of m X (≈ m S ) is thus determined from the Fermi GC fit. The observed spectral line width, which depends on the energy resolution of the instrument, is very sensitive to the Lorentz-boost from the S rest frame to the XX com frame; the result is relevant to the mass difference of m X and m S .
Using the low-velocity DM annihilation cross section obtained from the fit to the GC gamma-ray excess data, we can get the corresponding value of g dm in this secluded vector DM model.
As for a small mixing angle α 2 × 10 −6 , the hidden sector has been thermally decoupled from the bath before it becomes nonrelativistic, such that the resulting DM annihilation cross section that accounts for the correct relic density could be significantly boosted above the conventionally thermal WIMP value [17]. 3 Formulation of the gamma spectrum with prominent lines arising from one-step cascade DM annihilations 3.1 Basic formula of the differential gamma-ray flux originating from the DM annihilation The differential gamma-ray flux originating from the DM annihilation is given by Here, for the terms related to the particle physics, σv LV is the DM annihilation cross section into two hidden Higgs scalars in the low-velocity limit 3 (consistent with T → 0), and (dN γ /dE) X is the resulting photon spectrum produced per DM annihilation in the com frame of DM. On the other hand, the J-factor, related to the astrophysics, is the integral along the line of sight (l.o.s.) over the region of interest (ROI), which covers a rectangular solid angle ∆Ω with galactic latitude and longitude denoted by b and , respectively. For the GC data analysis, the l.o.s. described by the coordinate s is related to the distance to the GC by r = (s 2 + r 2 − 2r s cos cos b) 1/2 with r being the distance from the Sun to the GC. For the GC gamma-ray excess study, we adopt the generalized Navarro-Frenk-White (gNFW) profile [49,50] as a canonical DM density distribution in our Galaxy, where ρ is the local DM density corresponding to r = r .

JHEP07(2020)148
Below, we outline the calculation of the gamma-ray spectra generated from various S decays following the DM annihilation XX → SS.
3.2 Formulations of the spectra arising from various channels of one-step cascade annihilations in the DM rest frame For the process that DM annihilates into two on-shell mediators which subsequently decay through a small coupling into SM final states, the resulting photon spectrum (dN γ /dE) X defined in the com frame of DM can be written in terms of the spectrum (dN γ /dE ) S described in the S rest frame by considering a photon emitted at the angle θ measured from the boost axis along which we can boost the S rest frame by a relative velocity v = c 1 − m 2 S /m 2 X to obtain the result in the XX com frame [51]. The result is given by where x = E/m X , x = 2E /m S , and = m S /m X , the factor "2" on the right hand side (r.h.s.) results from the increased multiplicity due to the fact that each S decays to the SM final state f with a branching fraction Br(S → f ), and, after performing the angular integration, the second line is the convolution integral with the bounds for x , The kinematical range of the gamma-ray energy in the XX com frame satisfies (3.5)

Gamma-ray spectrum generated from S → SM SM
Considering the gamma spectrum arises from the S decay into a on-shell SM particle pair, we employ the PPPC4DMID package 4 [52,53] with the replacement of the DM mass there by m S /2 to generate the direct spectra (dN γ /dx ) S . This package, including the electroweak corrections, was obtained by using PYTHIA 8.135 [54]. Below, we consistently use the results generated by PPPC4DMID as the essential inputs to obtain spectra of the remaining channels.

Gamma-ray spectrum generated from the three-body decay
Below the V V kinematical threshold but m S > m V , with V ≡ W or Z, the hidden scalar decays into a gauge boson pair, of which one (V * ) is off-shell, resulting in S → V V * → 4 The package is also available from the website: "http://www.marcocirelli.net/PPPC4DMID.html".
In the S rest frame, the gamma-ray spectrum generated from S → V V * can be expressed as where, the former and latter terms on the r.h.s. describe gamma-ray spectra that are produced by the cascade decays of V and V * , respectively. Here, as before, we define x = 2E /m S with E the photon energy measured in the S rest frame. The gamma-ray spectrum generated from S → V V * can be obtained by convoluting the 3-body space with the spectrum arising from the cascade decay of V and V * . The relevant 3-body phase-space integral for the decay S → V V * → V f 1f2 is given by [55] The results will be briefly sketched as follows.
. For a photon emitted from the V cascade decay, the spectrum simply satisfies the relation, is the photon spectrum generated by the cascade decay of V which has energy E V with respect to the S rest frame. Changing variables, with m 1 and m 2 being the masses of f 1 andf 2 , respectively, we can recast the spectrum in the following form, where C V , which normalizes the spectrum, is given by

JHEP07(2020)148
We will simply take 1 = 2 = 0 and use (1/2) , where the subscript "PPPC" denotes the result generated by the PPPC4DMID package but with the DM mass replaced by E V ≡ ξm S . Here, the factor of "1/2" accounts for the fact that in PPPC4DMID the spectrum is generated by two gauge bosons, V . Note that the kinematic ranges ofx and x are given by 0 ≤x ≤ 1/2 and 0 ≤ x ≤ 1.
. For a photon emitted from the V * cascade decays, we have is the photon spectrum generated by the cascade decay of , we first consider the case with V * ≡ W * . Note that the charges of W have been summed in the width of S → W W * given in eq. (A.2). Above the thresholds of the following channels, the ratio of the W + * decays approximately follows + ν : where g V = T 3 /2 − Q f sin 2 θ W and g A = T 3 /2 with T 3 and Q f being the weak isospin and electric charge of f , respectively. For simplicity, we generically use V * → f 1,m f 2,m to denote the two-body decay of the virtual vector boson. We use the PPPC4DMID package to obtain the spectrum, is the gamma-ray spectrum arising from the cascade decays of the f i,m andf i,m pair in the PPPC4DMID package with the DM mass replaced by E V * /2 = (m S − E V )/2. Here, F m , depending on parameters such as N c , g V , g A and CKM matrix elements as shown above, is the relative fraction for each channel, m, which is above the threshold. Using the same notations as in eq. (3.9), -9 -

JHEP07(2020)148
we can rewrite this spectrum in the following form, where C V is the normalization factor of the spectrum as given in eq. (3.11).

Gamma-ray spectrum generated from S → γγ
We take the gamma line spectrum arising from the S → γγ decay as a simple δ-function form, in the rest frame of its parent S particle. Therefore, for the line spectrum in the XX com frame where the DM annihilates into two on-shell hidden scalars, each of which subsequently decays into two photons, the result can be written as When fitting the monochromatic line(s), which is likely much narrower than the experimental energy resolution, we account for the finite resolution of the instrument. The observed line spectrum measured by the detector at energy E(= x m X ) can be modeled by convolving the signal with a Gaussian energy dispersion, where σ is related to the detector energy resolution ξ as σ = ξ/(2 √ 2 ln 2) ξ/2.35, which is the ratio of the full peak width at half maximum to mean energy [56], and is the result given by eq. (3.18) but with x replaced by x 0 .

Gamma-ray spectrum generated from the decay S → Zγ
The decay S → Zγ exhibits a continuum spectrum plus a gamma line with a finite width, for which at the S rest frame, it has a central energy, -10 -

JHEP07(2020)148
depending on the mass of Z. The gamma line spectrum for this channel S → Zγ at the S rest frame can be expressed in terms of the decay width, Using the narrow width approximation for the resonance Z, the differential width can be written as which is a normalization factor corresponding to In the limit δ → m Z /Γ Z and Γ Z → 0, one has B = π. If taking δ = 2m Z /m S , our narrow width approximation is consistent with that used in refs. [57,58]. Numerically, we will use δ = 3. The result is insensitive to the value of δ 2, especially when δ 3. Changing the variable from M 2 to E , we obtain For S → Zγ, the continuum spectrum results from the cascade decay of Z. The energy of Z emitted from S → Zγ in the S rest frame are Thus, the continuum spectrum in the S rest frame can be written as In eq. (3.29), the factor of "1/2" is due to the fact that the PPPC4DMID spectrum is given by two Z bosons.
In summary, the photon spectrum of S → Zγ in the S rest frame is given by  [52,53], which was generated from Pythia 8.135 [54].
where x ≡ 2E /m S . As shown in eq. (3.20), we will further consider the energy resolution of the instrument for the gamma line signal by convolving the spectrum with a Gaussian kernel.

Gamma-ray spectrum for DM DM
In figure 3, using the above formulas, we show the gamma-ray spectra for DM DM → SS (blue curve) with m S = m h , in comparison with the case of DM DM → hh (red curve) obtained by the PPPC4DMID package [52,53], which was generated from Pythia 8.135 [54]. Physically, in the limit m S = m h , the produced spectrum, independent of the mixing angle α, should be the same for these two annihilation modes. For the case generating energetic S particles, our result is in good agreement with PPPC4DMID, while for the case of the final states S having a low kinetic energy, our result has a better resolution for the spectrum at energies about the S → γγ production (see eq. (3.18)).

Fits to the Galactic center excess spectrum
In order to satisfy the purpose of having a good fit to the Fermi GC excess spectrum and to show the spectral line structure, we take into account three cases: 18 GeV), for which the first two cases can figure out the boost dependence of the observed spectral line width due to the small mass difference of m X and m S , and the third case is expected to be consistent with the WIMP case dominated by XX → hh as it should be. We can use the third case to evaluate the validity of our calculation. We fit the DM mass m X and low-velocity annihilation cross section σv to the Fermi GC gamma-ray excess spectrum carried out by CCW [25]. CCW result covers the photon energy range 300 MeV-500 GeV, within ROI extended to | | ≤ 20 • and 2 • ≤ |b| ≤ 20 • . We JHEP07(2020)148 perform a χ 2 fit, given by [25] where the covariance matrix Σ contains statistical error, correlated empirical model systematics and correlated residual systematics, for which the latter two are non-diagonal. Here dΦ γ /dE i and dΦ obs γ /dE i respectively denote the model prediction and (CCW) central value of the observed flux in the ith energy bin with i ∈ [1,24] in the energy range.
For the gNFW halo profile, we use the scale radius r s = 20 kpc, r = 8.5 kpc, γ = 1.2 and ρ = 0.4 GeV/cm 3 as canonical inputs in the analysis. Because the CCW analysis was performed on the Fermi Pass 7 data, of which the energy resolution is about 10%, 5 [39], we thus use ξ = 0.1 for the line spectra, which are generated from S → γγ and Zγ, in the numerical fit.
In figure 4, Comparing with the GC gamma-ray excess data obtained by CCW, we show the spectrum in the middle panel of figure 4 using the best fit values of m X and σv . For the illustrative purposes, in the right panel, by multiplying the best fit values of m X and σv by 1.4, we draw the theoretical spectrum (the blue curve), where from up to down the p-values are respectively 0.12, 0.15, and 2 × 10 −6 , for which the last one is poor in the fit.
Using the canonical parameter set, for a nearly degenerate case with m X m S , we show the parameter space that provides a good-fit result (p-value ≥ 0.05) located in the range of m X ∈ [60, 132] GeV and σv ∈ [2.0, 6.8] × 10 −26 cm 3 /s, corresponding to the energy of the gamma-ray line ∈ [30,66] GeV. In this Higgs portal scenario, the gamma-ray line signal originating from S → Zγ is highly suppressed compared to the continuum signal, because of the smallness of its branching ratio for the value of m S preferred by the GC excess data 5 The energy resolution of Pass 8 (P8R3 SOURCE V2) has been improved to be about 6%-8.5% from 10 GeV to 200 GeV; see "http://www.slac.stanford.edu/exp/glast/groups/canda/lat Performance.htm".  (see figures 2 and 4). As shown in figure 4, the observed spectral line width is very sensitive to the boost of S. For further comparison with the effect due to the energy resolution of the instrument, we consider a higher energy resolution of ξ =2%, which would be achievable in the DAMPE [45] and GAMMA-400 [46,47] experiments, and show the results in figure 5.
One can thus expect that for the case m S 0.99 m X , the prominent gamma-ray line signal, generated from S → γγ, is distinguishable against the continuum spectrum. Once Fermi-   LAT can accumulate as much as 15 years of data [44], the gamma-ray signal with energy 40 GeV (corresponding to a larger Br(S → γγ)) predicted in the Higgs portal scenario is very likely to be directly examined in the near future.
We find that the best fit is m X m S 86 GeV, featuring a p-value of 0.42 (see table 1). In other words, the corresponding gamma-ray line peaks at 43 GeV. It is very interesting to note that Liang et al. recently found a line-like structure at ∼ 43 GeV with the significance ∼ 3.0σ after analyzing 85 month Pass 8 Fermi-LAT data (P8R2 ULTRACLEAN V6) in the directions of 16 Galaxy clusters which are expected to have large J factors [59]. Further extensive analyses of this line signal should be crucial for testing this scenario and identifying the nature of dark matter.

Constraints from other measurements
Here we present the parameter constraints from various measurements. For the analysis shown in figures 6 and 7, we adopt the fiducial value ρ = 0.4 GeV/cm 3 , and fix γ to be 1. of α, can be further constrained by the direct detection. As shown in the right panel of figure 6, a smaller mass difference of X and S requires a larger g dm to account for the GC data due to the fact that the phase space for XX → SS vanishes in the limit m S → m X . The detailed constraints from various measurements will be discussed as follows.
Fermi-LAT gamma-ray line search. The Fermi-LAT collaboration has recently placed constraints on the gamma-line signals [41,60]. The resultant limit depends on the mass and density profile of DM. We consider the Fermi R3 (ROI) fit from which the constraint, compared with other ROI results, is more restrictive. The R3 is defined to a very small circular regions of radius R GC = 3 • centered on the GC, and optimized for the contracted NFW 6 (NFWc) profile with γ = 1.3. Since this choice of ROI strongly depends on the value of γ (see the discussion in appendix B of ref. [60]), and since Fermi R3 line limit and CCW data correspond to different ROIs, we do not rescale the inner slope of the halo profile of the former one to γ = 1.2 to match each other. We remark that if the Fermi line data did not depend on γ, such a rescaling would weaken the constraint on the annihilation cross section by a factor of two. In figures 6 and 7, the 95% confidence level (C.L.) bound from the updated Fermi-LAT R3 (NFWc) gamma-ray line search (5.8 years of Pass 8 data) [41] is depicted by the solid magenta line corresponding to the use of ρ = 0.4 GeV/cm 3 . One should note that the Fermi gave the gamma-ray line limit for direct DM annihilation to the photon pair, while in our case four photons are produced per annihilation. Therefore, in our case the measured gamma-ray line energy is m X /2, and the limit for σv × Br(S → γγ) is equivalent to the value of 2 σv γγ given in ref. [41].
Fermi-LAT observations of dwarf spheroidal galaxies. We perform a combined likelihood analysis of 28 kinetically confirmed and 17 candidate dSphs with 6 years of the Fermi-LAT data 7 (passing the P8R2 SOURCE event class selections), where gamma-ray JHEP07(2020)148 energies are in the range from 500 MeV to 500 GeV [40]. As in ref. [40], we use the spectroscopically determined J-factors with errors for the confirmed dSphs, and adopt predicted values from the distance scaling relationship with a nominal uncertainty of 0.6 dex for the newly discovered candidates. We refer readers to ref. [16] for the detailed description of the likelihood analysis that we have used here.
The solid and long-dashed red lines shown in figure 6 represent the current dSph limit at 95% C.L. and Fermi-LAT projected sensitivity, respectively. Here we have assumed 60 dSphs (≡ N dSph ) observed and the 15 years of data (≡ N data ) collected for the projected sensitivity, which approximately rescales with N dSph × √ N data [61]. As shown in figures 6 and 7, the parameter space is much more restricted by the current dSphs constraint, compared with other measurements. The scenario that the hidden sector dark matter interacts with the SM through the Higgs portal can be further tested by the dSphs projection.
Planck cosmic microwave background. The CMB provides a probe into the DM annihilation at the epoch of recombination, and thus offers a complementary constraint compared with experiments of the gamma-ray observations. Planck sets a bound on the annihilation parameter, p ann , from TT, TE, EE+lowP (temperature and polarization) data combinations [2], where σv CMB at the epoch of recombination σv at the present day for s-wave DM annihilation (as the secluded vector dark model that we study in this paper), and the efficiency factor f eff is the fraction of the energy that is injected into the intergalactic medium from DM annihilations at redshift z. The efficiency factor depends on the spectra of e ± pairs and photons produced following DM annihilations, where we use f γ,e − eff (E) curve results suited for the "3 keV" baseline prescription obtained by Slatyer [62], and (dN e − /γ /dE) X is the electron/photon spectrum produced per DM annihilation in the DM rest frame. The calculation for (dN e − /dE) X , which originates from two-body and three-body S decays following the DM annihilation XX → SS, is completely the same as that for (dN γ /dE) X described in section 3.2, but using PPPC4DMID to generate the electron spectrum instead of the photon spectrum.
The Planck CMB 95% C.L. bound is sketched as the dot-dashed brown line in figures 6 and 7. The current Planck CMB limit seems to be considerably weaker than the Fermi-LAT dSphs limit.
Correct relic density. The thermodynamic evolution of the hidden sector (vector) dark matter interacting with the SM through the Higgs portal has been studied in ref. [17]. Here we will present the main properties, and refer readers to ref. [17] for the detailed results. The value of α is relevant to the coupling strengths of S to the SM particles, and thus determine -18 -JHEP07(2020)148 the decoupling temperature below which the hidden sector is kinetically decoupled from the SM bath.
For α 2 × 10 −6 , the correct relic density is set by the XX ↔ SS interaction, so that the DM particles can be in chemical and thermal equilibrium with S particles and with the SM bath (through S) before freeze-out. While this result is consistent with the conventional WIMP scenario, the annihilation cross section corresponding to the narrow gray range in figures 6 and 7 can account for the correct relic density. As for α less than 2 × 10 −6 , the dark sector has been kinetically decoupled from the thermal reservoir, before it becomes nonrelativistic. For this case, the DM annihilation cross section and coupling contant g dm , providing a correct relic density, could be boosted to the upper side of the gray range in the left panel and right panel of figures 6 and 7, respectively.
XENON1T result and LUX-ZEPLIN (LZ) projected sensitivity. Considering a specific DM model, one can set limits on (coupling) parameters from the direct detection experiments. For the secluded vector DM model shown in section 2, the elastic scattering cross section of X off a nucleon (N ), independent of the nuclear spin, is referred to as the "spin-independent cross section", which via the t-channel interactions with exchange of S and h is given by is the reduced mass of X and N , and f N = q N |qq|N · m q /m N 0.3 [63,64]. In the right panel of figure 6, using ρ = 0.4 GeV/cm 3 , the 95% C.L. bounds from XENON1T [42] and LUX-ZEPLIN (LZ) projected sensitivity [43] are respectively indicated by the gray and green solid contours, where the corresponding α values are denoted, and, in the region m X < m h (or m X > m h ), the r.h.s. (or left hand side (≡ l.h.s.)) of each line is allowed. We observe that XENON1T and LZ projected sensitivity are able to constrain the m X − g dm parameter space preferred by the GC gamma-ray excess data for α 0.01 and 0.002, respectively.
In figure 7, we do not consider direct detection limits on the case of m S = m h , for which there is no constraint for this perfect degenerate case. However, for a generic case of |m S − m h | < 4 GeV and |α| 0.17 (0.02), the region g dm < 1 evades the XENON1T searches (LZ projected sensitivity) for ρ ∈ [0.25, 0.6] GeV/cm 3 .
Big bang nucleosynthesis (BBN). The BBN measurement can set a lower bound on the mixing angle, α. When the scalar mediator S decays out-of-equilibrium into SM particles, the Universe becomes radiation-dominated again and experiences the reheating owing to large entropy injection. (See the discussions in ref. [17].) We can constrain α through the observation limit on the reheating temperature. If the reheating temperature T RH was on the order of the neutrino decoupling temperature, then the neutrinos would not be well thermalized [65]; if so, the relative rate of light element abundances would be changed, too. From the Y p + D/H analysis (with the helium nucleon fraction Y p ≡ -19 -

Mixing angle α constraint in the hidden Higgs portal dark matter model: from the thermodynamic point of view
If α is extremely small, the hidden sector can be decoupled from the SM bath at the very high temperature T m X,S . Like the hot dark matter case, after decoupling, the relativistic hidden sector almost maintains the same temperature as the bath, and its comoving number density is also approximately unchanged. Here we have neglected temperature variation of the SM bath due to the change of its relativistic degrees of freedom, when its temperature drops below m t or m h .
Therefore, for the case that the hidden sector was in thermal equilibrium with the bath in the earlier stage but later on was decoupled from the bath even at the very high temperature, once the decoupling has occurred, the relativistic hidden sector evolves with a temperature which is almost the same as the temperature of the bath. Moreover, for this case, when the hidden sector becomes nonrelativisitc, it gets hotter than the bath during the cannibal epoch [17]. In ref. [17], we have presented a comprehensive study on thermodynamic evolution of the hidden sector for this secluded vector DM model. Here we would like to discuss the minimum value of α for which the hidden sector was once in thermal equilibrium with the bath when T m S . This part did not mention in ref. [17].
In the following discussion, we assume that, before decoupling, number changing interactions among the dark sector particles are still active and guarantee their thermal equilibrium with zero chemical potential. We separately discuss the requirement of interactions, including (i) SS ↔ SM SM, (ii) S ↔ SM SM, and (iii) S SM ↔ S SM, that can account for the thermal equilibrium between S and the SM bath at a temperature T which is larger than m S .
If the hidden sector is in equilibrium with the thermal bath through the interaction SS ↔ SM SM at a temperature T m S , we need to have n eq S σv SS→SM SM H, which describes the S production rate from the inverse annihilation is larger than the cosmic dilution rate. Here n eq S is the equilibrium number density (with zero chemical potential) and H is the Hubble rate. In the limits of large energy and small α, because σv SS→SM SM ∝ α 2 /s 2 which is suppressed in the high temperature due to the fact that √ s ∝ T , therefore we can simply take T ≈ m S to obtain the lower bound of the mixing angle, α 10 −4 for this interaction. (See appendix B of ref. [17] for the exact form of the SS → SM SM amplitude.) If thermal equilibrium between S and the SM bath is due to S ↔ SM SM and hold at T m S , we have the inverse decay rate Γ S H. Since H ∝ T 2 , we can simply take T ≈ m S to get α 10 −5 (see also the result shown in the right panel of figure 1 in ref. [17]).

JHEP07(2020)148
As for the elastic scattering S SM ↔ S SM, we adopt the definition of temperature for the relativistic S, where the distribution function f S exp[−(E S − µ S )/T S ] with µ S the chemical potential, E S is the distribution function, g S = 1 is the internal degrees of freedom, n S is the number density, and p S are energy and 3-momentum of S, respectively. See further discussions for the definition of temperature in appendix D. Eq. (5.1) is a good approximation for the temperature definition even at the very high temperature, T S m S . Solving the Boltzmann moment equation, we obtain the temperature evolution of S for T S m S , where the elastic collision term is described by a semi-relativistic Fokker-Planck equation [66], Here the momentum relaxation rate is given by where f is the relevant relativistic SM species, and f f (T ) is its distribution function at temperature T , k is the 3-momentum of f , t is the momentum transfer squared between S and f . Note that this formula is a good approximation for a relativistic S under the condition E S √ t. Taking the limits E S m S and T (S) m S , we then obtain Note that in ref. [17], we consider the case satisfying T m S . If the elastic scattering S SM ↔ S SM can maintain S and the SM bath in thermal equilibrium at T m S , the -21 -JHEP07(2020)148 energy gained by S through the elastic scattering is larger in magnitude than the Hubble cooling rate,γ Because H ∝ T 2 andγ ∝ α 2 T 5 (from eq. (5.6)), we thus have α ∝ T −3/2 el with T el being the decoupling temperature for the elastic scattering interaction. A smaller mixing angle α will result in a higher T el . Using the result shown in figure 5 of ref. [17], from that we have α ∼ 10 −6 when T el ∼ m S , therefore for the case that the hidden sector is kinetically decoupled from the bath at T el ( m S ), the corresponding mixing angle reads As will be discussed below, the vacuum will become unstable when T 10 10 GeV for the secluded vector DM model. If T el is below this value, we shall need α 10 −18 .

Theoretical vacuum stability for the secluded vector dark matter model
Before concluding this paper, we study the scale-dependence of vacuum stability for the secluded vector dark matter model. The vacuum is required to be stable at the tree-level potential, i.e., the potential should be bounded from below and satisfies, Meanwhile, requiring det M 2 Higgs > 0, we also have λ 2 HS − λ H λ S < 0. From this and eq. (5.9) we thus get −|λ HS | + √ λ H λ S > 0. It was pointed out in ref. [67], where one-loop β-functions were considered, that the top quark can drive λ H to become negative at a certain higher scale, such that the electroweak vacuum is no longer the global minimum. To examine the vacuum stability, we study the renormalization group equations (RGEs) of the quartic scalar couplings, for which using SARAH [68][69][70][71] the β functions are calculated up to the two-loop level, and collected in appendix E. We find that, in the limit of m S → m X and α → 0, the scale-dependence of the quartic scalar couplings, highly insensitive to the values of m X and m S , depends only on the initial value of g dm . As seen in figure 8, for the m X ≈ m S case with a sizable mixing angle α, the stability condition is violated, i.e., λ H becomes negative, at the scale Q less than 10 10 GeV, but the violating scale will approach to 10 10 GeV in vanishing α limit. In other words, the present secluded vector DM model is an effective theory suitable for the scale below 10 10 GeV. Above this temperature, the Universe might experience the reheating and could be dominated by more massive particles during the oscillation epoch.

Conclusions
The gamma-ray line signal generated from the DM annihilation could be a clear signature which is distinguishable from astrophysical backgrounds and reveals the particle nature of DM. We are motivated by the recent studies that direct DM annihilation to two SM-like Higgses, produced close to rest, is capable of accounting for the GC gamma-ray excess data but with a little lower p-value 0.13, and also motived by the fact that the quality of fit can be significantly improved if DM mostly annihilates to a lighter Higgs pair, which soften the gamma-ray spectrum to have a better fit to the observation peaked at 1-3 GeV.
We therefore consider a Higgs portal DM model where the hidden scalar mediates the interaction of DM with the SM due to its mixing with the SM Higgs. In this model, the DM is secluded in the hidden sector and can annihilate directly to a pair of lighter scalar mediators, each of which, nearly degenerate with DM in mass, subsequently decays into the SM particles.
For the case of m X m S , we have obtained that the parameter region m X ∈ [60, 132] GeV can provide a good spectral fit to the Fermi GC gamma-ray excess data (pvalue ≥ 0.05), showing the energy of the gamma-ray line ∈ [30,66] GeV. The best fit to the data yields m X m S 86 GeV, featuring a p-value of 0.42, so that the corresponding gamma-ray line arising from S → γγ peaks at 43 GeV. The observed spectral line width, depending on the energy resolution of the detector (see figure 4 vs. 5), is sensitive to the Lorentz-boost from the mediator rest frame to the DM center of mass frame, that directly -23 -

JHEP07(2020)148
correlates with the mass difference of m X and m S . In the Higgs portal model, we expect that, for the secluded DM case 0.99 m S /m X < 1, a prominent gamma-ray line arising from S → γγ can be distinguished from the continuum spectrum, while the line signal originating from S → Zγ is highly suppressed.
The fitted value of the low-velocity DM annihilation cross section depends on the DM distribution. Adopting γ = 1.2 and ρ = 0.4 GeV/cm 3 , a good fit to the GC excess emission gives σv ∈ [2.0, 6.8] × 10 −26 cm 3 /s. Using a smaller (or larger) γ and/or ρ , the value of σv can be further raised (or lowered). We have derived constraints on the annihilation cross section from the Fermi-LAT gamma-ray line search, Fermi-LAT dSphs gamma-ray observations, and Planck CMB measurement. These detections can offer complementary probes. Currently, the dSphs constraint on the parameter space favored by the GC excess emission is more restrictive than that derived from other measurements. Considering the renormalizable secluded vector dark matter model, we have shown the results in the m X − g dm parameter space, where some regions favored by the GC excess emission can be excluded by XENON1T for α 0.01 and further by the LZ projected sensitivity for α 0.002.
The constraint from the Planck BBN measurement requires the S width Γ S (0.03 sec) −1 , from which we can put a lower bound on the mixing angle α 1 × 10 −10 . On the other hand, for a case with a small mixing angle α 2 × 10 −6 , the hidden sector has been kinetically decoupled from the bath before it becomes nonrelativistic [17]. As such, the correct relic density is described by a DM annihilation cross section which could be significantly boosted above the conventionally WIMP value [17]. In this paper, we have discussed the case with an extremely small value of α, for which the relativistic hidden sector can be decoupled from the SM bath at very high temperatures T m X,S , and, after decoupling, almost maintains the same temperature as the bath until T ∼ m S,X . Assuming that the number changing interactions among the dark sector particles can guarantee their thermal equilibrium with zero chemical potential before decoupling, we have shown that the mixing angle α and elastic decoupling temperature T el satisfy the relation: α ≈ 10 −6 (T el /m S ) −3/2 . We have also shown that in this present scenario the vacuum of the secluded vector DM model will become unstable when T 10 −10 . Thus, if T el 10 −10 is required for the secluded vector DM model, we obtain α 10 −18 .
For this Higgs portal (vector) DM model, the dSphs projected sensitivity can further probe the region of the annihilation cross section, which will be able to approach the conventional WIMP value. Furthermore, extensive analyses for the line signal with the energy in the range of 30-66 GeV should be crucially important in testing this Higgs portal scenario and identifying the nature of dark matter in the future.

B The annihilation cross section for XX → SS
As shown in figure 9, the cross section for the XX → SS in the laboratory frame, where one of the incoming particles is at rest with v lab being the relative velocity measured, is given by [17] σv lab = (σv lab ) 4v,s + (σv lab ) t,u + (σv lab ) int , with the center-of-mass energy of √ s 2m X for the low-velocity DM annihilation, and the triple SSS coupling being and v S = m X /g dm . Here (σv lab ) 4v,s is the cross section resulting from the 4-vertex and s-channel diagrams, (σv lab ) t,u is from the t-and u-channels, and (σv lab ) int is from the interference between (4-vertex, s) and (t, u). An interesting property is the DM annihilation amplitudes, Figure 9(a), 9(b), 9(c), 9(d) where i,µ is the polarization vector of the initial DM particle, and (B.7) -27 -

JHEP07(2020)148
In eq. (B.6) we have used Γ S ≈ 0 in the α → 0 limit, and P µν ≈ Q µν ≈ g µν , t = (p 1 −p 3 ) 2 ≈ 0 and u = (p 1 − p 4 ) 2 ≈ 0 in the limits of s → 4m 2 X and m X → m S . Therefore, in this limit, the total amplitude can be well approximated as −4ig 2 dm 1,µ µ 2 . In other words, for a nearly degenerate case of X and S with a small mixing angle α, we have σv lab ≈ (σv lab ) 4v,s , which is numerically confirmed. Note that we have neglected the diagram XX → h * → SS, which, corresponding to figure 9(b) but with the propagator replaced by h, is further suppressed by sin 2 α in the amplitude level, because the coupling of the X − X − h vertex is "is α g dm m X ", while the coupling of the h − S − S vertex is , in the small α limit. (B.9) For most of the GC favored regions in the Higgs portal model, we find that m X < m h (= 125.18 GeV), i.e. XX → hh is kinematically forbidden. Nevertheless, as shown in figure 4, if S and h are degenerate in mass, only a very small GC region, corresponding to p-value 0.09 and m X ∈ [m h , 128 GeV], is allowed; in this region, the DM annihilation is still dominated by XX → SS, while the XX → hh amplitude, for which the 4-vertex, u-, and t-channels ∝ sin 2 α and s-channel ∝ sin α, is relatively suppressed by sin α.
C The dependence of the allowed parameter space on the value of ρ The uncertainties of the dark matter distribution near the Galactic center and its local density are still large. CCW have analyzed the GC inner slope for 60 Galactic diffuse emission models, and found γ = 1.2 ± 0.1 preferred within a ROI: | | ≤ 20 • and 2 • ≤ |b| ≤ 20 • . To further illustrate the dependence on variation of ρ ∈ [0.25, 0.6] GeV/cm 3 for the GC excess favored region compared with other constraints, in figures 10 and 11, we employ three values of ρ = 0.25, 0.4 and 0.6 GeV/cm 3 . On the other hand, if using a smaller (or larger) γ, we can naïvely expect from the change of the J-factor that the value of σv is further raised (or lowered) for the GC favored region. This appendix is a complement to section 4.2.
The bounds of the Femi gamma-ray line search [41] and direct detections [42,43] also depend on the value of ρ , which are depicted in figures 10 and 11 by the solid, dotted, and dashed magenta lines, corresponding to the use of ρ = 0.4, 0.25, and 0.6 GeV/cm 3 , respectively. For the right panel of figure 11, the direct detection does not set the bound on this perfect degenerate case. Note again that for ρ ∈ [0.25, 0.6] GeV/cm 3 , the parameter space, where g dm < 1, |m S − m h | < 4 GeV and |α| 0.17 (0.02), can evade the bound from the XENON1T measurement (LZ projected sensitivity).

D The definition of temperature for S
In the high temperature limit E S m S , the temperature of S satisfies the relation, 2.

16.
m X (GeV) 〈σ v〉  where f S = [exp((E − µ S )/T S ) − 1] −1 with µ S the chemical potential, g S = 1 is the internal degrees of freedom, and n S is the number density. Assuming that S is in chemical equilibrium with zero chemical potential, we can approximate the r.h.s. of eq. (D.1) as (3) 6ζ (3) T S (0.90 + 0.10) , (D. 2) which shows that the term ∝ f 2 S on the r.h.s. of eq. (D.1) gives about 10% correction in amount. On the other hand, in high temperatures, we can expect that the average number of particles in each state of the phase space is much less than 1, i.e., 1 + f S 1, and thus approximate the S distribution as Using the approximate distribution, we obtain g S n S (T S ) which is a good approximation for the temperature of S. Here the approximation of the thermal equilibrium number density (with µ S = 0) is less than the exact value by a factor of 17%. We thus use eq. (D.4) as the benchmark to derive the Boltzmann moment equation of the hidden scalar temperature, which is suitable at high temperatures ( m S ).