Primordial black holes and scalar-induced gravitational waves in radiative hybrid inflation

Abstract

We study the possibility that primordial black holes (PBHs) can be formed from large curvature perturbations generated during the waterfall phase transition due to the effects of one-loop radiative corrections of Yukawa couplings between the inflaton and a dark fermion in a non-supersymmetric hybrid inflationary model. We obtain a spectral index $n_s$, and a tensor-to-scalar ratio r, consistent with the current Planck data. Our findings show that the abundance of PBHs is correlated to the dark fermion mass $m_N$ and peak in the GW spectrum. We identify parameter space where PBHs can be the entire dark matter (DM) candidate of the universe or a fraction of it. Our predictions are consistent with any existing constraints of PBH from microlensing, BBN, CMB, etc. Moreover, the scenario is also testable via induced gravitational waves (GWs) from first-order scalar perturbations detectable in future observatories such as LISA and ET. For instance, with inflaton mass $m \sim 2\times 10^{12}$ GeV, $m_N \sim 5.4\times 10^{15}$ GeV, we obtain PBHs of around $10^{-13}\, M_\odot$ mass that can explain the entire abundance of DM and predict GWs with amplitude $\Omega _\textrm{GW}h^2$ $\sim 10^{-9}$ with peak frequency f $\sim$ 0.1 Hz in LISA. By explicitly estimating fine-tuning we show that the model has very mild tuning. We discuss successful reheating at the end of the inflationary phase via the conversion of the waterfall field into standard model (SM) particles. We also briefly speculate a scenario where the dark fermion can be a possible heavy right-handed neutrino (RHN) which is responsible for generating the SM neutrino masses via the seesaw mechanism. The RHN can be produced due to waterfall field decay and its subsequent decay may also explain the observed baryon asymmetry in the universe via leptogenesis. We find the reheat temperature $T_R\lesssim 5\times 10^9$ GeV that explains the matter-anti-matter asymmetry of the universe.

1 Introduction

The seeds of the structure of the universe can be generated by quantum fluctuations of inflaton or curvaton during inflation. The amplitude of curvature perturbations is of the order of $10^{-5}$ at the cosmic microwave background (CMB) scale [1], whereas larger curvature perturbations may be generated at a smaller scale [2,3,4,5,6,7]. Observations of the supermassive  [8, 9] and stellar-mass black hole (BH) merger events by gravitational wave (GW) detectors [10, 11] imply the existence of primordial black holes (PBHs). These PBHs are generated via the collapse of high-dense regions [12,13,14].¹ For other mechanism of PBH and GW production see Refs.  [26,27,28]. The PBH is also a candidate for dark matter (DM) [29] if its mass is within $10^{17\text {--}23}$g (see, e.g., Refs. [30,31,32]). Such large curvature perturbations can be generated if the inflaton or a spectator field goes through a very flat potential during inflation (see Ref. [33] for a recent review on PBHs).

However, going beyond the single-field inflationary scenarios which involve some degrees of fine-tuning [34], one may also generate large curvature perturbations in the hybrid inflation models. Here, inflation ends with a waterfall phase transition [35] with the natural possibility of being embedded in grand unification schemes. Particularly, it is well known that the Grand Unified Theory (GUT) scale $M_\textrm{GUT}$ for such waterfall phase transition corresponds to the observed amplitude of primordial scalar fluctuations as seen in the CMB [36].

Well studied in the literature, the waterfall field potential may be flat leading to the generation of large curvature perturbations at small scales, since the waterfall transition happens at the later stages of inflation [37,38,39]. These large density perturbations which subsequently lead to large curvature perturbations may collapse into a PBH and give us observable effects. Particularly following Refs. [40, 41], one may understand that more often than not it led to the overproduction of PBHs of astrophysical size.² However, with a slightly modified waterfall field potential, this can be avoided and PBH abundance can be controlled with detectable scalar-induced GW signals in the next generation GW observatories [45,46,47]. The models discussed in  [40, 41, 48] also suffer with the problem of initial conditions that need to be set to satisfy current constraints. This nonetheless can be avoided by considering $\alpha$-attractor models [49]. Since there are several large uncertainties present in the critical process of PBH formation estimates, we study it in less detail. On the other hand, we study in detail the characterization of the stochastic gravitational wave background (SGWB) induced at the second-order by the large curvature perturbations (at small length scales) during horizon re-entry in the radiation dominated era. These GWs will be used to test the model involving dark fermions and seesaw physics involving radiative corrections to the hybrid inflation potential. Careful attention has been paid to the consistency with large-scale measurements of the CMB radiation anisotropy from the latest Planck/BICEP/Keck Array release [50]. This involves producing a peak at the small scale, but keeping the amplitude of the power spectrum small at the CMB scale. A best Planck fit spectral index is obtained without the problem of fine-tuning, which is usually the case for single-field inflationary scenarios. We show from our fine-tuning estimates involving the parameter space of dark fermion masses and quartic couplings that the fine-tuning in the model is mild and largely reduced from that of the single-field scenarios by explicitly estimating the fine-tuning quotient.

We investigate a non-supersymmetric particle physics framework where the inflaton field is responsible for dark fermion mass generation [51]. Because of the Yukawa coupling between the two, the inflaton field receives radiative quantum corrections which plays a vital role in setting the inflaton dynamics such that the power spectrum gets enhanced at smaller scales and sets the amplitude of the curvature spectrum at the CMB scale. However, in the hybrid inflationary scenario, there is further growth of perturbation during the waterfall transition which leads to the formation of PBH that can be a dark matter candidate and also induces tensor perturbations at second order propagating as SGWB. We show the second-order tensor perturbations propagating as GWs with amplitude $\Omega _\textrm{GW}h^2$ $\sim 10^{-9}$ and peak frequency f $\sim$ 0.1 Hz for LISA and $\Omega _\textrm{GW}h^2 \sim 10^{-11}$ and peak frequency of $\sim$ 10 Hz in ET. Production of PBH of mass around $10^{-13} M_\odot$ explains the entire abundance of DM in the universe that corresponds to the inflaton mass $m \sim 2\times 10^{12}$ GeV and dark fermion mass $m_N \sim 5.4\times 10^{15}$ GeV. This novel DM candidate is also a signature of the scale of dark fermion physics involving inflationary cosmology. We furthermore speculate that the dark fermions can be the possible candidate for a right-handed neutrino (RHN) responsible for generating SM neutrino mass via the seesaw mechanism and leaving imprints in the scalar power spectrum and GWs. This RHN may also explain the baryon asymmetry of the universe through baryogenesis via leptogenesis.

This paper is structured as follows: In Sect. 2, we present $\alpha$-attractor radiative hybrid inflation with dark fermions. In Sect. 3, we present the power spectrum of scalar perturbations and formation mechanism for PBHs. In Sect. 4, we present the power spectrum for scalar-induced GWs and calculate the signal-to-noise ratio for future planned detectors and correlation between the model parameters. In Sect. 5, we give the fine-tuning estimate and compare it with a single field and standard hybrid inflation and discuss the reheating. In Sect. 6, we explain that our proposed dark fermion can also be a RHN, we discuss the seesaw mechanism, reheating and non-thermal leptogenesis and conclude in Sect. 7. The paper also contains Appendix A where we briefly discuss our model in the Pulsar Timing Arrays (PTAs) band.

2 $\alpha$-Attractor radiative hybrid inflation with dark fermions

