1 Introduction↩︎

The visible matter content in the present Universe is highly asymmetric, referred to as the baryon asymmetry of the Universe (BAU) [1], [2]. Starting with a baryon symmetric Universe, such an asymmetry can be generated dynamically if Sakharov’s conditions are satisfied [3]. However, the standard model (SM) of particle physics fails to satisfy them in required amount, making the origin of BAU a longstanding puzzle in particle physics and cosmology. This has led to several beyond standard model (BSM) frameworks like baryogenesis [4], [5] as well as leptogenesis [6]. In leptogenesis, a non-zero lepton asymmetry is generated first which later gets converted into baryon asymmetry by electroweak sphalerons [7]. One interesting aspect of leptogenesis is its connection to the origin of light neutrino masses within canonical seesaw mechanisms like type-I seesaw [8]–[12]. In such a scenario, the out-of-equilibrium CP violating decay of a heavy right handed neutrino (RHN) is responsible for generating the required lepton asymmetry. In order to generate non-zero CP asymmetry while being consistent with neutrino oscillation data, at least two copies of such RHNs are required. Additionally, with hierarchical RHN spectrum, there exists a lower bound on the scale of thermal leptogenesis, known as the Davidson-Ibarra (DI) bound \(m_{N_{1}} \gtrsim 10^9\) GeV [13]. This keeps canonical leptogenesis models out of reach from terrestrial experiments.

Due to the absence of direct signatures, there have been attempts to probe high scale seesaw or leptogenesis scenarios via stochastic gravitational wave (GW) observations [14]–[28] with the GW source being cosmic strings [14]–[17], [24], domain walls [18], [25], [27]–[36] or bubbles generated at first order phase transition (FOPT) [19]–[23], [26]. However, all these probes require extensions of the minimal seesaw model and the origin of GW is not directly connected to the origin of lepton asymmetry. Recently, there have been attempts to find GW signature of such leptogenesis scenarios directly from the decay of heavy RHN or similar fields involved in other seesaw-leptogenesis models [37]–[41]. This involves gravitons emitted via bremsstrahlung-like process during RHN decay into lepton and Higgs. As pointed out in some of these works, the GW amplitude is enhanced with some hope for future detection only for very heavy RHN having mass close to the grand unified theory (GUT) scale \(\sim 10^{16}\) GeV. For thermal leptogenesis, this requires a very high reheat temperature, saturating the upper bound from cosmic microwave background (CMB) observations [42]. On the other hand, thermal leptogenesis with very heavy RHN mass suffers from the strong washout problem due to fast \(\Delta L=2\) processes [43], [44]. Finally, the peak of the GW spectrum from RHN decay remains buried under the GW background spectrum from thermal bath [45]–[47], making it difficult for experimental verification of the contribution from RHN decay alone.

Motivated by these, we consider a non-thermal leptogenesis scenario to avoid the above-mentioned issues related to thermal leptogenesis. While this can occur due to non-thermal production of RHNs during reheating after inflation [48]–[67], we consider the production of RHN from the evaporation of ultra-light primordial black holes (PBH) [68]–[70] with initial mass \(M_\text{in}\lesssim 3.4\times 10^8\,\text{g}\) which evaporate before the epoch of the big bang nucleosynthesis (BBN). Graviton bremsstrahlung from RHNs produced from reheating has already been studied in [40]. While the authors of [38] considered a setup with ultra-light PBH and type-I seesaw, they considered leptogenesis from thermal RHNs while graviton bremsstrahlung occurs dominantly from an ultra-long-lived RHN produced from PBH evaporation. Graviton bremsstrahlung from heavy scalar particle produced from PBH evaporation was studied by [71] but without any connection to the origin of neutrino mass or leptogenesis. In this work, we consider a type-I seesaw leptogenesis scenario with a PBH dominating epoch such that heavy RHNs are produced non-thermally and hence does not suffer from the problem of strong \(\Delta L=2\) washout [72]. Due to the early matter dominated phase of ultra-light PBH, it also separates the RHN decay originated GW spectra from the thermal background resulting in two distinct peaks. The PBH dominated epoch also sources other GW background PBH density fluctuations [73]–[77] as well as Hawking evaporation into gravitons [78]. For the most optimistic scenario with the largest possible GW amplitude from RHN decay, we find a high-frequency GW spectrum with three distinct peaks with the chosen parameter space being consistent with the observed neutrino oscillation data and baryon asymmetry via leptogenesis.

This paper is organized as follows. In section 2, we summarize the basics of PBH-assisted leptogenesis. In section 3, we discuss different GW sources in this scenario and finally conclude in section 4

2 PBH-assisted leptogenesis↩︎

2.1 Primordial Black Holes↩︎

Primordial black holes, originally proposed in [79] and later by Hawking [68], [69] have gaine significant attention due to their cosmological consequences [80], [81]. A recent review of PBH can be found in [70]. In the early universe, PBH can be formed in a variety of ways like, from inflationary perturbations [82]–[86], first-order phase transition (FOPT) [87]–[90], the collapse of topological defects [91], [92] etc. We remain agnostic about their specific production mechanism but assume their formation in a radiation-dominated epoch after the end of inflation. We also assume a monochromatic mass spectrum for PBH with initial mass \(M_{\rm in}\), zero spin and vanishing electric charge. Apart from initial mass, another key parameter related to PBH is the initial energy fraction \[\begin{align} \beta \equiv \frac{\rho_{\rm BH}(T_{\rm in})}{\rho_{\rm R}(T_{\rm in})}, \end{align}\] where \(T_{\rm in}\) denotes the plasma temperature at the time of PBH formation, while \(\rho_{\rm BH}\) and \(\rho_{\rm R}\) represent the PBH and radiation energy densities, respectively.

For PBH formed during radiation domination, the initial PBH mass is typically a fraction of the horizon mass at formation and can be expressed as \[M_{\rm in} = \gamma \frac{4\pi}{3\,\mathcal{H}^3(T_{\rm in})} \,\rho_{\rm R}(T_{\rm in}),\] where \(\mathcal{H}\) is the Hubble expansion rate and \(\gamma\) accounts for the efficiency of gravitational collapse. Numerical studies suggest a representative value \(\gamma \simeq 0.2\) [83]. The corresponding formation time and temperature can be obtained as \[\begin{align} t_{\rm in} &=& \frac{M_{\rm in}}{8\, \pi \gamma M_{P}^2}, \\ T_\text{in} &=& \Biggl(\frac{1440\,\gamma^2}{g_*\left(T_\text{in}\right)}\Biggr)^{1/4}\,\sqrt{\frac{M_P}{M_\text{in}}}\,M_P\,, \end{align}\] with \(M_{P}\) denoting the reduced Planck mass. The upper bound on the inflationary scale [42], [93] provides a lower bound on the initial PBH mass \(M_{\rm in} \gtrsim 0.1 \text{g}\).

