# Introducing the Voxel Interactive Contour Tool for Online Radiation
Intensity Analytics (VICTORIA)

Elise Badun
University of Victoria
Victoria, BC, Canada
<EMAIL_ADDRESS>Frédéric Tessier
National Research Council Canada
Ottawa, ON, Canada
<EMAIL_ADDRESS>Reid Townson
National Research Council Canada
Ottawa, ON, Canada
<EMAIL_ADDRESS>Ernesto Mainegra-Hing
National Research Council Canada
Ottawa, ON, Canada
<EMAIL_ADDRESS>Margaret-Anne Storey
University of Victoria
Victoria, BC, Canada
<EMAIL_ADDRESS>Magdalena Bazalova-Carter
University of Victoria
Victoria, BC, Canada
<EMAIL_ADDRESS>
###### Abstract

In this paper, the Voxel Interactive Contour Tool for Online Radiation
Intensity Analytics (VICTORIA) web viewer is presented as a solution to the
issue of inaccessible, expensive, and insecure Digital Imaging and
Communications in Medicine (DICOM) and Monte Carlo (MC) dose viewers. VICTORIA
is a tool for the visualization and comparison of dose distributions with
underlying anatomy that is accessible, free, designed with security in mind,
and available on the web111http://web.uvic.ca/~bazalova/VICTORIA/. The code is
open-source and under the GNU General Public License, and it can be found at
the EGSnrc Github repo in the DICOM-viewer
branch222https://github.com/ellieb/EGSnrc/tree/DICOM-
viewer/HEN_HOUSE/gui/dose-viewer. The viewer should be useful for researchers
using EGSnrc file types .egsphant and .3ddose from the EGSnrc toolkit,
patients who want a simple tool to view files, or researchers using DICOM
Computed Tomography (CT) and DICOM Radiotherapy (RT) Dose files in low- or
middle-income countries without access to treatment planning systems.

##  1 Introduction

Many popular software tools exist to visualize computed tomography (CT)
images. Programs like Osirix, 3D Slicer, ImageJ, Dicompyler, and Horos are
used to analyze and compare medical images. As important as these tools are,
almost none of them (with the exception of ImageJ) are free and browser-based.
If using software-specific formats such as .egsphant and .3ddose files, online
visualization programs are challenging to find. The goal of VICTORIA is to be
an easy tool to view and compare different dose and density files in a
standard web browser, requiring no additional software download and
installation. VICTORIA was initially developed as a web tool that could easily
view and compare phantom and dose files from the EGSnrc toolkit, an Electron
Gamma Shower software package by the National Research Council of Canada
(NRC). The EGSnrc Monte Carlo (MC) toolkit is a well-known medical physics
package that simulates the transport of ionizing radiation through matter [1].
It can model the transport of electrons, positrons, and photons with energies
ranging from 1 keV to 10 GeV, in user-defined geometries. The EGSnrc package
uses custom file types .egsphant and .3ddose, which are formatted plain text
files containing phantom and dose information.

The toolkit comes with legacy tools to visualize these file types however,
they are outdated and can be hard to install on some operating systems. For
example, a user who wants to run the EGSnrc dosxyz_show program on a Windows
computer would have to install a virtual machine to run the program on Linux.
These types of restrictions prompted the collaboration between the NRC and the
University of Victoria to build a web-based .egsphant and .3ddose viewer. Work
began on the project under the supervision of Magdalena Bazalova-Carter in
conjunction with Frédéric Tessier and Reid Townson and Ernesto Mainegra-Hing
from the NRC. Later, in collaboration with Margaret-Anne Storey, the ability
to load and view DICOM files was added, as requested by users of the site.

##  2 Materials and Methods

VICTORIA is a web-based viewer coded in JavaScript, HTML, and CSS. Most of the
data visualization capabilities come from D3, a JavaScript library made for
producing dynamic data-based displays by binding input to the Document Object
Model (DOM) and applying transformations [2]. D3 provides many user interface
functions in VICTORIA, particularly data plotting, axes support, click events,
zooming, and other changes to the DOM. The dicomParser library from
Cornerstone is used to read in and parse DICOM files [3]. The other libraries
in use are d3-legend, which creates colour-based legends, as well as filesaver
and dom-to-image, which assists in saving the volume viewer to PNG [4, 5, 6].
The code for this project can be found in the EGSnrc GitHub repository
333https://github.com/ellieb/EGSnrc/tree/DICOM-viewer/HEN_HOUSE/gui/dose-
viewer.

Due to the sensitive nature of medical imaging documents, all processing is
done by the browser. This means that any data loaded into the viewer remains
in the client context, so no files are uploaded over a network to a database,
improving file security. To improve site performance, Web Workers are used to
cache density slices. Web Workers allow for code execution in the background
without interfering with the user interface. Since VICTORIA is browser-based,
the only limitations are the browser and browser version a user has. The site
has been tested to work on Chrome versions 85+, Firefox versions 77+, and MS
Edge versions 89+, as of May 10, 2021.

##  3 Results

###  3.1 Supported file formats

The viewer accepts file types of .egsphant, .3ddose, DICOM CT Image, DICOM
Radiotherapy (RT) Dose, and DICOM RT Structure Set [7]. The files can be
loaded with drag and drop, or by opening a file manager to navigate through
the user directory. Alternatively, there are test files loaded on the site
that consist of a .egsphant density file, an RT Dose file, and a structure set
file of a pediatric patient, for users to test out the website. A timed upload
of a 21.7 MB .egsphant file with 161$\times$113$\times$85 voxels took 1.4 s to
load, while a batch of 171 DICOM CT files with a total size of 90.4 MB and
512$\times$512$\times$171 voxels took 13.2 s to load, using Google Chrome
version 90 on OS X El Capitan with a 2.7 GHz processor speed and 2 cores.

###  3.2 User interface

Figure 1: A volume viewer with a CT Image file and RT Dose file loaded and
displayed. The window and level of the density file has been adjusted. The
dose profile along the red cross-hairs has been plotted. The ROI contours are
shown and the corresponding DVH is plotted beneath the dose profiles. The dose
has been normalized to 4500 cGy at 90%, and this is reflected in the DVH
x-axis. A few ROI contours have been toggled off.

The user interface of VICTORIA needs to display multiple 3D volumes on screen,
in each of the three axes. This was accomplished by dividing the website into
panes called volume viewers, where each volume viewer can display a density
file (either .egsphant or CT Image), a dose file (.3ddose, RT Dose, or a dose
difference), and region of interest (ROI) contours in three panels for the
sagittal, coronal, and axial views. The density and dose can be selected using
drop down menus at the top of each volume viewer, and the ROI contours are
matched by the study unique identifier (UID) or position.

###  3.3 Navigation

The volumes need to be easy and quick to navigate, without noticeable lag for
the best user experience. Each panel can traverse through a view of the volume
by moving though its slices. Since the navigation of the slices is based on
position rather than voxel number, files that have incompatible voxel
dimensions but the same position can still be viewed together. To increase the
speed of navigating through a volume, density caching with Web Workers was
implemented. Using Google Chrome, navigating each slice for a 21.7 MB
.egsphant file took 9.2 ms before caching, and 3.1 ms after. For a batch of
171 DICOM CT files with a total size of 90.4 MB, navigating a slice took 45.5
ms before and 14.07 ms after caching, respectively.

###  3.4 Data display

The ability to view the difference between two dose files is a useful feature.
This can allow users to compare two different patient plans generated by
either a commercial treatment planning system (TPS) or by MC codes (i.e.
EGSnrc or TOPAS), to compare doses calculated by a MC code against a
commercial TPS, to compare dose distributions calculated by two difference MC
codes, or to compare MC dose distributions calculated with different transport
parameters. The dose difference is calculated by dividing each dose matrix by
its maximum dose value, then subtracting the second dose matrix from the
first, using tri-linear interpolation if the voxel dimensions are different.
The panels are navigable by using sliders to move through slice by slice, or
by clicking on a point to update the other two panels to the same location.
The panels also support panning and zooming.

The density files are displayed as an image, whereas the dose files are
displayed as isodose contours. The dose contours are calculated using the D3
contours function, which returns an array of polygons where the input values
are greater than the given thresholds, without any smoothing or interpolation
[8]. The default dose contours are from 10% to 100% in 10% increments. An
input box at the bottom of the dose contour legend gives the user the option
to input custom dose contours that are not already displayed (e.g., 5%). A
toggle function is also available to show/hide individual dose contours by
clicking the corresponding legend cell. The option to plot the dose profiles
is also available when a dose file is loaded. The dose profiles contain the
dose data along the cross-hairs where the user clicks in the z, x, and y
dimensions.

Similar functionality to the dose contours is available for the ROI contours.
The structure sets are matched to the volume by the existing density and dose
files, rather than the user selecting which structure set to display. They are
matched if they have the same study instance UID or if the extent of the
density/dose overlap with the ROI contours. The user can toggle on and off
which ROI contours are visible in the volume viewer. If there is a dose file
loaded, there is also an option to plot a dose volume histogram (DVH). This
will go through each ROI contour in the structure set, determine the dose
enclosed by each, and calculate the cumulative dose delivered in each portion
of the ROI. Again, toggling the ROI legend will toggle the data shown in the
DVH plot.

There are various options the user can make with the check-boxes at the top of
each volume viewer. The user can normalize the dose files from cGy to percent,
or vice versa. This changes the scale for the DVH. The normalization from
percent to cGy is mandatory if the user would like to plot a DVH while using a
dose file with relative doses. The voxel information of the position of the
panel click can also be displayed. This will show the voxel number, absolute
position, density, dose, and material of the position where the user clicks.
Another checkbox can plot the dose profile along the cross-hairs of where a
user clicks if a dose volume is currently displayed in the viewer. If a set of
ROI contours has been uploaded, a checkbox will be available to show or hide
them. If ROI contours and dose data are plotted, an option to plot a DVH is
also available.

The sliders on the upper right of each volume viewer control different display
options. The density sliders can be used to adjust the window and level of the
density plots. The max dose slider adjusts the maximum dose for plotting
purposes (i.e. if the max dose is 100 cGy, and the max dose slider is set to
80%, the plot will now be adjusted so the dose at 100% is 80 cGy, adjusting
the dose contours levels). If a dose comparison is plotted, the dose
comparison normalization slider will behave similarly, adjusting the
normalization factor of the second dose volume that is being compared, from
0.5 to 1.5.

###  3.5 Export formats

VICTORIA offers a range of data export formats for plot data and
visualizations. Exporting plot data to CSV allows for analysis in spreadsheet
software. If the dose profiles are displayed, the data can be downloaded by
selecting the ’Export dose profiles to CSV’ button. The positions and dose
values of the three dose profiles are then downloaded as a CSV file. The
’Export DVH to CSV’ button works similarly, exporting the percent of the ROI
volume with the dose value for each of the ROIs. To export the visualization
as an image, the ’Export visualization to PNG’ button can be used to download
a screenshot of the current viewer.

##  4 Discussion

VICTORIA was evaluated against some existing free web-based DICOM viewers
based on format support and performance. Viewers were only tested if they
required no downloads, required no account creation, did all processing
client-side, and were free to use. The Imaios DICOM Viewer was evaluated [9].
It was very fast and easy to use but did not support the display of DICOM RT
Dose nor DICOM RS Structure Set files. Similarly, both the U Dicom Viewer and
the Web Dicom Viewer were evaluated as quick and simple DICOM viewing tools
that can view remote files via URL, however, both viewers only had support for
DICOM CT Image files [10, 11]. The Rocket Viewer was another medical image
viewing website, with support for a large variety of file types however, it
did not support the DICOM RT Dose and RT Structure Set file types, and slices
had to be clicked through to navigate [12]. The majority of the online DICOM
viewers were fast and easy to use, and many had features such as the ability
to draw on lines, make measurements, and import remote files. However, none of
the tested viewers could upload all of DICOM CT Image, RT Dose, and RT
Structure Set files. Additionally, websites that could upload both CT Image
and RT Dose files could not overlay them. Furthermore, they did not have
features such as dose profiles nor DVH plotting. Currently, VICTORIA is the
only online viewer able to read and display the .egsphant and .3ddose file
types.

The VICTORIA viewer is in the early stages of development and evaluation and
could be further enhanced. There are many data manipulation features and user
interface design improvements to turn it into a robust viewer. Inspired by
other viewers, a desirable feature would be the ability to import remote files
hosted on DropBox or GoogleDrive. A library to convert a list of 2D polygons
to a 3D object should be included to build smoother ROI contours. Another
important feature would be speeding up the upload of large files. Some smaller
improvements to consider adding to VICTORIA would include a drop down menu to
change the colour mapping for the dose contours, a dose contour opacity
slider, the ability to close a volume viewer, allowing folder uploads, and an
option to plot the density in the dose profile plots. In you are interested in
contributing to this project, please head over to our github page
444https://github.com/ellieb/EGSnrc/tree/DICOM-viewer/HEN_HOUSE/gui/dose-
viewer and submit a pull request. Any questions about the project can be
directed to<EMAIL_ADDRESS>
##  5 Conclusion