We explore the radiative hybrid inflation in the context of an $\alpha$-attractor model where the kinetic part of the inflaton field $\phi$ is modified.³ Working with Einstein–Hilbert gravity action, the relevant terms for the inflaton field $\phi$, the waterfall field $\psi$ and a dark fermion N in the Lagrangian density can be written as,

$$\begin{aligned} \mathscr {L}\simeq&\dfrac{\left( \partial ^\mu \phi \right) ^2}{2\left( 1-\dfrac{\phi ^2}{6\,\alpha }\right) ^2}+\dfrac{\left( \partial ^\mu \psi \right) ^2}{2}+\dfrac{i}{2}\bar{N}\gamma ^\mu \partial _\mu N\nonumber \\&\quad -\kappa ^2\left( M^2-\dfrac{\psi ^2}{4}\right) ^2-\dfrac{1}{2}\,m^2\,\phi ^2\nonumber \\&-\dfrac{\lambda ^2}{4}\,\phi ^2\,\psi ^2 -\dfrac{1}{2}y\,\phi \,\bar{N}\,N-\dfrac{1}{2}Y\,\psi \,\bar{N}\,N-\dfrac{1}{2}m_N\,\bar{N}\,N. \end{aligned}$$

(2.1)

The Yukawa interaction between $\phi$ and N results in one-loop radiative corrections. The hybrid inflationary potential including one-loop radiative corrections and a linear term arising from the interaction between $\psi$ and N can be written as:

$$\begin{aligned} V(\psi ,\phi )&= \kappa ^2\left( M^2-\dfrac{\psi ^2}{4}\right) ^2+\dfrac{1}{2}\,m^2\,\phi ^2+\dfrac{\lambda ^2}{4}\,\phi ^2\,\psi ^2\nonumber \\&\quad +b^3\, \psi + V_\text {loop}(\phi ), \end{aligned}$$

(2.2)

where the one-loop radiative correction is given according to Ref. [51],

$$\begin{aligned} V_\text {loop}(\phi )&=\dfrac{1}{64\,\pi ^2}\left[ m^4\,\text {ln}\left( \dfrac{m^2}{\mu ^2}\right) +\dfrac{\lambda ^4}{4}\left( \phi ^2-\phi _c^2\right) ^2\,\right. \nonumber \\&\quad \times \left. \text {ln}\left( \dfrac{\lambda ^2|\phi -\phi _c|^2}{2\,\mu ^2}\right) \right. \nonumber \\&\quad \left. -2\,\left( m_N+y\,\phi \right) ^4\,\text {ln}\left( \dfrac{m_N+y\,\phi }{\mu }\right) ^2\right] . \end{aligned}$$

(2.3)

Here, M, m (inflaton mass) and b, are the dimensional mass parameters, while $\kappa$, $\lambda$ and y are dimensionless couplings. The parameters $\mu$ and $m_N$ are the cutoff scale and the dark fermion mass, respectively. The linear term in $\psi$ in the potential Eq. (2.2) can arise in many instances, e.g., from the non-perturbative generation of condensates with respect to some dark sector interactions [52,53,54,55,56], or due to gravitational interactions [57, 58], or from specific SUSY related motivations [36, 59,60,61,62,63] and from several other scenarios involving instanton physics [64,65,66]. Here, we do not go into such detailed scenarios but instead, use it as a phenomenological term since our analysis remains applicable to all such scenarios. Moreover, the linear term does not affect the PBH and GW cosmology, except to not overproduce the PBH. We will see later that it can be tuned accordingly to choose the relevant benchmark points (BPs) that we study.

Along the valley, the waterfall field stabilizes at $\psi =0$ as long as $\phi$ is larger than the critical field value $\phi _c$, see Fig. 1. After $\phi$ crosses $\phi _c$ the field $\psi$ falls into one of the two minima of the potential at $\psi \simeq \pm M$ depending on the sign of the coefficient of the linear term b. The radiative hybrid potential Eq. (2.2), in terms of the canonically normalized inflaton field [49], $\phi \rightarrow \sqrt{6\,\alpha } \text {Tanh}\left( \varphi /\sqrt{6\,\alpha }\right)$, can now be written as,

$$\begin{aligned} V(\psi ,\varphi )&= \kappa ^2\left( M^2-\dfrac{\psi ^2}{4}\right) ^2+\dfrac{1}{2}\,\left( m^2+\dfrac{\lambda ^2}{2}\psi ^2\right) \,\nonumber \\&\quad \times \left( \sqrt{6\,\alpha }\, \text {Tanh}\left( \dfrac{\varphi }{\sqrt{6\,\alpha }}\right) \right) ^2 +b^3\, \psi \nonumber \\&\quad +V_\text {loop}\left( \sqrt{6\,\alpha }\, \text {Tanh}\left( \varphi /\sqrt{6\,\alpha }\right) \right) . \end{aligned}$$

(2.4)

Fig. 1

A schematic picture of radiative hybrid inflation potential. The light orange bullet points show the inflationary trajectory of the inflaton $\varphi$, the white bullet is the critical point and the dark orange bullets show the waterfall regime. Inflation continues for some e-folds after the critical point during which very large curvature perturbations are generated

Table 1 Benchmark points (BPs) for model parameters

The schematic view of hybrid potential is shown in Fig. 1. The mass squared of the waterfall field at $\psi =0$ is,

$$\begin{aligned} M_\psi ^2=\left( -\kappa ^2\,M^2+\dfrac{1}{2}\left( \lambda \,\sqrt{6\,\alpha }\,\text {Tanh}\left( \dfrac{\varphi }{\sqrt{6\,\alpha }}\right) \right) ^2\right) . \end{aligned}$$

(2.5)

In this paper, we assume that $(\lambda \,\sqrt{6\,\alpha })^2/2 > \kappa ^2\,M^2$ such that $M_\psi ^2>0$ at large $\varphi >\varphi _c$ to stabilize the inflationary trajectory at $\psi = 0$, where,

$$\begin{aligned} \text {Tanh}^2\left( \dfrac{\varphi _c}{\sqrt{6\,\alpha }}\right) =\dfrac{\kappa ^2\,M^2}{3\,\alpha \,\lambda ^2}. \end{aligned}$$

(2.6)

During inflation, as long as $\varphi \lesssim \varphi _c$, the effective mass square of $\psi$ becomes negative which gives rise to tachyonic instability that will grow the curvature perturbations. These growing perturbations will enhance the scalar power spectrum at small scales and upon horizon re-entry the collapse of large density fluctuations produces the PBHs. Due to the linear term in the potential Eq. (2.4) the field $\psi$ will not relax exactly at $\psi =0$ but will be displaced depending upon the sign. of the coefficient of the linear term b. In this way, we can inflate the unnecessary topological defects and control the peak of the power spectrum at small scales to avoid PBH overproduction.

The slow-roll parameters are given by, [67],

$$\begin{aligned} \epsilon _V&= \dfrac{m_\text {Pl}^2}{2}\left( \dfrac{\partial _\xi V}{V}\right) ^2,\,\,\,\,\,\eta _V = m_\text {Pl}^2\left( \dfrac{\partial _\xi ^2 V}{V}\right) , \nonumber \\ \delta _V^2&= m_\text {Pl}^4 \left( \dfrac{\partial _\xi V\,\partial _\xi ^3 V}{V^2}\right) ,\,\,\,\,\, \sigma _V^3= m_\text {Pl}^6 \left( \dfrac{(\partial _\xi V)^2\,\partial _\xi ^4 V}{V^3}\right) . \end{aligned}$$

