sigway.kernels¶

The kernel is the time-integrated transfer function \(\overline{I^2}\) — how efficiently a scalar configuration sources gravitational waves in a given expansion era. Swap the kernel to change the cosmology.

Kernel classes¶

Base class¶

sigway.kernels.Kernel ¶

Kernel(norm=None)

Abstract base class for scalar-induced GW kernels \(\overline{I^2}\).

A kernel encapsulates \(\overline{I^2}(t, s, k, \ldots)\): the oscillation-averaged, time-integrated transfer function squared that measures how efficiently a pair of scalar curvature perturbations with momentum ratios \(u = p/k\) and \(v = q/k\) sources gravitational waves at wave-number \(k\) during a particular cosmological expansion era. The energy density spectrum is

\[\Omega_{\rm GW}(k) = \mathcal{N}(k) \int_0^\infty dt \int_{-1}^{1} ds\; \overline{I^2}(t, s, k)\; \mathcal{P}_\zeta(k\,u)\, \mathcal{P}_\zeta(k\,v),\]

where \(\mathcal{N}(k)\) is the normalisation returned by norm and \(\mathcal{P}_\zeta\) is the dimensionless scalar power spectrum.

Concrete subclasses implement overline_Isq (the smooth part of the kernel) and, when the kernel has a narrow resonant feature, overline_Isq_resonant. They also set the class attributes k_dependent, param_names, and resonant_t to tell the integrator how to call them.

Parameters:

Name	Type	Description	Default
`norm`	`str, float, or callable`	Overall normalisation \(\mathcal{N}(k)\) applied to \(\Omega_{\rm GW}\). A string must be one of the preset keys (`'RD'`, `'CT'`, `'bare'`); a float is used as a constant; a callable must accept a wave-number array \(k\) and return the prefactor. The `'RD'` preset gives today's astrophysical value \(c_g\,\Omega_{r,0}/12\); `'CT'` and `'bare'` give the dimensionless \(1/12\). Defaults to the subclass `_default_norm`.	`None`

Attributes:

Name	Type	Description
`k_dependent`	`bool`	`True` if \(\overline{I^2}\) depends explicitly on \(k\), requiring a separate kernel evaluation for each wave-number.
`param_names`	`tuple of str`	Names of any additional physical parameters the kernel requires beyond \((t, s, k)\), e.g. `('etaR',)` for the transition time.
`resonant_t`	`tuple of float`	Fixed \(t\) values at which the kernel has a resonant feature narrow enough to require a dedicated integration slice.
`nonsmooth_params`	`tuple of str`	Kernel parameters whose gradient requires finite differences rather than JAX auto-differentiation.
`norm_spec`	`str, float, or callable`	The normalisation specification as passed at construction.

Raises:

Type	Description
`ValueError`	If norm is a string that is not a recognised preset.

norm ¶

norm(k)

Return the \(\Omega_{\rm GW}\) overall prefactor evaluated at \(k\).

For the 'RD' preset this equals \(c_g\,\Omega_{r,0}/12\), giving \(\Omega_{\rm GW} h^2\) directly in today's radiation background. For 'CT' / 'bare' the prefactor is \(1/12\), yielding a dimensionless result normalised to the radiation density.

Parameters:

Name	Type	Description	Default
`k`	`array - like`	GW wave-number array (Mpc\(^{-1}\)).	required

Returns:

Type	Description
`float or Array`	Normalisation value(s) with the same shape as \(k\).

overline_Isq ¶

overline_Isq(t, s, k, *kparams)

Oscillation-averaged kernel \(\overline{I^2}\) (smooth part).