VICTORIA is a convenient web-based Monte Carlo and DICOM dose visualization
and comparison tool. It is unique in the fact that it is a web-based tool that
supports .3ddose and .egsphant file types. The viewer can be accessed at
http://web.uvic.ca/~bazalova/VICTORIA/ [13].

## References

  * [1] I Kawrakow, DWO Rogers, E Mainegra-Hing, F Tessier, RW Townson, and BRB Walters.  EGSnrc toolkit for Monte Carlo simulation of ionizing radiation transport.  doi:10.4224/40001303 [release v2021], 2000.
  * [2] Mike Bostock.  D3 data-driven documents.  https://d3js.org/, 2020.
  * [3] CornerstoneJS.  dicomparser.  https://github.com/cornerstonejs/dicomParser, 2020.
  * [4] Susie Lu.  d3-legend.  https://github.com/susielu/d3-legend, 2019.
  * [5] Eli Grey.  Filesaver.js.  https://github.com/eligrey/FileSaver.js/, 2020.
  * [6] Anatolii Saienko.  dom-to-image.  https://github.com/tsayen/dom-to-image, 2017.
  * [7] NEMA PS3 / ISO 12052.  Digital Imaging and Communications in Medicine (DICOM) Standard.  http://www.dicomstandard.org/.
  * [8] D3.  d3-contour.  https://github.com/d3/d3-contour, 2021.
  * [9] Imaios Dicom Viewer.  https://www.imaios.com/en/Imaios-Dicom-Viewer.
  * [10] webnamics.  U Dicom Viewer 0.8.1.  https://webnamics.github.io/u-dicom-viewer/, 2019.
  * [11] ivmartel.  dwv-jqui 0.5.0.  https://ivmartel.github.io/dwv-jqui/demo/stable/, 2021.
  * [12] BCN MedTech.  Rocket Viewer 1.0.  http://rkt-viewer.surge.sh/, 2018.
  * [13] VICTORIA.  http://web.uvic.ca/~bazalova/VICTORIA/, 2021.


# Non-Gaussianity consistency relations
and their consequences for the peaks

Mohammad Hossein Namjoo  Bahar Nikbakht<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Strong deviations from scale invariance and the appearance of high peaks in
the primordial power spectrum have been extensively studied for generating
primordial black holes (PBHs) or gravitational waves (GWs). It is also well-
known that the effect of non-linearities can be significant in both phenomena.
In this paper, we advocate the existence of a general single-field consistency
relation that relates the amplitude of non-Gaussianity in the squeezed limit
$f_{\text{NL}}$ to the power spectrum and remains valid when almost all other
consistency relations are violated. In particular, it is suitable for studying
scenarios where scale invariance is strongly violated. We discuss the general
and model-independent consequences of the consistency relation on the behavior
of $f_{\text{NL}}$ at different scales. Specifically, we study the size, sign
and slope of $f_{\text{NL}}$ at the scales where the power spectrum peaks and
argue that generally the peaks of $f_{\text{NL}}$ and the power spectrum occur
at different scales. As an implication of our results, we argue that non-
linearities can shift or extend the range of scales responsible for the
production of PBHs or GWs, relative to the window as determined by the largest
peak of the power spectrum, and may also open up new windows for both
phenomena.

††institutetext: School of Astronomy, Institute for Research in Fundamental
Sciences (IPM),
Tehran, Iran, P.O. Box 19395-5531

##  1 Introduction

The possibility of large enhancement of power at sub-CMB scales has been
widely studied due to its necessity for the abundant formation of the
primordial black holes (PBHs) Ozsoy:2023ryl ; Hertzberg:2017dkh ; Garcia-
Bellido:2017mdw  and the generation of induced gravitational waves (GWs)
Domenech:2021ztg . Such large enhancements in the power spectrum can occur if
scale invariance is broken during inflation as a result of phase transitions
and temporary phases of non-attractor evolution. It is also well-known that
the PBH formation and the generation of GWs are both sensitive to the non-
linearities that may be represented by higher order correlation functions
Bullock:1996at ; Cai:2018dig ; Young:2015kda ; Atal:2018neu ;
Franciolini:2023pbf . Thus, it is crucial to simultaneously study the power
spectrum and non-linearities at the scales where the enhancements occur. This
paper provides one step toward that direction by paying particular attention
to the scale-dependence of the two-point and three-point correlation
functions. More explicitly, we study the squeezed limit bispectrum which
corresponds to the correlation between the modes that experience large
enhancements and a mode that has a much larger wavelength.

A useful tool that may enable us to simultaneously study the power spectrum
and the squeezed limit bispectrum is the non-Gaussianity consistency
relations. However, the aforementioned reasons for the enhancement of the
power spectrum — i.e., the presence of phase transitions, the break-down of
scale invariance and the non-attractor phases of inflation — are in fact
different reasons for the violation of almost all non-Gaussianity consistency
relations, if the long mode is affected by these features Finelli:2017fml ;
Hinterbichler:2012nm ; Bravo:2017wyw ; Assassi:2012zq ; Cheung:2007sv ;
Kundu:2015xta . An exception is the consistency condition that is advocated in
Ref. Namjoo_2023  which is claimed to hold more generally because it does not
rely on specific symmetries or restrictive assumptions about the Lagrangian
governing inflation. In this paper, we confirm that claim by customizing the
consistency relation of Ref. Namjoo_2023  specifically for the situation of
interest, which is the single-field inflationary scenarios with possibly
multiple phases of non-attractor evolution and sharp phase transitions that
lead to the large enhancement of the power spectrum and the appearance of
several peaks and troughs (Sec. 2).

If the long mode is not affected by the aforementioned features in the
inflationary background — which happens if the scale-invariant-breaking phases
are sufficiently short — the slow-roll consistency relation Maldacena_2003
holds. This relation implies that the size of the squeezed limit bispectrum is
determined by the degree to which the power spectrum deviates from scale
invariance. Thus, there can be a strong correlation between the long mode and
the short mode due to the strong violation of scale invariance. We will show
that our general consistency relation reduces to the standard one in the same
limit (Sec. 2.2).

Since our consistency relation is generic, we can make general and model-
independent statements about the behavior of the squeezed limit bispectrum as
a function of the wavelength of the short modes. In particular, we study the
size, sign and slope of the bispectrum at the scales where the power spectrum
peaks (Sec. 3). Moreover, we show that the peaks of the power spectrum and the
squeezed limit bispectrum generally occur at different scales. Crucially, this
implies that the largeness of the power spectrum and the significance of non-
linearities may contribute to the PBH formation separately at different
scales, changing the mass window of PBHs relative to the predictions of the
linear theory. Likewise, we argue that the scales where non-linearities can
have a notable impact on the generation of GWs may differ from the peak scale
of the power spectrum. We confirm the validity of the consistency relation and
its implications using a concrete example (Sec. 4).

In the rest of the introduction, we provide some definitions and clarify our
notation. We work with the curvature perturbation on comoving hypersurfaces,
denoted by $\mathcal{R}$. Generally, its two-point, unequal-time correlation
function may be quantified by

$\left\langle\mathcal{R}_{\mathbf{k}}(\tau)\mathcal{R}_{\mathbf{k^{\prime}}}(\tau^{\prime})\right\rangle=(2\pi)^{3}\delta^{(3)}\left(\mathbf{k}+\mathbf{k^{\prime}}\right)P_{\mathcal{R}_{k}}(\tau,\tau^{\prime}).$
(1)

The observable correlation function is the equal-time at the boundary (i.e.,
at the end of inflation, $\tau\to 0$) which we denote by $P_{\mathcal{R}_{k}}$
(without an explicit time argument). Note that,
$P_{\mathcal{R}_{k}}(\tau,\tau^{\prime})$ is not necessarily positive or even
real, unless if $\tau^{\prime}=\tau$. However, we will need this notation only
for the long-wavelength, super-horizon mode which is classical, resulting in
the reality of that function. It will be useful to express the results in
terms of the dimensionless power spectrum
$\mathcal{P}_{\mathcal{R}_{k}}(\tau,\tau^{\prime})\equiv\frac{k^{3}}{2\pi^{2}}P_{\mathcal{R}_{k}}(\tau,\tau^{\prime})$
and
$\mathcal{P}_{\mathcal{R}_{k}}\equiv\frac{k^{3}}{2\pi^{2}}P_{\mathcal{R}_{k}}$
if $\tau=\tau^{\prime}=0$.

The three-point function at the boundary, the bispectrum, may be quantified by

$\left\langle\mathcal{R}_{\mathbf{k}_{\mathbf{1}}}\mathcal{R}_{\mathbf{k}_{\mathbf{2}}}\mathcal{R}_{\mathbf{k}_{\mathbf{3}}}\right\rangle\equiv(2\pi)^{3}\delta^{(3)}\left(\mathbf{k}_{\mathbf{1}}+\mathbf{k}_{\mathbf{2}}+\mathbf{k}_{\mathbf{3}}\right)B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right).$
(2)

We will be interested in the squeezed limit bispectrum where the momentum of
one mode (which we denote by $k_{\ell}$) is much smaller than the momentum of
the other two (which we denote by $k_{s}$). We quantify the bispectrum in this
limit by

$B_{\mathcal{R}}\left(k_{\ell},k_{s},k_{s}\right)\to(2\pi)^{4}\frac{1}{k_{s}^{3}k_{\ell}^{3}}\mathcal{P}_{\mathcal{R}_{\ell}}\mathcal{P}_{\mathcal{R}_{s}}\,\frac{3}{5}{f_{\text{NL}}}(k_{s})\,,\qquad\text{for
$k_{\ell}\ll k_{s}$}\ ,$ (3)

where we used $\mathcal{R}_{\ell}$ and $\mathcal{R}_{s}$ as shorthands for
$\mathcal{R}_{k_{\ell}}$ and $\mathcal{R}_{k_{s}}$, respectively. Eq. (3) also
defines ${f_{\text{NL}}}$ as the amplitude of non-Gaussianity which in general
depends on the wavenumber $k_{s}$. When a relation applies to both long and
short modes, we denote the corresponding wavenumber by $k$. Occasionally, we
will encounter the ratio of the power spectra of the long mode at different
times. Thus, we find it useful to also define

$\displaystyle
r(\tau)\equiv{\cal{P}}_{{\cal{R}}_{\ell}}(0,\tau)/{\cal{P}}_{{\cal{R}}_{\ell}}.$
(4)

To simplify the language, when confusion is not possible, we refer to the
scales where the power spectrum peaks as the peak scale or the peak. Likewise,
we refer to the scale where there is a trough in the power spectrum as the
trough scale or the trough. We may also be interested in the peaks of other
quantities but we refer to them explicitly to avoid ambiguities. Finally,
throughout this paper, we set $M_{\rm Pl}\equiv(8\pi G)^{-1/2}=1$.

##  2 A general consistency relation

A general consistency relation for single-field inflation is presented in Ref.
Namjoo_2023  where it is shown that it reproduces the results of many single-
field inflationary models including the ones that contain a non-attractor
phase of inflation. Here, we review the consistency condition of Ref.
Namjoo_2023  and customize it to be particularly convenient for the analysis
of the models of interest where some modes experience phases of non-attractor
evolution leading to the strong violation of scale invariance.

The first step in the derivation of the consistency relation is to note that
the correlation between the long mode (which is assumed to be a super-horizon,
classical fluctuation) and the short modes can be written by Creminelli:2004yq

$\displaystyle\langle{\cal{R}}_{\ell}(\tau){\cal{R}}_{s}(\tau){\cal{R}}_{s}(\tau)\rangle\simeq\langle{\cal{R}}_{\ell}(\tau)\,\langle{\cal{R}}_{s}(\tau){\cal{R}}_{s}(\tau)\rangle_{{\cal{R}}_{\ell}}\rangle\
,$ (5)

where $\langle{\cal{R}}_{s}^{2}\rangle_{{\cal{R}}_{\ell}}$ denotes the two-
point correlation function of the short modes under the influence of the long
mode and $\tau$ is the conformal time at which the correlations are computed
which, from now on, we set $\tau\to 0$. According to Eq. (5), the long-short
correlation arises from the short modes perceiving the long mode as a change
in their background evolution. This intuition is consistent with the separate
universe picture but note that we need not to assume that ${\cal{R}}_{s}$ is a
super-horizon mode when it starts correlating with the long mode. Assuming
that the effect of the long mode is small, Eq. (5) can be Taylor expanded.
Since we do not necessarily require ${\cal{R}}_{\ell}$ to be a constant in
time, its effect must be taken into account over the entire history of the
short modes. Thus, we perform a functional Taylor expansion as follows

$\displaystyle\langle{\cal{R}}_{s}^{2}(0)\rangle_{{\cal{R}}_{\ell}}\simeq\langle{\cal{R}}_{s}^{2}(0)\rangle_{0}+\int_{\tau_{\ell}}^{0}{\cal{R}}_{\ell}(\tilde{\tau})\dfrac{\delta}{\delta{\cal{R}}_{\ell}(\tilde{\tau})}\langle{\cal{R}}_{s}^{2}(0)\rangle\,d\tilde{\tau}\,,$
(6)