(2.7)

Here, $m_\text {Pl}\simeq 2.43 \times 10^{18}$ GeV is the reduced Planck mass and the subscript $\xi =\{\varphi ,\psi \}$ indicates the field derivative. In the slow roll limit, along the valley at the pivot scale, the spectral index $n_s$, its running and running of the running, are given by [67],

$$\begin{aligned} n_s&= 1-6\,\epsilon _V + 2\,\eta _V,\,\,\,\,\, \text {d}n_s/\text {d}\ln k=16\,\epsilon _V\,\eta _V\nonumber \\&\quad -24\,\epsilon _V^2-2\,\delta _V^2,\nonumber \\ \text {d}^2n_s/\text {d}\ln k^2&=-192\,\epsilon _V^3+192\,\epsilon _V^2\,\eta _V-32\,\epsilon _V\,\eta _V^2\nonumber \\&\quad -24\,\epsilon _V\delta _V^2+2\,\eta _V\delta _V^2+2\sigma _V^3. \end{aligned}$$

(2.8)

Fig. 2

Tensor-to-scalar ratio r vs. scalar spectral index $n_s$ for the corresponding parameter sets given in Table 1. The solid contours are the current Planck bounds [68], Planck/BICEP [1, 50, 69] and the dashed shaded region indicates the future proposed experiments (LiteBIRD, CMB S4-Euclid, Simons Observatory (SO)) [70,71,72]

The central measurements by Planck 2018 [1] in the $\Lambda$CDM model are; $n_s = 0.9647 \pm 0.012$, $\text {d}n_s/\text {d}\ln k = 0.0011 \pm 0.0099$ and $\text {d}^2n_s/\text {d}\ln k^2 = 0.009 \pm 0.012$. The tensor to scalar ratio $r=16\,\epsilon _V < 0.036$ at $95\%$ C.L. All these values are measured at the pivot scale, $k_\star =0.05\,\text {Mpc}^{-1}$. For BP-1 in Table 1, we find $\text {d}n_s/\text {d}\ln k \simeq -0.0001373$ and $\text {d}^2n_s/\text {d}\ln k^2\simeq 0.00002961$. The prediction of $n_s$ and r given in Table 1 is shown in Fig. 2 within $1\sigma$ bound of recent Planck results.

To avoid eternal inflation, M should be O(1) [45]. The parameter m along with the loop corrections controls the amplitude of the plateau in the valley. We choose $y > rsim \lambda /\sqrt{2}$ and set $m_N$ as a free parameter such that the dominant contribution to the potential comes from the last term in Eq. (2.3) thus making it radiatively generated. We set the normalization scale $\mu$ so the log term remains positive. The coupling $\lambda$ determines the number of e-folds in the waterfall regime and $\kappa$ will fix the amplitude of the power spectrum at the pivot scale, that is, $A_s(k_\star )\simeq 2.24\times 10^{-9}$. Taking into account all these constraints we define the BPs in Table 1 for the potential in Eq. (2.4) and study the production of PBH and the scalar-induced SGWB in the subsequent sections. The subscripts i and e in Table 1 refer to the values at the start and end of inflation. The presence of the linear term makes the potential unbounded close to the critical point since the coefficient b is very small and for larger field values it does not play a significant role.

3 Scalar perturbations and primordial black hole formation

Let us explain the generation and evolution of scalar perturbations during inflation.

3.1 Scalar spectra

The background equations of motion in the number of e-fold times are given by [73],

