Probing Gravitational Dark Matter

So far all evidences of dark matter (DM) come from astrophysical and cosmological observations, due to gravitational interactions of the DM. It is possible that the true DM particle in the universe joins gravitational interactions only, but nothing else. Such a Gravitational DM (GDM) acts as a weakly interacting massive particle (WIMP), which is conceptually simple and attractive. In this work, we explore this direction by constructing the simplest scalar GDM particle $\chi_s$. It is a $Z_2$ odd singlet under the standard model (SM) gauge group, and naturally joins the unique dimension-4 interaction with Ricci curvature, $\xi_s \chi_s^2 R$, where $\xi_s$ is the dimensionless nonminimal coupling. We demonstrate that this gravitational interaction $\xi_s \chi_s^2 R$, together with Higgs-curvature nonminimal coupling term $\xi_h H^\dag H R$, induces effective couplings between $\chi_s^2$ and SM fields which can account for the observed DM thermal relic abundance. We analyze the annihilation cross sections of GDM particles and derive the viable parameter space for realizing the DM thermal relic density. We further study the direct/indirect detections and the collider signatures of such a scalar GDM. These turn out to be highly predictive and testable.


Introduction
All evidences of dark matter (DM) come from astrophysical and cosmological observations so far, due to the gravitational interactions of the DM. It is possible that Nature may have designed the DM particle to join gravitational interactions only, but nothing else. Such a Gravitational DM (GDM) acts as a weakly interacting massive particle (WIMP), which is conceptually simple and attractive.
The standard model (SM) of particle physics successfully describes the electromagnetic, weak and strong forces in nature, while the gravitation is best theorized by Einstein general relativity (GR). It is apparent that the world is described by the joint effective theory [1] of the SM and GR, which could be valid up to high scales below the Planck mass. We are well motivated to study the intersection between the SM and GR within this effective theory. In this work, we construct the simplest scalar GDM particle χ s , which is a Z 2 -odd singlet under the SM gauge group, and joins gravitational interaction only. As such, there is a unique dimension-4 operator prescribing the interaction between the GDM χ s and the Ricci curvature R, where ξ s is the corresponding dimensionless nonminimal coupling. Since all SM particles enjoy gravitational interaction, gravity can serve as the natural messenger between the GDM and SM particles via (1.1). In the present work, we systematically study the constraints and tests of such a GDM for a variety of dark matter phenomenologies. In passing, we also note that a recent different study considered a gravity-mediated (composite) dark matter model -1 -

JCAP03(2015)052
in the context of warped extra-dimensions, where the radion and massive KK gravitons serve as the mediator [2]. This paper is organized as follows. In section 2, we present the minimal construction of GDM in both Jordan and Einstein frames. Then, in section 3, we analyze the GDM as a WIMP dark matter candidate and identify its viable parameter space for generating the observed dark matter relic abundance. Section 4 is devoted to the systematical analysis of (in)direct searches of the GDM, and the probe of the GDM at high energy hadron colliders. We finally conclude in section 5. Appendix A will present the formulas of radiative loop factors as needed for the physical applications in sections 3-4. In appendix B, we calculate the threshold and resonance effects for dark matter annihilations, which are needed for the thermal relic density analysis in section 3.

Minimal gravitational dark matter
In this section, we present the formulation of the scalar GDM χ s and derive its induced interactions with the SM particles. We first consider the GDM in Jordan frame, where the nonminimal coupling (1.1) is manifest. Then, we make the Weyl transformation on the metric and convert the action into Einstein frame, in which the nonminimal term (1.1) is fully transformed away and result in a new set of effective operators. With these, we will systematically derive the relevant Feynman vertices for χ s in Einstein frame.

Minimal GDM in Jordan frame
Within the joint effective theory of the SM + GR, we can write down the effective action by including this scalar GDM field χ s , where g (J) µν and R (J) denote the Jordan frame metric and Ricci scalar, respectively. 1 The Lagrangian term L F represents the fermion sector of the SM. In eq. (2.1), we define the gauge field strength, F a jµν = (G a µν , W a µν , B µν ), as well as the Higgs doublet field, H = π + , 1 where v EW 246GeV is the vacuum expectation value (VEV) of the SM Higgs at electroweak vacuum. The second line of (2.1) contains the nonminimal coupling terms for the GDM χ s and the Higgs doublet H. According to the constraints from the current LHC Higgs data [4] and from the perturbative unitarity [5,6], the nonminimal coupling ξ h receives an upper limit around O (10 15 ). In the electroweak vacuum, the Higgs non-minimal coupling term makes a contribution to the Einstein-Hilbert action: . Hence, we can identify M 2 + ξ h v 2 EW = M 2 Pl , where M Pl = (8πG) −1/2 2.44 × 10 18 GeV is the 1 We note that ref. [3] considered a real scalar serving as both the DM particle and the inflaton in the early universe. With ξ s ξ h , inflation occurs along the real scalar direction. To account for the cosmic fluctuation strength, the demanded ξ s is far below 10 15 , and thus fully differs from the relevant parameter range of ξ s in the present paper (cf. our figure 3). The nonminimal coupling in ref. [3] is negligible for low energy DM phenomenology, and its DM interacts with SM particles mainly via the conventional Higgs portal coupling λ hχ [cf. (2.2)]. Hence, our current GDM construction realizes a different DM mechanism from ref. [3] and invokes different parameter space of the DM nonminimal coupling. Pl , and thus M M Pl holds to good accuracy. We may also add the higher curvature terms R 2 and R µν R µν to the effective action (2.1) as well. But they do not affect the leading-order graviton contributions, and are irrelevant to the present analysis of scalar GDM.
In eq. (2.1), V (H, χ s ) is the general potential including the SM Higgs doublet H and the scalar GDM χ s . We construct χ s as a Z 2 -odd real singlet, which has vanishing VEV, χ s = 0. Then, we deduce the gauge-invariant scalar potential with CP and Z 2 symmetries as follows, Since we consider that the GDM field χ s joins gravitational interactions only, the χ s has no direct coupling with the SM particles except coupling to gravity and itself. So we will set the Higgs portal coupling λ hχ = 0, or be negligible for the current study. Note that if λ hχ = 0 holds at tree-level, we might expect it to be reinduced via nonminimal couplings (ξ h , ξ s ) due to graviton-exchange. But such graviton-exchanges just induce a new dimension-6 effective operator (H † H)∂ 2 χ 2 s in Einstein frame [cf. eq. (2.6)], 2 which differs from the λ hχ term of dimension-4. We also note that the graviton-loop may induce Higgs portal term. This should have a coefficient proportional to ξ s ξ h Λ 4 /M 4 Pl , where Λ is the UV cutoff for loop integration. Setting Λ as the unitarity bound Λ ∼ M Pl / |ξ s ξ h | (cf. section 2.3), we can estimate the graviton-loop-induced Higgs portal coupling λ hχ ∝ |ξ s ξ h | −1 1, which is negligible for |ξ s ξ h | 1 in the present study (cf. section 3). In practice, we only need to mildly set λ hχ O(10 −2 ) for our construction, which has negligible contribution to the DM thermal relic density. We also note that the Higgs portal term λ hχ H † Hχ 2 s was extensively studied in the literature [7,8] for realizing χ s as a DM. It induces interactions of DM with other SM particles and may provide the DM relic density if the coupling is sizable, λ hχ = O(0.1 − 1). In an extended scheme, we may consider both couplings λ hχ and (ξ s , ξ h ) to give comparable contributions to the DM relic density. But we will focus on the minimal GDM construction for the present study, where the DM interacts with SM particles only via gravity-induced interactions.
From the Jordan frame action (2.1), the dominant interactions for the GDM arise from its nonminimal coupling with the Ricci curvature. We perturb the metric under flat background, g (J) µν = η µν + κĥ µν , where κ ≡ √ 2/M Pl andĥ µν denotes graviton. Since χ s = 0 , there is no mixing betweenĥ µν and χ s . Then, we derive the Feynman vertex for gravity induced triple coupling where the first term comes from nonminimal coupling and is proportional to ξ s . The SM particles couple to gravity minimally through their energy-momentum tensor. The cubic JCAP03(2015)052 couplings of a pair of SM particles with h µν are suppressed by M −1 Pl . In the parameter region of |ξ s | 1 , the interactions between GDM and SM particles induced by graviton-exchange are largely enhanced. Furthermore, since Higgs field mixes with graviton via kinetic term, χ s can communicate with SM particles via Higgs-exchange. Such contributions are proportional to ξ s ξ h , which will be much more enhanced when both |ξ s |, |ξ h | 1. The explicit momentum structures of these interactions are determined by the complicated tensor structure of graviton propagator and the related vertices. We note that the analysis will be much simplified by transformation into Einstein frame. In the following, we will explicitly derive the new set of effective Feynman vertices involving the GDM interactions with the SM particles in Einstein frame.

