Probing superheavy dark matter with gravitational waves

We study the superheavy dark matter (DM) scenario in an extended $B-L$ model, where one generation of right-handed neutrino $\nu_R$ is the DM candidate. If there is a new lighter sterile neutrino that co-annihilate with the DM candidate, then the annihilation rate is exponentially enhanced, allowing a DM mass much heavier than the Griest-Kamionkowski bound ($\sim10^5$ GeV). We demonstrate that a DM mass $M_{\nu_R}\gtrsim10^{13}$ GeV can be achieved. Although beyond the scale of any traditional DM searching strategy, this scenario is testable via gravitational waves (GWs) emitted by the cosmic strings from the $U(1)_{B-L}$ breaking. Quantitative calculations show that the DM mass $\mathcal{O}(10^9-10^{13}~{\rm GeV})$ can be probed by future GW detectors.


Introduction
The freeze-out of weakly interacting massive particles (WIMPs) [1] has been the most popular explanation for the particle origin of dark matter (DM) for decades. In this paradigm, the DM relic abundance is determined by [2][3][4] (1.1) with σ ann. being the pair annihilation cross section of DM particles to the Standard Model (SM) particles, v the relative velocity, and ... the thermal average. In the second approximate equality we have used σ ann. v ∼ α 2 DM /M 2 DM where α DM is the finite structure constant of the coupling between the dark and SM sectors. Eq. (1.1) shows that if the dark matter is at electroweak (EW) scale and its coupling is of the order of the EW coupling, then the freeze-out mechanism can explain the observed DM density Ω DM h 2 ≈ 0.12 [5,6]. This is the so-called "WIMP miracle", which has motivated enormous efforts to search for EW scale WIMPs via direct [7], indirect [8] and collider [9] experiments.
However, as proposed in Refs. [34,35], DM mass beyond the GK bound is also possible within the thermal freeze-out framework, as long as there is a lighter unstable species that co-annihilates with the DM candidate (dubbed as the "zombie collision" [35]) and exponentially enhances the interaction rate. In other words, for the same coupling, such an annihilation allows an exponentially heavier DM mass compared to the normal DM pair annihilation scenario. The extra entropy produced from the late time decay of the lighter species further dilutes the DM density, leading to an even higher DM mass upper limit that can reach 10 16 GeV [34].
While the above zombie annihilation mechanism is theoretically appealing, it is very challenging to probe such superheavy DM via the traditional direct, indirect or collider experiments. In this article, we propose a zombie annihilation mechanism that is associated with the breaking of a U (1) symmetry, which leads to the formation of cosmic strings that can be detected via the gravitational wave (GW) signals at current or future GW detectors. As a benchmark, we study an extended B − L model, where one generation of the right-handed neutrino (RHN) is the DM candidate, and a new Dirac sterile neutrino serves as the lighter species for co-annihilation. Section 2 introduces the model, while Section 3 is devoted to the calculation of thermal freeze-out and DM relic abundance. We then investigate the corresponding cosmic strings and GW signals in Section 4, where the correlation to the DM scenario is obtained, and the recent NANOGrav excess is also commented. The conclusion will be given in Section 5.

