# Constraints on Non-Gaussian primordial curvature perturbation from the LIGO-
Virgo-KAGRA third observing run

Ryoto Inui  Santiago Jaraba  Sachiko Kuroyanagi  Shuichiro Yokoyama

###### Abstract

The scalar-induced gravitational waves, arising from large amplitude
primordial density fluctuations, provide a unique observational test for
directly probing the epoch of inflation. In this work, we provide constraints
on the SIGW background by taking into account the non-Gaussianity in the
primordial density fluctuations, using the third observing run (O3) data of
the LIGO-Virgo-KAGRA collaboration. We find that the non-Gaussianity gives a
non-negligible effect on the GW energy density spectrum and starts to affect
the analysis of the O3 data when the non-Gaussianity parameter is
$F_{\mathrm{NL}}>3.55$. Furthermore, the constraints exhibit asymptotic
behavior given by $F_{\mathrm{NL}}A_{g}=\rm{const.}$ at large
$F_{\mathrm{NL}}$ limit, where $A_{g}$ denotes the amplitude of the curvature
perturbations. In the limit of large $F_{\mathrm{NL}}$, we placed a 95%
confidence level upper limit $F_{\mathrm{NL}}A_{g}\leq 0.13,0.09,0.10$ at
fixed scales of $10^{16},10^{16.5},10^{17}~{}\mathrm{Mpc}^{-1}$, respectively.

##  1 Introduction

The first detection of gravitational waves from a compact binary coalescence
(CBC): GW150914 opened a new window to explore the universe [1]. Subsequent to
this landmark detection, the Laser Interferometer Gravitational Wave
Observatory (LIGO) and Virgo collaborations have identified approximately 90
GW events from CBCs up until the third observing run (O3). More recently,
multiple Pulsar Timing Array (PTA) experiments reported possible evidence of
the stochastic gravitational wave backgrounds, potentially originating from
the superposition of GWs from supermassive black hole binary systems [2, 3, 4,
5]. Nonetheless, ongoing discussions continue to search for other
astrophysical and cosmological sources of the SGWBs (see, e.g., Ref. [6]).

With such substantial observational progress, the cosmological SGWBs have
garnered significant attention as a means to explore the universe even before
the cosmic microwave background (CMB). Various mechanisms for generating such
SGWBs in the early universe prior to the CMB have been extensively studied.
These mechanisms include the production of primordial GWs originating from
inflation and a preheating phase [7, 8, 9, 10, 11], cosmic strings [12, 13,
14], the bubble collisions at the first-order phase transition [15, 16, 17,
18], and more. One of the most interesting mechanisms is the scalar-induced
gravitational waves (SIGWs), which can be generated from the primordial
curvature perturbations through the tensor-scalar coupling predicted in the
second-order cosmological perturbation theory (see, e.g., Refs. [19, 20, 21,
22]).

In general, it is challenging to directly probe small-scale primordial scalar
fluctuations. However, by investigating such SIGWs, we can indirectly explore
such primordial fluctuations even on much smaller scales than the large-scale
structure in our Universe (see, e.g., Ref. [23]). Probing such small-scale
primordial fluctuations should be valuable not only for testing inflationary
models but also for indirectly searching for primordial black holes [21].

In Refs. [24, 25], the upper limits on the amplitude of primordial curvature
perturbations have been placed by searching for SIGWs using the LIGO-Virgo-
KAGRAs O2 and O3 data.111The search for SIGWs in the recent pulsar timing
array (PTA) experiments has been frequently discussed (see, e.g., Refs. [6,
26, 27, 28, 29]). While the previous studies assumed a Gaussian distribution
of the primordial curvature perturbations, in this paper, we investigate the
impact of non-Gaussianity on the SIGW signal and provide new constraints on
non-Gaussian primordial curvature perturbations using the LVKs O3 data.

In fact, recent theoretical studies have shown that several inflationary
models can generate enhanced primordial curvature perturbations with non-
Gaussian distributions (see, e.g., Refs. [30, 31, 32, 33, 34, 35, 36, 37, 38,
39]). Such non-Gaussianity not only affects the abundance of PBHs (see, e.g.,
Refs. [40, 41]) but also the spectrum of the SIGWs [42, 43, 44, 45, 46].
Therefore, the observational constraints must be adjusted when considering
such inflationary models predicting non-Gaussian primordial curvature
perturbations. While the form of non-Gaussianity can vary significantly
depending on the inflationary model, in this paper, as an initial approach, we
employ the simplest non-linearity parameter, $F_{\rm NL}$, to characterize the
primordial non-Gaussianity, and the monochromatic power spectrum for the
primordial Gaussian fluctuations.

This paper is organized as follows. In Sec. 2, we give a brief review of SIGWs
in the presence of the primordial non-Gaussianity, based on Refs [47, 46, 45].
Then, we briefly describe the parameter estimation method for a SGWB searched
by the cross-correlation analysis in Sec. 3. Then we place observational
constraints using LVK O3 data in Sec. 4. The conclusion is given in Sec. 5.

##  2 Gravitational waves sourced by scalar perturbations

Here, we briefly review the SIGWs with primordial non-Gaussianity, based on
Ref. [46]. The perturbed metric in the conformal Newtonian gauge can be
written as

$\displaystyle\differential{s}^{2}=a(\tau)^{2}\left[-(1+2\Phi)\differential{\tau}^{2}+\left((1-2\Phi)\delta_{ij}+\frac{1}{2}h_{ij}\right)\differential{x}^{i}\differential{x}^{j}\right],$
(2.1)

where $\tau$ represents the conformal time, $\Phi$ is the curvature
perturbation in the Newtonian gauge, and $h_{ij}$ is the transverse traceless
tensor perturbation. Here, we have neglected the vector perturbations and the
anisotropic stress.

The tensor perturbation can be expanded by the Fourier modes
$h_{\lambda}(\tau,\mathbf{k})$

$\displaystyle
h_{ij}(\tau,\mathbf{x})=\sum_{\lambda=+,\times}\int\frac{\differential^{3}k}{(2\pi)^{3}}\mathrm{e}^{i\mathbf{k}\cdot\mathbf{x}}e^{\lambda}_{ij}h_{\lambda}(\tau,\mathbf{k}),$
(2.2)

where $e^{\lambda}_{ij}(k)$ denotes the transverse traceless polarization
tensors that satisfy
$e^{\lambda}_{ij}e^{\lambda^{\prime},ij}=\delta^{\lambda\lambda^{\prime}}$.

###  2.1 Scalar induced gravitational waves

At the second order, the Fourier mode of GWs follows the equation of motion
(EoM)

$\displaystyle
h^{\prime\prime}_{\lambda}(\tau,\mathbf{k})+2\mathcal{H}h^{\prime}_{\lambda}(\tau,\mathbf{k})+k^{2}h_{\lambda}(\tau,\mathbf{k})=4S_{\lambda}(\tau,\mathbf{k})\,,$
(2.3)

where $S_{\lambda}(\tau,\mathbf{k})$ is the source term which is proportional
to the quadratic terms of $\Phi$, $\mathcal{H}=a(\tau)H(\tau)$ is the
conformal Hubble parameter, and the prime denotes the partial derivative with
respect to $\tau$. In terms of the primordial curvature perturbation in the
uniform energy density slice on super-horizon scales, $\zeta$, the source term
can be expressed as

$\displaystyle
S_{\lambda}(\tau,\mathbf{k})=\int\frac{\differential^{3}q}{(2\pi)^{3}}Q_{\lambda}(\mathbf{k},\mathbf{q})f(\absolutevalue{\mathbf{k}-\mathbf{q}},q,\tau)\zeta(\mathbf{q})\zeta(\mathbf{k}-\mathbf{q})\,.$
(2.4)

Here, the projection factor $Q_{\lambda}(\mathbf{k},\mathbf{q})$ is given by

$\displaystyle
Q_{\lambda}(\mathbf{k},\mathbf{q})=e^{\lambda}_{ij}(\mathbf{k})q^{i}q^{j}=\frac{q^{2}}{\sqrt{2}}\sin^{2}\theta\times\begin{cases}\cos(2\psi)&(\lambda=+)\\\
\sin(2\psi)&(\lambda=\times)\end{cases}\,,$ (2.5)

where we work in spherical coordinates and set $\mathbf{k}$ in the
$z$-direction. The source factor $f(p,q,\tau)$ in the radiation-dominated (RD)
universe is given by

$\displaystyle f(p,q,\tau)$
$\displaystyle=3\phi(p\tau)\phi(q\tau)+\frac{\differential\phi(p\tau)}{\differential\ln
p\tau}\frac{\differential\phi(q\tau)}{\differential\ln
q\tau}+\pqty{\phi(p\tau)\frac{\differential\phi(q\tau)}{\differential\ln
q\tau}+\frac{\differential\phi(p\tau)}{\differential\ln p\tau}\phi(q\tau)},$
(2.6)

where $\phi(x)$ is the linear transfer function between $\Phi$ in the
Newtonian gauge and $\zeta$

$\displaystyle\Phi(\tau,\mathbf{k})=\phi(k\tau)\zeta(\mathbf{k})\,,$ (2.7)

with222Here we consider the adiabatic scalar perturbations. In the case of the
isocurvature perturbations, see e.g. Ref. [48] about the transfer function.

$\displaystyle\phi(x)=-\frac{2}{3}\frac{9}{x^{2}}\pqty{\frac{\sin(x/\sqrt{3})}{x/\sqrt{3}}-\cos(x/\sqrt{3})}\,,$
(2.8)

in the RD era. The EoM can be solved by the Green’s function method and the
particular solution is obtained as

$\displaystyle
h_{\lambda}(\tau,\mathbf{k})=\frac{4}{a(\tau)}\int^{\tau}\differential{\tilde{\tau}}G_{\mathbf{k}}(\tau,\tilde{\tau})a(\tilde{\tau})S_{\lambda}(\tilde{\tau},\mathbf{k}),$
(2.9)

where the Green’s function $G_{\mathbf{k}}(\tau,\tilde{\tau})$ is given by

$\displaystyle G_{\mathbf{k}}(\tau,\tilde{\tau})=\frac{\sin
k(\tau-\tilde{\tau})}{k}\,.$ (2.10)

Substituting Eq. (2.4) into Eq. (2.9), the two-point function of induced GWs
is

$\braket{h_{\lambda_{1}}(\tau,\mathbf{k}_{1})h_{\lambda_{2}}(\tau,\mathbf{k}_{2})}=\int\frac{\differential^{3}q_{1}}{(2\pi)^{3}}\int\frac{\differential^{3}q_{2}}{(2\pi)^{3}}Q_{\lambda_{1}}(\mathbf{k}_{1},\mathbf{q}_{1})Q_{\lambda_{2}}(\mathbf{k}_{2},\mathbf{q}_{2})\\\
\times
I_{k}(\absolutevalue{\mathbf{k}_{1}-\mathbf{q}_{1}},q_{1},\tau)I_{k}(\absolutevalue{\mathbf{k}_{2}-\mathbf{q}_{2}},q_{2},\tau)\braket{\zeta(\mathbf{q}_{1})\zeta(\mathbf{k}_{1}-\mathbf{q}_{1})\zeta(\mathbf{q}_{2})\zeta(\mathbf{k}_{2}-\mathbf{q}_{2})},$
(2.11)

with the kernel function

$\displaystyle
I_{k}(p,q,\tau)=4\int^{\tau}\differential{\tilde{\tau}}G_{\mathbf{k}}(\tau,\tilde{\tau})\frac{a(\tilde{\tau})}{a(\tau)}f(p,q,\tilde{\tau}).$
(2.12)

###  2.2 Gravitational waves induced by non-Gaussian curvature perturbations

Let us now consider SIGWs in the presence of non-Gaussianity in the primordial
curvature perturbations. While various types of non-Gaussianities can be
considered depending on inflationary models, in this work, we focus on the
simplest and well-studied type known as quadratic local-type non-Gaussianity
[42, 43, 45]. This type can be described as

$\displaystyle\zeta(\mathbf{x})=\zeta_{g}(\mathbf{x})+F_{\mathrm{NL}}\zeta^{2}_{g}(\mathbf{x}),$
(2.13)

where $F_{\mathrm{NL}}$ is the non-linearity parameter and $\zeta_{g}$ denotes
the Gaussian curvature perturbation. In the following analysis, we assume that
the power spectrum of the Gaussian curvature perturbation is monochromatic.
Once we substitute the above ansatz (Eq. (2.13)) for the non-Gaussian
curvature perturbation into Eq. (2.11), we can obtain seven contributions in
total, which have different momentum integrals [45, 46].

By introducing the tensor (dimensionless) power spectrum
$\mathcal{P}_{\lambda\lambda^{\prime}}(\tau,k)$ which is defined as

$\displaystyle\braket{h_{\lambda}(\tau,\mathbf{k})h_{\lambda^{\prime}}(\tau,\mathbf{k}^{\prime})}=(2\pi)^{3}\delta^{3}(\mathbf{k}+\mathbf{k}^{\prime})\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\lambda\lambda^{\prime}}(\tau,k)\,,$
(2.14)

we can describe

$\displaystyle\mathcal{P}_{\lambda\lambda^{\prime}}(\tau,k)=\sum_{n=1}^{7}\mathcal{P}_{\lambda\lambda}^{(n)}(\tau,k)\delta_{\lambda\lambda^{\prime}}\,,$
(2.15)

and each contribution can be calculated as [46] 333 The coefficients differ
slightly from the expressions in Ref. [46], but the difference is due to the
inclusion of the “deformation” factor.

$\displaystyle\mathcal{P}_{\lambda\lambda}^{(1)}(\tau,k)=2\,\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q},\mathbf{q}\mid\mathbf{q},\mathbf{k}-\mathbf{q}),$
(2.16) $\displaystyle\mathcal{P}^{(2)}_{\lambda\lambda}(\tau,k)$
$\displaystyle=2^{2}\,(2!F_{\mathrm{NL}})^{2}\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{2}\mid\mathbf{q}_{2},\mathbf{k}-\mathbf{q}_{2},\mathbf{q}_{1}-\mathbf{q}_{2}),$
(2.17) $\displaystyle\mathcal{P}^{(3)}_{\lambda\lambda}(\tau,k)$
$\displaystyle=2^{2}\,(2!F_{\mathrm{NL}})^{2}\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{2}\mid\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{k}-\mathbf{q}_{1}-\mathbf{q}_{2}),$
$\displaystyle\mathcal{P}_{\lambda\lambda}^{(4)}(\tau,k)$
$\displaystyle=2^{2}\,\frac{(2!F_{\mathrm{NL}})^{2}}{2!}\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{1}\mid\mathbf{k}-\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{1}-\mathbf{q}_{2}),$
$\displaystyle\mathcal{P}^{(5)}_{\lambda\lambda}(\tau,k)$
$\displaystyle=2\,\frac{(2!F_{\mathrm{NL}})^{4}}{(2!)^{2}}\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{1}\mid\mathbf{q}_{1}-\mathbf{k}+\mathbf{q}_{3},\mathbf{q}_{3},\mathbf{q}_{2},\mathbf{q}_{2}-\mathbf{q}_{1}),$
(2.18) $\displaystyle\mathcal{P}^{(6)}_{\lambda\lambda}(\tau,k)$
$\displaystyle=2\,(2!F_{\mathrm{NL}})^{4}\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{2}\mid\mathbf{q}_{1}-\mathbf{q}_{3},\mathbf{q}_{2}-\mathbf{q}_{3},\mathbf{q}_{3},\mathbf{q}_{3}-\mathbf{k}),$
$\displaystyle\mathcal{P}^{(7)}_{\lambda\lambda}(\tau,k)$
$\displaystyle=(2!F_{\mathrm{NL}})^{4}\mathcal{I}_{\lambda\lambda}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{2}\mid\mathbf{q}_{1}-\mathbf{k}+\mathbf{q}_{2}-\mathbf{q}_{3},\mathbf{q}_{1}-\mathbf{q}_{3},\mathbf{q}_{2}-\mathbf{q}_{3},\mathbf{q}_{3})\,,$

with

$\mathcal{I}_{\lambda\lambda^{\prime}}(\tau,\mathbf{k}\mid\mathbf{q}_{1},\mathbf{q}_{2}\mid\mathbf{k}_{1},\mathbf{k}_{2},\cdots)\\\
\coloneqq\frac{k^{3}}{2\pi^{2}}\int\frac{\differential[3]{q_{1}}}{(2\pi)^{3}}\frac{\differential[3]{q_{2}}}{(2\pi)^{3}}\cdots
Q_{\lambda}(\mathbf{k},\mathbf{q}_{1})Q_{\lambda^{\prime}}(\mathbf{k},\mathbf{q}_{2})I_{k}(\absolutevalue{\mathbf{k}-\mathbf{q}_{1}},q_{1},\tau)I_{k}(\absolutevalue{\mathbf{k}-\mathbf{q}_{2}},q_{2},\tau)\\\
\times P_{g}(k_{1})P_{g}(k_{2})\cdots.$ (2.19)

Here, $P_{g}(k)$ is the power spectrum of the Gaussian curvature perturbation
$\zeta_{g}$ defined as

$\displaystyle\braket{\zeta_{g}(\mathbf{k})\zeta_{g}(\mathbf{k}^{\prime})}=(2\pi)^{3}\delta^{3}(\mathbf{k}+\mathbf{k}^{\prime})P_{g}(k),$
(2.20)

and for the monochromatic case, it can be described as

$\displaystyle
P_{g}(k)=\frac{2\pi^{2}}{k^{3}}A_{g}\delta(\ln{k}-\ln{k_{*}})\,,$ (2.21)

where $A_{g}$ is the power spectrum amplitude and $k_{*}$ is the peak scale.
Note that the constant parameters $A_{g}$ and $F_{\mathrm{NL}}$ can be
separated out from the integrals in $\mathcal{I}_{\lambda\lambda^{\prime}}$.
Then, we can easily see the order of each contribution as
$\mathcal{P}_{\lambda\lambda}^{(1)}=O(F_{\mathrm{NL}}^{0}A_{g}^{2})$,
$\mathcal{P}_{\lambda\lambda}^{(2,3,4)}=O(F_{\mathrm{NL}}^{2}A_{g}^{3})$, and
$\mathcal{P}_{\lambda\lambda}^{(5,6,7)}=O(F_{\mathrm{NL}}^{4}A_{g}^{4})$.

Figure 1: GW energy spectrum for different $F_{\mathrm{NL}}$ values. We set
$A_{g}=0.001$ as a reference value.

To investigate the observational constraint on the SGWBs, the density
parameter per logarithmic wave number
$\Omega_{\mathrm{GW}}\equiv\frac{1}{\rho_{\mathrm{c}}}\frac{\differential{\rho_{\mathrm{GW}}}}{\differential{\log{k}}}$
is commonly used, where $\rho_{\mathrm{c}}$ and $\rho_{\mathrm{GW}}$
respectively denote the critical energy density and the energy density of GW.
For the SIGWs, it is given by

$\displaystyle\Omega_{\mathrm{GW}}(\tau,k)=\frac{1}{48}\left(\frac{k}{\mathcal{H}(\tau)}\right)^{2}\sum_{\lambda,\lambda^{\prime}=+,\times}\overline{\mathcal{P}_{\lambda\lambda^{\prime}}(\tau,k)}\,,$
(2.22)

where the overline means the oscillation average. In the RD era, the source
term of GWs given by Eq. (2.4) becomes negligible soon after the curvature
perturbation re-enters the horizon, and the SIGWs behave as the radiation
without the source. We denote a constant GW density parameter in a deep
subhorizon limit during the RD era as $\Omega_{\mathrm{GW}}^{\mathrm{RD}}(k)$,
and then the present density parameter can be described as (see, e.g., Ref.
[49])

$\displaystyle\Omega_{\mathrm{GW}}(\tau_{0},k)=\Omega_{\mathrm{r},0}\frac{g^{*}(\tau_{i})}{g^{*}_{,0}}\left(\frac{g^{*}_{s,0}}{g^{*}_{s}(\tau_{i})}\right)^{4/3}\Omega_{\mathrm{GW}}^{\mathrm{RD}}(k)\,,$
(2.23)

where the subscript “$0$” denotes the value at present, $\Omega_{\mathrm{r}}$
is the density parameter of the radiation component, $g^{*}$ and $g^{*}_{s}$
respectively represent the effective number of relativistic degrees of freedom
contributing to the energy and entropy densities, and $\tau_{i}$ is the
conformal time when the SIGWs start to behave as radiation. Here, we assume
$g^{*}(\tau_{i})=g^{*}_{s}(\tau_{i})=106.75$, which corresponds to the epoch
in which the relevant scales for LVK band re-entered the horizon. For more
detailed calculations on $\Omega_{\mathrm{GW}}^{\mathrm{RD}}(k)$, one can
refer to Refs. [45, 46]. Finally, in Fig. 1, we plot the GW energy spectrum,
Eq. (2.23), for different $F_{\mathrm{NL}}$ values. The relation between the
gravitational waves frequency $f$ and the scale $k$ is given by

$\displaystyle f=\frac{k}{2\pi}=25\left(\frac{k}{1.67\times 10^{16}~{}\rm
Mpc^{-1}}\right)\rm Hz.$ (2.24)

##  3 Observational constraint by LVK detectors

In this section, we begin with a brief review of the detection method for
SGWBs and then outline our parameter estimation setup.

###  3.1 Basic formulas

Let us consider the time series strain data of detector $I$ positioned at
$\bm{x}_{I}$; $s_{I}(t)=H_{I}(t)+n_{I}(t)$, where $H_{I}(t)$ and $n_{I}(t)$
denote the GW signal measured at the detector and the intrinsic noise of the
detector, respectively. The time series strain data can be transformed in the
Fourier space as

$\displaystyle\tilde{s}_{I}(f)=\tilde{H}_{I}(f)+\tilde{n}_{I}(f)\,.$ (3.1)

The measured signal $\tilde{H}_{I}(f)$ is described by using the Fourier
amplitude of the GW signal $\tilde{h}_{\lambda}(f,\hat{\Omega})$ and the
antenna pattern function $F^{\lambda}_{I}(f,\hat{\Omega})$, which
characterizes the response to a plane wave coming from the sky position
$\hat{\Omega}$, as

$\displaystyle\tilde{H}_{I}(f)=\sum_{\lambda}\int\differential^{2}{\hat{\Omega}}F^{\lambda}_{I}(f,\hat{\Omega})\tilde{h}_{\lambda}(f,\hat{\Omega})e^{-i2\pi
f\bm{x}_{I}\cdot\hat{\Omega}}\,.$ (3.2)

To extract the GW signal, the strain data are cross-correlated with a detector
pair denoted as $I$ and $J$

$\displaystyle\braket{\tilde{s}^{*}_{I}(f)\tilde{s}_{J}(f^{\prime})}=\braket{\tilde{H}^{*}_{I}(f)\tilde{H}_{J}(f^{\prime})}+\braket{\tilde{H}^{*}_{I}(f)\tilde{n}_{J}(f^{\prime})}+\braket{\tilde{n}^{*}_{I}(f)\tilde{H}_{J}(f^{\prime})}+\braket{\tilde{n}^{*}_{I}(f)\tilde{n}_{J}(f^{\prime})}.$
(3.3)

Assuming uncorrelated noise and no correlation between the GWs signal and the
detector noise, the noise terms get suppressed 444If there is a correlated
noise, the last term does not vanish. Schumann Resonances is one of the
sources that induces such a correlated noise [50, 51, 52, 53]. So far, the O3
data shows no evidence of the correlated magnetic noise [54], and many
statistical checks guarantee the stability and Gaussian nature of the noise,
which is assumed in this paper. while the first term remains and can be
enhanced through the accumulation of observation time [55]. Assuming an
isotropic, Gaussian, stationary, and unpolarized SGWB signal, the spectral
density $S_{h}(f)$ can be defined by the ensemble average of the Fourier
amplitude as 555 In Sec. 3, we use the definition of the polarization tensor,
$\mathrm{e}^{\lambda}_{ij}\mathrm{e}^{\lambda^{\prime},ij}=2\delta^{\lambda\lambda^{\prime}}$,
following the common notation of the GWs data analysis. Since $\Omega_{\rm
GW}$ is defined to account for the difference in polarization definitions, we
can confidently use $\Omega_{\rm GW}$ for comparing the theory and data, even
in the presence of notation differences between Sec. 2 and Sec. 3.

$\displaystyle\braket{\tilde{h}^{*}_{\lambda}(f,\hat{\Omega})\tilde{h}_{\lambda^{\prime}}(f^{\prime},\hat{\Omega}^{\prime})}=\frac{1}{2}\delta(f-f^{\prime})\frac{5}{8\pi}S_{h}(f)\delta^{(2)}(\hat{\Omega},\hat{\Omega}^{\prime})\delta_{\lambda\lambda^{\prime}}~{},$
(3.4)

where $\delta(f-f^{\prime})$ denotes the finite-time width delta function,
which typically can be treated as the ideal widthless delta function, and
$\delta^{(2)}(\hat{\Omega},\hat{\Omega}^{\prime})=\delta(\psi-\psi^{\prime})\delta(\cos{\theta}-\cos{\theta^{\prime}})$
is a Dirac delta function on the two dimension sphere. Then the first term in
Eq. (3.3) becomes

$\displaystyle\braket{\tilde{H}^{*}_{I}(f)\tilde{H}_{J}(f^{\prime})}=\frac{1}{2}\delta(f-f^{\prime})\gamma_{IJ}(f)S_{h}(f)\,.$
(3.5)

The function $\gamma_{IJ}(f)$ is the overlap reduction function, which is
determined by the relative position and orientation between the detectors and
given by

$\displaystyle\gamma_{IJ}(f)=\frac{5}{8\pi}\sum_{\lambda}\int\differential^{2}{\hat{\Omega}}F^{\lambda}_{I}(f,\hat{\Omega})F^{\lambda}_{J}(f,\hat{\Omega})e^{-i2\pi
f(\bm{x}_{I}-\bm{x}_{J})\cdot\hat{\Omega}},$ (3.6)

The spectral density can be related to the GWs density parameter Eq. (2.22) as

$\displaystyle
S_{h}(f)=\frac{3H_{0}^{2}}{10\pi^{2}}\frac{1}{f^{3}}\Omega_{\mathrm{GW}}(f)\,,$
(3.7)

where $H_{0}$ is the Hubble parameter.

###  3.2 Parameter estimation

We use the LVK O3 data [56] and perform the multi baseline study summing the
corresponding log-likelihoods for individual pairs of detectors with the
Python package pygwb [57, 58]. The optimal cross-correlation estimator between
times $t$ and $t+T$ is computed from data of detectors $I$ and $J$ as [59, 51,
54]

$\displaystyle
C_{IJ}(f;t)=\frac{2}{T}\frac{\mathrm{Re}\left[\tilde{s}^{*}_{I}(f;t)\tilde{s}_{J}(f;t)\right]}{\gamma_{IJ}(f)S_{0}(f)},$
(3.8)

where $S_{0}(f)=3H_{0}^{2}/(10\pi^{2}f^{3})$ is a factor that converts the GW
strain power spectrum into the fractional energy density. The variance of
$C_{IJ}$ can be written in terms of the one-sided power spectral density of
detector $I$, $P_{I}(f;t)$, and the frequency resolution $\Delta f$ as

$\displaystyle\sigma^{2}_{IJ}(f;t)=\frac{1}{2\Delta
fT}\frac{P_{I}(f;t)P_{J}(f;t)}{\gamma^{2}_{IJ}(f)S_{0}^{2}(f)}.$ (3.9)

In practice, the data is divided into short time segments of $T=192{\rm s}$,
so the estimator Eq. (3.8) and its variance Eq. (3.9) is evaluated at each
time bin. These segments are optimally combined through the following
procedure

$\displaystyle\hat{C}_{IJ}(f)$
$\displaystyle=\frac{\sum_{K}C_{IJ}(f;t_{K})\sigma^{-2}_{IJ}(f;t_{K})}{\sum_{K}\sigma^{-2}_{IJ}(f;t_{K})},$
(3.10) $\displaystyle\hat{\sigma}^{2}_{IJ}(f)$
$\displaystyle=\left(\sum_{K}\sigma^{-2}_{IJ}(f;t_{K})\right)^{-1},$ (3.11)

where $K$ is the index of the time segments and runs from $1$ to $N$ for an
$N$ set of time segments.

In the following parameter estimation, we assume the Gaussian log-likelihood
for each detector pair defined by [60, 51, 25]

$\displaystyle\ln{p(\hat{C}_{IJ}|\bm{\Theta},\xi)}=-\sum_{i}\frac{\left[\hat{C}_{IJ}(f_{i})-\xi\,\Omega_{\mathrm{GW}}(f_{i},\bm{\Theta})\right]^{2}}{2\hat{\sigma}^{2}_{IJ}},$
(3.12)

where the labels $I$ and $J$ run for different combinations of detectors {H,
L, V} and $\xi$ is a factor that accounts for potential calibration
uncertainties [61, 62, 63]. In pygwb, it is modeled as a normal distribution
centered at 1 with a variance $\epsilon^{2}$, and the more conservative O3b
value is taken here for the entire period of O3: $\epsilon=0.148$, 0.123 and
0.108 for the HL, HV and LV baselines, respectively [64]. In order to obtain a
conservative bound, we include a contribution from the CBC background, which
is characterized by the power law spectrum
$\Omega_{\mathrm{CBC},0}(f)=\Omega_{\mathrm{CBC}}(f/f_{\mathrm{ref}})^{2/3}$
with $f_{\mathrm{ref}}=25~{}\mathrm{Hz}$ [54]. The model
$\Omega_{\mathrm{GW}}(f;\bm{\Theta})=\Omega_{\rm
SIGW}(f)+\Omega_{\mathrm{CBC},0}(f)$ contains the fitting parameters
$\bm{\Theta}=(A_{g},F_{\mathrm{NL}},k_{*},\Omega_{\mathrm{CBC}})$ and we
display the priors used for each variable in Table 1. We have selected prior
ranges for the amplitude parameter $A_{g}$ and the peak scale $k_{*}$ such
that the peak of the GW spectrum falls within the LIGO/Virgo sensitivity band.
The range of $\Omega_{\mathrm{CBC}}$ is consistent with that used in the
previous work [25]. The upper bound on the prior range of $F_{\mathrm{NL}}$ is
chosen so that the tail of the posterior distribution can be observed. From a
theoretical perspective, one might consider $F_{\mathrm{NL}}=10^{4}$ to be
relatively large. However, as the upper bound on $A_{g}$ decreases for larger
values of $F_{\mathrm{NL}}$, we have confirmed that the curvature perturbation
$\zeta$ itself does not exceed $1$. Thus, the assumption of cosmological
perturbation expansion remains valid within the parameter space we are
exploring.

Parameters |  Prior
---|---
$\Omega_{\mathrm{CBC}}$ |  $\text{Log-Uniform}[10^{-10},~{}10^{-7}]$
$A_{g}$ |  $\text{Log-Uniform}[10^{-3.5},~{}10^{1}]$
$F_{\mathrm{NL}}$ |  $\text{Log-Uniform}[10^{-1},~{}10^{4}]$
$k_{*}/\mathrm{Mpc}^{-1}$ |  $\text{Log-Uniform}[10^{15.5},~{}10^{18.5}]$
Table 1: Prior distributions used for the parameter estimation.
$\Omega_{\mathrm{CBC}}$ is the amplitude of the CBC power spectrum at $25$Hz,
$A_{g}$ is the primordial scalar power spectrum amplitude, $F_{\mathrm{NL}}$
is the non-Gaussian parameter of primordial scalar perturbations, and $k_{*}$
is the peak position of the primordial scalar power spectrum.

##  4 Results

Figure 2: Posterior distributions for the combined SIGWs and CBC search. The
blue and light blue contours correspond to the $68$ and $95$% CL allowed
regions. For reference, we plot the red line corresponding to
$F_{\mathrm{NL}}^{2}A_{g}=1$, which is the boundary where the non-Gaussian
term starts to dominate over the Gaussian contribution in the GW power
spectrum. The orange line corresponds to the asymptotic behavior at large
$F_{\mathrm{NL}}$ values, which is characterized as
$F_{\mathrm{NL}}A_{g}=0.355$.

The resultant posterior distributions with different parameters are shown in
Fig. 2. The orange line in the $(\log_{10}A_{g}-\log_{10}F_{\mathrm{NL}})$
plane represents the asymptotic behavior of $95$% 666We would like to remark
that the $68$ and $95$% CL contours do not correspond to the notion of 1 and
2$\sigma$ regions in 2 dimensions, which is the default output of pygwb
provided through the corner package [65]. This can be a common misconception.
See the corner documentation [66] for a detailed explanation. confidence level
(CL) curve characterized by $F_{\mathrm{NL}}A_{g}=0.355$ at large
$F_{\mathrm{NL}}$ values. In this regime, the curvature perturbations
asymptotically coincide with
$\zeta=\zeta_{g}+F_{\mathrm{NL}}\zeta^{2}_{g}\rightarrow\zeta\sim
F_{\mathrm{NL}}\zeta^{2}_{g}$. The red line corresponds to the boundary where
non-Gaussianities begin to dominate in the SGWB power spectrum expression.
Considering that we had contributions of
$\mathcal{P}_{\lambda\lambda}^{(1)}=O(F_{\mathrm{NL}}^{0}A_{g}^{2})$,
$\mathcal{P}_{\lambda\lambda}^{(2,3,4)}=O(F_{\mathrm{NL}}^{2}A_{g}^{3})$, and
$\mathcal{P}_{\lambda\lambda}^{(5,6,7)}=O(F_{\mathrm{NL}}^{4}A_{g}^{4})$, we
roughly set this boundary as $F^{2}_{NL}A_{g}=1$ by taking the ratio between
the $\mathcal{P}_{\lambda\lambda}^{(1)}$ and
$\mathcal{P}_{\lambda\lambda}^{(5,6,7)}$ terms. Thus, in the region above the
red line, non-Gaussianity gives the non-negligible effect on the GWs power
spectrum. We observe a transition at $F_{\mathrm{NL}}\simeq 3.55$, which
corresponds to the intersection between the red line and the $95$% CL curve.
This is because, for $F_{\mathrm{NL}}>3.55$, we expect the constraint is set
by the $\mathcal{P}_{\lambda\lambda}^{(5,6,7)}$ terms and it should follow the
asymptotic behavior, which is shown by the orange line,
$F_{\mathrm{NL}}A_{g}=0.355$. On the other hand, when $F_{\mathrm{NL}}<3.55$,
the non-Gaussian contribution is comparable or subdominant, and thus the upper
bound on $A_{g}$ has mild dependence on $F_{\mathrm{NL}}$.

From this result, we obtain the marginalized constraint on the primordial
amplitude $A_{g}<0.011$ at a $95\%$ CL near the best sensitive band of LIGO-
Virgo detectors $k_{*}/\mathrm{Mpc}^{-1}\simeq 2.04\times 10^{16}$. However,
we must note that this bound depends on the chosen prior range for
$F_{\mathrm{NL}}$ and $A_{g}$. In general, when dealing with marginalized
constraints derived from posterior distributions that don’t have tails
converging to zero, the constraints show a strong dependency on the selected
prior range. For this reason, in Fig. 3, we present the $95$% CL curves in the
$(A_{g}-k_{*})$ plane while keeping the value of $F_{\mathrm{NL}}$ fixed. From
this figure, it becomes evident that the upper bound on $A_{g}$ decreases for
larger values of $F_{\mathrm{NL}}$ because it enhances the SGWB spectral
amplitude. It is also worth noting that the scale dependence of $A_{g}$ varies
with different values of $F_{\mathrm{NL}}$, and the most stringent exclusions
are obtained at the scales corresponding to LIGO-Virgo’s most sensitive
frequency range.

In Fig. 4, we present the $95\%$ CL curves in the $(A_{g}-F_{\mathrm{NL}})$
plane by fixing the peak scale $k_{*}$. We show three cases where
$k_{*}/\mathrm{Mpc}^{-1}=\\{10^{16},10^{16.5},10^{17}\\}$ and all three curves
exhibit an asymptotic behavior characterized by
$F_{\mathrm{NL}}A_{g}=\text{const.}$ in the limit of large $F_{\mathrm{NL}}$.
The 95% CL upper limits of this quantity are

$\displaystyle
F_{\mathrm{NL}}A_{g}\leq\left\\{\begin{array}[]{ll}0.13&~{}(k_{*}/\mathrm{Mpc}^{-1}=10^{16})\\\
0.09&~{}(k_{*}/\mathrm{Mpc}^{-1}=10^{16.5})\\\
0.10&~{}(k_{*}/\mathrm{Mpc}^{-1}=10^{17}).\end{array}\right.$ (4.4)

The reason why the constraint at $k_{*}/\text{Mpc}^{-1}=10^{16.5}$ is stronger
is simply that it exhibits a peak in the GW spectrum within the most sensitive
range of the LVK detectors. For small $F_{\mathrm{NL}}$, we observe different
behavior between the cases of $k_{*}/\text{Mpc}^{-1}=10^{16}$ and
$k_{*}/\text{Mpc}^{-1}=10^{16.5}$ or $10^{17}$. Below the red dashed line, the
non-Gaussian contribution to the SGWB spectrum is comparable or subdominant.
As depicted in Fig. 1, the non-Gaussian contribution tends to appear at
frequencies higher than $f_{*}$. In the former case, the upper bound on
$A_{g}$ has a mild dependence on $F_{\mathrm{NL}}$ because the non-Gaussian
effect emerges around the most sensitive frequencies of the LVK detectors. On
the other hand, in the latter case, we find no dependence on
$F_{\mathrm{NL}}$, which is because the non-Gaussian effect occurs outside the
sensitive band. In this figure, we have included dotted lines representing the
combinations of $F_{\mathrm{NL}}$ and $A_{g}$ that would yield a primordial
black hole (PBH) abundance sufficient to constitute $100$% of the dark matter,
denoted as $f_{\rm PBH}=1$. We analytically calculated this curve using the
Peak-theory [67, 68] for estimating PBH abundance (for its application to
local-type non-Gaussianity, refer to [40, 41]). We can see that our current
constraints are weaker than those imposed by the overproduction of PBHs.

Figure 3:  The 95% CL upper bounds on the curvature perturbation amplitude
$A_{g}$ plotted as a function of the peak scale $k_{*}$. These bounds are
obtained by fixing $F_{\mathrm{NL}}$, and different colors represent various
values of $F_{\mathrm{NL}}$. The region below each curve corresponds to the
allowed parameter space.  Figure 4:  The $95\%$ CL counters in the
$A_{g}-F_{\mathrm{NL}}$ plane, obtained by fixing the value of the peak scale.
We present three cases with $k_{*}/\mathrm{Mpc}^{-1}=10^{16}$ (blue),
$10^{16.5}$ (black), $10^{17}$ (light blue). Additionally, we include an
orange line denoting the asymptotic behavior for large non-Gaussianity, where
$F_{\mathrm{NL}}A_{g}=\text{const.}$, and a red line representing the boundary
where the non-Gaussian term dominates over the Gaussian contribution in the GW
power spectrum. Dotted lines correspond to analytically derived relations
between $F_{\mathrm{NL}}$ and $A_{g}$, indicating the point at which PBHs make
up 100% of the dark matter. These lines are plotted for different values of
$k_{*}$ but the difference is tiny.

##  5 Discussion and Conclusion

In this paper, we have investigated the constraint on the non-Gaussian
primordial curvature perturbations, using the most recent data from the LVK
O3. By considering the second-order cosmological perturbations, the primordial
curvature perturbations can induce the GWs called SIGWs, and such SIGWs are a
unique probe for studying inflation physics, as the scales probed by GW
interferometer experiments are much smaller than those of the CMB.

Our main result is in Sec. 4 where we performed the parameter estimation for
the combined model of SIGWs and CBC. We put an upper limit on the amplitude of
the primordial curvature power spectrum $A_{g}$ and the strength of non-
Gaussianity $F_{\mathrm{NL}}$. From the margnialized constraint, shown in Fig.
2, the most stringent limit on the amplitude, $A_{g}<0.011$ at the $95$% CL,
was obtained at the peak scale $k_{*}/\mathrm{Mpc}^{-1}\sim 2.04\times
10^{16}$. Additionally, we observed that the influence of non-Gaussianities on
the GW power spectrum becomes non-negligible when $F_{\mathrm{NL}}\simeq
3.55$. Since marginalized constraints inevitably depend on the range of
priors, we have also presented the analysis by fixing either the
$F_{\mathrm{NL}}$ parameter or the peak scale $k_{*}$ (Figs. 3 and 4).
Although the current constraints appear to be weaker than those imposed by the
overproduction of primordial black holes (PBHs), it is anticipated that these
constraints will improve with the third-generation GW detectors such as the
Einstein Telescope (ET) and Cosmic Explorer (CE).

While we have simply assumed the quadratic order local type non-Gaussianity
(the simple $F_{\mathrm{NL}}$ parametrization) to impose the form of non-
Gaussianity, it is well known that non-Gaussian tails can vary significantly
depending on the mechanisms generating large curvature perturbations. One of
the primary theoretical challenges lies in developing a calculation method for
the SIGW spectrum capable of handling arbitrary shapes of non-Gaussianity. We
also need to establish a robust framework for providing valuable observational
constraints for a wide range of primordial models from the data, given the
anticipated increase in observational opportunities for testing SGWBs in the
near future.

## Acknowledgments

We acknowledge the i-LINK 2021 grant LINKA20416 by the Spanish National
Research Council (CSIC), which initiated and facilitated this collaboration.
The authors are grateful to Alba Romero-Rodríguez for a very useful
discussion. R.I. is supported by JST SPRING, Grant Number JPMJSP2125 and the
“Interdisciplinary Frontier Next-Generation Researcher Program of the Tokai
Higher Education and Research System”. S.J. and S.K. acknowledge support from
the research project PID2021-123012NB-C43 and the Spanish Research Agency
(Agencia Estatal de Investigación) through the Grant IFT Centro de Excelencia
Severo Ochoa No CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033.
S.J. is supported by the FPI grant PRE2019-088741 funded by
MCIN/AEI/10.13039/501100011033. S.K. is supported by the Spanish Atracción de
Talento contract no. 2019-T1/TIC-13177 granted by Comunidad de Madrid, the I+D
grant PID2020-118159GA-C42 funded by MCIN/AEI/10.13039/501100011033, the
Consolidación Investigadora 2022 grant CNS2022-135211 of CSIC, and Japan
Society for the Promotion of Science (JSPS) KAKENHI Grant no. 20H01899,
20H05853, and 23H00110. S. Y. is supported by JSPS KAKENHI Grant no. 20K03968,
and 23H00108. The MCMCs have been run in the Hydra HPC cluster at the IFT.
This research has made use of data or software obtained from the Gravitational
Wave Open Science Center (gwosc.org), a service of the LIGO Scientific
Collaboration, the Virgo Collaboration, and KAGRA. This material is based upon
work supported by NSF’s LIGO Laboratory which is a major facility fully funded
by the National Science Foundation, as well as the Science and Technology
Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS),
and the State of Niedersachsen/Germany for support of the construction of
Advanced LIGO and construction and operation of the GEO600 detector.
Additional support for Advanced LIGO was provided by the Australian Research
Council. Virgo is funded, through the European Gravitational Observatory
(EGO), by the French Centre National de Recherche Scientifique (CNRS), the
Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef,
with contributions by institutions from Belgium, Germany, Greece, Hungary,
Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by the
Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan
Society for the Promotion of Science (JSPS) in Japan; National Research
Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia
Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.

## References

  * [1] LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * [2] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8 [2306.16213].
  * [3] J. Antoniadis et al., _The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals_ , 2306.16214.
  * [4] D.J. Reardon et al., _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L6 [2306.16215].
  * [5] H. Xu et al., _Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I_ , _Res. Astron. Astrophys._ 23 (2023) 075024 [2306.16216].
  * [6] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Search for Signals from New Physics_ , _Astrophys. J. Lett._ 951 (2023) L11 [2306.16219].
  * [7] A.A. Starobinsky, _Spectrum of relict gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682.
  * [8] M.S. Turner, _Detectability of inflation produced gravitational waves_ , _Phys. Rev. D_ 55 (1997) R435 [astro-ph/9607066].
  * [9] S.Y. Khlebnikov and I.I. Tkachev, _Relic gravitational waves produced after preheating_ , _Phys. Rev. D_ 56 (1997) 653 [hep-ph/9701423].
  * [10] J. Garcia-Bellido, D.G. Figueroa and A. Sastre, _A Gravitational Wave Background from Reheating after Hybrid Inflation_ , _Phys. Rev. D_ 77 (2008) 043517 [0707.0839].
  * [11] N. Aggarwal et al., _Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies_ , _Living Rev. Rel._ 24 (2021) 4 [2011.12414].
  * [12] LIGO Scientific, Virgo collaboration, _Constraints on cosmic strings using data from the first Advanced LIGO observing run_ , _Phys. Rev. D_ 97 (2018) 102002 [1712.01168].
  * [13] LIGO Scientific, Virgo, KAGRA collaboration, _Constraints on Cosmic Strings Using Data from the Third Advanced LIGO–Virgo Observing Run_ , _Phys. Rev. Lett._ 126 (2021) 241102 [2101.12248].
  * [14] N. van Remortel, K. Janssens and K. Turbang, _Stochastic gravitational wave background: Methods and implications_ , _Prog. Part. Nucl. Phys._ 128 (2023) 104003 [2210.00761].
  * [15] E. Witten, _Cosmic Separation of Phases_ , _Phys. Rev. D_ 30 (1984) 272.
  * [16] A. Kosowsky, M.S. Turner and R. Watkins, _Gravitational radiation from colliding vacuum bubbles_ , _Phys. Rev. D_ 45 (1992) 4514.
  * [17] A. Kosowsky and M.S. Turner, _Gravitational radiation from colliding vacuum bubbles: envelope approximation to many bubble collisions_ , _Phys. Rev. D_ 47 (1993) 4372 [astro-ph/9211004].
  * [18] C. Caprini and D.G. Figueroa, _Cosmological Backgrounds of Gravitational Waves_ , _Class. Quant. Grav._ 35 (2018) 163001 [1801.04268].
  * [19] K.N. Ananda, C. Clarkson and D. Wands, _The Cosmological gravitational wave background from primordial density perturbations_ , _Phys. Rev. D_ 75 (2007) 123518 [gr-qc/0612013].
  * [20] D. Baumann, P.J. Steinhardt, K. Takahashi and K. Ichiki, _Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations_ , _Phys. Rev. D_ 76 (2007) 084019 [hep-th/0703290].
  * [21] R. Saito and J. Yokoyama, _Gravitational wave background as a probe of the primordial black hole abundance_ , _Phys. Rev. Lett._ 102 (2009) 161101 [0812.4339].
  * [22] G. Domènech, _Scalar Induced Gravitational Waves Review_ , _Universe_ 7 (2021) 398 [2109.01398].
  * [23] K. Inomata and T. Nakama, _Gravitational waves induced by scalar perturbations as probes of the small-scale primordial spectrum_ , _Phys. Rev. D_ 99 (2019) 043511 [1812.00674].
  * [24] S.J. Kapadia, K. Lal Pandey, T. Suyama, S. Kandhasamy and P. Ajith, _Search for the Stochastic Gravitational-wave Background Induced by Primordial Curvature Perturbations in LIGO’s Second Observing Run_ , _Astrophys. J. Lett._ 910 (2021) L4 [2009.05514].
  * [25] A. Romero-Rodriguez, M. Martinez, O. Pujolàs, M. Sakellariadou and V. Vaskonen, _Search for a Scalar Induced Stochastic Gravitational Wave Background in the Third LIGO-Virgo Observing Run_ , _Phys. Rev. Lett._ 128 (2022) 051301 [2107.11660].
  * [26] G. Franciolini, A. Iovino, Junior., V. Vaskonen and H. Veermae, _The recent gravitational wave observation by pulsar timing arrays and primordial black holes: the importance of non-gaussianities_ , 2306.17149.
  * [27] L. Liu, Z.-C. Chen and Q.-G. Huang, _Implications for the non-Gaussianity of curvature perturbation from pulsar timing arrays_ , 2307.01102.
  * [28] K.T. Abe and Y. Tada, _Translating nano-Hertz gravitational wave background into primordial perturbations taking account of the cosmological QCD phase transition_ , 2307.01653.
  * [29] Y.-F. Cai, X.-C. He, X. Ma, S.-F. Yan and G.-W. Yuan, _Limits on scalar-induced gravitational waves from the stochastic background by pulsar timing array observations_ , 2306.17822.
  * [30] J. Garcia-Bellido, M. Peloso and C. Unal, _Gravitational waves at interferometer scales and primordial black holes in axion inflation_ , _JCAP_ 12 (2016) 031 [1610.03763].
  * [31] J. Garcia-Bellido, M. Peloso and C. Unal, _Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter_ , _JCAP_ 09 (2017) 013 [1707.02441].
  * [32] Y.-F. Cai, X. Chen, M.H. Namjoo, M. Sasaki, D.-G. Wang and Z. Wang, _Revisiting non-Gaussianity from non-attractor inflation models_ , _JCAP_ 05 (2018) 012 [1712.09998].
  * [33] V. Atal, J. Cid, A. Escrivà and J. Garriga, _PBH in single field inflation: the effect of shape dispersion and non-Gaussianities_ , _JCAP_ 05 (2020) 022 [1908.11357].
  * [34] 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].
  * [35] S. Pi and M. Sasaki, _Primordial Black Hole Formation in Non-Minimal Curvaton Scenario_ , 2112.12680.
  * [36] 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].
  * [37] LISA Cosmology Working Group collaboration, _Cosmology with the Laser Interferometer Space Antenna_ , _Living Rev. Rel._ 26 (2023) 5 [2204.05434].
  * [38] J.M. Ezquiaga, J. García-Bellido and V. Vennin, _Massive Galaxy Clusters Like El Gordo Hint at Primordial Quantum Diffusion_ , _Phys. Rev. Lett._ 130 (2023) 121003 [2207.06317].
  * [39] S. Pi and M. Sasaki, _Logarithmic Duality of the Curvature Perturbation_ , _Phys. Rev. Lett._ 131 (2023) 011002 [2211.13932].
  * [40] C.-M. Yoo, J.-O. Gong and S. Yokoyama, _Abundance of primordial black holes with local non-Gaussianity in peak theory_ , _JCAP_ 09 (2019) 033 [1906.06790].
  * [41] N. Kitajima, Y. Tada, S. Yokoyama and C.-M. Yoo, _Primordial black holes in peak theory with a non-Gaussian tail_ , _JCAP_ 10 (2021) 053 [2109.00791].
  * [42] R.-g. Cai, S. Pi and M. Sasaki, _Gravitational Waves Induced by non-Gaussian Scalar Perturbations_ , _Phys. Rev. Lett._ 122 (2019) 201101 [1810.11000].
  * [43] C. Unal, _Imprints of Primordial Non-Gaussianity on Gravitational Wave Spectrum_ , _Phys. Rev. D_ 99 (2019) 041301 [1811.09151].
  * [44] C. Yuan and Q.-G. Huang, _Gravitational waves induced by the local-type non-Gaussian curvature perturbations_ , _Phys. Lett. B_ 821 (2021) 136606 [2007.10686].
  * [45] P. Adshead, K.D. Lozanov and Z.J. Weiner, _Non-Gaussianity and the induced gravitational wave background_ , _JCAP_ 10 (2021) 080 [2105.01659].
  * [46] K.T. Abe, R. Inui, Y. Tada and S. Yokoyama, _Primordial black holes and gravitational waves induced by exponential-tailed perturbations_ , 2209.13891.
  * [47] K. Kohri and T. Terada, _Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations_ , _Phys. Rev. D_ 97 (2018) 123532 [1804.08577].
  * [48] G. Domènech, S. Passaglia and S. Renaux-Petel, _Gravitational waves from dark matter isocurvature_ , _JCAP_ 03 (2022) 023 [2112.10163].
  * [49] K. Ando, K. Inomata and M. Kawasaki, _Primordial black holes and uncertainties in the choice of the window function_ , _Phys. Rev. D_ 97 (2018) 103528 [1802.06393].
  * [50] E. Thrane, N. Christensen, R.M.S. Schofield and A. Effler, _Correlated noise in networks of gravitational-wave detectors: subtraction and mitigation_ , _Phys. Rev. D_ 90 (2014) 023013 [1406.2367].
  * [51] P.M. Meyers, K. Martinovic, N. Christensen and M. Sakellariadou, _Detecting a stochastic gravitational-wave background in the presence of correlated magnetic noise_ , _Phys. Rev. D_ 102 (2020) 102005 [2008.00789].
  * [52] E. Thrane, N. Christensen and R. Schofield, _Correlated magnetic noise in global networks of gravitational-wave interferometers: observations and implications_ , _Phys. Rev. D_ 87 (2013) 123009 [1303.2613].
  * [53] Y. Himemoto and A. Taruya, _Impact of correlated magnetic noise on the detection of stochastic gravitational waves: Estimation based on a simple analytical model_ , _Phys. Rev. D_ 96 (2017) 022004 [1704.07084].
  * [54] KAGRA, Virgo, LIGO Scientific collaboration, _Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run_ , _Phys. Rev. D_ 104 (2021) 022004 [2101.12130].
  * [55] B. Allen and J.D. Romano, _Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities_ , _Phys. Rev. D_ 59 (1999) 102001 [gr-qc/9710117].
  * [56] KAGRA, VIRGO, LIGO Scientific collaboration, _Open Data from the Third Observing Run of LIGO, Virgo, KAGRA, and GEO_ , _Astrophys. J. Suppl._ 267 (2023) 29 [2302.03676].
  * [57] A.I. Renzini et al., _pygwb: A Python-based Library for Gravitational-wave Background Searches_ , _Astrophys. J._ 952 (2023) 25 [2303.15696].
  * [58] https://pypi.org/project/pygwb.
  * [59] LIGO Scientific, VIRGO collaboration, _Searching for stochastic gravitational waves using data from the two colocated LIGO Hanford detectors_ , _Phys. Rev. D_ 91 (2015) 022003 [1410.6211].
  * [60] V. Mandic, E. Thrane, S. Giampanis and T. Regimbau, _Parameter Estimation in Searches for the Stochastic Gravitational-Wave Background_ , _Phys. Rev. Lett._ 109 (2012) 171102 [1209.3847].
  * [61] L. Sun et al., _Characterization of systematic error in Advanced LIGO calibration_ , _Class. Quant. Grav._ 37 (2020) 225008 [2005.02531].
  * [62] L. Sun et al., _Characterization of systematic error in Advanced LIGO calibration in the second half of O3_ , 2107.00129.
  * [63] J.T. Whelan, E.L. Robinson, J.D. Romano and E.H. Thrane, _Treatment of Calibration Uncertainty in Multi-Baseline Cross-Correlation Searches for Gravitational Waves_ , _J. Phys. Conf. Ser._ 484 (2014) 012027 [1205.3112].
  * [64] “Parameter estimation demo in pygwb documentation.” https://pygwb.docs.ligo.org/pygwb/run_pe.html.
  * [65] D. Foreman-Mackey, _corner.py: Scatterplot matrices in python_ , _The Journal of Open Source Software_ 1 (2016) 24.
  * [66] “Note about sigmas in corner python package documentation.” https://corner.readthedocs.io/en/latest/pages/sigmas/.
  * [67] J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, _The Statistics of Peaks of Gaussian Random Fields_ , _Astrophys. J._ 304 (1986) 15.
  * [68] C.-M. Yoo, T. Harada, J. Garriga and K. Kohri, _Primordial black hole abundance from random Gaussian curvature perturbations and a local density threshold_ , _PTEP_ 2018 (2018) 123E01 [1805.03946].

  *[SIGW]: scalar-induced gravitational wave
  *[CBC]: compact binary coalescence
  *[LIGO]: Laser Interferometer Gravitational Wave Observatory
  *[SGWBs]: stochastic gravitational wave background
  *[CMB]: cosmic microwave background
  *[GWs]: gravitational wave
  *[SIGWs]: scalar-induced gravitational wave
  *[PTA]: pulsar timing array
  *[LVKs]: LIGO-Virgo-KAGRA
  *[PBHs]: primordial black hole
  *[SGWB]: stochastic gravitational wave background
  *[EoM]: equation of motion
  *[RD]: radiation-dominated
  *[GW]: gravitational wave
  *[LVK]: LIGO-Virgo-KAGRA
  *[CL]: confidence level