Minimal GDM in Einstein frame
The Einstein frame is defined by the conventional metric that satisfies Einstein equation. This is achieved by eliminating non-minimal coupling terms via the Weyl transformation. For notational convenience, we will suppress the superscript "(E)" for geometric quantities in Einstein frame. The Weyl transformation is defined as, g µν = Ω 2 g (J) µν , and the factor Ω 2 is given by where |π| 2 = 2π + π − + (π 0 ) 2 . Accordingly, the Weyl transformation of Ricci scalar takes the following form, Substituting this into (2.1), we derive the Einstein frame action for bosonic sector, For nonzero ξ h , the higher dimensional operator in the first line of (2.6) yields additional contribution to the Higgs kinetic term. Together with the original one, we have the following kinetic term for the Higgs and Goldstone boson fields, Hence, we can normalize the kinematic term of Higgs boson by a field redefinition,φ = ζφ, with the rescaling factor, Then, the canonical field φ is identified as the 125 GeV Higgs boson, which was recently discovered at the LHC [9][10][11][12]. We note that this rescaling only applies to the Higgs fieldφ, but does not affect its constant vacuum expectation value v EW . The same operator in (2.6) also induces self-interactions for scalars. For the fermionic sector, we write down the pure kinetic term and mass-term for a generic Dirac spinor f (quark or lepton) in Jordan frame, where e q ν and ω µ mn denote the vierbein and spin-connection, and σ mn = i 2 [γ m , γ n ]. Setting the flat background in Einstein frame, we deduce the metric in Jordan frame, g (J) µν = Ω −2 η µν . Thus, we can express the vierbein and spin-connection in Jordan frame as functions of Ω, With these, we can explicitly write down the kinetic term and mass-term for the SM fermions (quarks or leptons) in the Einstein frame [5,6], In the following, we summarize the vertices relevant for DM annihilation processes, by expanding Ω at the leading order of 1/M 2 Pl .
The only triple coupling relevant to the following analysis comes from the vertex χ s −χ s − φ. It induces Higgs invisible decay when M χ < 1 2 m φ , and also generates interactions between the GDM and SM particles by exchanging the Higgs boson. We derive the corresponding Feynman vertex at the leading order, where all momenta flow inwards. Then, we deduce the quartic couplings between the GDM χ s and Higgs/Goldstone bosons at the leading order, where q = p 1 + p 2 . They also contribute to the dark matter annihilations in early universe and today. To obtain leading order contributions form Higgs exchange to these vertices, we need the following triple couplings, (2.15) Then, we deduce the quartic coupling for the vertex χ s (p 1 )−χ s (p 2 )−π +,0 (p 3 )−π −,0 (p 4 ) at the leading order, where Γ φ stands for the Higgs boson width. For the quartic coupling with Higgs bosons, t(u)-channel exchange of χ s also contribute, and the vertex where t = (p 1 −p 3 ) 2 and u = (p 1 −p 4 ) 2 . The quartic couplings for the 4φ and 4χ s vertices as well as for the Higgs-Goldstone interactions receive quite similar contributions. They will be included in our coupled channels analysis of perturbative unitarity.
• GDM Interactions with Weak Gauge Bosons: Under Weyl transformation the gauge boson kinetic terms remain intact as in eq. (2.6). We note that the tree-level interactions between the GDM and massive gauge bosons arise from the gauge boson mass-term. For weak gauge bosons, this is associated with the Higgs kinetic term in eq. (2.6). Thus, we derive the interaction term, − Pl δ V V µ V µ χ 2 s , with the notation V ∈ (W, Z) and coefficients (δ W , δ Z ) = (2, 1). Hence, we infer the Feynman vertex of gravity-induced contact interaction for the GDM and weak bosons,

