A Possible Explanation of Low Energy $\gamma$-ray Excess from Galactic Centre and Fermi Bubble by a Dark Matter Model with Two Real Scalars

We promote the idea of multi-component Dark Matter (DM) to explain results from both direct and indirect detection experiments. In these models as contribution of each DM candidate to relic abundance is summed up to meet WMAP/Planck measurements of $\Omega_{\rm DM}$, these candidates have larger annihilation cross-sections compared to the single-component DM models. This results in larger $\gamma$-ray flux in indirect detection experiments of DM. We illustrate this fact by introducing an extra scalar to the popular single real scalar DM model. We also present detailed calculations for the vacuum stability bounds, perturbative unitarity and triviality constraints on this model. As direct detection experimental results still show some conflict, we kept our options open, discussing different scenarios with different DM mass zones. In the framework of our model we make an interesting observation: The existing direct detection experiments like CDMS II, CoGeNT, CRESST II, XENON 100 or LUX together with the observation of excess low energy $\gamma$-ray from Galactic Centre and Fermi Bubble by FGST already have the capability to distinguish between different DM halo profiles.


Introduction
The overwhelming cosmological and astrophysical evidences have now established the existence of an unknown non-luminous matter present in the universe in enormous amount, namely the dark matter (DM). Experiments like Wilkinson Microwave Anisotropy Probe (WMAP) [1], BOSS [2] or more recently Planck [3] measure the baryonic fraction precisely to consolidate the fact that this non-baryonic DM constitutes around ∼ 26.5% of the content of the universe. The particle nature of DM candidate is still unknown. The relic density of dark matter deduced from cosmological observations mentioned above tends to suggest that most of the DM could be made of weakly interacting massive particles (WIMPs) [4,5,6,7] and they are nonrelativistic or cold in nature. This calls for an extension of the standard model (SM) of particle physics. Many such extensions have been suggested in the literature in the framework of supersymmetry, extra dimensions, axion etc. Models such as Kaluza Klein [8], inert triplet [9] or supersymmetry breaking models like mAMSB [10] predict very massive DM whereas models like SMSSM [11], axion [12] predict DM of lower mass. Phenomenology of simpler extensions of SM like fermionic DM model [13] or inert doublet model [14] has been elaborately studied. Amongst all such options, extending the scalar sector is particularly interesting because of its simplicity.
The minimal extension with a single gauge singlet real scalar stabilised by a Z 2 symmetry in the context of dark matter was proposed by Silveira and Zee in Ref. [15] and then it was extensively studied in the literature [16] - [37]. In Ref. [38] the singlet scalar DM model has been discussed with a global U(1) symmetry.
Amongst the non-minimal extensions, a DM model where SM is extended by a complex singlet scalar has been considered in Refs. [39,40,41]. A DM model with two real scalars has been discussed in Refs. [42,43], where one scalar is protected by a Z 2 symmetry, but the Z 2 symmetry protecting the other one spontaneously breaks. In all these non-minimally extended models there is, however, only one DM candidate.
A pertinent question to ask at this stage is that, if our visible sector is enriched with so many particles, why the DM should be composed of only one component? We therefore intend to discuss in this paper a model with two DM candidates. In some earlier works [44] the idea of multicomponent dark matter has been discussed in details. We extend the SM with two gauge singlet real scalars protected by unbroken Z 2 symmetries. As mentioned in the abstract, the advantage of such a multi-component DM model is that the DM annihilation can be enhanced so that one can expect spectacular signals in the indirect detection experiments. Because of this simple fact, the thermal averaged annihilation cross-sections in this model can enjoy enhancement upto two orders of magnitude compared to the models with one real scalar. As a result, in the present model with two real scalars we can make an attempt to explain both direct and indirect detection DM experimental observations, which was not possible with the model with a single real scalar.
Direct detection DM experiments can detect DM by measuring the recoil energy of a target nucleon of detecting material in case a DM particle happens to scatter off such nucleons. Experiments like CDMS [45,46], DAMA [47], CoGeNT [48] or CRESST [49] present their results indicating allowed zones in the scattering crosssection -DM mass plane. These experiments seem to prefer low dark matter masses ∼ 10 GeV. Some earlier works on ∼ 10 GeV DM mass have been done [50,51]. XENON 100 [52,53], however, did not observe any potential DM event contradicting claims of the earlier experiments and has presented an upper bound on DM-nucleon scattering cross-section for various DM masses. Recent findings by LUX [54] have fortified claims by XENON 100 collaboration.
The indirect detection of DM involves detecting the particles and their subsequent decay products, produced due to DM annihilations. Huge concentration of DM are expected at the centre of gravitating bodies such as the Sun or the galactic centre (GC) as they can capture DM particles over time.
The region in and around the GC are looked for detecting the dark matter annihilation products such as γ, ν etc. Fermi Gamma Space Telescope (FGST), operated from mid of '08, have been looking for the gamma ray from the GC [55]. The low energy gamma ray from GC shows some bumpy structures around a few GeV which cannot be properly explained by known astrophysics. A plausible explanation of such a non-power law spectrum is provided by DM annihilations [56,57].
The emission of gamma rays from Fermi Bubble may also be partially caused by DM annihilations. The Fermi Bubble is a lobular structure of gamma ray emission zone both upward and downward from the galactic plane has been discovered recently by Fermi's Large Area Telescope [58]. The lobes spread up to a few kpc above and below the galactic plane and emit gamma ray with energy extending from a few GeV to about a hundred GeV. The gamma emission is supposed to be produced from the inverse Compton scattering (ICS) of cosmic ray electrons. But more involved study of this emission reveals that while the spectra from the high galactic latitude region can be explained by ICS taking into consideration cosmic electron distribution, it cannot satisfactorily explain the emission from the lower latitudes. The γ-ray flux from possible DM annihilation in the galactic halo may help explain this apparent anomaly [59,60,61].
As mentioned earlier, in this work, the proposed model with two gauge singlet real scalars protected by a Z 2 ×Z 2 symmetry is confronted with the experimental findings of both direct and indirect DM experiments. As direct DM experimental results are contradictory, we keep an open mind in considering them. We find constraints on the parameter space from these results and then try to constrain them further imposing the requirement of producing the relic density of the DM candidates consistent with Planck observations. We then choose benchmark points from this constrained parameter space to explain results from indirect detection experiments from DM annihilation taking into consideration different DM halo profiles.
The paper is organised as follows. In Sec. 2 we discuss the theoretical framework of our proposed model. The theoretical constraints from vacuum stability, perturbative unitarity, triviality and experimental constraints from the invisible branching ratio of the Higgs boson have been discussed in Sec. 3. The next section contains the relevant relic density calculations. The model is confronted with direct detection experiments and Planck observations in Sec. 5. Explanation of the observed excess of γ-ray from GC and Fermi bubble by our model is studied in Sec. 6. We conclude in Sec. 7.