where higher order terms in the expansion are neglected. $\tau_{\ell}$ is
around the horizon-crossing time of ${\cal{R}}_{\ell}$ and we assumed that all
correlations between the long mode and the short modes occur after
$\tau_{\ell}$. This is justified as long as the sub-horizon short modes around
$\tau_{\ell}$ are in the Bunch-Davies vacuum. Note that, the right hand side
of Eq. (6) being a Taylor expansion, we must set ${\cal{R}}_{\ell}=0$ in
$\delta\langle{\cal{R}}_{s}^{2}\rangle/\delta{\cal{R}}_{\ell}$ after the
functional derivative is taken. Substituting this relation to Eq. (5) results
in

$\displaystyle
B(k_{s},k_{s},k_{\ell})\simeq\int_{\tau_{\ell}}^{0}\langle{\cal{R}}_{\ell}(0){\cal{R}}_{\ell}(\tilde{\tau})\rangle^{\prime}\dfrac{\delta}{\delta{\cal{R}}_{\ell}(\tilde{\tau})}\langle{\cal{R}}_{s}^{2}(0)\rangle^{\prime}\,d\tilde{\tau}\,,$
(7)

where $\langle.\rangle^{\prime}$ denotes $\langle.\rangle$ with the factor
$(2\pi)^{3}$ and the momentum conserving delta function dropped. The notation
introduced in Sec. 1 (in particular, the definition $r(\tau)$ according to Eq.
(4)) allows us to also write

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}=\int_{\tau_{\ell}}^{0}r(\tilde{\tau})\,\dfrac{\delta\ln{\cal{P}}_{{\cal{R}}_{s}}}{\delta{\cal{R}}_{\ell}(\tilde{\tau})}\,d\tilde{\tau}.$
(8)

This is the key relation that is shown in Ref. Namjoo_2023  that can reproduce
the squeezed limit bispectra from many single-field models of inflation. We
refer the readers to Ref. Namjoo_2023  for a summary of the assumptions behind
this relation and for its application to some non-trivial examples where many
other consistency relations are violated. Here we only mention that, besides
the conditions stated above, this relation assumes that the long mode can be
described by a single independent solution. That is why the functional
derivatives are taken with respect to ${\cal{R}}_{\ell}$ and not its conjugate
momentum. In most situations, this assumption is satisfied (even if
${\cal{R}}_{\ell}$ is time dependent) but the extension of Eq. (7) to also
include the effect of the conjugate momentum is straightforward. To utilize
Eq. (8), one needs to inspect how the late time power spectrum of the short
mode ${\cal{P}}_{{\cal{R}}_{s}}$ depends on the history of the evolution and
then study how this history may be affected by a long wavelength mode. In
practice, ${\cal{P}}_{{\cal{R}}_{s}}$ depends on a few instances of time in
its history leading to a significant simplification.

###  2.1 General consistency relation for non-scale-invariant correlations

One can stop with the abstract formula Eq. (8) and apply it to specific
examples. However, to reach some generic conclusions for the scenarios that
are of interest in this paper, we make further — but still generic —
assumptions to simplify Eq. (8). What we have in mind is inflationary
scenarios with transient phases of non-attractor evolution that result in a
large enhancement of the power spectrum at sub-CMB scales. The non-attractor
phases can happen e.g., due to the presence of constant segments in the
potential (the so-called ultra-slow-roll phase) or breaks in the potential or
its derivative. The transitions from a short non-attractor phase to the
subsequent slow-roll phase must be sharp enough so that the effect of the non-
attractor evolution on the power spectrum is not washed out Cai_2018 . Then
one expects to have a large peak followed by several (usually smaller) ones in
the power spectrum.111Obviously, when there are peaks, troughs would also
exist but peaks are of more interest due to their dominant contribution to the
PBH or GW production. Our goal is to understand how the squeezed limit
bispectrum behaves at the scales where the power spectrum deviates
significantly from scale invariance and, specifically, at the peaks.

Assume that ${\cal{P}}_{{\cal{R}}_{s}}$ is an explicit function of several
instances of time and the background initial conditions, i.e.,
${\cal{P}}_{{\cal{R}}_{s}}={\cal{P}}_{{\cal{R}}_{s}}(k_{s},\phi_{0},\tau_{*},\tau_{0},\tau_{1},...,\tau_{n})$,
where $\tau_{0}$ is some initial conformal time222$\tau_{0}$ must be early
enough so that the effect of even earlier history on
${\cal{P}}_{{\cal{R}}_{s}}$ can be neglected but it also needs to be late
enough so that the long mode has already left the horizon (i.e.,
$|\tau_{0}|\lesssim|\tau_{\ell}|$)., $\phi_{0}$ depicts collectively the
dependence on the initial conditions (i.e., the value of the inflaton field
and its conjugate momentum at $\tau_{0}$) and we have also taken into account
a possible dependence on the horizon-crossing time $\tau_{*}=-1/k_{s}$.333We
assume that all other conformal times in the argument of
${\cal{P}}_{{\cal{R}}_{s}}$ are independent of $k_{s}$. We do not make
explicit assumptions about the other instants $\tau_{j}$ for $j=1,...,n$ but
they can appear e.g., when $n$ distinct phase transitions occur, responsible
for the large enhancement of the short mode power spectrum. To derive a more
explicit expression for ${f_{\text{NL}}}$ from Eq. (8), we need to inspect how
the presence of the long mode affects ${\cal{P}}_{{\cal{R}}_{s}}$ due to its
effect on the history of ${\cal{R}}_{s}$. First, note that the long mode
induces an initial perturbation to the inflaton field (and thus to its
conjugate momentum) via
$\Delta\phi_{0}=-(\dot{\phi}_{0}/H){\cal{R}}_{\ell}(\tau_{0})$. Furthermore,
the line element as seen by the short mode in the presence of the long mode is
given by $ds^{2}\simeq-dt^{2}+a^{2}e^{2{\cal{R}}_{\ell}}d{\bf{x}}^{2}$.
Therefore, the effect of ${\cal{R}}_{\ell}$ on different instances of time can
be described by a rescaling of the scale factor,
$a\to\hat{a}=ae^{{\cal{R}}_{\ell}}$. A positive ${\cal{R}}_{\ell}$ leads to a
larger scale factor, which implies that the conformal time of each instant
which ${\cal{P}}_{{\cal{R}}_{s}}$ depends on moves forward. One exception is
the horizon-crossing conformal time which must remain fixed for a fixed
comoving wavenumber $k_{s}$ which implies that a positive ${\cal{R}}_{\ell}$
shifts $\tau_{*}$ backward. Taking all these effects into account results in

$\displaystyle\dfrac{\delta\ln{\cal{P}}_{{\cal{R}}_{s}}}{\delta{\cal{R}}_{\ell}(\tilde{\tau})}=\dfrac{\partial\phi_{0}}{\partial{\cal{R}}_{\ell}}\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial\phi_{0}}\delta(\tilde{\tau}-\tau_{0})+\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln\tau_{*}}\delta(\tilde{\tau}-\tau_{*})-\sum_{j=0}^{n}\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln\tau_{j}}\delta(\tilde{\tau}-\tau_{j})\
,$ (9)

where the first term is schematic and takes into account all effects of
${\cal{R}}_{\ell}$ on the initial conditions. Plugging this into Eq. (8)
results in

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}=r(\tau_{0})\,\dfrac{\partial\phi_{0}}{\partial{\cal{R}}_{\ell}}\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial\phi_{0}}+r(\tau_{*})\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln\tau_{*}}-\sum_{j=0}^{n}r(\tau_{j})\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln\tau_{j}}.$
(10)

To attain some insight from this key result, we note that the characteristic
comoving scales are only meaningful relative to each other. Thus, defining
$x_{k_{s}}=-k_{s}\tau_{n}$, and $x_{j}=\tau_{j}/\tau_{n}$ for $0\leq j<n$, the
power spectrum can be described by
${\cal{P}}_{{\cal{R}}_{s}}={\cal{P}}_{{\cal{R}}_{s}}(\phi_{0},\tau_{*},x_{k_{s}},x_{1},...,x_{n-1})$
where we kept the dependence on $\tau_{*}$ explicit due to its different role
and meaning. This allows us to rewrite Eq. (10) as

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}$ $\displaystyle=$ $\displaystyle
r(\tau_{0})\,\dfrac{\partial\phi_{0}}{\partial{\cal{R}}_{\ell}}\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial\phi_{0}}+\sum_{j=0}^{n-1}\big{(}r(\tau_{n})-r(\tau_{j})\big{)}\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial\ln
x_{j}}$ (11) $\displaystyle+$
$\displaystyle\big{(}r(\tau_{*})-r(\tau_{n})\big{)}\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln\tau_{*}}-r(\tau_{n})(n_{s}-1)\,,$

where $(n_{s}-1)$ is the spectral index that takes into account both the
explicit and implicit dependence on $k_{s}$:

$\displaystyle n_{s}(k_{s})-1\equiv\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln
k_{s}}=x_{k_{s}}\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial
x_{k_{s}}}-\tau_{*}\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial\tau_{*}}.$
(12)

Eq. (10) (or Eq. (11)) is the final form of ${f_{\text{NL}}}$ that, we claim,
reproduces the bispectrum of all examples for which our assumptions are
satisfied (including many models studied in the literature for the purpose of
generating PBHs or GWs). This expression may be further simplified when
further assumptions are made or a more explicit form of
${\cal{P}}_{{\cal{R}}_{s}}$ is considered. In the rest of this section and in
Sec. 3, we discuss the consequences of this generic result and in Sec. 4 we
study a concrete example to examine our consistency relation.

###  2.2 Conserved long mode and the transient non-attractor evolution

First, we consider a situation where the long mode ${\cal{R}}_{\ell}$ has left
the horizon during the slow-roll inflation prior to any non-attractor phase
and is completely frozen when the short modes experience non-attractor
evolution. The attractor behavior is expected to be reached quickly due to the
rapid Hubble expansion, making the non-attractor phases of inflation very
short. Therefore, this situation indeed covers most realistic scenarios that
generate large peaks in the power spectrum. In this case, the initial time
$\tau_{0}$ is in the slow-roll phase well before the onset of the non-
attractor evolution. This implies that the power spectrum is insensitive to
the initial conditions and the first term of Eq. (11) can be neglected.
Moreover, ${\cal{R}}_{\ell}$ is conserved during the entire evolution that
lays down the final form of ${\cal{P}}_{{\cal{R}}_{s}}$. That is, we have
$r(\tau)\equiv{\cal{P}}_{{\cal{R}}_{\ell}}(0,\tau)/{\cal{P}}_{{\cal{R}}_{\ell}}=1$
for all $\tau>\tau_{\ell}$ which significantly simplifies Eq. (11) and reduces
it to

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}=-(n_{s}-1)\ .$ (13)

Remarkably, this is precisely the standard, slow-roll consistency relation
(but holds even if the short modes — but not the long mode — experience non-
attractor evolution) Maldacena_2003 . This result is in fact expected since
this consistency relation only relies on the conservation of the long mode and
does not require restrictions on the evolution of the short modes.

It is well-known that the consistency relation Eq. (13) can be violated if the
long mode is not conserved which can happen if the long mode also experiences
non-attractor phases of inflation Namjoo:2012aa . Our general consistency
relation Eq. (11) holds even in this case since the conservation of the long
mode is not required.

###  2.3  Bispectrum at large and short scales

If the mode of interest ${\cal{R}}_{s}$ experiences non-attractor phases when
its wavelength is larger than the horizon, the resulting changes to the
correlation functions will be scale-invariant. This is the situation where the
methods that rely on the gradient expansion for all modes, such as the $\delta
N$-formalism, are expected to hold. Because of the scale invariance, the third
and the last terms in Eq. (11) vanish and we are left with the dependence on
the initial conditions and the relative time scales of the non-attractor
phases. We will see in Sec. 4 that our consistency relation reproduces the
results from the $\delta N$-formalism at these scales (while our result is
more generic and works for all scales).

Eq. (11) allows us to also make generic statements about ${f_{\text{NL}}}$ at
very short scales $k_{s}\to\infty$. As discussed in Sec. 2, the scenarios of
interest are expected to have sharp transitions. Consequently, at short scales
the power spectrum oscillates rapidly. Because of the strong violation of
scale invariance and the rapid oscillations of the power spectrum, in Eq. (11)
the term proportional to the spectral index dominates and we simply have
${f_{\text{NL}}}\propto(n_{s}-1)$ in this limit.444Note that in realistic
scenarios, the transitions are not infinitely sharp. Consequently, this
argument fails if one goes to extremely short scales where the corresponding
modes perceive the transition as smooth. The proportionality constant is
model-dependent and generally differs from the one in the slow-roll
consistency relation Eq. (13). This result will have important consequences
for the peak scales which will be discussed in Sec. 3. In the explicit example
of Sec. 4.2 we will see that ${f_{\text{NL}}}$ converges quickly to this
result as one goes to short scales and strong deviations from it only occur at
the first few peaks.

##  3 Consequences of the consistency relation for the peaks