Following their formation, PBH lose mass through Hawking evaporation by emitting particles into the surrounding thermal bath [68], [69]. Here, we use \(M_{\rm in}\) to denote the initial mass of the black hole at formation, while \(M_{\rm BH}\) denotes the black hole mass at a later stage of its evolution. At a particular instance of evolution, the Hawking temperature and Bekenstein–Hawking entropy possessed by PBH are given by \[\begin{align} T_{\rm BH} &=& \frac{M_P^2}{M_{\rm BH}}, \\[2mm] S &=& \frac{1}{2} \left( \frac{M_{\rm BH}}{M_P} \right)^2 = \frac{1}{2} \left( \frac{M_P}{T_{\rm BH}} \right)^2 , \label{eq:entropy} \end{align}\tag{1}\] respectively.

The net mass loss rate of a PBH due to Hawking evaporation can be expressed as \[\begin{align} \label{eq:massloss} \frac{dM_{\rm BH}}{dt} = \dot{M}_{\rm BH} \simeq - \epsilon \frac{M^4_{P}}{M^2_{\rm BH}}, \end{align}\tag{2}\] where \(\epsilon\) carries the information of the emitted particles and is given as \[\epsilon = \frac{ 3.8\, \pi g_{*, H}(T_{\rm BH})}{480}, \,\, g_{*, H}(T_{\rm BH}) = \sum_i \omega_i g_{i,H}, \,\,g_{i,H}= \begin{cases} 1.82 &\text{for }s=0\,,\\ 1.0 &\text{for }s=1/2\,,\\ 0.41 &\text{for }s=1\,,\\ 0.05 &\text{for }s=2\,,\\ \end{cases} \label{eq:gsh}\tag{3}\] with \(\omega_i=2s_i+1\) for massive particles of spin \(s_i\), \(\omega_i=2\) for massless species with \(s_i>0\), and \(\omega_i=1\) for spinless species. Integrating Eq. 2 and assuming radiation domination throughout the evolution of a PBH, its corresponding evaporation temperature in terms of initial mass can be calculated as \[\begin{align} \label{eq:Tev} T_{\rm ev} = \left(\frac{40}{\pi^2 g_{*}(T_{\rm ev})}\right)^{1/4} \left(\frac{ 3 \, \epsilon\, M^{5}_{\rm P}}{M^{3}_{\rm in}}\right)^{1/2}. \end{align}\tag{4}\] However, if PBH dominate the energy density of the Universe at some point in their evolution, the bath temperature after their complete evaporation becomes \(\bar{T}_{\rm ev} = \frac{2}{\sqrt{3}} T_{\rm ev}\). Assuming PBH domination, an upper bound on PBH mass \(M_{\rm in} \lesssim 4\times 10^{8} \text{g}\) is obtained by demanding that PBH evaporate before the onset of Big Bang Nucleosynthesis (BBN), \(\bar{T}_{\rm ev} \gtrsim T_{\rm BBN} \sim 4 \,\text{MeV}\). Hence, throughout this work, we restrict our analysis to primordial black holes in the mass range: \(0.1\, \text{g} \lesssim M_{\rm in} \lesssim 4\times 10^{8} \, \text{g}\) with an initial abundance (\(\beta\)) large enough for PBH to dominate the energy density of the Universe prior to their evaporation.

2.2 Non-thermal leptogenesis↩︎

As mentioned earlier, we consider a type-I seesaw mechanism where three gauge singlet RHNs \(N_i, i\in (1-3)\) are added to the SM. The relevant part of the Lagrangian is given by \[-\mathcal{L}_{N} \supset (Y_\nu)_{\alpha i} \overline{\ell_{L\alpha}} \tilde{H} N_i + \frac{1}{2} m_{N_i} \overline{N^c_i} N_i + {\rm h.c.} \label{eq:10}\tag{5}\] where \(\ell_L, H\) denote lepton and Higgs doublet of the SM respectively. Once the electroweak symmetry gets spontaneously broken after the Higgs acquires a non-zero vacuum expectation value (VEV) \(v\), light neutrino masses arise in the seesaw limit \(M_D = Y_\nu v/\sqrt{2} \ll m_N\) as \[m_\nu = -M_D m^{-1}_N M^T_D = -Y_\nu m^{-1}_N Y^T_\nu \frac{v^2}{2}.\]

In thermal leptogenesis, RHNs are produced from the thermal bath due to the Yukawa interactions. In non-thermal leptogenesis setup we are considering here, RHNs are produced from PBH evaporation. There have been several works in the literature where such role of PBH evaporation on generating baryon asymmetry of the Universe has been dstudied [60], [68], [72], [81], [94]–[104]. Due to gravitational interactions, the evaporation of PBH produces all particles including the RHNs. Here, we consider a scenario where the RHNs produced from PBH evaporation never reaches thermal equilibrium such that we remain in the non-thermal leptogenesis regime. Subsequently, the CP violating decay of these RHNs generate a non-zero lepton asymmetry which gets converted into baryon asymmetry via \(B+L\)-violating electroweak sphaleron transitions [7].

In order to compute the PBH-assisted baryon asymmetry, we need to track the energy densities of PBH (\(\rho_{\rm BH}\)), radiation (\(\rho_{\rm R}\)), RHNs (\(\rho_{\rm N}\)), as well as the bath temperature (\(T\)) throughout the PBH evolution. The corresponding Boltzmann Equations (BEs) are \[\begin{align} \label{eq:BE95PBH} \frac{d \rho_{\rm BH}}{dt} &+& 3\,\mathcal{H} \rho_{\rm BH} = \frac{\dot{M}_{\rm BH}}{M_{BH}} \rho_{\rm BH} \nonumber \\ \frac{d \rho_{\rm R}}{dt} &+& 4\,\mathcal{H} \rho_{\rm R} = - \frac{\epsilon_{\rm SM}}{\epsilon_{\rm SM}+\epsilon_{\rm N}} \frac{\dot{M}_{\rm BH}}{M_{BH}} \rho_{\rm BH} \nonumber \\ \frac{d \rho_{\rm N}}{dt} &+& 3\,\mathcal{H} \rho_{\rm N} = - \frac{\epsilon_{\rm N}}{\epsilon_{\rm SM}+\epsilon_{\rm N}} \frac{\dot{M}_{\rm BH}}{M_{BH}} \rho_{\rm BH} \nonumber \\ \mathcal{H}^{2} &=& \frac{\rho_{\rm BH} + \rho_{\rm R} + \rho_{\rm N}}{3 M^2_{\rm P}} \nonumber \\ \frac{\dot{T}}{T} &=& - \frac{1}{\Delta} \left(\mathcal{H} + \frac{\epsilon_{\rm SM}}{4 (\epsilon_{\rm SM}+\epsilon_{\rm N})} \frac{\rho_{\rm BH}}{\rho_{\rm R}} \frac{\dot{M}}{M} \frac{g_{*}(T)}{g_{*,s}(T)}\right). \end{align}\tag{6}\] Here \(^{.}\) denotes derivatives with respect to cosmic time and \(\Delta = 1 + \frac{T}{3} \frac{d\, g_{*,s}(T)/dT}{g_{*,s}(T)}\). The quantities \(\epsilon_{\rm SM}\) and \(\epsilon_{\rm N}\) carry the information of PBH emitted particles for SM and RHN as described in Eq. 3 . By solving the BEs given in Eqs. 6 , the number density of \(N_{1}\) is obtained at the time of PBH evaporation. Finally, the total baryon asymmetry today from PBH-assisted leptogenesis can be expressed as \[\begin{align} \label{eq:asymmetry1} Y^{\Delta L} = \epsilon^{\Delta L}_{1} a_{\rm sph} \frac{n_{\rm N_{1}}(T_{\rm ev})}{s (T_{\rm ev})}, \end{align}\tag{7}\] where \(a_{\rm sph}= 28/79\) is the sphaleron conversion factor. Here, the CP asymmetry parameter is : \[\begin{align} \epsilon^{\Delta L}_{1} = \frac{1}{8\pi} \sum_{j=2,3} \frac{\text{Im}[(Y^{\dagger}_{\nu}Y_{\nu})^2_{1j}]}{(Y^{\dagger}_{\nu}Y_{\nu})_{11}}\left\{{f_{V}\left(\frac{m^2_{N_{j}}}{m^2_{N_{1}}}\right)+f_{S}\left(\frac{m^2_{N_{j}}}{m^2_{N_{1}}}\right)}\right\}, \end{align}\] where \(Y_{\nu}\) is the Yukawa matrix corresponding to the term \((Y_{\nu})_{\alpha i} \bar{\ell}_{L\alpha}\tilde{H} N_{i}\) and \(f_{V}(x) = \sqrt{x}\left[1-(1+x)\text{ln}\left(\frac{1+x}{x}\right)\right]\), and \(f_{S}(x) = \frac{\sqrt{x}}{1-x}\) denote the contributions of vertex and self-energy respectively. In the hierarchical RHN mass limit of \(m_{N_{1}} \ll m_{N_{2,3}}\), \[\begin{align} f_{V}(x) + f_{S}(x) \simeq - \frac{3}{2\sqrt{x}}, \,\,\,\,\,\,\, \text{for}\,\,\,\, x\gg 1. \end{align}\]

