Filtered-Squared Bispectrum

• Filtered-squared bispectrum is a technique that converts three-point statistical information into a two-point cross-spectrum using quadratic fields from filtered maps.
• It targets specific triangle configurations in the bispectrum, efficiently isolating gravitational, bias-induced, and near-isosceles contributions.
• The method streamlines covariance estimation and computational scaling by leveraging FFT and pseudo-Cℓ pipelines in both 3D large-scale structure and full-sky analyses.

The filtered-squared bispectrum (FSB) is a class of bispectrum estimators that compresses three-point information into a two-point statistic, typically the cross-power spectrum between a quadratic field built from a filtered map and a linear field. In three-dimensional large-scale structure this compression appears as density–quadratic-field cross-spectra, while in projected and full-sky analyses it is the cross-correlation between the square of a field filtered on a range of scales and the original field or a second field. In both settings, the estimator is a weighted integral of the underlying bispectrum, retains direct sensitivity to selected triangle families, and reuses much of the computational and covariance machinery developed for power-spectrum analysis (Schmittfull et al., 2014, Harscouet et al., 2024).

1. Definition and compression principle

The defining operation is to filter a field, form a quadratic composite, and cross-correlate it with a linear field. On the flat sky, for a scalar field $X(\theta)$ with filtered Fourier modes $X_{f,\ell}=W(\ell)X_\ell$, one defines the squared map

$S(\theta)=[X_f(\theta)]^2,$

and the filtered-square cross-spectrum

$$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$

With the bispectrum defined by $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$, this becomes

$$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$

so the FSB is an angle-averaged line integral through bispectrum configuration space (Harscouet et al., 2024).

The same compression structure appears in three-dimensional large-scale structure. If $Q$ is a quadratic functional of the density field $\delta$, the estimator is the cross-spectrum

$$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$

where the superscript $R$ denotes prior smoothing. Its expectation value is a bispectrum-weighted convolution over the kernel defining $X_{f,\ell}=W(\ell)X_\ell$0, with one filtered linear leg and two filtered quadratic legs (Schmittfull et al., 2014).

A central property of the FSB is that it turns a three-variable statistic into a one-variable data vector. This compression is not arbitrary: in the projected formulation, the estimator is most sensitive to triangles with two sides in the filter support and the third at the cross-spectrum multipole, so diagonal FSBs emphasize near-isosceles configurations, while generalized filtered-multiply constructions extend the coverage to neighboring or asymmetric bands (Harscouet et al., 2024). In three-dimensional perturbation theory, the analogous compression can be chosen to isolate the Legendre components of the tree-level bispectrum kernel, making the estimator directly matched to the dominant gravitational and bias-induced bispectrum contributions (Schmittfull et al., 2014).

2. Three-dimensional large-scale-structure formulation

For the matter or galaxy density contrast $X_{f,\ell}=W(\ell)X_\ell$1, with Fourier transform $X_{f,\ell}=W(\ell)X_\ell$2, the leading-order matter bispectrum in standard Eulerian perturbation theory is

$X_{f,\ell}=W(\ell)X_\ell$3

with

$X_{f,\ell}=W(\ell)X_\ell$4

Using the Legendre basis $X_{f,\ell}=W(\ell)X_\ell$5, $X_{f,\ell}=W(\ell)X_\ell$6, and $X_{f,\ell}=W(\ell)X_\ell$7, this kernel decomposes into monopole, dipole, and quadrupole pieces with coefficients $X_{f,\ell}=W(\ell)X_\ell$8, $X_{f,\ell}=W(\ell)X_\ell$9, and $S(\theta)=[X_f(\theta)]^2,$0 (Schmittfull et al., 2014).

The FSB construction then introduces quadratic fields whose kernels match these angular structures. Writing

$S(\theta)=[X_f(\theta)]^2,$1

the three basic quadratic fields are the squared density,

$S(\theta)=[X_f(\theta)]^2,$2