Having equipped with the general consistency relation, in this section we make
general statements about the size, sign and slope of the squeezed limit
bispectrum at the peak scales. What we obtain here is of significant
importance concerning PBHs or GWs production.

###  3.1 Size of bispectrum

Concerning the size of ${f_{\text{NL}}}$, the most realistic situation — i.e.,
when the non-attractor phase is short — is indeed the simplest one to analyze
because of the simplicity of the consistency relation Eq. (13). The simple
conclusion is that ${f_{\text{NL}}}$ vanishes at the peaks, where
$d{\cal{P}}_{{\cal{R}}_{s}}/dk\propto(n_{s}-1)=0$.

When Eq. (13) is violated, the more general form of consistency relation Eq.
(11) holds which allows for large values of ${f_{\text{NL}}}$ at the peaks. We
will confirm this in an explicit example in Sec. 4.2. However, we note that
when we go to scales shorter than the first peak — as discussed in Sec. 2.3 —
${f_{\text{NL}}}$ will be dominated by a term proportional to $(n_{s}-1)$.
This implies that scales away from the peaks may have much larger values of
${f_{\text{NL}}}$. This may cause the realization of large fluctuations to
become more probable at non-peak scales which will be explored in the explicit
example of Sec. 4 in Secs. 4.3 and 4.4.

###  3.2 Sign of bispectrum

The sign of the squeezed limit bispectrum is also of interest since the rate
of PBH formation is sensitive to it Young:2015kda . When the slow-roll
consistency relation Eq. (13) holds, ${f_{\text{NL}}}$ vanishes at the peaks
so this question is not relevant. Therefore, we focus on the cases where Eq.
(13) is violated.

The sign of ${f_{\text{NL}}}$ is studied in Ref. firouzjahi2023sign  which
claims that ${f_{\text{NL}}}$ at the peaks is always positive. We will see in
the explicit example of Sec. 4.2 that ${f_{\text{NL}}}$ at the peak is
different from what is obtained in Ref. firouzjahi2023sign  and argue that the
result of Ref. firouzjahi2023sign  can indeed become negative. However, our
consistency relation Eq. (11) is sufficiently powerful to allow us to extract
the assumptions required for the positivity of ${f_{\text{NL}}}$ at the peaks
(besides the one required for the validity of Eq. (11) itself). Recall that
$\tau_{n}$ was an arbitrary reference time that ${\cal{P}}_{{\cal{R}}_{s}}$
depends on. To simplify the analysis, without loss of generality, we assume
that $\tau_{n}$ is chosen so that $r(\tau_{n})\geq r(\tau_{j})$ for all $j\neq
n$. Then it is easy to see that, according to Eq. (11), ${f_{\text{NL}}}$ is
positive at the peaks if (i) the dependence of ${\cal{P}}_{{\cal{R}}_{s}}$ on
the horizon crossing time $\tau_{*}$ can be neglected, (ii) the dependence of
${\cal{P}}_{{\cal{R}}_{s}}$ on the initial conditions can be neglected, (iii)
$r(\tau_{n})$ (hence ${\cal{P}}_{{\cal{R}}_{\ell}}(0,\tau_{n})$) is positive
and (iv) a longer separation between $\tau_{n}$ and $\tau_{j}$ for all $j\neq
n$ leads to a larger power spectrum at the peaks.555Note that a simultaneous
negation of conditions (iii) and (iv) also leads to the same conclusion, i.e.,
the positivity of ${f_{\text{NL}}}$ at the peak. However, this does not happen
— at least — for the simplest non-attractor models of inflation. We do not
claim that when any of the above assumptions is violated ${f_{\text{NL}}}$
will or even can be negative. If there is a more general theorem stating the
positivity of ${f_{\text{NL}}}$ beyond these assumptions, it remains to be
proven. Assumptions (i) and (ii) are satisfied in many scenarios that predict
large peaks in the power spectrum.666See, however, the discussion in Ref.
Namjoo_2023  where it is shown that the example of Ref. Chen:2013aj  violates
condition (ii). Condition (iii) holds if the mode function of
${\cal{R}}_{\ell}$ does not change sign during which it affects the short
modes (which is often the case for non-attractor evolutions with sharp
transitions). (iv) is also likely to be satisfied in most scenarios since
$\tau_{n}$ is likely to be the conformal time at the final phase transition
after which the attractor behavior is reached and then condition (iv)
corresponds to requiring a larger ${\cal{P}}_{{\cal{R}}_{s}}$ at the peak for
a longer non-attractor period.

We conclude that despite the restrictive appearance of the above assumptions,
they are indeed satisfied in many scale-invariance-breaking models of
inflation that lead to large enhancement of the power spectrum, hence
${f_{\text{NL}}}$ is likely to be positive at the peaks. The complex example
that will be studied in Sec. 4 confirms that conclusion in the entire
parameter space. It would be interesting to explore the possibility of
violations of these assumptions and whether that can lead to a negative
${f_{\text{NL}}}$ at the peaks. In particular, the sensitivity to the initial
conditions can happen in models like the one in Ref. Chen:2013aj . Also,
conditions (iii) and (iv) may be violated if a smooth phase transition occurs
or if multiple phases of non-attractor evolution with sufficiently complex
dynamics exist.

###  3.3 Slope of bispectrum

We can also study the slope of ${f_{\text{NL}}}$ around the peak scales using
the consistency relations. In the scenarios where the slow-roll consistency
relation Eq. (13) holds, it is obvious that we have
${f^{\prime}_{\text{NL}}}>0$ at the peaks. When Eq. (13) is violated, from Eq.
(11) the same conclusion can be deduced if conditions (i) and (iii) of Sec.
3.2 are satisfied. Since these assumptions typically hold in models of
interest, we conclude that the positivity of ${f^{\prime}_{\text{NL}}}$ at the
peaks is generic. This result may have notable implications for the PBHs and
GWs. In particular, together with the positivity of ${f_{\text{NL}}}$, it
shows that ${f_{\text{NL}}}$ is larger at shorter-than-the-peak scales. Thus,
if non-linearities are strong enough to shift the window of scales responsible
for the generation of PBHs or GWs, the shift ought to be toward the shorter-
than-the-peak scales. (See Secs. 4.3 and 4.4 for related illustrative
examples.)

We end this section with a few remarks on the troughs. A similar analysis can
be done to investigate the size, sign and slope of the bispectrum at the
trough scales. In particular, one can argue for the non-positivity of
${f_{\text{NL}}}$ and negativity of ${f^{\prime}_{\text{NL}}}$ at the troughs,
under the similar conditions discussed above. Crucially, in different examples
that we have studied, when condition (iv) holds and the power spectrum at the
peak scales increases for a longer period of non-attractor evolution, the
power spectrum at the trough scales decreases. This is the reason for the
negativity of ${f_{\text{NL}}}$. However, it is worth noting that — depending
on the details of the model — the power spectrum at the trough scale may
vanish which is the minimum value that it can take. In this case,
${f_{\text{NL}}}$ also vanishes (rather than being negative) according to the
consistency condition Eq. (11). Thus, in this sense, we are in a similar
situation as where the slow-roll consistency relation Eq. (13) holds, even
though it may be strongly violated.

##  4 A concrete example: Ultra-slow-roll inflation followed by a step

In this section, we verify the validity of the consistency relation Eq. (11)
in a concrete example and make explicit its consequences for the peaks. We
consider a transient ultra-slow-roll (USR) phase sandwiched between two slow-
roll (SR) phases and also consider an infinitely sharp step in the potential
at the end of the USR phase. Such a sharp step is unrealistic but we take it
into account to examine our consistency relation in a more complex but yet
tractable example. The USR phase occurs when there is a flat segment of the
potential. We will have two transitions (from SR to USR and vice versa) that
we assume are both instantaneous and model the potential as follows

$V(\phi)=V_{0}\left[1+\left(\kappa+\sqrt{2\epsilon_{V}}(\phi-\phi_{e})\right)\,\theta(\phi_{e}-\phi)+\sqrt{2\epsilon_{i}}(\phi-\phi_{i})\,\theta(\phi-\phi_{i})\right],$
(14)

where $V_{0}$ is a constant (and the dominant term), $\epsilon_{i}$ and
$\epsilon_{V}$ control the constant slopes of the slow-roll pieces (and
coincide with the attractor slow-roll parameter at each piece) and $\phi_{i}$
and $\phi_{e}$ are the field values at the SR-USR and USR-SR transitions,
respectively. $\kappa$ determines the size of the step after the USR phase; if
$\kappa>0$ ($\kappa<0$) we have an upward (a downward) step. We assume
$|\kappa|\ll 1$. We denote the conformal time at the two transitions by
$\tau_{i}$ and $\tau_{e}$, respectively. We also denote the slow-roll
parameter $\epsilon\equiv-\dot{H}/H^{2}$ at the end of the USR phase by
$\epsilon_{e}$ and right after the step by $\epsilon_{s}$. To better express
our final results we follow Refs. Cai:2022erk ; Cai_2018  and also define

$\displaystyle h\equiv-6\sqrt{\frac{\epsilon_{V}}{\epsilon_{s}}}\ ,\qquad
g\equiv\sqrt{\frac{\epsilon_{s}}{\epsilon_{e}}}.$ (15)

$h$ determines the relaxation time to the final slow-roll phase and $g$
quantifies the abrupt change of the slow-roll parameter due to the step. By
definition, we always have $h<0$; larger values of $|h|$ correspond to the
faster relaxation (after which the super-horizon ${\cal{R}}$ freezes).
Furthermore, we always have $g>0$; $0<g<1$ for an upward and $g>1$ for a
downward step. The step disappears when $g=1$. Since we assume that the step
is instantly taken by inflaton, from the energy conservation, we have
$\epsilon_{s}=\epsilon_{e}-3\kappa$. During the USR phase the slow-roll
parameter decays like $\epsilon\sim a^{-6}$ leading to the rapid growth of the
curvature perturbation like $a^{3}$.

One advantage of this model is that it allows us to also study a situation
where the long mode experiences the non-attractor phase (by extending the USR
phase, i.e., taking the limit $\tau_{i}\to-\infty$). In this case, the
standard consistency relation Eq. (13) is violated but, we claim, the more
general one Eq. (11) holds.

This model has already been discussed in Cai:2022erk  but the bispectrum is
only obtained using the $\delta N$-formalism. We will obtain
${f_{\text{NL}}}(k_{s})$ that is valid for all values of $k_{s}$ including the
peaks where we show that it generically differs from the results of the
$\delta N$-formalism. The model without the step is also extensively discussed
in the literature, see e.g., Chen:2013aj . The full expression for
${f_{\text{NL}}}(k_{s})$ is studied in this simplified situation in
tada2023cancellation  but not for the case of extended USR. Below, we only
outline the existing analysis and refer the readers to e.g., Cai:2022erk  and
tada2023cancellation  for further details but also provide some details in
Appendix A. While we keep the analytic calculations generic throughout this
section, when we present numerical results, we choose $g=1$ which corresponds
to the case without the step. This is because, as mentioned earlier, an
infinitely sharp step is unrealistic and one expects the slope of the step to
be observable to sufficiently short modes. The choice of other parameters will
be reported in the figures.

###  4.1 Transient ultra-slow-roll inflation

We start with analyzing the situation where the USR phase is transient, so we
have two phase transitions to take into account. The mode function of the
curvature perturbation can be expressed generally by

$\mathcal{R}_{k}(\tau)=\frac{H}{\sqrt{4\epsilon(\tau)k^{3}}}\left[\alpha_{k}(1+ik\tau)e^{-ik\tau}+\beta_{k}(1-ik\tau)e^{ik\tau}\right],$
(16)

where $H$ is the Hubble parameter and is assumed to be a constant while the
time dependence of $\epsilon$ is arbitrary. $\alpha_{k}$ and $\beta_{k}$ are
the Bogolyubov coefficients that may be determined at each phase by matching
the mode function to the one from the previous phase at the transition time,
using the continuity of ${\cal{R}}$ and $\epsilon{\cal{R}}^{\prime}$
Namjoo_2012 . The coefficients at the first SR phase are also fixed by the
Bunch-Davies vacuum. We leave the details of the derivation of the Bogolyubov
coefficients to the Appendix. A.2 and only present here the late time power
spectrum which is given by

$\mathcal{P}_{\mathcal{R}}(k)=\dfrac{H^{2}}{8\pi^{2}\epsilon_{V}}|\alpha_{k}^{(3)}+\beta_{k}^{(3)}|^{2}\
,$ (17)

where $\alpha^{(3)}_{k}$ and $\beta^{(3)}_{k}$ are the Bogoliubov coefficients
at the third stage of the evolution (i.e., the final SR phase) which are
presented in Eq. (39) of Appendix. A.2. Note that the dependence of
$\mathcal{P}_{\mathcal{R}}$ on the history of the evolution is via the
explicit appearance of $\tau_{i}$ and $\tau_{e}$ in the Bogoliubov
coefficients as well as the dependence of $h$ and $g$ on
$\epsilon_{e}=\epsilon_{i}(\tau_{i}/\tau_{e})^{6}$. It is clear from Eqs. (38)
and (39) and from the functional form of $\epsilon_{e}$ that the power
spectrum can be entirely written in terms of relative scales, as claimed in
Sec. 2. The left panel of Fig. 1 shows the behavior of the power spectrum.