Moreover, the upper bound on the CP asymmetry parameter can be connected with the low-energy neutrino parameter as \[\begin{align} |\epsilon^{\Delta L}_{1}| \lesssim \frac{3 m_{N_{1}} \sqrt{(\Delta m_{\rm atm})^2}}{8 \pi v^{2}}. \end{align}\]

2.2.1 Analytical estimation of Baryon Asymmetry↩︎

In terms of the number of the lightest RHN \(N_{1}\) species, \(N_{N_{1}}\), produced from the complete evaporation of a PBH, the Eq. 7 can be expressed as \[\begin{align} Y^{\Delta L} = \epsilon^{\Delta L}_{1} a_{\rm sph} \frac{N_{\rm N_{1}} n_{\rm BH}(a_{\rm ev})}{s (a_{\rm ev})}. \end{align}\] The expression for \(N_{N_{1}}\) is obtained from the generic expression of the number of a species with mass \(m_{X}\) and internal degrees of freedom \(g_{X}\) produced from the complete evaporation of a PBH, which is given as \[\begin{align} \label{eq:Nx} N_{X} \simeq \frac{27}{128} \frac{\xi g_{j} \zeta(3)}{\pi^3 \epsilon} \begin{cases} \left(\frac{M_{\rm in}}{M_P}\right)^{2},& \text{if m_{X} < T^{\rm in}_{\rm BH}} \\ \left(\frac{M_P}{m_{X}}\right)^{2},& \text{if m_{X} > T^{\rm in}_{\rm BH}}. \end{cases} \end{align}\tag{8}\] Here \(\xi\) takes the value of \(1\) for bosons and \(3/4\) for fermions. Using the relation \(\rho_{\rm BH}(a_{\rm ev}) = \rho_{\rm R}(a_{\rm ev}) \simeq n_{\rm BH}(a_{\rm eq})\times M_{\rm in}\), the number density of PBH at the end of evaporation is found to be \[\begin{align} \label{eq:naev} n_{\rm BH} (a_{\rm ev}) = 12\, M^3_{\rm P}\, \epsilon^2 \left(\frac{M_{P}}{M_{\rm in}}\right)^7. \end{align}\tag{9}\] Combining the expression of Eqs. 8 , 9 and 4 in Eq. 7 , the PBH assisted baryon asymmetry can be approximately written as \[\begin{align} \label{eq:lepto95analy} Y^{\Delta L} \simeq 1.04\times10^{-2}\, \frac{g^{3/4}_{*}(T_{\rm ev})}{g_{*,s}(T_{\rm ev})}\, \epsilon^{\Delta L}_{1}\, a_{\rm sph}\, \begin{dcases} \left(\frac{M_{\rm P}}{M_{\rm in}}\right)^{1/2} ,& \text{if m_{N_{1}} < \,T_{\rm in}} \\ \left(\frac{M^{9}_{\rm P}}{M^{5}_{\rm in} m^{4}_{N_{1}}}\right)^{1/2} ,& \text{if m_{N_{1}} > \,T^{\rm in}_{\rm BH}}. \end{dcases} \end{align}\tag{10}\]

The left panel of Fig. 1 shows the viable parameter space for a hierarchical RHN mass spectrum in the \(M_{\rm in}\) – \(m_{N_{1}}\) plane. The observed behavior can be understood from Eq. 10 . In the regime \(m_{N_{1}} < T^{\rm in}_{\rm BH}\), the generated baryon asymmetry scales as \(Y^{\Delta L} \propto \frac{m_{N_{1}}}{\sqrt{M_{\rm in}}},\) implying a positive correlation between \(M_{\rm in}\) and \(m_{N_{1}}\) along the lower boundary of the allowed region. In contrast, for \(m_{N_{1}} > T^{\rm in}_{\rm BH}\), the asymmetry scales as \(Y^{\Delta L} \propto \frac{1}{M^{5/2}_{\rm in} m_{N_{1}}},\) leading to an anti-correlation between \(M_{\rm in}\) and \(m_{N_{1}}\) along the upper boundary. These scaling behaviors explain the distinct slopes of the lower and upper limits of the viable parameter space. We find that the allowed PBH mass range in the hierarchical scenario is approximately \(0.1~{\rm g} \lesssim M_{\rm in} \lesssim 10~{\rm g}\). Additionally, the bottom gray-shaded region denotes the parameter space incompatible with non-thermal leptogenesis, obtained by approximately imposing the condition \(m_{N_{1}}< T_{\rm ev}\).

The generated lepton asymmetry can be substantially enhanced when two RHNs are nearly degenerate in mass, i.e., \(\Delta m \equiv m_{N_2}-m_{N_1} \ll \bar{m} \equiv (m_{N_1}+m_{N_2})/2\), giving rise to the resonant leptogenesis scenario [105]. In particular, when the mass splitting becomes comparable to the RHN decay width, \(\Delta m \simeq \Gamma\), the CP asymmetry parameter can be resonantly amplified and reach values of \(\mathcal{O}(1)\). As a consequence, successful leptogenesis can be achieved over a significantly larger region of parameter space, as illustrated in the right panel of Fig. 1 where we consider the CP asymmetry to be \(0.1\). The vertical dashed lines indicate the condition \(T_{\rm ev}=T_{\rm sph}\), corresponding to the epoch of sphaleron decoupling.