JCAP03(2015)052
Besides, the GDM can interact with weak bosons by exchanging the Higgs boson. With the gravity-induced χ s − χ s − φ vertex in eq. (2.13) and the V µ −V ν −φ vertex from the SM, we derive the following contribution via Higgs-exchange at the leading order, where q = p 1 + p 2 . Then, we deduce an effective (nonlocal) vertex for (2.20) For |ξ h |, |ξ s | 1, it is the quadratic term of ξ h ξ s in eq. (2.19) that will make dominant contribution. For the scattering process V L V L → χ s χ s , the non-renormalizable gravityinduced interactions will contribute a net E 2 -dependence in the amplitude, and cause perturbative unitarity violation at high energies. Furthermore, this vertex will lead to the GDM pair-productions via weak boson scattering V V → χ s χ s at the LHC and future high energy pp colliders.
• GDM Interactions with Fermions: According to eq. (2.11), the dark matter can interact with fermions via their kinetic terms or mass-terms. Intuitively, the kinetic terms in the parentheses seem to induce momentum-dependent higher dimensional operators with Ω −4 1 − 2χ 2 s /M 2 * . But, for on-shell fermions, the contributions from kinetic terms share the same structure as that from mass-terms. The total contribution to the contact interaction is s . In addition, the Higgs-exchange induces a nonlocal contribution to the same vertex at the leading order. Thus, we explicitly derive the Feynman vertex We note that the terms in the brackets of (2.21) and (2.20) take the same form. At high energies, the scattering amplitude of χ s χ s →f f contains non-canceled leading E 1 terms, which will eventually violate perturbative unitarity as the energy E increases [13][14][15].
• GDM Interactions with Massless Gauge Bosons: As shown in eq. (2.6), the gauge boson kinetic terms remain intact under Weyl transformation. So there is no contact interaction of the GDM with massless gauge bosons (gluons or photons) at the leading order. Nevertheless, there are loop-induced higher dimensional operators. For instance, the dimension-6 operator χ 2 s G a µν G aµν can be generated by the top quark triangle-loop in figure 1, where the diagram (a) involves leading order χ s χ sf f contact interaction, and the diagram (b) includes Higgs-exchange with Higgs effective coupling to gluons. It will initiate gluon-fusion production of χ s χ s at the LHC and the future high energy hadron colliders. In parallel, the dimension-6 operators χ 2 s A µν A µν and χ 2 s A µν Z µν can be generated from both W ± loop and fermion loop, and are relevant to the indirect detections of dark matter. Inspecting (2.20) and (2.21), we note the similarity between the contact vertex for χ 2 s V V (χ 2 sf f ) and the corresponding φV V ( φf f ) vertex. With the same structure, couplings of the former can be reproduced from the latter by the substitution v EW → −M 2 Pl /ξ s . Hence, we can directly infer the form of these one-loop generated vertices from the conventional results for the SM Higgs boson [16,17] as follows, where the form factors (C g , C γ , C γz ) are energy-dependent, The explicit expressions of A V,F (τ ) and B V,F (τ, η) are given in appendix A. For analysis of non-relativistic dark matter annihilations in section 3, we will have q 2 ≈ 4M 2 χ .

Perturbative unitarity
In this subsection, we derive perturbative unitarity bound from high energy scattering processes involving the GDM, as induced by the non-renormalizable gravitational interactions in section 2.2. For gauge bosons and the Goldstone bosons, the leading order amplitudes in the high energy limit are given by O(E 2 ) terms. Thus, we derive the following amplitudes, where E = √ s is the center of mass energy of the scattering. Here, we keep the leading order . We compute the scattering amplitudes for both the longitudinal gauge boson final state V a L V a L and the corresponding Goldstone boson final state π a π a . This verifies the equivalence theorem [18][19][20][21][22][23] at high energies and serves as nontrivial consistency checks of our analysis. For the Higgs final state, we find that the leading amplitude is given by T [χ s χ s → φφ] T [χ s χ s → π a π a ], in the high energy regime, which arises from the contact interaction (2.14) at O(1/M 2 Pl ). To derive the optimal perturbative unitarity constraint, we further perform the coupled channel analysis for the normalized two-body scalar states, |π 2 |χ s χ s , |π 0 χ s and |φχ s . The partial wave amplitude is given by We inspect the leading contributions at O(E 2 /M 2 Pl ). The leading amplitudes without involving χ s were derived before in ref. [5,6]. Combined these with the amplitudes of (2.24) and related results, we deduce the full s-wave amplitude in matrix form, The submatrices in (2.26) take the following form, where we only keep the leading E 2 -terms under the limit |ξ h |, |ξ s | 1. In the above formulas (2.27b), Here, for the convenience of applying the unitarity conditions below, we have included proper kinematical phase factor of each scattering channel (such as η 12 and η 33 ), which were generally defined in appendix B of ref. [14]. In the present unitarity analysis, it suffices to keep only the mass M χ (which could reach TeV scale) and ignore other small masses of weak bosons and Higgs boson in comparison with the large scattering energy E. Thus, in eq. (2.27), only the scattering channels involving external χ s state have nontrivial phase factor η ij = 1. After diagonalization, we deduce the eigenvalue amplitudes, The s-wave amplitude should obey the unitarity condition |â 0 | < 1 (or, |Reâ 0 | < 1/2) [14]. Imposing condition |â 0 | < 1 on the maximal eigenvalue, we derive the unitarity bound Λ U = E max , . Defining the coupling ratio, r ≡ |ξ s /ξ h |, we can reexpress (2.29) as an upper bound on |ξ h ξ s | for each given energy E, Here, the strongest limit corresponds to E = E max = Λ U , which serves as an ultraviolet (UV) cutoff of this effective theory. In our present study, we will set up the parameter space |ξ s | > |ξ h | 1, where typically we take the coupling ratio r = 5 − 30. It is clear that for the range of r = 5 − 30, we have r −2 1 and thus the bound (2.30) is not so sensitive to the ratio r.
Besides, since ξ 2 s > 0 and ξ 2 h > 0 in (2.29), we can always derive an upper bound on ξ h alone (for each given energy scale E ), We further note that the scattering amplitude of χ s χ s → χ s χ s vanishes at O(E 2 ) due to the crossing symmetry, we may further consider its subleading terms at O(E 0 ), which is still enhanced by ξ 2 s . From the Lagrangian (2.12), we derive the following amplitude, which leads to the s-wave amplitude, for ξ s 1. Thus, imposing the unitarity condition |η 12 a 0 | < 1 [14], we deduce the following bound for ξ s 1, , we see that the bound (2.35) is independent of the scattering energy E, but inversely suppressed by the dark matter mass M χ .
The second class of processes involves a pair of fermions, φφ → ff or χ s χ s → ff . At high energies, their amplitudes are dominated by O(E 1 ) terms. The amplitudes φφ → ff and χ s χ s → ff are enhanced by ξ h and ξ s , respectively. For |ξ h |, |ξ s | 1, the unitarity bounds (2.29) and (2.34) from pure scalar scatterings are much stronger than these processes with a fermion pair.