The extended B − L model
We begin with the B − L model [36][37][38][39], which gauges the U (1) B−L group, with B and L being the baryon and lepton number, respectively. In the B − L model, three generations of Majorana RHNs ν i R (with B − L = −1) are introduced for gauge anomaly cancellation. The new gauge boson of U (1) B−L is denoted as Z , and a complex scalar field Φ (with B − L = 2) is introduced to break the U (1) B−L and provide mass for Z . The relevant Lagrangian reads where i, j are generation indices and D µ = ∂ µ −ig B−L XZ µ is the gauge covariant derivative with X being the corresponding B −L number. The scalar potential triggers a spontaneous symmetry breaking at |Φ| = v φ / √ 2. If we parametrize Φ as (φ + iη)/ √ 2, then after the U (1) B−L breaking η is absorbed as the longitudinal mode of Z , and the particles obtain the following masses The Lagrangian (2.1) can elegantly explain the extremely small mass of the SM neutrinos and the matter-antimatter asymmetry of the Universe via Type-I seesaw [40] and leptogenesis [41], respectively. Especially, the explanation of neutrino mass ∼ 0.1 eV [5,42] prefers superheavy RHNs, since the seesaw mechanism requires M ν R ∼ λ D λ † D × 10 14 GeV. For the sake of a DM candidate, we adjust Eq. (2.1) by assigning a Z 2 symmetry, under which the third generation RHN is odd while all other particles are even. Consequently, λ i3 D = 0 for i = 1, 2, 3 and λ j3 R = 0 for j = 1, 2, making ν 3 R free from the ν 3 R → H decay and hence can be a DM candidate. This setup is generally adopted in the ν R -DM scenarios [43][44][45][46][47][48]. 2 To realize the zombie co-annihilation, we further extend the model by introducing a Dirac sterile neutrino ψ (with B − L = −1) and a gauge singlet real scalar mediator S. 3 The quantum numbers of the particles beyond the SM (BSM) are summarized in Table 1. Since hereafter we only discuss the DM candidate ν 3 R , for simplicity we will denote ν 3 R as ν R and λ 33 R as λ R . The Lagrangian relevant for DM reads which provides the zombie collisions and their charge conjugations via the exchange of an S mediator. If we name ν R as a "survivor" and ψ as a "zombie", then Eq. (2.4) is infecting a survivor to a zombie, and that is why it is called a zombie collison [35]. As we will see in the next section, for M ψ < M ν R = λ R v φ / √ 2, the annihilation rate is exponentially enhanced and M ν R can reach 10 13 GeV while still generating the correct DM relic abundance.
We work in the parameter space so that the DM decay channels ν R → ψS/ψS or ν R → ψψψ are kinematically forbidden, while the mediator S can decay to pairs ofψν R , ψν R , and ψψ, with the width

6)
2 Two generations of Z2-even RHNs are already sufficient to explain the neutrino oscillation data and realize leptogenesis, see the recent review [49] and the references therein. 3 See Ref. [50] for the low-scale phase transition study in the complex singlet extended U (1)B−L model. and The final essential ingredient of the zombie mechanism is an appropriate decay channel for the Dirac sterile neutrino ψ. After the freeze-out of ν R , ψ needs to decay into the SM particles, otherwise itself is a WIMP DM candidate that satisfies the GK bound, and hence ν R cannot exceed the GK bound due to the M ν R < 3M ψ constraint. Since ψ is odd under the Z 2 while the SM particles are even, such ψ decay must involve some Z 2 -breaking interactions, which can also result in the decay of the DM candidate ν R . To make ν R sufficiently long-lived (with a lifetime longer than 10 27 s, as required by the diffuse gammaray spectrum [51][52][53]), the breaking of Z 2 should be mediated by either extremely small couplings or extremely high scales with some level of fine-tuning. 4 For the former case, we can have an interacting vertex like y ψ¯ L Hψ R with |y ψ | 1; while for the latter case, a dimension-6 operator can be induced by exchanging a new heavy color triplet scalar where i, j and j are generation indices and j = j . The color indices of quarks are implicitly summed up in a completely asymmetric way via the Levi-Civita symbol to yield a color singlet. We will take Eq. (2.8) as an example for the further study, but keep in mind that the zombie DM mechanism holds as long as ψ can decay via a tiny Z 2 -breaking interaction. According to Eq. (2.8), the decay width of ψ is assuming an exclusive uds-quark final state. Note that Eq. (2.8) also allows the decay ν R → ψψ ( * )ψ( * ) → 9j via one or two off-shell ψ's, and hence a high Λ is required to keep ν R sufficiently long-lived. On the other hand, the same high Λ also results in the late time decay of ψ, which can produce entropy to dilute the ν R abundance.
3 The superheavy ν R dark matter