3 Gravitational wave signatures↩︎

3.1 Gravitational Bremsstrahlung↩︎

As mentioned earlier, graviton bremsstrahlung from the decay of RHN [37]–[40] offers an interesting GW probe directly connected to the processes responsible for generating lepton asymmetry. The interaction between the SM and BSM sectors is described by the action \[\begin{align} S \supset \int d^{4}x \sqrt{-g} \left[ \frac{M_{\rm P}^{2}}{2} R + \mathcal{L}_{\rm SM} + \mathcal{L}_{N} \right], \end{align}\] where \(R\) denotes the Ricci scalar and \(g \equiv \det(g_{\mu\nu})\) is the determinant of the spacetime metric. In the weak-field limit, the metric can be expanded around the flat Minkowski background as \[\begin{align} g_{\mu\nu} \simeq \eta_{\mu\nu} + \frac{2}{M_{\rm P}} h_{\mu\nu}, \end{align}\] where \(h_{\mu\nu}\) represents the canonically normalized graviton field. The corresponding interaction between the graviton and the energy-momentum tensor \(T_i^{\mu\nu}\) of a generic particle species \(i\) is given by \[\begin{align} \mathcal{L}_{\rm int}^{\rm gravity} \supset -\frac{2}{M_{\rm P}} h_{\mu\nu} T_i^{\mu\nu}. \end{align}\] The leading-order processes contributing to the GW spectrum through gravitational bremsstrahlung accompanying the decay of RHNs are depicted in Fig. 2. As the graviton is emitted from one of the external legs of the same RHN decay process relevant for leptogenesis, this offers a powerful and unavoidable probe of leptogenesis.

In order to obtain the spectrum of GW produced via bremsstrahlung, we first need to solve the BE for the energy density of graviton, \(\rho_{\rm gw}\), which is given as [37], [39]

\[\begin{align} \frac{d \rho_{\rm gw}}{dt} + 4 \mathcal{H} \rho_{\rm gw} = \int\frac{d^3 k}{(2\pi)^{3}} E_{k} \frac{1}{2 E_{k}} \int\frac{d^3 p}{(2\pi)^{3} \,2 E_{p}} \int\frac{d^3 q}{(2\pi)^{3}\, 2 E_{q}} \int\frac{d^3 r}{(2\pi)^{3} \, 2E_{r}} |\mathcal{M}|^{2} f_{N}(p) \nonumber \\ \times (2\pi)^4 \delta^{(4)}(P-Q-K-R), \end{align}\] where \(P,\, Q,\, R,\, K\) denote the four momentum of RHN, lepton, Higgs, and graviton, respectively. \(f_{N}\) is the distribution function of RHN and \(|\mathcal{M}|^{2}\) represents the spin and polarization averaged matrix element squared. The Eq. 11 can be written as \[\begin{align} \frac{d \rho_{\rm gw}}{dt} + 4 \mathcal{H} \rho_{\rm gw} &=& \int\frac{d^3 p}{(2\pi)^{3}} f_{N}(p)\, \frac{m_{N_{1}}}{E_{N_{1}}}\int \frac{d\Gamma^{1\to3}}{dE_{k}} E_{k}\, dE_{k} \end{align}\] which in the non-relativistic limit can be written as \[\begin{align} \label{eq:BE95gw} \frac{d \rho_{\rm gw}}{dt} + 4 \mathcal{H} \rho_{\rm gw} &=& \left[\int \frac{d \Gamma ^{1 \to 3}}{d E_{\rm k}}\left(\frac{E_{\rm k}}{m_{\rm N_{1}}}\right)dE_{\rm k}\right] \rho_{N_{1}}. \end{align}\tag{11}\] Here, \[\begin{align} \label{eq:395decay95width} \frac{d\Gamma^{1\to 3}}{dE_{\rm k}}= \frac{(Y_{\nu}^{\dagger}Y_{\nu})_{11}}{256 \pi^{3}} \frac{m^2_{N_{1}}}{M^2_{P}} \frac{(x-1)^2(1-2x)}{x}, \end{align}\tag{12}\] with \(x = E_{k}/m_{N_{1}}\). Calculating \(\rho_{\rm gw}\) requires solving Eqs. 11 and 12 along with the coupled BEs in Eqs. 6 . We solve these equations numerically. Below, we also provide an analytical estimation of the calculation of \(\Omega_{\rm GW}h^2\).

3.1.1 Analytical estimation of GW via bremsstrahlung↩︎