Theoretical Framework
We propose a model where two real scalar singlets (S and S ) are added to the standard model.
The general form of the renormalisable scalar potential is then given by, where H is the ordinary (SM) Higgs doublet. In the above δ's denote the couplings between the singlets and the Higgs and k's are the couplings between these singlets themselves.
The stability of DM particles is achieved by imposing a discreet symmetry Z 2 onto the Lagrangian. Depending on whether S and S are odd under the same Z 2 or not, we discuss two scenarios for completeness.

Lagrangian Invariant under
If only S and S are odd under the same Z 2 , and the rest of the particles are even, some parameters of the potential vanish: so that the scalar potential (1) reduces to the following, After the spontaneous symmetry breaking the mass matrix for S and S is given by v √ 2 denotes the vacuum expectation value of the Higgs. After diagonalisation the masses of the physical eigenstates S 1 and S 2 are given by where tan 2θ = 2M 12 2.2 Lagrangian Invariant under Z 2 × Z 2 If S and S are stabilised by different discrete symmetries, so that the scalar potential (4) further reduces to After spontaneous symmetry breaking the respective masses of S and S are given by The four beyond SM parameters determining the masses of the scalars are k 2 , k 2 , δ 2 and δ 2 .
In both Z 2 × Z 2 and Z 2 × Z 2 cases if SS ↔ S S scattering processes can be avoided, the model can give rise to a two-component DM scenario. However, as the later case has fewer number of beyond SM parameters, in the following we will restrict ourselves only to the Z 2 × Z 2 invariant Lagrangian.

Constraints on Model Parameters
The extra scalars present in the model modifies the scalar potential. Hence it is prudent to revisit constraints emanating from vacuum stability conditions and triviality of the Higgs potential. Perturbative unitarity can also get affected by these scalars. Limits on the invisible decay width of Higgs from LHC severely restricts such models. In the following we elaborate on these constraints.

Vacuum Stability Conditions
Calculating the exact vacuum stability conditions for any model is generally difficult. However, for many dark matter models the quartic part of the scalar potential can be expressed as quadratic form (λ ab S 2 a S 2 b ) with the squares of real fields as single entity. Lagrangian respecting Z 2 symmetry which ensures the stability of scalar dark matter have the terms which can be expressed like that. The scalar potential of our proposed model can also be expressed in a similar form as above because of preservation of Z 2 symmetry. The criteria for copositivity allow one to derive properly the analytic vacuum stability conditions for the such matrix λ ab from which sufficient conditions for vacuum stability can be obtained. 4 The necessary conditions for a symmetric matrix A of order 3 to be copositive are given by [62,63,64,65], a 11 0, a 22 0, a 33 0, a 12 = a 12 + √ a 11 a 22 0, a 13 = a 13 + √ a 11 a 33 0, and √ a 11 a 22 a 33 + a 12 √ a 33 + a 13 √ a 22 + a 23 √ a 11 + √ 2ā 12ā13ā23 0.
The last criterion given in Eq. (13) is a simplified form of the two conditions (Eqs. (14) and (15)) below √ a 11 a 22 a 33 + a 12 √ a 33 + a 13 √ a 22 + a 23 √ a 11 0, det A = a 11 a 22 a 33 − (a 2 12 a 33 + a 2 13 a 22 + a 11 a 2 23 ) + 2a 12 a 13 a 23 0, where one or the other inequality has to be satisfied [63] 5 . The conditions Eq. (12) impose that the three 2 × 2 principal submatrices of A are copositive. The matrix of quartic couplings Λ in the (h 2 , S 2 , S 2 ) basis for the potential Eq. (9) is given by Copositivity criteria of Eqs. (12) and (13) yield the necessary and sufficient vacuum stability conditions, λ 0, k 4 0, k 4 0, and The conditions of Eqs. (17) and (18) simply determine the vacuum stability bounds on our model. We restrict the parameter space by these conditions for later calculation.

