# Parity-Odd Power Spectra: Concise Statistics for Cosmological Parity Violation

Drew Jamieson, Angelo Caravano, Jiamin Hou, Zachary Slepian, Eiichiro Komatsu

► To cite this version:

Drew Jamieson, Angelo Caravano, Jiamin Hou, Zachary Slepian, Eiichiro Komatsu. Parity-Odd Power Spectra: Concise Statistics for Cosmological Parity Violation. Monthly Notices of the Royal Astronomical Society, 2024, 533 (3), pp.2582-2598. 10.1093/mnras/stae1924. hal-04635826

HAL Id: hal-04635826

https://hal.science/hal-04635826v1

Submitted on 21 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Parity-odd power spectra: concise statistics for cosmological parity violation

Drew Jamieson<sup>1,1*</sup> Angelo Caravano,<sup>2,★</sup> Jiamin Hou,<sup>3,4</sup> Zachary Slepian<sup>3</sup> and Eiichiro Komatsu<sup>1,5,6</sup>

Accepted 2024 August 5. Received 2024 August 1; in original form 2024 July 15

## ABSTRACT

We introduce the parity-odd power (POP) spectra, a novel set of observables for probing parity violation in cosmological N-point statistics. POP spectra are derived from composite fields obtained by applying non-linear transformations, involving also gradients, curls, and filtering functions, to a scalar field. This compresses the parity-odd trispectrum into a power spectrum. These new statistics offer several advantages: they are computationally fast to construct, estimating their covariance is less demanding compared to estimating that of the full parity-odd trispectrum, and they are simple to model theoretically. We measure the POP spectra on simulations of a scalar field with a specific parity-odd trispectrum shape. We compare these measurements to semi-analytic theoretical calculations and find agreement. We also explore extensions and generalizations of these parity-odd observables.

**Key words:** methods: analytical – methods: data analysis – methods: statistical – inflation – large-scale structure of Universe.

## 1 INTRODUCTION

Modern physics has advanced considerably by proposing theoretical symmetries and testing whether or not those symmetries are respected or violated in nature. A crucial breakthrough for the standard model of particle physics occurred with the surprising discovery of weak nuclear parity violation in 1957 (Lee & Yang 1956; Garwin, Lederman & Weinrich 1957; Wu et al. 1957). Parity is the negation of spatial coordinates: $x \to -x$, and its violation implies that the laws of physics distinguish between right-handed and left-handed chiralities. While parity is known to be maximally violated on small scales by weak nuclear interactions, the status of parity on cosmological scales remains an open question.

Recent investigations using the method proposed in Cahn, Slepian & Hou (2023) and Cahn & Slepian (2023) have uncovered intriguing evidence of cosmological parity violation in the large-scale structure (LSS) of the Universe. Notably, analyses of the BOSS galaxy four-point correlation function detected parity violation with significance of up to $7\sigma$ (Philcox 2022; Hou, Slepian & Cahn 2023b). Confirmation of parity violation in the galaxy clustering would have profound implications. If the signal is primordial in origin, it is a form of primordial non-Gaussianity and would greatly inform our early-Universe inflationary models (Bartolo et al. 2004). Current cosmic microwave background (CMB) observations are consistent with parity conservation (Philcox 2023; Philcox & Shiraishi 2024). However, since the CMB is measured on a 2D surface and is sensitive to different scales than the LSS, it is still unclear what this implies for the BOSS signal. The central challenge in the current BOSS analysis lies in robustly estimating the covariance (Cahn et al. 2023), which can be influenced by both instrumental systematics and observational effects (Hou et al. 2023b) as well as the type of mocks used to calibrate or compute the covariance (Hou et al. 2022; Philcox & Ereza 2024).

Next-generation 3D spectroscopic surveys such as DESI (Aghamousa et al. 2016; Adame et al. 2024), Euclid (Amendola et al. 2018), Subaru PFS (Takada et al. 2014), Roman (Wang et al. 2022), and SPHEREX (Doré et al. 2014), as well as photometric efforts such as LSST (Ivezić et al. 2019), will provide significantly more data than is currently available. To constrain large-scale parity violation robustly with this wealth of data, we need analysis tools that facilitate accurate covariance estimation and characterization of observational systematics.

For scalar fields like the primordial curvature perturbation or the matter density contrast, N-point statistics are sensitive to parity at orders four and above (Shiraishi 2016). Measurements of these correlators have an enormous number of degrees of freedom [e.g. the BOSS analysis of Hou et al. (2023b) had 18 000]. Covariance estimation for such large data vectors requires a large number of simulated mock data sets, often

<footnote n="1"><sup>1</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, UMR 7095 du CNRS et de Sorbonne Université, 98 bis Bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="4"><sup>4</sup>Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="5"><sup>5</sup>Ludwig-Maximilians-Universität München, Schellingstr. 4, D-80799 München, Germany</footnote>

<footnote n="6"><sup>6</sup>Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), University of Tokyo, Chiba 277-8582, Japan</footnote>

<footnote>*E-mails: jamieson@mpa-garching.mpg.de (DJ); caravano@iap.fr (AC); jiamin.hou@ufl.edu (JH)</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Here we display a Fourier-space tetrahedron representing the trispectrum configuration $T(k_1, k_2, k_3)$. The diagonals (in dashed magenta) are defined as $K = k_1 + k_2$ and $\bar{K} = k_1 + k_4$. The tetrahedron is formed by joining two triangles, $\{k_1, k_2, -K\}$ on the left and $\{k_3, k_4, K\}$ on the right, along their shared edge.</figcaption></figure>

combined with approximate analytic covariance templates [e.g. Hou et al. (2022); also Slepian & Eisenstein (2015) and Xu et al. (2012)]. Even with statistically accurate covariance estimation, systematics in observational data can bias both the detected signal and its significance.

A complementary approach to high-order N-point statistics is lower-order statistics, such as power spectra, that are computed on composite fields. These composite fields are created by applying non-linear transformations to the original scalar field. Composite-field statistics average or compress the high-dimensional correlators down to low-dimensional correlators. As a result, the covariances of these composite-field power spectra are less demanding to estimate, due to having fewer degrees of freedom, while still carrying higher-order information. The trade-off involved in this approach is the loss of information through compression. Examples of weighted, compressed analyses for bispectra and parity-even trispectra have been done before and are known as skew-spectra and kurt-spectra (Schmittfull, Baldauf & Seljak 2015; Moradinezhad Dizgah et al. 2020; Schmittfull & Moradinezhad Dizgah 2021; Munshi et al. 2022; Hou et al. 2023a), as well as works on the unweighted integrated three-point correlation function (Slepian et al. 2017), trispectrum (Gualdi, Gil-Marín & Verde 2021; Gualdi & Verde 2022), and four-point function (Sabiu et al. 2019); the clustering fossil estimators (Jeong & Kamionkowski 2012); and a Fourier-transform-based algorithm for the full four-point function is Sunseri et al. (2023). An estimator of the CMB lensing power spectrum can also be constructed in this manner (Hu 2001).

In the current work, we introduce new methods for detecting cosmological parity violation by constructing power-spectrum-like statistics sensitive to parity. We call our new observables parity-odd power spectra, or POP spectra. The POP spectra are compressions of the 6D four-point statistics down to 1D power spectra. They are computationally efficient to construct and their lower dimensionality alleviates the burden of full four-point covariance estimation.

Importantly, POP spectra are sensitive to soft limits of parity-violating trispectrum shapes, which can help to place strong constraints on primordial parity violation. Moreover, being formulated in Fourier space, they are easy to interpret and can facilitate theoretical understanding of the signal, which is crucial to distinguish it from potential observational systematics. It is also possible to construct real-space and spherical-harmonic-space versions of these statistics. As we will demonstrate below, constructing parity-odd statistics is non-trivial, requiring more sophisticated techniques than those required for parity-even compressed estimators.

This work is structured as follows. In Section 2, we discuss the primordial parity-odd trispectrum. In Section 3, we present two POP spectra constructions based on composite vector and scalar fields. In Section 4, we validate the estimator on simulated data, including comparisons with semi-analytical calculations. In Section 5, we explore extensions and generalizations of our parity-odd observables. Finally, we conclude and discuss prospects for measuring these parity-odd statistics on observational data in Section 6.

## 2 PARITY-VIOLATING TRISPECTRUM

Under statistical homogeneity and isotropy, the four-point correlation function is the lowest-order parity-sensitive statistic for a single scalar field. In Fourier space, this corresponds to the trispectrum, $T(k_1, k_2, k_3)$, defined through the four-point correlator,

$$ \langle \Phi(k_1) \Phi(k_2) \Phi(k_3) \Phi(k_4) \rangle = (2\pi)^3 \delta_D^{(3)}(k_1 + k_2 + k_3 + k_4) T(k_1, k_2, k_3). \quad (1) $$

The Dirac delta function imposes translational invariance by requiring that the sum of the four wave vectors vanishes. Thus, they form a closed loop in Fourier space, and define a tetrahedron, as displayed in Fig. 1. More precisely, the four wave vectors define an equivalence class of tetrahedra, consisting of the 24 tetrahedra produced by permuting them. These are not all geometrically distinct, since a cyclic permutation of wave vectors within a given tetrahedron yields a congruent tetrahedron. The three independent wave vectors, which we take to be $k_1, k_2$, and $k_3$, uniquely determine these tetrahedra. Statistical isotropy, however, reduces the dimensionality of the trispectrum from nine to six continuous degrees of freedom. We will choose the wave vector magnitudes $k_1, k_2, k_3$, and $k_4$ as four of these continuous parameters, and refer to them as sides of the tetrahedron. Then we take $K = |k_1 + k_2| = |k_3 + k_4|$ and $\bar{K} = |k_1 + k_4| = |k_2 + k_3|$ as the other two parameters and refer to these as diagonals to distinguish them from the wave vector magnitudes.
Among the set of 24 tetrahedra, there is an additional discrete degree of freedom that splits them into two groups of 12. These are distinguished by the sign of the vector triple product $k_1 \cdot (k_2 \times k_3)$. This can be thought of as the handedness or helicity of the tetrahedron. The triple product is positive for a right-handed configuration and negative for a left-handed configuration. Under a parity transform, the sign of the triple product changes, so a parity transformation interchanges the two helicities.

We can isolate the trispectrum for only right-handed configurations, and obtain the right-handed trispectrum $T_R(k_1, k_2, k_3)$, and similarly for the left-handed trispectrum $T_L(k_1, k_2, k_3)$. Let $\pi_{123} \equiv \text{sgn}(k_1 \cdot (k_2 \times k_3))$, then

$$T(k_1, k_2, k_3) = \frac{1 + \pi_{123}}{2} T_R(k_1, k_2, k_3) + \frac{1 - \pi_{123}}{2} T_L(k_1, k_2, k_3). \qquad (2)$$

The right- and left-handed parts of the trispectrum are not necessarily equal. If they are different, this indicates that the trispectrum has a parity-odd component. For the modes of a real field, a parity transformation is complex conjugation,

$$\Phi(-k) = \Phi(k)^*, \qquad (3)$$

so the parity-even component of the trispectrum is its real part and is proportional to the sum of $T_R + T_L$. The parity-odd component is its imaginary part and is proportional to the difference $T_R - T_L$:

