title: PBHs and secondary GWs from ultra slow roll and punctuated inflation
authors: H. V. Ragavendra Pankaj Saha L. Sriramkumar Joseph Silk
text: 
  
  
  
  E-mail: ragavendra@physics.iitm.ac.in 
  ^1Department of Physics, Indian Institute of Technology Madras, 
  Chennai 600036, India
  
  E-mail: pankaj@physics.iitm.ac.in
  ^1Department of Physics, Indian Institute of Technology Madras, 
  Chennai 600036, India
  
  E-mail: sriram@physics.iitm.ac.in
  ^1Department of Physics, Indian Institute of Technology Madras, 
  Chennai 600036, India
  
  E-mail: silk@iap.fr
  ^2Institut d'Astrophysique de Paris, UMR 7095, CNRS/UPMC Université 
  Paris 6, Sorbonne Universités, 98 bis boulevard Arago, F-75014 Paris, France
  ^3Institut Lagrange de Paris, Sorbonne Universités, 
  98 bis Boulevard Arago, 75014 Paris, France
  ^4Department of Physics and Astronomy, The Johns Hopkins University, 
  3400 N. Charles Street, Baltimore, MD 21218, U.S.A.
  ^5Beecroft Institute for Cosmology and Particle Astrophysics,
  University of Oxford, Keble Road, Oxford OX1 3RH, U.K.
  
  
  
  The primordial scalar power spectrum is well constrained by the cosmological 
  data on large scales, primarily from the observations of the anisotropies in 
  the cosmic microwave background (CMB).
  Over the last few years, it has been recognized that a sharp rise in power on 
  small scales will lead to enhanced formation of primordial black holes (PBHs) 
  and also generate secondary gravitational waves (GWs) of higher and, possibly,
  detectable amplitudes.
  It is well understood that scalar power spectra with COBE normalized amplitude 
  on the CMB scales and enhanced amplitudes on smaller scales can be generated due 
  to deviations from slow roll in single, canonical scalar field models of 
  inflation.
  In fact, an epoch of so-called ultra slow roll inflation can lead to the
  desired amplification.
  We find that scenarios that lead to ultra slow roll can be broadly classified
  into two types, one wherein there is a brief departure from inflation (a scenario
  referred to as punctuated inflation) and another wherein such a departure does 
  not arise.
  In this work, we consider a set of single field inflationary models involving
  the canonical scalar field that lead to ultra slow roll and punctuated inflation
  and examine the formation of PBHs as well as the generation of secondary GWs in
  these models.
  Apart from considering specific models, we reconstruct potentials from 
  certain functional choices of the first slow roll parameter leading to
  ultra slow roll and punctuated inflation and investigate their observational
  signatures.
  In addition to the secondary tensor power spectrum, we calculate the secondary 
  tensor bispectrum in the equilateral limit in these scenarios.
  Moreover, we calculate the inflationary scalar bispectrum that arises in all the 
  cases and discuss the imprints of the scalar non-Gaussianities on the extent of 
  PBHs formed and the amplitude of the secondary GWs generated. 
  We conclude with a discussion on the wider implications of our results.
  
  
  PBHs and secondary GWs from ultra slow roll and punctuated inflation
      Joseph Silk^2,3,4,5
      December 30, 2023
  ====================================================================
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  With the recent observations of gravitational waves (GWs) from merging 
  binary black holes involving a few to tens of solar 
  masses <cit.>, there has been a considerable
  interest in examining whether such black holes could have a primordial 
  origin <cit.>.
  The most popular mechanism to generate primordial black holes (PBHs) is the 
  inflationary scenario (for earlier discussions, see, for example,
  refs. <cit.>; also see the recent 
  reviews <cit.>).
  PBHs are formed when the curvature perturbations generated during inflation
  reenter the Hubble radius during the radiation and matter dominated epochs.
  However, most inflationary models permit only slow roll inflation and, in 
  such cases, the extent of PBHs produced proves to be considerably smaller
  than required for any astrophysical implications (see, for example,
  ref. <cit.>).
  Recall that, on large scales, the primordial scalar power spectrum is strongly
  constrained by the increasingly precise observations of the anisotropies in 
  the cosmic microwave background (CMB) (for recent constraints from Planck,
  see refs. <cit.>).
  In order to lead to a significant amount of PBHs, the scalar power spectrum 
  on small scales should be considerably enhanced from the COBE normalized 
  values over the CMB scales (for an early discussion in this context, see, 
  for instance, ref. <cit.>).
  In inflation, this is possible only when there are strong departures from
  slow roll.
  It boils down to identifying inflationary potentials that permit slow roll
  initially and then violating it for a certain period of time, before restoring 
  it again until close to the termination of inflation.
  
  
  
  In models of inflation driven by a single, canonical scalar field, the 
  so-called ultra slow scenario has turned out to be the most popular 
  mechanism in the literature to enhance scalar power on small scales.
  This scenario involves a period during inflation wherein the first slow 
  roll parameter turns very small (for the initial discussions, see 
  refs. <cit.>; 
  in this context, also see, for instance,
  refs. <cit.>).
  In fact, one finds that the scenario can be further divided into two types, 
  those which admit a brief period of departure from inflation and another 
  wherein no such departure arises.
  The scenario wherein inflation is interrupted briefly is referred to as 
  punctuated inflation (for the original discussions, see 
  refs. <cit.>;
  for later and recent efforts, see 
  refs. <cit.>; for a 
  discussion in the context of PBHs, see refs. <cit.>).
  Interestingly, in such scenarios, the interruption of inflation is inevitably 
  followed by an epoch of ultra slow roll which aids in boosting the power on 
  small scales.
  While, in the case of punctuated inflation, all the slow roll parameters 
  (including the first) turn large briefly, in ultra slow roll inflation, the 
  first slow parameter remains small until the very end of inflation and slow 
  roll is said to be violated due to the large values achieved by the second and 
  higher slow roll parameters.
  
  
  
  Often, the above-mentioned scenarios are achieved with the aid of potentials 
  which contain a point of inflection <cit.>.
  The inflection point seems to play a crucial role in these scenarios in inducing
  a period of ultra slow roll after the short epoch of deviation from slow roll.
  The two stages of slow roll and ultra slow roll lead to either a step or a 
  bump-like feature in the resulting inflationary scalar power spectrum, depending 
  on the details of the intermediate departure from slow roll.
  The lower level of the step is associated with the large scale modes that leave 
  the Hubble radius during the first epoch of slow roll and the power is enhanced 
  on small scales corresponding to modes that leave the Hubble radius during the 
  later epoch of ultra slow roll.
  We should mention here that the punctuated inflationary scenario has been 
  considered to explain the lower power observed at the small multipoles in the 
  CMB data.
  If one chooses the drop in power to occur at scales roughly corresponding 
  to the Hubble radius today, one finds that the resulting power spectrum
  can improve the fit to the CMB data to a certain extent (for an earlier
  analysis, see ref. <cit.>; for a recent discussion, see 
  ref. <cit.>).
  
  
  
  We mentioned above that both ultra slow roll inflation and punctuated inflation 
  can lead to a sharp rise in power on small scales. 
  Evidently, if one chooses the rise to occur at suitable scales, one can utilize 
  these power spectra to lead to enhanced formation of PBHs.
  As has been established, such an enhanced amplitude for the scalar power spectrum
  can induce secondary GWs when these modes reenter the Hubble radius at later 
  times during the radiation dominated epoch (for the original discussions, see,
  for example, refs. <cit.>;
  for recent discussions 
  in this context, see refs. <cit.>).
  These secondary GWs with boosted amplitudes can, in principle, be detected 
  by current and forthcoming observatories such as 
  LIGO/Virgo <cit.>, Pulsar Timing 
  Arrays (PTA) <cit.>,
  the Laser Interferometer Space Antenna 
  (LISA) <cit.>,
  the Big Bang Observer (BBO) <cit.>,
  the Deci-hertz Interferometer Gravitational wave 
  Observatory (DECIGO) <cit.> and
  the Einstein Telescope (ET) <cit.>.
  Moreover, the deviations from slow roll inflation, even as they boost the scalar 
  power spectrum on small scales, also lead to larger levels of scalar 
  non-Gaussianities on these scales (in this context, see, for example,
  refs. <cit.>).
  These non-Gaussianities can, in principle, further increase the extent of PBH
  formation (for early discussions, see, for example,
  refs. <cit.>; for recent 
  discussions, see refs. <cit.>)
  as well as the strength of the secondary GWs 
  (see refs. <cit.>; for a
  very recent discussion, also see ref. <cit.>).
  In this work, we examine the enhanced formation of PBHs and the generation of 
  secondary GWs in ultra slow roll and punctuated inflation.
  We also numerically evaluate the inflationary scalar bispectrum generated on 
  small scales in these scenarios and utilize the results to discuss the 
  corresponding imprints on the extent of PBHs formed and the amplitude of 
  secondary GWs.
  In addition to considering specific potentials that lead to the scenarios of 
  our interest, we choose functional forms for the first slow roll parameter 
  leading to ultra slow roll and punctuated inflation, reverse engineer potentials
  and examine the observational implications (for other efforts in these 
  directions, see, for instance, Refs. <cit.>).
  Interestingly, such an exercise also confirms the understanding that, in models 
  of inflation involving a single, canonical scalar field, a point of inflection in 
  the potential seems essential to lead to ultra slow roll or punctuated inflation.
  
  
  
  This paper is organized as follows.
  In the following section, we shall introduce the different models of our interest
  which lead to ultra slow roll and punctuated inflation.
  In section <ref>, we shall discuss the power spectra that arise in these
  models and illustrate how the intrinsic entropy perturbation associated with the 
  scalar field proves to be responsible for enhancing the amplitude of the curvature
  perturbations.
  In this section, we shall also highlight some of the challenges that one encounters
  in constructing viable models of ultra slow roll and punctuated inflation.
  In section <ref>, we shall consider specific forms for the first slow roll
  parameter leading to ultra slow roll and punctuated inflation, and reverse 
  engineer the potentials that lead to such scenarios.
  We shall also discuss the power spectra that arise in these cases.
  In sections <ref> and <ref>, we shall discuss extent of PBHs 
  formed and calculate the dimensionless parameters characterizing the power as 
  well as bispectra of secondary GWs generated in the models and scenarios of 
  interest.
  We shall also compare our results with the constraints from observations.
  In section <ref>, we shall calculate the dimensionless 
  non-Gaussianity parameter  associated with the scalar bispectrum in 
  all the different cases.
  We shall highlight some of the properties of the non-Gaussianity parameter  
  and then go on to discuss the imprints of the scalar non-Gaussianities on the
  formation of PBHs and the generation of secondary GWs.
  In section <ref>, we shall conclude with a summary of the main results.  
  We shall relegate some of the related discussions to six appendices.
  
  
  
  A few remarks on our conventions and notations are in order at this stage 
  of our discussion. 
  We shall work with natural units such that ħ=c=1 and set the reduced 
  Planck mass to be =ł(8 π  G)̊^-1/2.
  We shall adopt the signature of the metric to be (-,+,+,+).
  Note that Latin indices shall represent the spatial coordinates, except 
  for k which shall be reserved for denoting the wave number. 
  We shall assume the background to be the spatially flat
  Friedmann-Lemaître-Robertson-Walker (FLRW) 
  line element described by the scale factor a and the Hubble parameter H.
  Also, an overdot and an overprime shall denote differentiation with 
  respect to the cosmic time t and the conformal time η, respectively.
  Moreover, N shall denote the number of e-folds.
  
  
  
  
  
  § MODELS OF ULTRA SLOW ROLL AND PUNCTUATED INFLATION
  
  
  In this section, we shall briefly describe the specific models of interest 
  that lead to ultra slow roll and punctuated inflation.
  We should mention that all the five models that we shall discuss in the following
  two subsections contain a point of inflection.
  Recall that, the first slow roll parameter is defined as ϵ_1=-Ḣ/H^2. 
  The higher order slow roll parameters are defined in terms of the first slow
  roll parameter ϵ_1 through the relations 
  
      ϵ_n+1=ln ϵ_nṆ
  
  for n ≥ 1.
  As it is the first three slow roll parameters, viz. ϵ_1,  ϵ_2, 
  and ϵ_3, that determine the amplitude and shape of the power spectrum
  as well as the bispectrum, we shall illustrate the behavior of these slow roll 
  parameters in the models of interest.
  
  
  
  
  
   §.§ Potentials leading to ultra slow roll inflation
  
  
  We shall consider two specific models that permit ultra slow roll inflation.
  The first potential we shall consider which leads to a period of ultra slow 
  roll inflation is often written in the following form (see, for instance,
  ref. <cit.>):
  
      V(ϕ) = V_0 6 x^2 - 4 α x^3 + 3 x^4(1 + β x^2)^2,
  
  where x = ϕ/v, with v being a constant rescaling factor.
  We shall work with the following choices of the parameters involved: V_0/^4 
  = 4×10^-10, v/ = √(0.108), α = 1 and β = 1.4349.
  For these choices of parameters, the inflection point, say, ϕ_0, is located 
  at 0.39.
  We find that, if we choose the initial value of the field to be 
  ϕ_i=3.614, then inflation lasts for about 63 e-folds 
  in the model.
  For convenience, we shall hereafter refer to the potential (<ref>), 
  along with the above-mentioned set of parameters, as USR1.
  
  
  
  The second potential that we shall consider is given by <cit.>
  
      V(ϕ) = V_0 ł{tanhł(ϕ√(6) )̊ 
      + A sinł[tanhł[ϕ/ł(√(6) )̊]̊f_ϕ]̊}̊^2,
  
  and we shall work with the following values of the parameters involved: V_0/^4 
  = 2×10^-10, A = 0.130383 and f_ϕ = 0.129576.
  We find that, for these values of the parameters, the inflection point occurs at 
  ϕ_0 = 1.05.
  For the initial value of the field ϕ_i=6.1, we obtain about
  66 e-folds of inflation in the model.
  We shall refer to the potential (<ref>) and the above set of parameters 
  as USR2.
  
  
  
  As we mentioned, the background dynamics driven by these potentials can be 
  well captured by the behavior of the first three slow roll parameters
  ϵ_1, ϵ_2 and ϵ_3.
  We have plotted the evolution of these quantities as a function of e-folds N
  in figure <ref>. 
  
  
  
  It is clear from the behavior of ϵ_1 that these models permit two 
  different regimes of slow roll, separated by a short phase of departure 
  from slow roll. 
  Note that the value of ϵ_1 during the second regime of slow roll is
  a few orders of magnitude smaller than its value during the initial regime,
  thereby leading to the nomenclature of ultra slow roll inflation.
  We should point out that there is no deviation from inflation in these models, 
  as the first slow roll parameter always remains smaller than unity until the 
  very end of inflation. 
  The transition from slow roll to ultra slow roll is rather rapid and this 
  aspect is reflected by the sharp rise and fall in the amplitude of the 
  second and third slow roll parameters within a short period. 
  It should also be highlighted that the second slow roll parameter ϵ_2 
  is large and negative (about -6 and -7 in USR1 and USR2) during the ultra 
  slow phase when the first slow roll parameter ϵ_1 is rapidly decreasing. 
  The parameter ϵ_2 changes sign when ϵ_1 begins to rise 
  as the field crosses the point of inflection and rolls down towards the 
  minimum of the potential.
  But, ϵ_2 continues to remain relatively large (it is about 0.2 
  and 0.9 in the cases of USR1 and USR2) even during this latter phase, 
  when compared to the typical slow roll values encountered, say, at early 
  times before the transition to the epoch of ultra slow roll.
  
  
  
  To gain a better understanding of the dynamics involved, in figure <ref>, 
  we have also plotted the evolution of the scalar field in phase space for the 
  case of USR2.
  
  
  
  Evidently, trajectories from different initial conditions eventually merge 
  with the primary trajectory of interest.
  The transition to the ultra slow roll regime corresponds to the sharp upward 
  turn in the phase space trajectory when the velocity of the field decreases
  as it nears the point of inflection.
  It is interesting to note that the solution obtained in the slow roll
  approximation closely follows the primary trajectory even during the ultra 
  slow roll regime. 
  The field crosses the point of inflection, eventually emerging from the ultra 
  slow regime, and inflation ends as the field approaches the minimum of the
  potential.
  
  
  
  
  
   §.§ Potentials permitting punctuated inflation
  
  
  As we have discussed, punctuated inflation corresponds to a scenario wherein 
  a short period of departure from inflation is sandwiched between two epochs
  of slow roll.
  With the help of specific examples, we shall illustrate that the period of 
  departure from inflation is inevitably followed by an epoch of ultra slow roll 
  inflation.
  
  
  
  A simple model that has been examined in the early literature which permits 
  interrupted inflation is described by the potential (see ref. <cit.>;
  also see refs. <cit.>)
  
      V(ϕ) = V_0 ł(1 + B ϕ^4)̊.
  
  It should be evident that the inflection point for this model is located
  at ϕ = 0.
  For B/^4 = 0.5520, one finds that the model leads to two epochs of 
  inflation separated by a brief interruption of inflation.
  In fact, around the interruption, the first slow roll parameter rises above 
  unity and quickly falls to very small values, resulting in a period of 
  ultra slow roll.
  It is easy to argue that such a behavior arises due to the constant 
  term V_0 in the potential <cit.>.
  But, the presence of the constant term simultaneously leads to an 
  important drawback of the model.
  Once inflation is restored after the interruption, it is found that 
  the eventual slow roll regime lasts forever.
  There is no conventional termination of inflation as the constant 
  term V_0 sustains slow roll evolution even when the field has 
  reached the bottom of the potential.
  So, one is either forced to terminate inflation by hand or invoke an 
  additional source to end inflation.
  Despite these drawbacks, we shall nevertheless briefly discuss the model 
  due to its simplicity.
  We shall work with the above-mentioned value for the parameter B and 
  choose V_0/^4 = 8 × 10^-13.
  We shall set the initial value of the field to be ϕ_i= 17, and we shall assume that inflation ends after 70 e-folds.
  We shall hereafter refer to this model as PI1.
  
  
  
  The second potential that we shall consider can be expressed as (see, for 
  instance, refs. <cit.>)
  
      V(ϕ) = m^2/2 ϕ^2
      - ł(√(2 λ (n-1)) mn)̊ ϕ^n 
      + λ4 ϕ^2 (n-1),
  
  where n is an integer.
  These potentials contain a point of inflection at  
  
      ϕ_0=ł[2 m^2λ (n-1)]̊^1/[2 (n-2)].
  
  We shall focus on the case n=3, wherein the potential above reduces to
  
      V(ϕ) = m^22 ϕ^2 - 2 m^23 ϕ_0 ϕ^3 
      + m^24 ϕ_0^2 ϕ^4,
  
  and we shall work with the following values of the parameters: m/= 1.8 
  × 10^-6 and ϕ_0/ = 1.9777.
  As we shall soon discuss, these choice of parameters indeed admit 
  punctuated inflation.
  However, one finds, as in the case of PI1, the above potential (for 
  the parameters mentioned) does not naturally result in an end of 
  inflation.
  Despite this limitation, we shall discuss the model, since, it should
  be clear that, modulo the denominator, the potential describing 
  USR1 [cf. eq. (<ref>)] is essentially the same as 
  the potential (<ref>).
  We shall choose the initial value of the field to be ϕ_i
  = 20, and we shall again assume that inflation ends after 70 e-folds.
  We shall refer to this model as PI2.
  
  
  
  Another model we shall consider that permits punctuated inflation is
  motivated by supergravity.
  It is described by the potential (see ref. <cit.>; 
  for a very recent discussion, also see ref. <cit.>)
  
      V(ϕ) = V_0 ł[c_0 + c_1 tanh ł(ϕ√(6 α))̊ 
      + c_2 tanh^2ł(ϕ√(6 α))̊
      + c_3 tanh^3ł(ϕ√(6 α))̊]̊^2,
  
  and we shall work with the following values for the parameters involved: 
  V_0/^4=2.1 × 10^-10, c_0=0.16401, c_1=0.3, c_2=-1.426,
  c_3=2.20313 and α=1.
  This model too contains a point of inflection and, for the above values for
  the parameters, the inflection point is located at ϕ_0 = 0.53.
  If we choose the initial value of the field to be ϕ_i = 7.4, 
  we find that inflation ends after about 68 e-folds.
  We shall refer to this model as PI3. 
  For the above choice of the parameters, apart from a plateau for large field 
  values, the potential admits a second plateau at smaller values of the field.
  As we shall see soon, it is these aspects of the potential that permits 
  punctuated inflation and thereby aids in boosting the scalar power spectrum 
  at small scales.
  
  
  
  As in the case of the ultra slow roll models we had discussed in the previous
  sub-section, we have plotted the first three slow roll parameters ϵ_1,
  ϵ_2 and ϵ_3 for the models PI1, PI2 and PI3 in 
  figure <ref>.
  It is easy to see from the plots that the behavior of the three slow roll 
  parameters are very similar across the models and they differ only in their 
  location of the departures from slow roll.
  Evidently, after an initial slow roll regime, a brief departure from inflation 
  occurs with ϵ_1 growing above unity.
  The interruption of inflation is immediately followed by a period of ultra slow
  roll with ϵ_1 falling to a value that is considerably smaller than its
  value during the initial slow roll regime. 
  Moreover, other than PI3, the models have no definite end of inflation since
  ϵ_1 does not rise to unity once the ultra slow roll regime has 
  begun. 
  Further, note that, when the epoch of ultra slow roll sets in, as in USR1 and 
  USR2, the second slow roll parameter ϵ_2 turns large and negative in 
  all the cases of PI1, PI2 and PI3.   
  The parameter ϵ_2 eventually approaches zero in the cases of PI1 and 
  PI2, since the first slow roll parameter never rises from its very low values
  in these models. 
  However, in PI3, since ϵ_1 rises ultimately leading to the end of
  inflation, the second slow roll parameter ϵ_2 eventually turns 
  positive (from nearly -7) and attains a large value (around 1.2), in 
  very much the same manner it had in USR2.
  As with USR2, we have plotted the behavior of the field in phase space for the
  case of PI3 in figure <ref>.
  It should be clear from the figure that the velocity of the field reaches 
  larger values in the case of PI3 than in the case of USR2 prior to entering
  the ultra slow roll regime.
  Evidently, it is this behavior that is responsible for the brief interruption
  of inflation. 
  
  
  
  
  
  § EVOLUTION OF THE CURVATURE PERTURBATION AND POWER SPECTRA
  
  
  In this section, we shall discuss the scalar and tensor power spectra that 
  arise in the models permitting ultra slow roll and punctuated inflation
  we had introduced in the previous section.
  However, before we go on to discuss the power spectra, we shall illustrate
  the behavior of the curvature perturbations during the period of deviation 
  from slow roll.
  Specifically, we shall highlight the role played by the intrinsic entropy
  perturbations in the enhancement of the amplitude of the curvature 
  perturbations over wave numbers that leave the Hubble radius either immediately
  prior to or during the departure from slow roll.
  
  
  
  
  
   §.§ Scalar and tensor modes, and power spectra
  
  
  Let  and γ_ij denote the curvature and the tensor perturbations 
  at the first order, respectively.
  Also, let _ and γ_ij^ denote the Fourier modes associated
  with these perturbations. 
  Recall that the modes _ and γ_ij^ satisfy the differential
  equations
  
  
      _”+2 z'z _' + k^2 _   =    0,
      γ_ij^”+2 a'a γ_ij^' + k^2 γ_ij^   =    0,
  
  
  where z=√(2 ϵ_1)  a, with ϵ_1 being the first 
  slow roll parameter.
  Moreover, note that, if _ and γ̂_ij^ denote the 
  operators associated with the scalar and tensor modes on quantization, the 
  scalar and tensor power spectra (k) and (k) are defined in terms of
  these operators through the relations
  
  
      ⟨_() _'()⟩   =   2 π^2k^3  P__ S(k) δ^(3)ł(+')̊,
      ⟨ γ̂_ij^() γ̂^ij_'() ⟩   =   2 π^2k^3 𝒫__ T(k) δ^(3)(+'),
  
  
  where  is the conformal time at late times, close to the end of 
  inflation.
  We should mention that, in the above expressions, the expectation values 
  on the left hand side are to be evaluated in the specified initial quantum 
  state, which we shall assume to be the Bunch-Davies vacuum.
  Let f_k and g_k denote the positive frequency modes (associated with
  the Bunch-Davies vacuum) in terms of which the operators _
  and γ̂_ij^ are decomposed. 
  Then, in terms of the quantities f_k and g_k, the power spectra (k) 
  and (k) can be expressed as
  
  
      (k)
         =   k^32 π^2 | f_k()|^2,
      (k)
         =    8 k^32 π^2 | g_k()|^2.
  
  
  
  
  
  
  
   §.§ Role of the intrinsic entropy perturbation
  
  
  Often the evolution of the curvature perturbations in 
  non-trivial scenarios involving departures from slow roll inflation are 
  examined in terms of the behavior of the quantity z (see, for instance,
  Refs. <cit.>).
  We find that it proves to be instructive to understand this 
  aspect from the behavior of the intrinsic entropy 
  perturbations <cit.>. 
  It is well known that, in contrast to perfect fluids, scalar fields, 
  in general, possess non-vanishing non-adiabatic pressure perturbation 
  δ p__NA or, equivalently, the intrinsic entropy
  perturbation , which are related through the expression (in this 
  context, see, for example, refs. <cit.>)
  
      δ p__NA=p'ℋ ,
  
  where p denotes the pressure associated with the background and ℋ
  =a H is the conformal Hubble parameter.
  In the case of inflation driven by a single, canonical scalar field, one can
  show that the intrinsic entropy perturbation _k associated with a given
  mode of the field can be expressed in terms of the corresponding curvature
  perturbation, say, _k, as follows <cit.>:
  
      _k'=-ł[2 a^2 p'^2  (ℋ'-ℋ^2)]̊ ł(11-c__A^2)̊ _k,
  
  where c__A=√(p'/ρ') is adiabatic speed of the scalar 
  perturbations, with ρ being the background energy density.
  It is easy to show using the equation of motion (<ref>) describing 
  the curvature perturbation that, in the super Hubble limit, the intrinsic 
  entropy perturbation _k decays as e^-2 N.
  However, it is found that, during deviations from slow roll, for modes which 
  are either about to leave or have just left the Hubble radius, the amplitude
  of the intrinsic entropy perturbation briefly increases, sourcing the
  curvature perturbation <cit.>.
  This, in turn, alters the amplitude of the curvature perturbation for modes 
  which cross the Hubble radius just before or during the departure from slow 
  roll.
  
  
  
  To demonstrate these effects, in figure <ref>, we have 
  plotted the evolution of the curvature and the intrinsic entropy perturbations 
  in the inflationary models USR2 and PI3. 
  
  
  
  In order to highlight the differences in the behavior of the modes, we have 
  plotted the evolution of the amplitudes for three modes which leave the Hubble 
  radius just prior to the start of the departure from slow roll inflation, 
  immediately after start of the period of transition, and during the middle 
  of the transition.
  We should point out that we have plotted the imaginary parts of _k
  and _k since they dominate at late times.
  Moreover, they allow us to highlight the oscillations in the sub-Hubble regime.
  The time when these oscillations cease is an indication that the modes have
  crossed the Hubble radius.
  Evidently, there is a sharp rise in the amplitude of the intrinsic entropy 
  perturbation for all the modes during the departure from slow roll inflation.
  We should add here that the corresponding real parts of _k and _k 
  behave in a roughly similar manner.
  It is the sharp rise in _k that is responsible for either an enhancement 
  or a suppression in the asymptotic (i.e. late time) amplitude of the curvature 
  perturbation, thereby leading to features in the power spectrum (for 
  related discussions in this context, also see, for instance,  
  refs. <cit.>).
  In contrast, we find that there is relatively little effect of the deviation 
  from slow roll on the evolution of the amplitude of the tensor perturbations.
  Due to this reason, the tensor power spectrum exhibits far less sharper 
  features than the scalar power spectrum.
  
  
  
  
  
   §.§ Scalar and tensor power spectra
  
  
  We shall now turn to the scalar and tensor power spectra that arise in the 
  ultra slow roll and punctuated inflationary scenarios we had discussed in 
  the last section.
  Barring the brief rise of ϵ_1 above unity in the models of punctuated
  inflation and the location of the deviations from slow roll inflation, we had 
  seen that the behavior of the first three slow roll parameters were very similar 
  in the different models of our interest (cf. figure <ref>).
  We can expect these features to be reflected in the corresponding power spectra.
  In figure <ref>, we have plotted the power spectra arising in
  all the five models, viz. USR1, USR2, PI1, PI2 and PI3. 
  
  
  
  
  
  
  We shall first point out the features in the scalar power spectra that are 
  common to all the models.
  All the models exhibit a rise in scalar power on small scales corresponding 
  to modes that leave the Hubble radius during the second stage of slow roll.
  Moreover, the location of the rise in power is determined by the time when 
  the deviation from slow roll occurs.
  This is due to the fact that, as we discussed in the previous subsection,
  it is the amplitude of the modes which exit the Hubble radius during the 
  phase of departure from slow roll that are enhanced compared to the amplitudes 
  of modes which leave during the initial phase of slow roll.
  Further, the modes that exit the Hubble radius during the epoch of ultra 
  slow roll carry the imprints of the extremely small values of the first 
  slow roll parameter and hence exhibit higher amplitudes.
  
  
  
  Let us now consider the power spectra in the models USR1 and USR2.
  The location of features in the spectra is determined by the finely 
  tuned values of parameters of the potential and the time when the modes 
  leave the Hubble radius.
  Note that both USR1 and USR2 have a definite end of inflation.
  Let us say that the pivot scale k_∗=0.05 Mpc^-1 
  leaves the Hubble radius N_∗ number of e-folds prior to the 
  end of inflation.
  For USR1 and USR2, to arrive at the power spectra
  plotted in figure <ref>, we have assumed that 
  N_∗=(50.0,56.2).
  The occurrence of a peak in the scalar power spectra 
  at small scales in these models can be easily understood if we recall 
  the behavior of the slow roll parameters in these cases.
  Note that, in slow roll inflation, the scalar spectral index  is 
  given in terms of the first two slow roll parameters as = 1 - 2 ϵ_1 - ϵ_2. 
  Though the regime of our interest does not strictly correspond to slow 
  roll dynamics, we can utilize this relation to roughly understand the 
  rise and fall of the scalar power spectra. 
  We had earlier mentioned that, as ϵ_1 decreases rapidly during
  the epoch of ultra slow roll and eventually rises from its very small 
  values, ϵ_2 changes from relatively large negative values to 
  positive values in USR1 and USR2.
  Since ϵ_1 is very small during the ultra slow roll regime, for
  modes which leave around this epoch, the spectral index  mimics the 
  behavior of -ϵ_2, changing from large positive values (corresponding 
  to an initially blue spectrum) to negative values (corresponding to a 
  red spectrum on smaller scales), leading to a peak in the power spectra.
  Clearly, we also require that the power spectra at large scales 
  are consistent with the current constraints on the scalar spectral 
  index  and the tensor-to-scalar ratio r from the CMB 
  data <cit.>.
  We find that the models USR1 and USR2 lead to (, r)= (0.945,0.015) 
  and (0.946,0.007) at the pivot scale.
  We should add a word of caution in this regard.
  The above values for  and r lie barely within the 2-σ 
  limits on the respective parameters according to the latest constraints 
  from Planck <cit.>.
  Importantly, if one were to even slightly change the values of the model 
  parameters, the features in the power spectra get considerably altered.
  In other words, there is a severe fine tuning involved in arriving at the
  desired power spectra, an aspect which is well known and has been 
  highlighted earlier (in this regard, see, for instance, 
  ref. <cit.>).
  
  
  
  Let us now turn to the power spectra arising in the punctuated inflationary
  models.
  Once again, we can understand the behavior of the spectra at small scales in 
  these cases from the relation between the scalar spectral index and the slow 
  roll parameters.
  Recall that, while PI3 has a finite duration of inflation, there exists the 
  problem of termination of inflation in the models PI1 and PI2.
  Due to this reason, as should be evident from the power spectra plotted in
  figure <ref>, the power never comes down in PI1 and PI2 because
  the eventual slow roll regime lasts for a long duration.
  However, since the evolution of the slow roll parameters in PI3 mimic their
  behavior in USR1 and USR2, the resulting scalar power spectrum exhibits a peak 
  for the same reason that we discussed above, viz. the relatively large values
  and the change in the sign of the second slow roll parameter ϵ_2.
  For the three models of PI1, PI2 and PI3, we have set N_∗ = (60.0, 60.0,
  54.5) to arrive at their respective spectra presented 
  in figure <ref>.
  We find that, for the choice of parameters that lead to COBE normalized
  scalar amplitude on large scales, the scalar spectral index and the 
  tensor-to-scalar ratio at the pivot scale prove to be 
  (, r)= (0.885, 0.580), (0.909, 0.461) and (0.944, 0.009) in 
  PI1, PI2 and PI3, respectively.
  Evidently, PI1 and PI2 are ruled out due to the large tensor-to-scalar ratio 
  (beyond the upper limits from Planck) generated on the CMB scales in these 
  models. 
  In contrast, PI3 leads to a rather small tensor-to-scalar ratio that is 
  consistent with the bounds from the Planck data and also comes close to 
  satisfying the constraints on  <cit.>. 
  As far as the extent of boosting the power on small scales and the tunability 
  of the model parameters are concerned, PI3 seems to require the same extent of 
  fine-tuning as USR1 and USR2.
  In contrast to PI3, we find that it is easier to achieve sustained amplification 
  of power over a wider range of scales in PI1 and PI2.
  But, obviously, it is achieved at the high cost that inflation does not end within 
  the desired duration, essentially making them unviable.
  Nevertheless, we believe that there are lessons to be learnt from the simpler
  models PI1 and PI2 and we will exploit the main features of these models to
  reverse engineer desired potentials in the following section.
  
  
  
  Lastly, let us make a few remarks on the tensor power spectra that we obtain 
  in the various models. 
  Note that the tensor power spectra also exhibit a step-like feature in all the 
  models, but the step is in the opposite direction as compared to the scalars,
  with the amplitude of tensors at small scales being a few orders of magnitude
  smaller than their amplitude over large scales <cit.>.
  This can be attributed to the fact that after the period of deviation from slow 
  roll, the inflaton evolves over smaller values of the field and hence smaller 
  values of the potential. 
  
  
  
  
  
   §.§ Challenges in constructing viable models
  
  
  With the experience of examining a handful of inflationary models, let us 
  briefly summarize the challenges in constructing viable and well motivated 
  models that lead to enhanced power on small scales.
  
  
  
  To begin with, we need to ensure that the scalar spectral index  and 
  the tensor-to-scalar ratio r are consistent with the cosmological data 
  over the CMB scales.
  Moreover, in order to boost the extent of PBHs formed and the amplitude of 
  the secondary GWs, we require enhanced power on small scales.
  Simultaneously, we need to make sure that inflation ends in a reasonable 
  number of (say, about 65) e-folds.
  It is found that, as one attempts to resolve one issue, say, reduce the level
  of fine tuning or permit room to shift the location of the features in the 
  scalar power spectrum, another difficulty, such as the prolonged duration of
  inflation, creeps in.
  
  
  
  We should point out here that, a given potential which admits ultra slow 
  roll inflation for a set of values of the parameters involved may permit 
  punctuated inflation for another set (in this context, see 
  appendix <ref>).
  For that reason, we should stress that the potentials themselves cannot 
  always be classified as ultra slow roll or punctuated inflationary models.
  Hence, the dichotomy of ultra slow roll and punctuated inflationary scenarios
  that we have created may be considered somewhat artificial.
  However, we find it intriguing that whenever a potential admits restoration
  of inflation after a brief interruption, it seems to naturally result in a 
  regime of ultra slow roll inflation.
  We believe that this aspect ought to be exploited to construct well motivated 
  and viable canonical, single field inflationary models that also lead to 
  enhanced PBH formation and generate secondary GWs of significant amplitudes.
  
  
  
  With the eventual aim of overcoming these difficulties in single, canonical
  scalar field models of inflation, we shall now attempt to reconstruct
  potentials that possess the desired features.
  
  
  
  
  
  § REVERSE ENGINEERING POTENTIALS ADMITTING ULTRA SLOW ROLL
  AND PUNCTUATED INFLATION
  
  
  In this section, we shall assume specific time-dependence for the first
  slow roll parameter ϵ_1 so that it leads to ultra slow roll or 
  punctuated inflation.
  With the functional form of ϵ_1(N) at hand, we shall reconstruct 
  the potentials using the equations of motion for the background and evaluate
  the resulting scalar and tensor power spectra that arise in the different 
  scenarios <cit.>.
  
  
  
  
  
   §.§ Choices of ϵ_1(N)
  
  
  We shall consider the following two forms for ϵ_1(N) which 
  lead to ultra slow roll or punctuated inflation for suitable choice 
  of the parameters involved:
  
  
  
      ϵ_1^I(N) 
         =   ł[ϵ_1a ł(1+ϵ_2a N)̊]̊ ł[1 - tanhł(N - N_1Δ N_1)̊]̊ 
      + ϵ_1b + expł(N - N_2Δ N_2)̊,
      ϵ_1^II(N) 
         =   ϵ_1^I(N) 
      + cosh^-2ł(N - N_1Δ N_1)̊.
  
  
  We find that considering a parametrization of the first 
  slow roll parameter rather than the quantity z or the scale factor a 
  proves to be much more convenient and easy to model the scenarios of our
  interest (in this context, see the recent 
  efforts <cit.>).
  The approach we adopt also allows us to easily ensure that the CMB constraints 
  on large scales are satisfied.
  The above forms of ϵ_1(N) are supposed to represent the ultra 
  slow roll and the punctuated inflationary scenarios we had discussed 
  earlier.
  For convenience, we shall hereafter refer to the reconstructed inflationary
  scenarios arising from the forms of ϵ_1(N) in eqs. (<ref>)
  and (<ref>) as RS1 and RS2, respectively.
  We shall now highlight a few points concerning the above constructions before 
  proceeding to calculate the resulting power spectra. 
  
  
  
  Consider RS1 described by ϵ_1(N) in eq. (<ref>).
  Note that the functional form contains seven parameters, viz. ϵ_1a, 
  ϵ_1b, ϵ_2a, N_1, N_2, Δ N_1 and Δ N_2.
  For suitable choices of these parameters, this form of ϵ_1(N) leads to 
  a period of slow roll followed by an epoch of ultra slow roll, before inflation 
  eventually ends, as encountered in the ultra slow models USR1 and USR2 we had
  discussed in the last section.
  While ϵ_1a and ϵ_1b determine the values of the first slow 
  roll parameter during slow roll and ultra slow roll, the parameters N_1 and N_2 
  determine the duration of these two phases.  
  Note that the first term in the functional form (<ref>) is expressed
  as a product of two parts.
  The first part involving the parameter ϵ_2a induces a small time
  dependence during the early stages.
  Such a time dependence is necessary to achieve slow roll inflation which leads 
  to scalar and tensor power spectra that are consistent with the CMB data.
  Recall that, in slow roll inflation, the scalar spectral index and the 
  tensor-to-scalar ratio are given by = 1-2 ϵ_1-ϵ_2
  and r = 16 ϵ_1, with the slow roll parameters evaluated at 
  the time when the modes cross the Hubble radius.
  For suitable choices of ϵ_1a and ϵ_2a, we find that we 
  can arrive at spectra that are consistent with the constraints on  
  and r from CMB, viz. = 0.9649 ± 0.0042 and r < 0.056 at the pivot
  scale <cit.>.
  The second part of the first term containing the hyperbolic tangent function
  aids in the transition from the slow roll to the ultra slow roll phase around 
  the e-fold N_1.
  We need to set N_1 so that all the large scale modes leave the Hubble radius
  during the first slow roll phase.
  
  
  
  The second term ϵ_1b in equation (<ref>) essentially 
  prevents the first slow parameter ϵ_1 from reducing to zero 
  beyond N_1. 
  Since ϵ_1b defines the ultra slow roll phase of the model, 
  we shall choose the parameter to be much smaller than ϵ_1a.
  The last term involving the exponential factor has been included to
  essentially ensure that ϵ_1 rapidly rises at later times, 
  crossing unity at N_2, resulting in the termination of inflation. 
  Lastly, the rapidity of the transitions from slow roll to ultra slow roll 
  and from ultra slow roll to the end of inflation are determined by the 
  parameters Δ N_1 and Δ N_2, respectively.
  In summary, since  ϵ_1a and  ϵ_2a are constrained
  by the CMB data on large scales, we have five free parameters, viz.
  ϵ_1b, N_1, N_2, Δ N_1 and Δ N_2, to construct 
  the features we desire in the scalar power spectra over small scales.
  
  
  
  Let us now turn to RS2 with ϵ_1(N) described by eq. (<ref>).
  In this case, evidently, the term involving the hyperbolic cosine function 
  has been added to the form of ϵ_1(N) in RS1. 
  This additional terms leads to a brief interruption of inflation around the 
  e-fold N_1, as is encountered in the punctuated inflationary models PI1, 
  PI2, and PI3 discussed earlier.
  
  
  
  Both the constructions of ϵ_1 above have been motivated to simplify
  the study of models containing an epoch of ultra slow roll with or without 
  punctuation and thus producing inflationary spectra with either extended or 
  localized features on small scales. 
  The advantage of these constructions is that the parameters are easy to tune,
  which allows us to directly infer the corresponding effects on the background 
  dynamics and importantly on the power spectra, unlike the specific inflationary
  models examined earlier.
  Of course, this has been possible due to the fact the reconstructions involve
  more parameters than the potentials we have considered.
  
  
  
  
  
   §.§ Reconstructed potentials and the corresponding scalar and tensor 
  power spectra
  
  
  Using the Friedmann equations and the equation of motion governing the
  inflaton, it is straightforward to show that the time evolution of the scalar 
  field ϕ(N) and the Hubble parameter H(N) can be expressed in terms 
  of the slow roll parameter ϵ_1(N) as follows:
  
  
      ϕ(N)    =   ϕ_i 
      - ∫^N_N_iṆ √(2 ϵ_1(N)),
      
      H(N)    =    H_i expł[-∫^N_N_iṆ ϵ_1(N)]̊,
  
  
  where ϕ_i and H_i are the values of the scalar field
  and the Hubble parameter at some initial e-fold N_i.
  We can use the above relations to arrive at the required background quantities
  given a functional form for ϵ_1(N).
  These background quantities can then be utilized to evaluate the resulting scalar 
  and tensor power spectra.
  It is useful to note that the potential V(N) can be expressed in terms of 
  the Hubble parameter and the first slow roll parameter as
  
      V(N) = ^2 H^2(N) ł[3-ϵ_1(N)]̊.
  
  Having obtained ϕ(N) and V(N), clearly, we can construct V(ϕ) 
  parametrically.
  
  
  
  In figure <ref>, we have plotted the two choices (<ref>)
  for ϵ_1(N) and the corresponding potentials for a small range of
  the parameter Δ N_1 that determines the duration of the transition 
  from slow roll to ultra slow roll.
  
  
  
  The parameters we have worked with in the case of the reconstructed 
  scenario RS1 are as follows: ϵ_1a=10^-4, 
  ϵ_2a=5× 10^-2,
  ϵ_1b = 10^-10, N_1 = 42, N_2 = 72 and Δ N_2 = 1.1. 
  We have varied the parameter Δ N_1 over the range (0.3345,0.7) 
  to obtain the bands of ϵ_1 and the corresponding potential in 
  the figure.
  Similarly, in the case of RS2, the parameters we have chosen to work
  with are as follows: ϵ_1a=8× 10^-5, 
  ϵ_2a=6.25× 10^-2,
  ϵ_1b= 10^-10, N_1 = 48, N_2 = 72
  and Δ N_2 = 0.8. 
  The parameter Δ N_1 has been varied over the range (0.3847,0.5) 
  to arrive at the bands of ϵ_1 and the corresponding potential. 
  We should note that the band describing the potential is more pronounced 
  in the case of RS2 than in RS1. 
  The choices for ϵ_1a and ϵ_2a have been made so that
  the resulting power spectra are consistent with the Planck constraints on
  the scalar spectral index  and the tensor-to-scalar ratio r at the 
  pivot scale that we mentioned earlier.
  For comparison, in the figure, we have also included the behavior 
  of the first slow parameter as well as the form of the potential 
  in the models USR2 and PI3.
  It should be clear that, for suitable values of the parameters, our 
  functional forms for ϵ_1(N) closely mimic the corresponding 
  behavior in these models.
  Moreover, from the parametric forms of V(ϕ) constructed numerically, 
  we have been able to determine if the reconstructed potentials in the cases 
  of RS1 and RS2 contain a point of inflection.
  At an accuracy of 0.1%, we find that the reconstructed potentials 
  indeed contain an inflection point.
  
  
  
  With the background quantities at hand, it is straightforward to compute the
  power spectra by integrating the differential equations (<ref>) for 
  the curvature and the tensor perturbations.
  In figure <ref>, we have plotted the power spectra that arise 
  in the scenarios RS1 and RS2.
  
  
  
  We have also compared the power spectra in these cases with the spectra in 
  USR2 and PI3.
  It is clear that, while the scalar power spectra from the reconstructed 
  potentials are indeed very similar to the power spectra from USR2 and PI3,
  the corresponding tensor power spectra exhibit some differences. 
  Since we shall be focusing on the observational imprints of the scalar 
  perturbations generated during inflation, we shall ignore these differences 
  for now.
  We shall make a few clarifying remarks regarding this point in the 
  concluding section.
  
  
  
  Earlier, we had emphasized the point that the models USR2
  and PI3 are highly fine-tuned and that it is difficult to move the locations 
  of the peaks in the scalar power spectra substantially without either 
  considerably affecting the duration of inflation or the spectra over the 
  CMB scales.
  In contrast, because of the presence of the additional parameters, the scenarios
  RS1 and RS2 are easier to tune and, as a result, we find that we can shift the
  location of the peak as well as broaden its width. 
  In figure <ref>, apart from the spectra in RS1 and RS2 which
  closely mimic the scalar spectra that arise in USR2 and PI3, we have 
  plotted the power spectra for two other sets of parameters which lead to 
  peaks at different locations and also exhibit a broader peak.
  These spectra have been achieved by choosing different values for the 
  parameter N_1, while keeping the other parameters fixed at the values 
  mentioned earlier. 
  To arrive at the spectra with the broader peaks in figure <ref>,
  we have set N_1=34 and 26 in the case of RS1 and N_1 = 40 and 32 in 
  the case of RS2.
  We should mention that a smaller choice of N_1 leads to a peak at a smaller
  wave number.
  Moreover, the bands associated with these two spectra correspond to the 
  variation of the parameter Δ N_1 over the domain we had mentioned 
  before.
  
  
  
  In the next two sections, we shall study the imprints of the various power 
  spectra on the formation of PBHs and the generation of secondary GWs.
  
  
  
  
  
  § FORMATION OF PBHS
  
  
  Let us begin by recalling a few essentials.
  Scales with wave numbers greater than k≃ 10^-2 Mpc^-1 
  renter the Hubble radius during the radiation dominated epoch.
  When these modes reenter the Hubble radius, the perturbations in the matter
  density at the corresponding scales collapse to form structures.
  We shall assume that the density contrast in matter characterized by the
  quantity δ is a Gaussian random variable described by the 
  probability density
  
      P(δ) 
      = 1√(2 π σ^2) expł(-δ^22 σ^2)̊,
  
  where σ^2 is the variance of the spatial density fluctuations.
  Let us assume that perturbations with a density contrast beyond a certain 
  threshold, say, δ_c, are responsible for the formation 
  of PBHs.
  In such a case, the fraction, say, β, of the density fluctuations 
  that collapse to form PBHs is described by the integral (in this context, 
  see the reviews <cit.>)
  
      β = ∫^1__δ_cδ̣ (δ)
      ≃12 ł[1-erfł(δ_c√(2 σ^2))̊]̊,
  
  where erf(z) denotes the error function.
  Note that the lower limit of the above integral is the threshold value of 
  the density contrast beyond which matter is expected to collapse to form PBHs.
  We should clarify here that the value of δ_c 
  is not unique and it is expected to depend on the amplitude of the perturbation 
  at a given scale (see refs. <cit.>;
  in this context, also see the recent discussions <cit.>).
  The choice of δ_c becomes important for the reason that the 
  extent of PBHs formed is exponentially sensitive to its value.
  In order to calculate the extent of PBHs formed, we shall work with the 
  following values of δ_c: 1/3, 0.35 and 0.4.
  
  
  
  During the radiation dominated epoch, the matter power spectrum P_δ(k) 
  and the inflationary scalar power spectrum (k) are related through the
  expression 
  
      P_δ(k) 
      = 1681 ł(kaH)̊^4 (k).
  
  The variance in the spatial density fluctuations σ^2, which determines
  the fraction β of PBHs formed [cf. eq. (<ref>)], can be expressed
  as an integral over the matter power spectrum P_δ(k).
  In order to introduce a length scale, say, R, the variance is smoothened over 
  the scale with the aid of a window function W(k R).
  The variance σ^2(R) can then be written as
  
      σ^2(R) 
      = ∫_0^∞ḳk P_δ(k)  W^2(k R),
  
  and we shall work with a Gaussian window function of the form
  W(k R) = e^-(k^2 R^2)/2.
  
  
  
  There remains the task of relating the scale R to the mass, say, M, 
  of the PBHs formed.
  Let M__H denote the mass within the Hubble radius H^-1 
  at a given time.
  It is reasonable to suppose that a certain fraction of the total mass 
  within the Hubble radius, say, M=γ  M_H, goes on to
  form PBHs when a mode with wave number k reenters the Hubble radius. 
  The quantity γ that has been introduced reflects the efficiency 
  of the collapse.
  In the absence of any other scale, it seems natural to choose k=R^-1, 
  and make use of the fact that k=a H when the modes reenter the Hubble 
  radius, to finally obtain the relation between R and M.
  One can show that R and M are related as follows:
  
      R=2^1/4γ^1/2 ł(g_∗,kg_∗,eq)̊^1/12 ł(1k_eq)̊ ł(MM_eq)̊^1/2,
  
  where k_eq is the wave number that reenters the Hubble radius 
  at the epoch of radiation-matter equality, and M_eq denotes 
  the mass within the Hubble radius at equality.
  Also, the quantities g_∗,k and g_∗,eq represent
  the number of relativistic degrees of freedom at the times of PBH 
  formation and radiation-matter equality, respectively.
  It can be easily determined that M_eq = 5.83× 10^47 kg, 
  so that we can express the above relation between R and M in terms of the
  solar mass M_⊙ as follows:
  
      R=4.72×10^-7 ł(γ0.2)̊^-1/2 ł(g_∗,kg_∗,eq)̊^1/12 ł(MM_⊙)̊^1/2 Mpc.
  
  
  
  
  On using the above arguments, we can arrive at the fraction of PBHs, say,
  , that contribute to the dark matter density today.
  The quantity (M) can be expressed as  
  
      (M) = 2^1/4 γ^3/2 β(M) ł(Ω_m h^2Ω_c h^2)̊ ł(g_∗,kg_∗,eq)̊^-1/4 ł(MM_eq)̊^-1/2,
  
  where Ω_m and Ω_c are the dimensionless 
  parameters describing the matter and cold matter densities, with the Hubble 
  parameter, as usual, expressed as H_0=100 h km sec^-1 Mpc^-1.
  In our calculations, we shall choose γ = 0.2, g_∗,k 
  = 106.75 and g_∗,eq = 3.36 and set Ω_m h^2
  = 0.14, Ω_c h^2 = 0.12, with the last two being the best
  fit values from the recent Planck data <cit.>.
  On substituting these values, one can arrive at the following expression 
  for (M):
  
      (M) 
      =ł(γ0.2)̊^3/2 ł(β(M)1.46× 10^-8)̊ ł(g_∗,kg_∗,eq)̊^-1/4 ł(MM_⊙)̊^-1/2.
  
  
  
  
  
  
  Given a primordial power spectrum (k), we can utilize the
  relations (<ref>) and (<ref>) to arrive at 
  the quantity σ^2(R).
  Then, using the relation (<ref>), we can determine σ^2
  as a function of M and utilize the result (<ref>) to 
  obtain β(M).
  With β(M) at hand, we can use the relation (<ref>) to
  finally arrive at (M) for a given inflationary scalar power 
  spectrum.
  In figure <ref>, we have plotted (M) for the models of 
  USR2, PI3, RS1, and RS2.
  
  
  
  In the figure, we have also indicated the constraints from the various 
  observations such as constraints from gravitational 
  lensing <cit.>, 
  constraints due to the limits on extragalactic background 
  photons from PBH evaporation <cit.>, constraints 
  from microlensing searches by Kepler <cit.>,
  MACHO <cit.>, EROS <cit.> 
  and OGLE <cit.>, 
  constraints from the large scale structure <cit.>,
  constraints from the CMB anisotropies due to accretion onto 
  PBHs (FIRAS and WMAP3) <cit.> and, finally, constraints
  from the dynamics of ultra-faint dwarf galaxies <cit.>. 
  (For the latest and comprehensive list of these constraints 
  and a detailed discussion, see refs. <cit.>.
  For related discussions in these contexts, also see 
  refs. <cit.>.)
  We find that, in the cases of USR2 and RS1, where the location of the peaks 
  in the scalar power spectra approximately match, the maximum values of 
  achieved are 1.5×10^-2 and 0.10, respectively. 
  For the models PI3 and RS2, when the peaks are located at roughly the same
  wave number, we similarly obtain  to be 3×10^-3 and 0.11 at their 
  respective maxima. 
  In these cases, the maxmima in (M) are located over the domain M
  ≃ 10^-16–10^-12  M_⊙.
  For peaks in the scalar power spectra that occur at smaller wave numbers in 
  the cases of RS1 and RS2, as expected, the locations of the maxima in (M) 
  shift towards larger masses of PBHs.
  Interestingly, for the power spectra in RS1 and RS2 which exhibit a broad peak 
  beginning at k≃ 10^6 Mpc^-1, there arise maxima in  
  at tens of solar masses.
  However, the corresponding maximum value of  at M≃ 10 M_⊙ is a 
  few orders of magnitude smaller than the maximum values we discussed above at 
  smaller masses. 
  This arises despite the fact the amplitude of the scalar power spectra at their
  peak is the same in all these cases.
  We believe that this result can be attributed to the dependence of  on M
  as M^-1/2 [cf. eq. (<ref>)].
  We should point out here that the shaded bands corresponding to RS1 and RS2 
  in figure <ref> indicate the range of  that can be generated 
  by varying the parameter Δ N_1 in the functional forms of ϵ_1(N) 
  [cf. eqs. (<ref>)].
  The intersection of the shaded bands with the constraints readily translate 
  to the limits on this parameter in our reconstructions RS1 and RS2. 
  We find that a smaller Δ N_1 leads to a steeper growth of power and 
  hence to a higher fraction of PBHs. 
  Therefore, for a fixed set of values for the other parameters, the constraints
  essentially restrict the rapidity of the transition of inflation from slow roll 
  to ultra slow roll epoch in our reconstructions.
  
  
  
  
  
  § GENERATION OF SECONDARY GWS
  
  
  In this section, we shall calculate the secondary power and bispectrum of 
  GWs induced by the scalar perturbations at the second order.
  
  
  
  
  
   §.§ The secondary tensor power spectrum
  
  
  Earlier, we had described the scalar and tensor perturbations at first order 
  in terms of the curvature perturbation  and the quantity γ_ij 
  (cf. subsection <ref>).
  It is well known that, at the linear order, the scalar and tensor perturbations 
  evolve independently, with their evolution being governed by the corresponding
  equations of motion, viz. eqs. (<ref>).
  However, one finds that, at the second order, the tensor perturbations are 
  sourced by quadratic terms involving the first order scalar perturbations 
  (for early discussions in this context, see for instance, 
  refs. <cit.>).
  These contributions due to the scalar perturbations become important particularly
  when the amplitude of the scalar power spectrum is boosted over small scales such
  as in the situations leading to enhanced formation of PBHs.
  In this subsection, we shall calculate the dimensionless density parameter
  associated with the GWs, say, , generated due to the scalar perturbations 
  in the different models and scenarios of interest.
  
  
  
  Let us begin by outlining the primary steps towards the calculation of (f), 
  where f is the frequency associated with the wave number k.
  We shall start with the following perturbed metric:
  
      ṣ^2 = a^2(η) ł{-ł(1+2 Φ)̊ η̣^2
      + ł[ł(1-2 Ψ)̊ δ_ij+12 h_ij]̊ x̣^i x̣^j}̊,
  
  where Φ and Ψ are the Bardeen potentials describing the scalar 
  perturbations at the first order, while the quantity h_ij represents 
  the second order tensor perturbations.
  We should clarify that we have denoted the second order tensor perturbation
  as h_ij in order to distinguish them from the first order tensor
  perturbations γ_ij which we had introduced earlier.
  The transverse and traceless nature of the tensor perturbations implies that 
  ^i h_ij=0 and h^i_i=0.
  In our discussion below, we shall assume that anisotropic stresses are absent
  so that Φ=Ψ.
  
  
  
  The tensor perturbations h_ij can be decomposed in terms of the Fourier modes, 
  say, h_ k, as 
  
      h_ij(η,)
      =∫^̣3(2 π)^3/2 ł[e_ij^+()  h_^+(η)
      +e_ij^×()  h_^×(η)]̊  e^i ·,
  
  where e_ij^+() and e_ij^×() denote the polarization tensors
  which have non-zero components in the plane perpendicular to the direction of 
  propagation, viz. .
  The polarization tensors e_ij^+() and e_ij^×() can be 
  expressed in terms of the set of orthogonal unit vectors (e(),
  e̅(),) in the following manner (see, for instance, the
  review <cit.>):
  
  
      e_ij^+()
         =   1√(2) ł[e_i() e_j()
      -e̅_i() e̅_j()]̊,
      
      e_ij^×()
         =   1√(2) ł[e_i() e̅_j()
      +e̅_i()  e_j()]̊.
  
  
  The orthonormal nature of the vectors e() and e̅() lead to the 
  normalization condition: e_ij^λ()  e^λ',ij()
  =δ^λλ', where λ and λ' can be 
  either + or ×.
  
  
  
  The equation of motion governing the Fourier modes h_ can be arrived at using 
  the second order Einstein equations describing the tensor perturbation h_ij
  and the Bardeen equation describing the scalar 
  perturbation Ψ at the first order (see, for example, 
  refs. <cit.>; for recent discussions, see 
  refs. <cit.>).
  One finds that the equation governing h_ can be written as
  
      h_^λ”+ 2 ℋ h_^λ'+ k^2  h_^λ
      =S_^λ
  
  with the source term S_^λ being given by
  
      S_^λ(η)
         =    4 ∫^̣3 /(2 π)^3/2  e^λ(,) {2 Ψ_(η) Ψ_-(η)
            + 43  (1+w) ℋ^2  ł[Ψ_'(η)+ H Ψ_(η)]̊ ł[Ψ_-'(η)+ H Ψ_-(η)]̊ },
  
  where, evidently, Ψ_ represents the Fourier modes of the Bardeen potential,
  while H and w denote the conformal Hubble parameter and the equation of
  state parameter describing the universe at the conformal time η.
  Also, for convenience, we have defined the quantity e^λ(,)=
  e^λ_ij() p^i p^j.
  While discussing the formation of PBHs earlier, we had assumed that the scales of 
  our interest reenter the Hubble radius during the epoch of radiation domination.
  In such a case, we have w=1/3 and H=1/η.
  Moreover, during radiation domination, it is well known that we can express the 
  Fourier modes Ψ_ of the Bardeen potential in terms of the inflationary 
  Fourier modes _ of the curvature perturbations generated during inflation
  through the relation
  
      Ψ_(η)=23 (k η) _,
  
  where (k η) is the transfer function given by
  
      (k η)=9ł(k η)̊^2 ł[ sinł(k η/√(3))̊k η/√(3)
      - cosł(k η/√(3))̊]̊.
  
  Utilizing the Green's function corresponding to the tensor modes during 
  radiation domination, we can express the inhomogeneous contribution to 
  h_^λ as <cit.>
  
      h_^λ(η) 
         =   4/9 k^3 η ∫^̣3 (2 π)^3/2  
      e^λ(,) _ _- ł[_cł(pk,|-|k)̊  cosł(k η)̊
      +_sł(pk,|-|k)̊  sinł(k η)̊]̊,
  
  where the quantities _c(v,u) and _s(v,u) are described by the integrals
  
  
      _c(v,u)   =   -4 ∫_0^∞ τ̣ τ  sin τ {2 (v τ) (u τ)
      +ł[(v τ)+v τ _vτ(v τ)]̊ ł[(u τ)+u τ _uτ(u τ)]̊},
      _s(v,u)   =    4 ∫_0^∞ τ̣ τ  cos τ {2 (v τ) (u τ)
      +ł[(v τ)+v τ _vτ(v τ)]̊ ł[(u τ)+u τ _uτ(u τ)]̊},
  
  
  with _z=/ẓ.
  The above integrals can be carried out analytically and they are given by
  
  
      _c(v,u)    =    -27 π4 v^3 u^3 Θł(v+u-√(3))̊  (v^2+u^2-3)^2,
      _s(v,u)    =    -274 v^3 u^3  (v^2+u^2-3) ł[4 v u+ (v^2+u^2-3)  log |3-(v-u)^23-(v+u)^2|]̊,
  
  
  where Θ(z) denotes the theta function.
  It is useful to note that _c,s(v,u) =_c,s(u,v).
  
  
  
  The power spectrum of the secondary GWs, say, (k,η), generated
  due to the second order scalar perturbations can be defined as follows:
  
      ⟨ h_^λ(η) h_'^λ'(η)⟩
      =2 π^2k^3 _h(k,η) δ^(3)(+') δ^λλ'.
  
  Note that h_^λ involves products of the Fourier modes _ 
  and _- of the curvature perturbations generated during inflation 
  [cf. eq. (<ref>)].
  Evidently, the power spectrum (k) of the secondary GWs will involve 
  products of four such variables.
  Since, the quantity _ is a Gaussian random variable, we can express
  the four-point function in terms of the two-point functions or, equivalently,
  the inflationary scalar power spectrum (k) [cf. eq. (<ref>)]
  as
  
      (k,η)
         =   481 k^2 η^2∫_0^∞ṿ ∫_| 1-v|^1+vụ ł[4 v^2-(1+v^2-u^2)^24 u v]̊^2 (k v) (k u)
            × ł[_c(u,v)  cosł(k η)̊+_s(u,v)  sinł(k η)̊]̊^2.
  
  We shall now choose to average (k,η) over small time scales so that 
  the trigonometric functions in the above expressions are replaced by their
  average over a time period.
  In such a case, only the overall time dependence remains, leading 
  to <cit.>
  
      (k,η)   =   281 k^2 η^2∫_0^∞ṿ ∫_| 1-v|^1+vụ ł[4 v^2-(1+v^2-u^2)^24 u v]̊^2 (k v) (k u)
            × ł[_c^2(u,v)+_s^2(u,v)]̊,
  
  where the line over (k,η) implies that we have averaged over 
  small time scales.
  The energy density of GWs associated with a Fourier mode corresponding to 
  the wave number k at a time η is given by <cit.>
  
      ρ__GW(k,η)
      = ^28 ł(k/a)̊^2 (k,η).
  
  The corresponding dimensionless density parameter (k,η) can
  be defined in terms of the critical density ρ_cr(η) 
  as <cit.>
  
      (k,η)
      =ρ__GW(k,η)ρ_cr(η)
      =124 ł(kℋ)̊^2 (k,η).
  
  
  
  
  Note that the dimensionless density parameter (k,η) above has been 
  evaluated during the radiation dominated epoch.
  Once the modes are inside the Hubble radius, the energy density of GWs 
  decay just as the energy density of radiation does.
  Upon utilizing this point, we can express (k) today in terms of the 
  above (k,η) as follows:
  
      h^2 (k)
         =   ł(g_∗,kg_∗,0)̊^-1/3 Ω_r h^2 (k,η)
         ≃    
      1.38×10^-5 ł(g_∗,k106.75)̊^-1/3 ł(Ω_r h^24.16×10^-5)̊ (k,η),
  
  where Ω_r and g_∗,0 denote the dimensionless energy 
  density of radiation and the number of relativistic degrees of freedom today.
  We should point out here that, since ℋ∝η^-1 during 
  radiation domination and (k,η)∝η^-2, the quantity 
  (k,η) in the expression (<ref>) is actually independent 
  of time.
  Moreover, the observable parameter today is usually expressed as a function 
  of the frequency, say, f, which is related to the wave number k as
  
      f = k/2 π 
      = 1.55×10^-15 ł(k1 Mpc^-1)̊ Hz.
  
  
  
  
  In figure <ref>, we have plotted the quantity (f) arising
  in the models USR2 and PI3 as well as the reconstructed scenarios RS1 
  and RS2.
  
  
  
  In the figure, we have also included the sensitivity curves associated  
  with the various current and forthcoming observatories, viz. PTA and 
  the Square Kilometre Array (SKA) <cit.>, LISA <cit.>, 
  MAGIS-100 <cit.>, 
  BBO <cit.>, 
  DECIGO <cit.>,
  ET <cit.>, 
  advanced LIGO + Virgo <cit.>
  and CE <cit.>. 
  (For a summary of the sensitivity curves and their updated versions, see
  ref. <cit.> and the associated web-page.)  
  We should mention here that the estimated sensitivity curves have been arrived
  at assuming a power law spectrum (the so-called `power-law integrated curves') 
  over the bands of interest. 
  These sensitivities are expected to be achieved by integrating over frequency 
  in addition to integrating over time <cit.>.
  It should be evident from the figure that the strength of the GWs generated 
  in the models and scenarios we have examined here is significant enough to
  be detectable by one or more of these observatories.
  Recall that, spectra arising in the scenarios RS1 and RS2
  with broad peaks starting from a wave number of about 10^6 Mpc^-1 
  had led to PBHs with tens of solar masses.
  It should be clear from figure <ref> that the constraints from PTA 
  on  already rule out such spectra for certain values of Δ N_1.
  
  
  
  
  
   §.§ The secondary tensor bispectrum
  
  
  In this section, we shall evaluate the secondary tensor bispectrum generated 
  in the inflationary models and scenarios of our interest.
  The secondary tensor bispectrum, say, 
  ^λ_1λ_2λ_3_h(,,) is defined as
  
      ł⟨ h_^λ_1(η)  h_^λ_2(η) 
      h_^λ_3(η)⟩̊=(2 π)^3 _h^λ_1λ_2λ_3(,,,η) δ^(3)(++).
  
  We can evaluate the above tensor bispectrum during the radiation dominated 
  era by using the expression (<ref>) for h^λ_(η).
  As we had discussed, h^λ_(η) is quadratic in the 
  Gaussian variables _.
  Therefore, obviously, the bispectrum 
  _h^λ_1λ_2λ_3(,,,η) will involve 
  six of these variables.
  Upon utilizing Wick's theorem applicable to Gaussian random variables, one 
  can show that the tensor bispectrum consists of eight terms all of which lead 
  to the same contribution <cit.>.
  For convenience, we shall define 
  G_h^λ_1λ_2λ_3(,,,η)
  =(2 π)^-9/2 _h^λ_1λ_2λ_3(,,,η) 
  and hereafter refer to G_h^λ_1λ_2λ_3(,,,η) 
  as the secondary tensor bispectrum.
  We find that the secondary tensor bispectrum can be expressed as
  
      G_h^λ_1λ_2λ_3(,,,η)
         =   ł(8 π/9)̊^3 1(k_1 k_2 k_3 η)^3
            × ∫^̣3 _1  e^λ_1(_1,_1) 
      e^λ_2(_2,_2) e^λ_3(_3,_3) (p_1)p^3_1 (p_2)p^3_2 (p_3)p^3_3
            × Jł(p_1k_1,p_2k_1,η)̊ 
      Jł(p_2k_2,p_3k_2,η)̊ 
      Jł(p_3k_3,p_1k_3,η)̊,
  
  where _2 = _1-, _3 = _1 + and, for convenience, we 
  have set
  
      Jł(p_1k_1,p_2k_1,η)̊
      =_cł(p_1k_1,p_2k_1)̊  cosł(k_1 η)̊
      +_sł(p_1k_1,p_2k_1)̊  sinł(k_1 η)̊,
  
  with _c(v,u) and _s(v,u) given by eqs. (<ref>).
  In a manner partly similar to the case of the secondary tensor 
  power spectrum, we shall replace the trigonometric functions 
  by their averages so that the function J(x,y,η) is instead 
  given by
  
      J̅(v,u) =1√(2) ł[_c^2(v,u)+_s^2(v,u)]̊^1/2.
  
  
  
  
  Our aim in this work is to understand the amplitude of the secondary tensor 
  bispectrum generated due to the scalar perturbations for modes that reenter 
  the Hubble radius during the radiation dominated era. 
  For simplicity, we shall restrict our analysis to the equilateral limit of 
  the bispectrum so that k_1=k_2=k_3=k.
  In order to determine the integrals involved in the expression (<ref>),
  we shall choose a specific configuration for the vectors ,  and .
  We shall assume that the vectors lie in the x-y-plane with  oriented
  along the negative x-direction.
  In such a case, we find that the vectors (, , ) in the equilateral 
  limit are given by
  
      =ł(k/2, √(3) k/2, 0)̊,   = ł(k/2, -√(3) k/2, 0)̊,   = (-k, 0, 0).
  
  We shall choose _1 =(p_1x, p_1y, p_1z) so that, since 
  _2 = _1- and _3 = _1 +, we have
  
      _2 = ł(p_1x-k/2, p_1y-√(3) k/2, p_1z)̊,  _3 = (p_1x-k, p_1y, p_1z).
  
  We find that such a choice of Cartesian coordinates
  proves to be convenient to carry out the integrals involved than the 
  cylindrical polar coordinates that have been adopted
  earlier <cit.>.
  Therefore, the tensor bispectrum in the equilateral limit 
  G_h^λ_1λ_2λ_3(k) can be written as
  
      k^6 G_h^λ_1λ_2λ_3(k,η) 
         =   ł(8 π/9 √(2))̊^3 1(k η)^3
            × ∫_-∞^∞p̣_1x ∫_-∞^∞p̣_1y ∫_-∞^∞p̣_1z e^λ_1(_1,_1) 
      e^λ_2(_2,_2) e^λ_3(_3,_3)
            × (p_1)p^3_1 (p_2)p^3_2 (p_3)/p^3_3 J̅ł(p_1/k,p_2/k)̊ J̅ł(p_2/k,p_3/k)̊ J̅ł(p_3/k,p_1/k)̊.
  
  The factors e^λ(,) involving the polarization tensor can 
  be readily evaluated for our configurations of (,,) and 
  (_1,_2,_3) (for details, see appendix <ref>).
  Since λ can be + or ×, clearly, the tensor bispectrum
  G_h^λ_1λ_2λ_3(k,η) has eight components.
  However, we find that e^×(,) is odd in p_1z 
  [cf. eqs. (<ref>)].
  As a result, the tensor bispectrum proves to be non-zero only for the 
  following combinations of (λ_1λ_2λ_3): (+++), 
  (+××), (×+×) and (××+).
  Also, note that the integral above describing the tensor bispectrum in 
  the equilateral limit is symmetric under the simultaneous interchange 
  of λ_1↔λ_2, ↔ and
  _1↔_2.
  This implies that, in the equilateral limit of interest, the tensor 
  bispectrum for the three components (+××), (×+×) 
  and (××+) are equal.
  Hence, we are left with only G_h^+++(k,η) and, say, 
  G_h^+××(k,η) to evaluate.
  
  
  
  We proceed to numerically evaluate G_h^+++(k) and G_h^+××(k)
  in the situations of our interest, viz. namely USR2, PI3, RS1, and RS2. 
  Because the scalar power spectra in these cases exhibit a localized maxima, 
  we restrict our evaluation of the tensor spectrum to the range of wave numbers
  around the peak.
  We find that the integrand in eq. (<ref>) exhibits a maximum around
  |_1|≃ k and, beyond that, it quickly decreases in all the
  three directions of integration.
  In fact, the contributions to the integral prove to be negligible 
  for |_1|≳ 100 k.
  So, we choose the limits for our integrals over p_1x, p_1y 
  and p_1z to be (-10^3 k, 10^3 k).
  
  
  
  In order to understand the behavior of the tensor bispectrum, we shall calculate 
  the dimensionless quantity referred to the shape function, say, _h(k), which 
  is defined as <cit.> 
  
      S_h^λ_1λ_2λ_3(k) 
      = k^6 G^λ_1λ_2λ_3_h(k,η)√(_h^3(k,η)).
  
  Note that, in this expression, both the quantities 
  k^6 G^λ_1λ_2λ_3_h(k) and _h^3(k) are 
  dimensionless.
  Moreover, the overall dependence on time cancels leading to a shape function that 
  is time-independent.
  In figure <ref>, we have plotted the shape functions S_h^+++(k) and
  S_h^+××(k) for the four cases of interest, viz. USR2, PI3, RS1
  and RS2.
  
  
  
  We find that the amplitude of S_h(k) for a given model or scenario is maximum 
  around the wave number where the scalar power spectrum exhibits a peak.
  This is true for both the cases of  S_h^+++(k) and S_h^+××(k)
  though there is a certain asymmetry in the behavior of the functions about the peak.
  Note that the amplitude of S_h(k) remains large over large wave numbers, while 
  it quickly reduces to small values at smaller wave numbers.
  In fact, this behavior should not come as a surprise since such a behavior was also
  encountered in the case of (f) (cf. figure <ref>).
  It is interesting to note that S_h^+++(k) and S_h^+××(k) 
  settle down to about 10 and -250, respectively, at large wave numbers.
  Recall that the secondary tensor bispectra and hence the shape functions we have 
  illustrated in figure <ref> have been evaluated during the radiation 
  dominated epoch, when the modes are well inside the Hubble radius. 
  They will have to be evolved until today to examine the
  corresponding observational imprints which may possibly be detected by 
  upcoming missions such as, say, LISA and PTA (in this context, see 
  ref. <cit.>; also see
  refs. <cit.>).
  
  
  
  
  
  § CONTRIBUTIONS TO PBH FORMATION AND SECONDARY GWS FROM 
  SCALAR NON-GAUSSIANITIES
  
  
  Until now, we have focused on the imprints of the scalar power spectrum 
  on the extent of PBHs formed and the generation of secondary GWs.
  Clearly, if the scalar non-Gaussianities prove to be large in a given 
  inflationary model, it seems plausible that they would significantly alter
  the observables ,  and S_h <cit.>.
  To understand the possible effects of non-Gaussianities on ,  
  as well as S_h, in this section, we shall first calculate the scalar 
  bispectrum and thereby the corresponding non-Gaussianity parameter  
  in the two inflationary models USR2 and PI3 and the reconstructed scenarios 
  RS1 and RS2.
  We shall then discuss the corresponding contributions from the scalar 
  bispectrum to ,  and S_h.
  
  
  
  
  
   §.§ Evaluating the scalar bispectrum
  
  
  The scalar bispectrum is the three point function of the curvature 
  perturbation in Fourier space, and it is defined in terms of the
  operator _ that we had introduced earlier as 
  follows <cit.>:
  
      ⟨_() _() _()⟩ 
      =(2 π)^3 __ S(,,) δ^(3)(++).
  
  Recall that,  is a time close to the end of inflation and, in this expression,
  the expectation value on the left hand side is to evaluated in the perturbative 
  vacuum <cit.>.
  Note that the three wave vectors (,,) form the edges of a triangle.
  For convenience, we shall hereafter set
  
      __ S(,,)=
      (2 π)^-9/2  G(,,)
  
  and refer to G(,,) as the scalar bispectrum.
  
  
  
  The so-called Maldacena formalism is the most complete approach to 
  evaluate the scalar bispectrum in a given 
  inflationary model <cit.>.
  In this approach, one first obtains the third order action governing 
  the curvature perturbation.
  With the third order action at hand, the scalar bispectrum is
  evaluated using the standard rules of perturbative quantum field 
  theory.
  For the case of inflation driven by a single, canonical scalar field,
  the third order action is found to consist of six bulk terms terms,
  apart from the boundary terms <cit.>.
  One can show that the scalar bispectrum G(,,) generated
  by such an action can be expressed as follows (see, for instance, 
  refs. <cit.>; in this context, also see
  ref. <cit.>):
  
      G(,,) 
         =   ∑_C=1^7  G__C(,,)
         =   ^2 ∑_C=1^6 [f_k_1()  f_k_2() f_k_3() __C(,,)+complex conjugate]
      + G_7(,,),
  
  where, as we discussed earlier, f_k are the positive frequency Fourier
  modes of the curvature perturbation.
  Amongst the seven terms in the above expression for the scalar bispectrum, 
  the first six correspond to the bulk terms in the third order action, whereas
  the seventh arises due to a boundary term, and it is usually absorbed through
  a field redefinition <cit.>.
  The quantities __C(,,), with C=(1,6), are integrals  
  associated with the bulk terms in the action and, as one can expect, apart 
  from the background quantities, they involve the modes f_k and its 
  derivative f_k'.
  (We have listed these integrals explicitly in appendix <ref>.)
  The seventh term G_7(,,) that arises due to the contribution 
  from a boundary term can be expressed as <cit.>
  
      G_7(,,)
         =    -i ^2 ł[f_k_1() f_k_2() f_k_3()]̊
            × [a^2ϵ_1ϵ_2 
      f_k_1^∗(η) f_k_2^∗(η) f_k_3'^∗(η) 
      +two permutations]_η_i^ 
      + complex conjugate,
   
  where  is the time when the initial conditions are imposed on the
  scalar perturbations.
  We should mention that the remaining boundary terms do not contribute 
  in the scenarios of our interest.
  
  
  
  As in the case of the scalar power spectrum,  due to the deviation
  from slow roll, it proves to be difficult to evaluate the scalar 
  bispectrum analytically in the inflationary models of interest.
  Therefore, we resort to numerics.
  There now exists a standard procedure to numerically compute the scalar 
  bispectrum in inflationary models involving a single, canonical 
  scalar field <cit.>.
  Recall that, in the case of the power spectrum, it is adequate to 
  impose the Bunch-Davies initial conditions on the modes when they 
  are sufficiently inside the Hubble radius.
  Apart some special situations wherein the boundary conditions 
  may need to be imposed deeper inside the Hubble, one often imposes 
  the conditions when k/(a H)≃ 10^2.
  Since the amplitude of the scalar as well as tensor perturbations 
  freeze when they are adequately outside the Hubble radius, say, 
  when k/(a H)≃ 10^-5, one can evaluate the power spectra
  at such a time for the different modes.
  Note that, in order to arrive at the bispectrum we need to carry
  out integrals which involve the background quantities, the scalar 
  modes f_k and its time derivative f_k' 
  [cf eqs. (<ref>) and (<ref>)].
  These integrals need to be carried out from a time  when the 
  initial conditions are imposed on the modes until the late 
  time  towards the end of inflation.
  We had mentioned that the amplitudes of the modes freeze soon after
  they leave the Hubble radius.
  Due to this reason, one finds that, the super-Hubble contributions
  to the scalar bispectrum prove to be negligible <cit.>.
  Therefore, one can carry out the integrals from the time when 
  k/(a H)≃ 10^2 to the time when k/(a H)≃ 10^-5.
  However, since the bispectrum involves three modes, in general, one 
  needs to integrate from the time when the smallest of the three wave 
  numbers is well inside the Hubble radius to the time until when the 
  largest of the wave numbers is sufficiently outside.
  Moreover, in order to choose the correct perturbative vacuum, one
  has to impose a cut-off in the sub-Hubble regime <cit.>.
  We impose a democratic (in wave number) cut-off of the 
  form exp-ł[κ (k_1+k_2+k_3)/(3 a H)]̊, 
  where κ is a positive definite and small
  quantity <cit.>.
  In fact, such a cut-off aids in the efficient numerical computation 
  of the integrals involved.
  One can choose a suitable value of κ depending on how deep
  from inside the Hubble radius the integrals are to be carried out.
  
  
  
  
  
   §.§ Amplitude and shape of 
  
  
  The non-Gaussianity parameter, say, (,,), corresponding 
  to the scalar bispectrum is defined as (see, for 
  instance, Refs. <cit.>)
  
      (,,)
          =   -10/3 1/ł(2 π)̊^4 k_1^3  k_2^3  k_3^3 
      G(,,) [k_1^3 (k_2) (k_3) 
      + two permutations]^-1,
  
  where (k) denotes the scalar power spectrum [cf. eq. (<ref>)].
  With the scalar power and bispectra at hand, evidently, it is straightforward 
  to arrive the non-Gaussianity parameter  for a given model.
  
  
  
  Based on prior experience, we would like to emphasize a few points 
  concerning the expected shape and amplitude of the scalar 
  bispectrum before we go on to present the results for  in the 
  different models and scenarios we have introduced earlier.
  As is well known, in slow roll inflationary models involving a single,
  canonical scalar field, the scalar non-Gaussianity parameter  
  proves to be of the order of the first slow roll 
  parameter ϵ_1 <cit.>.
  In other words, the parameter  is typically of the order of 10^-2 
  or smaller in such situations.
  Moreover, the bispectrum is found to have an equilateral shape, with 
  the  parameter slightly peaking when k_1=k_2=k_3 (in this 
  context, see, for instance, ref. <cit.>).
  However, when departures from slow roll occur, the non-Gaussianity
  parameter  can be expected to be of the order of unity or larger,
  depending on the details of the background dynamics.
  Further, in contrast to the slow roll case, wherein there is only a 
  weak dependence of the parameter  on scale, when departures from
  slow roll occur, the parameter turns out to be strongly scale dependent.
  Needless to say, we can expect that the non-Gaussianity parameter 
  to be relatively large as well as strongly scale dependent in the
  situations of our interest.
  
  
  
  Let us now discuss the results we obtain in the different models we 
  have introduced.
  In order to illustrate the complete shape of the bispectrum, the 
  non-Gaussianity parameter  is usually presented as a density plot in, 
  say, the (k_3/k_1)-(k_2/k_1)-plane <cit.>.
  It proves to be a bit of a numerical challenge to compute the complete 
  shape of the bispectrum across the wide range of wave numbers over which 
  we have evaluated the power spectra. 
  As a result, we shall focus on the amplitude of  in the equilateral 
  and the squeezed limits, i.e. when k_1=k_2=k_3=k and when k_1→ 0, 
  k_2≃ k_3=k, respectively.
  It is easier to calculate the scalar bispectrum in the equilateral
  limit as we just need to follow the evolution of one mode at a time.
  To arrive at the scalar bispectrum in the squeezed limit, we shall
  set k_2=k_3=k and choose k_1=10^-3 k.
  We have confirmed that our results are robust against choosing a smaller
  value of k_1.
  Before we go to illustrate the amplitude and shape of the non-Gaussianity 
  parameter , let us understand the behavior of the scalar 
  bispectrum G(,,) itself.
  In figure <ref>, we have plotted the scalar bispectra that arise
  in the equilateral and squeezed limits in the models of USR2 and PI3.
  
  
  
  We would like to highlight a few aspects regarding the amplitude and shape of 
  the bispectra.
  Note that the scalar bispectra have roughly the same shape in the equilateral
  and squeezed limits.
  Also, they closely resemble the corresponding scalar power spectra and, in 
  particular, they exhibit a dip and a peak around the same locations
  (cf. figure <ref>).
  Moreover, at small scales, the scalar bispectra have a larger amplitude in the
  equilateral limit than in the squeezed limit.
  Further, in the equilateral limit, the scalar bispectra have almost the same
  amplitude as the power spectra near the peak.
  
  
  
  Let us now understand the behavior of the non-Gaussianity parameter .
  In figures <ref> and <ref>, we have plotted 
  the behavior of the  parameter in the equilateral and squeezed limits 
  over a wide range of wave numbers in the models USR2 and PI3 as well as the 
  scenarios RS1 and RS2.
  
  
  
  
  
  The following points are evident from the two figures.
  Firstly, in the equilateral limit, the non-Gaussianity parameter  proves 
  to be fairly large (of the order of 10^1–10^4) over a small range of wave 
  numbers.
  In fact, the  exhibit an upward spike in their amplitude around exactly 
  the same wave numbers wherein the scalar power spectra exhibit a 
  downward spike (cf. figures <ref> and <ref>).
  Since the definition of the parameter  [cf. eq. (<ref>)] contains
  the scalar power spectrum in the denominator, the upward spike can be
  partly attributed to the downward spike in the power spectrum.
  If we ignore the large spike, we find that ≃ 1–10 around these 
  wave numbers.
  It is worth noting that these wave numbers correspond to those modes which
  leave the Hubble radius just prior to or during the transition from the 
  slow roll to the ultra slow roll regime.
  In contrast, the non-Gaussianity parameter  proves to be relatively 
  small (at most of order unity) over wave numbers where the scalar power 
  spectra exhibit their peak.
  However, we should clarify that, though the value of  is smaller than unity 
  around this domain, it is considerably larger than its typical value in slow roll
  inflation (of about 10^-2, such as over the CMB scales in our models).
  For instance, in USR2 and PI3, we find that, in the equilateral limit,  is
  about -0.37 and -0.44, respectively, near the locations of the peak in the 
  power spectra.
  This can be attributed to the large value of ϵ_2 during the ultra slow 
  roll regime.
  Secondly, in the squeezed limit, the scalar bispectrum is expected to 
  satisfy the so called consistency condition wherein it can be completely
  expressed in terms of the scalar power spectrum <cit.>.
  This translates to the condition ^_CR(k)=(5/12)  [(k)
  -1] in the squeezed limit, where (k)-1=ln (k)/ln  k is the scalar spectral index.
  In figures <ref> and <ref>, apart from 
  plotting  in the squeezed limit, we have also plotted the quantity
  ^_CR obtained from the scalar spectral index.
  We should add that we have also examined the validity of
  the consistency relation more closely by working with a smaller k_1.
  We find that the consistency condition is indeed satisfied even when there 
  arise strong features in the scalar power spectrum in all the scenarios of 
  our interest (in this context, however, see appendix <ref>).
  Therefore, in the squeezed limit, we find that  is at most of order
  unity around the peaks of the scalar power spectra.
  
  
  
  It seems important that we clarify a point regarding the 
  validity of the consistency condition at this stage of our discussion. 
  One may be concerned if the period of ultra slow roll, with its large value 
  of ϵ_2, could lead to a violation of the consistency condition over 
  wave numbers that leave the Hubble radius during this epoch (in this context, 
  see refs. <cit.>).
  Recall that the amplitude of scalar modes over a certain range of wave numbers 
  are modified to some extent during the transition from slow roll to ultra 
  slow roll (cf. figure <ref>).
  However, since, in the cases of our interest, the epoch of ultra slow roll
  ends leading to the eventual termination of inflation, the amplitude of the
  scalar modes asymptotically freeze at sufficiently late times (for further 
  details, see appendix <ref>; in this context, also see 
  refs. <cit.>).
  Due to this asymptotic behavior of the scalar modes, it should not come as a 
  surprise that the consistency condition is satisfied in the models and scenarios 
  of our interest despite the phase of ultra slow roll 
  (for very recent 
  discussions in this context, see Refs. <cit.>).
  
  
  
  
  
   §.§ Imprints of  on  and 
  
  
  Recall that the observationally relevant dimensionless, scalar non-Gaussianity
  parameter  is usually introduced through the following relation (see
  ref. <cit.>; also see refs. <cit.>):
  
      (η, )=^G(η, )
      -35  ł[^_G(η, )]̊^2
  
  where ^G denotes the Gaussian contribution. 
  In Fourier space, this relation can be written as (see, for instance,
  ref. <cit.>)
  
      _
      =^G_
      -35  ∫^̣3(2 π)^3/2 ^G_ ^G_-.
  
  If one uses this expression for _ and evaluates the corresponding
  two-point correlation function in Fourier space, one obtains 
  that <cit.> 
  
      ⟨_ _'⟩
      =2 π^2k^3 δ^(3)( +') ł[(k)+ł(35)̊^2 k^32 π ^2 ∫ d^3 (p)p^3 ł(|-|)̊|-|^3]̊,
  
  where (k) is the original scalar power spectrum defined in the Gaussian 
  limit [cf. eq. (<ref>)], while the second term represents the 
  leading non-Gaussian correction.
  We find that we can write the non-Gaussian correction to the scalar power 
  spectrum, say, (k), as follows:
  
      (k)   =   ł(35)̊^2 ^2 ∫_0^∞ṿ∫_| 1-v|^1+vụv^2 u^2 (k v) (k u)
         =   ł(125)̊^2 ^2 ∫_0^∞ṣ∫_0^1ḍ(s^2-d^2)^2 [k (s+d)/2] [k (s-d)/2].
  
  
  
  
  Since we have evaluated the scalar non-Gaussianity parameter in the inflationary
  models of our interest, we can now calculate the non-Gaussian corrections (k) 
  to the scalar power spectrum and the corresponding modifications to , 
  and S_h.
  However, before we do so, we need to clarify an important point.
  In introducing the scalar non-Gaussianity parameter through the
  relation (<ref>), it has been assumed that  is local, i.e. 
  it is independent of the wave number <cit.>.
  In contrast, the parameter  proves to be strongly scale dependent in all 
  the situations we have considered.
  In order to be consistent with the fact that the  in eq. (<ref>)
  is local, we shall consider the squeezed limit of the parameter (in this 
  context, also see the discussions in ref. <cit.>).
  Moreover, in the expression (<ref>) for (k), we shall assume that 
   is dependent on the wave number k, with k_2=k_3≃ k and k_1≪ k 
  to be consistent with the squeezed limit. 
  In figure <ref>, we have plotted the original Gaussian power spectrum
  as well the modified power spectrum including the non-Gaussian corrections (k).
  
  
  
  Recall that the non-Gaussianity parameter  had contained sharp spikes 
  around the wave numbers where the Gaussian scalar power spectra had
  exhibited a downward spike (cf. figures <ref> 
  and <ref>). 
  While evaluating the modified power spectra, we have regulated the maximum value
  of these spikes to be ||≃ 100.
  Evidently, the non-Gaussian corrections to the scalar power spectrum are 
  insignificant.
  This can be attributed to the fact that the peaks in the original power 
  spectrum (k) and the non-Gaussianity parameter  are located 
  at different wave numbers.
  Therefore, we find the corresponding modifications to ,  and 
  S_h are insignificant as well.
  This conclusion can also be understood from the fact the amplitude of the
  dimensionless bispectrum in the squeezed limit is considerably smaller 
  than the amplitude of the scalar power spectrum around its peak 
  (cf. figure <ref>).
  
  We should clarify a particular point regarding the non-Gaussian corrections 
  we have calculated in this section.
  Note that we have calculated the cubic order non-Gaussian corrections to the 
  power spectrum.
  This method proves to be adequate to examine the imprints of non-Gaussianities 
  on the dimensionless energy density  describing the secondary GWs. 
  However, the approach does not completely account for the effects of 
  non-Gaussianities on the fraction  of PBHs produced (for an early
  discussion on the topic, see ref. <cit.>; for recent 
  discussions, see refs. <cit.>).
  In the context of PBHs, the non-Gaussianities also change the shape of the 
  probability distribution characterizing the over-densities at the time of 
  their formation, which we have assumed to be a Gaussian [cf. eq. (<ref>)]
  These effects due to the non-Gaussianities are expected to be larger (than the 
  corrections to the power spectrum we have calculated), and they need to be taken 
  into account to arrive at the modified  <cit.>.
  
  
  
  
  
  § CONCLUSIONS
  
  
  In this work, we had considered models involving a single, canonical scalar 
  field that lead to ultra slow roll or punctuated inflation.
  All these models had contained a point of inflection, which seems essential 
  to achieve the epoch of ultra slow roll required to enhance scalar power on 
  small scales.
  We had also examined the extent of PBHs formed and the secondary GWs generated 
  in these models and had compared them with the constraints on the corresponding
  observables  and .
  These models require a considerable extent of fine tuning in order to 
  lead to the desirable duration of inflation (of say, 60–70 e-folds), 
  be consistent with the constraints from the CMB on large scales, and 
  simultaneously exhibit higher scalar power on small scales.
  
  
  
  In order to explore the possibilities in single field models further, 
  we had also considered scenarios wherein the functional forms for the 
  first slow roll parameter closely mimic the typical behavior in ultra 
  slow roll and punctuated inflation.
  We had reconstructed the potentials associated with these scenarios,
  evaluated the resulting scalar and tensor power spectra as well as
  the corresponding imprints on ,  and S_h.
  The presence of extra parameters in the choices for ϵ_1(N)
  had allowed us to construct the required scenarios rather easily.
  Interestingly, we had found that the reconstructed potentials too 
  contain a point of inflection as the original models do.
  This lends further credence to the notion that a point of inflection 
  is essential to achieve ultra slow roll or punctuated inflation.
  However, we should add a note of caution that, while we were able to 
  broadly capture the expected shape of the scalar power spectra in the 
  reconstructed scenarios, there were some differences in the tensor 
  power spectra in these scenarios and the original models.
  Moreover, we find that these reconstructed scenarios allow us to easily
  examine the rate of growth of the scalar power from the CMB scales to 
  small scales 
  (for a discussion in this context, 
  see refs. <cit.>).
  While the steepest growth possible in the reconstructed scenario RS1 has 
  -1≃ 4, we find that the growth is non-uniform but faster in 
  RS2 with -1 between 4 and 6 over the relevant range of wave
  numbers (for details, see appendix <ref>).
  Further, though we have been able to reconstruct the potentials numerically
  in the scenarios RS1 and RS2, it would be worthwhile to arrive at analytical
  forms of these potentials <cit.>.
  
  
  
  We had also computed the scalar bispectrum and the associated non-Gaussianity
  parameter  is these models and scenarios.
  We had found that the parameter  is strongly scale dependent in all
  the cases.  
  Also, the non-Gaussianities had turned out to be fairly large (with, say, 
  >10 over a range of wave numbers) in the equilateral limit.
  Moreover, we had found that the consistency condition governing the 
  non-Gaussianity parameter is always satisfied, despite the period of
  sharp departure from slow roll, implying that the non-Gaussianity 
  parameter in the squeezed limit is at most of order unity around the 
  domain where the scalar power spectra exhibit their peak.
  Due to this reason, we had found that the non-Gaussian corrections to 
  power spectra were negligible leading to insignificant modifications to
  the observables ,  and S_h on small scales.
  However, we should point out that the effects of non-Gaussianities on 
   and  have been included in a simple fashion and a more
  detailed approach seems required to account for the complicated scale 
  dependence of  <cit.>.
  It has recently been argued that, in the squeezed limit
  of the bispectrum, the part satisfying the consistency relation should 
  be subtracted away as it cannot be observed (in this context, see 
  refs. <cit.>; however also
  see Ref. <cit.>).
  If this is indeed so, since the scalar bispectrum satisfies the consistency 
  condition in the squeezed limit in the models and scenarios we have examined, 
  the cubic order non-Gaussian corrections to the power spectrum would then 
  identically vanish.
  
  
  
  Moreover, we had calculated the secondary tensor bispectrum 
  generated in the different inflationary models of interest during 
  the radiation dominated epoch.
  Interestingly, we had found that the shape function characterizing the 
  tensor bispectrum has an amplitude of about 10–250 at small wave 
  numbers in all the models and scenarios of interest.
  It seems important to evolve the shape function until today and examine 
  the possibility of observing its imprints in ongoing efforts such as
  PTA <cit.> and forthcoming missions such as 
  LISA <cit.>.
  We are currently investigating these issues
  in a variety of single and two field models of inflation <cit.>.
  
  
  
  
  
  The authors wish to thank Dhiraj Hazra, Rajeev Jain and Subodh Patil for
  discussions and detailed comments on the manuscript. 
  HVR and LS also wish to thank Arindam Chatterjee, Arul Lakshminarayan and 
  Jérôme Martin for related discussions.
  HVR and PS would like to thank the Indian Institute of Technology Madras (IIT
  Madras), Chennai, India, for support through the Half-Time Research
  Assistantship and the Institute Postdoctoral Fellowship, respectively.
  The authors wish to acknowledge use of the cluster computing facilities at 
  IIT Madras, where some of the numerical computations were carried out.
  LS also wishes to acknowledge support from the Science and Engineering Research  
  Board, Department of Science and Technology, Government of India, through the 
  Core Research Grant CRG/2018/002200.
  
  
  
  
  
  
  
  § THE DICHOTOMY OF ULTRA SLOW ROLL AND PUNCTUATED 
  INFLATION
  
  
  With the help of an example, in this appendix, we shall illustrate that a 
  given inflationary potential can permit ultra slow roll as well as punctuated 
  inflation for different sets of parameters.
  The potential that we shall consider, when expressed in terms of the 
  quantity x = ϕ/v that we had introduced in the context of USR1, is
  given by <cit.>
  
      V(ϕ) = V_0 α x^2-β x^4+γ x^6/(1+δ 
      x^2)^2.
  
  In figure <ref>, we have plotted the evolution of the 
  first slow roll parameter ϵ_1 in the above potential for the 
  following two sets of parameters:
  V_0/^4=1.3253×10^-9, γ=1, δ= 1.5092 and 
  (v/,α,β)=(4.3411,8.522× 10^-2,0.469)
  and (10,8.53× 10^-2, 0.458).
  
  
  
  We obtain about 75 e-folds of inflation in these cases for
  ϕ_i=17.245 and ϕ_i=13.4.
  It is clear from the figure that, while the first set of parameters lead to 
  punctuated inflation, the second set does not permit an interruption of 
  inflation until the very end.
  This example illustrates the point that a potential itself cannot be classified
  as an ultra slow roll or a punctuated inflationary model.
  
  
  
  
  
  § THE FUNCTIONAL FORMS OF THE POLARIZATION FACTORS
  
  
  Recall that, e^λ(,)=e_ij^λ() p^i p^j.
  For our choice of (,,) and (_1,_2,_3) [cf.
  eqs. (<ref>) and (<ref>)], we find that 
  e^λ(,) can be evaluated to be
  
  
  
      e^+(,_1) 
         =   14 √(2) ł(3 p_1x^2 + p_1y^2 
      - 2 √(3) p_1x p_1y - 4 p_1z^2)̊,
      
      e^+(,_2) 
         =   14 √(2) (3 p_1x^2 + 3 k^2 + p_1y^2
      + 2 √(3) p_1x p_1y- 6 k p_1x
      - 2 √(3) k p_1y - 4 p_1z^2),
      
      e^+(,_3) 
         =   1√(2) ł(p_1y^2-p_1z^2)̊,
      
      e^×(,_1) 
         =    -1√(2) ł(√(3) p_1x-p_1y)̊ p_1z,
      
      e^×(,_2) 
         =   1√(2) ł[√(3) (p_1x - k) + p_1y]̊ p_1z,
      
      e^×(,_3)    =    -√(2) p_1y p_1z.
  
  
  
  
  
  
  
  § INTEGRALS DETERMINING THE SCALAR BISPECTRUM
  
  
  The quantities __C(,,) appearing in the 
  expression (<ref>) for the scalar bispectrum represent 
  six integrals that involve the scale factor, the slow roll 
  parameters, the modes f_k and their time derivatives f_k'.
  They correspond to the six bulk terms appearing in the cubic 
  order action governing the curvature perturbation, and they are 
  described by the following 
  expressions <cit.>: 
  
  
      _1(,,)
         =    2 i ∫_^η̣  a^2 ϵ_1^2 ł(f_k_1^∗ f_k_2'^∗ 
      f_k_3'^∗ + two permutations)̊,
      _2(,,)
         =   -2 i ł(· + two permutations)̊ ∫_^η̣  a^2 ϵ_1^2  f_k_1^∗ f_k_2^∗ 
      f_k_3^∗,
      _3(,,)
         =   -2 i ∫_^η̣  a^2 ϵ_1^2 ł(·k_2^2 
      f_k_1^∗ f_k_2'^∗  f_k_3'^∗
      + five permutations)̊,
      _4(,,)
         =    i ∫_^η̣  a^2 ϵ_1 ϵ_2' ł(f_k_1^∗ f_k_2^∗ f_k_3'^∗
      +two permutations)̊,
      _5(,,)
         =   i/2 ∫_^η̣  
      a^2 ϵ_1^3 ł(·k_2^2 
      f_k_1^∗ f_k_2'^∗  f_k_3'^∗ 
      + five permutations)̊,
      _6(,,) 
         =   i/2 ∫_^η̣  a^2 ϵ_1^3 ł(k_1^2 ł(·)̊k_2^2 k_3^2  
      f_k_1^∗  f_k_2'^∗  f_k_3'^∗ 
      + two permutations)̊.
  
  
  These integrals are to be evaluated from a sufficiently early time, say, , 
  when all the modes are well inside the Hubble radius, until suitably late times,
  which can be conveniently chosen to be a time close to the end of 
  inflation, say, .
  
  
  
  
  
  § A CLOSER EXAMINATION OF THE CONSISTENCY RELATION
  
  
  We had pointed out that, in the squeezed limit, i.e. when k_2≃ k_3
  =k and k_1→ 0, the non-Gaussianity parameter  is expected to satisfy 
  the consistency condition ^_CR(k)=(5/12)  [(k)-1], 
  where (k)-1=ln (k)/ln  k is the scalar 
  spectral index.
  In the results presented earlier (in figures <ref> and
  <ref>), we had worked with k_1=10^-3  k to arrive
  at  in the squeezed limit.
  While we find that the consistency condition is satisfied to better than
  5% over a wide range of scales, we notice that there is some departure 
  around wave numbers corresponding to the peak in the scalar power spectrum.
  To investigate this point more closely, in figure <ref>, we have 
  plotted the numerical results around the peak in the scalar power spectrum
  for the original choice of k_1 as well as for k_1 = 10^-1 k and 
  k_1 = 10^-5 k in the case of the model PI3.
  
  
  
  We have considered the case of k_1=10^-1 k since we find that roughly a 
  decade of modes exit the Hubble radius during the ultra slow roll phase.
  Evidently, such a value of k_1 would be insufficient for it to be considered
  a squeezed mode.
  We find that the value of  remains of order unity even when we confine 
  to modes which leave the Hubble radius during the period of ultra slow roll.
  Also, as one would expect, we find that the consistency condition is satisfied 
  better and better as we work with a smaller value of k_1.
  We should clarify that adequate care needs to be taken while evaluating the 
  integrals involved in the calculation of the bispectrum during the ultra 
  slow roll regime. 
  Since there occur rapid changes in the slow roll parameters during this epoch,
  we should regulate the integrals with an appropriate choice for the cut-off 
  parameter κ, especially for the dominant contribution G_4(,,) 
  [cf. appendix <ref>].
  With an appropriate cut-off and with smaller values for the squeezed mode k_1, 
  we find that the match between  and ^_CR indeed improves.
  Nevertheless, even with a smaller of choice of k_1, we still notice some 
  difference near the peak in the power spectrum.
  We feel that this is an artefact and we believe that the difference can be 
  overcome with a further smaller value for k_1.
  However, working with a very small k_1 poses certain numerical challenges, 
  and we will leave it for future investigation.
  We should mention that this an independent issue and stress that it does not 
  affect our main conclusions related to PBHs and GWs. 
  
  
  
  
  
  § ASYMPTOTIC BEHAVIOR OF THE CURVATURE PERTURBATIONS
  
  
  As we mentioned, it has been shown that an indefinite ultra slow 
  roll regime of inflation leads to the violation of the consistency 
  condition <cit.>.
  Since all the models of our interest contain an ultra slow roll 
  phase, one may wonder if a violation of the consistency condition 
  would occur in these cases.
  As we have seen, the consistency condition is satisfied in all the
  cases we have considered.
  This is primarily due to the fact that the ultra slow roll phase lasts 
  only for a finite duration in our models, permitting the eventual 
  freezing of the amplitude of the curvature perturbations.
  
  
  
  In this appendix, we shall illustrate this point with the aid of a 
  truncated version of the scenario RS1.
  We shall consider the following two functional forms for ϵ_1(N):
  
      ϵ_1^III(N) 
         =   ł[ϵ_1a ł(1+ϵ_2a N)̊]̊ ł[1 - tanhł(N - N_1Δ N_1)̊]̊ ,
      ϵ_1^IV(N) 
         =   ł[ϵ_1a ł(1+ϵ_2a N)̊]̊ ł[1 - tanhł(N - N_1Δ N_1)̊]̊ + ϵ_1b.
  
  Evidently, while the first choice lead to an indefinite period of ultra slow 
  roll beyond the e-fold N_1, the second choice restores slow roll 
  when ϵ_1(N) attains the value of ϵ_1b.
  In figure <ref>, we have plotted the behavior of these slow roll
  parameters as well as the evolution of the curvature perturbation for three
  modes which leave the Hubble radius just prior to and after the onset of 
  the ultra slow roll phase.
  
  
  
  We have worked with the following values for parameters involved in plotting
  the figure: ϵ_1a = 10^-4, ϵ_2a = 0.05, N_1=42, 
  Δ N_1 = 0.5 and ϵ_1b = 10^-10.
  It should be clear that, while the amplitude of the curvature perturbation 
  grows indefinitely when the ultra slow roll continues, the amplitude 
  freezes when slow roll inflation is restored.
  
  
  
  
  § THE STEEPEST GROWTH OF THE SCALAR POWER SPECTRUM
  
  
  In models of ultra slow roll and punctuated inflation, we have seen that 
  the scalar power grows rapidly from its COBE normalized values on the CMB 
  scales to higher values at smaller scales over wave numbers that leave 
  the Hubble radius during the transition from slow roll to ultra slow roll. 
  An interesting issue that is worth understanding is the steepest such growth 
  that is possible in models of inflation driven by a single, canonical scalar 
  field.
  It has been argued that the fastest growth will
  have -1≃ 4 over this range of wave numbers (in this 
  context, see ref. <cit.>; also 
  see ref. <cit.>).
  We find that the reconstructed scenarios RS1 and RS2 easily permit us to 
  examine this issue.  
  Recall that, in these scenarios, the parameter Δ N_1 determines the 
  rapidity of the transition from the slow roll to the ultra slow roll regime
  [cf. eqs. (<ref>)].
  We find that it is this parameter that dictates the steepness of the growth
  in the corresponding scalar power spectra, with smaller Δ N_1 producing 
  a faster rise.
  We have examined the rate of growth in the cases of RS1 and RS2 by varying 
  Δ N_1 over a certain range, while keeping the other parameters fixed.
  In figure <ref>, we have illustrated the spectra for four values 
  of Δ N_1 which are relatively smaller than those we had used for the 
  reconstructions discussed earlier.
  
  
  
  It should be clear from the figure that, in the case of RS1, the rise is 
  fairly steady as the value of Δ N_1 is made smaller, with -1 
  ≃ 4 over the growing regime. 
  In the case of RS2, we find that -1 varies between 4 and 6 over 
  the growing regime and therefore corresponds to a steeper but non-uniform 
  growth of the spectra.
  
  
  apsrev4-2
  
  
  