Probing superheavy dark matter with gravitational waves

The freeze-out of weakly interacting massive particles (WIMPs) Lee:1977ua 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 Kolb:1990vq ; Bertone:2004pz ; Lisanti:2016jxe

with $\sigma_{\rm ann.}$ being the pair annihilation cross section of DM particles to the Standard Model (SM) particles, $v$ the relative velocity, and $\left\langle...\right\rangle$ the thermal average. In the second approximate equality we have used $\left\langle\sigma_{\rm ann.}v\right\rangle\sim\alpha_{\rm DM}^{2}/M_{\rm DM}^{2}$ where $\alpha_{\rm DM}$ is the finite structure constant of the coupling between the dark and SM sectors. Eq. (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 $\Omega_{\rm DM}h^{2}\approx 0.12$ Planck:2018vyg ; ParticleDataGroup:2020ssz . This is the so-called “WIMP miracle”, which has motivated enormous efforts to search for EW scale WIMPs via direct Schumann:2019eaa , indirect Gaskins:2016cha and collider Kitano:2010fa experiments.

However, WIMP mass deviating from EW scale is possible. For example if $\alpha_{\rm DM}\gg 0.01$, then $M_{\rm DM}\gg 100$ GeV can yield the correct DM abundance. Since $\alpha_{\rm DM}$ has an upper limit $\sim 4\pi$ set by the unitarity bound, $M_{\rm DM}$ also has an upper limit, which is $\sim 300$ TeV derived from the partial wave analysis, known as the Griest-Kamionkowski (GK) bound Griest:1989wd .¹¹1This bound applies to the elementary particle DM. If the DM is a composite object, another bound $R_{\rm DM}\gtrsim 7.5\times 10^{-7}$ fm applies to the DM size Griest:1989wd . It is known that DM can be heavier than the GK bound in case of non-thermal dynamics, non-standard cosmological history Kolb:1998ki ; Hui:1998dc ; Chung:1998rq ; Chung:2001cb ; Harigaya:2014waa ; Davoudiasl:2015vba ; Randall:2015xza ; Harigaya:2016vda ; Berlin:2016vnh ; Berlin:2016gtr ; Bramante:2017obj ; Hamdan:2017psw ; Cirelli:2018iax ; Babichev:2018mtd ; Hashiba:2018tbu ; Hooper:2019gtx ; Davoudiasl:2019xeb ; Chanda:2019xyl or the first-order cosmic phase transition Baker:2019ndr ; Chway:2019kft ; Marfatia:2020bcs ; Baldes:2020kam ; Azatov:2021ifm . However, as proposed in Refs. Berlin:2017ife ; Kramer:2020sbb , 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” Kramer:2020sbb ) 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 Berlin:2017ife .

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.

We begin with the $B-L$ model Davidson:1978pm ; Marshak:1979fm ; Mohapatra:1980qe ; Davidson:1987mh , 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 $\nu_{R}^{i}$ (with $B-L=-1$) are introduced for gauge anomaly cancellation. The new gauge boson of $U(1)_{B-L}$ is denoted as $Z^{\prime}$, and a complex scalar field $\Phi$ (with $B-L=2$) is introduced to break the $U(1)_{B-L}$ and provide mass for $Z^{\prime}$. The relevant Lagrangian reads

where $i$, $j$ are generation indices and $D_{\mu}=\partial_{\mu}-ig_{B-L}XZ_{\mu}^{\prime}$ is the gauge covariant derivative with $X$ being the corresponding $B-L$ number. The scalar potential triggers a spontaneous symmetry breaking at $\left\langle|\Phi|\right\rangle=v_{\phi}/\sqrt{2}$. If we parametrize $\Phi$ as $(\phi+i\eta)/\sqrt{2}$, then after the $U(1)_{B-L}$ breaking $\eta$ is absorbed as the longitudinal mode of $Z^{\prime}$, and the particles obtain the following masses

The Lagrangian (2) can elegantly explain the extremely small mass of the SM neutrinos and the matter-antimatter asymmetry of the Universe via Type-I seesaw Minkowski:1977sc and leptogenesis Fukugita:1986hr , respectively. Especially, the explanation of neutrino mass $\sim 0.1$ eV Planck:2018vyg ; KATRIN:2019yun prefers superheavy RHNs, since the seesaw mechanism requires $M_{\nu_{R}}\sim\lambda_{D}\lambda_{D}^{\dagger}\times 10^{14}~{}{\rm GeV}$.