$$T(k_1, k_2, k_3) = \frac{1}{2} (T_R(k_1, k_2, k_3) + T_L(k_1, k_2, k_3)) + \frac{\pi_{123}}{2} (T_R(k_1, k_2, k_3) - T_L(k_1, k_2, k_3)) \qquad (4)$$

$$= T_+(k_1, k_2, k_3) + i T_-(k_1, k_2, k_3). \qquad (5)$$

The relationship between the right-/left-trispectra and the parity-even/parity-odd trispectra are,

$$T_{R/L} = T_+ \pm i\pi_{123}T_-, \qquad (6)$$

$$T_{\pm} = \frac{1}{2}(-i\pi_{123})^{(1\mp 1)/2}(T_R \pm T_L). \qquad (7)$$

Due to isotropy, the only vectors we can use to form a parity-odd structure are $k_1, k_2$, and $k_3$, so $T_-(k_1, k_2, k_3)$ must be proportional to the triple product $k_1 \cdot (k_2 \times k_3)$. We parametrize the shape of the imaginary trispectrum as (Coulton, Philcox & Villaescusa-Navarro 2024)

$$T_-(k_1, k_2, k_3) = k_1 \cdot (k_2 \times k_3) \tau_-(k_1, k_2, k_3, k_4, K, \tilde{K}). \qquad (8)$$

The left-hand side of equation (1) is totally symmetric under the interchange of any of the wave vectors, and the triple product in equation (8) is totally antisymmetric. These symmetries require that $\tau_-(k_1, k_2, k_3, k_4, K, \tilde{K})$ is totally antisymmetric under the interchange of any two of the four wave vectors. The trispectrum cannot depend on one diagonal and not the other since $K \to \tilde{K}$ under the interchange of $k_2$ and $k_4$. Also, under the interchange of $k_2$ and $k_3$, $K$ becomes $|k_1 + k_3|$, which is not independent, since

$$|k_1 + k_3|^2 = -K^2 - \tilde{K}^2 + k_1^2 + k_2^2 + k_3^2 + k_4^2. \qquad (9)$$

This parametrization is not unique. We could choose any six geometrically independent quantities of the tetrahedra to parametrize the continuous degrees of freedom. For example, we could reparametrize by replacing $K$ with the angle between $k_1$ and $k_2$. Similarly, we can choose either parity or helicity to parametrize the discrete degree of freedom.

## 3 POP SPECTRA

CMB and LSS observables are consistent with cosmic structure originating from a single primordial scalar potential, $\Phi(x)$. In this section, we construct POP-spectrum-like observables, which we call POP spectra. These are compressions of the 6D parity-odd trispectrum down to a 1D power spectrum of composite fields. First, we construct the vector POP spectrum, defined below in equation (14), by cross-correlating a vector field with a pseudovector field, both derived as composite fields quadratic in $\Phi$. Next, we construct the scalar POP spectrum, defined below in equation (24), by cross-correlating a scalar field with a pseudoscalar field. The latter is formed as a composite field that is cubic in $\Phi$.

### 3.1 Vector POP spectrum

Our goal in this subsection is to construct composite vector and pseudovector fields out of a scalar field such that the cross-power spectrum between them is a compressed estimator of the parity-odd trispectrum. Since a vector is parity odd and a pseudovector is parity even, their cross-correlation will be a POP spectrum. Obtaining a vector field from a scalar is trivial: take its gradient, $\nabla\Phi(x)$. Obtaining a pseudovector is non-trivial. One possibility is to convolve the field with a filtering or smoothing function $f_a(x)$:

$$\Phi_a(x) = \int d^3y f_a(x - y)\Phi(y) \qquad (10)$$

$$= \int_k f_a(k) \Phi(k) e^{ik \cdot x}, \qquad (11)$$

where the integral notation is defined in Table 1. For two distinct filtering functions, $f_a(x)$ and $f_b(x)$, the composite vector field $V_{ab}(x) \equiv$
<table><caption>**Table 1.** Here, we list our Fourier-space integral notation and conventions.</caption><thead><tr><th>Notation</th><th>Definition</th></tr></thead><tbody><tr><td>∫<sub>k</sub></td><td>∫ <sup>d<sup>3</sup>k</sup> / (2π)<sup>3</sup></td></tr><tr><td>∫<sub>q<sub>1</sub>...q<sub>n</sub></sub></td><td>∫<sub>q<sub>1</sub></sub> ... ∫<sub>q<sub>n</sub></sub></td></tr><tr><td>∫<sub>q</sub></td><td>∫<sub>0</sub><sup>∞</sup> <sup>q<sup>2</sup> dq</sup> / (2π<sup>2</sup>)</td></tr><tr><td>∫<sub>q<sub>1</sub>...q<sub>n</sub></sub></td><td>∫<sub>q<sub>1</sub></sub> ... ∫<sub>q<sub>n</sub></sub></td></tr></tbody></table>

The bold symbols indicate wave vectors and the integral is over all of Fourier space. The non-bold symbols indicate the wave vector magnitude.

$\Phi_a(\mathbf{x}) \nabla \Phi_b(\mathbf{x})$ is not the gradient of a scalar, so it also has a curl component. Therefore, we can carry out a Helmholtz decomposition,

$$V_{ab}(\mathbf{x}) = \nabla \phi_{ab}(\mathbf{x}) - \nabla \times A_{ab}(\mathbf{x}). \qquad (12)$$

Since $V_{ab}$ and the spatial derivatives are parity odd, $A_{ab}$ is a pseudovector. Notice that if the two fields, $\Phi_a$ and $\Phi_b$, were identical, $A_{ab}$ would vanish. Although one of these can be the original, unfiltered field. We impose that its divergence vanishes because any possible divergence would not contribute to $V_{ab}(\mathbf{x})$. Then we solve for its modes by taking the curl of $V_{ab}(\mathbf{x})$ and solving the resulting Poisson equation,

$$A_{ab}(\mathbf{k}) = -\frac{1}{k^2} i \mathbf{k} \times V_{ab}(\mathbf{k}). \qquad (13)$$

In general, we are free to choose four distinct filter functions: $f_a(\mathbf{x})$, $f_b(\mathbf{x})$, $f_c(\mathbf{x})$, and $f_d(\mathbf{x})$. From the first two, we define the composite vector field $V_{ab}(\mathbf{x})$. From the second two, we define the composite vector field $V_{cd}(\mathbf{x})$ and then isolate the pseudovector field $A_{cd}(\mathbf{x})$ from it. The dot product $V_{ab}(\mathbf{x}) \cdot A_{cd}(\mathbf{x})$ is parity odd. As a result, we can define our first POP spectrum as

$$(V_{ab}(\mathbf{k}) \cdot A_{cd}(\mathbf{k}')) = (2\pi)^3 \delta_D^{(3)}(\mathbf{k} + \mathbf{k}') P_{\text{vector}}(\mathbf{k}). \qquad (14)$$

Explicitly, we can express the power spectrum as

$$P_{\text{vector}}(k) = \frac{1}{\mathcal{N}_{\text{vector}}(k)} \int_{q_1, q_2, q_3, q_4} \left[ (2\pi)^6 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 - \mathbf{k}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 + \mathbf{k}) f_a(\mathbf{q}_1) f_b(\mathbf{q}_2) f_c(\mathbf{q}_3) f_d(\mathbf{q}_4) \times [\mathbf{k} \cdot (\mathbf{q}_1 \times \mathbf{q}_3)]^2 \tau_-(q_1, q_2, q_3, q_4, k, |\mathbf{q}_1 + \mathbf{q}_4|) \right]. \qquad (15)$$

The wave vector integrals are defined in Table 1. The normalization factor $\mathcal{N}_{\text{vector}}(k)$ is somewhat arbitrary, but if chosen poorly, the estimator will depend strongly on the geometry and resolution of the region where we sample $\Phi(\mathbf{x})$. A sensible choice is

$$\mathcal{N}_{\text{vector}}(k) = \int_{q_1, q_2, q_3, q_4} (2\pi)^6 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 - \mathbf{k}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 + \mathbf{k}) f_a(\mathbf{q}_1) f_b(\mathbf{q}_2) f_c(\mathbf{q}_3) f_d(\mathbf{q}_4) [\mathbf{k} \cdot (\mathbf{q}_1 \times \mathbf{q}_3)]^2, \qquad (16)$$

so that equation (15) is a weighted average over the parity-odd trispectrum shape $\tau_-(q_1, q_2, q_3, q_4, k, |\mathbf{q}_1 + \mathbf{q}_4|)$.

If we choose isotropic filtering functions that depend only on the magnitudes of the $\mathbf{q}_i$, the vector POP spectrum is a real-valued quantity that is non-vanishing only if the parity-odd trispectrum is non-vanishing. With isotropic filtering functions, the parity-even part of the trispectrum does not contribute to the imaginary part of the POP spectrum. To see this, we note that the parity-even trispectrum is invariant under $\mathbf{q}_i \rightarrow -\mathbf{q}_i$ while the triple product receives a minus sign. Therefore, the integrand is antisymmetric and the integral vanishes.<fnref n="1" />

The vector-pseudovector construction compresses the trispectrum down to a power spectrum by correlating two fields of order $\Phi^2$. The modes of these fields represent triangles of wave vectors satisfying $\mathbf{k} + \mathbf{q}_1 + \mathbf{q}_2 = 0$ and $\mathbf{k}' + \mathbf{q}_3 + \mathbf{q}_4 = 0$. Computing the power spectrum enforces that the two triangles connect along a side of equal length, $k$, which imposed $\mathbf{k} + \mathbf{k}' = 0$. Thus, the two triangles form a tetrahedron (shown in Fig. 1) as the trispectrum requires. The vector POP spectrum is a weighted average over all trispectrum configurations while holding fixed the diagonal side length $k$, which is the argument of this POP spectrum. For this reason, the vector POP carries crucial information about the soft limit where the diagonal approaches zero.

If the trispectrum has no explicit dependence on the diagonals and the filtering functions depend only on the magnitudes of the wave vectors in Fourier space, we can compute the angular integrals in equations (15) and (16) analytically. Then the expression for the vector POP

<footnote n="1">¹The sign of the $\mathbf{k}$ in the Dirac delta functions does not matter, since $P_{\text{vector}}(\mathbf{k})$ is a function of only the magnitude of $\mathbf{k}$.</footnote>
spectrum simplifies (as detailed in Appendix A) to

$$P_{\text{vector}}(k) = \frac{\pi^4}{2\mathcal{N}_{\text{vector}}(k)} \int_{q_1, q_2, q_3, q_4} \left[ \frac{q_1}{q_2} \Theta(\sin^2\theta_{12}) \sin^2\theta_{12} \frac{q_3}{q_4} \Theta(\sin^2\theta_{34}) \sin^2\theta_{34} f_a(q_1) f_b(q_2) f_c(q_3) f_d(q_4) \times \tau_-(q_1, q_2, q_3, q_4) \right], \quad (17)$$

and

$$\mathcal{N}_{\text{vector}}(k) = \frac{\pi^4}{2} \int_{q_1, q_2, q_3, q_4} \frac{q_1}{q_2} \Theta(\sin^2\theta_{12}) \sin^2\theta_{12} \frac{q_3}{q_4} \Theta(\sin^2\theta_{34}) \sin^2\theta_{34} f_a(q_1) f_b(q_2) f_c(q_3) f_d(q_4). \quad (18)$$