In the transient USR scenario, the long mode leaves the horizon before the
first phase transition and is completely frozen during the non-attractor
evolution, hence the standard consistency relation Eq. (13) holds. The
validity Eq. (13) is already confirmed in Ref. tada2023cancellation  in a
similar but simpler setup (no step is considered). So we do not show the
explicit comparison that we have made (including the step) between the
consistency relation and the in-in calculations of this model (but we will do
so for the more complex situation of Sec. 4.2). The right panel of Fig. 1
shows the behavior of ${f_{\text{NL}}}$ according to Eq. (13), where we have
also marked the extrema of the power spectrum. Clearly, ${f_{\text{NL}}}$
vanishes at the peaks and troughs.

Figure 1: The power spectrum (left) and the squeezed limit bispectrum (right)
for the transient USR model as functions of $x_{k_{s}}=-k_{s}\tau_{e}$. The
shaded regions indicate the modes that exit the horizon during the USR phase.
On the right panel, we also mark the peak and trough scales of the power
spectrum. For this figure, we set $g=1$, $h=-10$ and $\tau_{i}/\tau_{e}=-10$.

###  4.2 Extended ultra-slow-roll inflation

A more non-trivial exercise for the application of our consistency relation is
the case where the USR phase is extended so that the long mode leaves the
horizon in the USR phase. In this case, we have only one phase transition
(from USR to SR) and the power spectrum can be obtained by taking the limit
$\tau_{i}\to-\infty$ of Eq. (17) which results in

$\displaystyle\mathcal{P}_{\mathcal{R}}=\frac{H^{2}}{32\pi^{2}g^{2}x_{k}^{6}\epsilon_{V}}|\Gamma|^{2}\
,$ (18)

where $x_{k}=-k\tau_{e}$ and

$\displaystyle\Gamma$ $\displaystyle\equiv$
$\displaystyle\left(6-g^{2}(6+h)(1-ix_{k})-2x_{k}(3i+x_{k})\right)\left(\sin(x_{k})-x_{k}\cos(x_{k})\right)$
$\displaystyle+2g^{2}x_{k}^{2}(1-ix_{k})\sin(x_{k}).$

The left panel of Fig. 2 visualizes this power spectrum.

Let us now study the squeezed limit bispectrum by applying our consistency
relation Eqs. (10) or (11). We need to inspect how the power spectrum
${\cal{P}}_{{\cal{R}}_{s}}$ (Eq. (18) with the replacement $k\to k_{s}$)
depends on the history of its evolution. It is clear from Eq. (18) that the
final power spectrum is sensitive to $\tau_{e}$ so the effect of the long mode
on the transition time must be taken into account. Clearly, there is no
dependence on the horizon-crossing time. Moreover, one can argue that the
contribution from the initial conditions can also be neglected: The power
spectrum Eq. (18) is affected by the initial conditions via the slow-roll
parameter $\epsilon=(\dot{\phi}_{0}^{2}/2H^{2})(\frac{\tau}{\tau_{0}})^{6}$
where $\dot{\phi}_{0}$ is the initial inflaton’s velocity and $\tau_{0}$ is
some initial conformal time (i.e., we assume
$|\tau_{\ell}|<|\tau_{0}|\ll|\tau_{e}|$).777As discussed in Sec. 2 the choice
of $\tau_{0}$ is rather arbitrary and we only require it to be early enough so
that an even earlier effect of the long mode can be neglected. Let us stress
that this expression for $\epsilon$ is different from the one for the
transient USR scenario where we had $\epsilon=\epsilon_{i}(\tau/\tau_{i})^{6}$
in the USR phase where $\tau_{i}$ is the conformal time at the onset of the
USR phase and $\epsilon_{i}$ is the attractor slow-roll parameter prior to the
USR phase. In principle, ${\cal{R}}_{\ell}$ can affect both $\dot{\phi}_{0}$
and $\tau_{0}$. The former is negligible since the induced field
$\Delta\phi=-(\dot{\phi}/H){\cal{R}}_{\ell}$ is almost conserved even in the
USR regime so that it cannot contribute significantly to the inflaton’s
velocity. The effect of ${\cal{R}}_{\ell}$ on $\tau_{0}$ can also be neglected
since $r(\tau)$ at $\tau_{0}$ is much smaller than that at later times due to
the rapid growth of the power spectrum.

We are thus left with the only contribution from the impact of
${\cal{R}}_{\ell}$ on $\tau_{e}$. We need to know $r(\tau_{e})$ (which is
${\cal{R}}_{\ell}(\tau)$ at $\tau_{e}$ relative to that at $\tau\to 0$ ) which
can be deduced from the expansion of the mode function Eq. (16) for small $k$.
This yields $r(\tau_{e})=\frac{g^{2}h}{g^{2}h-6}$. Then Eqs. (10) and (11)
reduce, respectively, to

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}(k_{s})=\big{(}\frac{g^{2}h}{6-g^{2}h}\big{)}\dfrac{d\ln{\cal{P}}_{{\cal{R}}_{s}}}{d\ln\tau_{e}}=\big{(}\frac{g^{2}h}{g^{2}h-6}\big{)}\left[\dfrac{\partial\ln{\cal{P}}_{{\cal{R}}_{s}}}{\partial\ln
x_{0}}-(n_{s}-1)\right]\ ,$ (20)

where, for the last equality, we noted that $\tau_{e}$ appears in the
formulation of Eq. (11) via $x_{k_{s}}=-k_{s}\tau_{e}$ and
$x_{0}=\tau_{0}/\tau_{e}$ and we neglected $r(\tau_{0})$ relative to
$r(\tau_{e})$ in the last term of Eq. (11). Note that $g$ and $h$ in Eq. (18)
also depend on $\tau_{e}$ through $\epsilon_{e}\propto\tau_{e}^{6}$ (and
recall that $\epsilon_{s}=\epsilon_{e}-3\kappa$). Remarkably, this result
coincides with the one from explicit in-in calculations. The in-in
calculations are rather tedious but we outline the main steps in Appendix. A
where we also present the explicit form of ${f_{\text{NL}}}$ (Eq. (51)). The
right panel of Fig. 2 shows the behavior of ${f_{\text{NL}}}$ as a function of
$k_{s}$ where we also mark the location of peaks and troughs.

Let us study the analytic form of ${f_{\text{NL}}}$ in the two limits of
interest. We first take the limit $k_{s}\to 0$ which corresponds to the modes
that are super-horizon when the transition to the SR phase occurs. We have

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}(k_{s})\to\dfrac{6h(g^{4}h+6g^{2}-6)}{(g^{2}h-6)^{2}}\qquad\text{for
$k_{s}\to 0$}\,.$ (21)

Note that, as stated in Sec. 2.3, this is the limit where the $\delta
N$-formalism is also valid and we reproduced precisely that result Cai:2022erk
. The same result is also obtained in Ref. firouzjahi2023sign  but it is
claimed that it holds for the peak. Our expression for the peak would be
different since $k_{s}\to 0$ is not a valid limit for the first peak mode
which indeed leaves the horizon after the phase transition (i.e.,
$x_{k_{s}}>1$ for the first peak; see the left panel of Fig. 2). Crucially,
contrary to the claim of Ref. firouzjahi2023sign , it is possible to have a
negative ${f_{\text{NL}}}$ according to Eq. (21) (for an appropriate range of
$g$ when $|h|<6$) whereas our full result Eq. (20) implies that at the peaks
${f_{\text{NL}}}$ always remains positive (see the right panel of Fig.
2).888We note that Ref. firouzjahi2023sign  assumes continuity of
$\dot{\cal{R}}$ at the phase transitions for derivation of a consistency
relation which does not hold in this setup. More generally, and for the
scenario we consider here (which is the same as the one considered in Ref.
firouzjahi2023sign ), $\epsilon\dot{\cal{R}}$ is the continuous quantity at
the transitions. This is because all the assumptions discussed in Sec. 3.2 for
the positivity of ${f_{\text{NL}}}$ at the peaks are satisfied in this
example.

In the opposite limit (i.e., $k_{s}\to\infty$) , Eq. (20) (or, more
explicitly, Eq. (51)) reduces to

$\displaystyle\dfrac{12}{5}{f_{\text{NL}}}(k_{s})\to\big{(}\frac{g^{2}h}{6-g^{2}h}\big{)}(n_{s}-1)\
,\qquad\text{for $k_{s}\to\infty$}\ ,$ (22)

where in the same limit we have

$\displaystyle(n_{s}-1)\to\begin{cases}\frac{2g^{2}\left(g^{4}-1\right)hx_{k_{s}}\sin(2x_{k_{s}})}{\left(g^{2}h-6\right)\big{(}\left(g^{4}-1\right)\cos(2x_{k_{s}})-g^{4}-1\big{)}}\,,\quad\text{for
$g\neq 1$}\,,\\\ \\\ \frac{h^{2}\cos(2x_{k_{s}})}{h-6}\,,\hskip
93.89418pt\text{for $g=1$}\,.\end{cases}$ (23)

Despite the significant difference in the behavior of $(n_{s}-1)$ for $g=1$
(no step) and $g\neq 1$, the relation ${f_{\text{NL}}}\propto(n_{s}-1)$
remains universal, consistent with our generic argument in Sec. 2.3.

From the short wavelength expansion of Eq. (20) and keeping next-to-leading
order terms, one can also deduce that
${f_{\text{NL}}}(k_{\text{peaks}})\simeq-{f_{\text{NL}}}(k_{\text{troughs}})\to
5h|g^{2}-1|/2(g^{2}h-6)$ which is a constant (and zero for $g=1$) and the same
for all peaks and troughs in this limit. We see once again that
${f_{\text{NL}}}$ is non-negative at the peaks and non-positive at the
troughs. Finally, note that it is clear from Eq. (20) that
${f^{\prime}_{\text{NL}}}>0$ at all peaks and ${f^{\prime}_{\text{NL}}}<0$ at
all troughs as claimed in Sec. 3.3.

We have accomplished our main objective in this paper which is to derive and
confirm a generalized consistency relation, applicable to the single-field but
not necessarily scale-invariant models of inflation and studying its
consequences for the peaks. In the next two subsections, we make some side
remarks about the implications of the consistency relation on the predictions
of the PBH and GW production and leave a thorough analysis for future works.

Figure 2: The power spectrum (left) and the squeezed limit bispectrum (right)
for the extended USR model ($\tau_{i}\to-\infty$). On the right panel, we also
mark the peak and trough scales of the power spectrum. For this figure we used
$g=1$ and $h=-10$ (similar to Fig. 1).

###  4.3 Probability of large fluctuations

We have seen that ${\cal{P}}_{{\cal{R}}_{s}}$ and ${f_{\text{NL}}}$ peak at
different scales. This raises the question of which scales exhibit the highest
probability of the realization of large fluctuations which is of great
importance for the estimation of the mass and abundance of PBHs. The full
answer to that question requires non-perturbative knowledge of the probability
distribution of smoothed density contrast in real space Hooshangi:2021ubn ;
Ezquiaga:2019ftu ; Biagetti:2021eep ; Musco:2018rwt . In contrast, our
approach only offers partial and perturbative access to the correlation
functions in Fourier space. In principle, it is possible to obtain real space
(but still perturbative) results if the full bispectrum is calculated but it
will be a complex task due to the high scale dependence and is beyond the
scope of this paper. In this section, we only try to make a rough estimate of
the effect of non-linearities by studying approximations to the full
probability distribution. Our results should be considered as speculations,
motivating a more careful analysis.

We start with defining a quantity $\beta(k_{s})$ as follows

$\displaystyle\beta(k_{s})\equiv\int_{1}^{\infty}\rho_{\hat{\cal{R}}}(k_{s})\,d\hat{\cal{R}}\
,$ (24)

where $\hat{\cal{R}}$ is a random variable (to be defined) with the
probability distribution function $\rho_{\hat{\cal{R}}}$. $k_{s}$ may be
understood as the variable that determines the smoothing comoving scale of
fluctuations in real space $R_{\text{smooth}}\simeq 1/k_{s}$. According to
that definition, $\beta(k_{s})$ is the probability of realization of the
random variable $\hat{\cal{R}}$ larger than a certain threshold (which is
chosen to be 1) at scale $R_{\text{smooth}}$.

Now we need to determine a probability distribution for $\hat{\cal{R}}$. A
natural relation for the random variable $\hat{{\cal{R}}}$ that takes into
account the leading order non-linearities may be expressed by

$\displaystyle\hat{\cal{R}}\equiv{\cal{R}}_{g}+\frac{3}{5}f^{\text{eff}}_{\text{NL}}(k_{s})\left({\cal{R}}^{2}_{g}-\sigma^{2}(k_{s})\right).$
(25)