Perturbative Unitarity Bounds
The potential of the Z 2 × Z 2 model is bounded from below if Eq. (17) and Eq. (18) are simultaneously satisfied. Then, λ, k 4 > 0 and δ 2 > 0 or The Higgs mechanism generates a mass of M 2 H = 1 2 λv 2 for the Higgs and also contributes to the mass of the S particle For H = (0, v/ √ 2) and S = 0, S = 0 be a local minimum we should have M 2 S > 0 and M 2 S > 0 are required. This is also a global minimum as long as 17]. The potential of the scalar sector after electroweak symmetry breaking in the unitary gauge can be written as, After that tree-level perturbative unitarity [66] to scalar elastic scattering processes have been applied in this model (Eq. (22)). The zeroth partial wave amplitude, must satisfy the condition |Re(a 0 )| ≤ 1 2 [67]. In the above, s is the centre of mass (CM) energy, p CM i,f are the initial and final momenta in CM system and T 2→2 denotes the matrix element for 2 → 2 processes with θ being the incident angle between two incoming particles.
The possible two particle states are HH, HS, SS, S S , S S, HS and the scattering processes include many possible diagrams such as HH → HH, SS → SS, HS → HS, HH → SS, SS → HH, S S → HH, S S → SS, HS → HS , S S → S S, HH → S S , SS → S S , S S → S S . The matrix elements (T 2→2 ) for the above 2 → 2 processes are calculated from the tree level Feynman diagrams for corresponding scattering and given by, Now using Eq. (23), we have calculated the partial wave amplitude for each of the scattering processes and the coupled amplitude can be written as a matrix form, Requiring |Re(a 0 )| ≤ 1 2 for each individual process above we obtain for HS → HS and HH → SS for S S → SS and S S → S S |k a 4 | ≤ 8π , for HS → HS and HH → S S |δ 2 | ≤ 8π .

Triviality Bound
The requirement for 'triviality bound' on any model is guaranteed by one of the conditions that the renormalization group evolution should not push the quartic coupling constant of such models (say, λ) to infinite value up to the ultraviolet cutoff scale Λ of the model. This requires that Landau pole of the Higgs boson should be in higher scale than Λ. Therefore to check the triviality in our model, namely the two scalar singlet model with Z 2 × Z 2 symmetry, we have to solve the renormalization group (RG) evolution equations for all the running parameters of this model. We have chosen only one-loop contribution in determining the beta functions for our model. The RG equations for the couplings in the model, namely, λ, k 2 , k 2 , k 4 , k 4 , k a 4 , δ 2 , δ 2 are thus obtained at one-loop level as where t = log(µ/M ). In Eq. (40) − Eq. (47) µ denotes the renormalization scale and M is an arbitrary scale. Here γ h = −(9/4)g 2 2 − (3/4)g 2 1 + 3y 2 t and g 1 , g 2 , y t are U(1) Y , SU(2) L gauge couplings and top Yukawa coupling, respectively. In our calculation, the RG equations for gauge and top Yukawa couplings are also taken into account. We have taken the initial condition, We have solved all the RG equations given above and checked the consistency of all the quartic couplings within the suitably chosen scale of the theory. For the initialisation, we have taken λ init corresponding to the recent value of Higgs mass, 126 GeV. The other initial values of parameters in our model should have been chosen from the allowed region of parameter space. In Fig. 1 the variation of different parameters (k 4 , k 4 , k a 4 , k 2 , k 2 , δ 2 , δ 2 , λ) with scale used in this model is shown. The benchmark point 4A of Table 4 has been chosen for assigning initial values in the evaluation of the running of various couplings. The variation of mass of each scalar S (or S ) with energy scale can be obtained from the plots as it is determined by couplings k 2 and δ 2 (or k 2 and δ 2 ). Although we have solved the RG equation for δ 2 and δ 2 , the influence of δ 2 and δ 2 on λ is very small as we can see from Eq. (42) that the RG equation of λ is deviated from the SM RG equation only by the almost smooth term (2δ 2 2 + 2δ 2 2 ). But the allowed region for 'triviality bound' for a given Higgs mass shrinks as the term, (2δ 2 2 + 2δ 2 2 ) starts growing.

Constraints from Invisible Higgs Decay Width
If kinematically allowed, Higgs boson can decay to SS or S S . Such invisible decay channels are severely restricted by the present data from Large Hadron Collider (LHC). The branching fraction is bounded at 95% CL to be less than 19% by the global fits to the Higgs data keeping Higgs to fermion couplings fixed to their SM values. If such a restriction is lifted and additional particles are allowed in the loops the bound get relaxed to B(H → inv) < 38% [69]. Γ SM denotes the SM Higgs decay width and Γ inv is the invisible Higgs decay width, which in our model is given by [36], The benchmark point 4A in Table 4, consistent with the XENON 100 direct detection results, gives B(H → inv) ∼ 0.26 which at present is allowed at 95% CL [70,71,72,73,74]. However as we intend to interpret the low mass regions of dark matter claimed to be probed by several other dark matter direct search experiments (CDMS II, CRESST, CoGeNT etc.) along with indirect searches (low energy γ-ray from Fermi Bubble and Galactic Centre), in some cases M S , M S ≤ M H /2 and the Higgs boson decays invisibly to SS or S S , with a B(H → inv) disfavoured by the LHC observations. This is a well known problem with all such models, where DM annihilation is mediated by the SM Higgs. It can be circumvented by adding extra degrees of freedom lighter than DM particles [75]. In that case it is possible to delineate the dominant Higgs-mediated annihilation channels from the diagrams contributing to the invisible Higgs decay. Such an elaborate model building issue is out of the scope of the present work as we feel given the non-observation of low mass DM events by XENON 100 and LUX, it is still premature to rely on CDMS II and other experiments betting on the low mass DM. For quantitative estimations we have chosen the Higgs to be the 126 GeV SM Higgs as a benchmark.