The radial Fourier integrals are defined in Table 1. Here, $\Theta(x)$ is the Heaviside function,

$$\Theta(x) = \begin{cases} 0, & \text{if } x < 0 \\ \frac{1}{2}, & \text{if } x = 0, \\ 1, & \text{if } x > 0 \end{cases}, \quad (19)$$

and $\theta_{12}$ the angle between $q_1$ and $k$ in the first triangle, so

$$\sin^2 \theta_{12} = 1 - \left( \frac{q_1^2 + k^2 - q_2^2}{2q_1 k} \right)^2, \quad (20)$$

and similarly for $\theta_{34}$, with $1, 2 \to 3, 4$. Notice that if $q_i, q_j$, and $k$ fail to satisfy the triangle inequalities $|q_j - q_i| \le k \le q_i + q_j$, then the right-hand sides in equation (20) would be negative and the corresponding sine would be purely imaginary. Thus, the Heaviside functions in equations (17) and (18) enforce the triangle inequalities.

As a final remark, note that the pseudovector defined in equation (13) is not unique. For example, one could alternatively take the following curl,

$$\mathbf{B}_{ab}(\mathbf{x}) = \nabla \Phi_a(\mathbf{x}) \times \nabla \Phi_b(\mathbf{x}). \quad (21)$$

Using this instead of $\mathbf{A}_{ab}$ simply multiplies the vector POP spectrum $P_{\text{vector}}(k)$ by a factor of $-k^2$. Instead of applying separate filtering functions to each field to the cross-product does not vanish, we could define an antisymmetric convolution kernel. Using convolution kernels is similar to the skew spectrum (Schmittfull et al. 2015), and is the method used for clustering fossil estimators (Jeong & Kamionkowski 2012).<fnref n="2" /> These are better suited to model-specific optimal estimators, whereas in this work we are taking a more model-independent approach.

## 3.2 Scalar POP spectrum

We can construct another POP spectrum by cross-correlating a scalar field with a pseudoscalar field. The construction of a pseudoscalar, similar to that of the pseudovector from the previous subsection, is non-trivial and requires convolving the scalar field with filtering functions. In this case, we first construct the pseudovector given in equation (21). We have the freedom to further smooth this pseudovector by another filtering function, with modes

$$\mathbf{B}_{bcd}(\mathbf{Q}) = -f_d(\mathbf{Q}) \int_{q_3:q_4} (2\pi)^3 \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathbf{Q}) f_b(\mathbf{q}_3) f_c(\mathbf{q}_4) (\mathbf{q}_3 \times \mathbf{q}_4) \Phi(\mathbf{q}_3) \Phi(\mathbf{q}_4). \quad (22)$$

Then, we construct the triple product field

$$\Psi_{abcd}(\mathbf{x}) = \nabla \Phi_a(\mathbf{x}) \cdot \mathbf{B}_{bcd}(\mathbf{x}), \quad (23)$$

which is a pseudoscalar. Its cross-power spectrum with $\Phi(\mathbf{k})$ is a POP spectrum,

$$\langle \Phi(\mathbf{k}) \Psi_{abcd}(\mathbf{k}') \rangle = (2\pi)^3 \delta_D^{(3)}(\mathbf{k} + \mathbf{k}') P_{\text{scalar}}(\mathbf{k}). \quad (24)$$

The scalar POP spectrum is given by

$$P_{\text{scalar}}(k) = \frac{1}{\mathcal{N}_{\text{scalar}}(k)} \int_{q_2, q_3, q_4, Q} \left[ (2\pi)^6 \delta_D^{(3)}(\mathbf{k} + \mathbf{q}_2 + \mathbf{Q}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathbf{Q}) f_a(\mathbf{q}_2) f_b(\mathbf{q}_3) f_c(\mathbf{q}_4) f_d(\mathbf{Q}) \times [\mathbf{Q} \cdot (\mathbf{q}_2 \times \mathbf{q}_4)]^2 \tau_-(\mathbf{k}, \mathbf{q}_2, \mathbf{q}_3, \mathbf{q}_4, \mathbf{Q}, |\mathbf{q}_2 + \mathbf{q}_3|) \right], \quad (25)$$

<footnote n="2">²There are conditions under which the vector-mode clustering fossil estimator and our vector POP spectrum exactly coincide (up to the choice of normalization and inverse-variance-weighting). Referring to equation (1) of Jeong & Kamionkowski (2012), if their mode coupling kernels for right-/left-handed vector modes are the identical and have the form $f_{R/L}(\mathbf{k}_1, \mathbf{k}_2) = f_a(k_1) f_b(k_2) - f_b(k_1) f_a(k_2)$, then this is the same construction as our vector POP spectrum in the case where $f_c(k) = f_a(k)$ and $f_d(k) = f_b(k)$. There is not yet an equivalent or analogous estimator in the clustering fossil approach for the scalar POP spectrum that we present in the next section.</footnote>
with normalization

$$N_{\text{scalar}}(k) = \int_{q_2, q_3, q_4, Q} (2\pi)^6 \delta_D^{(3)}(\mathbf{k} + \mathbf{q}_2 + \mathbf{Q}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathbf{Q}) f_a(\mathbf{q}_2) f_b(\mathbf{q}_3) f_c(\mathbf{q}_4) f_d(\mathbf{Q}) [\mathbf{Q} \cdot (\mathbf{q}_2 \times \mathbf{q}_4)]^2. \quad (26)$$