Returns the main, smoothly varying piece of the transfer function squared, integrated over the bulk of the \((t, s)\) domain. Narrow resonant features at fixed \(t\) values are handled separately by overline_Isq_resonant.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1 \geq 0\).	required
`s`	`Array`	Dimensionless combination \(s = u - v \in [-1, 1]\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`*kparams`		Additional physical parameters listed in `param_names` (e.g. \(\eta_R\) for InstantEMDKernel).	`()`

Returns:

Type	Description
`Array`	\(\overline{I^2}\) values.

Raises:

Type	Description
`NotImplementedError`	Must be implemented by subclasses.

overline_Isq_resonant ¶

overline_Isq_resonant(t, s, k, *kparams)

Oscillation-averaged kernel \(\overline{I^2_{\rm res}}\) at the resonant slice.

Called at each fixed \(t\) value listed in resonant_t, where the integrand has a narrow peak that cannot be resolved by the smooth quadrature grid. Kernels without a resonance (e.g. RadiationKernel) do not need to override this method.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\), pinned to a value in `resonant_t` (e.g. \(\sqrt{3} - 1\) for the sound-horizon resonance).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`*kparams`		Additional physical parameters listed in `param_names`.	`()`

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm res}}\) values.

Raises:

Type	Description
`NotImplementedError`	Must be implemented by subclasses that declare `resonant_t`.

Radiation domination¶

sigway.kernels.RadiationKernel ¶

RadiationKernel(norm=None)

Bases: Kernel

Kernel for scalar-induced GWs produced entirely during radiation domination.

Use this kernel when all relevant scalar modes re-enter the Hubble radius during a standard radiation-dominated era. The underlying physics is that the GW source is active from Hubble re-entry until the present; once averaged over many oscillation cycles, the kernel takes the closed-form expression in Eqs. (4.21)–(4.22) of arXiv:2501.11320. It has two pieces: a smooth logarithmic term that is always present, and a \(\pi^2\) resonant contribution that switches on when \(u + v > \sqrt{3}\) (equivalently \(t > \sqrt{3} - 1\)), where the combined scalar momentum equals the sound horizon.

Because the Green's function for radiation domination is oscillation-averaged analytically, the kernel is independent of \(k\) — all \(k\) values share the same \((t, s)\) integrand. There is no separate resonant slice to integrate.

The default normalisation preset 'RD' gives \(\Omega_{\rm GW} h^2\) in terms of today's radiation density via \(c_g\,\Omega_{r,0}/12\).

Parameters:

Name	Type	Description	Default
`norm`	`str, float, or callable`	Normalisation override; see Kernel. Defaults to `'RD'`.	`None`

Examples:

>>> from sigway.kernels import RadiationKernel
>>> from sigway.spectrum import OmegaGW   # kernel is the 2nd argument

overline_Isq ¶

overline_Isq(t, s, k, *kparams)

Evaluate the radiation-domination kernel \(\overline{I^2_{\rm RD}}\).

Delegates to I_sq_RD, the numerically stable \((t, s)\) form of the oscillation-averaged kernel for pure radiation domination.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)); unused, accepted for a uniform call signature.	required
`*kparams`		Ignored (this kernel has no extra parameters).	`()`

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm RD}}\) values.

Early matter domination¶

sigway.kernels.InstantEMDKernel ¶

InstantEMDKernel(norm=None)

Bases: Kernel

Kernel for scalar-induced GWs produced in an early matter era with a sudden reheating.

Models a universe that starts in an early matter-dominated (EMD) era (e.g. dominated by a pressureless oscillating field) and transitions instantaneously to radiation domination at conformal time \(\eta_R\). Scalar modes that re-enter the Hubble radius before \(\eta_R\) experience enhanced growth inside the horizon during the matter era; after the transition the sourcing of GWs continues in the radiation era.

The oscillation-averaged kernel receives two distinct contributions that must be integrated separately:

Smooth (large-\(V\)) part — the bulk contribution from modes with \(u \sim v \gg 1\) (large \(t\)), evaluated via I_sq_IRD_LV in overline_Isq. It depends on \(k\) and \(\eta_R\) through \(x_R = k\,\eta_R\).
Resonant slice at \(t = \sqrt{3} - 1\) (i.e. \(u + v = 1/c_s = \sqrt{3}\)) — a narrow peak from the sound-horizon resonance, evaluated via I_sq_IRD_res in overline_Isq_resonant and declared in resonant_t.

The cutoff wave-number \(k_{\max}\), which marks the end of the enhanced scalar power spectrum, is a property of the perturbation object (typically a heaviside ScalarPerturbations) rather than of the kernel itself; it sets the upper boundary of the \((t, s)\) integration domain externally. Accordingly, the underlying numeric cores receive kmax=0.0 here and domain clipping is handled by the integrator.

