Primordial black holes in nonminimal derivative coupling inflation with quartic potential and reheating consideration

We investigate the generation of primordial black holes (PBHs) with the aid of gravitationally increased friction mechanism originated from the nonminimal field derivative coupling (NMDC) to gravity framework, with the quartic potential. Applying the coupling parameter as a two-parted function of inflaton field and fine-tuning of five parameter assortments we can acquire ultra slow-roll phase to slow down the inflaton field due to high friction. This enables us to achieve enough enhancement in the amplitude of curvature perturbations power spectra to generate PBHs with different masses. The reheating stage is considered to obtain criteria for PBHs generation during radiation dominated era. We demonstrate that three cases of asteroid mass PBHs (10-12M⊙\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-12}M_{\odot }$$\end{document}, 10-13M⊙\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-13}M_{\odot }$$\end{document}, and 10-15M⊙)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-15}M_{\odot })$$\end{document} can be very interesting candidates for comprising 100%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$100\%$$\end{document} , 98.3%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$98.3\%$$\end{document} and 99.1%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$99.1\%$$\end{document} of the total dark matter (DM) content of the universe. Moreover, we analyse the production of induced Gravitational Waves (GWs), and illustrate that their spectra of current density parameter (ΩGW0)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\Omega _{{\mathrm{GW}}_\mathrm{0}}})$$\end{document} for all parameter Cases foretold by our model have climaxes which cut the sensitivity curves of GWs detectors, ergo the veracity of our outcomes can be tested in light of these detectors. At last, our numerical results exhibit that the spectra of ΩGW0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega _{{\mathrm{GW}}_\mathrm{0}}}$$\end{document} behave as a power-law function with respect to frequency, ΩGW0(f)∼(f/fc)n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega _{{\mathrm{GW}}_\mathrm{0}}} (f) \sim (f/f_c)^{n} $$\end{document}, in the vicinity of climaxes. Also, in the infrared regime f≪fc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\ll f_{c}$$\end{document}, the power index satisfies the relation n=3-2/ln(fc/f)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=3-2/\ln (f_c/f)$$\end{document}.