$$\begin{aligned}&\varphi ^{''}+\left( \dfrac{H^{'}}{H}+3\right) \varphi ^{'}+\dfrac{V_\varphi }{H^2}=0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\nonumber \\&\psi ^{''}+\left( \dfrac{H^{'}}{H}+3\right) \psi ^{'}+\dfrac{V_{\psi }}{H^2}=0, \end{aligned}$$

(3.1)

where H, the Hubble rate is defined to be, $H^2=2\,V/(6-\varphi ^{'2}-\psi ^{'2})$. Here, prime is the derivative with respect to the number of e-folds and $V_\xi =dV/d\xi$ where $\xi =\{\varphi ,\psi \}$. The evolution of the fields is shown in Fig. 3, where the background is calculated until the end of inflation, that is, when $\epsilon _H\equiv -H^{'}/H=1$. For different BPs given in Table 1, the critical point shifts and affects the number of e-folds in the waterfall regime. Usually, the evolution of quantum fluctuations is affected due to the back-reaction of long modes. Although, for detailed analysis one needs to perform the lattice simulations which is beyond the scope of the present work. As a naive criterium, which should anyway give a good estimate of when back-reaction is negligible, we may see that the background energy density; $H^2\,\phi ^{'2}/2+m^2\,\phi ^2/2\ll M^4$, for the BPs we chose in Table 1. This enables us to disregard the backreaction on metric as discussed in detail in [3, 45]

The scalar perturbations of the Friedmann–Lemaître–Robertson–Walker (FLRW) metric in longitudinal gauge can be expressed as [73];

$$\begin{aligned} ds^2&=a(\tau )^2\left[ (1+2\,\Phi _{\text {B}})d\tau ^2\right. \nonumber \\&\quad \left. +\left[ (1-2\,\Psi _{\text {B}})\delta _{ij}+\dfrac{h_{ij}}{2}\right] dx^idx^j\right] , \end{aligned}$$

(3.2)

where $\Phi _{\text {B}}$ and $\Psi _{\text {B}}$ are the Bardeen potentials, $h_{ij}$ is the transverse-traceless tensor metric perturbation i.e. $h_i^i= 0 = h^j _{i,j}$. The conformal time $\tau$, is related to cosmic time, $dt=a\,d\tau$, a being the scale factor. Working in the conformal Newtonian gauge we set, $\Phi _{\text {B}}=\Psi _{\text {B}}$. Following the dynamics given in [73, 74] to calculate the power spectrum, the Klein–Gordon equation to evaluate scalar perturbations is given by,

$$\begin{aligned}&\delta \xi _i^{''}+(3-\epsilon )\delta \xi _i^{'}+\sum _{j=1}^{2}\dfrac{1}{H^2}V_{\xi _i \xi _j}\delta \xi _j+\dfrac{k^2}{a^2H^2}\delta \xi _i\nonumber \\&\quad =4\Phi ^{'}_{\text {B}}\,\xi ^{'}_i-\dfrac{2\,\Phi _{\text {B}}}{H^2}V_{\xi _i}. \end{aligned}$$

(3.3)

Here $\xi$ with subscript (i, j) refers to the fields ($\varphi$, $\psi$), k is the comoving wave vector, the equation of motion for $\Phi _{\text {B}}$ is given by,

$$\begin{aligned} \Phi ^{''}_{\text {B}}+(7-\epsilon )\,\Phi ^{'}_{\text {B}}+\left( \dfrac{2\,V}{H^2}+\dfrac{k^2}{a^2H^2}\right) \Phi _{\text {B}}+\dfrac{V_{\xi _i}}{H^2}\,\delta \xi _i=0. \end{aligned}$$

(3.4)

Fig. 3

Field and the slow-roll parameter $\epsilon _H\equiv -H^{'}/H$ evolution with the number of e-folds from pivot scale to the end of inflation. We evaluate solving the full background Eq. (3.1) using potential in Eq. (2.4) for the BPs in Table 1

The initial conditions (i.c) for field perturbations in e-fold time are given as,

$$\begin{aligned} \delta \xi _{i,\text {i.c}}{=}\,\dfrac{1}{a_{\text {i.c}}\sqrt{2 k}},\,\,\,\,\,\,\,\,\,\,\,\,\,\, \delta \xi _{i,\text {i.c}}^{'}{=}\,{-}\dfrac{1}{a_{\text {i.c}}\sqrt{2 k}} \left( 1+\iota \dfrac{k}{a_{\text {i.c}}H_{\text {i.c}}}\right) . \end{aligned}$$

(3.5)

The initial conditions for the Bardeen potential and its derivative are given by,

$$\begin{aligned} \Phi _{\text {B,i.c}}&{=} \sum _{i=1}^{2}\dfrac{\left( H_{\text {i.c}}^2 \xi _{i,\text {i.c}}^{'}\delta \xi _{i,\text {i.c}}^{'}{+}\left( 3H_{\text {i.c}}^2 \xi _{i,i.c}^{'}{+}V_{\xi _i,\text {i.c}}\right) \delta \xi _{i,\text {i.c}}\right) }{2H_{\text {i.c}}^2\left( \epsilon _{\text {i.c}}{-}\dfrac{k^2}{a_{\text {i.c}}^2H^2_{\text {i.c}}}\right) },\,\,\,\,\nonumber \\ \Phi ^{'}_{\text {B,i.c}}&=\sum _{i=1}^{2}\dfrac{\xi _{i,\text {i.c}}^{'}\delta \xi _{i,\text {i.c}}}{2}-\Phi _{\text {B,i.c}}. \end{aligned}$$

(3.6)

The scalar power spectrum $P_s(k)$ is defined as,

$$\begin{aligned} P_s(k)=\dfrac{k^3}{2\pi ^2}\left| \zeta \right| ^2=\dfrac{k^3}{2\pi ^2}\left| \Phi _\mathrm{{B}}+\dfrac{\sum _{i=1}^{2}\xi ^{'}_i\delta \xi _i}{\sum _{j=1}^{2}\xi ^{'2}_j}\right| ^2. \end{aligned}$$

(3.7)

Fig. 4

Power spectrum from pivot scale to the end of inflation by solving the exact scalar perturbation equations for BP given in Table 1. The shaded region corresponds to the constraints from the present (solid) and future (dashed) experiments, see the main text for the details. The orange line in PBH overproduction bound represents the window where PBHs can explain the DM in totality

In this paper, we present an exact power spectrum (numerically solved) for the potential Eq. (2.4) and the the resulting spectrum is shown in Fig. 4 from the pivot scale to the end of inflation. One can see from Fig. 4 that for a certain choice of the inflaton mass which is large (see BP-1 in Table 1), it ruins the power spectrum predictions and runs into severe constraints at the CMB measurement. However, for this same choice, if the inflaton is coupled to the dark fermion N, then the quantum radiative loop corrections rectify the power spectrum and give a Planck consistent $n_s\simeq 0.964$ at the CMB scale. For BP-2, without loop corrections, $n_s\simeq 0.974$ but with loop corrections, we obtain $n_s=0.969$ which is consistent with the measurements, see the inset of the bottom right in Fig. 4. Note that the duration of the waterfall is controlled by vev M. To avoid an abrupt end or a longer period of inflation, M should be of the order 1 in Planckian units. As the value of M increases, we expect the duration of the waterfall transition to increase, which enhances the curvature perturbation $P_{s}(k)$. As in Fig. 4, the top right panel clearly depicts that the larger M enhances the waterfall phase and the peak of the spectrum is shifted towards the left and the smaller value reduces the waterfall phase and the peak of the spectrum is moved towards the right.

We illustrate for comparison the current bounds from Planck [68], Lyman-alpha [75], PIXIE [76], COBE/Firas [77], PTA [78], the PBH overproduction that is related to the critical density $\delta _c$ which we discuss later and a window where PBH as the entire DM candidate of the universe (see orange line in the PBH overproduction bound in Fig. 4). At scales $10^{-4}\lesssim k/\text {Mpc}^{-1}\lesssim 1$, the power spectrum is constrained by the angular resolution of the current CMB measurements. However, inhomogeneities at these scales result in isotropic deviations from the usual black body spectrum and are known as CMB spectral distortions [79]. These distortions are usually categorized into two: $\mu$-distortions, associated with chemical potential that occur at early times, and Compton y-distortions, generated at redshifts $z \lesssim 5 \times 10^4$. A $\mu$-distortion is associated with a Bose-Einstein (BE) distribution with $\mu \ne 0$. The most stringent constraints on spectral distortions currently come from observations made by the COBE/FIRAS experiment, which restricts $|\mu |\lesssim 9.0 \times 10^{-5}$ and $|y |\lesssim 1.5\times 10^{-5}$ at $95\%$ CL [77]. In future, a PIXIE-like detector can investigate distortions with magnitudes $\mu \lesssim 2\times 10^{-8}$ and $y\lesssim 4 \times 10^{-9}$ [76]. We also present acoustic reheating (AR) constraints on the spectrum [80]. Note the solid lines are present and dashed are the future experiments. We find that the model parameters as shown in Table 1 satisfy the constraints of COBE/FIRAS. However, for the broad-width power spectrum, to explain the NANOGrav signal [81, 82], one not only needs to assume a diversion from the $\Lambda$CDM model but also the power spectrum conflicts with COBE/Firas. For further discussion on this issue, see Appendix A.

3.2 Primordial black hole formation

The mass of the PBH (in solar mass units $M_\odot$) formation is associated with a wave vector k and is given by [83],

$$\begin{aligned} M_{\text {PBH}}=3.68\left( \dfrac{\gamma _c}{0.2}\right) \left( \dfrac{g_{*}(T_f)}{106.75}\right) ^{-1/6}\left( \dfrac{10^6\,\text {Mpc}^{-1}}{k}\right) ^2\,M_{\odot }. \end{aligned}$$

(3.8)

The fractional abundance of PBHs, $\Omega _{\text {PBH}}/\Omega _{\text {DM}}\equiv f_{\text {PBH}}$ is defined in terms of PBH mass [83],

$$\begin{aligned} f_{\text {PBH}}&= \dfrac{\beta (M_{\text {PBH}})}{3.94\times 10^{-9}}\left( \dfrac{g_{*}(T_f)}{106.75}\right) ^{-1/4}\left( \dfrac{\gamma _c}{0.2}\right) ^{1/2}\nonumber \\&\quad \times \left( \dfrac{0.12}{\Omega _{\text {DM}} h^2}\right) \left( \dfrac{M_{\text {PBH}}}{M_{\odot }}\right) ^{-1/2}, \end{aligned}$$

(3.9)

where $h^2\Omega _{\text {DM}}=0.12$ is the current energy density of DM, $\gamma _c=0.2$ is the factor that depends on the gravitation collapse.⁴ Using the Press Schechter approach,⁵$\beta$ the fractional energy density at the time of formation is given by [88],

$$\begin{aligned} \beta (M_{\text {PBH}})=\dfrac{1}{2\pi \sigma ^2(M_{\text {PBH}})}\int _{\delta _c}^{\infty }d\delta \,\,\, \text {exp}\left( -\dfrac{\delta ^2}{2\sigma ^2\left( M_{\text {PBH}}\right) }\right) . \end{aligned}$$

(3.10)

Fig. 5

PBH abundance as DM given in Eq. (3.9). The shaded regions represent the observational constraints on the PBH abundance from various experiments (solid lines present and dashed for future), see the main text for the details

Here $\delta _c$ is the critical density perturbation of PBH formation, and it takes the value between $0.4-0.6$ [89,90,91,92,93] that corresponds to $\sigma ^2(M_\text {PBH})\simeq 10^{-2}-10^{-3}$ to explain the full abundance of DM. The variance $\sigma (M_\text {PBH})$, of the curvature perturbations, is given by

$$\begin{aligned} \sigma ^2(M_\text {PBH}(k))=\dfrac{16}{81}\int \dfrac{dk^{'}}{k^{'}}(k^{'}/k)^4W^2(k^{'}/k) P_s(k^{'}),\nonumber \\ \end{aligned}$$

(3.11)

where $W(x)=\text {exp}(-x^2/2)$ is the Gaussian window function. The abundance of PBH Eq. (3.9) is shown in Fig. 5. For given parameter sets in Table 1, the PBHs explain the entire abundance of DM for BP-1. However, for BP-2, some fraction of DM can be explained due to different observational constraints from the experiments, as shown by the shaded region in Fig. 5. As from Eq. (3.8), $M_\text {PBH}\propto k^{-2}$, therefore to explain the PBHs in PTAs the power spectrum has to be enhanced at low k and large $M_\text {PBH}$ which are constrained by the Microlensing effects. Therefore, the entire abundance of DM from PBHs cannot be observed in PTAs.

In Fig. 5 we depict the constraints on $f_\textrm{PBH}$ (see Refs. [94,95,96,97] for details on constraints). Evaporation of PBH via Hawking radiation sets severe constraints in: CMB [98], EDGES [99], INTEGRAL [100, 101], Voyager [102], 511  keV [103], EGRB [104]; HSC (Hyper-Supreme Cam) [105], EROS [106], OGLE [107] and Icarus [108] which are all categorized as micro-lensing related observations, several constraints arising due to modification of the CMB spectrum which happens if PBHs accrete as investigated in Ref. [109] (see also Ref. [110]). Finally, the range around $M_{\odot }$ is constrained by LIGO-VIRGO-KAGRA (LVK) observations on PBH-PBH merger  [111,112,113,114,115,116,117]), while future GW interferometer based detectors like Cosmic Explorer (CE) and Einstein Telescope (ET) are expected to set limits on the PBH abundance as shown in Refs. [118,119,120,121,122,123], these are shown in dashed lines in the plot. We also show future sensitivity reaches of the Nancy Grace Roman Space Telescope (NGRST) from micro-lensing, see Ref. [124]. Our prediction in Table 1, for the BP-1 the PBHs explain the entire abundance of DM. While for BP-2 and 3 some fraction of DM can be explained but can be tested in future experiments like NGRST.

4 Scalar-induced gravitational waves

Assuming the formation of PBHs in the radiation dominated era, the present-day energy density of the GWs in terms of scalar power spectrum Eq. (3.7), is given by [125,126,127,128],

$$\begin{aligned} \Omega _{\text {GW}}(k)&=\dfrac{c_g\,\Omega _{r}}{6}\left( \dfrac{g_{*}(T_f)}{106.75}\right) \int _{-1}^{1}dd\int _{1}^{\infty }ds\,\,\nonumber \\&\quad \times P_s\left( k\dfrac{s-d}{2}\right) P_s\left( k\dfrac{s+d}{2}\right) I(d,s), \nonumber \\ I(d,s)&=\dfrac{288(d^2-1)^2(s^2-1)^2(s^2+d^2-6)^2}{(d-s)^8(d+s)^8}\nonumber \\&\quad \times \left\{ \left( d^2-s^2+\dfrac{d^2+s^2-6}{2}\text {ln}\left| \dfrac{s^2-3}{d^2-3}\right| \right) ^2\right. \nonumber \\&\quad \times \left. \dfrac{\pi ^2}{4}(d^2+s^2-6)^2\Theta (s-\sqrt{3}))\right\} . \end{aligned}$$

(4.1)

In the expression above, $\Omega _r=5.4\times 10^{-5}$ is the present day value of energy density, $c_g=0.4$ in the SM, $\Theta$ is the Heaviside function and the effective degrees of freedom $g_{*}(T_f)$ at the temperature of PBH formation $T_f$ is 106.75 for SM like spectrum. Furthermore, using $k=2\pi f$, $1\text {Mpc}^{-1}=0.97154\times 10^{-14}\,\text {s}^{-1}$ and $h=0.68$, we transformed into the $h^2\Omega _{\text {GW}} (f)-f$ plane. The GW spectra for the BPs in Table 1 are shown in Fig. 6 with the variation of different model parameters as labeled. The model can be tested in future planned experiments like SKA [129], THEIA [130], LISA [131], $\mu$-ARES [132], BBO [133], U-DECIGO [134, 135], CE [136] and ET [137] experiments presented by shaded region in Fig. 6. For the observation of GW spectra in the PTA band, see Appendix A.

Fig. 6

The energy density of GWs for Eq. (4.1) for the BPs given in Table 1. The colored shaded regions indicate the sensitivity curves of present (solid boundaries) LIGO O3 [138], NANOGrav [81] and future (dashed boundaries) LIGO O5, SKA [129], THEIA [130], LISA [131], $\mu$-ARES [132], BBO [133], U-DECIGO [134, 135], CE [136] and ET [137] experiments. The hatched region shows the astrophysical background [139]. See text for details

When looking for stochastic GW background of cosmic origin we expect several astrophysical sources of GW to be present which will be like a background, mainly in the form of LIGO/VIRGO observed binary black hole (BH–BH) merging events [140,141,142,143,144,145,146] and binary neutron star (NS–NS) events [147, 148]. Therefore to distinguish scalar-induced GWs of cosmic origin, the NS and BH foreground may be subtracted with sensitivities of BBO and ET or CE windows, particularly in the range $\Omega _\textrm{GW} \simeq 10^{-15}$ [149] and $\Omega _\textrm{GW} \simeq 10^{-13}$ [150]. In the LISA window, the binary white dwarf galactic and extra-galactic (WD–WD) may be of greater significance than the NS–NS and BH–BH foregrounds [151,152,153] and can be subtracted [154] with the expected sensitivity at $\Omega _\textrm{GW} \simeq 10^{-13}$ [155, 156]. This subtraction procedure alongside the fact that the GW spectrum generated by the astrophysical foreground increased with frequency as $f^{2/3}$ [157] is different from the GW spectrum generated by second-order GWs from radiative hybrid inflation ( at low frequency $f^{3/2}$, at higher frequency $f^{-3/2}$). Now, one will be able to pin down the GW signals from the scalar-induced gravitational wave source from the hybrid inflation as we predict.

4.1 Signal-to-noise ratio

The ability of an interferometer to quantify SGWB signal with an energy density $\Omega _\text {GW}(f)$ for observation time $T_\text {obs}$ is quantified as the signal-to-noise ratio (SNR). To estimate the SNR, we refer to the standard way of computation [158],

$$\begin{aligned} \text {SNR} \equiv \left[ T_\text {obs} \int _{f_\text {min}}^{f_\text {max}} \left( \frac{\Omega _{\text {GW}}(f)}{\Omega _{\text {Noise}}(f)} \right) ^2 d {f} \right] ^{1/2} \, , \end{aligned}$$

(4.2)

where $f_\text {min}$, $f_\text {max}$ denotes the detector bandwidth. We have utilized the noise curve $\Omega _{\text {Noise}}(f)$ for a given experiment and have assumed the duration of each GW mission to be $T_\text {obs}=4$ years. We present the SNR as a function of different model parameters in Fig. 7 and show SNR=10 as a reference with a gray dashed line. We find that most of the future planned experiments LISA, BBO, DECIGO, ET, CE, THEIA and $\mu$-ARES are all capable of testing the model with SNR $> 10$.

Fig. 7

Signal-to-noise ratio (SNR) as a function of various model parameters indicating projections of future sensitivity of different experiments as indicated by the color coding in the inset. The dashed gray line indicates SNR $=10$

Fig. 8

Parameter space in the $\lambda$ vs. $m_N$ with varying SNR are shown by the vertical bar legends for LISA, THEIA, BBO, CE, ET and $\mu$-ARES as labeled

The behavior of SNR can be comprehended from Fig. 6 where we vary one parameter and the rest of them are fixed according to BP-1 in Table 1. Variations of $\lambda$ in the top left plot in Fig. 6 do not have a large impact on the parameter space within the sensitivity curves like LISA, BBO, ET, CE and $\mu$-ARES, this gives a slight increase in SNR by increasing $\lambda$. There is a noticeable decrease in THEIA, since, by increasing $\lambda$, the power spectrum will span less space of the sensitivity curve for THEIA, as shown in the top left plot of Fig. 7. Variation of M in the top right of Fig. 6 predicts less sensitivity for ET, CE, BBO and more sensitivity for LISA $\mu$-ARES, THEIA with increasing M. This behavior is shown in the upper left panel of Fig. 7. In the middle left of Fig. 6, one can see increasing sensitivity for ET and CE and a decrease for LISA, BBO, THEIA and $\mu$-ARES with increasing $m_N$. This is consistent with the SNR shown in the middle left of Fig. 7. Parameter m affects the power spectrum in such a way that if it leaves the sensitivity curve at a lower frequency, it will span a slightly more parameter space at a higher frequency, therefore increasing SNR with increasing m is depicted in the middle right of Fig. 7 except for THEIA. Increasing the coefficient of the linear term b, reduces the peak of the spectrum and therefore SNR will decrease as shown in the third row of Fig. 7.

In Fig. 8, we present the correlation between the parameter space for the model parameters $\lambda$ vs. $m_N$ with the variation in SNR shown by a vertical bar for LISA, THEIA, BBO, CE, ET and $\mu$-ARES as indicated. In Fig. 8, the SNR for LISA and $\mu$-ARES increases with the decreasing $m_N$, and increasing $\lambda$, this feature is depicted in the top left and bottom right of Fig. 8. For ET and CE, increasing both $m_N$ and $\lambda$ increases the SNR; see the top right and bottom left of Fig. 8. For THEIA, the SNR increases with decreasing $m_N$ and $\lambda$, see the upper right of Fig. 8. For BBO, the behavior is different due to its sensitivity frequency range, SNR will increase with increasing $\lambda$ for a particular range of $m_N$, see the middle right of Fig. 8. The correlation parameter space for M vs. $m_N$ for the indicated experiments is given in Fig. 9. For LISA and BBO the SNR increases for some intermediate values of M and $m_N$, for the rest of the parameter space, it decreases. For THEIA, SNR increases with increasing M and for $\mu$-ARES the increasing M and decreasing $m_N$ gives the larger SNR. For ET and CE, the behavior is similar, small M and large $m_N$ increase the SNR; see Fig. 9.

Fig. 9

Parameter space in the M vs. $m_N$ with varying SNR are shown by the vertical bar legends for LISA, THEIA, BBO, CE, ET and $\mu$-ARES as labeled

5 Fine-tuning estimates

The single-field inflationary models, where the PBHs arise due to inflection points, etc. require a high level of fine-tuning of the parameters involved in the enhancement of the power spectrum at small scales [83, 159]. In this section, we show that, unlike single-field inflation, the amount of fine-tuning is much smaller in hybrid inflation. To estimate the required amount of fine-tuning, we calculate a quantity $\Delta _x$ given by

$$\begin{aligned} \Delta _x = \text {Max}\bigg |\dfrac{\partial \ln P_s^\text {Peak}}{\partial \ln x}\bigg |, \end{aligned}$$

(5.1)

where $x\in \{m,M,m_N,\lambda \}$ is the model parameter. Evaluating numerically, a fine-tuning estimate for theory parameters is given in Table 2. The larger the $\Delta _x$ is, the larger the amount of required fine-tuning. Note that we did not consider the case for b, since the variation of b does not keep the height of the power spectrum fixed and therefore affects the abundance of DM.

Table 2 Fine-tuning (FT) estimate of model parameters for BP$-1$ in Table 1 with a peak of the spectrum around $5\times 10^{-3}$

Fig. 10

In the left panel, variation of the field values with respect to the number of e-folds is shown for different initial conditions. In the right panel, we have shown the field value from the horizon exit till the end of inflation

The maximum fine-tuning we obtain is 15, which is almost five orders of magnitude smaller than single-field inflation [48] and one order of magnitude smaller than supersymmetric hybrid inflation [160].⁶

As can be seen in the right panel of Fig. 10, variation of initial conditions has no impact within our constrained parameter space. Due to attractor behavior, the power spectrum is not affected by the initial condition of $\phi$ within a constrained number of e-folds. In the right panel of Fig. 10, we show the field evolution for the total number of e-folds for different initial conditions. In the left panel of Fig. 10, we present the last 57 e-folds from horizon exit till the end of inflation, contributing as a BP-1. Therefore, the initial conditions do not contribute to the fine-tuning estimate within the constrained number of e-folds. For further details on the initial conditions see [45]. Also note that due to the presence of the linear term in Eq. (2.4), $\psi$ will not relax exactly at zero but is slightly displaced depending on the sign of the coefficient b. Therefore, even if $\psi _i=0$, the field will evolve with the number of e-folds and can be seen in Fig. 10.

6 Reheating

Let us now consider the SM Higgs h and its coupling with $\varphi$ and $\psi$ that is conducive to reheating the universe. The potential is written as,

$$\begin{aligned} V&=\lambda _h\left( h^2-\dfrac{v_h^2}{2}\right) ^2+2\,\lambda _{\phi \,h}\,\left( h^2-\dfrac{v_h^2}{2}\right) ^2\phi ^2\nonumber \\&\quad +2\,\lambda _{\psi h}\left( \psi ^2-\dfrac{M^2}{2}\right) \left( h^2-\dfrac{v_h^2}{2}\right) , \end{aligned}$$

(6.1)

where $v_h$ is the vacuum expectation value (VEV) of the SM Higgs, $\lambda _{\phi h}$, $\lambda _h$, $\lambda _{\psi h}$ are dimensionless couplings. For simplicity, we assume ${\lambda _{\phi h},\lambda _{\psi h}}$ to be very tiny such that the inflation and waterfall transition are not affected by the SM Higgs dynamics. However, $\lambda \gg \lambda _{\phi h}$ need not be exactly zero but small, as this will be responsible for the reheating.⁷ Now, once inflation ends, the inflaton will decay and reheat the universe. The inflaton predominantly decays into the SM Higgs and the corresponding reheat temperature is estimated assuming perturbative reheating by,

$$\begin{aligned} T_R\simeq \sqrt{\left( \dfrac{90}{\pi ^2\,g_\star }\right) ^{1/2}\Gamma _\varphi \,m_\text {Pl}}. \end{aligned}$$

(6.2)

Fig. 11

Reheat temperature as a function of $\lambda _{\varphi h}$ is shown in red color and the gray dashed lines are the upper and lower bounds, see the main text for an explanation

Here, $g_\star$ is the effective degrees of freedom and $\Gamma _\varphi (\varphi \rightarrow hh)$ is the decay rate of $\varphi$ given by,

$$\begin{aligned} \Gamma _\varphi \simeq \dfrac{\lambda _{\varphi h}^2}{8\,\pi }m. \end{aligned}$$

(6.3)

The coupling $\lambda _{\varphi h}$ is bounded below from the reheat temperature, that is $T_R > rsim 4$ MeV and above from the assumption $\lambda \gg \lambda _{\varphi h}$ therefore, $10^{-17}\lesssim \lambda _{\varphi h} \lesssim 10^{-9}$ and is shown in Fig. 11. The maximum reheat temperature we acquire is $2\times 10^5$ GeV. We see that even very tiny $\lambda _{\varphi h}$ leads to successful reheating of the universe and this is also consistent with our earlier assumption. Such choices of the very small Higgs-portal coupling allowed for SM Higgs dynamics not to let the SM Higgs acquire too large quantum fluctuation during inflation as explicitly studied in Ref. [161].⁸ For inflaton-Higgs coupling considered in our analysis $\lambda _{\varphi h} \ll 10^{-10}$, non-perturbative preheating particle production is negligible, therefore we consider simple perturbative reheating estimates from inflaton decay [163, 164].

7 Right-handed Neutrino, seesaw and non-thermal leptogenesis

In this section, we speculate a possibility that our proposed dark fermion could be the heavy right-handed neutrino (RHN) which via seesaw generates tiny masses for the SM neutrinos. Considering the relevant portion of Lagrangian, which is responsible for reheating and seesaw,

$$\begin{aligned} \mathscr {L}&\supset -\dfrac{1}{2}y_h\,h\,\bar{v}_L\,N-\dfrac{1}{2}y\,\phi \,\bar{N}\,N-\dfrac{1}{2}Y\,\psi \,\bar{N}\,N\nonumber \\&\quad -\dfrac{1}{2}m_N\,\bar{N}\,N. \end{aligned}$$

(6.4)

SM Higgs h, gives mass to the neutrino through Yukawa coupling $y_h$, with the SM neutrino $v_L$. Accordingly, the neutrinos acquire the mass via the seesaw mechanism [165],

$$\begin{aligned} m_v = \dfrac{y_h^2\,v_h^2}{m_N}, \end{aligned}$$

(6.5)

where $v_h = 174$ GeV is the VEV of the Higgs field.

Consider the BP-3 in Table 1, the effective mass of the waterfall field Eq. (2.5) at the end of inflation is $M_{\psi }=3.3\times 10^{-4}\, m_\text {Pl}> 2\,m_N$. This kinematic bound has to be satisfied for the successful decay of the waterfall field into the RHNs. The waterfall field after the end of inflation will decay into RHNs and the corresponding decay rate is given by,

$$\begin{aligned} \Gamma _\psi (\psi \rightarrow NN)\simeq \dfrac{Y^2}{8\,\pi }M_\psi . \end{aligned}$$

(6.6)

The reheating temperature is given by,

$$\begin{aligned} T_R\simeq \sqrt{\left( \dfrac{90}{\pi ^2\,g_\star }\right) ^{1/2}\Gamma _\psi \,m_\text {Pl}}. \end{aligned}$$

(6.7)

The decay channel $\Gamma _\psi$ is important for successful leptogenesis. Taking into account the non-thermal leptogenesis and assuming instantaneous decay of the waterfall field into RHN and subsequent decays of RHN to SM particles to generate the lepton asymmetry and transfer it to baryon asymmetry via sphaleron [166],

$$\begin{aligned} \dfrac{n_B}{s}&\simeq 4.43\times 10^{-17}\,\text {Br}_{\psi \rightarrow NN}\left( \dfrac{T_R}{1\,\text {GeV}}\right) \left( \dfrac{2\,m_N}{M_\psi }\right) \nonumber \\&\quad \left( \dfrac{m_{\nu _3}}{0.05\,\text {eV}}\right) \delta _\text {eff}. \end{aligned}$$

(6.8)

Here, $m_{\nu _3}\simeq 0.05$ eV and $\delta _\text {eff} \le 1$ is the CP violating phase factor. The branching ratio Br $=1$ if the dominant decay of the waterfall field is into RHNs. For the BP-3 in Table 1 and for observed baryon to entropy ratio; $n_B/s\simeq 8.7\times 10^{-11}$ gives the reheat temperature,

$$\begin{aligned} T_R\simeq 5 \times 10^9\,\, \text {GeV}, \end{aligned}$$

(6.9)

Fig. 12

Scalar power spectrum (left) and GW power spectrum (right) for BP-3 where the dark fermion is regarded as the RHN

to explain the entire baryon asymmetry of the universe via the leptogenesis mechanism. For BP-3, $m_N \gg T_R$ which ensures that RHN is not produced from the thermal radiation bath so the only contribution to RHN abundance comes from waterfall field decay, therefore the non-thermal nature of leptogenesis.⁹ The impact of the presence of RHN on the inflaton and waterfall field evolution modifying the critical waterfall transition point with respect to the number of e-folds is shown in Fig. 3 (see BP-3). Due to this modification, the predictions in the scalar and GW spectra get affected as shown in Fig. 12. PBHs here may explain some fraction of DM, see Fig. 5 and is testable in the upcoming Nancy Grace Roman Space Telescope.

8 Higgs vacuum stability and quantum fluctuations

In this section, we briefly discuss the fate of the universe due to the presence of the Higgs field during inflation. The electroweak vacuum stability or meta-stability demands that the Hubble rate H during inflation should be less than the Higgs instability scale $\Lambda$ that lies between $10^9\lesssim \Lambda /\text {GeV} \lesssim 10^{12}$ depending upon the precise measurement of the top quark mass [171,172,173]. The Higgs field may get pushed over the barrier due to quantum fluctuations that destabilize the Higgs field during inflation if the typical momentum $k\sim H$ is greater than the potential barrier leading to decay of the electroweak vacuum [174,175,176]. Usually, this problem is avoided by the presence of new physics at the instability scale $\Lambda$ or via Higgs-inflaton couplings [174, 177,178,179]. For our BPs in Table 1, $H \sim 10^{12}$ GeV we just give an example where the Higgs vacuum is still not destabilized however a careful analysis involving renormalization group equation (RGEs) following Ref. [161] could be needed to be studied which we plan to take up in future publication.

The effective mass $m_h$ of the Higgs field during inflation should be larger than H in order to ensure that the SM Higgs does not affect the hybrid inflationary scenario we illustrated. For this purpose, following Ref. [161] one may modify the potential Eq. (2.2), by introducing the SM singlet S with the following interactions terms,

$$\begin{aligned} V(S,h)&\simeq \left( -2\,\lambda _{\phi S}\,\phi ^2 + \lambda _S\left( S^2-\dfrac{v_S^2}{2}\right) \right. \nonumber \\&\quad \left. +2\,\lambda _{\psi S}\left( \psi ^2-\dfrac{M^2}{2}\right) \right. \nonumber \\&\quad \left. +2\,\lambda _{hS}\left( h^2-\dfrac{v_h^2}{2}\right) \right) \left( S^2-\dfrac{v_S^2}{2}\right) . \end{aligned}$$

(6.10)

where it can be assumed that the couplings $\lambda _{hS} \ll \lambda$ and $\lambda _{\phi S}$, $\lambda _{\psi S}, \lambda _{S}$ are negligibly small such that they do not affect the inflaton and waterfall dynamics. Due to the additional S field, $v_S$ can be tuned such that the SM Higgs effective mass can be large during inflation such that $m_h \gg H$ condition is satisfied and the Higgs vacuum remains stable [161]. Therefore, the Higgs fluctuations do not excite iso-curvature perturbations and the Higgs field is frozen during inflation for the Hubble rate around $10^{12}$ GeV for the BPs we considered in Table 1 at the end of inflation.

To explain the baryon asymmetry of the universe via leptogenesis, the heavy neutrino $m_N > rsim 10^9$ GeV [180, 181]. Moreover, to keep the model perturbative involving Yukawas for seesaw, there is an upper bound on the RHN that is $m_N\lesssim 10^{15}$ GeV see [180, 182] which is ensured in our choice in BP-3 as shown in Table 1. Lastly, Yukawa coupling $y_h$ being O(1) may affect the RGEs of the SM Higgs and make the vacuum unstable, which we also ensure in our choice of BP-3 in Table 1. For detailed analysis involving RGE, we refer the reader to Ref. [180].¹⁰

9 Discussion and conclusion

In summary, we presented a two-field inflationary model based on the original hybrid model and studied the generation of both gravitational waves and primordial black holes. The effective scalar potential derived by hybrid models has the advantage that they do not require a high level of fine-tuning of the parameters to describe an amplification in the scalar power spectrum. As the issue of fine-tuning is regarded as a major problematic feature in many proposed inflationary models, studying hybrid models should be a plausible scenario to describe enhancement in the scalar power spectrum. To predict acceptable values for the spectral index $n_s$ we introduced $\alpha$-attractor specific scenario involving a pole in the kinetic term following [45] where we evaluated the prediction for the cosmological constraints. In our proposed model we have achieved acceptable values for both spectral index ($n_s\simeq 0.964$) and tensor-to-scalar ratio ($r\simeq 0.0087$) satisfying PBH as dark matter and detectable GW signal. We summarize the main findings of our analysis below:

• Heavy neutrino-like dark fermions impact radiatively the hybrid inflation potential via quantum loop corrections to the inflaton field involved. Consequently, the predictions for CMB observables are affected due to the variation of the dark fermion mass m$_N$ (see Table 1). The prediction for tensor-to-scalar ratio $r\simeq 0.0087$ and scalar spectral index $n_s\simeq 0.964$, lies within current Planck bounds and testable in future CMB experiments like LiteBIRD, SO and CMB S4-Euclid.

• We estimate the power spectrum across all k-values which provides constraints on the dark fermion mass scale from the measurements of CMB spectral distortions (see Fig. 4). Along the valley, when the inflaton field becomes smaller than the critical value $\varphi _c$, the effective mass square of the waterfall field becomes negative. This gives rise to tachyonic instability and the power spectrum shoots up at small scales. The radiative corrections from the dark fermions along with the inflaton control the amplitude of the plateau and define the number of e-folds along the valley (see middle left of Fig. 4).

• Second-order tensor perturbations propagating as GWs that can be with amplitude $\Omega _\textrm{GW}h^2$ $\sim 10^{-9}$ and peak frequency f $\sim$ 0.1 Hz by LISA and $\Omega _\textrm{GW}h^2 \sim 10^{-11}$ and peak frequency of $\sim$ 10 Hz in ET in this model (see Fig. 6).

• Production of PBH of mass around $10^{-13} M_\odot$ as the sole DM candidate in the universe is proposed. This novel DM candidate is also a signature of the scale of dark fermion physics involving inflationary cosmology (see Fig. 5).

• We estimate fine-tuning in our model and found that it is around five orders of magnitude smaller than single field inflation and one order of magnitude smaller than other hybrid inflationary scenarios studied in the literature (see Table 2).

• For choice of benchmark point 3 (BP-3) in Table 1, one finds that dark fermions can be a possible candidate for being a RHN which is responsible for the generation of SM neutrino mass via the seesaw mechanism and leaves imprints in the power spectrum and GW. Not only it can be a fractional DM of the universe but also be tested in future probes of the Nancy Grace Roman Space Telescope (see Fig. 5). Reheating in such a case may proceed via a waterfall field decaying into RHN and is suitable for leptogenesis via subsequent decays of RHN. We find the reheating temperature $T_R\lesssim 5\times 10^{9}$ GeV that may explain matter–antimatter asymmetry leptogenesis, neutrino mass $m_N\simeq 8.28\times 10^{11}$ GeV and the corresponding PBHs are of $10^{-9}\,M_\odot$.

It has been well studied that the presence of non-Gaussianities affects the abundance of PBH formation quite a significant manner; see Refs. [5, 187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203] and along with PBH clustering, see Refs. [204,205,206,207,208,209,210]). It also affects the predictions of induced GW signals (see, e.g., Refs. [211,212,213,214,215,216,217,218] for a review on this topic). Since in our scenario the curvature perturbation is generated and enhanced during the waterfall transition, this non-Gaussianity may have an impact [41] but it is beyond the scope of the present work and will be taken up in future publication; for generic discussions of $f_\text {NL}$ on SIGW and PBH see [219,220,221]. Recently, the quantum loop correction from the PBH-scale perturbation to the CMB-scale is being actively debated for the single-field inflationary scenario [222,223,224,225,226,227,228,229,230,231]. However, the possible No-Go theorem proposed in Ref. [223] would also be interesting to understand for the loop correction in hybrid inflation. We envisage that GW astronomy with the planned global network of GW detectors can make the dream of testing high-scale and fundamental BSM scenarios like seesaw scale and neutrino physics involving UV-completion and inflationary cosmology a reality.

Appendix A: Gravitational wave signal in pulsar timing array

To test the model presented in this paper in the PTA band, the $n_s$ conflicts with the Planck 2018. Moreover, the power spectrum is also ruled out by COBE/Firas due to the steep enhancement at large scales, as shown in the left panel of Fig. 13 for the BP in Table 3. This conflict can be encountered by including the effect of non-Gaussanity parameter $f_\text {NL}$, we leave this for future work.¹¹

Fig. 13

Left panel shows the power spectrum by solving the exact scalar perturbation equations for the BP given in Table 3. The shaded region corresponds to the constraints from the present (solid) and future (dashed) experiments. For further details on experimental bounds, see the main text in Sect. 3. The right panel indicates the energy density of GWs for Eq. (4.1) for the BPs given in Table 3. The colored shaded regions indicate the sensitivity curves of present (solid boundaries) LIGO O3 [138], NANOGrav [81] and future (dashed boundaries) LIGO O5, SKA [129], THEIA [130], LISA [131], $\mu$-ARES [132], BBO [133], U-DECIGO [134, 135], CE [136] and ET [137] experiments. The hatched region shows the astrophysical background [139]

Table 3 Benchmark points for model parameters in the PTA