the shift term,

$S(\theta)=[X_f(\theta)]^2,$3

with displacement field $S(\theta)=[X_f(\theta)]^2,$4, and the tidal term,

$S(\theta)=[X_f(\theta)]^2,$5

with

$S(\theta)=[X_f(\theta)]^2,$6

These three operators form a complete basis for the angular dependence of the leading gravitational bispectrum and the leading nonlinear bias terms $S(\theta)=[X_f(\theta)]^2,$7 and $S(\theta)=[X_f(\theta)]^2,$8 (Schmittfull et al., 2014).

To suppress nonlinear small-scale mode coupling, the fields are smoothed first. The standard choice is a Gaussian filter

$S(\theta)=[X_f(\theta)]^2,$9

with $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$0 in the baseline implementation. The quadratic fields are then built from $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$1, so the cross-spectrum $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$2 carries three factors of the smoothing window in ensemble averages (Schmittfull et al., 2014).

This formulation has a direct physical interpretation. The shift cross-spectrum isolates the dipole component and, in the $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$3 limit and in the absence of velocity bias, has no independent nonlinear bias contribution. The squared-density and tidal cross-spectra isolate the monopole and quadrupole components associated with $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$4 and $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$5. In the large-scale limit the three cross-spectra become equal up to trivial window factors, but away from that limit they decorrelate and become complementary (Schmittfull et al., 2014).

3. Near-optimality, bias estimation, and large-scale limits

The three-dimensional FSB was developed as a near-optimal estimator for bias determination. With Eulerian bias expanded to second order,

$$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$6

the matter–matter–halo cross-spectra satisfy

$$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$7

Subtracting $$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$8 times the matter–matter result isolates the nonlinear bias terms,

$$P^{SX}(L)\equiv V_L^{-1}\int_{|\ell|\in L} d^2\ell\,\langle S_\ell X_{-\ell}\rangle.$$9

For halo–halo–halo cross-spectra, the $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$0 shift piece is purely $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$1-dependent, while the $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$2 and $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$3 pieces carry $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$4 and $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$5, so the three cross-spectra provide nearly orthogonal constraints on $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$6 (Schmittfull et al., 2014).

The near-optimality claim follows from a maximum-likelihood derivation. For any separable theoretical bispectrum template $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$7, the weakly non-Gaussian maximum-likelihood estimator for its amplitude reduces to an inverse-variance weighted integral over a single cross-spectrum of a quadratic filtered density with a linearly filtered density,

$$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$8

Applied to the $$\langle X_k X_q X_p\rangle=(2\pi)^2\delta_D(k+q+p)B(k,q,p)$$9 pieces of the tree-level gravitational bispectrum, this yields the three quadratic fields $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$0, $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$1, and $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$2 as matched filters. In practice, replacing exact inverse-variance weighting by Gaussian smoothing and using the theoretical covariance is sufficient to achieve near-optimality on large scales (Schmittfull et al., 2014).

The same structure explains why joint fits to the power spectrum and FSB break the $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$3–$$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$4 degeneracy. At fixed shape, the FSB amplitudes scale with the fluctuation amplitude through the linear power spectra inside the convolution integrals, while the three cross-spectra depend differently on $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$5, $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$6, and $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$7 than the power spectrum alone (Schmittfull et al., 2014).

The large-scale limit makes this decomposition especially transparent. As $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$8, the halo–halo–halo cross-spectra approach a common expression involving $$P^{SX}(L)=\int \frac{d^2\ell}{(2\pi)^2}W(\ell)W(|L-\ell|)B(\ell,|L-\ell|,L),$$9, $Q$0, and $Q$1, but the shift cross-spectrum remains purely proportional to $Q$2 at lowest order, whereas the $Q$3 and $Q$4 channels carry $Q$5 and $Q$6. A common misconception is therefore that the three estimators are redundant because they coincide on very large scales; in fact, their complementarity is realized on perturbative but nonzero $Q$7, where they decorrelate (Schmittfull et al., 2014).