Introduction
The concept of primordial black holes (PBHs) generation from the over-dense regions of the early universe should be addressed in the first instance to the recommendation by Zel'dovich and Novikov in 1967 [1], and subsequent progressive researches by Hawking and Carr in 1970s [2][3][4]. Lata e-mail: s.heydari@uok.ac.ir b e-mail: kkarami@uok.ac.ir (corresponding author) terly, the prosperous discovery of gravitational waves (GWs) arises from two coalescing black holes with masses around 30M (M denotes the solar mass) by LIGO-Virgo Collaboration [5][6][7][8][9], has attracted wide attention to PBHs as a possible source of the universe dark matter (DM) content and GWs [10][11][12][13][14][15]. Inasmuch as the PBHs generate from the gravitationally collapse of the primordial density perturbations, hence they are not bounded by Chandrasekhar mass limitation owing to their non-stellar origin, and have wide range of masses. In mass scale around O(10 −5 )M , PBHs locate in the permitted region by OGLE data [16]. Hence they can be considered as the origin of ultrashort-timescale microlensing events in the OGLE data with the fractional abundance about O(10 −2 ) of the total DM content. By reason of the invalidity of gravitational femtolensing of gammaray bursts [17] through the declined lensing effects by the wave effects [18,19], the Subaru Hyper Supreme-Cam (Subaru HSC) microlensing observations impose no restriction on PBHs with masses less than 10 −11 M . Thereupon, considering ineffectiveness of the restriction imposed by white dwarfs [20] on mass range (10 −14 −10 −13 )M due to numerical simulations in [21,22], the existence of PBHs in the mass window of asteroid masses (10 −16 −10 −11 )M as interesting candidates comprising all DM content of universe is possible [23].
Generally, producing a large enough amplitude of primordial curvature perturbations (R) during inflationary epoch is necessary to generate PBHs during radiation dominated (RD) era. Overdense regions can be formed when superhorizon scales associated with the large amplitude of R become subhorizon during RD era, and gravitationally collapse of these overdensities generate PBHs. Notwithstanding the restriction of the power spectrum of R at large scales by CMB anisotropies normalization to 2.1 × 10 −9 [24], PBHs generation requires an enhancement in the power spectrum of R to order O(10 −2 ) at scales smaller than CMB scales. Newly, different techniques for multiplying the amplitude of the power spectrum of R at small scale by seven orders of magnitude in comparison with CMB scales have been inspected in many investigation . One of the proper ways to achieve a rise in the scalar power spectrum is a brief period of ultra slow-roll (USR) inflation due to declining the speed of inflaton field via gravitationally enhanced friction. The framework of nonminimal derivative coupling to gravity (NMDC) beside the fine-tuning of the parameters of the model can give rise to increase friction gravitationally [50][51][52][53][54][55][56][57][58][59].
The nonminimal derivative coupling model is a subclass of a generic scalar-tensor theory with second-order equations of motion namely Horndeski theory [56][57][58][59][60], which prevents the model from negative energy and pertinent instabilities [61]. A characteristic of the nonminimal field derivative coupling to gravity is that the gravitationally increased friction mechanism can be applied for generic steep potentials such as quartic potential [57]. In view of the Planck 2018 TT,TE,EE+lowE+lensing+BK14+BAO data the quartic potential in the standard inflationary model is incompatible with the latest CMB data and completely ruled out [24]. As regards a lot of effort have been put into making the Higgs field consistent with inflation, One way is employing the nonminimal field derivative coupling to gravity [55] (see [62] for original study). With this in mind, we contemplate the quartic potential in the nonminimal derivative coupling setup and try to rectify its observational results at CMB scales and examine the feasibility of PBHs generation of the mentioned mass scales simultaneously.
It is worth noting that an analogous investigation has been done by Fu et al. [51] and Dalianis et al. [53]. In [51] for the chosen parameter the influence of the temporary nonminimal derivative coupling framework far from the peak position (φ = φ c ) is faded away and the standard slow-roll inflation is recovered. In the other word the NMDC framework just governs during USR phase. But in our work the authentication of the nonminimal derivative coupling framework throughout the inflationary epoch is respected because of the form of our suggested coupling function between field derivative and Einstein tensor.
Akin to [53] we choose quartic potential in NMDC framework with the same coupling function, but in our model inflaton is not Higgs and formation of multifarious PBHs in asteroid, earth, and stellar mass ranges thereto their concurrent scalar gravitational waves are predicted. Moreover reheating consideration for our model terminate in obtaining criteria for PBHs generation during RD era.
Generating of induced gravitational waves can be another effect of re-entry the large enough densities of primordial perturbations to the horizon during the RD epoch, concurrent with the PBHs formation [63][64][65][66][67][68][69][70][71][72][73][74]. Heretofore miscellaneous mechanisms of PBHs and detectable GWs generation owing to intensifying the small scale primordial perturbations have been suggested in many inflationary models. The popular one of them is the single-field model with a potential having the flat region pertinent to USR phase in the environs of an inflection-point which leads to a rise in the curvature perturbations and generate PBHs and GWs [75][76][77]. Another models with Gaussian form of scalar power spectrum [71,78,79], or multifield inflationary models such as double [80,81], and hybrid [82,83] inflation have been able to explicate the formation of PBHs and induced GWs via enhancement of scalar power spectrum during the inflationary era. In the following of this work, we evaluate the current fractional density parameter of induced GWs in our NMDC setup and check the rectitude of our numerical conclusion with the latest sensitivity curves of GWs detectors.
The reminder of this paper is organized as follows. In Sect. 2, we study the basis formulae of the nonminimal derivative coupling framework briefly. In Sect. 3, we explain the suitable mechanism of gravitationally increased friction to enhance the amplitude of the scalar power spectrum to around order O(10 −2 ) in our model. To specify when the superhorizon climax scales of perturbations re-enter the horizon, we take into account reheating consideration in Sect. 4. In the succeeding sections we calculate different masses and fractional abundances for PBHs in Sect. 5 and current energy spectra for induced GWs in Sect. 6. Ultimately, we epitomize our outcomes in Sect. 7.