Calculation of Relic Abundance
In order to calculate the relic abundance for the dark matter candidates in the present formalism, we need to solve the relevant Boltzmann equations.
In presence of one singlet scalar S, the Boltzmann eqn. is given by where n S and n Seq are the number density and equilibrium number density of the singlet scalar, S. σv is the thermal averaged annihilation cross-section of dark matter annihilating to SM particles. Defining dimensionless quantities Y = n e and x = M S T , where e is the total entropy density, Eq. (50) can be written in the form, where g * is the degrees of freedom. The relic density Y 0 (value of Y at T = T 0 ) is obtained by integrating Eq. (51) from initial x value, where T 0 and T f are the present photon temperature (2.726 o K) and freeze-out temperature respectively. The relic density of a dark matter candidate, S, in the units of critical energy density, ρ cr = 3H 2 /8πG, can be expressed as where e 0 is the present entropy density evaluated at T 0 . It follows that knowing Y 0 , we can compute the relic density of the dark matter candidate from the relation [76], In Eq. (53) h = 100 km sec −1 Mpc −1 is the Hubble constant. The Planck survey provides the constraints on the dark matter density Ω DM h 2 from precision measurements of anisotropy of cosmic microwave background radiation as consistent with the previous WMAP measurement 0.1093 < Ω DM h 2 < 0.1183. In our model with two real scalars, both S and S contribute to the relic density. Their individual contributions can be obtained by solving the following coupled Boltzmann equations, where X stands for a standard model particle. Thus, first term on the right hand side of Eqs. (55) and (56) take care of the contribution of annihilation to SM particles whereas the second terms of both the equations denote the contribution of the selfscattering of the two scalars. In the very early universe, both of the scalars are in thermal and chemical equilibrium. But as the universe expands, the temperature falls resulting some species to be decoupled from the universe plasma and contribute to the relic density. The heavier scalar decouples earlier than the lighter one. But today they are both freezed out giving rise to a total contribution in relic abundance that is probed by WMAP/Planck or other cosmological observations. Hence the total relic abundance (Ω DM ) due to both of the scalars is the sum of their individual contributions (Ω S and Ω S ), If the self-scattering cross-section between the two scalars is small compared to that of the scalars going to standard model particles, such that then the coupled Boltzmann equations in Eqs. (55) and (56) can be treated as two decoupled equations each one of which describes the evolution of each of the scalars independently. In the present work we have ensured this condition by taking the masses of S and S close enough 6 so that σv SS↔S S is negligible from phase space considerations. We have then used micrOMEGAs [77,78] to calculate Ω S and Ω S .
The thermally-averaged values of cross-section ( σv ) for the annihilation channels of dark matter to Standard Model particles (DM + DM → SM + SM ) can be expressed as ( [36], [76] denote the ith order modified Bessel function of second kind. σ(s) is the normalized annihilation cross-section of dark matter for DM + DM → SM + SM processes. For non-relativistic dark matter σv can be approximated as σv ∼σ(4M 2 S )/4M 2 S . The Feynman diagrams for singlet scalar (S or S ) pair annihilation into SM particles in the unitary gauge are shown in Fig. 2. The corresponding expressions for cross-sections can be found in Refs. [17,36]. zones in the mass of DM vs. DM-nucleon cross-section plane. If we can express these cross-sections in terms of the parameters of our model, we can translate the results obtained from direct detection experiments into some preferred or excluded region in the parameter space of our model comprising of δ 2 , M S , δ 2 and M S . The requirement of producing the right relic abundance will further restrict the allowed parameter space.
As presented in Appendix A, the expressions for singlet scalar-nucleon elastic scattering cross-section are rather involved. But for all practical purposes Eq. (64) can be approximated as [17] Similar expression works for S as well.
Empowering ourselves with the expression of cross-section in terms of model parameters, we will now go ahead in constraining the model parameter space from direct detection experimental results and relic density requirements from Planck survey. We have broadly explored three DM mass ranges. CDMS II and CoGeNT vouch for low (∼ 10 GeV) mass DM. CRESST II data prefer a relatively higher mass zone (∼ 25 GeV), in addition to the low zone. XENON 100 and LUX provide exclusion regions from non-observation of any interesting event. Only high mass DM (> 50 GeV) are consistent with these two experiments and Planck data. We will now elaborate more on them in the following.  This restricts Ω S from Planck constraints on relic abundance. Only a small part of the previously allowed region in the M S -δ 2 plane (shown as the blue shaded region) can now account for such Ω S . We then complete choosing our benchmark point (see

Constraints from CoGeNT and Planck data
CoGeNT dark matter direct detection experiment predicts dark matter particle with a mass roughly ∼ 7-11 GeV and elastic scattering cross-section with nucleon which is ∼ 10 −41 -10 −40 cm 2 . The other direct detection experiments like DAMA/LIBRA or CRESST II have also reported signals nearly in that zone which is not consistent with the known background sources. Also the spectrum of the events reported by experiments like CRESST II and CoGeNT are consistent with each other [79] and possibly attribute to dark matter of mass ∼ 10 GeV.
With CoGeNT preferred zone we can do similar analysis as we did with CDMS II. If we are to explain CoGeNT findings with either S (or S ), we find the olive hatched zone in the M S -δ 2 (or M S -δ 2 ) plane (see Fig. 5). We choose a point (M S = 7.8 GeV, δ 2 0.56) in the M S -δ 2 plane. This fixes Ω S . Planck results then restrict Ω S , which then can be reproduced by a tiny region in the M S -δ 2 plane (shown as blue shaded region). We then take M S = 8.2 GeV, δ 2 0.61 from this region to complete the CoGeNT benchmark point presented in Table 2.

Constraints from CRESST II and Planck data
We have analysed the 1σ contour of the CRESST II data. We could have chosen 2σ region of CRESST II data as well. But as the low mass part of the 2σ contour has crossover with CDMS II and CoGeNT low mass regions, the outcome is expected to be similar to these experiments. We rather prefer to work with DM of higher mass ∼ 25 GeV. CRESST II 1σ contour can similarly be translated to the hatched region in the M S -δ 2 (or M S -δ 2 ) plane (see Fig. 6). The best fit point for the CRESST II 1σ contour (mass of DM = 25.3 GeV, σ SI = 1.6×10 −42 cm 2 ) corresponds to the point (M S = 25.3 GeV, δ 2 0.36) in the M S -δ 2 plane. Choice of this point fixes Ω S . The relic abundance constraint from Planck survey then restricts Ω S , which in turn limits the allowed parameter space in the M S -δ 2 plane to the tip (shown as the blue zone in Fig. 6) of the CRESST II allowed zone. We then choose M S = 23.3 GeV, δ 2 0.47 from this region to complete the CoGeNT benchmark point presented in Table 3.

Constraints from XENON 100 and Planck data
XENON 100 collaboration did not observe any prospective signal of DM. From this non-observation they have set an upper bound on spin independent scattering crosssection σ SI for various dark matter masses. In the context of our model this translates to an allowed region in the M S −δ 2 (or M S −δ 2 ) plane, indicated as the olive region in Fig 7(a). Here we assume only S (or S ) participate to have a conservative estimate. We need to find the parameter space suitable for producing DM relic abundance compatible with Planck observations. Our strategy here is a bit different from earlier experiments which indicated some allowed zone in the scattering cross-section -mass of DM plane. For M S − M S = 2 GeV we scan the parameter space spanned by M S , δ 2 and δ 2 . The outcome is presented in Fig 7(b). We see that there is a small island of allowed parameter space at around M S ∼ 55 GeV and a continuum from ∼ 85 GeV onwards. From each of these two regions we will now choose a benchmark point.
First we choose (M S = 54 GeV, δ 2 0.022). This fixes Ω S . This in turn restricts Ω S from Planck data. To reproduce this Ω S window, we do a parameter scan in the M S − δ 2 plane imposing the constraint M S − M S ≤ 2 GeV. We thus get the zone (shown as blue dots in Fig 7(a)) consistent with both XENON 100 and Planck observations. We thus get benchmark point 4A as given in Table 4.
Then choosing (M S = 90 GeV, δ 2 0.05) and proceeding as above leads to the benchmark point 4B as shown in Table 4.

Constraints from LUX and Planck data
LUX collaboration has recently published their results which confirm to XENON 100 findings. As we have done for XENON 100, we show the allowed zones by LUX and Planck observations in Fig. 8. We see that similar to XENON 100, here also we get an island in the parameter space around M S ∼ 57 GeV. But the continuum starts around 135 GeV. For the M S ∼ 57 GeV point the phenomenology will be similar to XENON 100. However high DM masses ∼ 135 GeV can not reproduce the morphological features of indirect detection observations. So we will not discuss the LUX allowed parameter space any further in this work.
To compare with the literature which attempts explaining the experimental observations assuming certain branching fractions of Higgs to SM particles we have denoted in Tables 1-4 the relevant branching ratios for the chosen benchmark points. The main feature of our benchmark points is that in our case the braching ratios are determined from the model precisely whereas the previous analysis has been performed assuming certain branching fractions. Still the experimental data can be confronted remarkably well.

Confronting Indirect Dark Matter Detection Experiments
The region surrounding the Milky Way is rich in astrophysics and is assumed to have a high density of dark matter. This region is promising for better understanding of the properties of dark matter, as no other astrophysical source or region is as accessible as the galactic centre (GC). The Fermi Gamma Ray Space Telescope (FGST) has been employed to survey the high luminous gamma ray emission between ∼ 50 MeV to ∼ 100 GeV. DM distribution follows a density function, ρ( r), with r is the position vector from the centre of the galaxy. Several such DM halo profiles are available in the literature. We choose some representative cuspy to flat profiles for our numerical estimations. These are presented in the Appendix C. For a particular DM halo profile one can calculate the photon flux due to DM annihilation using Eq. (69) in Appendix B.
We now discuss observations of excess γ-ray emission from GC and low latitude of Fermi Bubble, which does not appear to have "standard" origins, but can be understood in the light of DM annihilation. In particular we show that the morphological features of these observations can be explained by our proposed model with parame-ter spaces consistent with DM direct detection experiments and DM relic abundance constraints from Planck survey. We will point out in subsequent discussions that the uncertainties involved in understanding the significance of such astrophysical observations are quite substantial. So at this point we do not intend to fit the data, but rather limit ourselves to reproduce the morphological features of the observations in terms of DM annihilation by our model.

Explanation of Excess Gamma Ray Emission from Galactic Centre
The low energy (few GeV) γ-ray data from galactic centre region observed by Fermi telescope give a hint of a low mass dark matter. In this section we have discussed the phenomenology of gamma ray from the annihilation of dark matter from galactic centre in this present formalism. We have computed γ-ray flux from the singlet scalars in this model constrained by CDMS II and CoGeNT experiments and finally comparison with the observed γ-ray flux has been done. Detailed studies on spectral and morphological features of the gamma rays from the galactic centre region has been studied in Refs. [79,56]. In Ref. [79], the spectrum of the gamma ray emission from the region that encompasses 5 o surrounding the galactic centre is studied after subtracting the known sources from the data of the Fermi Second Source Catalog [80] and disc emission template. The main reason of the disc template emission is the gamma ray produced from the neutral pion decay which is outcome of cosmic ray interaction with gas. Though inverse Compton and Bremsstrahlung can also contribute. Assuming the gas distribution to be of the following form, ρ gas ∝ e −|z|/zsc(R) e −R/Rsc , for R > 7 kpc, where (z, R) denotes the location relative to the GC in cylindrical coordinates. The chosen values of R sc = 3.15 kpc and z sc (R) = 0.1 + 0.00208 × (R/kpc) 2 [81,82] kpc as these values are best suitable to fit for the observational data of 21-cm H line surveys which is the conventional tool used to probe the density of neutral hydrogen. The flux of gamma rays from pion decay is estimated by integrating this distribution over the line-of-sight and it is found to be in good agreement wih the observed morphology of the diffuse emission The residual gamma ray spectrum is brighter for γ-energy range from 300 MeV to 10 GeV and drops by order of magnitude beyond 10 GeV. From the morphological characteristics of this residual gamma ray emission from the central region of our galaxy, it has been shown in [79,57], that below 300 MeV the residual gamma ray could originate from a point-like source but at higher energies it could originate from spatially extended components or may be from annihilating dark matter. Also, if the very high energetic portion of the residual gamma ray emission from the galactic centre is analysed, the spectral shape is found to match fairly well with the gamma emission from galactic ridge. The galactic ridge is an inner region of galaxy extending up to a width of 5 o galactic latitude and ±40 o galactic longitude containing huge amount of white dwarfs [83]. The standard convention is that the high energy cosmic nucleons interact with molecular cloud in the ridge and pions are produced in huge amount which subsequently decay to high energy gamma. Therefore, the residual emission from GC considered here can be assumed to contain low energy tail of ridge emission. Also very low energy part of the residual spectrum is supposed to be dominated by the point source [84,57] which loses its dominance at above GeV scale. However there are some astrophysical propositions that can morphologically explain the gamma-ray flux structure from the inner part of galactic centre. These include millisecond pulsar population [85], central supermassive black hole [56,84] etc. which can explain this spatially extended gamma ray distribution feature from GC.
Super-massive black holes can also accelerate both electrons and cosmic ray protons. These accelerated electrons then produce gamma ray from inverse Compton scattering and can be accounted for any unresolved gamma ray emission from galactic centre region. But these electrons produce γ-ray in TeV-scale [86] which may in principle explain the high energy gamma ray from galactic centre observed by different experiments like HESS, HAWC. Hence this type of mechanism cannot fully account for the FGST data for low energy gamma rays. But the cosmic ray protons accelerated by the black hole can produce pions through the interaction with interstellar gas. Decay of these pions yield gamma rays of lower energy. This scenario may partially explain the FGST residual emission feature as there appear a lot of astrophysical parameters like ISM gas distribution or unknown diffusion coefficient for proton propagation through ISM gas etc. which are not fully understood.
The other astrophysical objects, the gamma rays from which may yield spectra similar to that observed by FGST data are Millisecond pulsars. The spectra from the millisecond pulsars are hard in nature beyond a few GeV, i.e., it falls off with much rapidity after a few GeV. his tend to indicate that surrounding the galactic centre, there may be other millisecond pulsars in considerable numbers which are still to be probed experimentally. But there are discrepancies which immediately contradict this scenario. From the FGST's first pulsar catalog, the spectral index of gamma from pulsars is centred at 1.38 but a much harder spectrum for the average pulsar is required to match the observed gamma spectrum. Although a very few pulsars have certainly a very hard spectral index that can be accounted for the residual   Table 1). SS-annihilation is calculated for M S = 8.6 GeV and δ 2 = 0.45 and is denoted by the dotted violet line. For S S annihilation we use M S = 6.7 GeV and δ 2 = 0.82 and is represented by the dash-dotted cyan line. Total calculated residual γray flux is denoted by the solid black curve. Each sub-figure is calculated for different DM halo profiles, as indicated in their respective captions.
emission [85] below 10 GeV but to fit the gamma flux of the bumpy spectral shape one needs to have larger number of these types of pulsars which is not present in Fermi pulsar catalog. The globular clusters, rich in gamma pulsars have also been studied to measure the spectral index but here too, the data do not favour very hard spectral nature. In order to comply the angular distribution pattern of the emission, the pulsar density should decrease very rapidly along the outward radial distance but the significance of such rapidity has not been found from astrophysical data. From all the above discussions, one may conclude the fact that some different mechanism is required to explain this bumpy spectral shape of the residual emission from galactic centre observed by FGST.   Table 2). SS-annihilation is calculated for M S = 7.8 GeV and δ 2 = 0.56. For S S annihilation we use M S = 8.2 GeV and δ 2 = 0.61. Notations are same as in Fig. 9.
A potentially strong proposition about the nature of this bumpy feature of the residual gamma emission from GC is dark matter annihilation as indicated by Hooper et al. [87,88,89]. As in dark matter scenario, the angular width of the spectra is narrow since the astrophysical factor for flux calculation contains ρ 2 (r) which falls off very rapidly with radial distance from GC explaining the "bump". It also resolves the problem posed from pulsar explanation.
In Ref. [87] it has been argued that by considering few annihilating dark matter scenarios with some standard dark matter halo profiles, low mass dark matter can fit the spectrum with good statistics. Few benchmark cases, such as 10 GeV dark matter annihilating to leptonic channels [90] or 30 GeV dark matter annihilating to bb channel with NFW halo profile have been shown to fit data [87]. In order to   Table 3). get an idea of where our specific model fits in such discussion of generic models, in Tables 1-4 we have quoted branching ratios for different DM annihilation channels for different benchmark points. We see that although for 10 − 55 GeV DM the bb channel has a branching ratio ∼ 80%, but for higher masses, when the W + W − or ZZ channel opens up, it drastically changes. However such a DM candidate is also compatible with data.
In this section we extend the above discussion for the multi-component DM scenario as discussed in our model. As mentioned in the abstract, presence of more than one DM candidate helps enhance the total γ-ray emission due to DM annihilation. We work with benchmark points chosen from the model parameter space already constrained by direct detection experiments and Planck survey. For each such benchmark points we try to match the observed spectra from the theoretically calculated ones.   Table 4). SS-annihilation is calculated for M S = 54 GeV and δ 2 = 0.022. For S S annihilation we use M S = 56 GeV and δ 2 = 0.011. Notations are same as in Fig. 9.
We see that for low mass DM, the plots Fig. 9 and Fig. 10 corresponding to benchmark points 1 and 2 respectively, indicate that a flat DM halo profile like Isothermal profile offers a better agreement with the data. For benchmark point 3, Fig. 11 shows that Isothermal profile is still the promising one, whereas Moore profile overestimates the data. The XENON 100 benchmark point 4A is used for Fig. 12, where we see that for DM masses ∼ 55 GeV, all DM profiles other than the cuspy Moore profile, the DM annihilation contribution is rather small compared to contributions from point sources and galactic ridge. NFW profile works better for   Table 4). SS-annihilation is calculated for M S = 90 GeV and δ 2 = 0.05. For S S annihilation we use M S = 92 GeV and δ 2 = 0.045. Notations are same as in