Thermal freeze-out
Assume the reheating temperature after inflation is higher than M ν R and hence both ν R and ψ are originally in thermal equilibrium with each other and the SM particles, and the number densities are described by where α = ν R , ψ orψ, K 2 is the modified Bessel function of the second kind, and the factor 2 is for spin degeneracy. As the Universe expands and cools down, the reaction rate between ν R and ψ (see Eq. (2.4)) drops. When the interaction rate is lower than the Universe expansion rate, ν R deviates from the chemical equilibrium and eventually freeze-out to be the DM candidate. The freeze-out process can be characterized by a set of Boltzmann equations, which are discussed in detail below. Consider the radiation dominated era that the energy and entropy densities are respectively with g * being the number of relativistic degrees of freedom which we take to be the SM value 106.75. The Hubble constant can be solved as Finally, we define the particle abundance as Y α = n α /s, i.e. the ratio of number density to the entropy density, and hence the abundances of the equilibrium distributions are Under above conventions, the Boltzmann equations of ν R and ψ can be expressed as a set of ordinary differential equations, is the freeze-out relic abundance of the DM candidate ν R . Here the γ ab→cd 's are the corresponding interaction rates, γ ab→cd = n eq a n eq b σ ab→cd v , (3.8) and the detailed definitions can be found in Appendix A. Since M ν R > M ψ , the abundance of ψ is higher than ν R , i.e. n eq ψ ∼ (n eq ν R ) Mν R /M ψ , leading to an exponentially enhanced annihilation rate compared with the WIMP scenario. The factor "2" in front of γ ν R ψ→ψψ , γ ν Rψ →ψψ and γ ν R ν R →ψψ comes from charge conjugation (e.g. γ ν Rψ →ψψ ), and we have made use of Y ψ ≡ Yψ due to CP conservation. For simplicity, we assume M φ , M Z M ν R so that the B − L scalar φ and gauge boson Z do not participate in the Boltzmann equations explicitly, but Z contribute to γ ν R ν R →ff and γ ψψ→ff (where f denotes the SM fermions in equilibrium) via the off-shell s-channel diagrams. In principle, the decay ψ → jjj and scattering ψj → jj induced by the dimension-6 operator in Eq. (2.8) also affect the evolution of Y ψ , but the effect is completely negligible due to the large Λ. Instead, the late time ψ decay after ν R freeze-out plays an important role, as discussed in the next subsection.