The Eq. 11 can be re-written as \[\begin{align} \frac{d}{dt} \left(\frac{d\rho_{\rm gw}}{d E_{k}}\right) + 4 \mathcal{H} \frac{d\rho_{\rm gw}}{d E_{k}} = \frac{d \Gamma ^{1 \to 3}}{d E_{\rm k}} E_{\rm k}\, n_{N_{1}}. \end{align}\] In terms of scale factor \(a\), \[\begin{align} \frac{d}{da} \left(a^5 \frac{d\rho_{\rm gw}}{d \, \text{ln} \, E_{k}}\right) = \frac{a^4}{\mathcal{H}} \frac{d \Gamma ^{1 \to 3}}{d E_{\rm k}} E^2_{\rm k}\, n_{N_{1}}. \end{align}\] Integrating the above equation from the evaporation of PBH, \(a_{\rm ev}\) to the decay of RHN, \(a_{N}\), we get \[\begin{align} \frac{d\rho_{\rm gw}(a_{\rm N})}{d \, \text{ln}\, E_{k}} \simeq \frac{2}{3} \tau_{N}\, n_{N_{1}}(a_{\rm ev}) \left(\frac{a_{\rm ev}}{a_{N}}\right)^3 \left[1-\left(\frac{a_{\rm ev}}{a_{\rm N}}\right)^3\right]\left.\left(\frac{d \Gamma ^{1 \to 3}}{d E_{\rm k}}\right)\right|_{a_{N}} E^2_{\rm k} (a_{\rm N}), \end{align}\] where \(\tau_{N}\) is the decay lifetime of RHN. Hence, the GW amplitude today can be written as \[\begin{align} \Omega_{\rm GW}h^2 &=& \frac{h^2}{\rho_{\rm c, 0}} \frac{d\rho_{\rm gw}(a_{\rm N})}{d \, \text{ln}\, E_{k}} \left(\frac{a_{\rm N}}{a_{0}}\right)^{4} \\ \nonumber &\simeq& \frac{h^2}{\rho_{\rm c, 0}} \frac{2}{3} \tau_{N}\, n_{N_{1}}(a_{\rm ev}) \left.\left(\frac{d \Gamma ^{1 \to 3}}{d E_{\rm k}}\right)\right|_{a_{N}} E^2_{\rm k} (a_{\rm N}) \left(\frac{a_{\rm ev}}{a_{N}}\right)^3 \left[1-\left(\frac{a_{\rm ev}}{a_{\rm N}}\right)^3\right] \left(\frac{a_{\rm N}}{a_{0}}\right)^{4}, \end{align}\] with, \(\rho_{c,0}\) denoting the critical energy density today. The frequency today, \(f\), is related to the energy of the graviton at the time of RHN decay as \[\begin{align} E_{k}(a_{\rm N}) = 2 \pi f \left(\frac{a_{0}}{a_{\rm N}}\right). \end{align}\] Simplifying the expression for \(\Omega_{\rm GW}h^2\) in terms of \(f\) for \(m_{N_{1}}<T^{\rm in}_{\rm BH}\), we obtain \[\begin{align} \label{eq:GWs95amp} \Omega_{\rm GW}h^2 \simeq 2.7\times10^{-20} \left(\frac{f}{10^{10} \, \text{Hz}}\right) \left(\frac{m_{N_{1}}}{10^{15}\, \text{GeV}}\right)^2 \sqrt{\frac{1\, \text{g}}{M_{\rm in}}}. \end{align}\tag{13}\] For the scenario where \(m_{N_{1}}>T^{\rm in}_{\rm BH}\), the amplitude takes the form of \[\begin{align} \label{eq:GWs95amp951} \Omega_{\rm GW}h^2 \simeq 3\times 10^{-24} \left(\frac{f}{10^{10}\, \text{Hz}}\right) \left(\frac{1 \, \text{g}}{M_{\rm in}}\right)^{5/2}. \end{align}\tag{14}\] Moreover, as the graviton can carry the maximum energy of \(m_{N_{1}}/2\) at the time of formation, the peak GW frequency can be written as \[\begin{align} f^{\rm peak} &=& \frac{m_{N_{1}}}{4} \left(\frac{a_{N_{1}}}{a_{0}}\right) \\ \nonumber &=& \frac{m_{N_{1}}}{4} \left(\frac{T^3_0\, g_{*,s}(T_{0})}{T^3_N\, g_{*,s}(T_{N})}\right)^{1/3}, \end{align}\] where \(T_{N}\) denotes the temperature when \(N_{1}\) decays. Using the relation \(\mathcal{H}=\frac{1}{2 \tau_{N}}\), \(T_{N}\) can be determined as \[\begin{align} T_{N} = \left(\sqrt{\frac{90}{g_{*,s}(T_{N})}}\,\frac{(Y^{\dagger}_{\nu} Y_{\nu})_{11} m_{N_{1}} M_{\rm P}}{16\pi^2} \right)^{1/2}. \end{align}\] Consequently, the peak frequency can be expressed as \[\begin{align} \label{eq:peak95f} f^{\rm peak} = 1.2\times 10^{16}\, \text{Hz} \left(\frac{10^{-5}}{(Y^{\dagger}_{\nu} Y_{\nu})^{1/2}_{11}}\right) \left(\frac{m_{N_{1}}}{M_{P}}\right)^{\frac{1}{2}}. \end{align}\tag{15}\]

Figure 3 displays the GW spectra generated through bremsstrahlung processes for the three benchmark points (BPs) listed in Table 1. The solid curves correspond to the full numerical results, while the dashed curves represent the analytical approximations. As anticipated from the analytical estimates, the peak frequency shifts toward higher values with increasing \(m_{N_1}\) and decreasing Yukawa coupling. A comparison between BP1 and BP2 reveals a small difference. While the analytical approximation predicts nearly identical GW amplitudes at a fixed frequency for these two benchmark points, the numerical results exhibit a slight suppression of the amplitude for BP2. This discrepancy originates from the different evaluations of \(n_{N_1}(a_{\rm ev})\) in the analytical and numerical treatments. In particular, the analytical approximation becomes less reliable in the regime where \(T_{\rm BH}^{\rm in}\sim m_{N_1}\), resulting in a modest deviation from the full numerical calculation. In contrast, the GW amplitude is noticeably smaller for BP3 than for BP2 at a given frequency. This suppression is primarily due to the larger initial PBH mass, \(M_{\rm in}\), associated with BP3. In Fig. 3, we restrict our analysis to PBH masses up to \(M_{\rm in}=1~{\rm g}\). For heavier PBH, the GW spectrum at a fixed frequency would be further suppressed, as suggested by Eqs. 13 and 14 . Consequently, such spectra are expected to lie below the curves shown in the figure and are therefore not displayed. We stick to very low mass limits of ultra-light PBH such that the graviton bremsstrahlung spectrum remains in the most optimistic range for any possible future detections.

The values of the CP asymmetry parameter required to achieve successful leptogenesis for the benchmark points listed in Table 1 are \(\epsilon^{\Delta L}_{1}=10^{-3}\), \(\epsilon^{\Delta L}_{1}=10^{-5}\), and \(\epsilon^{\Delta L}_{1}=3\times10^{-3}\), respectively. In Appendix 5, we present the Yukawa coupling structures corresponding to specific choices of the lightest neutrino mass and the three complex rotation angles of the orthogonal matrix in Casas–Ibarra parametrisation [113]. These choices are consistent with the observed neutrino mass spectrum and mixing parameters, while simultaneously reproducing the required values of \(\epsilon^{\Delta L}_{1}\) and \((Y_\nu^\dagger Y_\nu)_{11}\) for each benchmark point.

For the resonant leptogenesis scenario, the GW spectrum exhibits the same qualitative features as in the non-resonant case. However, its amplitude is enhanced by approximately a factor of two due to the presence of two nearly degenerate right-handed neutrinos, both of which contribute to the GW production.

In Fig. 4, we present the peak GW amplitude, \(\Omega_{\rm GW}^{\rm peak}h^2\), and the corresponding peak frequency, \(f^{\rm peak}\), in the \(M_{\rm in}\)–\(m_{N_1}\) plane. Throughout this analysis, we fix the lifetime of \(N_1\) such that its decay lifetime is ten times the PBH lifetime. We also assume a hierarchical RHN mass spectrum with \(m_{N_2}=10\,m_{N_1}\) and \(m_{N_3}=50\,m_{N_1}\). The parameter space enclosed by the magenta contours yields successful leptogenesis with hierarchical RHNs, while the gray-shaded region outside the magenta contours is incompatible with the DI bound. Imposing the perturbativity condition \(\max \!\left[(Y_{\nu}^{\dagger}Y_{\nu})_{ii}\right] < 4\pi\), we obtain the corresponding upper bound on \(m_{N_{1}}\), shown by the horizontal cyan-shaded region. Analogous results for resonant leptogenesis are presented in Fig. 5 where we consider \(m_{N_{2}}\approx m_{N_{1}}\) and \(m_{N_{3}} = 50\, m_{N_{1}}\).

From Figs. 4 and 5, it is evident that, over most of the parameter space, the predicted GW peak frequency lies beyond the sensitivity reach of both current and planned GW observatories. Furthermore, even in the regions where the peak frequency falls within the projected sensitivity range of future detectors, the corresponding GW amplitude is too small to produce an observable signal.

3.2 GW from thermal bath↩︎