Analyzing thermal relic density of GDM
In this section, we study the property of the GDM χ s as a WIMP dark matter candidate. We explore the intriguing possibility that the GDM alone fully accounts for the observed thermal relic density. From this, we will analyze the viable parameter space for the GDM. We find three independent parameters involved for this analysis: the dark matter mass M χ and two nonminimal couplings (ξ h , ξ s ). As we will elaborate, in most of the parameter space, the prediction of thermal relic abundance is only sensitive to the product of two nonminimal couplings ξ h ξ s . Hence, our GDM construction is very economical and highly predictive.
In figure 2, we display all channels for the GDM annihilations into the two-body final states at leading order, 5 χ s χ s → V V, φφ, ff , where V = (W, Z) and f = (s, µ, c, τ, b, t). Given the present sensitivities to WIMP via various experimental searches, we will consider the GDM mass-range, O(1GeV) M χ O(1TeV). From eq. (2.17) and eqs. (2.20)-(2.21), we see that the gravity-induced couplings to fermion and gauge boson are proportional to their masses. Hence, for heavy mass-range O(100GeV) M χ O(1TeV), the dark matter annihilations are dominated by the channels χ s χ s → W + W − , ZZ,tt, φφ , while for light mass-range M χ O(10GeV) only the annihilation channels χ s χ s →bb,cc, τ τ,ss, µµ are allowed at the freeze-out temperature.
Before performing systematical numerical analyses, we may first estimate the required size of the nonminimal couplings (ξ h , ξ s ) for accommodating the DM thermal relic density (when M χ is away from any threshold or resonance). For s-wave annihilation, the cross section dictated by the relic density abundance is, (3.1) 4 Inspecting (2.6), we see that the interactions involving more than 4χ s need further expansion of 1/Ω 2 , which brings in additional ξ s /M 2 Pl for each pair of χ s . Thus the scattering channel χ s χ s → χ s χ s places the best unitarity constraint on ξ s . 5 For channels with final states containing more two particles, they may come from decays of off-shell heavy particles. For instance, we can estimate the size of χ s χ s → W W * → Wf f . In the intermediate mass-range and further suppressed by the on-shell W momentum for decreasing M χ . Thus, it is reasonable to just count on the leading two-body annihilation channels. Also, the loop-induced annihilations from the effective vertex (2.22) are negligible. For M χ much heavier than the weak scale, the following heavy modes dominate in final states. At leading order, the thermal averaged cross sections equal the zero-temperature expression with s 4M 2 χ , where we consider the parameter region of |ξ h |, |ξ s | 1. Note that the leading order contributions for gauge boson final states come from the longitudinal modes and coincide with that of the Higgs final state. For M χ = O(10 2 − 10 3 )GeV, the product of nonminimal couplings is required to have a size around |ξ h ξ s | ∼ O(10 14.5 ). When M χ 100GeV, the heavy final states W W, ZZ, hh, tt decouple in advance, and the cross section is dominated by the annihilations χ s χ s → bb, cc, τ τ .
To compute thermal relic density of dark matter with a wide mass range, we will take into account threshold and resonance effects in the thermal integration [24]. Compared with zero temperature case, annihilations into the final state slightly heavier than dark matter could be active, due to its Boltzman distribution at finite temperature. Also, when dark matter annihilate near the Higgs mass pole, i.e., 2M χ ∼ m φ , the cross section is largely enhanced over the case away from the pole. To properly treat the cross section around the pole, we will make thermal integration numerically without any expansion for velocity. Using the couplings of (2.17) and (2.20)-(2.21), we derive zero temperature annihilation cross JCAP03(2015)052 section for the relevant final states, where f = t, b, c, s, τ, µ and N c = 3(1) for quarks (leptons). We also denote V = W, Z and (δ W , δ Z ) = (2, 1). For the parameter range of interest, we only keep the leading order contributions under |ξ h |, |ξ s | 1. Thus, the cross sections are controlled by the DM mass M χ and the product of nonminimal couplings ξ h ξ s . The width of a SM Higgs boson with mass 125 GeV is rather small, Γ SM φ 4.03MeV. In our model, the Higgs total width Γ φ could deviate from the SM value only when the invisible decay channel φ → χ s χ s is open. Since the DM mass-range for active annihilation process χχ → φφ is far away from resonance region, we can safely neglect Γ φ in eq. (3.3c). We present the calculation of thermal averaged cross sections by including the threshold and resonance effects in appendix B.
Given the thermal average cross section σ A v as function of the DM mass M χ and coupling product ξ h ξ s , we will derive thermal relic abundance. It is convenient to define a ratio x ≡ M χ /T . Thus, the freeze-out temperature x f = M χ /T f can be derived from the following formula to a good accuracy [25], where σ A v f denotes the thermal averaged cross section at the freeze-out temperature T f . We denote the total effective relativistic degrees of freedom as g * and its counterpart for entropy as g * S . We further assume that all species in the universe have the same temperature, and g * g * S . The coefficient c is a free-parameter for fitting the numerical solution, and we use the conventional choice: c(c + 2) = 1 [25]. Then, the above equation is simplified as Integrating out the differential equation for number density per comoving volume, we can derive thermal relic abundance, where ρ c /h 2 = 1.88 × 10 −29 g cm −3 , and around the freeze-out, g * g * S = 106.75. In the context of ΛCDM scenario, the latest measurement from Planck satellite gives, Ω χ0 h 2 =   0.1199 ± 0.0027 [26]. In the minimum setup, the relic abundance is only sensitive to M χ and |ξ s ξ h |.
In figure 3, we present the contours of thermal relic density in M χ − |ξ h ξ s | plane, where the red solid curve corresponds to Ω χ0 h 2 0.12. For illustration, we focus on the GDM mass range above M χ 48GeV. This is the lower bound set by the Higgs invisible decay constraint (cf. figure 6), above which the narrow width assumption of Higgs holds well. With the increase of M χ , more heavy SM modes contribute to the annihilation cross sections, and thus the required |ξ h ξ s | becomes smaller. Around the Higgs mass-pole, 2M χ m φ , the required |ξ h ξ s | becomes almost one order of magnitude smaller due to the resonance enhancement. For M χ m W , the forbidden channel χ s χ s → W W contributes right below the threshold due to the thermal fluctuations shown in appendix B. Thus, the required nonminimal coupling for realizing the relic density could be smaller than what expected for zero temperature case. From eq. (3.6), we note that Ω χ0 h 2 is roughly proportional to (ξ h ξ s ) −2 via the thermal averaged cross section, i.e., |ξ h ξ s | ∝ (Ω χ0 h 2 ) −1/4 . This means that the solid curve in figure 3 is very insensitive to the experimental uncertainty of Ω χ0 h 2 . Hence, provided that χ s dominates in thermal relics, the constraint on |ξ h ξ s | is so robust that almost no visible variations from the solid curve are allowed.
The present analysis is fully based upon perturbative expansion in the effective theory formulation [1], which combines the SM with nonrenormalizable Einstein general relativity. Such an effective theory normally has an UV cutoff scale Λ U , above which the perturbative expansion breaks down and new physics is expected to show up. 6 For the validity of our JCAP03(2015)052 perturbative analysis, we should derive perturbative unitarity constraints on the parameter space. The shaded areas in the two plots of figure 3 depict the unitarity violation regions from two different unitarity conditions (2.30) and (2.35). The first condition (2.30) is derived from O(E 2 ) leading terms of the scattering amplitudes. It is insensitive to the DM mass M χ , as shown in figure 3(a). It is also insensitive to the coupling ratio r = |ξ s /ξ h | for r 5. So we will take a sample input of r = 5 for illustration. We present the unitarity violation regions in figure 3(a), for the sample UV cutoff (set by unitarity bound) Λ U = 5, 10 TeV, respectively. We see that the unitarity constraint is mild, and our effective theory remains perturbative for a wide range of M χ . The second unitarity condition (2.35) is derived from χ s χ s → χ s χ s channel at O(E 0 ). This is shown in figure 3(b), where the bounds are quickly enhanced with the increase of DM mass. Here, the shaded purple and pink regions correspond to the sample cutoff scale Λ U = 5, 10TeV, respectively. In each case, we have presented the unitarity constraints for the coupling ratio r = 5 (darker shaded area) and r = 30 (lighter shaded area). This plot shows that for a valid perturbation analysis, the ratio r cannot be too large, namely, the values of the two nonminimal couplings ξ s and ξ h should not have a large hierarchy. Finally, we note that the current LHC measurements on the Higgs signal rates put a mild constraint on the coupling ξ h via the kinetic rescaling factor ζ in (2.8). For instance, from the latest CMS (ATLAS) data [11,12], we can infer the 3σ upper bound, |ξ h | < 3.4×10 15 (|ξ h | < 2.3 × 10 15 ). In the following analysis, we will study various experimental searches of GDM within the viable parameter space that generates the DM thermal relic abundance and obeys the perturbative unitarity bounds (figure 3) as well as the current LHC bound on the nonminimal coupling ξ h .
As a final remark in this section, we note that even though the GDM can have a large nonminimal coupling ξ s , it would not cause any large effect on the long distance gravitational behavior. Since GDM field has no VEV, its nonminimal coupling (1.1) in Jordan frame does not contribute to the Planck mass. This nonminimal coupling term is fully transformed away in Einstein frame, and is replaced by a set of higher dimensional operators involving effective interactions between the GDM and SM fields (section 2.2). These dimension-6 operators are suppressed by (ξ 2 s , ξ 2 h , ξ s ξ h )(v 2 EW , E 2 )/M 2 Pl and do not cause any sizable effect at long distance (low energy). They are further constrained by perturbative unitarity bounds (section 2.3 and figure 3) at high energies. These GDM effective couplings properly generate the observed DM relic density. It means that our GDM belongs to a kind of WIMP dark matter and does not cause extra visible change at long distance.