4. Projected and full-sky generalizations

The projected FSB formalism recasts the same idea for maps on the plane or sphere. On the full sky, if $Q$8 are the spherical-harmonic coefficients of a scalar field and $Q$9, then the full-sky FSB is the cross-power between the squared filtered map and the original map,

$\delta$0

Its expectation value is a linear functional of the reduced bispectrum,

$\delta$1

with

$\delta$2

Thus the FSB is a bandpower-weighted compression of the full angular bispectrum (Harscouet et al., 2024).

Filter choice determines which triangles are emphasized. With a single filter used twice, the basic diagonal estimator weights triangles with two sides in the filter band and the third at the cross-spectrum multipole. Top-hat bands are simple and effective for wide large-scale-structure bands; smooth bands reduce ringing; matched filters can target specific shapes. A generalized filtered-multiply variant $\delta$3 targets squeezed shapes more explicitly, and adding products of neighboring bands substantially increases triangle coverage (Harscouet et al., 2024).

Under masking, the estimator is naturally implemented in the pseudo-$\delta$4 framework. The practical pipeline is to filter the observed map in harmonic space, square it in real space, optionally re-mask the squared map, compute the pseudo-spectrum between the squared map and the observed field, and deconvolve the mask with the standard pseudo-$\delta$5 mixing matrix. Mean subtraction of the squared map removes the induced monopole, and a linear term built from Gaussian simulations can reduce variance in anisotropic conditions (Harscouet et al., 2024).

A mathematically sharper full-sky formulation appears in the optimal binned bispectrum literature. In the mask-free ideal limit, the binned estimator is

$\delta$6

with filtered maps $\delta$7 defined by harmonic binning and weighting. When two bins are equal, this becomes the filtered-squared statistic $\delta$8. In the presence of masking, beam convolution, or inpainting, the estimator remains unbiased after Fisher deconvolution, and the linear term is required for optimality at low $\delta$9 (Philcox, 2023).

The projected formalism also extends directly to cross-bispectra of multiple fields. For maps $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$0, $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$1, and $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$2,

$$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$3

compresses the bispectrum $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$4 with the same Wigner-$$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$5 kernel structure. This generalization underlies recent projected galaxy–galaxy–convergence analyses (Harscouet et al., 10 Jul 2025).

5. Covariance, implementation, and validation

A major attraction of the FSB is that its covariance is much simpler than the covariance of the full bispectrum. In three-dimensional large-scale structure, the leading Gaussian covariance of two cross-spectra at the same wavenumber is

$$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$6

which implies the large-scale inverse-variance weight proportional to $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$7. For halo fields, Poisson stochasticity propagates into the cross-spectra, and phenomenological replacements $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$8 and $$\hat P_{\delta,Q}(k)\equiv \langle \delta^R(k)Q_R(-k)\rangle_{\rm angles},$$9 are used to model deviations from Poissonity (Schmittfull et al., 2014).

For projected fields, the disconnected covariance can be treated as the covariance of a cross-power spectrum between a field and its squared filtered counterpart. In the full sky,

$R$0

With masks, the improved Narrow-Kernel Approximation (iNKA) reuses pseudo-$R$1 machinery with measured mode-coupled spectra. The dominant off-diagonal corrections are $R$2 for FSB auto-covariances and $R$3 for FSB–power cross-covariances, and these terms can be built from measured spectra and generalized FSBs in a model-independent manner (Harscouet et al., 2024).

The computational scaling follows the same logic. In the three-dimensional FFT implementation, one filters $R$4, forms $R$5, $R$6, and $R$7, Fourier-transforms these quadratic fields, and measures three cross-spectra with $R$8. This reduces the task to three FFT-based power spectra rather than brute-force triangle counting with nominal cost $R$9 (Schmittfull et al., 2014). On the sphere, each filter requires one forward spherical-harmonic transform, one inverse transform, one squaring step, and one pseudo-$X_{f,\ell}=W(\ell)X_\ell$00 cross-spectrum, so the total cost scales as $X_{f,\ell}=W(\ell)X_\ell$01 for $X_{f,\ell}=W(\ell)X_\ell$02 filters (Harscouet et al., 2024).