Foundation of nonminimal derivative coupling framework
We initiate this section with the general action of the nonminimal derivative coupling (NMDC) model, in which the inflaton field is nonminimally coupled to the Einstein tensor via its derivative as follows [55][56][57][58][59] wherein g is the determinant of the metric g μν , R is the Ricci scalar, G μν is the Einstein tensor, ξ is the coupling parameter with dimension of (mass) −2 , and V (φ) is the potential of a scalar field φ. As we explained in the previous section the action (1) appertains to the Horndeski theory, which produces second order equations. The Lagrangian of this class of theories consists of the term where G 5 is a general function of φ and kinetic term X = − 1 2 g μν ∂ μ φ∂ ν φ. By taking G 5 = − (φ)/2, then integrating by parts, and defining ξ ≡ d /dφ, the mentioned Horndeski's Lagrangian leads to the NMDC action (1). In [55][56][57][58][59] the Case of constant ξ has been investigated. Moreover, in [84] another investigation has been done by considering ξ = β/φ n , where β is a parameter with dimension of (mass) n−2 and n is an integer, in order to adapt the observational results of steep potentials like Higgs with the present CMB data.
In this work, we consider a more general and appropriate function of φ such as ξ = θ(φ), so that not only the observational results of the quartic potential on large scales is reanimated, but also a tenable outline for PBHs and induced GWs production on small scales can be achieved. For this purpose, we start to study the background evolution of the homogeneous and isotropic universe with the flat Friedmann-Robertson-Walker (FRW) metric and mostly positive signature g μν = diag − 1, a 2 (t), a 2 (t), a 2 (t) , wherein a(t) is the scale factor and t is the cosmic time. We also adhere to convention of reduced Planck mass is identical to unity , over all this paper. In addition, we review the fundamental relations for produced perturbations during inflation epoch in our model described by action (1).
The Friedmann equations and the equation of motion for φ can be obtained via taking variation of action (1) with respect to g μν and φ as follows where H ≡ȧ/a is the Hubble parameter, the dot indicates derivative with respect to the cosmic time t, and moreover the notation (, φ) indicates d/dφ. In the NMDC setup, the slow-roll parameters are specified as [56,57] It is known that, the validity of the slow-roll inflation is Under the slow-roll approximation, the energy density of the inflaton is dominated by the potential energy and the Eqs. (2), (3) and (4) abbreviate to 3HφA wherein Also because of the ease of calculation, we take the below assumption throughout the slow-roll domain By utilizing the assumption (10), the Eqs. (6), (7) and (8) take the following simpler form as 2Ḣ + Aφ 2 0, 3HφA Applying Eqs. (11) and (13), we obtain where Note that ε ε V is concluded for A 1, and one can attain the restoration of the standard slow-roll inflationary formalism. As we can see from (14), A 1 (high friction domain) results in ε ε V , therefore inflaton embogs owing to gravitationally increased friction, and rolls more slowly in comparison to that in the standard slow-roll story line. On the other hand, we expect that the severe decline of ε in the mentioned high friction domain can give rise to an increase in the amplitude of curvature perturbations (R). With this aspect in mind, we should inspect the power spectrum of the curvature perturbations, which is calculated by [58] at the moment of the Hubble radius crossing by the comoving wavenumber (c s k = a H) as follows wherein and It is obvious, for defined ξ as a general function of φ, one can easily find the coincidence of the relations (17)- (22) and their counterparts in [57]. Applying the approximated background equations in the slow-roll domain (11)-(13), the power spectrum (16) can be written as The restricted amplitude of the scalar perturbations power spectrum evaluated by the Planck team from cosmic microwave background (CMB) anisotropy at the CMB pivot scale (k * = 0.05 Mpc −1 ) [24] is Since the scalar spectral index (n s ) can be calculated from the scaler power spectrum through the definition n s − 1 ≡ d ln P s /d ln k, therefore one can easily find the following relation between n s and slow-roll parameters in the NMDC framework [54] in which It is evident that for θ(φ) = 0 the coupling term disappears and the relations of standard slow-roll inflation are recovered. The slow-roll approximation of tensor power spectrum at c t k = a H and then the tensor-to-scalar ratio in NMDC setup have been evaluated in [58] as follows The observational restriction according to the Planck 2018 TT,TE,EE+lowE+lensing+BK14 +BAO data at the 68% CL on the scalar spectral index (n s ), and the upper bound on the tensor-to-scalar ratio (r ) at 95% CL [24] are as follows