As with the vector-pseudovector construction, we choose the normalization so that the POP spectrum is a weighted average over the imaginary trispectrum. In this case, we have constructed the compressed four-point statistic by cross-correlating a pseudoscalar of order $\Phi^3$ with the original potential $\Phi$. The resulting power spectrum is an average over the trispectrum tetrahedra, fixing the magnitude of one of the wave vector sides. The average integrates over the diagonal $Q$ this time (corresponding to the $k$ in equations (14)–(18), and the scalar POP spectrum is a function of only the fixed wave vector magnitude $k$, a side length of the tetrahedron. This property makes the scalar POP particularly useful for probing the soft limit where one of the side lengths approaches zero.

The side length that we label $k$ here was labelled $q_1$ for the vector POP spectrum in equations (15)–(18), so in the following equations a subscript 1 will refer to this side length. Again, if the trispectrum does not depend on the diagonals and the filtering functions depend only on the magnitudes of the wave vectors, we can compute the angular integrals analytically, and the expression simplifies (as detailed in Appendix A) to

$$P_{\text{scalar}}(k) = \frac{\pi^4}{2k N_{\text{scalar}}(k)} \int_{q_2, q_3, q_4, Q} \left[ q_2 \Theta(\sin^2\theta_{21}) \sin^2\theta_{21} \frac{q_4}{q_3} \Theta(\sin^2\theta_{43}) \sin^2\theta_{43} f_a(q_2) f_b(q_3) f_c(q_4) f_d(Q) \right. \\ \left. \times \tau_-(k, q_2, q_3, q_4) \right], \quad (27)$$

and

$$N_{\text{scalar}}(k) = \frac{\pi^4}{2k} \int_{q_2, q_3, q_4, Q} q_2 \Theta(\sin^2\theta_{21}) \sin^2\theta_{21} \frac{q_4}{q_3} \Theta(\sin^2\theta_{43}) \sin^2\theta_{43} f_a(q_2) f_b(q_3) f_c(q_4) f_d(Q). \quad (28)$$

Here, we evaluate $\sin^2 \theta_{21}$ and $\sin^2 \theta_{43}$ using

$$\sin^2 \theta_{21} = 1 - \left( \frac{q_2^2 + Q^2 - k^2}{2q_2 Q} \right)^2, \quad (29)$$

$$\sin^2 \theta_{43} = 1 - \left( \frac{q_4^2 + Q^2 - q_3^2}{2q_4 Q} \right)^2. \quad (30)$$

### 3.3 Constructing the estimators

To construct the vector POP spectrum from a scalar field $\Phi(\mathbf{x})$ on a discrete grid, we start by taking its Fourier transform, $\Phi(\mathbf{q})$. Then, we make four copies of these modes and rescale them by the filter functions {$f_a(\mathbf{q}), f_b(\mathbf{q}), f_c(\mathbf{q}), f_d(\mathbf{q})$}, obtaining {$\Phi_a(\mathbf{q}), \Phi_b(\mathbf{q}), \Phi_c(\mathbf{q}), \Phi_d(\mathbf{q})$}. Multiplying the modes of $\Phi_b(\mathbf{q})$ and $\Phi_d(\mathbf{q})$ by $i\mathbf{q}$ and then inverse Fourier transforming yields the modes of $\nabla\Phi_b(\mathbf{x})$ and $\nabla\Phi_d(\mathbf{x})$. We also inverse Fourier transform the other two fields, $\Phi_a$ and $\Phi_c$ and, now in position space, form the products $V_{ab}(\mathbf{x}) = \Phi_a(\mathbf{x})\nabla\Phi_b(\mathbf{x})$ and $V_{cd}(\mathbf{x}) = \Phi_c(\mathbf{x})\nabla\Phi_d(\mathbf{x})$. Next, we Fourier transform both of these and, from $V_{cd}(\mathbf{q})$, solve for the modes of $A_{cd}(\mathbf{q})$ using equation (13). The unbinned vector POP spectrum estimator is

$$P_{\text{vector}}(k) = \frac{1}{V_{\text{box}} N_{\text{vector}}(k)} \sum_{|\mathbf{q}|=k} V_{ab}(\mathbf{q}) \cdot A_{cd}^*(\mathbf{q}), \quad (31)$$

where the sum is over all modes with equal wave vector magnitudes on the grid.

The normalization is computed by initializing four grids of modes corresponding to the filter functions: $f_a(\mathbf{q}), f_b(\mathbf{q}), f_c(\mathbf{q}),$ and $f_d(\mathbf{q})$. From here, obtaining the normalization corresponding to equation (16) is somewhat complicated. When constructing the power spectrum, we get one triple product from the estimator and another from the parity-odd trispectrum. Overall, the power spectrum involves the square of the triple product. To get a squared triple product in the normalization factor, we need to compute the Hessian of the field $f_a(\mathbf{x})$ and construct the composite tensor field with components

$$F_{ab}^{ij}(\mathbf{x}) = f_a(\mathbf{x}) \nabla^i \nabla^j f_b(\mathbf{x}). \quad (32)$$

After Fourier transforming the Hessian of the filter function fields, we take the curl of one of the two components,

$$J_{ab}^{ij}(\mathbf{q}) = \epsilon^{imn} i k^m F_{ab}^{nj}(\mathbf{q}), \quad (33)$$

where repeated indices imply summation. Following the same steps for $f_c(\mathbf{x})$ and $f_d(\mathbf{x})$, the desired normalization is

$$N_{\text{vector}}(k) = - \sum_{|\mathbf{q}|=k} J_{ab}^{ij}(\mathbf{q}) [J_{cd}^{ij}(\mathbf{q})]^{\dagger}. \quad (34)$$
The dagger † denotes the Hermitian conjugate of the matrix. We then bin the estimator in wave number bins of width Δk,

$$P_{\text{vector}}(k|\Delta k) = \frac{1}{N_{\text{vector}}(k|\Delta k)} \sum_{q \in k} P_{\text{vector}}(q) N_{\text{vector}}(q), \qquad (35)$$

where

$$N_{\text{vector}}(k|\Delta k) = \sum_{q \in k} N_{\text{vector}}(q). \qquad (36)$$

Here, the sum over $q \in k$ means summing over all $q \in [k - \Delta k/2, k + \Delta k/2)$.

The computation of the scalar POP spectrum follows similar steps. In this case, we compute the gradients $\nabla\Phi_a(\mathbf{x})$, $\nabla\Phi_b(\mathbf{x})$ and $\nabla\Phi_c(\mathbf{x})$. We use the latter two to form $B_{bc}(\mathbf{x})$ according to equation (21). We Fourier transform this and then multiply by the filtering function modes $f_d(\mathbf{q})$, resulting in the modes of $B_{bcd}(\mathbf{q})$. We then inverse Fourier transform this, and its dot product with $\nabla\Phi_a(\mathbf{x})$ yields $\Psi_{abcd}(\mathbf{x})$. The unbinned scalar POP spectrum is

$$P_{\text{scalar}}(k) = \frac{1}{V_{\text{box}} N_{\text{scalar}}(k)} \sum_{|\mathbf{q}|=k} \Phi(\mathbf{q}) [\Psi_{abcd}(\mathbf{q})]^* . \qquad (37)$$

Constructing the normalization is again non-trivial due to the squared triple product factor. In this case, we require the Hessians of the field $f_a(\mathbf{x})$, $f_b(\mathbf{x})$, and $f_c(\mathbf{x})$. From the latter two, we construct the tensor field $F_{bc}^{ij}(\mathbf{x})$ as in equation (32). By Fourier transforming all components and multiplying the modes by the filter function $f_d(\mathbf{q})$, we obtain $F_{bcd}^{ij}(\mathbf{q}) \equiv f_d(\mathbf{q})F_{bc}^{ij}(\mathbf{q})$. Inverse Fourier transforming this field and multiplying by the Hessian of $f_a(\mathbf{x})$ gives

$$G_{abcd}^{ij}(\mathbf{x}) = \epsilon^{imn} \epsilon^{jpq} F_{bcd}^{mp}(\mathbf{x}) \nabla^n \nabla^q f_a(\mathbf{x}). \qquad (38)$$

Fourier transforming this, the scalar POP spectrum normalization is

$$N_{\text{scalar}}(k) = - \sum_{|\mathbf{q}|=k} k^i k^j G_{abcd}^{ij}(\mathbf{q}). \qquad (39)$$

The binned estimator is

$$P_{\text{scalar}}(k|\Delta k) = \frac{1}{N_{\text{scalar}}(k|\Delta k)} \sum_{q \in k} P_{\text{scalar}}(q) N_{\text{scalar}}(q), \qquad (40)$$

where

$$N_{\text{scalar}}(k|\Delta k) = \sum_{q \in k} N_{\text{scalar}}(q). \qquad (41)$$

If we choose filter functions that only select wave vector magnitude shells of width Δk, then both POP spectra coincide and equal the binned trispectrum estimator in equation (A30). This fact illustrates the computational advantage of our POP spectra. Evaluating all the bins of one POP spectrum is computationally equivalent to constructing a 1D slice of the full trispectrum. The tradeoff is that we have significantly compressed the information of the full trispectrum. However, the presence of two estimators, each carrying complementary information from the parity-odd trispectrum, partially compensates for this reduction in information. We also have the freedom to choose the filter functions, which can be optimized. The extent of information loss could be further mitigated by considering extensions and generalizations of the POP spectra, which we explore in Section 5.

The algorithms that we presented above for constructing these estimators are efficient since the most costly operations are the Discrete Fourier Transforms (DFTs), which scale as $N \log N$, with N the number of grid points. All other operations are point-wise multiplications, which scale as N. Ignoring the normalization factors, the vector POP spectrum requires 15 and the scalar POP spectrum requires 17 3D DFTs. Constructing the estimators without the normalization factor takes about 6 s for the vector and 8 s for the scalar POP spectrum running on a single node with 128 cores. The normalization factors are more expensive to compute, but these only need to be computed once and can be saved and reused since they are not data-dependent.

# 4 VALIDATION OF THE ESTIMATORS

## 4.1 Random parity-violating realizations

To validate our POP estimators, we generate a set of random primordial curvature fields that include a specific parity-violating trispectrum. We start by generating the modes of a Gaussian random field primordial curvature perturbation $\zeta_G(\mathbf{k})$,

$$\zeta_G(\mathbf{k}) = \sqrt{\frac{V_{\text{Box}} \pi^2 A_s}{k^3}} \left(\frac{k}{k_p}\right)^{n_s-1} [N_1(\mathbf{k}) + i N_2(\mathbf{k})], \qquad (42)$$

where $N_1$ and $N_2$ are drawn from a standard normal distribution, $A_s$ is the amplitude of the primordial power spectrum, $n_s$ is the spectral tilt, and $k_p = 0.05 \text{ Mpc}^{-1}$ is the pivot scale. We then transform this Gaussian random field into a non-Gaussian field in real space as (Coulton et al.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Here we show the trispectrum in the thin bin limit (top panel) and the binning scheme (lower six panels). The shape of the trispectrum template from equation (45) appears in the top panel. The bottom six panels indicate the binning scheme with $n_i = k_i/k_F$, $N = K/k_F$, and $\tilde{N} = \tilde{K}/k_F$. Here, $k_F$ is the fundamental mode of a box interpreted as having length $4 \text{ Gpc h}^{-1}$. We impose $k_1 < k_2 < k_3 < k_4$. We have used 10 k-bins of width $k_F$.</figcaption></figure>

2024)

$$ \zeta(\mathbf{x}) = \zeta_G(\mathbf{x}) + g_- \nabla \zeta_G^{[\alpha]}(\mathbf{x}) \cdot [\nabla \zeta_G^{[\beta]}(\mathbf{x}) \times \nabla \zeta_G^{[\gamma]}(\mathbf{x})]. \quad (43) $$

The coefficient $g_-$ controls the amplitude of the primordial parity-odd trispectrum, and the modes of the fields in the triple product are

$$ \zeta_G^{[\alpha]}(\mathbf{k}) = k^\alpha \zeta_G(\mathbf{k}). \quad (44) $$

The leading order imaginary part of the trispectrum for this template is given by

$$ T_-(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = g_- \mathbf{k}_1 \cdot (\mathbf{k}_2 \times \mathbf{k}_3) (2\pi^2 A_s)^3 (k_1^{\alpha-4+n_s} k_2^{\beta-4+n_s} k_3^{\gamma-4+n_s} k_4^0 \mp 23 \text{ signed permutations}). \quad (45) $$

The 24 terms in the parenthesis are the permutations of {$k_1, k_2, k_3, k_4$}. Even permutations (with an even number of transpositions) get a positive sign, and odd permutations get a negative sign in the sum. We have included the $k_4^0$ factor to emphasize that every term has one of the four wave vector magnitudes raised to the exponent zero.

The trispectrum is scale invariant if $\alpha + \beta + \gamma = -3$. To preserve the large-scale power spectrum we also restrict these exponents to be less than or equal to zero. We choose $\alpha = -2$, $\beta = -1$, and $\gamma = 0$, which is the same template used in Coulton et al. (2024). For simplicity, we choose $A_s = 2 \times 10^{-9}$, $n_s = 1$, and $g_- = \pm 10^6$. In Fig. 2, we display the shape of this trispectrum.

Choosing a binning scheme where $k_1 < k_2 < k_3 < k_4$, the trispectrum shape function is dominated by the terms

$$ \tau_-(k_1, k_2, k_3, k_4) \approx \frac{g_-(2\pi^2 A_s)^3}{k_1^5 k_2^4} \left( \frac{1}{k_3^3} - \frac{1}{k_4^3} \right). \quad (46) $$

This shape peaks when the wave vectors $\mathbf{k}_1, \mathbf{k}_2$, and $\mathbf{k}_3$ are parallel so that $k_4 = k_1 + k_2 + k_3$. The diagonals for this configuration are $K = k_1 + k_2$ and $\tilde{K} = k_4 - k_1$. However, this configuration is insensitive to parity, since the triple product vanishes when any two wave vectors are colinear. The trispectrum peaks on a tetrahedron shape that deviates from colinearity with diagonals $K < k_1 + k_2$ and $\tilde{K} > k_4 - k_1$ that maximize the product of the right-hand side of equation (46) and the triple product $\mathbf{k}_1 \cdot (\mathbf{k}_2 \times \mathbf{k}_3)$. This peak analysis is true for any scale-invariant trispectrum of the form in equation (45), although the particular diagonals that maximize the trispectrum depend on the values of the template exponents. These parity-odd trispectrum templates diverge in the soft limit $k_1 \to 0$.

Since we are dealing with a scale-invariant primordial potential, the box length is arbitrary, but we interpret it as $L_{\text{Box}} = 4 \text{ Gpc h}^{-1}$ so that the scales involved are relevant for CMB and LSS. We generate 64 random pairs of these non-Gaussian potential fields on a grid of size $N_{\text{grid}} = 512^3$. The pairs have the same underlying Gaussian realization but opposite signs for their non-linear terms. Subtracting the POP spectra measured on these pairs suppresses cosmic variance from the purely Gaussian part of the fields. As can be seen in Fig. 3, the cosmic
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The vector POP spectrum of the primordial potential defined in equation (43). In the left panels, the amplitude of the primordial trispectrum is $g_{-} = 10^{6}$, while in the middle panels $g_{-} = -10^{6}$. We used the same underlying Gaussian realizations for both signs of $g_{-}$. The right panels show the variance-suppressed estimator, taking the difference between results with positive and negative $g_{-}$. The dashed curves are the expected signal, computed from equations (17) and (18). The bottom row shows fractional residuals of the simulated data with respect to the semi-analytical calculation. Square data points show the estimator on a Gaussian random field. These scatter around zero and are consistent with a vanishing parity-violating signal. All data points have error bars obtained by bootstrap averaging over the 64 simulations, resampling with replacement $10^{5}$ times. In the right-most panel, the sample variance cancellation is not perfect. This is because the POP spectrum involves four weakly non-Gaussian fields.</figcaption></figure>

variance is not perfectly cancelled. This is because the POP spectrum is a compression of the four-point correlator of a weakly non-Gaussian field.

There is a considerable amount of freedom in our estimators since we can choose any set of filtering functions. For the vector POP spectrum, $f_a(k)$ must be distinct from $f_b(k)$, and $f_c(k)$ must be distinct from $f_d(k)$, otherwise $P_{\text{vector}}(k)$ vanishes identically. For the scalar POP spectrum, $f_a(k)$, $f_b(k)$, and $f_c(k)$ must all be distinct for $P_{\text{scalar}}(k)$ not to vanish. Here, we choose

$$f_a(k) = k^2 \Theta(k - k_{\text{min}}) \Theta(k_{\text{max}} - k), \qquad (47)$$

$$f_c(k) = k^{-2} \Theta(k - k_{\text{min}}) \Theta(k_{\text{max}} - k), \qquad (48)$$

$$f_b(k) = f_d(k) = \Theta(k - k_{\text{min}}) \Theta(k_{\text{max}} - k). \qquad (49)$$

Each filtering function involves the same scale cuts, $k_{\text{min}} = 5 \times 10^{-3} \text{ Mpc}^{-1} \text{ h}$ and $k_{\text{max}} = 2 \times 10^{-1} \text{ Mpc}^{-1} \text{ h}$, which selects a thick spherical shell of modes. As long as the fundamental mode in the box is less than $k_{\text{min}}$, and the Nyquist mode is greater than $k_{\text{max}}$, the box contains the full range of modes and the estimator should be independent of the box size and resolution. However, the discreteness of the grid is noticeable on large scales if $k_{\text{min}}$ is very close to the fundamental mode.

The autocorrelation of the non-linear term in equation (43) will introduce corrections to the power spectrum. These corrections are of the order $\sim g_{-}^2 A_s^3$ and are subdominant compared to the Gaussian power spectrum for our template at the scales we consider. The non-Gaussian corrections to the power spectrum grow at small scales and increase the power of the Nyquist mode, $k = 0.4 \text{ Mpc}^{-1} \text{ h}$, in our simulated boxes by about 1 per cent. After imposing the scale cuts in our POP spectra the corrections to the power spectrum are far below percent level for the modes we analyse. These corrections to the power spectrum can be removed by further rescaling the non-Gaussian field (Coulton et al. 2024). However, corrections to the power spectrum could be a physical effect. The (inflationary) mechanisms that generate the primordial trispectrum can also affect the shape of the power spectrum. This template does not introduce a primordial bispectrum because all contributions would involve an odd number of Gaussian fields, resulting in a vanishing expectation value.

We have omitted many complications of observational survey data in this analysis. Cosmological surveys target observables that are biased tracers of the underlying matter field, which has undergone non-linear evolution. Accurate modelling will require the linear matter transfer function, the bias expansion, and possibly the non-linearity of gravitational clustering depending on the scales considered. Surveys
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Fig. 3 but for the scalar POP spectrum. In this case, the signal grows in the limit $k \to 0$, showing sensitivity to the soft limit of the trispectrum as one of the tetrahedron wave vector side lengths becomes small. In contrast, Fig. 3 demonstrates decreasing sensitivity in the diagonal soft limit of the vector POP spectrum because the specific trispectrum template that we simulate does not peak in this limit. The behaviour of the residual in the right-most panel is again (as in Fig. 3) due to the additional sample variance contributed by the weakly non-Gaussian nature of the fields.</figcaption></figure>

observe objects in redshift space where isotropy does not hold. Observational analyses must also model the survey mask, selection function, and shot noise. We will investigate the impact of these aspects of observational data on the POP spectra in future work.

## 4.2 COMPARISON WITH ANALYTICAL COMPUTATION

In Fig. 3, we display the vector POP spectrum compared with the theoretical expectation obtained by numerically evaluating the integrals in equations (17) and (18). The left panel shows the results for 64 simulations with $g_{-} = +10^{6}$. The middle panel shows results for the same underlying Gaussian random fields but with $g_{-} = -10^{6}$. We significantly reduce cosmic variance by computing the difference between these two estimators, as illustrated in the right panel of Fig. 3. We estimated the error bars by bootstrap averaging over the estimators constructed for the 64 realizations, resampling $10^{5}$ times with replacement.

We obtained the semi-analytic result in Fig. 3 by evaluating equation (17) on a discretized lattice for the radial integrals over all $q_i$ and $k = |q_1 + q_2|$. Specifically, for a given bin, we construct a 3D grid of values for {$q_1, q_2, k$}, with $q_1$ and $q_2$ spaced uniformly between $k_{min}$ and $k_{max}$ and $k$ spanning the width of the power spectrum bin. We then mask out regions that violate the triangle inequalities. Next, we construct a second grid for {$q_3, q_4, k$}, masking the non-triangular regions. Finally, we integrate equation (17) term by term. For each permutation term in the trispectrum of equation (45), we construct the integrand for $q_1$ and $q_2$, and compute the discretized integral using Romberg integration with $2^8 + 1$ sample points. We do the same for $q_3$ and $q_4$, except we include explicit factors of $k$ appearing in equation (45) in this grid. We then multiply the results from the two grids. Finally, we compute the discretized integral of the diagonal $k$ using Romberg integration with $2^7 + 1$ sample points. The procedure for integrating equation (18) is the same, omitting the factors from the primordial trispectrum template shape. We parallelize this calculation by computing multiple bins simultaneously. From Fig. 3, we see excellent agreement with the semi-analytic calculation of the expected signal. This agreement demonstrates that we accurately recover the expected signal from the trispectrum shape injected into the data.

In Fig. 4, we display similar results for the scalar POP spectrum. The theoretical calculation of equations (27) and (28) differs from the vector case in that $k$ in the triangle {$k, q_2, Q$} is integrated over the estimator's bin width spanning the range between $k_{min}$ and $k_{max}$ and the filter function $f_d(Q)$ applies to the diagonal $Q$. Again, we find excellent agreement between our simulated data and the expected shape of $P_{\text{scalar}}(k)$. The bias at low $k$ is due to the discreteness of the grid, compared with the continuous-limit integrals computed in the semi-analytic prediction. This bias is not noticeable for the vector estimator in Fig. 3 due to the signal vanishing at low $k$.

As a final remark note that, unlike the vector POP spectrum, the scalar POP spectrum increases as $k$ approaches zero. This occurs because the argument $k$ of the scalar POP spectrum corresponds to the side length of the tetrahedron. Our specific trispectrum template, equation (45), diverges as one side length goes to zero, so the scalar POP spectrum also diverges in this limit. If we chose a template that peaks in the soft limit where a diagonal approaches zero, the vector POP spectrum would increase as its argument approaches zero. We could also construct a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** On the left panel, we display the correlation matrix, equation (50), for the vector-pseudovector (upper triangle) and scalar-pseudoscalar (lower triangle) estimators. The correlation matrix is dominated by the diagonal, with off-diagonal contributions less than 10 per cent. On the right panels, we display histograms of the $\chi^2$ values from 64 random non-Gaussian fields compared with the distribution from Gaussian realizations. The peak of the dashed analytical $\chi^2$ curve coincides with the number of bins analysed.</figcaption></figure>

template that diverges in the soft side length and soft diagonal limits. Then the vector and scalar POP spectra would increase as their arguments go to zero. This demonstrates that the $k \to 0$ limits of our POP spectra encode crucial and complementary information about the soft limits of the trispectrum.

## 4.3 Correlation matrix and $\chi^2$ test

The POP spectra integrate over many configurations to obtain a power spectrum from a trispectrum. For this reason, we may worry that the covariance matrices of the POP spectra have complicated structures with large off-diagonal contributions. To test this, we generate 1024 realizations of random Gaussian fields, construct the vector and scalar POP spectra estimators, and compute the covariance matrix for each,

$$C(k_i, k_j) = \langle P(k_i) P(k_j) \rangle - \langle P(k_i) \rangle \langle P(k_j) \rangle, \qquad (50)$$

where $k_i$ indicates the $i^{th}$ wave number bin, and $P$ can be either $P_{\text{vector}}(k)$ or $P_{\text{scalar}}(k)$. In Fig. 5 we display the correlation matrices,

$$R(k_i, k_j) = \frac{C(k_i, k_j)}{\sqrt{C(k_i, k_i)C(k_j, k_j)}}, \qquad (51)$$

for the vector POP spectrum in the upper-right triangle and the scalar POP spectrum in the lower-left triangle. We find the correlation matrix is dominated by its diagonal for both estimators, with off-diagonal contributions less than 10 per cent.

We assess the sensitivity of our compressed estimators by computing the $\chi^2$ from the inverse of the covariance matrix,

$$\chi^2 = \sum_{ij} P(k_i) C^{-1}(k_i, k_j) P(k_j). \qquad (52)$$

We compute the $\chi^2$ values for all 64 of our parity-violating non-Gaussian fields with positive amplitude for the primordial trispectrum and their underlying Gaussian fields. Fig. 5 displays histograms of the $\chi^2$ values. The scalar POP spectrum exhibits greater sensitivity for detecting primordial parity violation than the vector POP spectrum. We can understand this from the form of primordial non-Gaussianity in equation (43) that leads to the shape of the trispectrum in equation (45). This shape peaks as $k_1$, the argument of the scalar POP spectrum, goes to zero. We see this in the low-$k$ behaviour in Fig. 4.

The sensitivity depends on the choice of filtering kernels and the shape of the primordial signal. Different parity-odd trispectrum templates will generally yield different $\chi^2$ distributions with more or less sensitivity. In principle, we can find the optimal filter functions to maximize the signal-to-noise ratio of the POP power spectra for a given shape of the underlying primordial trispectrum, as was done for the CMB lensing estimator (Hu 2001), but we do not explore that in this work. Since the vector and scalar POP spectra contain complementary information, it would be best to do a joint analysis of both simultaneously. This joint analysis requires the joint covariance, which will have off-diagonal contributions between the vector and scalar POP spectra bins. We leave the analysis of the joint covariance matrix to future work.
# 5 EXTENSIONS AND GENERALIZATIONS

The examples of POP spectra explored above are the simplest examples of a more general class of constructions. In this section, we outline some of these generalizations.

## 5.1 Connection to the polarization basis

Before exploring some generalizations of our POP spectra, we first write down an alternative and equivalent definition of the vector POP spectrum based on polarization vectors, which will be useful later. As described in Section 2, parity violation is closely related to right- and left-handed helicities. Here we will demonstrate the relationship between our vector POP spectrum and these helicities. To do so, let us take two vector fields $V_{ab}$ and $V_{cd}$. From these, we can cross-correlate their Cartesian components as

$$ \langle V_{ab}^i(\mathbf{k}) V_{cd}^j(\mathbf{k}') \rangle = (2\pi)^3 \delta_D^{(3)}(\mathbf{k} + \mathbf{k}') P^{ij}(\mathbf{k}). \quad (53) $$

Assuming isotropy, we can decompose this power spectrum matrix into its trace, its traceless symmetric, and its antisymmetric parts,

$$ P^{ij}(\mathbf{k}) = \frac{1}{3} \delta^{ij} P_{\parallel}(\mathbf{k}) + \left( \frac{k^i k^j}{k^2} - \frac{1}{3} \delta^{ij} \right) P_{\perp}(\mathbf{k}) + i \epsilon^{ijk} \frac{k^k}{k} P_{-}(\mathbf{k}). \quad (54) $$

Assuming parity invariance, only the first two terms are present. These are the longitudinal and transverse parts of the parity-even power spectrum matrix. The final term is the parity-odd part of the power spectrum matrix.

One can choose left-/right-handed polarization vectors and form a complete basis {$e_\parallel(\mathbf{k})$, $e_R(\mathbf{k})$, $e_L(\mathbf{k})$} in Fourier space,

$$ V_{ab}(\mathbf{k}) = V_{ab,\parallel}(\mathbf{k}) e_\parallel(\mathbf{k}) + V_{ab,R}(\mathbf{k}) e_R(\mathbf{k}) + V_{ab,L}(\mathbf{k}) e_L(\mathbf{k}). \quad (55) $$

The transverse polarization vectors are defined as eigenvectors of the curl operator,

$$ ik \times e_{R/L}(\mathbf{k}) = \pm k e_{R/L}(\mathbf{k}), \quad (56) $$

and $e_\parallel(\mathbf{k}) = ik/k$. These quantities transform into one another after a parity transformation, as $e_{R/L}(-\mathbf{k}) = \pm ie_{L/R}(\mathbf{k})$. Therefore, one can define the eigenvectors of the parity operator as

$$ e_{\pm}(\mathbf{k}) = \frac{1}{\sqrt{2}} (e_{R}(\mathbf{k}) \pm i e_{L}(\mathbf{k})). \quad (57) $$

Using this alternative basis, the vector field is decomposed as

$$ V_{ab}(\mathbf{k}) = V_{ab,\parallel}(\mathbf{k}) e_{\parallel}(\mathbf{k}) + V_{ab,-}(\mathbf{k}) e_{-}(\mathbf{k}) + V_{ab,+}(\mathbf{k}) e_{+}(\mathbf{k}). \quad (58) $$

The $\pm$ labels of the coefficients are the opposite of the basis labels because the vector must be parity odd, so the coefficients transform with the opposite sign compared to the transverse basis vectors. We can then project the power spectrum matrix onto the polarization vectors,

$$ P_{RR/LL}(\mathbf{k}) = P^{ij}(\mathbf{k}) e_{R/L}^j(\mathbf{k}) e_{R/L}^j(-\mathbf{k}), \quad (59) $$

or

$$ P_{+-}(\mathbf{k}) = P^{ij}(\mathbf{k}) e_{+}^i(\mathbf{k}) e_{-}^j(-\mathbf{k}). \quad (60) $$

It is straightforward to show that

$$ P_{-}(\mathbf{k}) = -\frac{k}{2} P_{\text{vector}}(\mathbf{k}) \quad (61) $$

$$ = \frac{1}{2} (P_{RR}(\mathbf{k}) - P_{LL}(\mathbf{k})) \quad (62) $$

$$ = P_{+-}(\mathbf{k}), \quad (63) $$

where $P_{\text{vector}}(\mathbf{k})$ is the vector POP spectrum.

## 5.2 Higher-order derivatives: angular dependence

The vector and scalar POP spectra that we have considered thus far are minimal in the sense that they involve the minimum number of spatial derivatives required to form the parity-odd triple product. We are free to add any number of spatial derivatives to these constructions, forming higher-rank tensor POP spectra. For example, we consider adding an additional derivative to the vector and pseudovector fields from the vector POP spectrum, as

$$ H_{\text{lab}}^{ij}(\mathbf{x}) = \nabla^i \Phi_a(\mathbf{x}) \nabla^j \Phi_b(\mathbf{x}), \quad (64) $$

$$ H_{2cd}^{ij}(\mathbf{x}) = \epsilon^{ikl} \nabla^j \nabla^k \Phi_c(\mathbf{x}) \nabla^l \Phi_d(\mathbf{x}). \quad (65) $$
The POP spectrum formed by contracting the indices of these two fields $\langle H_{1ab}^{ij}(k) H_{2cd}^{ij}(k') \rangle$ has a similar form to the vector POP spectrum in equation (15), the difference being that the integrand would now contain an additional factor of $q_1 \cdot q_3$. Therefore, this POP tensor spectrum probes a different form of angular dependence in the parity-odd trispectrum.

We can also define higher-rank tensor POP spectra by contracting $\langle H_{1ab}^{li}(k) H_{2cd}^{lj}(k') \rangle$. We can then project this onto the basis tensors

$$e_{R/L}^{ij}(k) = e_{R/L}^{i}(k) e_{R/L}^{j}(k), \qquad (66)$$

where $e_{R/L}^{i}(k)$ are components of the polarization vectors, defined in equation (56). This decomposition gives the right-handed and left-handed components of the transverse tensor spectra, and the difference between these is a POP spectrum. This kind of construction can be easily generalized to higher-derivative and higher-rank tensor POP spectra. Some specific examples of these constructions can be found in Jeong & Kamionkowski (2012) and Shim et al. (2024). Similar constructions for the scalar POP spectra are also possible.

## 5.3 Higher-order constructions: parity-odd bispectra and beyond

The scalar bispectrum is defined by the three-point correlation function in Fourier space through

$$\langle \Phi(k_1) \Phi(k_2) \Phi(k_3) \rangle = (2\pi)^3 \delta_D^{(3)}(k_1 + k_2 + k_3) B(k_1, k_2, k_3). \qquad (67)$$

The Dirac delta function enforces statistical translational invariance and makes the bispectrum a function over 2D, closed triangles. Since in three dimensions, we can rotate a parity-transformed 2D triangle back to its initial configuration, the bispectrum is insensitive to parity under statistical isotropy.

We construct the parity-odd bispectrum that compresses the parity-odd four-point function in a similar way to the vector POP spectra,

$$\langle A_{ab}(k_1) \cdot i k_2 \Phi(k_2) \Phi(k_3) \rangle = (2\pi)^3 \delta_D^{(3)}(k_1 + k_2 + k_3) B_{\text{vector}}(k_1, k_2, k_3). \qquad (68)$$

In this case, we have used the gradient of the scalar rather than the composite vector field. With the composite vector field, this construction would involve five scalar fields corresponding to a compression of the parity-odd five-point function. Similarly, the pseudoscalar requires a product of three scalar fields for its constructions, so the parity-odd bispectrum of one composite pseudoscalar with two scalar fields compresses the parity-odd five-point statistics. Additionally, we can correlate the vector-pseudovector inner product with two pseudoscalar fields or correlate three pseudoscalar fields. These constructions are of order six in the underlying scalar field, compressing the parity-odd six-point function. Thus, equation (68) is the only form for a parity-odd compression of the trispectrum to a bispectrum.

## 6 CONCLUSION

Exploring parity violation on large cosmological scales will provide crucial insights into the fundamental symmetries governing the origins and evolution of the Universe. The search for a parity-odd signal within galaxy surveys presents a key challenge in observational cosmology, especially in light of recent hints of a parity-odd signal in the four-point correlation function of galaxies (Philcox 2022; Hou et al. 2023b), both of which used the isotropic basis functions of Cahn & Slepian (2023), the method proposed in Cahn et al. (2023), and the covariance template of Hou et al. (2022). Direct measurements of four-point correlations pose several technical challenges related to the high dimensionality of the data vector and the resulting size of the covariance matrix, making it difficult to interpret the signal and assess its statistical significance. Motivated by these challenges, we developed a novel set of statistics for parity violation.

The new observables are two-point correlators calculated between composite fields derived via non-linear transformations of the original scalar field. These constructions compress the trispectrum information into a power spectrum, resulting in a significant simplification compared to the full four-point correlation function. The advantages of this approach are two-fold. First, the computation of a two-point function is faster because of the significant dimensionality reduction of the data vector; our computation on a 512³ grid takes fewer than 10 s for each spectrum on 128 cores. This reduction facilitates a more efficient and interpretable analysis of the underlying parity-odd signal and its covariance. Secondly, formulating the estimators in Fourier space enables efficient characterization of the scale dependence of the signal and its soft limits.

We defined two sets of estimators: one based on correlating a vector and a pseudovector field and another correlating a scalar and a pseudoscalar field, both derived from the original scalar field. These two estimators are distinct compressions of the trispectrum and thus carry complementary information on the parity-odd signal. In particular, they are sensitive to different soft limits of the trispectrum. To empirically validate the effectiveness of these estimators, we tested them on a set of mock simulations with an injected parity-odd four-point correlation function. On these mocks, we compared the result of the estimator with semi-analytical theoretical calculations. The agreement between the data and our theoretical calculations demonstrates that our estimators robustly capture the injected parity-odd signal.

Beyond validation, we investigated the sensitivity of our new estimators by computing the correlation matrix on Gaussian realizations. We found that different bins of our estimator are predominantly uncorrelated. Furthermore, a $\chi^2$-squared test on simulated parity-violating mocks highlighted the substantial statistical significance of both estimators, with the scalar variant exhibiting enhanced sensitivity compared to its vector counterpart. This sensitivity is specific to the choice of the primordial parity-violating trispectrum template and filtering functions. In particular, we did not attempt to optimize the filtering functions by maximizing signal-to-noise. Such an optimization could dramatically increase the significance. We plan to explore this in future work.
Applying our estimator to realistic survey data requires several developments. Surveys observe galaxies within a specific survey mask and selection function. We must incorporate these into our estimator to remove the spurious clustering signals they would otherwise produce. Moreover, galaxies are discrete tracers, and this results in shot noise that should be characterized and subtracted from the POP spectra. We leave these developments for future work. After overcoming these challenges, we expect the POP estimators will emerge as an essential tool in the search for parity violation in current and future CMB and galaxy surveys.

## ACKNOWLEDGEMENTS

We thank Toshiki Kurita, Fabian Schmidt, Azadeh Moradinezhad Dizgah, Robert Cahn, Donghui Jeong, and Ue-Li Pen for helpful discussions. AC acknowledges funding support from the Initiative Physique des Infinis (IPI), a research training program of the Idex SUPER at Sorbonne Université. JH has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 101025187. This work was supported in part by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy: Grant No. EXC-2094–390783311. This work has also received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 101007633. The Kavli IPMU is supported, in part, by World Premier International Research Center Initiative (WPI), MEXT, Japan. ZS thanks the Max Planck Institute for Astrophysics for hospitality during some of the period during which this work was performed.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The software used to generate and analyse the data, and the analysis results data in Figs 3, 4, and 5 are publicly available at https://github.com/dsjamieson/parity_violatingTemplates.

## REFERENCES

* Adame A. G. et al., 2024, preprint (arXiv:2404.03002)

* Aghamousa A. et al., 2016, preprint (arXiv:1611.00036)

* Amendola L. et al., 2018, Living Rev. Rel., 21, 2

* Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Phys. Rept., 402, 103

* Cahn R. N., Slepian Z., 2023, J. Phys. A, 56, 325204

* Cahn R. N., Slepian Z., Hou J., 2023, Phys. Rev. Lett., 130, 201002

* Coulton W. R., Philcox O. H. E., Villaescusa-Navarro F., 2024, Phys. Rev. D, 109, 23531

* Doré O. et al., 2014, preprint (arXiv:1412.4872)

* Garwin R. L., Lederman L. M., Weinrich M., 1957, Phys. Rev., 105, 1415

* Gualdi D., Verde L., 2022, J. Cosmol. Astropart. Phys., 9, 50

* Gualdi D., Gil-Marín H. e., Verde L., 2021, J. Cosmol. Astropart. Phys., 7, 8

* Hou J., Cahn R. N., Philcox O. H. E., Slepian Z., 2022, Phys. Rev. D, 106, 43515

* Hou J., Moradinezhad Dizgah A., Hahn C., Massara E., 2023a, J. Cosmol. Astropart. Phys., 3, 45

* Hou J., Slepian Z., Cahn R. N., 2023b, MNRAS, 522, 5701

* Hu W., 2001, Phys. Rev. D, 64, 83005

* Ivezić v. et al., 2019, ApJ, 873, 111

* Jeong D., Kamionkowski M., 2012, Phys. Rev. Lett., 108, 251301

* Lee T. D., Yang C.-N., 1956, Phys. Rev., 104, 254

* Moradinezhad Dizgah A., Lee H., Schmittfull M., Dvorkin C., 2020, J. Cosmol. Astropart. Phys., 4, 11

* Munshi D., Lee H., Dvorkin C., McEwen J. D., 2022, J. Cosmol. Astropart. Phys., 11, 20

* Philcox O. H. E., 2022, Phys. Rev. D, 106, 63501

* Philcox O. H. E., 2023, Phys. Rev. Lett., 131, 181001

* Philcox O. H. E., Ereza J., 2024, preprint (arXiv:2401.09523)

* Philcox O. H. E., Shiraishi M., 2024, Phys. Rev. D, 109, 63522

* Sabiu C. G., Hoyle B., Kim J., Li X.-D., 2019, ApJS, 242, 29

* Schmittfull M., Baldauf T., Seljak U., 2015, Phys. Rev. D, 91, 43530

* Schmittfull M., Moradinezhad Dizgah A., 2021, J. Cosmol. Astropart. Phys., 3, 20

* Shim J., Pen U.-L., Yu H.-R., Okumura T., 2024, preprint (arXiv:2406.06080)

* Shiraishi M., 2016, Phys. Rev. D, 94, 83503

* Slepian Z., Eisenstein D. J., 2015, MNRAS, 454, 4142

* Slepian Z. et al., 2017, MNRAS, 468, 1070

* Sunseri J., Slepian Z., Portillo S., Hou J., Kahraman S., Finkbeiner D. P., 2023, RAS Techniques Instruments, 2, 62

* Takada M. et al., 2014, Publ. Astron. Soc. Jap., 66, R1

* Wang Y. et al., 2022, ApJ, 928, 1

* Wu C. S., Ambler E., Hayward R. W., Hoppes D. D., Hudson R. P., 1957, Phys. Rev., 105, 1413

* Xu X., Padmanabhan N., Eisenstein D. J., Mehta K. T., Cuesta A. J., 2012, MNRAS, 427, 2146
APPENDIX A: THE BINDED PARITY-ODD TRISPECTRUM AND THE BINDED POP SPECTRA

Measuring the parity-odd trispectrum from data requires choosing a binning scheme for the tetrahedral configurations of wave vectors. To form a bin of nearby tetrahedra, we specify the side lengths $k_1, k_2, k_3k_4$, and the diagonal lengths $K$ and $\tilde{K}$. We then select the set of wave vectors satisfying

$$k_1 + k_2 + K = 0, \qquad (A1)$$

$$k_3 + k_4 - K = 0, \qquad (A2)$$

$$k_1 + k_4 + \tilde{K} = 0. \qquad (A3)$$

We define the 6D binning scheme through the integral

$$\int_{\text{bin}} \equiv \int_{\mathbf{q}_1 \in k_1} \int_{\mathbf{q}_2 \in k_2} \int_{\mathbf{q}_3 \in k_3} \int_{\mathbf{q}_4 \in k_4} \Theta(4(|\mathbf{q}_1 + \mathbf{q}_2| - K)^2 - \Delta k^2) \Theta(4(|\mathbf{q}_1 + \mathbf{q}_4| - \tilde{K})^2 - \Delta k^2). \qquad (A4)$$

The integral notation is defined in Table A1. The Heaviside functions enforce $K - \Delta k/2 \le |\mathbf{q}_1 + \mathbf{q}_2| \le K + \Delta k/2$ and $\tilde{K} - \Delta k/2 \le |\mathbf{q}_1 + \mathbf{q}_4| \le \tilde{K} + \Delta k/2$, which defines the binning scheme for the diagonals.

We isolate the binned parity-odd trispectrum by integrating the product of the triple product and the four-point expectation value over the bin,

$$\bar{\tau}_{-}(k_1, k_2, k_3, k_4, K, \tilde{K}) = \mathcal{N}_{\tau_{-}}^{-1} \int_{\text{bin}} -i \mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3) \langle \Phi(\mathbf{q}_1) \Phi(\mathbf{q}_2) \Phi(\mathbf{q}_3) \Phi(\mathbf{q}_4) \rangle \qquad (A5)$$

$$= \mathcal{N}_{\tau_{-}}^{-1} \int_{\text{bin}} (2\pi)^3 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{q}_3 + \mathbf{q}_4) [\mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3)]^2 \tau_{-}(q_1, q_2, q_3, q_4, |\mathbf{q}_1 + \mathbf{q}_2|, |\mathbf{q}_1 + \mathbf{q}_4|) . \qquad (A6)$$

Here, $\mathcal{N}_{\tau_{-}}$ is a normalization factor, which we will define later. The minus sign in the first line is for convenience, and we could absorb it into the normalization factor. We split the Dirac delta function as follows:

$$(2\pi)^3 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{q}_3 + \mathbf{q}_4) = (2\pi)^9 \int_{\mathcal{Q}, \tilde{\mathcal{Q}}} \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 + \mathcal{Q}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathcal{Q}) \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_4 + \tilde{\mathcal{Q}}) \times [\mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3)]^2 \tau_{-}(q_1, q_2, q_3, q_4, \mathcal{Q}, \tilde{\mathcal{Q}}) . \qquad (A7)$$