For the sake of a DM candidate, we adjust Eq. (2) by assigning a $\mathbb{Z}_{2}$ symmetry, under which the third generation RHN is odd while all other particles are even. Consequently, $\lambda_{D}^{i3}=0$ for $i=1$, 2, 3 and $\lambda_{R}^{j3}=0$ for $j=1$, 2, making $\nu_{R}^{3}$ free from the $\nu_{R}^{3}\to\ell H$ decay and hence can be a DM candidate. This setup is generally adopted in the $\nu_{R}$-DM scenarios Okada:2010wd ; Okada:2012sg ; Okada:2016tci ; Okada:2016gsh ; Okada:2018ktp ; Borah:2018smz .²²2Two generations of $\mathbb{Z}_{2}$-even RHNs are already sufficient to explain the neutrino oscillation data and realize leptogenesis, see the recent review Xing:2020ald and the references therein. To realize the zombie co-annihilation, we further extend the model by introducing a Dirac sterile neutrino $\psi$ (with $B-L=-1$) and a gauge singlet real scalar mediator $S$.³³3See Ref. Bian:2019szo for the low-scale phase transition study in the complex singlet extended $U(1)_{B-L}$ model. The quantum numbers of the particles beyond the SM (BSM) are summarized in Table 1. Since hereafter we only discuss the DM candidate $\nu_{R}^{3}$, for simplicity we will denote $\nu_{R}^{3}$ as $\nu_{R}$ and $\lambda_{R}^{33}$ as $\lambda_{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 $\nu_{R}$ as a “survivor” and $\psi$ as a “zombie”, then Eq. (5) is infecting a survivor to a zombie, and that is why it is called a zombie collison Kramer:2020sbb . As we will see in the next section, for $M_{\psi}<M_{\nu_{R}}=\lambda_{R}v_{\phi}/\sqrt{2}$, the annihilation rate is exponentially enhanced and $M_{\nu_{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 $\nu_{R}\to\psi S/\bar{\psi}S$ or $\nu_{R}\to\psi\psi\bar{\psi}$ are kinematically forbidden, while the mediator $S$ can decay to pairs of $\bar{\psi}\nu_{R},\psi\nu_{R}$, and $\psi\bar{\psi}$, with the width

where

and $a_{R}\equiv M_{\nu_{R}}^{2}/M_{S}^{2}$, $a_{\psi}\equiv M_{\psi}^{2}/M_{S}^{2}$.

The final essential ingredient of the zombie mechanism is an appropriate decay channel for the Dirac sterile neutrino $\psi$. After the freeze-out of $\nu_{R}$, $\psi$ needs to decay into the SM particles, otherwise itself is a WIMP DM candidate that satisfies the GK bound, and hence $\nu_{R}$ cannot exceed the GK bound due to the $M_{\nu_{R}}<3M_{\psi}$ constraint. Since $\psi$ is odd under the $\mathbb{Z}_{2}$ while the SM particles are even, such $\psi$ decay must involve some $\mathbb{Z}_{2}$-breaking interactions, which can also result in the decay of the DM candidate $\nu_{R}$. To make $\nu_{R}$ sufficiently long-lived (with a lifetime longer than $10^{27}$ s, as required by the diffuse gamma-ray spectrum Cirelli:2012ut ; Essig:2013goa ; Blanco:2018esa ), the breaking of $\mathbb{Z}_{2}$ should be mediated by either extremely small couplings or extremely high scales with some level of fine-tuning.⁴⁴4The word “fine-tuning” is also in the sense that we assume only $\psi$ directly participates in the $\mathbb{Z}_{2}$-breaking interactions. For the former case, we can have an interacting vertex like $y_{\psi}\bar{\ell}_{L}H\psi_{R}$ with $|y_{\psi}|\ll 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^{\prime}$ are generation indices and $j\neq j^{\prime}$. 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. (9) as an example for the further study, but keep in mind that the zombie DM mechanism holds as long as $\psi$ can decay via a tiny $\mathbb{Z}_{2}$-breaking interaction. According to Eq. (9), the decay width of $\psi$ is

assuming an exclusive $uds$-quark final state. Note that Eq. (9) also allows the decay $\nu_{R}\to\psi\psi^{(*)}\bar{\psi}^{(*)}\to 9j$ via one or two off-shell $\psi$’s, and hence a high $\Lambda$ is required to keep $\nu_{R}$ sufficiently long-lived. On the other hand, the same high $\Lambda$ also results in the late time decay of $\psi$, which can produce entropy to dilute the $\nu_{R}$ abundance.

The BPs obtained in the last section are for $M_{\nu_{R}}\gtrsim 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_{\nu_{R}}$ is associated with a high-scale $U(1)_{B-L}$ breaking, which forms cosmic strings that can generate detectable stochastic GW signals Buchmuller:2013lra ; Dror:2019syi ; Auclair:2019wcv ; Fornal:2020esl ; Samanta:2020cdk ; Masoud:2021prr ; Buchmuller:2021mbb .

Cosmic strings are one dimensional topological defects form during a spontaneous symmetry breaking if the topology of the vacuum is not simply connected Nielsen:1973cs . 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 $\mu\sim v_{\phi}^{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_{\rm CS}(\ell,t)$ is the number density of sub-horizon string loops with invariant length $\ell$ at cosmic time $t$, while $P_{\rm GW}(f,\ell)$ is the loop power spectrum describing the power of GW with frequency $f$ emitted from a loop with length $\ell$. We follow the method described in Ref. Auclair:2019wcv to calculate the GWs from the Nambu-Goto string Vachaspati:1984gt by transforming Eq. (25) into Auclair:2019wcv ; Blanco-Pillado:2017oxo

where $n=1$, 2, …, labels the radiation frequencies $\omega_{n}=2\pi n/(\ell/2)$, and $P_{n}$ is the corresponding average loop power spectrum which we use the numerical results from Ref. Blanco-Pillado:2017oxo , and

where we have transformed the integral variable from time $t$ to the redshift $z$.

To evaluate the integration in Eq. (27), the cosmic time and Hubble constant should be expressed as functions of redshift

with the abundance Planck:2018vyg

and the function

enfolds the change of relativistic degrees of freedom at $e^{+}e^{-}$ annihilation (200 keV, $z=10^{9}$) and QCD phase transition (200 MeV, $z=2\times 10^{12}$) Binetruy:2012ze . The cosmic string number density is Blanco-Pillado:2013qja

for the loops in radiation dominated era, and

for the loops produced in radiation dominated era but survive until matter domination, where $\Gamma=50$ Blanco-Pillado:2017oxo , and $t_{\rm eq}=2.25\times 10^{36}~{}{\rm GeV}^{-1}$ is the matter-radiation equality time. In the matter dominated era, the loop number density is

Up to now, the GW spectrum is calculable for a given $G\mu$.

Using Eq. (24), we evaluate the GW signals for the four BPs in Table 2, and plot them as red, green, blue and purple lines in Fig. 2. Due to the continuous GW emission, the signal spectra are flat in a large frequency range. The sensitivity curves of current and future GW detectors are also plotted in the figure, including the pulsar timing arrays (PTAs) NANOGrav McLaughlin:2013ira ; NANOGRAV:2018hou ; Aggarwal:2018mgp ; Brazier:2019mmu , PPTA Manchester:2012za ; Shannon:2015ect , EPTA Kramer:2013kea ; Lentati:2015qwp ; Babak:2015lua , IPTA Hobbs:2009yy ; Manchester:2013ndt ; Verbiest:2016vem ; Hazboun:2018wpv and SKA Carilli:2004nx ; Janssen:2014dka ; Weltman:2018zrl , the space-based laser interferometers LISA LISA:2017pwj , BBO Crowder:2005nr , TianQin TianQin:2015yph ; Hu:2017yoc ; TianQin:2020hid and Taiji Hu:2017mde ; Ruan:2018tsw , and the ground-based interferometers LIGO LIGOScientific:2014qfs ; LIGOScientific:2019vic , CE Reitze:2019iox and ET Punturo:2010zz ; Hild:2010id ; Sathyaprakash:2012jk .⁵⁵5The sensitivity curves shown here are strain noise spectra. For the corresponding power-law-integrated sensitivity curves and peak-integrated sensitivity curves, see Ref. Schmitz:2020syl and the references therein. The existing constraints are shown as shaded regions, while the projected sensitivities are plotted as dashed lines.

The current searches for the SGWB of PTAs constrains $G\mu<10^{-11}$ Ringeval:2017eww ; Blanco-Pillado:2017rnf which requires the $M_{\nu_{R}}/\lambda_{R}<3.2\times 10^{13}$ GeV according to Eq. (24). In Fig. 2, we can see that BP1 (with $M_{\nu_{R}}=1.5\times 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 NANOGrav:2020bcs , which might be a hint for the stochastic GW background from cosmic string Ellis:2020ena ; Bian:2020urb ; Blasi:2020mfx . 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 ground-based 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_{\nu_{R}}=3.0\times 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 $\sim 10^{9}$ GeV and $\sim 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.⁶⁶6For example, the LIGO-Virgo O3 data can exclude $G\mu\gtrsim 4\times 10^{-15}$ for some specific string network models LIGOScientific:2021nrg , much stronger than the results presented in Fig. 2.

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 $\mathcal{O}(10^{9}~{}{\rm GeV})$ and $\mathcal{O}(10^{13}~{}{\rm GeV})$ can be probed by future GW experiments; especially, for DM with mass $\sim 10^{13}$ GeV, the signal might already be reached by the recent NANOGrav results, but more data are needed to clarify this.