GDM detections and collider searches
In this section, we explore various searches of the GDM. With the gravity-induced interactions between the GDM and SM particles, we find it is quite difficult to probe the GDM by direct detections due to the small-momentum suppression. On the other hand, indirect detections can be promising to reach the parameter space that successfully accounts for thermal relic abundance. Finally, we study the collider searches of the GDM from the Higgs invisible decay, and further discuss the probe of a heavier GDM χ s at the LHC (14 TeV) and future high energy hadron colliders. O(10TeV) will reveal its Kaluza-Klein modes at energies above this scale.  The region below black short-dashed curve denotes the parameter space sensitive to neutrino background. The red solid curve is our prediction which accounts for the DM thermal relic abundance Ω χ0 h 2 = 0.12.

Direct detection of GDM
In the present model, the GDM scalar χ s interacts with nucleons via the gravity-induced interactions between light fermions and χ s . As a scalar dark matter, the GDM-nucleon interaction is spin-independent. From eq. (2.21), we derive the GDM-nucleon scattering cross section under the limit |ξ s |, |ξ h | 1, where f N is the effective form factor, which can be estimated from the QCD chiral perturbative theory, the pion-nucleon scattering and the lattice simulations. We will use f N = 0.345 for the following estimate [27]. For the nucleon mass, we input the averaged mass of proton and neutron, m N = 0.939GeV. The typical scale of momentum-exchange is around 100MeV, i.e., t ≈ −(100MeV) 2 . In the parameter space of interest, the second term in the parentheses of eq. (4.1) is dominant. 7 Comparing with annihilation cross section at √ s 2M χ , we see that the GDM-nucleon cross section has a momentum-suppression factor (t/m φ ) 2 , due to the energy-dependent structure of χ s − χ s − φ vertex. In figure 4, using the data of relevant direct detection experiments, we present a summary of their exclusions (at 90% C.L.) in M χ − |ξ h ξ s | plane. For M χ 10GeV, the strongest constraint comes from LUX experiment [28], as depicted by blue solid curve. The 7 This approximation numerically holds well for ξ h within the perturbatively unitary range. future reach of Xenon1T projection [29] is represented by the blue dashed curve. Due to the low-momentum suppression of σ SI in eq. (4.1), the magnitude of |ξ h ξ s | as dictated by the observed DM thermal relic abundance (red solid curve) is even below the required sensitivity to the neutrino background and out of the reach of direct detection. Hence, the GDM that can fully account for the thermal relic abundance is unlikely to be detected through the nucleus recoil. This is a typical feature of the GDM in contrast to other WIMP DM candidates.

Indirect detection of GDM
Many astrophysical experiments aim at finding indirect evidences of dark matter annihilations in the sky. In the present model, we find that searching for gamma ray signals from target with high dark matter density is most promising. 8 There are two types of gamma ray signals. One is monochromatic photon "line" arising from the dark matter annihilation χ s χ s → γX, where X denotes any other possible SM bosons. The other one is a diffuse continuum spectrum from secondary production of photons from primary dark matter annihilation χ s χ s → W + W − , ZZ, bb, τ + τ − , µ + µ − . The secondary photon may be initiated from final state radiation or hadronization with decays π 0 → γγ. In the following, we study the impacts of these measurements on our model in turn.
As discussed in section 2.2, the effective operators for χ s χ s → γX can be induced from gravitational interactions at one-loop order. From (2.22), we infer the zero temperature cross section for χ s χ s → γγ and χ s χ s → γZ , where α 1/128 is the fine structure constant, and C γ , C γz are energy dependent loopfactors defined in (2.23). The annihilation channel χ s χ s → γZ is active for M χ > m Z . In comparison with tree-level processes, these two channels are suppressed by a loop factor. Nevertheless, since the "line" shape search features a better sensitivity than that of the continuum spectrum, the constraint from monochromatic spectrum would be potentially important. Provided that the DM annihilations into γγ and γZ are the only sources to generate gamma ray line, it is possible to extract an upper bound on the quantity 2(σ A v) γγ + (σ A v) γZ from galactic center γ-ray line search [36], i.e., Fermi-LAT in low photon energy range [37] and H.E.S.S in high energy range [38]. The limits depend on the DM halo profiles as well as the signal region of interest (selected by the experimental group for analyses).
To demonstrate the potential of these experiments for testing our model, we present the strongest constraint from the gamma-ray line search [36] for illustration. Figure 5 Figure 5. Constraints in M χ − |ξ h ξ s | plane from indirect DM detections. (a). 95% C.L. exclusions from gamma ray line search. The light brown curve depicts the strongest limit from FermiLAT in the low photon energy region, and the dark brown curve denotes the limit from H.E.S.S [36]. The areas above these curves are excluded. (b). Exclusions from gamma ray continuum spectrum. The shaded regions are excluded at 95% C.L. The blue, light green and dark green curves represent the bounds from detections via three primary annihilation channels W + W − (ZZ), bb and τ + τ − , respectively. The solid and dashed curves denote bounds from Fermi-LAT and CAT (projection), respectively. In each plot, the red solid curve gives the prediction by realizing the GDM thermal relic density Ω χ0 h 2 = 0.12.