The right-hand side of this expression imposes the condition that the tetrahedron consists of two triangles, $\{\mathbf{q}_1, \mathbf{q}_2, -\mathcal{Q}\}$ and $\{\mathbf{q}_3, \mathbf{q}_4, \mathcal{Q}\}$, which are joined along their common side of length $\mathcal{Q}$. The vectors $\mathbf{q}_1$ and $\mathbf{q}_4$ are additionally constrained to form a triangle $\{\mathbf{q}_1, \mathbf{q}_4, -\tilde{\mathcal{Q}}\}$. After substituting equation (A7) into equation (A6), the Heaviside functions limit the integrations over $\mathcal{Q}$ and $\tilde{\mathcal{Q}}$. Then, with the integral notation in Table A1, the binned parity-odd trispectrum is given by

$$\bar{\tau}_{-}(k_1, k_2, k_3, k_4, K, \tilde{K}) = \mathcal{N}_{\tau_{-}}^{-1} \int_{6D-\text{bin}} [(2\pi)^9 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 + \mathcal{Q}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathcal{Q}) \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_4 + \tilde{\mathcal{Q}}) \times [\mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3)]^2 \tau_{-}(q_1, q_2, q_3, q_4, \mathcal{Q}, \tilde{\mathcal{Q}})] . \qquad (A8)$$