The default normalisation preset 'CT' gives the dimensionless ratio \(\Omega_{\rm GW} / \Omega_r\) via the prefactor \(1/12\).

Parameters:

Name	Type	Description	Default
`norm`	`str, float, or callable`	Normalisation override; see Kernel. Defaults to `'CT'`.	`None`

Attributes:

Name	Type	Description
`k_dependent`	`bool`	Always `True`; the kernel depends on \(k\) through \(x_R = k\,\eta_R\).
`param_names`	`tuple of str`	`('etaR',)` — the conformal time of the EMD → RD transition must be supplied as a kernel parameter when calling `overline_Isq` and `overline_Isq_resonant`.
`resonant_t`	`tuple of float`	`(sqrt(3) - 1,)` — the single resonant integration slice at \(t = \sqrt{3} - 1\).

overline_Isq ¶

overline_Isq(t, s, k, etaR)

Evaluate the smooth (large-\(V\)) EMD → RD kernel \(\overline{I^2_{\rm IRD,LV}}\).

Delegates to I_sq_IRD_LV with kmax=0.0: domain clipping based on \(k_{\max}\) is applied externally by the integrator.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`etaR`	`float`	Conformal time at the EMD → RD transition (Mpc).	required

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm IRD,LV}}\) values.

overline_Isq_resonant ¶

overline_Isq_resonant(t, s, k, etaR)

Evaluate the resonant EMD → RD kernel \(\overline{I^2_{\rm IRD,res}}\).

Called at the resonant slice \(t = \sqrt{3} - 1\) declared in resonant_t, where \(u + v = \sqrt{3} = 1/c_s\). Delegates to I_sq_IRD_res with kmax=0.0.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\), evaluated at \(\sqrt{3} - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`etaR`	`float`	Conformal time at the EMD → RD transition (Mpc).	required

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm IRD,res}}\) values.

Helper functions¶

These are used internally by the kernels and the integrator; most users never call them directly.

Internal-momentum geometry¶

The change of variables \(u=(t+s+1)/2\), \(v=(t-s+1)/2\) and the geometric factor.

sigway.kernels.get_u ¶

get_u(t, s)

Recover the momentum ratio \(u = p/k\) from the integration variables \((t, s)\).

Inverts the change of variables in Eq. (4.19) of arXiv:2501.11320: \(u = (t + s + 1) / 2\), where \(t = u + v - 1\) and \(s = u - v\).

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1 \geq 0\), integrated over \([0, \infty)\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\), integrated over \([-1, 1]\).	required

Returns:

Type	Description
`Array`	The momentum ratio \(u = p/k\), with the same shape as the broadcast of \(t\) and \(s\).

sigway.kernels.get_v ¶

get_v(t, s)

Recover the momentum ratio \(v = q/k\) from the integration variables \((t, s)\).

Inverts the change of variables in Eq. (4.19) of arXiv:2501.11320: \(v = (t - s + 1) / 2\), where \(t = u + v - 1\) and \(s = u - v\).

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1 \geq 0\), integrated over \([0, \infty)\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\), integrated over \([-1, 1]\).	required

Returns:

Type	Description
`Array`	The momentum ratio \(v = q/k\), with the same shape as the broadcast of \(t\) and \(s\).

sigway.kernels.polynomial ¶

polynomial(t, s)

Geometric projection factor for tensor-mode sourcing (Eq. (4.20) of 2501.11320).

This \(k\)-independent factor arises from contracting the GW polarisation tensor with the stress-energy quadrupole of the two scalar modes. In terms of the integration variables \((t, s)\) it reads

\[2\left[\frac{t(2+t)(s^2-1)}{(1-s+t)(1+s+t)}\right]^2.\]

It vanishes when \(s = \pm 1\) (i.e. \(u = 0\) or \(v = 0\), collinear configuration), and is largest near \(u \approx v \approx 1/\sqrt{3}\) (the resonance).

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1 \geq 0\).	required
`s`	`Array`	Dimensionless combination \(s = u - v \in [-1, 1]\).	required

Returns:

Type	Description
`Array`	Polynomial factor values with the same shape as the broadcast of \(t\) and \(s\).

Kernel cores (per era)¶

The closed-form transfer functions each kernel evaluates.

sigway.kernels.I_sq_RD ¶

I_sq_RD(t, s, k)

Oscillation-averaged radiation-domination kernel, expressed in \((t, s)\) variables.

Evaluates \(\overline{I^2_{\rm RD}}\) for a universe in pure radiation domination, using the numerically stable \((t, s)\) form of Eqs. (4.21)–(4.22) of arXiv:2501.11320. The kernel contains two pieces: a smooth logarithmic term (present for all \(t > 0\)) and a \(\pi^2\) resonant piece that switches on when \(u + v > \sqrt{3}\), i.e. when \(1 + t > \sqrt{3}\), via a Heaviside step. The result is independent of \(k\) because the Green's function for radiation domination has been oscillation-averaged analytically.

This is the production kernel used by RadiationKernel.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)); unused, present for a uniform call signature across all kernel cores.	required

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm RD}}\) values, same shape as the broadcast of \(t\) and \(s\).

sigway.kernels.I_sq_RD_uv ¶

I_sq_RD_uv(t, s, k)

Oscillation-averaged radiation-domination kernel, expressed in \((u, v)\) variables.

Evaluates \(\overline{I^2_{\rm RD}}\) for a universe that remains in radiation domination from horizon re-entry to the present, using the \((u, v)\) form of Eqs. (4.21)–(4.22) of arXiv:2501.11320. Because the Green's function for a radiation-dominated universe has been oscillation-averaged analytically, the result is independent of \(k\).

This function is kept for cross-validation against I_sq_RD; at large \(t\) it is numerically less stable than the \((t, s)\) form. All production calculations use I_sq_RD instead.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)); unused, present for a uniform call signature across all kernel cores.	required

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm RD}}\) values, same shape as the broadcast of \(t\) and \(s\).

sigway.kernels.I_sq_MD ¶

I_sq_MD(t, s, k)

Oscillation-averaged kernel for a universe in pure matter domination.

Returns the constant \(\overline{I^2_{\rm MD}} = 18/25\), independent of \(t\), \(s\), and \(k\). This is the analytic result for a universe that stays matter-dominated forever: the matter-era Green's function grows as a power law, and after oscillation-averaging the transfer function saturates at \(18/25\).

This is an unphysical limiting case — a realistic early-matter-dominated era eventually transitions to radiation domination. For that physical scenario use I_sq_IRD_LV and I_sq_IRD_res. This function is retained for reference and is no longer called in production.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\) (unused).	required
`s`	`Array`	Dimensionless combination \(s = u - v\) (unused).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)); unused.	required

Returns:

Type	Description
`Array`	Constant \(18/25\), broadcast to the shape of \(t\).

sigway.kernels.I_sq_IRD_LV ¶

I_sq_IRD_LV(t, s, k, kmax, etaR)

Smooth (large-\(V\)) contribution to the instantaneous EMD → RD kernel.

Evaluates \(\overline{I^2_{\rm IRD,LV}}\), the dominant bulk piece of the oscillation-averaged transfer function for modes that re-enter the Hubble radius during the early matter-dominated era and are then amplified by the abrupt transition to radiation domination at conformal time \(\eta_R\).

This part of the kernel corresponds to the regime \(u \sim v \gg 1\) (large \(t\)), where the mode functions have undergone many oscillations inside the Hubble radius before the transition. The result depends on \(k\) and \(\eta_R\) only through the dimensionless combination \(x_R = k\,\eta_R\), and is proportional to the Si/Ci combination _sici_precomp. The integration domain in \((t, s)\) is bounded by \(k_{\max}\) via \(x_{\max,R} = k_{\max}\,\eta_R\).

Used by InstantEMDKernel.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`kmax`	`float`	Cutoff wave-number of the scalar power spectrum (Mpc\(^{-1}\)); sets the upper boundary of the \((t, s)\) integration domain via \(x_{\max,R} = k_{\max}\,\eta_R\).	required
`etaR`	`float`	Conformal time at the EMD → RD transition (Mpc).	required

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm IRD,LV}}\) values, same shape as the broadcast of \(t\), \(s\), and \(k\).