Decay of ψ and the dark matter relic abundance
Due to its long life time, ψ can be treated as a stable particle during the freeze-out of ν R , and its "relic abundance" is given by . After freeze-out, the energy density of ψ scales as a −3 where a(t) is the Friedmann-Lemaitre-Robertson-Walker (FLRW) scale factor, while the radiation energy density scales as a −4 . As a result, if ψ is sufficiently long-lived, it dominates the Universe energy at the late time, and its decay would generate significant entropy that further suppresses the ν R relic abundance. The dilution factor can be estimated quickly by the following considerations: ψ decays at T = T ψ when Γ ψ ∼ H, i.e.
If ψ decays to jjj very promptly at t ψ = 1/Γ ψ and reheats the Universe up to T ψ , by energy conservation we have Combining these two equations, one obtains the entropy enhancement factor (or equivalently, the DM density dilution factor) 11) and the final DM relic abundance is which will be adopted in our numerical study. Now we are ready to calculate the DM relic abundance by solving the Boltzmann equations (3.6) and (3.7) to get Y ∞ ν R , and combining the dilution factor ∆ ψ . For the numerical study, we fix g B−L = 1, λ R = 0.1, M ν R = 1.9M ψ , M S = 2M ν R and λ 1 = λ 2 , and vary M ν R , λ 1 and Λ to find four benchmark points (BPs) listed in Table 2. All BPs can  Table 2. The BPs for the superheavy RHN DM scenario. g B−L = 1, λ R = 0.1, M ν R = 1.9M ψ , M S = 2M ν R and λ 1 = λ 2 are fixed, and the first three columns M ν R (RHN mass), λ 1 (ν R -ψ-S coupling) and Λ (Z 2 breaking scale) are free parameters. The ψ decay temperature T ψ and reheating temperature T ψ , dilution factor ∆ ψ and the DM life time τ ν R are derived, see the text.  Table 2. Y ν R and Y eq ν R are shown in blue and green curves, respectively. For comparison, we also show in orange curves the Y ν R evolution without the zombie collision. yield the correct relic abundance where H 0 = 100 h km/s/Mpc (with h = 0.674) and s 0 = 2891.2 cm −3 are current Hubble constant and entropy density, respectively [6]. As we can see, due to the enhanced cross section from the light sterile neutrino (i.e. zombie collision) and the dilution factor from ψ decay, a 1.5 × 10 13 GeV RHN can still yield the correct DM abundance with a moderate coupling λ 1 = λ 2 = 1.5. Even for a small coupling λ 1 = λ 2 = 0.30, the RHN DM can be as heavy as 3.0 × 10 9 GeV, well above the GK bound. In our mass setup, ν R can decay to 9 jets via 2 off-shell ψ's, and the life time is calculated by the FeynRules [55] and MadGraph5 aMC@NLO [56] packages. τ ν R 10 27 s is satisfied in all BPs, as required by the diffuse gamma-ray spectrum [51][52][53].
To distinguish and compare the contributions from "zombie collision" and "ψ decay dilution" to the DM relic density, we plot the Y ν R evolution as blue curves in the BPs in Fig. 1. The equilibrium distribution Y eq ν R is shown in green curves for reference, such that one can clearly see the ν R deviates from the thermal equilibrium and freeze-out to a constant abundance Y ∞ ν R . The late time decay of ψ will dilute the freeze-out abundance to Y ∞ ν R /∆ ψ , which is eventually Y DM ≈ 0.8 eV/M ν R , as shown in red dashed straight lines. To manifest the importance of the zombie collision, we also plot the Y ν R evolution assuming λ 1 = λ 2 = 0 in orange curves. In that case, ν R freeze-out in a very high abundance and hence cannot provide the correct DM density even with the help of the ψ decay dilution. By comparing the curves of Y ν R , Y ν R (without zombie) and Y DM , one can see that the contributions from zombie collision and ψ decay dilution are comparable.