Validation results show that this compression is not merely formal. In the original three-dimensional study, with $X_{f,\ell}=W(\ell)X_\ell$03 at $X_{f,\ell}=W(\ell)X_\ell$04, leading-order predictions agree with $X_{f,\ell}=W(\ell)X_\ell$05-body results for matter–matter–matter and matter–matter–halo cross-spectra at the $X_{f,\ell}=W(\ell)X_\ell$06 level up to $X_{f,\ell}=W(\ell)X_\ell$07, while halo–halo–halo requires additional stochasticity corrections over the same range (Schmittfull et al., 2014). In projected analyses, 3D LPT full-sky simulations (100) and $X_{f,\ell}=W(\ell)X_\ell$08-body GLAM lightcones (600) give FSB measurements unbiased within $X_{f,\ell}=W(\ell)X_\ell$09, and 2D LPT full-sky simulations (6000) show that the disconnected iNKA covariance underestimates the variance near filtered scales by $X_{f,\ell}=W(\ell)X_\ell$10, while adding $X_{f,\ell}=W(\ell)X_\ell$11 and $X_{f,\ell}=W(\ell)X_\ell$12 reproduces the full simulation covariance with $X_{f,\ell}=W(\ell)X_\ell$13 means differing from the simulation-based truth by $X_{f,\ell}=W(\ell)X_\ell$14 (Harscouet et al., 2024).

The multi-field projected extension has also been validated on real survey geometry. In the DESI luminous red galaxy and Planck lensing analysis, the covariance combines Gaussian pseudo-$X_{f,\ell}=W(\ell)X_\ell$15 terms with multi-field generalizations of $X_{f,\ell}=W(\ell)X_\ell$16 and $X_{f,\ell}=W(\ell)X_\ell$17, rescaled by effective $X_{f,\ell}=W(\ell)X_\ell$18 factors calibrated with Gaussian simulations (Harscouet et al., 10 Jul 2025).

6. Relations, applications, and limitations

The FSB sits at the intersection of several pre-existing bispectrum compressions. In redshift-space-distortion studies, the skew-spectrum is the cross-power between a quadratic field and the density field; with the quadratic field taken as the square of a filtered density, it is precisely a filtered-squared bispectrum compression. The integrated bispectrum, by contrast, is a position-dependent power spectrum and correlates local small-scale power with the mean overdensity of a subvolume, so it targets the squeezed limit specifically. The line correlation function is different again: it uses phase information rather than amplitude-based quadratic filtering (Regan, 2017, Chiang et al., 2014).

This relation is especially clear in the soft-limit literature. The position-dependent power spectrum implements a filtered–squared-field estimator of the squeezed-limit bispectrum, and in Fourier space the cross-spectrum of a squared small-scale filtered field with a long-wavelength field is

$X_{f,\ell}=W(\ell)X_\ell$19

which is exactly the filtered-squared structure. In this sense, the integrated bispectrum, skew-spectrum, and FSB become equivalent weighted measures of squeezed configurations, differing mainly by the window and normalization conventions (Chiang et al., 2014, Munshi et al., 2016).

In CMB-style analyses, the FSB is closely analogous to KSW-like quadratic estimators: optimal or separable bispectrum estimators reduce to cross-correlating a linearly filtered field with a quadratic transform of the field. The full-sky optimal binned bispectrum estimator makes this correspondence exact, and the filtered-squared specialization $X_{f,\ell}=W(\ell)X_\ell$20 is simply the case with two equal bins (Philcox, 2023).