JCAP03(2015)052
The light brown curve is extracted from the searches of FermiLAT [37]. In the intermediate mass range 80 − 160 GeV, since the line shape is sensitive to relative strength of the two processes, no reliable model-independent limit could be inferred [36]. The dark brown curve is extracted from H.E.S.S [38]. As before, the red solid curve is our GDM prediction by accommodating the DM thermal relic abundance. We see that the GDM with mass between 60 − 80GeV is already excluded by FermiLAT. In low mass range below m φ /2, due to the resonance enhancement from thermal integration, i.e., (σ A v) σ A v , the GDM prediction is still viable. For the GDM mass M χ > 80 GeV, the relic density is dominated by the tree-level annihilation into heavier final states. In this mass range, our prediction is significantly below the reach of the gamma ray "line" searches.
Next, we study constraints on dark matter annihilation cross sections from diffuse continuum spectrum. The latest results come from the 4-years data of Fermi-LAT observation of 15 Milky Way dwarf spheroidal satellite galaxies [39]. In the future, the next generation experiments with better angular resolution (such as CTA [40]) will largely improve the sensitivity over a wider mass range. Normally, the upper limit on DM annihilation cross sections is extracted by assuming 100% branching fraction for each primary annihilation channel. Since these limits are sensitive to the spectrum shape for each channel, they could not be straightforwardly mapped to a given model where all annihilation channels contribute in a certain pattern. Nevertheless, following ref. [36], we may estimate the conservative constraint by taking into account the fraction of each channel in the total annihilation cross section. The bound is derived as follows, jj is the experimental upper bound. BR jj ≡ (σ A v) ii /(σ A v) tot , with (σ A v) jj defined in eq. (3.3a) and (σ A v) tot summing over all these channels. Note that BR jj is only a function of mass M χ , and is insensitive to (ξ h , ξ s ). We then deduce the lower bound on |ξ h ξ s | from (σ A v) 95%,res jj for W + W − (ZZ), bb, τ + τ − , respectively. 9 The constraints are summarized in figure 5(b), where the blue, light green and dark green curves represent three primary annihilation channels W + W − (ZZ), bb and τ + τ − , respectively. The shaded regions above solid curves are excluded by Fermi-LAT experiment at 95% C.L., and the dash curves present the sensitivity of CAT. In the mass range M χ 100GeV, the strongest constraint on our model comes from measurements of W + W − (ZZ) channels. Our prediction from the GDM thermal relic abundance is within reach of the future indirect detection experiments. Recently, some studies suggested the gamma ray excess from Galactic Center [41], which can be interpreted as a signal predicted by a 31 − 40GeV dark matter annihilating mostly into bb final state with cross section (σ A v) = (1.4 − 2.0) × 10 −26 cm 3 s −1 . In the present model, the dominant annihilation in this intermediate mass range is indeed the bb channel, but the required parameter range in M χ − |ξ h ξ s | plane is already excluded by Higgs invisible decays.