<table><caption>**Table A1.** Fourier-space integral and harmonic-space summation notation.</caption><thead><tr><th>Notation</th><th>Definition</th></tr></thead><tbody><tr><td>∫<sub>q∈k</sub></td><td>∫<sub>k-Δk/2</sub><sup>k+Δk/2</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">k<sup>2</sup>dK</div><div style="padding: 0 0.2em;">2π<sup>2</sup></div></div></td></tr><tr><td>∫<sub>q̃</sub></td><td>∫ <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">dΩ<sub>q̃</sub></div><div style="padding: 0 0.2em;">4π</div></div></td></tr><tr><td>∫<sub>q∈k</sub></td><td>∫<sub>q̃</sub> ∫<sub>q∈k</sub></td></tr><tr><td>∫<sub>6D-bin</sub></td><td>∫<sub>q<sub>1</sub>∈k<sub>1</sub></sub> ∫<sub>q<sub>2</sub>∈k<sub>2</sub></sub> ∫<sub>q<sub>3</sub>∈k<sub>2</sub></sub> ∫<sub>q<sub>4</sub>∈k<sub>4</sub></sub> ∫<sub>Q∈K</sub> ∫<sub>Q̃∈K</sub></td></tr><tr><td>∫<sub>5D-bin</sub></td><td>∫<sub>q<sub>1</sub>∈k<sub>1</sub></sub> ∫<sub>q<sub>2</sub>∈k<sub>2</sub></sub> ∫<sub>q<sub>3</sub>∈k<sub>3</sub></sub> ∫<sub>q<sub>4</sub>∈k<sub>4</sub></sub> ∫<sub>Q∈K</sub></td></tr><tr><td>∫<sub>5D-bin</sub></td><td>∫<sub>q<sub>1</sub>∈k<sub>1</sub></sub> ∫<sub>q<sub>2</sub>∈k<sub>2</sub></sub> ∫<sub>q<sub>3</sub>∈k<sub>3</sub></sub> ∫<sub>q<sub>4</sub>∈k<sub>4</sub></sub> ∫<sub>Q∈K</sub></td></tr><tr><td>∑<sub>L<sub>i</sub>M<sub>i</sub></sub></td><td>∞ ∑<sub>L<sub>1</sub>=0</sub> ∑<sub>M<sub>1</sub>=-L<sub>1</sub></sub> ∞ ∑<sub>L<sub>2</sub>=0</sub> ∑<sub>M<sub>2</sub>=-L<sub>2</sub></sub> ∞ ∑<sub>L<sub>3</sub>=0</sub> ∑<sub>M<sub>3</sub>=-L<sub>3</sub></sub></td></tr></tbody></table>