Applications now span several observables. For the thermal Sunyaev–Zel'dovich effect, a filtered map $X_{f,\ell}=W(\ell)X_\ell$21 and its square $X_{f,\ell}=W(\ell)X_\ell$22 yield either the scalar estimator $X_{f,\ell}=W(\ell)X_\ell$23 or the cross-spectrum $X_{f,\ell}=W(\ell)X_\ell$24, whose expectation value is a Wigner-$X_{f,\ell}=W(\ell)X_\ell$25-weighted sum over the tSZ bispectrum. In that context the bispectrum amplitude scales as $X_{f,\ell}=W(\ell)X_\ell$26, is principally sourced by massive clusters at redshifts around $X_{f,\ell}=W(\ell)X_\ell$27, and is less sensitive to astrophysical uncertainties than the tSZ power spectrum at $X_{f,\ell}=W(\ell)X_\ell$28 (Bhattacharya et al., 2012).

For projected galaxy clustering and CMB lensing tomography, the multi-field FSB has enabled a combined power-spectrum and bispectrum analysis using DESI luminous red galaxies and Planck lensing. In that study the projected galaxy bispectrum is detected at very high significance, above $X_{f,\ell}=W(\ell)X_\ell$29 in all redshift bins, and the galaxy–galaxy–convergence bispectrum is detected above $X_{f,\ell}=W(\ell)X_\ell$30 in the three highest-redshift bins. The combination of $X_{f,\ell}=W(\ell)X_\ell$31 with power spectra improves the $X_{f,\ell}=W(\ell)X_\ell$32 precision by $X_{f,\ell}=W(\ell)X_\ell$33, while giving constraints consistent with the traditional $X_{f,\ell}=W(\ell)X_\ell$34-point combination (Harscouet et al., 10 Jul 2025).

A more specialized extension appears in pulsar timing arrays. There the primary estimator is a filtered cubic statistic selecting a fixed frequency triangle $X_{f,\ell}=W(\ell)X_\ell$35, but in the squeezed limit $X_{f,\ell}=W(\ell)X_\ell$36 it approaches a filtered-squared-times-residual form, correlating a long-wavelength residual with a filtered square of short-wavelength residuals. The paper states that in the regime where the bispectrum is sharply peaked in the squeezed limit and the long mode is nearly constant across the observation, this filtered-squared estimator is equivalent up to normalization to the full bispectrum estimator with $X_{f,\ell}=W(\ell)X_\ell$37 (Tsuneto et al., 2018).

The main limitations are structural rather than conceptual. The basic diagonal FSB does not measure the full bispectrum shape space uniformly; it emphasizes near-isosceles configurations, and full shape recovery requires neighboring-band products or a filtered-multiply generalization. In three-dimensional perturbation theory, the choice of smoothing scale $X_{f,\ell}=W(\ell)X_\ell$38 trades signal against perturbative control: too small $X_{f,\ell}=W(\ell)X_\ell$39 increases signal but strains perturbation theory, while too large $X_{f,\ell}=W(\ell)X_\ell$40 suppresses signal-to-noise. The statement that the shift term has no independent nonlinear bias holds only in the $X_{f,\ell}=W(\ell)X_\ell$41 limit and in the absence of velocity bias. Real-data applications must additionally handle survey windows, redshift-space distortions, spin fields, or highly irregular sampling, and the current projected cosmology analyses therefore adopt conservative scale cuts and simple tree-level bias prescriptions (Schmittfull et al., 2014, Harscouet et al., 2024, Harscouet et al., 10 Jul 2025).

Taken together, these formulations define the filtered-squared bispectrum as a general strategy for compressing bispectrum information into power-spectrum-like observables. Its distinguishing features are separable map-level construction, compatibility with FFT and pseudo-$X_{f,\ell}=W(\ell)X_\ell$42 pipelines, tractable covariance modeling, and the ability to target physically meaningful triangle families without reverting to full triangle-by-triangle estimation (Schmittfull et al., 2014, Harscouet et al., 2024).