Explanation of Gamma Ray Bump from Fermi Bubble's Low Galactic Latitude
From the Fermi Gamma-Ray Space Telescope (FGST) data a pair of bilateral lobular structure that contain large amount of gamma-ray had been found in the upper and lower regions of galactic centre. These lobes, known as Fermi Bubbles, emit γ rays between ∼ few GeV to ∼ 100 GeV range and they are extended almost ∼ 50 o (r = ±10 kpc) up and down from the galactic plane. In Ref. [95], the bubble emission has been studied as an extension of WMAP haze [96] which is the non-thermal, microwave emission from the inner part of the galaxy confirmed from data of different ongoing experiments worldwide such as Planck [97] and ROSAT [98] X-ray emission data. Evidences [97] show that near the galactic plane, the γ-ray bubbles and the haze can have a strong correlation that attribute to the fact that they might have been a common origin. When we move far from the galactic plane along the Fermi bubble the gamma ray spectrum follows a power law, E −α with spectral index, α = 2 over all the energy range observed by the FGST. This type of gamma ray spectrum can be well explained by approximate power law spectrum of electron distribution with spectral index, 3, i.e., E −3 where the inverse Compton scattering (ICS) is the mechanism of production of these types of gamma rays. Also, similar distribution can produce radio emission in the galaxy [95,99] due to the synchrotron radiation effect with the interaction of microgauss galactic magnetic field. The picture of the gamma ray emission from Fermi Bubble from the low galactic latitude is somewhat different. In this case the γ-spectrum has a peak at a few GeV energy range. This cannot be explained by the known astrophysical processes like inverse Compton scattering of light sourceby cosmic electrons in steady-state. This non-ICS nature of the spectrum has generated some interest in astrophysics community and different origins for this bumpy nature of γ-ray spectrum from bubble very near the galactic plane have been proposed. Population of millisecond pulsars [100], cosmic ray interaction with gas [101] or an annihilating dark matter scenario [56,87,102,103] have been studied in great details. A detailed study on the morphology and spectral signature of Fermi Bubble is given in Hooper et al. [61].
The explanation of this low latitude excess γ-ray emission as given by the diffuse emission mechanism is due to the fact that the cosmic ray protons are scattered with the gas present in the Milky Way region. But the explanation cannot fully provide the observed phenomenology as the gas distribution is merely correlated with the morphological structure of the γ-emission and also the spectrum of the cosmic ray protons should follow a bumpy nature around ∼ 25 GeV or so to provide a good description of the observed emission. But this type of peak in the cosmic proton spectral feature is not fully understood from the known astrophysical observation. Another possibility of this γ-ray excess can be attributed to the excess population of millisecond pulsars which have the advantage of producing γ-ray emission with very high luminosity over another types of pulsars. But the nature of γ-ray spectrum generated by such objects is not very well understood as a very few pulsars of this type have been discovered. Also the distribution of such objects outside the galactic plane as proposed is much more constrained from various astrophysical observations.
On the other hand the spectral nature of excess gamma emission from the lower latitude of Fermi Bubble may be consistent with the gamma spectrum calculated from the annihilating dark matter scenario at the galactic halo region.
As the DM density is expected to be high for regions close to the GC, we con-  centrate on DM annihilation from low latitude | b | = 1 o -10 o zone of Fermi Bubble. Like we did for GC, we work with benchmark points consistent with direct detection experiments and Planck survey. However, rather than exploring for all DM halo profiles, we present the plot for that DM profile for which the GC data was better explained. These plots are presented in Fig. 14. Here the observed flux is shown after deducting inverse Compton scattering contribution of best-fit steady state electron spectrum of the Bubble. We see that for very low DM mass ∼ 7 − 11 GeV, as preferred by CDMS II or CoGeNT, the spectrum peaks at a lower energy than that obtained from the data. The higher mass zone ∼ 25 GeV, as preferred by CRESST II or ∼ 55 GeV, as allowed by XENON 100 works better. For very high DM masses ∼ 90 GeV allowed by XENON 100, the calculated spectra tend to peak at a bit higher energy than the observed spectrum. But overall we can conclude that the model holds some promise to explain the morphological feature of the Fermi Bubble low latitude γ-ray excess.