sigway.kernels.d_I_sq_IRD_LV ¶

d_I_sq_IRD_LV(index, t, s, k, kmax, etaR)

Analytic gradient of the smooth EMD → RD kernel with respect to \(k_{\max}\) or \(\eta_R\).

Returns the partial derivative of I_sq_IRD_LV selected by index. The gradient with respect to \(k_{\max}\) is identically zero inside the kernel body because \(k_{\max}\) only shifts the integration domain boundaries; the corresponding boundary term is handled separately by the integrator. The gradient with respect to \(\eta_R\) follows from the chain rule on \(x_R = k\,\eta_R\).

Parameters:

Name	Type	Description	Default
`index`	`int`	Selects the differentiation target: `0` for \(k_{\max}\), `1` for \(\eta_R\).	required
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`kmax`	`float`	Cutoff wave-number of the scalar power spectrum (Mpc\(^{-1}\)).	required
`etaR`	`float`	Conformal time at the EMD → RD transition (Mpc).	required

Returns:

Type	Description
`Array`	Gradient values, same shape as the broadcast of \(t\), \(s\), and \(k\).

sigway.kernels.I_sq_IRD_res ¶

I_sq_IRD_res(t, s, k, kmax, etaR)

Resonant contribution to the instantaneous EMD → RD kernel.

Evaluates \(\overline{I^2_{\rm IRD,res}}\), the sharply peaked piece of the transitioning kernel that arises when \(u + v = 1/c_s = \sqrt{3}\), i.e. when the combined momentum of the two scalar modes equals the sound horizon at the transition. In integration variables this resonance sits at the fixed slice \(t = \sqrt{3} - 1\), where the integrand is not smooth and must be treated separately.

Near the resonance the \(\mathrm{Ci}\) function diverges logarithmically; at the transition scale \(x_R = k\,\eta_R\) it is approximated by \(\mathrm{Ci}(x_R/2) \approx 7.97727 / x_R\), which captures the dominant behaviour for the \(k\) values of interest.

Used by InstantEMDKernel as the resonant integration slice declared in resonant_t.

Parameters:

Name	Type	Description	Default
`t`	`Array`	Dimensionless combination \(t = u + v - 1\); this function is evaluated at \(t = \sqrt{3} - 1\) to capture the resonance.	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`kmax`	`float`	Cutoff wave-number (Mpc\(^{-1}\)); unused here, present for a uniform call signature.	required
`etaR`	`float`	Conformal time at the EMD → RD transition (Mpc).	required

Returns:

Type	Description
`Array`	\(\overline{I^2_{\rm IRD,res}}\) values, same shape as the broadcast of \(t\), \(s\), and \(k\).

sigway.kernels.d_I_sq_IRD_res ¶

d_I_sq_IRD_res(index, t, s, k, kmax, etaR)

Analytic gradient of the resonant EMD → RD kernel with respect to \(k_{\max}\) or \(\eta_R\).

Returns the partial derivative of I_sq_IRD_res selected by index. The gradient with respect to \(k_{\max}\) is identically zero (the resonant integrand does not depend on \(k_{\max}\) directly). The gradient with respect to \(\eta_R\) follows from the power-law dependence of the kernel on \(x_R = k\,\eta_R\): since the resonant kernel scales as \(x_R^7\), the derivative is \(7 / \eta_R\) times the kernel value.

Parameters:

Name	Type	Description	Default
`index`	`int`	Selects the differentiation target: `0` for \(k_{\max}\), `1` for \(\eta_R\).	required
`t`	`Array`	Dimensionless combination \(t = u + v - 1\).	required
`s`	`Array`	Dimensionless combination \(s = u - v\).	required
`k`	`Array`	GW wave-number (Mpc\(^{-1}\)).	required
`kmax`	`float`	Cutoff wave-number of the scalar power spectrum (Mpc\(^{-1}\)).	required
`etaR`	`float`	Conformal time at the EMD → RD transition (Mpc).	required

Returns:

Type	Description
`Array`	Gradient values, same shape as the broadcast of \(t\), \(s\), and \(k\).