Collider searches for GDM
The GDM may be produced at hadron colliders in several ways. For a light GDM with mass M χ < m φ /2 , it can be produced via invisible decays of the 125 GeV Higgs boson, φ → χ s χ s , JCAP03(2015)052 due to the cubic vertex χ s − χ s − φ in (2.13). For |ξ h |, |ξ s | 1 , we deduce the invisible decay width, Accordingly, its invisible decay branching fraction is given by where Γ SM φ 4.3MeV denotes the Higgs decay width from the SM contributions alone. Currently, LHC searches invisible Higgs decays via the vector boson associated production, vector boson fusion, and top associated production. By assuming the SM production rate, the best upper limit on the invisible branching fraction comes from combining all existing measurements at the LHC, BR inv < 40% at 95% C.L. [42]. Setting BR inv = BR χχ , we can translate this limit into a constraint on our GDM parameter space. We present this constraint in the M χ − |ξ h ξ s | plane, as shown in figure 6. Since the invisible width Γ(φ → χ s χ s ) is proportional to |ξ h ξ s | 4 , the constraint on |ξ h ξ s | is insensitive to either M χ or the kinetic rescaling factor ζ for φZZ vertex (given that ξ h itself obeys the LHC bound). It is also rather insensitive to the experimental limit on the invisible decay branching fraction around BR inv = O(0.1). Since the same cubic vertex φχ s χ s determines both the Higgs invisible decays and the GDM thermal relic abundance, we find that the LHC bound on Higgs invisible decays puts a nontrivial constraint on the required coupling |ξ h ξ s | for generating the observed thermal relic density. As shown in figure 6, we deduce that the GDM with mass M χ < 48GeV is excluded at 95% C.L. The GDM effective interaction to light fermions in eq. (2.21) initiates the χ s χ s production via quark annihilation at hadron colliders. With a mono-jet, photon and W/Z radiation from the initial state quarks qq , the / E T may be observed. This type of processes has been extensively studied in literature [43][44][45]. Given the null result, we can infer an upper bound on |ξ h ξ s | as function of the GDM mass M χ . In the high energy regime, q 2 m 2 φ , the interaction (2.21) amounts to an effective operator χ s χ sf f . Constraints on the cutoff scale of various effective operators were derived from combining the results of different initial states measured by ATLAS and CMS at LHC (7 TeV) [46]. For scalar type operators in our model, lower bound on M Pl / |ξ h ξ s | is around O(10)GeV. The improvement at the LHC (8 TeV) [47] and the sensitivity estimated for the LHC (14 TeV) search are fairly mild [43][44][45]. We find that these bounds are quite weak as compared to the constraint from Higgs invisible decays in figure 6. The case is further studied for the high luminosity LHC and future pp colliders [48], but the limit is improved by no more than a factor of 10. Thus, it appears uneasy to probe the effective interactions between the GDM and light fermions at hadron colliders.
The GDM may be produced by gluon fusions as well, via loop-induced effective operator χ 2 s G aµν G a µν in eq. (2.22). As we learn from the SM Higgs production, the suppression from one-loop factor can be compensated by the large gluon parton distribution function in high energy pp collisions. So the gluon fusions provide the most significant production at the LHC (14TeV). However, in comparison with the Higgs production, the loop-factor C g in eq. (2.23) for the GDM production is energy-dependent and diminishes for √ s m f . Thus, the production cross section in high energy pp collisions becomes much smaller than what is expected for the SM Higgs production. The gluon fusion production of DM as induced by  Figure 6. Constraint from searching for Higgs invisible decays at the LHC [42]. The yellow region is excluded at 95% C.L. The red solid curve presents our prediction by generating the observed DM thermal relic density Ω χ0 h 2 = 0.12. the top-DM effective operator has been analyzed by using mono-jet searches at the LHC [49]. This greatly improves the sensitivity over the standard search based on light fermion effective operators. But, due to the loop-factor suppression, the lower bound on cutoff scale is around 100 GeV, which is still weaker than that derived from generating the thermal relic abundance by the GDM.
In our model, the GDM has much larger couplings to heavy particles. Thus, we can effectively produce the GDM pair through its interactions with the third generation quarks (t, b) or the vector bosons (W, Z). In figure 7, we present two additional production mechanisms for probing the GDM particles, i.e., the top pair (bottom pair) associated production [diagram (a)] and the vector boson fusions [diagram (b)]. The black dots denote the effective interactions (2.21) and (2.20). In the high energy regime q 2 m 2 φ , the propagator suppression for the dominant Higgs exchange diagram is compensated by energy enhancement in the χ s −χ s −φ vertex, and the effective interactions become contact. The top pair associated DM production can effectively probe scalar-type interactions between the DM and quarks [50]. This is a typical feature of our GDM in the present model. Recently, CMS presented the analysis of this process in di-lepton final states for Dirac DM [51]. The lower bound on the cutoff scale is around 100 GeV for M χ 100 GeV, and decreases in higher mass range. So far, this limit is still too weak to constrain the GDM prediction. Given the heaviness of top quark, we expect a significant improvement of sensitivity in this channel at future circular pp colliders (50-100 TeV) [52]. 10 The DM pair production via vector-boson-fusions (VBF) was studied lately at the LHC (14 TeV) in the context of SUSY models [53]. But, the effective JCAP03(2015)052 vertex χ 2 s V µ V µ is not yet studied. 11 It is encouraging to perform systematical Monte Carlo simulations of this type of operator at the upcoming LHC runs with 13 − 14 TeV collision energy. Further studies at the future pp colliders (50-100 TeV) [52] 10 should effectively probe the heavier mass range of the GDM. This is fully beyond the current scope and will be considered elsewhere.

Conclusions
All the astrophysical and cosmological evidences of dark matter (DM) so far have demonstrated the role of its gravitational interactions only. An intriguing possibility is that the DM communicates with our visible world only via gravitation. In this work, we presented a minimal construction of such a gravitational dark matter (GDM), where a scalar GDM particle χ s couples to the SM through the unique dimension-4 operator (1.1) which contains the fields χ 2 s and Ricci curvature R. In section 2, we formulated this minimal GDM in both Jordan frame and Einstein frame. The GDM χ s is a real singlet scalar and odd under the Z 2 symmetry, which may serve as a WIMP DM candidate. In Jordan frame, both the dark matter particle χ s and the SM Higgs boson φ have gravitational interactions (2.1) with nonminimal couplings ξ s and ξ h , respectively. Due to the graviton-exchange and the graviton-Higgs kinetic mixing, the interactions between the dark matter χ s and SM particles will be enhanced by the coupling product |ξ s ξ h | 1, besides the suppression factor (v 2 EW , E 2 )/M 2 Pl . In Einstein frame, these effective interactions become manifest, as shown in eqs. (2.6) and (2.11). Our model only invokes three key parameters in the DM phenomenology: the GDM mass M χ and the nonminimal couplings (ξ s , ξ h ). For convenience of physical analysis, we derived all relevant Feynman vertices for the GDM in Einstein frame. We also derived the perturbative unitarity constraints on the new couplings (ξ s , ξ h ), and identified the valid perturbative parameter space in figure 3, which justifies our leading order analysis for the GDM. 11 We note that Ref. [54] had recast the CMS search for invisible Higgs decays in the VBF production and converted the CMS results into a bound on the DM mass with the effective operator approach. But, it focused on the gauge-invariant operator χ 2 s F µνa i F a µνi , which could be generated only at one-loop in our GDM model. Under the limit q 2 m 2 φ , our tree-level effective interaction (2.20) induces a contact operator χ 2 s V µ V µ , which was not studied before.

JCAP03(2015)052
In section 3, we systematically analyzed the GDM thermal relic density. For the viable parameter space, we found that the Higgs-exchange contributions dominate the dark matter annihilation cross sections, where only the dark matter mass M χ and coupling product ξ h ξ s are relevant. Since the leading order interactions between GDM and SM fields are proportional to the corresponding SM particle masses, the DM annihilations into heavy mode dominates in large M χ range. In figure 3, the compatibility of the predicted GDM thermal relic abundance with Planck data was demonstrated in the M χ − |ξ h ξ s | plane for a wide range of χ s mass. The red solid curve of figure 3 corresponds to the central value of Ω χ0 h 2 0.12. Since the GDM cross sections are proportional to ( |ξ h ξ s |) 4 , the predicted parameter space of |ξ h ξ s | in figure 3 has little sensitivity to the experimental uncertainty of Ω χ0 h 2 .
In section 4, we further studied possible direct and indirect detections of GDM, as well as discussing its collider searches. Direct detection of GDM relies on its effective interactions with light fermions. In comparison with the GDM annihilation cross sections, the GDMnucleon scattering is largely suppressed for the small momentum exchange. As shown in figure 4, the required range of |ξ h ξ s | for accommodating the thermal relic abundance predicts signals to be even lower than the general neutrino background, so it is out of reach of the current direct detection technique. For indirect detections, we mainly studied constraints from the observation of gamma ray spectrum which is most promising for searching the GDM. The line spectrum arises from direct annihilations of dark matter into γ's. These operators are generated at one-loop level for the GDM. The diffuse continuum spectrum reflects the secondary photons produced from primary dark matter annihilations into massive gauge bosons, quarks or leptons. We summarized the constraints for these processes in the GDM parameter space, as shown in figure 5. In contrast to the direct detection, we found that gamma ray searches are promising, and have higher sensitivity to the heavier GDM particles. For M χ O(100)GeV, the prediction of our model is within the reach of future gamma ray searches of diffused spectrum. For collider searches, we first studied the constraint from measuring Higgs invisible decays at the LHC. We derived a nontrivial upper bound on |ξ h ξ s | for low M χ region, which excludes the GDM with mass M χ < 48GeV at 95% C.L. Finally, we discussed the searches of χ s at hadron colliders with different production channels. The production of GDM particles in association with top pair (bottom pair) or from the vector-boson-fusions (figure 7) can be probed at the upcoming LHC runs and the future high energy pp colliders [52] 10 .
Note added at the proofreading stage. To probe the scalar-type interactions between DM and fermions, the bottom quark associated production can become important. The latest ATLAS search [55] analyzed the DM production in association with a single bottom quark [cf. figure 8(a)], in addition to the DM production associated with a pair of top or bottom quarks [cf. figure 7(a)]. It was found that the top associated DM production has higher sensitivity than the bottom associated DM production for scalar-type operators, where the bottom associated production has larger phase space for the final state, but not enough to compensate the suppression effect of bottom Yukawa coupling (relative to the top Yukawa coupling). The lower bound on the cutoff scale for the scalar DM is around O(10) GeV, which is generally weaker than that of the Dirac DM [51]. Another way to probe the interactions of DM with weak gauge bosons is to use the vector boson associated production, as depicted in figure 8(b), which contributes to mono-W/Z signals. In the limit q 2 m 2 φ , our tree-level effective interaction (2.20) will induce a contact operator χ 2 s V µ V µ , which enters figure 8 and was not studied before. 12 Hence, it is motivated to further probe the GDM particles via production processes of figure 8 at the upcoming LHC runs and future high energy pp colliders.