Intensification of curvature perturbations power spectrum
It is well known that an increase in the amplitude of curvature perturbations (R) on scales smaller than the CMB scales can lead to detectable generation of PBHs and GWs. It is understood that, such a rise in scalar power spectrum can be achieved with the aid of a span of USR inflation. Suitable choices of parameters and coupling function between the gravity and derivatives of the scalar field in the NMDC setup can establish such an USR phase due to gravitationally enhanced friction [50][51][52][53][54]. With this regard, in this section we define two-parted form of coupling function θ(φ), to have a correspondence with the available data on large scales (CMB) and a peak in the curvature perturbations on smaller scales, as follows in which The first part of our coupling function (32) is a general form of the ones taken in [56,84]. Another part (33) is taken by Fu et al. [51] in order to investigate PBHs and GWs formation in NMDC setup for the potential V (φ) ∝ φ 2/5 , but we consider a combination of these two functions pursuant to [53]. By changing α to (α − 1) in (32), regardless of α coefficient, the coupling function of [53] is retrieved. Concerning (33), one can easily find that θ I I (φ) has a climax at φ = φ c of the height and width that are defined by ω and σ respectively. It is obvious that for the far-off field values of φ c , the function θ I I (φ) evanesces. In our model θ I (φ) is necessary to make the quantum fluctuations on the CMB scales compatible with the latest data. On the other hand, the quantity θ I (φ)θ I I (φ) depending on the suitable choices for the parameters makes the scalar power spectrum of the curvature perturbations on small scales enlarge to get the enough amplitude to generate the detectable PBHs and GWs. Note that α and ω are the dimensionless parameters, while φ c , σ , and M are parameters with the dimensions of mass. Note that the quartic potential is entirely ruled out by Planck 2018 TT,TE,EE+lowE+len-sing+BK14+BAO in the standard inflationary groundwork [24]. Therefore, this persuades us to study the quartic potential in the nonminimal derivative coupling framework with the aim of getting compatible observational result with Planck 2018. Hence, we keep on our investigation to see if the quartic inflationary model can derive the observable inflation and generate detectable PBHs and GWs in NMDC framework simultaneously. The quartic potential is given as follows where λ is the dimensionless constant, and fixed by the scalar power spectrum restriction (24) in our investigation. From our numerical results listed in Tables 1 and 2, one can see for reduced λ up to order O(10 −9 ) in comparison with λ 0.1 appertain to Higgs field in [53], the quartic potential can drive the viable inflation and detectable enhanced PBHs and GWs in the nonminimal derivative coupling groundwork. Now we want to estimate the order of needed parameters to increase the amplitude of the scalar power spectrum (P s ) to order O(10 −2 ), which is enough value to produce PBHs. Hence, we derive an approximated relation between P s at φ = φ c and P s at φ = φ * , noting that φ * is the field value at the time of horizon crossing by k * . So applying Eqs. (23) and (31) where furthermore we suppose that eventually we conclude Clearly, regarding (24) and (38), the climax of the power spectrum at φ = φ c will be of order O(10 −2 ) provided by ω ∼ O(10 7 ). We know that in the vicinity of peak position the term of θ I (φ)θ I I (φ) is dominated in (31), ergo the power spectrum (23) around the peak position abbreviates as the following form where It is obvious from (39) that, the width of P s around the peak position is determined by parameter σ . In view of these estimations and expressed conditions in (37), we set q = 20 and contemplate five different Cases of parameters which are represented in Table 1. It is worth noting that an assortment consists of seven free parameters {α, M, λ, q, ω, φ c , σ } identifies our model, in which the parameters α, M, q, and λ are correlated together via Eq. (36). The calculated values of inflationary observables n s , r , and quantities related to generated PBHs are represented in Table 2.
It is known that the observable inflation takes place in about 50-60 e-folds number from the horizon crossing moment by k * to the end of inflation. So regarding the reheating process duration, we set the horizon crossing e-fold number (N * ) as 53.62 for Case A, 55 for Case B, 54.8 for Case C, 55.26 for Case D, and 53.35 for Case F. From Fig. 2, one can see that the end of inflation epoch is specified by resolving ε = 1 at the e-fold number N end = 0 for each Case. We solve the background Eqs. (2)-(4), using the quartic potential (34) and the NMDC coupling function (31) In this figure, one can see a plane zone lasting about 20 efolds around φ = φ c in all field excursion. In this zone the scalar field revolves very slowly owing to the high friction, and the slow-roll condition is broken, thus one has an Ultra Slow-Roll (USR) stage. Figure 3 shows that, the curvature perturbations power spectrum will be intensified during this USR stage. Figure 2 shows the evolution of the first slow-roll parameter (ε) (left plot), the second slow-roll parameter (δ φ ), and θ ,φ Hφ/A (right plot) with respect to e-fold number (N ) during the observable inflation ( N = N * − N end ) for Cases listed in Table 1. In this figure a severe reduction in the value of ε in the USR stage for each Case up to the order O(10 −10 ) can be seen, which gives rise to a large amplification in the amplitude of the curvature power spectrum. Also one can see δ φ violates the slow-roll condition due to exceed unity in this stage temporarily. As regards the evolution of the slow-roll parameters at horizon crossing e-fold number (N * ), one can see that the previously mentioned slow-roll conditions are respected for all Cases, therefore we can use Eqs. (25) and (28) to evaluate the scalar spectral index (n s ) and tensor-to-scalar ratio (r ) in our model. Regarding the numerical outcomes listed in Table 2, we can realize that, the values of r for all Cases and n s for Cases B and F are consonant with the 68% CL of the Planck 2018 TT,TE,EE+lowE+lensing+BK14+BAO data, and the values of n s for Cases A, C, and D are consistent with 95% CL [24]. Consequently, we have been able to revive the ruled out quartic potential in the standard inflationary scenario by means of selecting the nonminimal derivative coupling framework. It is worth mentioning that, the recent BICEP/Keck collaboration's constraints suggest that the tensor-scalar ratio r < 0.036 [85], which thereunder as we  (35) is acquired under the slow-roll conditions with the proviso (10), applying that to calculate the scalar power spectrum in the USR stage is impermissible owing to infraction of the slow-roll limitation in this stage. Inevitably, to calculate the exact scalar power spectrum, numerical analysis of the following Mukhanov-Sasaki (MS) equation is needed, for all Fourier modes The MS Eq. (40) represents the evolution of curvature perturbations (R) in the Fourier space (υ k ) during inflation epoch, from sub-horizon scale (c s k a H) until exiting the horizon and getting the constant value at super-horizon scales (c s k a H). In (40) the prime indicates a derivative with respect to the conformal time η = a −1 dt, and We consider the following Fourier transform of the Bunch-Davies vacuum state as the initial condition for the subhorizon scale (c s k a H) [59] υ k → e −ic s kη √ 2c s k .
After finding the numerical solutions of the Mukhanov-Sasaki Eq. (40), the exact scalar power-spectrum for each mode υ k is obtained as The numerical results for the exact value of the climax of the scalar power spectrum (P peak s ), and the related comoving wavenumber (k peak ) for each Case of Table 1 are represented  in Table 2. Also, in Fig. 3 the exact power spectra for all Cases of Table 1 as a function of comoving wavenumber (k), and the current observational restrictions are illustrated. The purple, green, red, blue, and brown lines are related to the actual P s for Cases A, B, C, D, and F, respectively. This figure exhibits that, during the slow-roll inflationary domain on large scales around the CMB scale (k ∼ 0.05 Mpc −1 ), the power spectra for all Cases are of order O(10 −9 ) which are consonant with the observational data (24). Furthermore, in the USR phase on smaller scales, one can see the multiplication of the power spectra to order O(10 −2 ), which is sufficient amount to generate PBHs. Into the bargain, we have inferred that the scaling behaviour of P s in the vicinity of peak position can be appraised as power-law function versus wavenumber P s ∼ k n . Thereby for Case D, P s ∼ k 2.58 for k < k peak = 3.94 × 10 13 Mpc −1 , and P s ∼ k −1.31 for k > k peak = 3.94 × 10 13 Mpc −1 have been delineated in Fig. 3.