Discussion and Conclusion
We have explored a DM model adding two real scalar gauge singlets to the SM. The stability is ensured by a Z 2 × Z 2 symmetry. We keep the symmetry unbroken to get a two component DM scenario. The annihilation of the heavier DM to the lighter ones is suppressed by considering a scenario where the DM candidates S and S are almost degenerate in mass. Such a two component real scalar DM model is better suited to explain both direct and indirect DM experiments compared to the DM models containing a single real scalar.
Detailed calculations for the vacuum stability, perturbative unitarity and triviality constraints on the model has been presented, which forms an integral part of the paper. The parameter space used to explain experimental results do respect these constraints. The model however brings forth the traits of any higgs-portal DM model. For low DM masses, the model predicts unacceptably high invisible Higgs decay width, which calls for adding lighter degrees of freedom to the model. In the present work we did not take up that exercise as we feel that till the conflict of CDMS II or CoGeNT observations with the XENON 100 or LUX observations are settled there is no pressing argument to believe that 7 − 11 GeV DM do exist. So we considered 126 GeV Higgs mediated DM annihilation processes throughout to reproduce indirect DM experimental results.
Guided by the direct detection experiments we considered three DM mass zones. The "low" zone of 7 − 11 GeV is indicated by CDMS II, CoGeNT and CRESST II experiments. CRESST II also favours a "mid" zone ∼ 25 GeV. As XENON 100 and LUX seem to rule out these zones, the only DM masses consistent with both XENON 100 or LUX and Planck observations belong to a "high" mass zone > 50 GeV. The advantage of dealing with this zone is that they do not give rise to unacceptable invisible branching ratio for Higgs. But a too high DM mass > 100 GeV predicts a photon flux from DM annihilations peaked at higher energies than what has been observed in the indirect detection experiments. This high DM mass zone will be probed by future XENON 1T [104] and LUX measurements.
We have chosen some representative "benchmark points" from the parameter space allowed by the direct detection experiments and Planck data. Now the question is how robust are they? DM annihilation cross-sections do depend on M S (or M S ) and δ 2 (or δ 2 ). But as it is proportional to δ 2 2 (or δ 2 2 ), it is quite sensitive to the choice of δ 2 (or δ 2 ). For this reason we choose to show the allowed zones in the plots, so that the reader can roughly estimate the changes in the photon flux in the indirect detection experiments when we choose different benchmark points within the allowed parameter space.
There is some advantage of addressing both direct and indirect detection experiments. The allowed model parameter space is rather restricted by the direct DM detection experiments and relic density constraints as imposed by Planck. This makes indirect DM detection predictions quite sensitive to the assumed DM halo profiles. We would like to make a point that once some agreement in the direct DM sector is established and the background effects in the indirect detection experiments are better understood to delineate DM annihilation effects, in the framework of a given model, the experiments with the existing precision show some promise to identify the right DM halo profile. We have illustrated this with our proposed DM model.
In conclusion, with the proposed model we do not intend to show that it can explain all the experimental results, which sometimes are contradictory in nature. Rather, in the framework of the model, we wanted to exploit the advantage of having a multi-component DM model satisfying both direct and indirect DM experiments, and in this process comment on the viability to choose the right DM halo profile. For completeness we also presented detailed calculations for theoretical constraints on this model as mentioned earlier.

C DM Halo Profiles
The dark matter distribution is usually parametrised as a spherically symmetric profile, ρ(r) = ρ 0 F halo (r) = ρ 0 (r/r c ) γ [1 + (r/r c ) γ ] (β−γ)/α , where α, β, γ and r c are the parameters that represent some particular halo profile listed in Table 5. ρ 0 is the local dark matter halo density at solar location (ρ(r )) taken to be 0.4 GeV/cm 3 with r is the distance between sun to the galactic centre (∼ 8.5 kpc). Another halo profile, namely Einasto profile has also been involved for our study. A different kind of parametric form is adopted in this halo profile [93] which can be written as, whereα is a parameter of the halo profile. In our work value ofα is chosen to be 0.17.