Apart from the production of GW from the decay of heavy RHNs, a stochastic GW background can also be generated through bremsstrahlung processes in the thermal plasma. This background is commonly referred to as the Cosmic Gravitational Microwave Background (CGMB) [39], [45]–[47]. The production of CGMB starts at very high temperatures, where graviton emission from the thermal bath is efficient. The evolution of the graviton energy density is governed by the Boltzmann equation \[\begin{align} \frac{d \rho_{\rm gw}}{dt} + 4 \mathcal{H} \rho_{\rm gw} &=& \int \frac{d^3k}{(2\pi)^3} \frac{4}{M^2_{\rm P}} T^4 \,\hat{\eta}\left(T,\frac{k}{T}\right), \end{align}\] where \(\hat{\eta}\) denotes the dimensionless source function defined in Ref. [47].

Expressing the above equation in terms of the entropy density \(s\) and neglecting the variation of \(g_{*,s}\) at temperatures well above the electroweak scale, one obtains \[\begin{align} \frac{d}{dt} \left(\frac{\rho_{\rm gw}}{s^{4/3}}\right) = \frac{4}{2\pi^2 M^2_{\rm P}} \frac{1}{s^{4/3}} \int^{\infty}_{0} k^2 dk \, T^4\, \hat{\eta}\left(T,\frac{k}{T}\right). \end{align}\] The dominant production of thermal gravitons takes place between the end of reheating and the onset of PBH domination. Integrating over this period yields \[\begin{align} \frac{\rho_{\rm gw}(T_{\rm E})}{s^{4/3} (T_{\rm E})} = \frac{12 \sqrt{10}}{2\pi^3 M_{\rm P}} \frac{1}{\left(\frac{2\pi^2}{45}\right)^{4/3}} \int^{T_{\rm RH}}_{T_{\rm E}} \frac{dT}{T^3 g^{4/3}_{*,s}(T) g^{1/2}_{*,s}} \int^{\infty}_{0} k^2 dk \, \hat{\eta}\left(T,\frac{k}{T}\right), \end{align}\] where \(T_{\rm RH}\) and \(T_{\rm E}\) denote the reheating temperature and the temperature at the onset of PBH domination, respectively. In terms of per logarithimic wavenumber interval and using the fact that the \(\hat{\eta}\) is weakly dependent on temperature, the above equation becomes \[\begin{align} \frac{d \rho_{\rm gw}}{d\, \text{ln}\, k_{\rm E}} (T_{\rm E}, k_{\rm E}) &=& \frac{12 \sqrt{10}}{2\pi^3 M_{\rm P}} g^{1/3}_{*,s} (T_{\rm E}) T^4_{\rm E} \frac{k^3_{\rm E}}{T^3_{\rm E}}\int^{T_{\rm RH}}_{T_{\rm E}} \frac{dT}{g^{4/3}_{*,s}(T) g^{1/2}_{*,s}(T)} \hat{\eta}\left(T,\frac{k_{\rm E}}{T_{\rm E}} \left(\frac{g_{*,s}(T)}{g_{*,s}(T_{\rm E})}\right)^{1/3} \right)\nonumber \\ &\simeq& \frac{12 \sqrt{10}}{2\pi^3 M_{\rm P}} T^4_{\rm E} \frac{k^3_{\rm E}}{T^3_{\rm E}}\, \frac{T_{\rm RH}}{ g^{1/2}_{*,s}(T_{\rm RH})} \hat{\eta}\left(T_{\rm RH},\frac{k_{\rm E}}{T_{\rm E}} \left(\frac{g_{*,s}(T)}{g_{*,s}(T_{\rm E})}\right)^{1/3} \right). \end{align}\] Redshifting this contribution to the present epoch gives \[\begin{align} \Omega^{\rm CGMB}_{\rm GW}h^2 (k_{\rm E}) \simeq \frac{h^2}{\rho_{\rm c, 0}} \left(\frac{T_{0}}{T_{\rm E}}\right)^4 \left(\frac{g_{*,s}(T_{0})}{g_{*,s}(T_{\rm E})}\right)^{4/3} \frac{1}{\alpha^{4/3}} \frac{d\rho_{\rm gw}}{d\, \text{ln}\, k_{\rm E}} (T_{\rm E}, k_{\rm E}). \end{align}\] Here we use the quantity \(\alpha\) to take into account the dilution of CGMB due to PBH domination era. In terms of frequency today, \(f= \frac{k_{\rm E}}{2\pi} \left(\frac{T_{0}}{T_{\rm E}}\right) \left(\frac{g_{*,s}(T_{0})}{g_{*,s}(T_{\rm E})}\right)^{1/3} \frac{1}{\alpha^{1/3}}\), the GW spectrum today can be expressed as \[\begin{align} \Omega^{\rm CGMB}_{\rm GW}h^2 (f) &\simeq& \frac{h^2\, T_{0}}{\rho_{\rm c, 0}} \left(\frac{g_{*,s}(T_{0})}{g_{*,s}(T_{\rm E})}\right)^{\frac{1}{3}} \frac{1}{\alpha^{1/3}} \frac{48 \sqrt{10}}{g^{1/2}_{*}(T_{\rm RH})} \frac{T_{\rm RH}}{M_{\rm P}} f^3\, \hat{\eta} \left(T_{\rm RH},\, \frac{2\pi f}{T_{0} } \left(\frac{g_{*,s}(T_{\rm E})}{g_{*,s}(T_{0})}\right)^{\frac{1}{3}} \alpha^{\frac{1}{3}}\right) \nonumber \\ &\simeq& 4.69\times10^{-12} \,\frac{1}{\alpha^{1/3}} \, \left(\frac{f}{10^{9} \, \text{Hz}}\right)^3 \left(\frac{T_{\rm RH}}{M_{\rm P}}\right) \hat{\eta} \left(T_{\rm RH},\, 0.056\times\left(\frac{f}{10^{9} \, \text{Hz}}\right) \alpha^{\frac{1}{3}}\right).\nonumber \\ \end{align}\]

The entropy dilution factor arising from PBH evaporation is given by [72] \[\begin{align} \alpha \equiv \frac{s(\tilde{T})}{s(T_{\rm ev})} = \left(1+ \beta \frac{T_{\rm in}}{T_{\rm ev}}\right)^{3/4}, \end{align}\] where \(\tilde{T}\) is the temperature of the SM plasma after PBH evaporation occurs.

Fig. 6 shows the CGMB for the benchmark points shown in Table 1, considering the maximum allowed reheating temperature \(10^{16}\) GeV. Comparing these results with those in Fig. 3, we find that the CGMB spectrum exhibits a lower peak frequency and a larger peak amplitude.

3.3 Gravitational Waves from PBH↩︎

PBH can generate a stochastic GW background through a variety of physical mechanisms. These include: (i) the emission of gravitons via Hawking evaporation, producing GWs at extremely high frequencies [78], (ii) GWs sourced by density fluctuations associated with the spatial clustering and Poisson statistics of PBH after their formation [74], [76], [77], [114], and (iii) second-order GWs induced by the large primordial scalar perturbations responsible for PBH formation in the early Universe [115]. In this work, we focus on the GW backgrounds arising from Hawking evaporation and from density fluctuations associated with the PBH population.

3.3.1 From density fluctuations↩︎