The reheating stage
To keep on our study on the produced PBHs and GWs from the intensified scalar power spectrum, we should know when the super-horizon perturbation modes come back to the horizon after the inflation epoch to collapse. It is known that to shift from a supercooled universe after the inflation epoch to the hot RD epoch a thermalizing phase, so-called reheating stage, is required (see [89] for a review). In spit of uncertainty of the physics of reheating, one can use the simple canoni-  11) and (13) cal outline [90][91][92], whereby generated cold gases from the oscillations of the inflaton filed around the minimum of the potential transform to a plasma of relativistic particles, after that the standard RD area of the Hot Big Bang theory takes over. The effective equation-of-state parameter of the canonical reheating stage is ω reh = 0, and duration of this epoch is proportional to the inflaton decay-rate and reheating temperature (T reh ). One should attach the importance of the reheating consideration insomuch that prolonged reheating stage with the low temperature can result in return the climax scales of the curvature power spectrum to the horizon in the reheating stage and form the PBHs in this stage instead of RD epoch. Accordingly, in the following we try to find a criterion to show that in our investigation the mentioned scales come back to the horizon in RD epoch. In this regard, we pursue the used technique in [29,30].
Suppose that a scale like k −1 exits the horizon N k e-folds before the end of inflation, right at the moment it returns to the horizon we have wherein a k,re , a end , and w denote the scale factor at the horizon re-entry, the scale factor at the end of inflation, and the equation-of-state parameter, respectively. Then we define the post-inflationary e-folds number related to the scale k −1 from the time of re-entry to the end of inflation asÑ k ,  [24]. The yellow, cyan, and orange shadowy zones represent the restrictions from the PTA observations [86], the effect on the ratio between neutron and proton during the big bang nucleosynthesis (BBN) [87], and the μ-distortion of CMB [88] respectively. The power-law form of P s with respect to k is illustrated by black dashed line for Case D Thus, from Eqs. (44) and (45) we drive the following relatioñ Duration of the reheating stage is defined by the e-folds number from the end of reheating to the end of inflation as N reh ≡ ln (a reh /a end ), wherein a reh denotes the scale factor at the end of reheating. Applying the continuity equation, we can specify the e-folds numberÑ reh in terms of the energy density at the end of reheating stage (ρ reh ) as follows in which the energy density at the end of inflation is evaluated via ρ end = 3H 2 end . By making a comparison betweenÑ k and N reh one can specify when the horizon re-entry takes place for the scale k −1 associated with the PBH formation. In the way that, ForÑ k >Ñ reh , the horizon re-entry takes place during the RD era, and forÑ k <Ñ reh the mentioned scale comes back to the horizon during the reheating stage.
Respecting the canonical reheating (ω reh = 0), duration of the observable inflation is related to the post-inflationary expanse [93,94] through As we mentioned heretofore, N = N * − N end , which N end represents the e-fold number at the end of inflation, moreover N * , ε * , and V * denote e-fold number, the value of the first slow-roll parameter, and the potential value when the pivot scale k * = 0.05Mpc −1 leaves the horizon, respectively. Furthermore, regarding the relation between ρ reh and T reh (reheating temperature) in the form of ρ reh = (π 2 /30)g * T 4 reh , and applying Eq. (47) we can evaluate T reh as where g * designates the effective number of relativistic species upon thermalization, which is determined as g * = 106.75 in the Standard Model of particle physics at high temperature. It is worth noting that, reheating stage has lied between the end of inflation and BBN, which means the constraint 10 −2 GeV T reh 10 16 GeV can be accepted. Now we are able to characterize a criterion for specifying that the PBHs formation occurs during reheating or RD era. To this aim, in Fig. 4 the duration of the observable inflation ( N ) is schemed for all Cases of Table 1 Table 3 for defined Cases in Table 1. Here, k −1 reh specified as k reh = e −Ñ reh /2 k end [29], and k −1 end are the scales correspond to the horizon re-entry at the end of reheating and end of inflation, respectively. Using of these results we can establish the following critical conditions for PBHs formation during RD era In other words, the climax scales of the scalar power spectrum re-enter the horizon and collapse to form PBHs during   Table 3 with their critical counterparts shows that in our study the conditions (51) are respected, and the scale k −1 peak comes back to the horizon and collapses to generate PBH during RD epoch. In Fig. 4 the behavior of the exact scalar power spectra in terms of k for all Cases of Table 1 is schemed. The shadowy zones symbolize the re-entered scale to the horizon during reheating stage, and it is obvious from this figure that, the climax scales comes back to the horizon during RD epoch. Subsequently, applying the corresponding formulations with RD epoch to evaluate mass fraction of the PBHs and energy density of the GWs is permitted in the succeeding sections.

