Scalar-induced gravitational waves

Scalar and tensor perturbations decouple at first order, but not at second: curvature perturbations source gravitational waves through the quadratic terms of the Einstein equations. For a nearly scale-invariant spectrum the effect is negligible, but any model that enhances small-scale power — the inflationary feature scenarios that produce primordial black holes above all — radiates a stochastic background that pulsar-timing arrays and space interferometers can reach. The induced background is the unavoidable counterpart of PBH formation, which makes it both a discovery channel and a constraint.

Implementation

The calculation is the semianalytic Kohri–Terada one: gravitational waves generated during radiation domination, the source integrated analytically against the exact Green's function, leaving a two-dimensional integral of a universal kernel against the primordial spectrum,

\[\Omega_{\rm GW,RD}(k) \;=\; \frac{1}{12}\int_0^\infty \!dt \int_{-1}^{1}\!ds\; \mathcal{K}(t,s)\; \mathcal{P}_\zeta(k u)\,\mathcal{P}_\zeta(k v),\]

with the kernel carrying the resonance structure (a logarithmic amplification where the source oscillates at twice the potential's frequency) evaluated in closed form. Ω_gw_rd performs this integral for any spectrum callable — including an InflatonSpectrum, so the chain from a potential $V(\phi)$ to an observable $\Omega_{\rm GW}(f)$ has no hand-offs. The implementation reproduces the published anchors of the method exactly: the scale-invariant result $\Omega_{\rm GW} = 0.8222\,A_\zeta^2$, the power-law coefficients $Q(n_s)$, and the monochromatic closed form (sigw_monochromatic), against which the numerical kernel converges in the delta-function limit.

sigw_spectrum maps the radiation-era result to today, $\Omega_{\rm GW,0} = \Omega_{\rm GW,RD}\,\Omega_{r,0}$ with an optional $g_*$ correction hook, and returns frequencies in Hz ($f = 1.55\times10^{-15}\,k\,\mathrm{Mpc}$).

The integrated bound

A gravitational-wave background is radiation, and its integral over frequency counts toward $N_{\rm eff}$:

\[\Delta N_{\rm eff} = \frac{8}{7}\left(\frac{11}{4}\right)^{4/3} \frac{1}{\Omega_{\gamma,0}} \int \Omega_{\rm GW,0}\, d\ln k .\]

sigw_ΔNeff evaluates this, and because the BBN network lives in the same package, the bound can be applied self-consistently: feed the result back through cosmology(N_eff = 3.044 + ΔN) and check the helium and deuterium abundances directly, rather than quoting a detached limit. Spectra with large peaks that clear a detector's sensitivity curve but fail this integral constraint are not viable, and checking is one function call.

Reference

Cosmic.SIGWSpectrum — Type

SIGWSpectrum

Induced-GW background evaluated today: wavenumbers k (1/Mpc), frequencies f (Hz), and Ωh² = ΩGW,0(k)h². Built by [`sigwspectrum`](@ref).

source

Cosmic.sigw_C_connected — Method

sigw_C_connected(Pζ, k; f_NL, ksupp, N = 32_000_000, seed = 4242) -> (Ω_C, err)

The connected "C" contribution to the f_NL² induced-GW density during RD (arXiv:2308.07155 eqs. 36–37):

Ω_C = (2/3π)F_NL² ∫du dv du' dv' dφ₁ cos2φ₁ · I(u,v)I(u',v') · uv/(u'²v'²w³) · 𝒫_ζ(u'k)𝒫_ζ(v'k)𝒫_ζ(wk),

with w = w₁₂ = |q−q'|/k (their eq. 38). Here the free momentum q = uk carries no direct spectrum leg, so the spectrum-leg importance sampling used for sigw_Z_connected is not enough; instead the third leg w is drawn ∝ 𝒫ζ and the azimuth φ₁ is solved from its definition (|cosφ₁| ≤ 1 gates the sample), using |dw/dφ₁| = uu'sθsθ'|sinφ₁|/w. Sampling u', v', w ∝ 𝒫ζ and u, v uniformly,

Ω_C = (4/3π)F_NL² A³ · E[cos2φ₁ · I(u,v)I(u',v') · v W_u W_v / (u'²v' w sθ sθ' |sinφ₁|)],

the expectation over all trials (rejected samples count as 0). Returns the estimate and its 1σ Monte-Carlo error. Validated against the monochromatic collapse (independent 2D quadrature) as the spectrum narrows.

source

Cosmic.sigw_Z_connected — Method

sigw_Z_connected(Pζ, k; f_NL, ksupp, N = 8_000_000, seed = 12345) -> (Ω_Z, err)

The connected "Z" contribution to the f_NL² induced-GW density during RD (arXiv:2308.07155 eqs. 29–30):

Ω_Z = (2/3π)F_NL² ∫du dv du' dv' dφ₁ cos2φ₁ · I(u,v)I(u',v') · uu'/(v²v'²w³) · 𝒫_ζ(vk)𝒫_ζ(v'k)𝒫_ζ(wk),

with w = w₀₁₂ = |k−q−q'|/k (their eq. 32) and FNL = 3/5 fNL. This is a genuine five-dimensional integral with an oscillatory cos2φ₁ and no factorisation, so it is evaluated by importance-sampled Monte Carlo: the two spectrum legs v, v' are drawn ∝ 𝒫_ζ (cancelling the peaked factors), giving

Ω_Z = (4/3)F_NL² A² · E[cos2φ₁ · I(u,v)I(u',v') · uu'/(vv'w³) · 𝒫_ζ(wk) · W_u W_u'],

A = ∫𝒫ζ dln k. Returns the estimate and its 1σ Monte-Carlo error; use Threads.nthreads() > 1 for speed. This is one of the two connected fNL² diagrams — the companion "C" term (their eq. 37) is not yet included; its free momentum has no direct spectrum leg and needs an angle-solving sampler.

source

Cosmic.sigw_f4_reducible — Method

sigw_f4_reducible(Pζ, k; f_NL, ksupp, conv = nothing, rtol = 1e-4)

The disconnected ("reducible") fNL⁴ contribution to ΩGW during RD (arXiv:2308.07155 eq. 42): with both external spectrum legs replaced by the phase-space convolution 𝒞,

Ω_reducible = (F_NL⁴/3) ∫du dv I²(u,v)/(u²v²) 𝒞(uk) 𝒞(vk) = F_NL⁴ · Ω_ts[𝒞·𝒞],

FNL = 3/5 fNL. Like the hybrid it shares the Gaussian kernel and is evaluated through the fast (t,s) integrator, reusing the convolution table. This is the disconnected part of the f_NL⁴ correction; the connected "planar" (eq. 44) and "non-planar" (eq. 48) diagrams are not yet included (see the module note).

source

Cosmic.sigw_gnl_factor — Method

sigw_gnl_factor(Pζ; kmin, kmax, g_NL, rtol = 1e-6)

The linear-gNL correction is a pure rescaling of the Gaussian spectrum (arXiv:2308.07155 eq. 50): Ω{gNL}(k) = 12 GNL A · Ωg(k), with GNL = 9/25 gNL and A = ∫dln p 𝒫ζ(p) the variance of the Gaussian curvature. This returns the multiplicative factor 12 GNL A to apply to the Gaussian ΩGW.

source

Cosmic.sigw_gnl_reducible — Method

sigw_gnl_reducible(Pζ; kmin, kmax, g_NL, rtol = 1e-6)

The full reducible (loop) gNL rescaling of the Gaussian spectrum, summing the closed-form disconnected gNL diagrams (arXiv:2308.07155 eqs. 50, 54, 66):

Ω_g_NL,red(k) = (12 G A + 54 G²A² + 108 G³A³) · Ω_g(k),

with G = 9/25 gNL and A = ∫𝒫ζ dln k. Returns that multiplicative factor. The linear term alone is sigw_gnl_factor (what sigw_spectrum_ng uses); the connected gNL² (tri, ring) and gNL³ (ring3) diagrams are not included here.

source

Cosmic.sigw_hybrid — Method

sigw_hybrid(Pζ, k; f_NL, ksupp, conv = nothing, rtol = 1e-4)

The disconnected ("hybrid") fNL² contribution to ΩGW during RD (arXiv:2308.07155 eq. 28):

Ω_hybrid = (2/3)F_NL² ∫du∫dv I²(u,v)/(u²v²) 𝒫_ζ(vk) 𝒞(uk) = 2F_NL² · Ω_ts[𝒞·𝒫_ζ],

with FNL = 3/5 fNL and 𝒞 the phase-space convolution above. Because the term shares the Gaussian kernel, it is evaluated through the same fast (t,s) integrator (_sigw_omega_ts) rather than the stiff (u,v) form. ksupp = (kmin, kmax) is the curvature-spectrum support; pass a prebuilt conv closure (from _sigw_conv_table) to reuse the tabulation across wavenumbers. This is the disconnected part of the fNL² correction; the connected pieces are [`sigwZconnected](@ref) and [sigwCconnected`](@ref) — the three sum to the full fNL² term.

source

Cosmic.sigw_monochromatic — Method

sigw_monochromatic(k̃)

Kohri–Terada eq. (29): the exact closed-form ΩGW,RD(k)/Aζ² for a monochromatic spectrum Pζ = Aζ δ(ln k/k), as a function of k̃ = k/k. The log-resonance at k̃ = 2/√3 and the cutoff at k̃ = 2 are physical. Used as the analytic anchor for the numerical kernel.

source

Cosmic.sigw_omega_g_kernel — Method

sigw_omega_g_kernel(Pζ, k; ksupp, rtol = 1e-4)

The Gaussian induced-GW density during RD written in the (u,v) kernel variables (arXiv:2308.07155 eq. 24), Ω = (1/3)∫du∫dv I²(u,v)/(u²v²) 𝒫ζ(uk)𝒫ζ(vk). The integration is confined to the curvature-spectrum support ksupp = (kmin, kmax), where 𝒫ζ(uk)𝒫ζ(vk) ≠ 0. Numerically identical to Ω_gw_rd; kept as the cross-check that fixes the kernel normalisation the NG terms inherit.

source

Cosmic.sigw_spectrum — Method

sigw_spectrum(Pζ, c::Cosmology; kmin, kmax, nk = 60, g_factor = k -> 1.0)

Today's induced-GW spectrum ΩGW,0(k)h² for primordial spectrum Pζ (callable; pass `k -> primordialpower(c, k)to use the cosmology's own, or an [InflatonSpectrum](@ref) directly).gfactor(k)` is the relative g* correction (defaults to 1; ≈ 0.39 for modes entering above the electroweak scale — supply your own for precision).

source

Cosmic.sigw_spectrum_ng — Method

sigw_spectrum_ng(Pζ, c::Cosmology; kmin, kmax, nk = 40, f_NL = 0, g_NL = 0,
                 connected = false, f4 = false, gnl_loops = false,
                 N = 16_000_000, rtol = 1e-3, g_factor = k -> 1.0)

Today's induced-GW spectrum ΩGW,0(k)h² including the primordial-NG corrections. The default sums the Gaussian piece, the disconnected fNL² ("hybrid") term, and the linear-g_NL rescaling. The flags add the further validated contributions:

connected — the connected fNL² diagrams [`sigwZconnected](@ref) + [sigwCconnected](@ref) (Monte Carlo withN` samples each per k), completing the full fNL² correction;
f4 — the disconnected fNL⁴ [`sigwf4_reducible`](@ref);
gnl_loops — the full loop gNL tower [`sigwgnl_reducible`](@ref) (12GA+54G²A²+108G³A³) in place of the linear rescaling.

kmin/kmax are the curvature-spectrum support. Still not included: the connected higher-order diagrams (fNL⁴ planar/non-planar, gNL² tri/ring, gNL³ ring3 — see the module note). With `fNL = gNL = 0this returns exactly the Gaussian [sigwspectrum`](@ref).

source

Cosmic.sigw_ΔNeff — Method

sigw_ΔNeff(s::SIGWSpectrum, c::Cosmology)

The integrated background as extra relativistic species, ΔNeff = (8/7)(11/4)^{4/3} ∫ ΩGW,0 dln k / Ωγ,0. Check it against the BBN/CMB bound (≲ 0.3) before trusting large peaks — and for a real gate, feed it back into `cosmology(Neff = 3.044 + ΔN_eff)` and rerun BBN.

source

Cosmic.Ω_gw_rd — Method

Ω_gw_rd(Pζ, k; tmax = 500, rtol = 1e-5)

The induced-GW density parameter during radiation domination at wavenumber k (1/Mpc), for the dimensionless primordial spectrum Pζ(k) (any callable). Kohri–Terada eqs. (6)+(17)+(27):

Ω_GW,RD(k) = (1/12) ∫₀^∞ dt ∫₋₁¹ ds  proj(t,s) · x²Ī²(t,s) · P_ζ(ku)P_ζ(kv)

with u = (t+s+1)/2, v = (t−s+1)/2 (the 1/12 = (1/24)·2, the 2 from the (t,s) Jacobian form). The kernel's log singularity and the resonance step both sit at t = √3−1, which is passed to the quadrature as a breakpoint.

source