Cosmological Signatures of Superheavy Dark Matter

We discuss two possible scenarios, namely the curvaton mechanism and the dark matter density modulation, where non-Gaussianity signals of superheavy dark matter produced by gravity can be enhanced and observed. In both scenarios, superheavy dark matter couples to an additional light field as a mediator. In the case of derivative coupling, the resulting non-Gaussianities induced by the light field can be large, which can provide inflationary evidences for these superheavy dark matter scenarios.


Introduction
The nature of dark matter (DM) is a frontier in the modern cosmology. There have been plenty of astronomical and cosmological evidences for DM, such as the galaxy rotation curve, virial velocities of galaxy clusters, gravitational lensing, bullet clusters, supernovae, cosmic microwave background, existence of galaxies in lifetime of the universe and existence of galaxies on scale of milky way.
The production mechanism of the dark matter is not known so far. It can be produced during inflation or subsequent cosmological evolutions. If the dark matter is produced during inflation, it is hopeful to use inflation as an avenue to probe the nature of the dark matter. There are multiple possible mechanisms. One of them is that gravitational production [1] during inflation. Relevant models include Planckian Interacting Dark Matter (PIDM) [2][3][4][5][6], WIMPZILLA [7,8], SUPERWIMP [9], FIMP [10] and so on. They can be scalar [11,12], vector, fermion or spin-2 particles [13] from different types of beyond standard model physics. Such dark matter candidates are usually hard to probe using collider experiments due to the large mass and small coupling (the gravitationally produced dark matter even do not have electroweak interactions with standard model particles).
It is thus interesting to study the dark matter properties, such as mass, width, spin, parity and couplings from cosmology, using the method of cosmological collider physics [14][15][16][17], which has attracted much attention recently . From the frequency of the oscillation on the squeezed limit non-Gaussianity, one can directly obtain information about the mass of extra fields during inflation. We investigated the search for the dark matter mass in the context of cosmological collider physics in [52], targeting at the class known as superheavy dark matter (SHDM) [68][69][70][71][72] with mass m σ ≥ H. Since the dark matter is produced gravitationally, the signal is subject to a suppression of power (H/M pl ), which makes it hardly observable in the near future experiments.
In this paper, we would like to investigate scenarios with signals within the reach of the near future experiments on primordial non-Gaussianities, such as CMB-S4 [73], Simons Observatory [74], DESI [75], EUCLID [76], SPHEREx [77] and LSST [78]. The idea utilizes the curvaton scenario [79][80][81], where more than one light field is present during inflation. The curvatons are subdominant in the energy density during inflation. In the subsequent evolutions, they are the main source of the primordial curvature perturbations. The dark matter field, on the other hand, is another type of field which has mass m σ ≥ H. As we will see, this scenario can give larger non-Gaussianities (NG's) on the primordial curvature perturbations. Such possibility of large NG from the curvaton scenario has also been studied in [64,82]. Another possibility to imprint the dark matter information through a light field is modulated production of dark matter, where the dark matter mass, and thus the production rate, is modulated by a light scalar. This is related to the mechanism of modulated reheating [83,84].
This paper is organized as follows: in Sec. 2, we setup the model which consists of an inflaton, a massive scalar field and a light field. In Sec. 3, we discuss the production of dark matter during inflation. In Sec. 4, we calculate the primordial three-point and four-point correlation functions of the light field. In Sec. 5 and Sec. 6, we discuss the post-inflationary evolutions of the light field and their observational constraints. More specifically, Sec. 5 devotes to the case where the power spectrum is mostly generated by χ. Sec. 6 devotes to the case where the power spectrum is mostly generated by the inflaton field. We give a brief summary in the end in Sec. 7.

Model
We consider a cosmology model with an inflaton φ, a heavy field σ with mass m σ H inf , and a light field χ with mass m χ m φ with an action Both the inflaton field φ and the curvaton χ are light and have rolling background. The heavy field σ does not have a VEV. Also, the inflaton does not couple to either the heavy field σ or the light field χ while all fields couples to gravity minimally for simplicity. We can decompose the light field χ into a time-dependent background and perturbations as χ(t, x) =χ(t) + δχ(t, x). We further assume that the background of light field varies slowly during inflation, such that it approximates to a constant during inflationχ(t) χ. Aside from these requirements, the form of χ potential V (χ) is left free as long as ρ χ becomes insignificant compared to other components in the late time universe. We discuss two possibilities for the interaction term S int , namely direct coupling and derivative coupling. In the case of direct coupling, the first term and the third term contribute to the effective mass of the heavy field σ, the second term and third term are interaction terms between the heavy field and χ.
On the other hand, it is also well motivated to adopt the derivative σ − χ coupling of the form 1 4 λ 2 σ 2χ2 , leading to interaction the first term and the third term contribute to the effective mass of the heavy field σ, the second term and third term are interaction terms between the heavy field and χ.
For the χ field, the backgroundχ is governed by the following equation of motion which holds during the post inflation era as well. Admitting the χ equation of motion Eq. 2.5 and imposingχ ∼ 0, we can get the χ change rate during the primordial era as: The inflaton is the field that drives inflation. In the rest of this paper, we assume that the change of energy density contributed from χ subdominates during inflation, i.e. χ is a spectator field during inflation. In other words, the following constraint is satisfieḋ where ρ χinf and ρ φinf denote the energy density of the light field χ and the inflaton φ during inflation, respectively. The observable of interest is the primordial three-point correlation or four-point correlation of the primordial curvature perturbation ζ, defined by the following metric written in the ζ gauge where there is no energy density fluctuation, (2.8)

Superheavy DM Gravitational Production During Inflation
From our model (2.1) it is obvious that the heavy field σ is protected by the Z 2 symmetry and only interacts with the χ field. As discussed before, without further interactions σ is then a massive stable particle and interacting with other fields via the small coupling and χ as the mediator, making itself a good DM candidate. For its gravitational production during inflation, we will mostly follow the calculations in [52]. During inflation and in the presence of a background χ field ∼χ, the equation of motion for σ isσ where m σeff is given by (2.2). During the inflationary era, the field σ can be expressed in the standard formalism (see also [85]): where a k and a † −k are the annihilation and creation operators that satisfy the commutation relations [a k , a k ] = 0 and [a k , a . We then find Since we assume minimal coupling of σ, the ξR term is zero in the definition of µ. During slow roll inflation, the scale factor is approximately a(t) = e H inf t = −1/(H inf τ ), and R = 12H 2 inf andḢ inf vanishes. Hence we refer H inf as H in the following discussions. When inflation ends, the scale factor starts to evolve in a non-accelerated way. There are two classes of solutions to the massive field equation of motion, corresponding to "in" state and "out" state, respectively [14]: where H (1) ν (x) is the Hankel function of the first kind, and J ν (x) is the Bessel function. The mode functions are related via a Bogoliubov transformation as Inserting the explicit expressions for the "in" state mode function and "out" state mode function, we obtain the following expressions for the Bogoliubov coefficients with the outgoing amplitude |β k | 2 = 1/(e 2πµ − 1). The comoving number density of σ DM produced in the pure slow roll approximation then reads which gives the total number of particles from the past infinity to future infinity, as shown below (see also [86]). Since particles are produced when |ω k /ω 2 k | takes its maximum and thus τ ∼ −µ/k for each mode. The k integral is then rewritten as: The corresponding physical number density of particles at each moment can be evaluated as: .
which is constant in time. Assuming matter domination before reheating and radiation domination afterwards, the current DM abundance can be expressed numerically [52,72]: in the unit of GeV.
Here Ω R and T 0 are the radiation energy fraction and temperature today. To produce non-negligible DM comparable to the observed value Ω CDM h 2 0.12, from the form of (3.10) it is clear that the model disfavors large µ and low H or T RH . In the above slow roll approximation the inflation actually never ends, however, the DM density's exponential dependence on µ and thus on λχ 2 would be crucial for generating NG in the DM sector. To better estimate the σ relic density, we consider the case where inflation is smoothly connected to Minkowski spacetime. The particle production in this scenario would approximate the particle production in more realistic scenarios where inflation is connected to a stage of the universe with much lower Hubble scale (and thus approximately Minkowski), for example, radiation-dominated universe with transition time of order H −1 . By introducing the Stokes line method to compute the particle production beyond exponential precision, one obtains the relic density of σ today as which is also in the unit of GeV. The result is similar to (3.10) but with different numerical prefactors and different powers of µ. For moderate µ, the dependence of n σ shall be dominated by the exponent factor. There are other ways to create σ through gravity, such as a sudden transition from the inflation phase to the radiation domination phase, which may be the case for models such as brane inflation [87][88][89][90][91] or quintessential inflation [92]. The σ number density in this case can be much larger than the smooth transition case. However, the produced DM number density is proportional to H 3 and hence shows little correlation withχ. Without the µ andχ dependence, ρ σ fluctuation cannot carry information of χ field. Even though, the sudden transition mechanism might be useful in the curvaton case, which we will discuss in Sec. 5.
Since the precise Ω σ h 2 is dependent on models of reheating and requires numerical evaluations much more complicated than the semi-analytical approaches we adopt here, it is convenient to parameterize Ω σ h 2 for later use: where A is an constant between O(10 4−9 ) characterizing the DM production efficiency of different models and α describes the remnant power dependence on µ.
In the above discussions we assume that σ particle number is fixed once it is produced near the end of inflation era, one may concern if σ annihilates or produced during reheating and radiation domination. This is indeed a good approximation as long as λ is small and χ only interacts with other fields weakly. Starting from possible σ annihilation. For a weakly interacting χ, the process 2σ → 2χ dominates the annihilation, with s-wave thermal averaged σv ∼ λ 2 /m 2 σ . For heavy DM like σ, when the radiation temperature T m σ , the physical number density n σ n eq σ T 3 , hence annihilation cannot reduce n σ . The potential constraint on λ and other parameters will be discussed below. When T m σ and thus n eq σ T 3 e −mσ/T n σ , the annihilation can happen, but already freeze out n σ σv H. The minuscule annihilation rate also ensures σ DM is safe against indirect detection bounds, as current indirect detection constraints based on DM annihilation are unable to put constraints on weakly interacting DM ( σv 10 −26 cm 3 /s) with mass heavier than O(10) GeV [93].
On the contrary, during the reheating or early stage of radiation domination, the thermal bath can create more σ via the mediation of χ. As a result, the constraint on thermal interaction put an upper limit for λ if we want to ensure that most σ are created by gravity and henceforth δρ σ can carry information of χ fluctuation rather than that of radiation. In the minimal case, radiation are SM like and massless, created by inflaton decays and interact with χ, allowing the later to decay eventually. The χ−radiation interaction coupling (assuming massless fermions) is then ∼ O( Γ χ /m χ ) 1. The mass radiation then couples to σ in the presence ofχ via the effective coupling: As long as λ σr 1, the thermal production of σ would be suppressed, which is usually satisfied according to our numerical results. Moreover, there are other scenarios where the thermal production of σ is further suppressed 1 . For simplicity, here we assume that the constraint of λ σr 1 is always satisfied and gravity dominates σ production.

Isocurvature Fluctuations During Inflation
In this section, we consider the NG of δχ during inflation. We consider bispectrum in Section 4.1 and trispectrum in Section 4.2. We mainly focus on the clock signals which we can use to probe the mass of the superheavy dark matter.

Bispectrum
The second order action of the primordial curvature perturbation can be written down following the procedure in [94,95] Quantizing it in the following way where c † k , c k are the creation and annihilation operators satisfying the usual commutation . The mode function satisfies the following equation of motionζ To the lowest order in slow roll parameter, the solution is (4.4) Using the Schwinger-Keldysh formalism, the Feynman diagram on the left hand side of Figure 1 is evaluated as where the subscript (2) denotes the number of interaction vertices. This type of integral in the squeezed limit k 1 k 2 k 3 is well-known [17,49,52,59,96] and the details are collected in Appendix. A. Eq. (4.5) is then evaluated as where The right hand side of Figure 1 is only differed from the Feynman diagram on the left by a factor of 2λχ 2 /H 2 . So the total contribution to the bispectrum is We can replace the interaction vertices in Figure 1 to derivative couplings, resulting in a similar integration: where the subscript (2) denotes the number of interaction vertices. Similar to the calculations in the direct coupling case, in the squeezed limit k 1 k 2 k 3 , (4.9) and combine it with the integration of another diagram as: where the loop suppression factor for derivative coupling

Trispectrum
To calculate the trispectrum we consider the leading contributions from the 4-point diagrams plotted in Figure 2. The contribution from the first diagram reads: where we have defined k I = k 1 + k 2 . In the collapsed limit k 1 k 2 k 3 k 4 k I , it is where (4.14) Second and third Feynman diagram of Figure 2 is differed from the Feynman diagram on the leftmost by a factor of 2λχ 2 /H 2 and 2λ 2χ4 /H 4 , respectively. So the total contribution to the bispectrum is Based on the same method, the first diagram in the derivative coupling case can be evaluated as where we have defined k I = k 1 + k 2 . In the collapsed limit k 1 k 2 k 3 k 4 k I and summing up all 3 leading diagrams, the result reads: where the loop suppression factor for the derivative case

Curvaton Scenario: NG Introduced by Light Field Decay
In this section and the next section, we discuss the post-inflationary evolution of the isocurvature fluctuations. The discussion here is aligned with the so-called "curvaton scenario" [79][80][81]. The idea dates back to the earlier works [97,98]. After inflation, the subsequent evolution consists of three stages: curvaton field starting oscillation, oscillation persisting for many Hubble times and curvaton decaying. The curvaton field can either dominate the energy density just before the decay or do not dominate during the whole process. Defining the ratio of the curvaton energy density ρ χ and the radiation energy density ρ r to be r dec at the time of curvaton decay, we have Another parameter helpful in characterizing different scenarios is the ratio of the power spectrum of the primordial curvature perturbations generated by the χ field to that generated by the inflaton, which we denoted as In this section we focus on the scenario where the power spectrum is mostly generated by χ (R 1) and leave the scenario where the power spectrum is mostly generated by inflaton to the next section. The current measured value of primordial scalar power spectrum P (0) ζ = 2.1 × 10 −9 [99] demands thatχ has a fixed ratio with H: We will simply fix this ratio in the rest of this section. Also, since by definition the curvaton would decay to radiation, it naturally satisfies the constraint that ρ χ is negligible in the late time universe. We take the minimal form that V (χ) = mχχ 2 2 here. The post-inflationary evolution starts with a period when the curvaton field starts to oscillate around the minimum of the potential. During this time, the universe is dominated by radiation. The oscillation starts when the Hubble scale coincide with the mass of the light field H ∼ m.
The energy density of the χ field is Based on the inhomogeneity we can define the density contrast and curvature perturbations as: Observation requires that [100]χ * H * for curvature perturbations where χ * and H * denotes the field value of χ and Hubble are evaluated at horizon crossing.
The curvature perturbation in the radiation section is simply given by: The curvature perturbation is thus The energy density of the χ field can either dominate or sub-dominate the energy density before the decay. In the following, we discuss these two cases separately. If χ field dominates the energy density before the decay (r ∼ 1), (5.8) becomes χ field sub-dominates the energy density before the decay. On the other hand, the case that the energy of the curvaton is subdominant compared with radiation r 1, we will obtain ζ r 4 δ . (5.10) In this work we will focus on the first scenario where r ∼ 1, as for the later case the calculations and results are qualitatively the same and the fact that the later case is strongly constrained by data.

NG Signals
The χ bispectrum and trispectrum are related to ζ trispectrum and bispectrum in the following way We define the shape function to be (5.14) So the shape functions are in the direct coupling case: Where the normalization factor is the primordial scalar power: C bi (µ) and C tri (µ) are dimensionless functions of µ, as defined in (4.11) and (4.18). If we adopt the derivative coupling, following a similar approach the shape functions becomes: In Fig. 3 we show two typical numerical parameter space for this scenario. The parameter choice are labeled on each plot. We simply take r = 1 to fulfill our expectation that r ∼ O(1). In fact, the size of r doesn't affect the curvature NG clock signal significantly as it only depends on r 3 rather than exponentially. In both cases R is sufficiently larger than 1, consistent with our assumption for curvaton case. In most cases, both the size of bispectrum S and scaled DM relic density f σ are sensitive to µ. The model thus prefers smaller µ that produces enough σ DM and large S.

Observational Constraints
For the curvaton-like case, the power spectrum is dominated by curvaton decays. For formally one can define the ratio According to previous studies [101], the curvaton itself can introduce a non-trivial f NL : As long as R 1 as preferred by data, the observational constraint that f NL < O(10) then requires r to be of O(1).
Another important constraint comes from matter-photon isocurvature, which depends on the perturbation differences between the rest of CDM (c, other than σ), radiation (γ) and baryon (b). We define the gauge invariant entropy fluctuations S xγ ≡ 3(ζ x − ζ γ ), where x denotes b or c. If σ is the dominant component of CDM, since it is created before curvaton decays, its perturbation imprints the information P φ instead of P χ . In this model, the production of DM may even correlate to δχ negatively asχ raises DM effective mass, which may make things even worse. However, as ∂ ln ρσ ∂χ 1 we can take only the leading contribution to isocurvature, which is not related to λ.
The current measurement strongly disfavor cases that all CDM are created before curvaton χ decay [99,102] as this creates a large S mγ /ζ O(1). Therefore in the curvaton scenario we only consider if σ is a subdominant part of CDM, while the rest of CDM is created after or by χ decay: Aside from when CDM generation, baryon number B creation can also important, assuming they are created in the radiation, if it is created before/by/after χ decay makes difference. On the other hand, the lepton number L, as long as it is not created by curvaton χ decay, makes little difference in our calculation (R ν = 0). For more details, see [102]. Following previous discussion, definition of the CDM fraction f σ , and calculation in [102], one can calculate S mγ different cases. For example, if both the rest of CDM and B are created after χ decay (b af ter , c af ter ): Similarly, we have the case b by , c af ter : and b by , c by Notice that this immediately rules out the case that f σ is close to 1 for all different scenarios. Together with current data constraints |S mγ /ζ| O(0.1), we can work out the possible range of r and f σ based on the arguments is shown in Fig. 4.
We also need to take the tri-spectrum τ N L constraint into consideration. According to [101]: .
The current constraint on g NL is of O(10 4 ) [103] and τ N L < 2800 [104]. With the assumption that R 1, the constraints from the trispectrum measurements cannot provide meaningful information about this scenario as long as r is of O(1).

DM modulation Scenario: NG Introduced by DM Production
In this section, we will turn to another limit that χ barely leaves any curvature perturbation observable, in contrast to the curvaton case. Or equivalently we are discussing the R → 0, r → 0 limit using the notation we use in the curvaton case. As inflaton dominates both the curvature perturbation power spectrum and the energy density when χ decays, the NG imprints left by χ will be suppressed by r n and would not introduce visible f NL constraints according to (5.20). Nevertheless, in our model where χ directly couples to the dark matter σ, the information of χ primordial perturbation will also introduce isocurvature fluctuations in DM, which may become significant. The fluctuation in χ during inflation creates NG in terms matter-radiation isocurvature. Certainly, such a model must also be constrained by current observation of isocurvature power spectrum. In the follow section we will show that in order to get a large NG effect in terms of isocurvature, f σ need to be larger and eventually taken to be 1.

Review of this Scenario
According to [52] and discussions in Sec. 3, the heavy DM density produced by gravity is proportional to µ α e −2πµ . For superhorizon modes, the DM relic density before χ decays will pick up an non-zero fluctuation term due to the contribution from χ fluctuation: since during this era we have ρ σ ρ χ ρ r and ζ r ζ φ . Then, up to the leading order, we can expand the expression and get the following expression.
Note that inflaton also contributes to the curvature perturbation, however, it gives subdominant contributions. For more details, see [105].
Since ζ φ and ζ χ/σ is uncorrelated, the main contributor to the NG of σ density mostly stems from the isocurvature mode. This means the three-point function ∂ ln ρ σ ∂χ and for derivative coupling the above relation changes to: where in the last relation we assume α 0 in (3.12) for simplicity. The term grabs a negative sign compared to the χ induced clock signal. The calculation of the clock signal is then similar to the case in Sec. 5. However, to correctly compare the NG in isocurvature with the Gaussian part, the normalization factor changed to primordial isocurvature perturbation at large scale P (0) I deduced from the same mechanism instead. The current experiment only gives an upper limit of the size. In this mechanism, the primordial isocurvature perturbation can be introduced by modulated σ dark matter production similar to in, which will be discussed below.
Now we turn to estimate the primordial isocurvature perturbation power spectrum introduced by χ-σ interaction at large scales, assuming no other source of isocurvature. Similar to the curvaton case but with much smaller R and r, the ζ γ mainly follows inflaton fluctuation: while ζ c and ζ b is assumed to be created by radiation since χ vanishes. This is corresponding to the c after , b after scenario, leaving Resulting in the total matter non-adiabatic perturbation In the last relation we use the approximation that δχ ≈ H/2π and the fact that Ω b is only a small fraction of Ω m . According to the current isocurvature constraints [106], the value of λχ/Hµ must be small enough ( O(10 −5 )). Similarly, we can write down the S mγ for the derivative coupling case: which is solvable once we impose the relation betweenχ andχ such as in Eq. 2.5. The constraint from primordial isocurvature is shown in Figure 5 as the gray shades. It is obvious from the form of Eq. 6.7 that the dependence onχ 2 /χH ∝ m 4 χ is small since χ is a light field. With a reasonable choice of m χ , such as 0.01H we chose here, the constraint from primordial isocurvature are much weaker, allowing larger NG signals.
Since in the modulator scenario the light field χ no longer generates the primordial curvature spectrum, the ratioχ/H is no longer fixed and hence becomes a free parameter. In this work we are particularly interested in the region whereχ/H is larger than 1 to use Eq. 2.5 directly. Therefore the parameter space we shall have an additional dimension compared to the curvaton case. If we fix µ to be 1 which is required to have enough NG, we can scan through the parameter space. We show the result of scans in both direct coupling and derivative coupling in Fig. 5. In both panels, the whole allowed regions can provide enough σ DM (f σ = 1) with a proper choice of parameter A between 10 4 and 10 9 . For direct coupling case it would be very hard to give a large NG signal without extreme fine tunings. The reason is due to the fact that only 2 parameters: the coupling λ and the ratioχ/H controls the strength of NG signals. Instead, for the derivative coupling case, a much larger NG signal can be achieved without much fine tuning (µ 0.1).

Conclusion and Outlook
We propose a mechanism that a light spectator field can leave nontrivial observation signals via its interaction between gravitationally produced heavy DM. The minimal model consists of an inflaton φ, a massive field σ and a light field χ. The conclusion is largely inflation model independent as all non-trivial interactions happens inside the σ-χ sector and the inflaton φ only acts as a background field. Specifically we consider two types of wellmotivated couplings, direct coupling and derivative one. For both type of couplings the primordial three-point and four-point correlation functions of the light field χ mediated by the heavy DM field σ are calculated. The Z 2 symmetry that protects σ demands that such interactions can only present at loop level, indicating suppressed NG signals. However, large enough signal can still appear with proper m σ and derivative couplings.
The first case we propose is the well-known curvaton scenario, where the curvature perturbation we observe today are dominated by χ. The current matter isocurvature observations strongly constraints the σ heavy DM relic density. Therefore in this scenario we consider the curvature NG introduced by χ after it decays to radiation. The curvatonheavy DM derivative interaction still leads to measurable clock signals while the direct coupling cannot provide enough NG.
In another scenario in which χ energy density and power spectrum is always subdominant, it is difficult to measure χ properties via curvature perturbations. However, the σ heavy DM production is sensitive to χ field values during the inflationary era. Consequently, the matter isocurvature mode keeps the imprints of the primordial χ field configuration. The result of both types of couplings are shown. We will leave other intermediate cases to future studies.
In this work, we have focused on relatively simple cases. There can be many possibilities to further enhance the signal, including symmetry breaking, chemical potential, multiple dark matter particles with different masses, and so on. It would be interesting to explore these possibilities.