The notation $q \in k$ indicates a wave vector magnitude belonging to a bin centered on $k$ with width $\Delta k$. The integral over $d\Omega_{q\hat{q}}$ is the solid angle integral over the azimuthal and polar angles associated with the wave vector $q$ in polar coordinates.
Unfortunately, the full 6D trispectrum is difficult to analyse directly. The first two Dirac delta functions factor into integrals over two triangles, but the third couples these together, making the whole integral not factorizable. For this reason, we extend the $\tilde{Q}$ bin to be unrestricted, $\tilde{Q} \in [0, \infty)$, which fully averages over the second diagonal. An alternative approach is to implement isotropic basis functions, which discretize the angular dependence through spherical harmonic expansions (Cahn & Slepian 2023), but we will not pursue that here.

The 5D binned trispectrum, with notation from Table A1, is

$$ \bar{\tau}_{-}(k_1, k_2, k_3, k_4, K) = \mathcal{N}_{\tau_{-}}^{-1} \int_{5D\text{-bin}} (2\pi)^6 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{Q}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathbf{Q}) [\mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3)]^2 \tau_{-}(q_1, q_2, q_3, q_4, \mathcal{Q}, |\mathbf{q}_1 + \mathbf{q}_4|), \quad (A9) $$

The third Dirac delta function in equation (A8) enforced $\tilde{Q} = |\mathbf{q}_1 + \mathbf{q}_4|$. We choose the somewhat arbitrary normalization factor $\mathcal{N}_{\tau_{-}}$ to be

$$ \mathcal{N}_{\tau_{-}} = \int_{5D\text{-bin}} (2\pi)^9 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{Q}) \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 - \mathbf{Q}) [\mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3)]^2, \quad (A10) $$

so equation (A9) is the weighted average of $\tau_{-}(q_1, q_2, q_3, q_4, \mathcal{Q}, |\mathbf{q}_1 + \mathbf{q}_4|)$ within the bin.

Expressing both Dirac delta functions as the Fourier transform of plane waves enables us to compute most angular integrals analytically. The calculation simplifies greatly by noticing that we can rewrite the triple product as $\mathbf{q}_1 \cdot (\mathbf{q}_2 \times \mathbf{q}_3) = -\mathbf{Q} \cdot (\mathbf{q}_1 \times \mathbf{q}_3)$, which is more symmetric between the two triangles forming the tetrahedron. However, if the trispectrum depends on the diagonals, the angular dependence on $\mathbf{q}_1 \cdot \mathbf{q}_4$ may prevent us from analytically computing all angular integrals. From here on, we will restrict to the case where the trispectrum has no explicit dependence on the diagonals.

The spherical harmonic expansion of the squared triple product is