Cosmic strings and the gravitational wave signals
The BPs obtained in the last section are for M ν R 10 9 GeV, much higher than the available energy scale of any traditional DM direct, indirect or collider searches, making it almost hopeless to test superheavy DM scenario. However, the recently developed GW astronomy offers an unique opportunity to probe this scenario: as the RHN mass M ν R is associated with a high-scale U (1) B−L breaking, which forms cosmic strings that can generate detectable stochastic GW signals [57][58][59][60][61][62][63].
Cosmic strings are one dimensional topological defects form during a spontaneous symmetry breaking if the topology of the vacuum is not simply connected [64]. In our scenario, we consider Nambu-Goto cosmic strings that can form after the U (1) B−L breaking, and the energy density per unit length µ ∼ v 2 φ . A very important observable for the cosmic strings is the dimensionless combination where G is the Newton's constant of gravitation. After formation, the collisions and self-interactions of strings produce sub-horizon, non-self-interacting string loops, which emit GWs via cusp, kink and kink-kink collisions, and they produce GWs throughout the Universe history. The incoherent superposition of such continuous emission results in today's stochastic GW signals.
According to above physical picture, the GW spectrum today can be expressed as where a(t) is the FLRW scale factor and t 0 is the current cosmic time. n CS ( , t) is the number density of sub-horizon string loops with invariant length at cosmic time t, while P GW (f, ) is the loop power spectrum describing the power of GW with frequency f emitted from a loop with length . We follow the method described in Ref. [59] to calculate the GWs from the Nambu-Goto string [65] by transforming Eq. (4.2) into [59,66] where n = 1, 2, ..., labels the radiation frequencies ω n = 2πn/( /2), and P n is the corresponding average loop power spectrum which we use the numerical results from Ref. [66], and where we have transformed the integral variable from time t to the redshift z.
The current searches for the SGWB of PTAs constrains Gµ < 10 −11 [99,100] which requires the M ν R /λ R < 3.2 × 10 13 GeV according to Eq. (4.1). In Fig. 2, we can see that BP1 (with M ν R = 1.5×10 13 GeV) is already within the reach of the PPTA and NANOGrav experiments, and the null results of these observations suggest that the BP1 is excluded; but a DM with mass slightly lower than BP1 is still allowed and can be tested in the near future. Interestingly, recently the NANOGrav collaboration reported a common-spectrum process based on the 12.5-yr data set [101], which might be a hint for the stochastic GW background from cosmic string [102][103][104]. However, there is some discrepancy between the results of the NANOGrav and PPTA collaborations, and we still have to wait for the future data of the PTA experiments or even the crosscheck from the space-and groundbased detectors to reveal the origin of the excess. The GW signals from BP2 and BP3 can be accessed by a few future detectors, while the GWs from BP4 (with M ν R = 3.0 × 10 9 GeV) can be reached only by BBO and CE. Therefore, the cosmic strings induced GWs can probe our scenario for RHN DM mass between ∼ 10 9 GeV and ∼ 10 13 GeV. We also note that there are still considerable uncertainties in calculating the cosmic string GWs, and the expected reach for DM mass range might vary when more accurate treatments are available. 6

Conclusion
In this article, we have realized the zombie collision DM mechanism in an extended U (1) B−L model, showing it allows RHN DM mass up to 10 13 GeV with a moderate interaction coupling. Such a superheavy DM scenario benefits from both the exponentially annihilation rate due to the lighter sterile neutrino, and the late time entropy production from the sterile neutrino decay.
As the DM mass significantly exceeds the sensitivity regions of the traditional DM detection experiments, we propose to probe the zombie collision DM scenario via GW astronomy by detecting the GW signals from the cosmic strings formed when the U (1) B−L breaking. Calculations show that RHN DM mass between O(10 9 GeV) and O(10 13 GeV) can be probed by future GW experiments; especially, for DM with mass ∼ 10 13 GeV, the signal might already be reached by the recent NANOGrav results, but more data are needed to clarify this.

A The interaction rates
The interaction rates in Section 3 are defined as γ ab→cd ≡ σ ab→cd v n eq a n eq b = where K 1 is the modified Bessel function of the first kind, , and the reduced cross section iŝ where σ ab→cd is the Lorentz invariant cross section, which is a function of s. All initial and final state spin and internal degrees of freedom have been summed up. Our convention of γ ab→cd is the same as Ref. [106].
Since the explicit expressions for the interaction rates in Section 3 are too tedious to be shown in the article, below we only list the amplitude squares for the corresponding processes in our mechanism. With the amplitude squares in hand, one can easily derives the cross sections and interaction rates. Assume the collision to be a(p 1 )b(p 2 ) → c(p 3 )d(p 4 ), and define the Mandelstam variables as s = (p 1 + p 2 ) 2 and t = (p 3 − p 1 ) 2 , for the zombie collisions we have Note that in the process ν R ψ → ψψ, the mediator S can be in the s-channel, thus its width should be included to make the integral in Eq. (A.1) converge.
For the pair annihilation of ν R ν R → ψψ, there are both contributions from the Yukawa interactions and U (1) B−L gauge interactions, while for ν R ν R → ψψ, only Yukawa interactions contribute, Finally, the Z mediated annihilation to the SM fermions reads Here the factor 13/2 comes from the summation of SM fermionic degrees of freedom.