Primordial black holes generation
As we described heretofore, it is known that an enough increase in the amplitude of curvature perturbations on scales smaller than the CMB scales can lead to detectable production of PBHs and GWs. In this section, we will inquire PBHs generation in NMDC model. Pursuant to the discussions in the previous section, the re-entry of super-horizon perturbation modes generated in the inflation epoch, occurs during RD era after inflation. If the power spectrum of these modes have large enough amplitude, the gravity of the related overdense zones can overpower the radiation pressure and collapse to generate PBHs. The mass of generated PBH is associated with the horizon mass at the moment of re-entry through wherein γ is the efficiency of collapse, and in this paper we contemplate γ = ( 1 √ 3 ) 3 [4]. As be mentioned in the previous section g * = 106.75 for PBHs generation in RD epoch. On the presumption that the distribution function of perturbations abide by Gaussian statistics, and applying Press-Schechter theory, the generation rate for PBHs with mass M(k) is evaluated [95,96] as follows in which "erfc" designate the error function complementary, and δ c symbolizes the threshold value of the density perturbations for PBHs generation. We contemplate δ c = 0.4 in the way that is recommended by [97,98]. In addition σ 2 (M) denotes the coarse-grained density contrast with the smoothing scale k, that is specified as wherein P s is the power spectrum of curvature perturbations, and W denotes the window function which is specified as Gaussian window W (x) = exp −x 2 /2 in this paper. In the light of characterizing the abundance of PBHs, the present fractional density parameter of PBHs ( PBH ) to Dark Matter ( DM ) is evaluated as where DM h 2 0.12 is given by Planck 2018 data [24] for current density parameter of Dark Matter. At last we are able to calculate the PBHs abundance for all Cases of Table 1, using the exact power spectrum attained via resolving the Mukhanov-Sasaki equation and Eqs. (52)- (55). The obtained numerical and schemed results are represented in Table 2 and Fig. 5. It is revealed from our conclusion that, for Case A the climax of anticipated PBHs abundance is located at 28.23M with f peak PBH 0.0010, this sort of stellarmass PBHs can gratify the restriction from the upper limit on the LIGO merger rate, and they are suitable candidate to expound the GWs and the LIGO events. For Case B the PBHs mass spectrum has a peak at M peak PBH = 4.97× 10 −6 M and f peak PBH = 0.0434 which is placed in the permitted area by the ultrashort-timescale microlensing events in OGLE data, thus such a Case of PBHs with earth-mass can be useful to explain microlensing events. In Cases C, D, and F our model has foretold three PBHs abundance in asteroid-mass scales, It is worth noting that, recently Picker in [111] has demonstrated that cosmological PBHs foretold by Thakurta metric solution in asteroid mass range of (10 −17 − 10 −12 )M can be hotter than common Schwarzschild case and having completely evaporated by today. On the other hand, it has been shown in [112,113] that, the cosmological black holes and their mass growth cannot be described by Thakurta metric in the early universe. The authors of [112] showed that radiation, PBHs or particle dark matter cannot generate the requisite energy flux for providing the mass growth of Thakurta black holes. Moreover they verified that this solution is not applicable for binary black hole. The authors of [113] proved that the depicted space-time by Thakurta metric has not black hole event or entrapping horizon, ergo the null energy proviso as a solution of the Einstein equation is violated. Nevertheless, this issue is currently under dispute in literatures, as shown quite recently by Kobakhidze and Picker in [114] that, the Thakurta metric can indeed describe a cosmological black hole with apparent horizon.  [99]. The border of the red shaded domain depicts the upper bound on the PBH abundance ensued from the LIGO-VIRGO event consolidation rate [100][101][102][103]. The brown shadowy domain portrays the authorized region for PBH abundance owing to the ultrashort-timescale microlensing events in the OGLE data [16]. The green shaded area allots to constraints of microlensing events from cooperation between MACHO [104], EROS [105], Kepler [106], Icarus [107], OGLE [16], and Subaru-HSC [108]. The pink shadowy region delineates the constraints related to PBHs evaporation such as extragalactic γ -ray background [109], galactic center 511 keV γ -ray line (INTEGRAL) [23], and effects on CMB spectrum [110] 6 Induced gravitational waves in NMDC model In this section, we want to study the production of induced GWs in our model. This occurs due to re-entry the large enough densities of primordial perturbations to the horizon concurrent with the PBHs formation during the RD epoch. Considering the second order effects in perturbation theory, it is proven that the dynamics of tensor perturbations is originated from first order scalar perturbations. Choosing the conformal Newtonian gauge, the perturbed FRW metric takes the form [115] where η, , and h i j indicate the conformal time, the firstorder scalar perturbations, and the perturbation of the secondorder transverse-traceless tensor in turn by turn. As we explained in Sect. 4, during the reheating stage after the inflationary phase the inflaton field almost completely decay into light particles to thermalize the universe to initiate the RD era of HBB. In this regard, the inflaton field has a trivial influence on the cosmic evolution during RD era, and can be ignored. Accordingly, one can easily apply the standard Einstein equation to inquire the generation of the induced GWs in RD epoch. Subsequently, the second-order tensor perturbations h i j gratify the following equation of motion [115,116] where H ≡ a /a indicates the conformal Hubble parameter, and T lm i j is the transverse-traceless projection operator. The GW source term S i j is considered as The scalar metric perturbation in the Fourier space during the RD is given by [116] where k is the comoving wavenumber, and the primordial perturbation ψ k is associated with the power spectrum of the curvature perturbation through the two-point correlation function as follows At last, the energy density of induced GWs during the RD epoch can be obtained as [63] GW (η c , k) wherein indicates the Heaviside theta function, and η c denotes the time when the increasing of GW is halted. The current energy spectra of the induced GWs is associated with the energy spectra at η c through the following equation [86] GW 0 h 2 = 0.83 g * 10.75 wherein r 0 h 2 4.2 × 10 −5 is the present radiation density parameter, and g * 106.75 indicates the effective degrees of freedom in the energy density at η c . The association between present frequency and wavenumber can be evaluated as Now we are in a position to attain the current energy spectra of the scalar induced GWs related to PBHs, applying the actual scalar power spectrum owing to solving the MS equation, and Eqs. (61)-(63) for all Cases of Table 1. The schemed outcomes are represented in Fig. 6. Furthermore, to inquire the validity of anticipated outcomes of our model, the sensitiveness curves of GWs detectors comprising European PTA (EPTA) [117][118][119][120], the square kilometer array (SKA) [121], Advanced Laser Interferometer gravitational wave observatory (aLIGO) [122,123], laser interferometer space antenna (LISA) [124,125], TaiJi [126], and TianQin [127], are designated in the figure. Additionally 95% CL from the NANOGrav 12.5 yr experiment is delineated by the black shadowy area [128][129][130] Ultimately, after inspecting the scaling of GW 0 , we perceive that in the proximity of climaxes, the current density parameter spectra of induced GWs behave as a power-law function with respect to frequency ( GW 0 ( f ) ∼ f n ) [52]. In Fig. 6, this parametrization is plotted fore Case D with dashed black line. For this Case, we estimate the frequency of climax as f c = 0.074Hz and evaluate GW 0 ∼ f 1.9 for f < f c , and GW 0 ∼ f −2.51 for f > f c . Besides, for the infrared domain f f c , we obtain a log-contingent power index as n = 3 − 2/ ln( f c / f ) which is in constancy with the analytical consequences acquired in [131,132].