$$ [\mathbf{Q} \cdot (\mathbf{q}_1 \times \mathbf{q}_3)]^2 = -6(4\pi)^3 (Q q_1 q_3)^2 \sum_{L_i, M_i} D_{L_1 L_2 L_3}^{M_1 M_2 M_3} Y_{L_1 M_1}(\hat{q}_1) Y_{L_2 M_2}(\hat{q}_3) Y_{L_3 M_3}(\hat{Q}), \quad (A11) $$

where the coefficients are

$$ D_{L_1 L_2 L_3}^{M_1 M_2 M_3} = \begin{pmatrix} L_1 & L_2 & L_3 \\ M_1 & M_2 & M_3 \end{pmatrix} \begin{Bmatrix} L_1 & L_2 & L_3 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{Bmatrix} \prod_{i=1}^{3} \sqrt{\frac{2L_i + 1}{4\pi}} \begin{pmatrix} 1 & 1 & L_i \\ 0 & 0 & 0 \end{pmatrix}. \quad (A12) $$

Here, the 2 by 3 matrices are Wigner 3-j symbols, and the 3 by 3 matrix is a Wigner 9-j symbol. The form of the coefficients in the product restricts $L_i \in \{0, 2\}$.

The trispectrum estimator is an integral over tetrahedra formed by two triangles: $\{\mathbf{q}_1, \mathbf{q}_2, -\mathbf{Q}\}$ and $\{\mathbf{q}_3, \mathbf{q}_4, \mathcal{Q}\}$. By taking one vector from each triangle and the common vector shared by both triangles in the triple product, we preserve the symmetry between the two pairs $(\mathbf{q}_1, \mathbf{q}_2)$ and $(\mathbf{q}_3, \mathbf{q}_4)$, greatly simplifying the calculation. Each triangle contributes a Dirac delta function and a $Y_{LM}$ from a vector that is not $\pm \mathbf{Q}$. Expressing the Dirac delta functions as Fourier transforms of plane waves allows us to compute all angular integrals,

$$ \int_{\mathbf{q}_1 \in k_1, \mathbf{q}_2 \in k_2} (2\pi)^3 \delta_D^{(3)}(\mathbf{q}_1 + \mathbf{q}_2 - \mathbf{Q}) Y_{L_1 M_1}(\hat{q}_1) = Y_{L_1 M_1}(\hat{Q}) \int_{\mathbf{q}_1 \in k_1, \mathbf{q}_2 \in k_2} 4\pi \int_0^\infty dx \, x^2 j_{L_1}(q_1 x) j_{0}(q_2 x) j_{L_1}(Qx) \quad (A13) $$

$$ = Y_{L_1 M_1}(\hat{Q}) \int_{\mathbf{q}_1 \in k_1, \mathbf{q}_2 \in k_2} I_{L_1 0 L_1}(\mathbf{q}_1, \mathbf{q}_2, \mathcal{Q}), \quad (A14) $$

and similarly,

$$ \int_{\mathbf{q}_3 \in k_3, \mathbf{q}_4 \in k_4} (2\pi)^3 \delta_D^{(3)}(\mathbf{q}_3 + \mathbf{q}_4 + \mathbf{Q}) Y_{L_2 M_2}(\hat{q}_3) = (-1)^{L_2} Y_{L_2 M_2}(\hat{Q}) \int_{\mathbf{q}_3 \in k_3, \mathbf{q}_4 \in k_4} 4\pi \int_0^\infty dx \, x^2 j_{L_2}(q_3 x) j_{0}(q_4 x) j_{L_2}(Qx) \quad (A15) $$

$$ = (-1)^{L_2} Y_{L_2 M_2}(\hat{Q}) \int_{\mathbf{q}_3 \in k_3, \mathbf{q}_4 \in k_4} I_{L_2 0 L_2}(\mathbf{q}_3, \mathbf{q}_4, \mathcal{Q}). \quad (A16) $$

Since $L_2 \in \{0, 2\}$, we have $(-1)^{L_2} = 1$, so we can drop this factor. Both triangles contribute factors to the integrand that have the same form, illustrating the symmetric roles of the two triangles that form the tetrahedron. The triple-spherical Bessel integrals are given by

$$ I_{L_1 0 L_1}(\mathbf{q}_1, \mathbf{q}_2, \mathcal{Q}) = 4\pi \int_0^\infty dx \, x^2 j_{L_1}(\mathbf{q}_1 x) j_0(\mathbf{q}_2 x) j_{L_1}(Qx) \quad (A17) $$

$$ = \frac{\pi^2}{q_1 q_2 Q} \Theta(\sin^2 \theta_{12}) P_{L_1}(\cos \theta_{12}). \quad (A18) $$

Here, $\theta_{12}$ is the angle between $\mathbf{q}_1$ and $\mathbf{Q}$. As discussed in the main text, the Heaviside function enforces the triangle inequalities that ensure $\mathbf{q}_1$, $\mathbf{q}_2$, and $-\mathbf{Q}$ can form a closed triangle. $P_{L_1}$ is the $L_1^{th}$ Legendre polynomial. The only relevant ones will be

$$ P_0(x) = 1, \quad P_2(x) = \frac{1}{2}(3x^2 - 1). \quad (A19) $$

At this point, three $Y_{L_i M_i}$ products remain, all with the argument $\hat{Q}$. Integrating over $\hat{Q}$ gives Gaunt's integral (divided by $4\pi$),<fnref n="3" />

$$ G_{L_1 L_2 L_3}^{M_1 M_2 M_3} = \int_{\mathcal{Q}} Y_{L_1 M_1}(\hat{Q}) Y_{L_2 M_2}(\hat{Q}) Y_{L_3 M_3}(\hat{Q}) \quad (A20) $$

<footnote n="3">³https://dlmf.nist.gov/34.3, equation (34.3.22).</footnote>
$$ = \begin{pmatrix} L_1 L_2 L_3 \\ M_1 M_2 M_3 \end{pmatrix} \begin{pmatrix} L_1 L_2 L_3 \\ 000 \end{pmatrix} \prod_{i=1}^{3} \sqrt{\frac{2L_i + 1}{4\pi}}. \qquad (A21) $$

The form of the trispectrum after computing all angular integrals is

$$ \bar{\tau}_{-}(k_1, k_2, k_3, k_4, K) = -6 \mathcal{N}_{\tau_{-}}^{-1} \int_{5D-\text{bin}} (Q q_1 q_3)^2 \sum_{L_i M_i} D_{L_1 L_2 L_3}^{M_1 M_2 M_3} G_{L_1 L_2 L_3}^{M_1 M_2 M_3} I_{L_1 0 L_1} (q_1, q_2, Q) I_{L_2 0 L_2} (q_3, q_4, Q) \tau_{-}(q_1, q_2, q_3, q_4), \qquad (A22) $$

where the summation notation is defined in Table A1. Using the orthogonality relation for the 3-j symbol<fnref n="4" />

$$ \sum_{M_1=-L_1}^{L_1} \sum_{M_2=-L_2}^{L_2} \sum_{M_3=-L_3}^{L_3} \left( \begin{pmatrix} L_1 & L_2 & L_3 \\ M_1 & M_2 & M_3 \end{pmatrix} \right)^2 = 1, \qquad (A23) $$

we can perform the sums over all $M_i$. Since neither triple-spherical Bessel integral involves $L_3$, we can also sum over it, which gives

$$ \bar{\tau}_{-}(k_1, k_2, k_3, k_4, K) = 6 \mathcal{N}_{\tau_{-}}^{-1} \int_{5D-\text{bin}} (Q q_1 q_3)^2 \sum_{L_1, L_2} H_{L_1 L_2} I_{L_1 0 L_1} (q_1, q_2, Q) I_{L_2 0 L_2} (q_3, q_4, Q) \tau_{-}(q_1, q_2, q_3, q_4), \qquad (A24) $$

where the coefficients now simplify considerably,

$$ \begin{aligned} H_{L_1 L_2} &= -\sum_{L_3} \begin{pmatrix} L_1 & L_2 & L_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{Bmatrix} L_1 & L_2 & L_3 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{Bmatrix} \prod_{i=1}^{3} (2L_i + 1) \begin{pmatrix} 1 & 1 & L_i \\ 0 & 0 & 0 \end{pmatrix} \\ &= \frac{i^{L_1+L_2}}{27}. \end{aligned} \qquad (A25) $$

The two remaining sums factorize,

$$ \bar{\tau}_{-}(k_1, k_2, k_3, k_4, K) = \frac{2}{9} \mathcal{N}_{\tau_{-}}^{-1} \int_{5D-\text{bin}} (Q q_1 q_3)^2 \left[ \sum_{L_1} i^{L_1} I_{L_1 0 L_1} (q_1, q_2, Q) \right] \left[ \sum_{L_2} i^{L_2} I_{L_2 0 L_2} (q_3, q_4, Q) \right] \tau_{-}(q_1, q_2, q_3, q_4). \qquad (A27) $$

Each sum only has two terms, so we can sum them explicitly and expand the Legendre polynomials,

$$ \bar{\tau}_{-}(k_1, k_2, k_3, k_4, K) = \frac{\pi^4}{2} \mathcal{N}_{\tau_{-}}^{-1} \int_{5D-\text{bin}} \frac{q_1 q_3}{q_2 q_4} \Theta(\sin^2 \theta_{12}) \sin^2 \theta_{12} \Theta(\sin^2 \theta_{34}) \sin^2 \theta_{34} \tau_{-}(q_1, q_2, q_3, q_4), \qquad (A28) $$

where we set $\tau_{-} = 1$ to obtain the normalization factor,

$$ \mathcal{N}_{\tau_{-}} = \frac{\pi^4}{2} \int_{5D-\text{bin}} \frac{q_1 q_3}{q_2 q_4} \Theta(\sin^2 \theta_{12}) \sin^2 \theta_{12} \Theta(\sin^2 \theta_{34}) \sin^2 \theta_{34}. \qquad (A29) $$

Instead of restricting the integration bounds, we could define filter functions with Heaviside functions that impose the binning. Then, the binned parity-odd trispectrum takes the form,

$$ \bar{\tau}_{-}(k_1, k_2, k_3, k_4, K) = \frac{\pi^4}{2} \mathcal{N}_{\tau_{-}}^{-1} \int_{q_1, q_2, q_3, q_4, Q} \left[ \frac{q_1 q_3}{q_2 q_4} \Theta(\sin^2 \theta_{12}) \sin^2 \theta_{12} \Theta(\sin^2 \theta_{34}) \sin^2 \theta_{34} \right. \\ \left. \times f_a(q_1) f_b(q_2) f_c(q_3) f_d(q_4) f_e(Q) \tau_{-}(q_1, q_2, q_3, q_4) \right]. \qquad (A30) $$

By taking all but $f_e$ to be generic filter functions, rather than just Heaviside functions selecting spherical shells, and taking the thin-bin limit of $f_e$, this has the same form as the vector POP spectrum in equation (17). Similarly, by taking all but $f_a$ to be generic filter functions and taking the thin-bin limit of $f_a$, this has the same form as the scalar POP spectrum in equation (27). Thus, computing a POP spectrum is equivalent to computing a 1D subset of bins for the full 5D trispectrum.

<footnote n="4">⁴https://dlmf.nist.gov/34.3 equation (34.3.18).</footnote>

<footnote>This paper has been typeset from a TeX/LATEX file prepared by the author.</footnote>