The stochastic GW background generated by Poisson fluctuations in the PBH number density can be approximated by [74] \[\Omega_{\rm GW}^{\rm PBH}(f) \simeq \Omega_{\rm GW}^{\rm peak} \left(\frac{f}{f_{\rm peak}^{\rm PBH}}\right)^{11/3} \Theta\!\left(f_{\rm peak}^{\rm PBH}-f\right), \label{eqn:omgw}\tag{16}\] where \(\Omega_{\rm GW}^{\rm peak}\) denotes the maximum GW energy density fraction, given by \[\Omega_{\rm GW}^{\rm PBH,\, peak} \simeq 2\times10^{-6} \left(\frac{\beta}{10^{-8}}\right)^{16/3} \left(\frac{M_{\rm in}}{10^7\,{\rm g}}\right)^{34/9}. \label{eqn:omgpeak}\tag{17}\]

The above description relies on treating the PBH distribution as a continuous fluid. However, this approximation breaks down on scales smaller than the average separation between PBH. Consequently, the induced GW spectrum exhibits an ultraviolet cutoff, whose characteristic frequency corresponds to the comoving scale associated with the mean PBH separation. The resulting peak frequency can be estimated as \[f^{\rm PBH,\, peak} \simeq 1.7\times10^{3}\,{\rm Hz} \left(\frac{M_{\rm in}}{10^4\,{\rm g}}\right)^{-5/6}. \label{eqn:fpk}\tag{18}\]

Fig. 7 shows the GW from density fluctuations due to PBH. The spectra for BP1 and BP2 are exactly the same due to identical PBH parameters and are shown by the red curve.

3.3.2 From PBH evaporation↩︎

The production of gravitons from PBH evaporation has been investigated in several studies [71], [100], [116]–[118]. In this work, we follow Ref. [116] to estimate the GW spectrum generated from PBH evaporation.

The present-day GW energy density spectrum can be approximated as \[\begin{align} \Omega^{\rm PBH,ev}_{\rm GW} h^2(f) \simeq 1.36\times 10^{-21} \left(\frac{g_{*,H}}{100}\right)^2 \left(\frac{1~{\rm g}}{M_{\rm in}}\right)^2 \left(\frac{f}{10^9~{\rm Hz}}\right)^4 I\!\left(\frac{2\pi f}{T_{\rm BH}^0}\right), \end{align}\] where \(T_{\rm BH}^0\) denotes the PBH Hawking temperature redshifted to the present epoch and is given by \[\begin{align} T_{\rm BH}^0 \simeq 4.53\times10^{15}~{\rm Hz} \left(\frac{100}{g_{*,s}(T_{\rm ev})}\right)^{1/12} \left(\frac{100}{g_{*,H}}\right)^{1/2} \left(\frac{M_{\rm in}}{10^5~{\rm g}}\right)^{1/2}. \end{align}\] The function \(I\) is defined as \[\begin{align} I\!\left(\frac{2\pi f}{T_{\rm BH}^0}\right) = \int_0^{z_{\rm max}} \frac{\sqrt{1+z}\, dz}{\exp\!\left[(1+z)\,2\pi f/T_{\rm BH}^0\right]-1}, \end{align}\] where \[\begin{align} 1+z_{\rm max} \simeq 7.36\times10^7 \left(\frac{M_{\rm in}}{1~{\rm g}}\right)^{4/3} \beta^{1/3}. \end{align}\]

Fig. 8 shows the GW spectrum arising from PBH evaporation. Since the PBH parameters are identical for BP1 and BP2, the corresponding GW spectra completely overlap and are represented by the red curve.

Finally, in Fig. 9 we present the total GW spectrum obtained by combining all relevant sources. At low frequencies, the signal is dominated by GW induced by enhanced density perturbations. In the intermediate frequency range, \(10^{10}\,{\rm Hz} \lesssim f \lesssim 10^{14}\,{\rm Hz}\), two main contributions arise: the CGMB and graviton emission from PBH evaporation. Among these, the latter provides the dominant contribution and is clearly visible in the figure. At even higher frequencies, around \(f \sim 10^{15}\,{\rm Hz}\), the GW signal produced through RHN bremsstrahlung becomes dominant, carrying characteristic signatures of high-scale leptogenesis.

We note from Fig. 9 that, for the benchmark points listed in Table 1, the GW spectrum arising from RHN decay in hierarchical leptogenesis contributes significantly only in the ultra-high-frequency range, \(\sim 10^{14}-10^{15}~\mathrm{Hz}\). In this regime, distinguishing the GW spectrum from RHN decay may be challenging for some of the benchmark points due to the presence of other dominant GW sources, particularly the contribution from graviton emission during PBH evaporation. This is even more difficult for most of the benchmark points in case of thermal high-scale leptogenesis without a PBH-dominated epoch, as shown in Appendix 6. From Eq. 15 , it is evident that the peak frequency, and consequently the peak amplitude, can be enhanced by considering a smaller Yukawa coupling for \(N_1\). In the framework of hierarchical leptogenesis, however, reducing the Yukawa coupling simultaneously suppresses the CP asymmetry parameter, making it difficult to generate the observed baryon asymmetry of the Universe. This limitation can be circumvented in resonant leptogenesis, where the Yukawa couplings and the CP asymmetry parameter become effectively decoupled. As a result, one can consider significantly smaller Yukawa couplings while still obtaining the required baryon asymmetry through resonant enhancement.

To illustrate this possibility, we present three benchmark points for resonant leptogenesis in Table 2, featuring smaller Yukawa couplings than those in Table 1. Here, \(\Delta m\) denotes the mass splitting between \(m_{N_2}\) and \(m_{N_1}\), while the mass of \(N_1\) and the remaining PBH parameters are kept identical to those in Table 1. The resulting GW spectra are shown in Fig. 10. In contrast to the hierarchical leptogenesis case displayed in Fig. 9, the RHN bremsstrahlung contribution extends over a much broader frequency range and constitutes a substantial component of the total GW spectrum, thereby significantly improving its prospects for observational distinguishability.

Before concluding, we would like add a remark that the energy density of RHNs redshifts more slowly than that of the radiation bath after their production from PBH evaporation. Consequently, for sufficiently small Yukawa couplings, the RHNs can eventually dominate the energy density of the Universe before decaying. In this work, however, we restrict our analysis to benchmark points for which the RHNs remain subdominant throughout their evolution. If an RHN-dominated era were to occur, both the dynamics of PBH-assisted leptogenesis and the resulting GW spectrum would be modified. A dedicated investigation of this scenario is left for future study.

4 Conclusion↩︎