A Formulas for radiative loop factors
In this appendix, we summarize the exact expressions of the loop factors [57] for effective interactions between the dark matter and gauge bosons. The loop factors (C g , C γ ) in eq. (2.23) contain only a single mass-ratio τ j = E 2 /4m 2 j via functions (A V , A F ), where E denotes the center of mass energy. The functions (A V , A F ) are defined as follows, Ref. [56] used the ATLAS analysis of W + / E T to derive constraints on DM-vector-boson effective operators. But it focused on the gauge-invariant operator χ 2 s F µνa i F a µνi , which could be generated only at one-loop in our GDM model and thus is negligible here.

JCAP03(2015)052
Since Z is massive, the loop factor C γz in eq. (2.23) involves another mass ratio η j = m 2 Z /4m 2 j . The functions (B V , B F ) are defined as follows, where we have defined (s W , c W ) ≡ (sin θ W , cos θ W ), and t W ≡ tan θ W , with θ W denoting the weak mixing angle. Also, N C = 3(1) corresponds to the color factor of quarks (leptons).

B Threshold and resonance effects in thermal relic density analysis
In this appendix, we present the calculation of thermal relic density by including the threshold and resonance effects. Given the mass-spectrum of SM particles, we note that the two effects take place in difference mass-ranges and can be treated separately. We first consider the threshold effect. For a generic DM annihilation process χ s χ s → f j f j , the zero-temperature cross section can be parameterized as (σ A v) = (a + bv 2 )v n j , (B.1) where v is relative velocity of two dark matter particles, and v j is the final state velocity from phase space integration. The parameters a and b represent s-wave and p-wave contributions, respectively. For scalar dark matter, we have n = 1(3) for bosonic (fermionic) final states. Under non-relativistic approximation for the DM, we derive where z ≡ m j /M χ and µ 2 + ≡ (1−z 2 )/z 2 . Around freeze-out temperature T f , the DM particle is non-relativistic and the thermal average cross section can be derived by integrating over the relative velocity, where x ≡ M χ /T . For cold dark matter, we have x f 1. Substituting parametrization (B.1) into (B.3), we deduce the approximate thermal averaged cross section for M χ m j , In the kinematically forbidden case of M χ < m j , the nonzero velocity v in (B.1) could make the annihilation viable, which sets a lower bound of v in the thermal integration. For M χ < m j , we define µ 2 − ≡ −µ 2 + > 0 and v 2 > 4µ 2 − . Imposing v > 2µ − in (B.3) and making change of variables, we derive the cross section for M χ < m j , Note that eqs. (B.4) and (B.5) agree at the threshold, i.e., z 1 and µ + µ − 0. Since the gravity-induced interaction only generates s-wave contribution at leading order, we will set b = 0 afterwards and focus on the a term in eq. (B.1) for the following discussion. For each channel, we may infer a from cross section in (3.3a) divided by v n j at v = 0, i.e., (1 − z 2 ) n/2 . Then, the s-wave thermal averaged cross section can be parameterized as In figure 9(a), we depict (I A , I F ) as functions of z for different x and n. Red and blue curves denote the cases with bosonic and fermionic final states, respectively. The (dotted, solid, dashed) curves correspond to x = (∞, 25,10). It is clear that the higher temperature (i.e., smaller x) leads to more enhanced thermal integral from the threshold effect. Comparing the red and blues curves, we see that the annihilation cross section with bosonic final states is more enhanced. Next, we consider the resonance effect, which is important around M χ ∼ m φ /2. From (3.3a), we see a common factor from Higgs-exchange for both ff and V V final states. Taking into account the finite temperature effect, we have s expressions are modified in following way, where u ≡ 4M 2 χ /m 2 φ and ≡ Γ φ /m φ . For non-relativistic GDM, the thermal averaged integration over the propagator factor (B.7) defines the following function, As shown in ref. [24], for narrow resonance like the SM Higgs boson, i.e., ∼ 10 −5 , any expansion over v 2 may yield considerable error around the resonance pole. Thus, we will perform thermal integration numerically for computing σ A v . For M χ > m φ /2, i.e., u > 1, the width effect quickly becomes subdominant and negligible. In the light mass range, M χ m φ /2, the cross section is more enhanced due to finite temperature integration. Also, the Higgs invisible decay starts to open and the width depends on M χ and ξ h ξ s . Figure 9(b) depicts K R as a function of √ u = 2M χ /m φ with = Γ φ /m φ . For illustration, we choose |ξ h ξ s | = 10 15 , 10 16 as two benchmarks. The dotted curve corresponds to x = ∞ , where little difference can be seen for the two cases. The (red, blue) solid curves represent |ξ h ξ s | = (10 15 , 10 16 ) at x = 25 , respectively. For the case of |ξ h ξ s | = 10 15 , we see significant resonance enhancement for √ u 1, as compared with the zero-temperature estimate. The other case of |ξ h ξ s | = 10 16 corresponds to a much larger Higgs width and thus the ratio eq. (4.4)]. As shown in figure 9(b), the resonance effect is much smaller in this case.