where ${\cal{R}}_{g}$ is a Gaussian random variable with mean zero and
variance $\sigma^{2}(k_{s})={\cal{P}}_{{\cal{R}}_{s}}(k_{s})$ and
$f^{\text{eff}}_{\text{NL}}(k_{s})$ is an effective size of non-linearity
which shall be discussed shortly. This relation resembles the expression for
the local non-Gaussianity in real space Komatsu:2001rj  which studying its
effect on the probability distribution is straightforward Byrnes:2012yx . An
expression like this (and with $f^{\text{eff}}_{\text{NL}}={f_{\text{NL}}}$)
is justified in a scale-invariant scenario but it is also used in models that
predict strong violation of scale invariance Cai:2018dig . Since we have only
access to the squeezed limit bispectrum by the consistency relations and, in
the scenarios of interest, its scale dependence is very different from the
scale dependence of the power spectrum, it may be hard to justify Eq. (25).
However, we try to roughly compensate these defects by choosing an appropriate
expression for $f^{\text{eff}}_{\text{NL}}$. To achieve that, we require Eq.
(25) to roughly give the same skewness to the probability distribution of the
field as it is generated by a non-zero three-point correlation function. We
assume a sharp, top-hat window function in Fourier space which significantly
suppresses the contributions from the modes with momenta greater than $k_{s}$
LoVerde:2007ri ; Franciolini:2018vbk . In this case, the main contribution to
the skewness comes from the squeezed configuration but due to the large
difference between the long and the short power spectra,
$f^{\text{eff}}_{\text{NL}}$ receives an extra factor as follows

$\displaystyle
f^{\text{eff}}_{\text{NL}}(k_{s})\simeq\big{(}\dfrac{{\cal{P}}_{{\cal{R}}_{\ell}}}{{\cal{P}}_{{\cal{R}}_{s}}}\big{)}{f_{\text{NL}}}(k_{s}).$
(26)

In this relation, possible ${\cal O}(1)$ factors are omitted but the prefactor
is fixed by requiring the equality
$f^{\text{eff}}_{\text{NL}}={f_{\text{NL}}}$ in the scale-invariant limit. A
more accurate result consists of integrals over long and short modes which are
neglected here. This approximation can be roughly interpreted as considering
the effect of a narrow range of long modes on a narrow range of short modes
responsible for the PBH formation. The factor behind ${f_{\text{NL}}}$ is
reasonable because a long mode with a smaller amplitude is expected to have a
smaller effect on the short modes for a fixed strength of correlation.

According to Eq. (26), $f^{\text{eff}}_{\text{NL}}$ is suppressed
significantly compared to ${f_{\text{NL}}}$ in the transient USR model due to
the large enhancement of the short modes. This suppression makes the non-
linear term in Eq. (25) ineffective; the resulting $\beta$ would be very close
to the Gaussian case. On the other hand, in the extended model, the peak mode
and the long mode are both enhanced by the USR phase so that
$f^{\text{eff}}_{\text{NL}}\sim{f_{\text{NL}}}$ around the peak. Due to that
distinction, in this section, we only study the extended model which is more
interesting concerning non-linearities. (However, we also note that this is
only the partial result and it would be interesting to study the full effect
of non-linearities in both models.) Crucially, for the extended model,
$f^{\text{eff}}_{\text{NL}}$ becomes larger than ${f_{\text{NL}}}$ at scales
shorter than the first peak since
${\cal{P}}_{{\cal{R}}_{s}}<{\cal{P}}_{{\cal{R}}_{\ell}}$ at these scales (see
Fig. 2).

Fig. 3 depicts the behavior of $\beta_{{f_{\text{NL}}}}$ — which is $\beta$ as
defined in Eq. (24) using the relation (25) — as a function of $k_{s}$ for the
extended USR model. For comparison, we also show the behavior of $\beta_{G}$
which is the result of the linear theory (${f_{\text{NL}}}\to 0$). We observe
that the presence of non-linearities has three crucial effects: (i) it
increases the probability of large fluctuations, (ii) it shifts and broadens
the window of the PBH formation associated with the first peak and (iii) a new
window opens up due to the large $f^{\text{eff}}_{\text{NL}}$ at shorter
scales. While (i) is already a well-known effect, to our knowledge, (ii) and
(iii) seem to be overlooked in the previous literature. We may expect (ii) and
(iii) to lead to a non-monochromatic mass spectrum of PBHs and also cause the
mass of the most dominant PBHs to be different from what is predicted by the
linear theory. In Sec. 3.3 we have argued that ${f^{\prime}_{\text{NL}}}>0$ at
the peak. Thus, the size of ${f_{\text{NL}}}$ (which is positive according to
the discussion of Sec. 3.1) is larger after the peak than it is before the
peak. This implies that the presence of non-Gaussianity must shift the window
of large $\beta$ from around the peak of the power spectrum to the shorter
scales. This generic prediction of the consistency relation is clearly
confirmed in this example as is visualized in Fig. 3.999As another interesting
consequence of the consistency relations, consider the situation where the
slow-roll consistency relation Eq. (13) holds in which case ${f_{\text{NL}}}$
vanishes at the peak. Since ${f^{\prime}_{\text{NL}}}>0$ at the peak,
${f_{\text{NL}}}$ is negative slightly before and positive slightly after the
peak. This may lead to a suppression (enhancement) of $\beta$ at scales
slightly larger (shorter) than the peak. This is a situation where the
presence of non-linearities can have a dual effect; it causes the window of
large $\beta$ to shrink from the left and expand to the right.

As a further remark related to (iii), let us point out that it is also
possible that non-linearities open up new windows for PBH production which
might be missed in the perturbative analysis. As an example, consider large
scales ($k_{s}\to 0$) where scale invariance is restored and the validity of
the $\delta N$-formalism is justified. Using the $\delta N$-formalism in its
non-perturbative form, one can compute $\beta$ non-perturbatively
Hooshangi:2021ubn . (See Ref. Cai:2022erk  for the explicit $\delta N$
relations for the model under consideration.) For the same choice of
parameters as for Fig. 3, one obtains $\beta\simeq 2.5\times 10^{-8}$ from
this non-perturbative result, whereas at at the same scales we have
$\beta_{{f_{\text{NL}}}}\simeq 10^{-15}$ and $\beta_{G}\simeq 2.5\times
10^{-29}$. This shows the significance of non-perturbative effects.
Furthermore, comparing this non-perturbative result with the numbers in Fig. 3
suggests that an abundant production of PBHs is possible at large scales, far
away from the peak scale of the power spectrum.

We conclude that estimating the mass and abundance of PBHs from studying the
first peak of the power spectrum may not always be justified. However, as
stated before, we note that our analysis is oversimplified in various ways. A
thorough study of the PBH formation in such scenarios remains an interesting,
open problem.

Figure 3: The probability of larger-than-one fluctuations $\beta$ as a
function of scale for the extended USR model. Note that we multiplied
$\beta_{G}$ and $\beta_{{f_{\text{NL}}}}$ by different factors to visualize
them simultaneously. We used the same power spectrum as depicted in 2 (with
$g=1$ and $h=-10$) and used ${\cal{P}}_{{\cal{R}}_{\ell}}=8\times 10^{-3}$.
The small circle and triangle on the blue curve indicate the location of the
first peak of ${\cal{P}}_{{\cal{R}}_{s}}$ and ${f_{\text{NL}}}$, respectively.
(The peaks of ${f_{\text{NL}}}$ and $f^{\text{eff}}_{NL}$ occur at
$x_{k_{s}}\simeq 3.62$ and $x_{k_{s}}\simeq 3.77$, respectively.)

###  4.4 Connected and disconnected trispectra

It would also be interesting to study the effect of largely scale-dependent
power spectrum and bispectrum on the generation of stochastic GWs. We continue
to be speculative, make an order of magnitude estimate of the effect and leave
a complete analysis to future work.

The amplitude of the induced tensor mode is sourced by the second order scalar
perturbation. Therefore, the corresponding energy density of the induced GWs
is controlled by the four-point correlation function — i.e., the trispectrum —
of ${\cal{R}}_{s}$. There are two contributions to the four-point function
both of which may contribute to the induced GWs, namely, the disconnected and
connected trispectrum. The disconnected trispectrum originates from the linear
theory whereas the connected one requires non-linear interactions. If the two
contributions become comparable, the effect of non-linearities becomes non-
negligible Cai:2018dig . As an order of magnitude comparison, we consider
${\cal{P}}_{{\cal{R}}_{s}}^{2}$ as the size of the disconnected contribution
and ${\cal{P}}_{{\cal{R}}_{s}}^{2}{\cal{P}}_{{\cal{R}}_{\ell}}f^{2}_{NL}$ as
the size of the connected contribution. The latter estimation may be justified
by considering a connected diagram with four ${\cal{R}}_{s}$ as the external
legs and ${\cal{R}}_{\ell}$ as the mediator with a cubic long-short
interaction with strength ${f_{\text{NL}}}$. Comparing the two pieces, we
roughly conclude that if ${\cal{P}}_{{\cal{R}}_{\ell}}f^{2}_{NL}(k_{s})\gtrsim
1$, the effect of non-linearities cannot be neglected at the comoving scales
around $1/k_{s}$. Fig. 4 illustrates the behavior of this factor for a
specific choice of parameters. Notably, this figure suggests that, if the
squeezed limit bispectrum is to be important, it ought to be so at scales
different from the first peak of the power spectrum. We stress that the full
result beyond this order of magnitude estimation requires the full trispectrum
at all configurations (see e.g., Ragavendra:2021qdu ; Zhang:2021vak ), the
calculation of which is out of the scope of this paper.

Figure 4: The behavior of ${\cal{P}}_{{\cal{R}}_{\ell}}f^{2}_{NL}$ (the
contribution of non-linearities to the induced GWs, relative to the linear
contribution) as a function of scale for the extended USR model. For this
figure, we used $g=1$, $h=-10$ and ${\cal{P}}_{{\cal{R}}_{\ell}}=8\times
10^{-3}$.

##  5 Conclusion

In this paper, we adapted the general single-field consistency condition
advocated in Ref. Namjoo_2023  (Eq. (7)) for the purpose of analyzing
situations where significant deviation from scale invariance occurs. The main
motivation for such scenarios is to generate a large abundance of primordial
black holes (PBHs) or observable stochastic gravitational waves (GWs). Our
main relations are Eqs. (10) and (11) which have important implications for
both phenomena. Our consistency relation allowed us to make general and model-
independent statements about the squeezed limit bispectrum at different scales
and, in particular, about the size, sign and slope of ${f_{\text{NL}}}$ at the
peaks of the power spectrum.

We have argued that non-linearities are capable of shifting or expanding the
range of scales responsible for the GW and PBH production and may also open up
new windows for both phenomena. The full answer to the question of which range
of scales plays the main role in generating PBHs and GWs requires a more
complete analysis and remains an open problem. Specific statements on this
matter are likely to be model-dependent, thus studying explicit examples such
as the ones considered here or in Refs. Hooshangi:2023kss ; Domenech:2023dxx ;
Pi:2022ysn ; Cai:2022erk ; Cai:2021zsp ; Ragavendra:2020sop ; Lin:2020goi ;
Saito:2008em ; Gao:2020tsa  would be fruitful.

Our general consistency relation shows that the value of ${f_{\text{NL}}}$ at
the largest peak of the power spectrum can be different from the predictions
of the $\delta N$-formalism or any other methods that rely on the gradient
expansion for the peak scale firouzjahi2023sign ; Iacconi:2023ggt . Thus,
claims like that the behavior of the tail as predicted by the $\delta
N$-formalism applies to the peak may not be justified Cai:2022erk . This may
be expected because these methods do not carry scale-dependent information
while the situation under consideration is highly scale-dependent. The $\delta
N$-formalism may be improved by going beyond the leading order in the gradient
expansion Takamizu:2010xy ; Jackson:2023obv ; Ozsoy:2019lyy . This may allow
us to study non-linearities at the first peak but not at the next peaks that
appear at shorter scales. Crucially, our study of the explicit examples
suggests that the scales shorter than the first peak can be as important in
generating PBHs or GWs due to the enhancement of the tail of the distribution
by non-linearities.

It would be crucial to go beyond the three-point correlation function and it
would be interesting to explore the possibility of extending the consistency
relation (Eq. (7)) for that purpose. In the specific case of four-point
function in the collapse limit, one may expect to have
$\tau_{\text{NL}}\propto f^{2}_{\text{NL}}$ which allows us to make general
statements about $\tau_{\text{NL}}$ from what we concluded for
${f_{\text{NL}}}$ Ozsoy:2021pws . We leave these studies to future works.

###### Acknowledgements.

We thank Hassan Firouzjahi and Mehrdad Mirbabayi for fruitful discussions.
This work is supported by INSF of Iran under the grant number 4022911.

##  Appendix A Explicit in-in Calculation

In this appendix, we provide some details of the model studied in Sec. 4 and
also calculate the non-Gaussianity of the extended USR model of Sec. 4.2,
using the standard in-in formalism.

###  A.1 Background evolution

The background equations of motion in single-field inflation are

$\ddot{\phi}+3H\dot{\phi}+V_{,\phi}(\phi)=0,\quad H^{2}=\frac{1}{3M_{\rm
Pl}^{2}}\big{(}\frac{1}{2}\dot{\phi}^{2}+V(\phi)\big{)},$ (27)