We have revisited the graviton bremsstrahlung signature of leptogenesis by considering a non-thermal production of heavy right-handed neutrinos from evaporation of ultra-light primordial black holes. As it requires very heavy RHNs close to the GUT scale for the corresponding GW spectrum to lie in the most optimistic regime of future verifications, non-thermal origin of RHNs avoid the strong washout issues of thermal leptogenesis at very high scale. In addition, the PBH-dominated era also keeps this GW spectrum distinct from the one produced by the thermal background keeping future experimental verification more plausible. We compute the graviton bremsstrahlung contribution to the GW spectra in this setup and compare it with the GW spectrum from thermal background by considering a few benchmark choices of PBH and model parameters consistent with light neutrino mass and successful leptogenesis. Compared to the other avenues of non-thermal leptogenesis, the presence of ultra-light PBH also brings two additional sources of GW due to density perturbations and graviton emission via Hawking evaporation. The combined GW spectrum including all these sources bears a multi-peak feature at frequencies above the MHz regime. For resonantly enhanced CP asymmetry, the spectrum due to graviton bremsstrahlung becomes more prominent and further separated from the thermal contribution compared to the hierarchical leptogenesis scenario. While most of the near future GW prospects lie in the low frequency regime (\(f \lesssim \mathcal{O}(\rm kHz)\)), our work shows one possible example of high-frequency GW with unique features and connected to the origin of light neutrino mass and baryon asymmetry of the Universe. More avenues of high-frequency GW detection should be able to verify these scenarios to a greater accuracy.

The work of D.B. is supported by the Science and Engineering Research Board (SERB), Government of India grant CRG/2022/000603.

5 Details of the benchmark points↩︎

The Yukawa coupling matrix \(Y_\nu\) in Eq. 5 can be conveniently expressed using the Casas–Ibarra parametrization [113], \[\begin{align} \label{eq:Ynu} Y_{\nu}=-i\,\frac{\sqrt{2}}{v}\,\mathcal{U}\, D_{\sqrt{m}}\,\mathcal{R}^{T}\,D_{\sqrt{M}}\,, \end{align}\tag{19}\] where \(\mathcal{U}\) is the Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix [2], which relates the flavour and mass eigenstates of the light neutrinos. The diagonal matrices \[\begin{align} D_{\sqrt{m}} &=\text{diag}\!\left(\sqrt{m_1},\,\sqrt{m_2},\,\sqrt{m_3}\right),\\ D_{\sqrt{M}} &=\text{diag}\!\left(\sqrt{m_{N_{1}}},\,\sqrt{m_{N_{2}}},\,\sqrt{m_{N_{3}}}\right), \end{align}\] contain the square roots of the light-neutrino and RHN masses, respectively, where the RHNs are assumed to be in a mass-diagonal basis. The matrix \(\mathcal{R}\) is a complex orthogonal matrix satisfying \(\mathcal{R}^{T}\mathcal{R}=\mathbb{I}\) and can be parametrised in terms of three complex angles \(\omega_{1,2,3}\).

Below we list the details of the benchmark points BP1, BP2, BP3 used in our analysis.

1. BP1 : \(M_{\rm in}=0.1\) g, \(\beta=5\times 10^{-4}\), \(m_{N_1} = 10^{15}\) GeV, \(m_{N_2} = 7\times10^{15}\) GeV, \(m_{N_3} = 10^{16}\) GeV, \(m_{1}=10^{-10}\) eV, \(\omega_{1} = 0\), \(\omega_{2}= 0.00056 - i\,0.00039\), \(\omega_{3} = 0.0055 - i\,0.0011\).
   The corresponding Yukawa matrix is given by
   \(Y_{\nu} = \begin{pmatrix} 1.59\times10^{-3} + 4.31\times10^{-4}\,i & 7.89\times10^{-1} - 4.91\times10^{-6}\,i & -6.11\times10^{-1} + 1.56\times10^{-1}\,i \\[2mm] 2.46\times10^{-3} + 1.43\times10^{-4}\,i & 9.85\times10^{-1} + 1.91\times10^{-2}\,i & 2.49 + 2.23\times10^{-6}\,i \\[2mm] -7.32\times10^{-4} - 5.97\times10^{-4}\,i & -6.28\times10^{-1} + 2.41\times10^{-2}\,i & 3.14 - 5.32\times10^{-7}\,i \end{pmatrix},\)
   with \(\sqrt{(Y^{\dagger}_{\nu}Y_{\nu})_{11}} = 10^{-5}\). The CP asymmetry parameter is \(\epsilon^{\Delta L}_{1} =0.001\).

2. BP2 : \(M_{\rm in}=0.1\) g, \(\beta=5\times10^{-4}\), \(m_{N_1} = 10^{14}\) GeV, \(m_{N_2} = 7\times10^{15}\) GeV, \(m_{N_3} = 10^{16}\) GeV, \(m_{1}=10^{-6}\) eV, \(\omega_{1} = 0\), \(\omega_{2}= 0.004 - i\,10^{-4}\), \(\omega_{3} = 0.012 + i\,10^{-3}\).
   The corresponding Yukawa matrix is given by
   \(Y_{\nu} = \begin{pmatrix} 2.36\times10^{-3} + 1.63\times10^{-4}\,i & 7.89\times10^{-1} - 2.18\times10^{-5}\,i & -6.11\times10^{-1} + 1.56\times10^{-1}\,i \\[2mm] 1.74\times10^{-3} + 1.56\times10^{-4}\,i & 9.85\times10^{-1} + 1.91\times10^{-2}\,i & 2.49 - 5.50\times10^{-6}\,i \\[2mm] 1.17\times10^{-3} - 2.65\times10^{-5}\,i & -6.28\times10^{-1} + 2.41\times10^{-2}\,i & 3.14 + 9.74\times10^{-7}\,i \end{pmatrix},\)
   with \(\sqrt{(Y^{\dagger}_{\nu}Y_{\nu})_{11}} = 10^{-5}\). The CP asymmetry parameter is \(\epsilon^{\Delta L}_{1} =10^{-5}\).

3. BP3 : \(M_{\rm in}=1\) g, \(\beta=5\times10^{-4}\), \(m_{N_1} = 10^{14}\) GeV, \(m_{N_2} = 7\times10^{15}\) GeV, \(m_{N_3} = 10^{16}\) GeV, \(m_{1}=10^{-6}\) eV, \(\omega_{1} = 0\), \(\omega_{2}= 0.004 - i\,10^{-4}\), \(\omega_{3} = 0.012 + i\,10^{-3}\).
   The corresponding Yukawa matrix is given by
   \(Y_{\nu} = \begin{pmatrix} 3.54\times10^{-4} + 4.26\times10^{-4}\,i & 7.89\times10^{-1} - 1.80\times10^{-5}\,i & -6.11\times10^{-1} + 1.56\times10^{-1}\,i \\[2mm] 5.24\times10^{-4} + 8.48\times10^{-4}\,i & 9.85\times10^{-1} + 1.91\times10^{-2}\,i & 2.49 - 1.12\times10^{-5}\,i \\[2mm] 1.56\times10^{-4} - 4.74\times10^{-5}\,i & -6.28\times10^{-1} + 2.41\times10^{-2}\,i & 3.14 + 2.05\times10^{-6}\,i \end{pmatrix},\)
   with \(\sqrt{(Y^{\dagger}_{\nu}Y_{\nu})_{11}} = 10^{-6}\). The CP asymmetry parameter is \(\epsilon^{\Delta L}_{1} =3\times10^{-3}\).

6 GW spectra in the absence of PBH↩︎

Fig. 11 shows the graviton bremsstrahlung and CGMB spectra for a few benchmark points without considering a PBH-dominated phase. As mentioned before, for thermally produced RHNs, the peak of the graviton bremsstrahlung spectra remains underneath the one generated by the thermal background.

References↩︎