Conclusions
Here, we investigated the feasibility of PBHs generation in a single-field inflationary model, based on the framework of nonminimal field derivative coupling to the Einstein tensor pertaining to the Horndenski theory recounted by action (1). The nonminimal derivative coupling to gravity leads to a brief period of ultra slow-roll inflation through gravitationally enhanced friction for the scalar field. This attribute inspired us to revise a steep potential in this framework to comprehend a viable inflation. Thereupon, we contemplated quartic potential (34), which is completely ruled out by Planck 2018 [24] in the standard inflationary model, and tried to recover its observational results in the NMDC framework. By assigning the coupling parameter as a two-parted function of inflaton field (31)- (33), and fine-tuning of the five parameter Cases (A, B, C, D, and F) of the model (see Table 1), we could acquire an epoch of ultra slow-roll inflation on scales smaller than CMB scale making the inflaton slow down, sufficient to generate PBHs. In addition, compatibility of the observational prognostications of the model with Planck 2018 data on CMB scales is another result of these selections.
Utilizing the exact solution of the background equations, we illustrated the behavior of inflaton field (φ) and slowroll parameters (ε and δ φ ) in Figs. 1 and 2 with respect to e-fold number. It is obvious from Fig. 2 that, during the USR phase the slow-roll conditions are respected by the first parameter (ε 1), but infracted by the second one ( δ φ 1). Accordingly, we obtained the numerical solution of the Mukhanov-Sasaki equation to calculate the exact power spectra of R for all Cases of Table 1. It is clear from the exhibited results in Table 2 and Fig. 3 that, the obtained actual power spectra are compatible with the Planck 2018 data on large scales, and have climaxes of sufficient height to generate PBHs on small scales. Furthermore, according to the observational prognostications of our model for inflation era, the values of r for all Cases and n s for Cases B and F are consistent with the 68% CL of the Planck 2018 TT,TE,EE+lowE+lensing+BK14+BAO data, and the values of n s for Case A, Case C, and Case D gratify the 95% CL [24]. Consequently, we have been able to revive the ruled out quartic potential in the standard inflationary scenario by means of the nonminimal derivative coupling framework in our investigation.
In addition, we proved that the superhorizon climax scales of the power spectra related to PBHs formation reenter the horizon after the reheating stage during RD era (see Fig. 4). Therefore, applying the actual P s and Press-Schechter formalism we could attain detectable PBHs abundance with masses of order O (10) Our conclusions indicate that PBHs of Case A is suitable to describe GWs and LIGO events, Case B can be useful to expound microlensing events in OGLE data, and PBHs of Cases C, D, and F can be interesting candidates for composing around 98.32%, 99.11%, and 100% of DM content of the universe (see Table 2 and Fig. 5).
In the end, we inquired generation of the induced GWs subsequent to PBHs formation for all Cases of our model. Our calculation of current density parameter spectra ( GW 0 ) indicates that, all Cases have climaxes at contrasting frequencies with nearly identical heights of order 10 −8 . The climaxes of GW 0 for Cases A and B have placed at frequencies 10 −10 Hz and 10 −7 Hz, respectively, and both Cases can be traced via the SKA detector. Moreover, the spectra of GW 0 for Cases C, F, and D have climaxes localized at mHz and cHz bands which are tracked down by LISA, TaiJi, and TianQin (see Fig. 6). Hence, validity of our model can be assessed in view of GWs via the extricated data of these detectors. Also, we demonstrated that in the vicinity of climaxes, the spectra of GW 0 behave as a power-law function with respect to frequency ( GW 0 ( f ) ∼ f n ). We evaluated GW0 ∼ f 1.9 for f < f c and GW 0 ∼ f −2. 51 for f > f c . For the infrared domain f f c the power index has been ascertained as n = 3 − 2/ ln( f c / f ), which is in constancy with the analytical consequences acquired in [131,132].
Acknowledgements The authors thank the referee for his/her valuable comments.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The provided figures, and also the cited articles, include all relevant data.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .