Scalar-Induced Gravitational Waves

Scalar perturbations generate tensor perturbations at second order even when the primordial tensor mode vanishes. The effect is negligible on CMB-normalized scales but can become important where the primordial curvature spectrum is strongly enhanced. Its most useful feature is also its central limitation: the tensor spectrum is a convolution of the scalar statistics and cosmic transfer, so it is not a one-to-one image of either a curvature spectrum or a primordial-black-hole abundance.

Helpful background. The Einstein Equation supplies Einstein Equations and transverse-traceless projection; Relativistic Matter and Compact Objects supplies fluid stress-energy; Gravitational Radiation and Waveforms supplies tensor power and stochastic-energy conventions; Relativistic Cosmology supplies scalar perturbations and the expansion history; and Partial Differential Equations and Green Functions supplies the retarded solution.

A quadratic scalar source

On a spatially flat FLRW background, use the Chapter 5 Newtonian-gauge convention

$$ds^2=a^2\left[-(1+2\Phi)d\eta^2 +\left((1-2\Psi)\delta_{ij}+h_{ij}\right)dx^i dx^j\right],$$

with $\partial_i h_{ij}=h^i{}_i=0$. For a perfect fluid with negligible scalar anisotropic stress, $\Phi=\Psi$. This page uses the direct decomposition $h_{ij}=\sum_\lambda e_{ij}^{\lambda}h_\lambda$ and the norm-two basis $e_{ij}^{\lambda}e^{\lambda'ij}=2\delta_{\lambda\lambda'}$; extracting a mode therefore means contracting a tensor equation with $e_\lambda^{ij}/2$. Projecting the second-order spatial Einstein equation then gives

$$h_{\lambda,\mathbf k}''+2\mathcal Hh_{\lambda,\mathbf k}' +k^2h_{\lambda,\mathbf k}=S_{\lambda,\mathbf k},$$ $$S_{\lambda,\mathbf k}(\eta) =2\int\frac{d^3q}{(2\pi)^3} e_{\lambda}^{ij}(\widehat{\mathbf k})q_iq_j\, \mathcal F(q,p;\eta)\, \mathcal R_{\mathbf q}\mathcal R_{\mathbf p}, \qquad \mathbf p=\mathbf k-\mathbf q.$$

$\mathcal F$ contains the two linear scalar transfer functions and their velocity terms. The inverse factor $1/2$ is why the source coefficient is $2$, rather than the coefficient $4$ used with unit-normalized polarization tensors. Terms proportional to $\delta_{ij}$ or $k_i$ disappear under the transverse-traceless projection, leaving the angular quadrupole $e^{ij}q_iq_j$. The integrand is symmetric under $\mathbf q\leftrightarrow\mathbf p$ after the integration variable is changed.

For adiabatic radiation with sound speed $c_s=1/\sqrt3$, write

$$\Psi_{\mathbf k}(\eta) =-\frac23\,T_\Psi(k\eta/\sqrt3)\mathcal R_{\mathbf k}, \qquad T_\Psi(y)=3\frac{\sin y-y\cos y}{y^3}.$$

The minus sign follows the site’s $\mathcal R=-\psi+\mathcal H\delta q/(\rho+p)$ convention; sources using the opposite comoving-curvature sign omit it. It cancels from a Gaussian quadratic power calculation but matters when translating odd non-Gaussian statistics. In this era the dimensionless source kernel can be written

$$\mathcal F=\frac49\left[ T_qT_p+\frac12 \left(T_q+\frac{T_q'}{\mathcal H}\right) \left(T_p+\frac{T_p'}{\mathcal H}\right) \right],$$

where $T_q=T_\Psi(q\eta/\sqrt3)$ and a prime differentiates with respect to $\eta$. Domènech derives the action, perfect-fluid source, and Green-function solution in Domènech 2021, §§3.1–3.2, printed pp. 17–21, eqs. (3.1)–(3.27) (Open PDF). His unit-normalized polarization tensors, potential names, and spatial-potential sign have all been translated as one convention package here.

Green function and scalar convolution

Define the tensor operator and its retarded Green function by

$$\mathcal L_k\equiv\partial_\eta^2+2\mathcal H\partial_\eta+k^2, \qquad \mathcal L_k^{(\eta)}G_k(\eta,\widetilde\eta) =\delta(\eta-\widetilde\eta), \qquad G_k(\eta,\widetilde\eta)=0\quad\text{for }\eta<\widetilde\eta.$$

The derivative of $G_k$ has unit jump at $\eta=\widetilde\eta$, which fixes its normalization. In this branch the initial data satisfy $h_{\lambda,\mathbf k}(\eta_i)=0$ and $h'_{\lambda,\mathbf k}(\eta_i)=0$, so causality gives

$$h_{\lambda,\mathbf k}(\eta) =\int_{\eta_i}^{\eta}d\widetilde\eta\, G_k(\eta,\widetilde\eta)S_{\lambda,\mathbf k}(\widetilde\eta).$$

The corresponding two-time tensor covariance is therefore the retarded double integral

$$\begin{aligned} &\left\langle h_{\lambda,\mathbf k}(\eta) h_{\lambda',\mathbf k'}^*(\eta')\right\rangle \\ &\quad=\int_{\eta_i}^{\eta}d\eta_1 \int_{\eta_i}^{\eta'}d\eta_2\, G_k(\eta,\eta_1)G_{k'}^*(\eta',\eta_2) \left\langle S_{\lambda,\mathbf k}(\eta_1) S_{\lambda',\mathbf k'}^*(\eta_2)\right\rangle . \end{aligned}$$

For Gaussian scalar fluctuations, Wick’s theorem reduces this source unequal-time correlator—which is a scalar four-point function—to products of unequal-time scalar two-point functions. Define $v=q/k$, $u=p/k$, $x=k\eta$, and fix the dimensionless time kernel explicitly:

$$\mathcal I(u,v,x) \equiv k^2\int_{\eta_i}^{\eta}d\widetilde\eta\, G_k(\eta,\widetilde\eta)\, \mathcal F(kv,ku;\widetilde\eta).$$

For radiation, $a\propto\eta$ and the retarded Green function is $G_k(\eta,\widetilde\eta)=(\widetilde\eta/\eta) \sin[k(\eta-\widetilde\eta)]/k$ for $\eta\ge\widetilde\eta$. Consequently, with $x_i=k\eta_i$,

$$\mathcal I(u,v,x) =\int_{x_i}^{x}d\widetilde x\, \frac{\widetilde x}{x}\sin(x-\widetilde x)\, \mathcal F(kv,ku;\widetilde x/k).$$

This definition includes the $4/9$ in $\mathcal F$; it must not be multiplied by another curvature-to-potential conversion. The factor $k^2$ in the first form compensates the dimensions of $d\widetilde\eta\,G_k$. Define the induced spectrum by the same tensor contraction as elsewhere in the chapter,

$$\mathcal P_t^{\rm ind}\equiv 2\sum_\lambda\mathcal P_{h_\lambda}^{\rm ind}.$$

For Gaussian scalars it is

$$\overline{\mathcal P_t^{\rm ind}(k,\eta)} =8\int_0^\infty dv \int_{\lvert1-v\rvert}^{1+v}du\, \mathcal K^2(u,v)\, \overline{\mathcal I^2(u,v,x)} \mathcal P_{\mathcal R}(ku)\mathcal P_{\mathcal R}(kv),$$ $$\mathcal K(u,v) =\frac{4v^2-(1-u^2+v^2)^2}{4uv}.$$

After the polarization contraction produces $\mathcal K^2$, the numerical coefficient $8$ contains the squared site source coefficient $2^2=4$ and the two Gaussian Wick pairings. It is not a stand-alone count of polarizations. The projected source coefficient $2$ and the complete contracted-power kernel must be translated together. The overbar is the late-time oscillation average taken after the scalar source has shut off; the ensemble average is already encoded in the power spectrum. It is not a second ensemble average. For example, take $\mathbf k$ along the $z$ direction and let $\theta,\varphi$ be the polar and azimuthal angles of $\mathbf q$. In the norm-two basis, $e_+^{ij}q_iq_j=q^2\sin^2\theta\cos2\varphi$ and $e_\times^{ij}q_iq_j=q^2\sin^2\theta\sin2\varphi$. The azimuthal integral for either squared projection is $\pi$ times $q^4\sin^4\theta$. With the two Wick pairings, each component therefore carries coefficient $2$ in the displayed $(u,v)$ integral. Summing both components and applying the contraction factor two gives $8$.

The triangular domain is momentum conservation. Both the domain and the full integrand are invariant under $u\leftrightarrow v$. Every integration variable, spectrum, and kernel is dimensionless, so $\mathcal P_t^{\rm ind}$ is dimensionless as required. For freely oscillating subhorizon modes, the site normalization then gives

$$\Omega_{\rm GW}(k,\eta) =\frac1{12}\left(\frac{k}{\mathcal H}\right)^2 \overline{\mathcal P_t^{\rm ind}(k,\eta)}.$$

The short-wave condition is essential. A second-order tensor outside the Hubble radius is not yet assigned the local effective wave energy derived in Isaacson 1968, Eq. (4.1), printed p. 1274. Domènech uses unit-normalized tensors. In symbols, $e_{ij}^{\lambda,\mathrm D}=e_{ij}^{\lambda,\mathrm site}/\sqrt2$; equality of the metric perturbation then requires $h_\lambda^{\mathrm D}=\sqrt2h_\lambda^{\mathrm site}$. Hence

$$\sum_\lambda\mathcal P_{h_\lambda}^{\mathrm D} =2\sum_\lambda\mathcal P_{h_\lambda}^{\mathrm site} =\mathcal P_t^{\rm ind}.$$

His projected source scales in the same way, so his source coefficient $4$ maps to the site’s coefficient $2$. With the Wick contractions and the definition of the time kernel translated at the same time, his complete summed-power coefficient $8$ maps unchanged to the site spectrum; that entire factor cannot be attributed only to polarization number or basis norm. The energy coefficient $1/12$ then applies to the contracted site power. Domènech gives the summed-polarization convolution and energy convention in Domènech 2021, §§2.1 and 3.2, printed p. 10 and pp. 20–21, eqs. (2.5)–(2.9) and (3.22)–(3.27) (Open PDF).

Narrow-peak support and acoustic resonance

Consider the idealized spectrum

$$\mathcal P_{\mathcal R}(k) =A_{\mathcal R}\,\delta\!\left(\ln\frac{k}{k_p}\right).$$

Both scalar legs then require $u=v=k_p/k$. The triangle condition becomes

$$0\lt k\lt2k_p.$$

The induced spectrum therefore has compact support above: no pair of scalar momenta of magnitude $k_p$ can add to a tensor momentum exceeding $2k_p$. In radiation, the time kernel resonates when the two scalar acoustic frequencies add to the tensor frequency,

$$c_s(u+v)=1 \quad\Longrightarrow\quad k_{\rm res}=2c_sk_p=\frac{2k_p}{\sqrt3}.$$

An exact delta peak produces an artificially sharp logarithmic resonance. A finite-width scalar spectrum, finite source duration, or nonlinear effects regularize it. Broad and broken-power-law scalar spectra smear the geometry over many $(u,v)$ pairs; a nonstandard equation of state changes both the sound speed and the Green kernel. These support and resonance checks are derived in Domènech 2021, §§4.2 and 5.2, printed pp. 28–30 and 34–37 (Open PDF).

Gaussianity, PBHs, and perturbative control

The Gaussian formula used a disconnected scalar four-point function. For non-Gaussian initial conditions, the connected trispectrum and loop corrections add momentum integrals and can change both amplitude and shape. Perna et al. enumerate local-type contributions through fifth order in the Gaussian seed in 2024, §2.3, printed pp. 8–15, and §§3.1–3.11, printed pp. 15–28 (Open PDF). Their expansion is a specific local model; it is not a proof that all higher orders are negligible for an arbitrary tail distribution.

The PBH connection is complementary rather than deterministic. The induced two-point function probes an integrated scalar four-point function, while PBH formation probes the rare tail of a smoothed density or compaction distribution and a collapse threshold. Two models with the same $\mathcal P_{\mathcal R}$ can therefore predict similar Gaussian induced spectra but exponentially different PBH abundances. Conversely, tuning non-Gaussianity to alter the PBH tail also changes the tensor trispectrum. Any joint claim must use the same window, curvature sign, non-Gaussian probability law, thermal history, and perturbative truncation on both sides; the PBH calculation is developed in Primordial Black Holes.

Second order is controlled only while omitted scalar and tensor corrections, background backreaction, and loop terms remain subdominant. A large formal $A_{\mathcal R}^2$ signal outside that parameter regime is not automatically a reliable prediction.

There is also a gauge boundary. A tensor perturbation is gauge invariant at first order, but a second-order change of coordinates contains products of first-order scalar shifts and can change the transverse-traceless metric piece. The equations here fix Newtonian gauge. Comparing another gauge requires translating its full second-order observable; the short-wave radiative part after the source decays is the relevant comparison, not an arbitrary nonoscillatory metric term. This is especially consequential in dust eras with persistent scalar sources; see Domènech 2021, §7, printed pp. 54–59 (Open PDF).

Results and limitations through 2026-08-19

Formulation or evidence	Result	Scope or limitation
Radiation-era Gaussian kernel	Gives a symmetric two-dimensional convolution with analytic narrow-peak support and resonance	Assumes Gaussian statistics, linear scalar transfer, radiation domination, and perturbative second order
Broad or broken scalar spectra	Produce smoother wave spectra whose slopes retain information about the scalar shape	Inversion is nonunique because the convolution loses phase and mixes scales
Perna et al. 2024 local non-Gaussian calculation	Exhibits connected and disconnected higher-order contributions with distinct shapes	Truncation and the local ansatz define the model class; forecast precision neglecting foregrounds is not an observational result
Dust or other nonstandard eras	Can enhance or reshape production and move resonances	Tensor energy and gauge choices are subtle when the source does not decay; matching and reheating are consequential
NANOGrav 15-year new-physics analysis	Tests an induced-wave interpretation of a measured low-frequency background	Priors, scalar-shape assumptions, and the competing SMBHB population prevent a unique primordial attribution
LVK O1–O4a early-universe search	Reports no evidence for the tested second-order scalar template	A null result in one frequency range is not a model-independent bound on every scalar spectrum

The observational model tests are documented in Afzal et al. 2023, §5.2, printed pp. 16–20, and §7, pp. 40–43 (Open PDF) and LIGO–Virgo–KAGRA 2025, §VII, printed pp. 15–17 (Open PDF). Detector likelihoods and mutable numerical limits remain within Gravitational-Wave Astronomy.

The unresolved inference problem is not whether the quadratic source exists; it is which scalar statistics, expansion history, and perturbative completion should be assigned to a measured spectrum.

What could change the conclusion. The interpretation could change after a stochastic-background detection resolves spectral structure, a non-Gaussian calculation changes the connected trispectrum or convergence domain, or a perturbativity, gauge, or backreaction analysis invalidates a commonly used scalar parameter regime.

Exercises

For the delta spectrum above, test tensor wavenumbers $k=1.5k_p$ and $k=2.1k_p$. Which is allowed by momentum conservation? Where does the radiation-era resonance lie, and why is its ideal divergence not a robust observable?

Solution

The support condition is $0\lt k\lt2k_p$. Thus $1.5k_p$ is allowed, while $2.1k_p$ is exactly absent in the delta-peak model. The acoustic resonance is

$$k_{\rm res}=\frac{2}{\sqrt3}k_p\simeq1.155k_p,$$

which lies inside the allowed interval. Its sharpness assumes an infinitely narrow scalar spectrum and indefinitely coherent linear transfer. Repeating the calculation with a finite-width peak spreads the $(u,v)$ support and regularizes the feature; non-Gaussian and nonlinear corrections can change it again. The support boundary is a momentum-geometry check, whereas the resonance height is formulation-dependent.

References

• Afzal, Adeela, et al. 2023. “The NANOGrav 15-Year Data Set: Search for Signals from New Physics.” Astrophysical Journal Letters 951, L11. DOI.
• Domènech, Guillem. 2021. “Scalar Induced Gravitational Waves Review.” Universe 7, 398. DOI.
• Isaacson, Richard A. 1968. “Gravitational Radiation in the Limit of High Frequency. II. Nonlinear Terms and the Effective Stress Tensor.” Physical Review 166, 1272–1280. DOI.
• LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration. 2025. “Cosmological and High Energy Physics Implications from Gravitational-Wave Background Searches in LIGO–Virgo–KAGRA’s O1–O4a Runs.” arXiv:2510.26848v2. Open preprint.
• Perna, Gabriele, Chiara Testini, Angelo Ricciardone, and Sabino Matarrese. 2024. “Fully Non-Gaussian Scalar-Induced Gravitational Waves.” arXiv:2403.06962v1. Open preprint.

Dated Research assessment. For open questions and recent results, follow the link below. Each assessment states the date through which it considers results and the limits of its conclusions.

• Scalar-Induced Gravitational Waves and PBHs