Bispectrum: Third-Order Spectral Analysis

• Bispectrum is a third-order spectral statistic that quantifies three-point phase correlations beyond conventional second-order methods.
• It detects nonlinear phase coupling and departures from Gaussianity, thus providing insight into reversibility and causality in signals.
• Applications of the bispectrum span turbulence, cosmological perturbations, and advanced signal processing techniques such as bicoherence analysis.

The bispectrum is a third-order spectral statistic whose precise definition depends on domain but whose central role is stable across applications: it encodes three-point structure that is invisible to second-order spectra. In stationary time-series analysis, the classical bispectrum is the two-dimensional Fourier transform of third-order joint cumulants, and a broader definition replaces absolute summability of cumulants by the requirement that the cumulant sequence be the inverse Fourier transform of an integrable function $B$ (Iglói et al., 2013). In turbulence and signal processing, the bispectrum is the Fourier three-point correlation function $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$ (O'Brien et al., 2022). In cosmology, it is the three-point function in Fourier or harmonic space, often written either for the primordial curvature perturbation, the matter density field, or angular fields on the sphere (Meerburg, 2010, Regan et al., 2014, Coulton et al., 2019). Across these settings, the bispectrum is the lowest-order statistic sensitive to departures from Gaussianity and to nonlinear phase coupling (Coulton et al., 2019, Maccarone et al., 2011).

1. Definitions, normalizations, and domain-specific conventions

For a stationary time series $\{X(t)\}$, the classical bispectrum is

$$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$

The generalized definition introduced by Iglói and Terdik instead specifies the bispectrum as any integrable function $B$ such that

$$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$

which extends the concept to settings where absolute summability of the cumulant sequence is too restrictive, including cases motivated by long memory (Iglói et al., 2013). For linear time series $X(t)=\sum_k c(k)Z(t-k)$ with transfer function $\beta(\omega)=\sum_k c(k)e^{-ik\omega}$, the bispectrum takes the explicit form

$$B(\omega_1,\omega_2)=\frac{\mathrm{cum}_3(Z(0))}{(2\pi)^2}\beta(\omega_1)\beta(\omega_2)\beta(-\omega_1-\omega_2)$$

almost everywhere (Iglói et al., 2013).

In signal-processing applications based on segmented Fourier transforms, the bispectrum is often estimated by an ensemble average over segments,

$$B(k,l)=\frac{1}{K}\sum_{i=0}^{K-1}X_i(k)X_i(l)X_i^*(k+l),$$

and the associated bicoherence is the normalized magnitude

$$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$0

with $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$1 indicating perfect nonlinear phase coupling and $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$2 indicating no consistent phase coupling (Maccarone et al., 2011). An analogous normalization is used in MHD turbulence, where the bicoherence isolates phase correlations independently of overall amplitude (O'Brien et al., 2022).

In cosmology, the bispectrum is commonly defined as the Fourier-space three-point function. For primordial curvature perturbations,

$$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$3

with $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$4 the model-dependent shape function (Meerburg, 2010). For matter perturbations sourced by cosmic strings,

$$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$5

and for large-scale structure more generally the reduced bispectrum

$$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$6

is used to factor out much of the overall power-spectrum dependence (Regan et al., 2014, Gil-Marín et al., 2011). On the sphere, the reduced bispectrum $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$7 appears after factoring out Wigner $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$8 symbols from $$B(\mathbf{k}_1,\mathbf{k}_2)=\langle \tilde\rho(\mathbf{k}_1)\tilde\rho(\mathbf{k}_2)\tilde\rho^*(\mathbf{k}_1+\mathbf{k}_2)\rangle$$9, which enforces rotational invariance and triangle conditions in harmonic space (Coulton et al., 2019).

2. Phase coupling, reversibility, and identifiability

In stationary time-series theory, real-valuedness of the bispectrum imposes strong structural constraints. If a stationary series is linear, has finite third absolute moment, and has almost everywhere positive spectrum, then a bispectrum that is real-valued and not almost everywhere zero implies reversibility; under the same positive-spectrum condition, a real-valued, nonzero bispectrum also rules out a causal linear representation (Iglói et al., 2013). The same paper shows that for linear processes with nonzero skewness and positive spectrum, reversibility in third order implies full reversibility (Iglói et al., 2013). These results connect a frequency-domain condition to symmetry of moving-average coefficients in the time domain.

In astrophysical timing, the bispectrum is used precisely because second-order spectra are degenerate. In observations of GRS 1915+105, nearly identical power spectra were accompanied by sharply different bicoherence morphologies, including “cross,” “hypotenuse,” and “web” patterns; strong bicoherence involving the QPO frequency showed that the QPOs are coupled to broadband noise rather than generated independently (Maccarone et al., 2011). This excludes models in which the variability components are independent and summed linearly, while supporting reservoir and nonlinear-oscillator interpretations at a qualitative level (Maccarone et al., 2011).

In interstellar-turbulence simulations, the bispectrum and bicoherence similarly expose couplings that the power spectrum cannot. The bicoherence identifies the turbulence driving scale in density and column-density fields, shows that the driving scale is phase-coupled to scales present in the turbulent cascade, and reveals enhanced large-scale phase coupling in the presence of an ordered magnetic field relative to the hydrodynamic case (O'Brien et al., 2022). Scrambling Fourier phases destroys the bicoherence signal, confirming that the statistic is genuinely phase sensitive (O'Brien et al., 2022). The paper therefore suggests the bispectrum and bicoherence as tools for searching for non-locality of wave interactions in MHD turbulence (O'Brien et al., 2022).

3. Primordial and inflationary bispectra

In inflationary cosmology, the bispectrum is the principal observable for leading-order primordial non-Gaussianity. Oscillatory bispectra arise in several classes of models: features in the inflaton potential generate sinusoidal modulations such as $\{X(t)\}$0; resonant models generate logarithmic oscillations such as $\{X(t)\}$1; and non-Bunch-Davies initial states produce more complicated rapidly varying shapes whose amplitudes can scale with frequency (Meerburg, 2010). Because standard smooth templates are inefficient for such signals, mode expansion in separable bases was proposed as a reconstruction strategy. Polynomial bases work well for smooth bispectra, whereas Fourier-based bases are much more efficient for resonant shapes and can reduce the number of required modes by a factor of $\{X(t)\}$2 (Meerburg, 2010).

For single-field inflation beyond slow roll, Adshead, Hu, and Miranda derived a Generalized Slow-Roll integral form for the bispectrum whose domain of validity includes models where the background is not slowly varying everywhere, while still preserving the squeezed-limit consistency relation (Adshead et al., 2013). Their representation reduces bispectrum evaluation for all triangle configurations to combinations of one-dimensional integrals, which expedites comparison of slow-roll-violating models with data (Adshead et al., 2013). In the example of a sharp step in the warped-brane tension of DBI inflation, the equilateral bispectrum associated with the power-spectrum oscillations favored by Planck is described as both extremely large and highly scale dependent (Adshead et al., 2013). The same calculation shows that the bispectrum distinguishes a step in the warp from a step in the potential in canonical single-field inflation, even when the corresponding power spectra are otherwise indistinguishable (Adshead et al., 2013).

The squeezed or soft limit also motivates integrated observables in late-time structure formation. The position-dependent power spectrum, or integrated bispectrum, measures the correlation between the local mean density and the local small-scale power spectrum and therefore captures squeezed bispectrum information (Munshi et al., 2016). In perturbative and hierarchical descriptions, this connects the bispectrum to skew-spectra, cumulant correlators, and higher-order position-dependent polyspectra, providing a unifying language for soft-limit mode coupling (Munshi et al., 2016).

4. Matter and galaxy bispectra in large-scale structure

In large-scale structure, the bispectrum is the leading non-Gaussian statistic and carries information complementary to the power spectrum (Baldauf et al., 2021). In modified-gravity $\{X(t)\}$3 simulations, the reduced dark-matter bispectrum shows deviations up to $\{X(t)\}$4–$\{X(t)\}$5 relative to $\{X(t)\}$6CDM when the initial power spectrum is matched and the final nonlinear power spectrum differs, with the largest effects occurring for squeezed configurations (Gil-Marín et al., 2011). However, after constructing $\{X(t)\}$7CDM models with nearly identical final nonlinear power spectra, the bispectrum differences shrink to $\{X(t)\}$8, often becoming comparable to residual power-spectrum mismatches (Gil-Marín et al., 2011). This indicates that, on mildly nonlinear scales, the bispectrum depends mainly on the power spectrum and less sensitively on the gravitational signatures of the $\{X(t)\}$9 model, which is why the reduced bispectrum remains useful for breaking galaxy-bias degeneracies even beyond general relativity (Gil-Marín et al., 2011).

For cosmic strings, the matter bispectrum has been computed by two distinct methods: a linear unequal-time-correlator approach based on a Gaussian string model, and a nonlinear wake model based on sheet-like overdensities (Regan et al., 2014). Agreement between the two requires compensation factors in both descriptions, and with that inclusion the two approaches give qualitatively and quantitatively consistent results (Regan et al., 2014). Even so, the string-induced matter bispectrum is found to be orders of magnitude below the standard gravitational bispectrum for string tensions compatible with current limits, implying that the matter bispectrum is unlikely to produce competitive constraints on a cosmic-string population (Regan et al., 2014).

Perturbative modeling has also advanced well beyond tree level. The two-loop EFT bispectrum uses kernels up to $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$0, shows that the four independent second-order EFT operators known from the one-loop bispectrum absorb the UV sensitivity of the double-hard region, and introduces one additional EFT parameter for the simplified treatment of the single-hard region (Baldauf et al., 2021). Relative to one-loop EFT, the range of wavenumbers with percent-level agreement with simulations, independently of triangle shape, extends from $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$1 to $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$2 at $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$3, and the two-loop terms begin to matter already at $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$4 (Baldauf et al., 2021).

A complementary route uses Lagrangian perturbation theory. In the Zeldovich approximation, linear displacements are captured to all orders in a manifestly infrared-safe way, and the resulting bispectrum calculation does not require an artificial split of the power spectrum into smooth and oscillatory components (Chen et al., 2024). The Eulerian prescription in which oscillatory components are Gaussian-damped appears as a saddle-point approximation of the LPT calculation, and at one loop the two IR-resummation schemes are in excellent agreement for the bispectrum; at tree level, however, resummed Eulerian perturbation theory captures the nonlinear damping of oscillations less well (Chen et al., 2024).

5. Estimation, compression, covariance, and survey geometry

Because the bispectrum depends on triangle size, shape, and in redshift space orientation, estimator design is central. An inventory of redshift-space bispectrum estimators distinguishes the full FFT-based estimator, modal decompositions, and compressed proxies such as the skew-spectrum, line correlation function, and integrated bispectrum (Regan, 2017). The standard FFT-based estimator measures all triangle configurations but yields a very high-dimensional data vector, whereas the modal estimator compresses the signal into a small number of basis coefficients and, in practice, as few as $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$5–$$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$6 modes can describe most of the relevant information (Regan, 2017). The compressed estimators are lower dimensional and useful as robustness checks, though they retain different subsets of the full three-point information (Regan, 2017).

A specific fast 3D estimator parameterizes triangles by $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$7, $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$8, and $$B(\omega_1,\omega_2)=\sum_{t_1=-\infty}^{\infty}\sum_{t_2=-\infty}^{\infty}\mathrm{cum}(X(0),X(t_1),X(t_2))\,e^{-i(t_1\omega_1+t_2\omega_2)}.$$9, thereby separating size and shape dependence (Shaw et al., 2021). Implemented with FFTs, its computational requirement scales as $B$0, and validation on a non-Gaussian test field showed estimated bispectrum values in good agreement with analytic predictions, with $B$1 deviation across much of the triangle-shape parameter space (Shaw et al., 2021). After introducing linear redshift-space distortion, the estimator remained in close agreement with the analytic monopole bispectrum in redshift space (Shaw et al., 2021).

For full-sky surveys, the Spherical Fourier-Bessel basis provides a natural representation because it diagonalizes the Laplacian in spherical coordinates and allows line-of-sight effects such as RSD and GR terms to be handled exactly rather than through plane-parallel expansions (Benabou et al., 2023). The first numerical calculation of the galaxy bispectrum in this basis was enabled by a Legendre decomposition of the redshift-space kernel $B$2, separate treatment of PNG and velocity-divergence terms, an identity reducing an integral over three spherical harmonics connected by a Dirac delta function to a simple sum, and a formalism for convolution with separable window functions (Benabou et al., 2023). The resulting implementation remains computationally challenging, but it establishes a path toward full large-scale information extraction using an SFB bispectrum (Benabou et al., 2023).

Accurate inference also requires covariance modeling. In redshift space, the Gaussian prediction for the halo bispectrum-monopole variance significantly underestimates numerical estimates for squeezed triangles, and a practical model in that limit is

$B$3

which captures both the variance and dominant off-diagonal terms when validated against a large suite of mock catalogs (Salvalaggio et al., 2024). Super-survey modes add supersample covariance as well; for the matter bispectrum this correction is roughly an order of magnitude smaller, relative to the small-scale covariance, than for the matter power spectrum, because the bispectrum’s own non-Gaussian covariance is already larger, while the power-spectrum–bispectrum cross covariance is as important as in the power-spectrum case (Chan et al., 2017). Geometric distortions introduce another layer: analytic Alcock–Paczynski mixing coefficients show that the leading bispectrum effect is a uniform dilation of all three wavevectors, with shape-dependent multipole mixing; the linear approximation is extremely accurate for the monopole, incurs sub-percent inaccuracies for the quadrupole, and fails for the hexadecapole (Khomeriki et al., 2023).

6. Astrophysical probes, inversion, and modern computational extensions

The bispectrum has become an observational diagnostic well beyond large-scale structure. For the thermal Sunyaev–Zel’dovich effect, the bispectrum amplitude scales as $B$4, compared to $B$5 for the power spectrum, and is sourced mainly by massive clusters at redshifts around $B$6 rather than by the lower-mass, higher-redshift systems that complicate the power spectrum (Bhattacharya et al., 2012). At $B$7, this makes the SZ bispectrum less sensitive to astrophysical uncertainties than the SZ power spectrum, and for an SPT-like survey the expected integrated signal-to-noise is $B$8 under the assumptions quoted in the study (Bhattacharya et al., 2012).

For polarized Galactic foregrounds, the bispectrum provides a direct measure of non-Gaussian contamination in CMB maps. Galactic dust has strong temperature bispectra and strong polarized bispectra that peak in squeezed configurations, including prominent parity-odd signals in BTT, BTE, and BEE, whereas after masking bright sources there is no evidence for polarized synchrotron bispectra and no evidence for dust–synchrotron cross bispectra (Coulton et al., 2019). Applied to component-separated Planck maps, the same tools found no evidence of residual foreground bispectra at Planck sensitivity, supporting their use as null tests for future experiments (Coulton et al., 2019).

More recently, bispectrum inversion has become an estimation primitive in inverse problems. In functional multi-target detection with continuous off-grid translations and correlated stationary Gaussian process noise, the signal bispectrum is first estimated from a debiased empirical third-order autocorrelation and then inverted either by a functional frequency-marching recursion or by a Kotlarski-type deconvolution formula (Little et al., 29 May 2026). Both algorithms come with non-asymptotic recovery guarantees for compactly supported signals without bandlimiting assumptions, and numerical experiments demonstrate accurate recovery in low-SNR regimes (Little et al., 29 May 2026).

In geometric deep learning, the $B$9-bispectrum has been developed as a complete invariant of a signal under a group action. A recent PyTorch library implements selective $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$0-bispectra for seven group actions as differentiable modules; for finite groups, selectivity reduces the coefficient count from $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$1 to $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$2, while for spherical $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$3 signals at band-limit $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$4 an augmented selective construction reduces the count from $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$5 to $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$6 coefficients (Mathe et al., 8 May 2026). The implementation achieves near-exact $$\mathrm{cum}(X(0),X(t_1),X(t_2))=\int_0^{2\pi}\int_0^{2\pi}e^{i(t_1\omega_1+t_2\omega_2)}B(\omega_1,\omega_2)\,d\omega_1\,d\omega_2,$$7-invariance with sub-millisecond GPU runtimes at commonly used bandlimits, and benchmark results show that bispectral pooling consistently outperforms alternative invariant pooling layers in the low-data, moderate-capacity regime (Mathe et al., 8 May 2026).

The resulting picture is that “bispectrum” names not a single formula but a family of third-order spectral constructions whose invariant content is tailored to the data geometry and the scientific question. In time-series theory it diagnoses reversibility and noncausality; in astrophysical variability and turbulence it exposes nonlinear phase coupling; in inflation and large-scale structure it is a primary carrier of non-Gaussian information; in survey analysis it drives estimator, covariance, and geometry-aware methodology; and in modern inverse problems and representation learning it functions as an identifiable invariant or a recoverable sufficient statistic under structured transformations (Iglói et al., 2013, Maccarone et al., 2011, Meerburg, 2010, Salvalaggio et al., 2024, Little et al., 29 May 2026, Mathe et al., 8 May 2026).