where the overdot denotes derivatives with respect to the cosmic time. Using
the number of e-folds ($N\propto\log(a)$) as the clock, the first equation can
be rewritten by

$\dfrac{\mathrm{d}^{2}\phi}{\mathrm{d}N^{2}}+3\dfrac{\mathrm{d}\phi}{\mathrm{d}N}+\dfrac{V^{\prime}(\phi)}{V(\phi)}\simeq
0\ ,$ (28)

where an $\epsilon$-suppressed term is neglected. The slow-roll parameters are
defined by

$\epsilon\equiv-\frac{\dot{H}}{H^{2}}=\frac{1}{2}\big{(}\dfrac{\mathrm{d}\phi}{\mathrm{d}N}\big{)}^{2},\quad\eta\equiv\frac{\dot{\epsilon}}{H\epsilon}.$
(29)

In the transient USR model with the potential as presented in Eq. (14), the
inflaton field begins its journey from an SR phase and reaches $\phi_{i}$, the
value at the first phase transition, at the conformal time $\tau_{i}$. Thus,
we have

$\text{For}\quad\tau\leq\tau_{i}:\quad\epsilon(\tau)\simeq\epsilon_{i},\quad\eta(\tau)\approx
0.$ (30)

Then it transits sharply to the USR phase after which $\eta\simeq-6$. Thus at
the transition, the SR parameters can be modeled as follows

$\text{For}\quad\tau_{i-}\leq\tau\leq\tau_{i+}:\quad\epsilon(\tau)=\epsilon_{i},\quad\eta(\tau)=-6\theta(\tau-\tau_{i}).$
(31)

In the USR phase, one can find the SR parameters by solving Eq. (28) for a
constant potential,

$\text{For}\quad\tau_{i}\leq\tau\leq\tau_{e}:\quad\epsilon(\tau)=\epsilon_{i}\big{(}\frac{\tau}{\tau_{i}}\big{)}^{6},\quad\eta(\tau)=-6,$
(32)

where $\tau_{e}$ is the comformal time when $\phi_{e}$ is reached, where the
USR phase terminates. Thus the SR parameters just at the end of the USR stage
would be
$\epsilon(\tau_{e-})=\epsilon_{i}\big{(}\frac{\tau_{e}}{\tau_{i}}\big{)}^{6}\equiv\epsilon_{e}$
and $\eta(\tau_{e-})=-6$.

After the USR phase, the inflation reaches a sharp step exactly located at
$\phi_{e}$. As a result, the inflaton speeds up or down depending on the sign
of $\kappa$ in Eq. (14). Denoting the SR parameter right after the step by
$\epsilon_{s}$, we can parameterize this change of speed with
$g\equiv\sqrt{\frac{\epsilon_{s}}{\epsilon_{e}}}$. Around the second phase
transition, the first slow-roll parameter is discontinuous resulting in

$\displaystyle\text{For}\quad\tau_{e-}\leq\tau\leq\tau_{e+}:$
$\displaystyle\epsilon(\tau)=\epsilon_{e}\theta(\tau_{e}-\tau)+\epsilon_{s}\theta(\tau-\tau_{e})\
,$ (33)
$\displaystyle\eta(\tau)=-6\theta(\tau_{e}-\tau)-(h+6)\theta(\tau-\tau_{e})-\frac{\tau}{\epsilon(\tau)}\delta(\tau-\tau_{e})\left(\epsilon_{s}-\epsilon_{e}\right),$

where $h\equiv-6\sqrt{\frac{\epsilon_{V}}{\epsilon_{s}}}$.

At the final SR phase, the inflaton’s trajectory relaxes to its new attractor.
The SR parameters during this relaxation time can be determined by solving the
Klein-Gordon equation Eq. (28) with the potential Eq. (14) for
$\phi<\phi_{e}$:

$\text{For}\quad\tau_{e}\leq\tau:\quad\epsilon(\tau)=\epsilon_{s}\left(\frac{h}{6}-(1+\frac{h}{6})\big{(}\frac{\tau}{\tau_{e}}\big{)}^{3}\right)^{2},\quad\eta(\tau)=-\frac{6(6+h)}{(6+h)-h\left(\frac{\tau_{e}}{\tau}\right)^{3}}.$
(34)

Now with all the SR parameters in hand, we can determine the mode functions
for each stage of inflation.

###  A.2 Mode Functions

The comoving curvature perturbation obeys the following equation

$\mathcal{R}_{k}^{\prime\prime}(\tau)-\dfrac{(2+\eta)}{\tau}\mathcal{R}_{k}^{\prime}(\tau)+k^{2}\mathcal{R}_{k}(\tau)=0.$
(35)

Using Eq. (35) along with the Bunch-Davies vacuum, one obtains the mode
function of the first slow-roll period as follows,

$\mathcal{R}^{(1)}_{k}(\tau)=\frac{H}{\sqrt{4\epsilon
k^{3}}}(1+ik\tau)e^{-ik\tau}.$ (36)

To obtain the mode functions for the USR and the second SR phases, we quantify
them using the Bogolyubov coefficients Cai_2018 ; Cai:2022erk ;
Firouzjahi:2023aum

$\mathcal{R}_{k}^{(2),(3)}(\tau)=\alpha^{(2),(3)}_{k}\frac{H}{\sqrt{4\epsilon
k^{3}}}(1+ik\tau)e^{-ik\tau}+\beta^{(2),(3)}_{k}\frac{H}{\sqrt{4\epsilon
k^{3}}}(1-ik\tau)e^{ik\tau},$ (37)

where the superscripts $(2)$ and $(3)$ represent the coefficients for the
second (USR) and third (SR) phases of the evolution, respectively. This mode
function is an approximate solution of Eq. (35) as long as $H$ is kept a
constant but $\epsilon$ can have an arbitrary time dependence Cai_2018 . The
Bogolyubov coefficients may be determined by the continuity of ${\cal{R}}$ and
$\epsilon{\cal{R}}^{\prime}$ at each phase transition Namjoo_2012  which
result in

$\alpha_{k}^{(2)}=1-\frac{3i}{2x_{k}^{3}x_{i}^{3}}\left(1+x_{k}^{2}x_{i}^{2}\right),\quad\beta_{k}^{(2)}=\frac{3i}{2x_{k}^{3}x_{i}^{3}}\left(1-ix_{k}x_{i}\right)^{2}e^{2ix_{k}x_{i}},$
(38)

and

$\begin{pmatrix}\alpha^{(3)}_{k}\\\ \\\
\beta^{(3)}_{k}\end{pmatrix}=\frac{1}{4gx_{k}^{3}}\begin{pmatrix}\gamma_{11}\quad\gamma_{12}\\\
\\\ \gamma_{21}\quad\gamma_{22}\end{pmatrix}\begin{pmatrix}\alpha^{(2)}_{k}\\\
\\\ \beta^{(2)}_{k}\end{pmatrix},$ (39)

where $x_{k}=-k\tau_{e}$, $x_{i}=\frac{\tau_{i}}{\tau_{e}}$ and

$\begin{gathered}\gamma_{11}=\gamma^{*}_{22}=-i\left[-6+g^{2}(6+h)+\left(-4+g^{2}(4+h)\right)x_{k}^{2}+2i(1+g^{2})x_{k}^{3}\right],\\\
\\\
\gamma_{12}=\gamma^{*}_{21}=e^{-2ix}\left(x_{k}-i\right)\left[-6+g^{2}(6+h)+i(-6+g^{2}(6+h))x_{k}+2(1-g^{2})x_{k}^{2}\right].\end{gathered}$
(40)

For $g=1$ i.e., when there is no step, $\alpha^{(3)}_{k}$ and
$\beta^{(3)}_{k}$ reduce to the Bogolyubov coefficients reported in
firouzjahi2023sign . The final power spectrum is given by

$\mathcal{P}_{\mathcal{R}}=\dfrac{H^{2}}{8\pi^{2}\epsilon_{V}}|\alpha_{k}^{(3)}+\beta_{k}^{(3)}|^{2}.$
(41)

The power spectrum for the extended USR phase can be obtained by taking the
limit $k\tau_{i}\to-\infty$ which results in $\alpha^{(2)}_{k}\to
1,\beta_{k}^{(2)}\to 0$ and

$\alpha^{(3)}_{k}=\frac{\gamma_{11}}{4gx_{k}^{3}},\quad\beta^{(3)}_{k}=\frac{\gamma_{21}}{4gx_{k}^{3}}.$
(42)

Hence, the final power spectrum in the extended USR model is

$\mathcal{P}_{\mathcal{R}}=\dfrac{H^{2}}{32\pi^{2}g^{2}x_{k}^{6}\epsilon_{V}}|\Gamma|^{2},$
(43)

where

$\Gamma=\frac{e^{-ix_{k}}}{2}\left(\gamma_{11}+\gamma_{21}\right).$ (44)

###  A.3 Non-Gaussianity

We are now prepared to calculate the three-point function for the extended USR
model. The leading contribution to the three-point function comes from the
following term in the cubic action Namjoo:2012aa

$S_{3}\supset\int d\tau
d^{3}x\frac{a^{2}\epsilon}{2}\eta^{\prime}\mathcal{R}^{2}\mathcal{R}^{\prime},$
(45)

which leads to the following expression for the bispectrum

$B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)=-\imaginary\mathcal{R}_{k_{1}}\left(\tau_{0}\right)\mathcal{R}_{k_{2}}\left(\tau_{0}\right)\mathcal{R}_{k_{3}}\left(\tau_{0}\right)\int_{-\infty}^{\tau_{0}}d\tau
a^{2}\epsilon\eta^{\prime}\left[\mathcal{R}_{k_{1}}^{*}(\tau)\mathcal{R}_{k_{2}}^{*}(\tau)\mathcal{R}_{k_{3}}^{*\prime}(\tau)+2\text{perm.}\right].$
(46)

We can break down the integral in Eq. (46) into different pieces for different
stages of inflation. In the USR phase, $\eta=-6$ is constant, and
$\eta^{\prime}=0$. Thus, $H_{\mathrm{int}}$ in Eq. (45) is negligible, and
there would be no significant non-Gaussianity from the USR phase. Thus, we are
left with the effect of the sharp USR-SR transition and the relaxation period
after the transition which we calculate their contribution separately.

To compute the bispectrum from the transition ($B^{T}_{\mathcal{R}}$) — i.e.,
the contribution to Eq. (46) from the range $\tau_{e-}\leq\tau\leq\tau_{e+}$ —
it is hard to work directly with $H_{\text{int}}$ as given in Eq. (45) mainly
due to the existence of terms proportional to $\delta^{\prime}$. To circumvent
this difficulty, we perform consecutive integrations by part and use the
equation of motion Eq. (35) to remove higher order time derivatives until no
term containing delta-functions or its derivatives remains. This procedure is
allowed as long as we keep track of the boundary terms Arroja:2011yj . Then,
all integrals become regular and vanish in the limit $\tau_{e-}\to\tau_{e+}$
and we are only left with the boundary terms. There are three boundary terms
as follow

$\int_{\tau_{e}-}^{\tau_{e}+}\frac{a^{2}\epsilon}{2}\eta^{\prime}\mathcal{R}^{2}\mathcal{R}^{\prime}=\frac{1}{2}\left[a^{2}\epsilon\left(\eta\mathcal{R}^{2}\mathcal{R}^{\prime}+\tau\mathcal{R}^{2}\nabla^{2}\mathcal{R}-2\tau\mathcal{R}\mathcal{R}^{\prime
2}\right)\right]_{\tau_{e-}}^{\tau_{e+}}\ ,$ (47)

which result in

$\displaystyle
B_{\mathcal{R}}^{T}=a^{2}\imaginary\\{\mathcal{R}^{*}_{k_{1}}(0)\mathcal{R}^{*}_{k_{2}}(0)\mathcal{R}^{*}_{k_{3}}(0)\big{[}\left(\epsilon\eta\mathcal{R}_{k_{1}}(\tau)\mathcal{R}_{k_{2}}(\tau)\mathcal{R}^{\prime}_{k_{3}}(\tau)\right.-\epsilon\tau
k_{1}^{2}\mathcal{R}_{k_{1}}(\tau)\mathcal{R}_{k_{2}}(\tau)\mathcal{R}_{k_{3}}(\tau)$
(48)
$\displaystyle\left.-2\epsilon\tau\mathcal{R}_{k_{1}}(\tau)\mathcal{R}^{\prime}_{k_{2}}(\tau)\mathcal{R}^{\prime}_{k_{3}}(\tau)\right)+2\mathrm{perms.}\big{]}_{\tau_{e-}}^{\tau_{e+}}\biggl{\\}}.$

Next, we need to perform the integration of Eq. (46) for the relaxation period
(i.e., for the range $\tau_{e}<\tau\leq 0$). Using Eq. (34) we obtain

$\displaystyle
B_{\mathcal{R}}^{R}=\imaginary\left[\frac{-729ie^{-i(k_{1}+k_{2}+k_{3})\tau}H^{4}\tau_{e}^{6}(6+h)}{g^{4}h^{2}k_{1}^{3}k_{2}^{3}k_{3}^{3}\epsilon_{e}\left((6+h)\tau^{3}-h\tau_{e}^{3}\right)^{3}}\prod_{i=1}^{3}\Lambda_{k_{i}}(\tau)\right]_{\tau=0}^{\tau=\tau_{e}}\
,$ (49)

where

$\displaystyle\Lambda_{k_{i}}(\tau)=(\alpha^{(3)}_{k_{i}}+\beta^{(3)}_{k_{i}})\left(\beta^{(3)*}_{k_{i}}(i-k_{i}\tau)+e^{2ik_{i}\tau}\alpha^{(3)*}_{k_{i}}(i+k_{i}\tau)\right).$
(50)

The final bispectrum is the sum of the two contributions
$B_{\mathcal{R}}=B^{T}_{\mathcal{R}}+B^{R}_{\mathcal{R}}$ which, after taking
the squeezed limit and using Eq. (3), results in

$\displaystyle{f_{\text{NL}}}=\Omega$
$\displaystyle\biggl{\\{}6\left(-6+g^{2}(6+h)\right)\left(-3-3g^{2}+g^{4}(6+h)\right)+12\left(-1+g^{2}\right)^{2}\left(1+g^{2}\right)x_{k}^{6}$
(51)
$\displaystyle+2\left(72+5g^{6}(4+h)(6+h)-3g^{2}(8+5h)-7g^{4}(24+5h)\right)x_{k}^{2}$
$\displaystyle+4\left(4-g^{2}(4+h)\right)\left(3+g^{2}-g^{4}(4+h)\right)x_{k}^{4}$
$\displaystyle+\biggl{[}3\left(-6+g^{2}(6+h)\right)\left(3+3g^{2}-g^{4}(6+h)\right)+2\left(3-g^{2}+7g^{4}-g^{6}(9+h)\right)x_{k}^{6}$
$\displaystyle\quad\,\,+\left(36-3g^{2}(-8+h)+g^{6}(6+h)(16+h)-g^{4}(156+19h)\right)x_{k}^{2}$
$\displaystyle\quad\,\,+\left(12+g^{2}(4-6h)-8g^{4}(19+3h)+g^{6}(136+3h(14+h))\right)x_{k}^{4}\biggl{]}2\cos(2x_{k})$
$\displaystyle+\biggl{[}\left(g^{4}(50h+204)+12g^{2}(2h+5)-72-\left(g^{6}(h+6)(9h+32)\right)\right)x_{k}^{2}$
$\displaystyle\,\,\,\,\,\,\,\,+\left(g^{6}(h(h+30)+140)-2g^{4}(7h+72)+28g^{2}-24\right)x_{k}^{4}-4g^{2}\left(g^{4}-1\right)x_{k}^{6}$
$\displaystyle\,\,\quad-12\left(g^{6}(h+6)^{2}-9g^{4}(h+6)-3g^{2}h+18\right)\biggl{]}x_{k}\sin(2x_{k})\biggl{\\}}\,,$

where

$\displaystyle\Omega\equiv\frac{5h}{12\left(-6+g^{2}h\right)|\Gamma|^{2}}.$
(52)

Remarkably, the above $f_{NL}$ matches the one we have derived via the
consistency relation Eq. (20).

## References

  * (1) O. Özsoy and G. Tasinato, _Inflation and Primordial Black Holes_ , _Universe_ 9 (2023) 203 [2301.03600].
  * (2) M.P. Hertzberg and M. Yamada, _Primordial Black Holes from Polynomial Potentials in Single Field Inflation_ , _Phys. Rev. D_ 97 (2018) 083509 [1712.09750].
  * (3) J. Garcia-Bellido and E. Ruiz Morales, _Primordial black holes from single field models of inflation_ , _Phys. Dark Univ._ 18 (2017) 47 [1702.03901].
  * (4) G. Domènech, _Scalar Induced Gravitational Waves Review_ , _Universe_ 7 (2021) 398 [2109.01398].
  * (5) J.S. Bullock and J.R. Primack, _NonGaussian fluctuations and primordial black holes from inflation_ , _Phys. Rev. D_ 55 (1997) 7423 [astro-ph/9611106].
  * (6) R.-g. Cai, S. Pi and M. Sasaki, _Gravitational Waves Induced by non-Gaussian Scalar Perturbations_ , _Phys. Rev. Lett._ 122 (2019) 201101 [1810.11000].
  * (7) S. Young and C.T. Byrnes, _Signatures of non-gaussianity in the isocurvature modes of primordial black hole dark matter_ , _JCAP_ 04 (2015) 034 [1503.01505].
  * (8) V. Atal and C. Germani, _The role of non-gaussianities in Primordial Black Hole formation_ , _Phys. Dark Univ._ 24 (2019) 100275 [1811.07857].
  * (9) G. Franciolini, A. Iovino, Junior., V. Vaskonen and H. Veermae, _Recent Gravitational Wave Observation by Pulsar Timing Arrays and Primordial Black Holes: The Importance of Non-Gaussianities_ , _Phys. Rev. Lett._ 131 (2023) 201401 [2306.17149].
  * (10) B. Finelli, G. Goon, E. Pajer and L. Santoni, _Soft Theorems For Shift-Symmetric Cosmologies_ , _Phys. Rev. D_ 97 (2018) 063531 [1711.03737].
  * (11) K. Hinterbichler, L. Hui and J. Khoury, _Conformal Symmetries of Adiabatic Modes in Cosmology_ , _JCAP_ 08 (2012) 017 [1203.6351].
  * (12) R. Bravo, S. Mooij, G.A. Palma and B. Pradenas, _A generalized non-Gaussian consistency relation for single field inflation_ , _JCAP_ 05 (2018) 024 [1711.02680].
  * (13) V. Assassi, D. Baumann and D. Green, _On Soft Limits of Inflationary Correlation Functions_ , _JCAP_ 11 (2012) 047 [1204.4207].
  * (14) C. Cheung, A.L. Fitzpatrick, J. Kaplan and L. Senatore, _On the consistency relation of the 3-point function in single field inflation_ , _JCAP_ 02 (2008) 021 [0709.0295].
  * (15) N. Kundu, A. Shukla and S.P. Trivedi, _Ward Identities for Scale and Special Conformal Transformations in Inflation_ , _JHEP_ 01 (2016) 046 [1507.06017].
  * (16) M.H. Namjoo, _One consistency relation for all single-field inflationary models_ , 2311.12777.
  * (17) J. Maldacena, _Non-gaussian features of primordial fluctuations in single field inflationary models_ , _Journal of High Energy Physics_ 2003 (2003) 013–013.
  * (18) P. Creminelli and M. Zaldarriaga, _Single field consistency relation for the 3-point function_ , _JCAP_ 10 (2004) 006 [astro-ph/0407059].
  * (19) Y.-F. Cai, X. Chen, M.H. Namjoo, M. Sasaki, D.-G. Wang and Z. Wang, _Revisiting non-gaussianity from non-attractor inflation models_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (2018) 012.
  * (20) M.H. Namjoo, H. Firouzjahi and M. Sasaki, _Violation of non-Gaussianity consistency relation in a single field inflationary model_ , _EPL_ 101 (2013) 39001 [1210.3692].
  * (21) H. Firouzjahi and A. Riotto, _The sign of non-gaussianity and the primordial black holes abundance_ , 2023.
  * (22) X. Chen, H. Firouzjahi, M.H. Namjoo and M. Sasaki, _A Single Field Inflation Model with Large Local Non-Gaussianity_ , _EPL_ 102 (2013) 59001 [1301.5699].
  * (23) Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang and Z. Zhou, _Highly non-Gaussian tails and primordial black holes from single-field inflation_ , _JCAP_ 12 (2022) 034 [2207.11910].
  * (24) Y. Tada, T. Terada and J. Tokuda, _Cancellation of quantum corrections on the soft curvature perturbations_ , 2023.
  * (25) M.H. Namjoo, H. Firouzjahi and M. Sasaki, _Multiple inflationary stages with varying equation of state_ , _Journal of Cosmology and Astroparticle Physics_ 2012 (2012) 018–018.
  * (26) S. Hooshangi, M.H. Namjoo and M. Noorbala, _Rare events are nonperturbative: Primordial black holes from heavy-tailed distributions_ , _Phys. Lett. B_ 834 (2022) 137400 [2112.04520].
  * (27) J.M. Ezquiaga, J. García-Bellido and V. Vennin, _The exponential tail of inflationary fluctuations: consequences for primordial black holes_ , _JCAP_ 03 (2020) 029 [1912.05399].
  * (28) M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, _The formation probability of primordial black holes_ , _Phys. Lett. B_ 820 (2021) 136602 [2105.07810].
  * (29) I. Musco, _Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations_ , _Phys. Rev. D_ 100 (2019) 123524 [1809.02127].
  * (30) E. Komatsu and D.N. Spergel, _Acoustic signatures in the primary microwave background bispectrum_ , _Phys. Rev. D_ 63 (2001) 063002 [astro-ph/0005036].
  * (31) C.T. Byrnes, E.J. Copeland and A.M. Green, _Primordial black holes as a tool for constraining non-Gaussianity_ , _Phys. Rev. D_ 86 (2012) 043512 [1206.4188].
  * (32) M. LoVerde, A. Miller, S. Shandera and L. Verde, _Effects of Scale-Dependent Non-Gaussianity on Cosmological Structures_ , _JCAP_ 04 (2008) 014 [0711.4126].
  * (33) G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, _Primordial Black Holes from Inflation and non-Gaussianity_ , _JCAP_ 03 (2018) 016 [1801.09415].
  * (34) H.V. Ragavendra, _Accounting for scalar non-Gaussianity in secondary gravitational waves_ , _Phys. Rev. D_ 105 (2022) 063533 [2108.04193].
  * (35) F. Zhang, J. Lin and Y. Lu, _Double-peaked inflation model: Scalar induced gravitational waves and primordial-black-hole suppression from primordial non-Gaussianity_ , _Phys. Rev. D_ 104 (2021) 063515 [2106.10792].
  * (36) S. Hooshangi, M.H. Namjoo and M. Noorbala, _Tail diversity from inflation_ , _JCAP_ 09 (2023) 023 [2305.19257].
  * (37) G. Domènech, G. Vargas and T. Vargas, _An exact model for enhancing/suppressing primordial fluctuations_ , 2309.05750.
  * (38) S. Pi and M. Sasaki, _Logarithmic Duality of the Curvature Perturbation_ , _Phys. Rev. Lett._ 131 (2023) 011002 [2211.13932].
  * (39) Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang and Z. Zhou, _One small step for an inflaton, one giant leap for inflation: A novel non-Gaussian tail and primordial black holes_ , _Phys. Lett. B_ 834 (2022) 137461 [2112.13836].
  * (40) H.V. Ragavendra, P. Saha, L. Sriramkumar and J. Silk, _Primordial black holes and secondary gravitational waves from ultraslow roll and punctuated inflation_ , _Phys. Rev. D_ 103 (2021) 083510 [2008.12202].
  * (41) J. Lin, Q. Gao, Y. Gong, Y. Lu, C. Zhang and F. Zhang, _Primordial black holes and secondary gravitational waves from $k$ and $G$ inflation_, _Phys. Rev. D_ 101 (2020) 103515 [2001.05909].
  * (42) R. Saito, J. Yokoyama and R. Nagata, _Single-field inflation, anomalous enhancement of superhorizon fluctuations, and non-Gaussianity in primordial black hole formation_ , _JCAP_ 06 (2008) 024 [0804.3470].
  * (43) Q. Gao, Y. Gong and Z. Yi, _Primordial black holes and secondary gravitational waves from natural inflation_ , _Nucl. Phys. B_ 969 (2021) 115480 [2012.03856].
  * (44) L. Iacconi, D. Mulryne and D. Seery, _Loop corrections in the separate universe picture_ , 2312.12424.
  * (45) Y.-i. Takamizu, S. Mukohyama, M. Sasaki and Y. Tanaka, _Non-Gaussianity of superhorizon curvature perturbations beyond $\delta$ N formalism_, _JCAP_ 06 (2010) 019 [1004.1870].
  * (46) J.H.P. Jackson, H. Assadullahi, A.D. Gow, K. Koyama, V. Vennin and D. Wands, _The separate-universe approach and sudden transitions during inflation_ , 2311.03281.
  * (47) O. Özsoy and G. Tasinato, _On the slope of the curvature power spectrum in non-attractor inflation_ , _JCAP_ 04 (2020) 048 [1912.01061].
  * (48) O. Özsoy and G. Tasinato, _Consistency conditions and primordial black holes in single field inflation_ , _Phys. Rev. D_ 105 (2022) 023524 [2111.02432].
  * (49) H. Firouzjahi, _One-loop corrections in power spectrum in single field inflation_ , _JCAP_ 10 (2023) 006 [2303.12025].
  * (50) F. Arroja and T. Tanaka, _A note on the role of the boundary terms for the non-Gaussianity in general k-inflation_ , _JCAP_ 05 (2011) 005 [1103.1102].