SARABANDE: 3/4 Point Correlation Functions with Fast Fourier
Transforms–Appendix

#  SARABANDE: 3/4 Point Correlation Functions with Fast Fourier Transforms

James Sunseri1 E-mail<EMAIL_ADDRESS>(JS) Zachary Slepian2 E-mail:
<EMAIL_ADDRESS>(ZS) Stephen Portillo3 E-mail<EMAIL_ADDRESS>(SP) Jiamin
Hou2,4 E-mail<EMAIL_ADDRESS>(JH) Sule Kahraman5 & Douglas P. Finkbeiner6

1Department of Astronomy E-mail<EMAIL_ADDRESS>(SK)E-mail:
<EMAIL_ADDRESS>(DPF) University of California  Berkeley  CA
94720-3411  USA
2Department of Astronomy  University of Florida  211 Bryant Space Science
Center  Gainesville  FL 32611  USA
3DIRAC Institute  Department of Astronomy  University of Washington  Seattle
WA 98195-1580  USA
4Max-Planck-Institut für Extraterrestrische Physik  Postfach 1312
Giessenbachstrasse 1  85748 Garching bei München  Germany
5Massachusetts Institute of Technology  Department of Electrical Engineering
and Computer Science  77 Massachusetts Avenue  Cambridge  MA 02139  USA
6Center for Astrophysics | Harvard & Smithsonian  60 Garden Street  Cambridge  MA 02138  USA

(Accepted XXX. Received YYY; in original form ZZZ; 2022)

###### Abstract

We present a new python package sarabande for measuring 3 & 4 Point
Correlation Functions (3/4 PCFs) in $\mathcal{O}(N_{\mathrm{g}}\log
N_{\mathrm{g}})$ time using Fast Fourier Transforms (FFTs), with
$N_{\mathrm{g}}$ the number of grid points used for the FFT. sarabande can
measure both projected and full 3 and 4 PCFs on gridded 2D and 3D datasets.
The general technique is to generate suitable angular basis functions on an
underlying grid, radially bin these to create kernels, and convolve these
kernels with the original gridded data to obtain expansion coefficients about
every point simultaneously. These coefficients are then combined to give us
the 3/4 PCF as expanded in our basis. We apply sarabande to simulations of the
Interstellar Medium (ISM) to show the results and scaling of calculating both
the full and projected 3/4 PCFs.

###### keywords:

methods: data analysis, statistical, numerical

##  1 Introduction

Extracting structural and physical information from density and intensity maps
is a common goal across multiple fields in astrophysics (Peebles & Ratra,
2003). Correlation functions are a natural means to characterize the
translation and rotation-invariant information in a density field. The 2-Point
Correlation Function (2PCF) or its Fourier-space analog the power spectrum are
ubiquitous, and quantify correlations between pairs of points as a function of
separation or Fourier-space wave-vector. However, for density fields beyond a
Gaussian random field, they do not capture all of the information (e.g.
Samushia et al. 2021).

To learn more about a given density field, one can look to higher-order
statistics (Peebles, 2001). The 3-Point Correlation Function (3PCF) and the
4-Point Correlation Function (4PCF) or their Fourier-space counterparts the
bispectrum and trispectrum are the next statistics to consider (Fry & Peebles,
1978; Collis et al., 1998; Burkhart & Lazarian, 2010; Koch et al., 2019; Sabiu
et al., 2019; Gualdi et al., 2021). These higher-order statistics measure
correlations between triplets and quadruplets of points as a function of
triangular and tetrahedral geometries in real or Fourier space respectively.
They are useful for tracing non-Gaussian signals as well as probing parity-
violation (Cahn et al., 2021; Hou et al., 2022b). The biggest hurdle for using
these higher-order statistics has been the computational complexity. A naive
calculation of an NPCF scales as $\mathcal{O}(N_{\mathrm{p}}^{\mathrm{N}})$
with $N_{\mathrm{p}}$ data points; this is infeasible for the ever-growing
scale of astrophysical datasets.

Recently a series of papers (Slepian & Eisenstein, 2015a, 2016, 2018; Philcox
et al., 2022; Philcox & Slepian, 2022; Hou et al., 2022a) have presented an
approach to compute the 3PCF, 4PCF, 5PCF, and 6PCF that formally scales as
$\mathcal{O}(N_{\mathrm{p}}^{2})$ for $N_{\mathrm{p}}$ particles. It can be
further accelerated using FFTs to $\mathcal{O}(N_{\rm g}\log N_{\rm g})$ with
$N_{\rm g}$ the number of grid points used for the FFTs. In this work we
present sarabande, a code package that computes the 3/4 PCFs using FFTs to
obtain the coefficients needed in our estimator at all points simultaneously
based on Slepian & Eisenstein (2015a) and Philcox et al. (2022).111Spherical-
harmonic-based methods like those presented here have also been developed for
the anisotropic 2-Point Correlation Function in Slepian & Eisenstein (2015b),
Philcox & Slepian (2021) and for the power spectrum in Hand et al. (2017); the
base 3PCF algorithm of Slepian & Eisenstein 2015a is also implemented in
python in nbodykit (Hand et al., 2018). This method scales as
$\mathcal{O}(N_{\rm g}\log N_{\rm g})$. sarabande is especially enabling for
analysis of simulations that are already on regular grids, as FFTs require
gridding but in this case, that will be lossless. In principle sarabande could
be used on any dataset as long as one is willing to grid it, and there are
well-known methods to correct any artifacts introduced by gridding (e.g. Jing
2005). We also present a version for 2D fields, suitable for computing the 3/4
PCF either on slices of gridded data or on a projected density field where one
of the three axes has been integrated out.

This paper is structured as follows: Section 2 describes the underlying
derivations and mathematical reasoning behind the FFT based approach to
measuring the 3/4 PCFs. Section 3 outlines how the code is implemented in
python and how to use it. Section 4 discusses the overall performance of the
algorithms and our process of validating the code. Section 5 offers potential
future applications of sarabande. Section 6 concludes. An Appendix describes
the algorithms used in sarabande.

##  2 The Algorithms

We here outline the derivations for the algorithms underlying sarabande. Some
of the algorithms have been presented in previous works (Portillo et al.,
2018; Philcox et al., 2022; Saydjari et al., 2021) but we present them here
for completeness.

###  2.1 Full 3PCF

Here we briefly recapitulate the mathematical structure of the full 3PCF
algorithm. More details of the base algorithm are in Slepian & Eisenstein
(2015a), and its implementation using FFTs in Slepian & Eisenstein (2016).

The full 3PCF can be parameterized by two sides of a triangle $r_{1}$ and
$r_{2}$ and the cosine of their enclosed angle, $\hat{r}_{1}\cdot\hat{r}_{2}$.
We can expand the full 3PCF, denoted $\zeta$, as a series of radial
coefficients $\zeta_{\ell}$ dependent on $r_{1}$ and $r_{2}$ times the
“isotropic” basis functions $\mathcal{P}_{\ell}(\hat{r}_{1},\hat{r}_{2})$ of
Cahn & Slepian (2020), which are orthonormal and up to a phase and rescaling
correspond to Legendre
polynomials.222$\mathcal{P}_{\ell}(x)=(-1)^{\ell}\sqrt{2\ell+1}/(4\pi)\mathcal{L}_{\ell}(x)$,
with $\mathcal{L}_{\ell}$ the Legendre polynomial of order $\ell$. We have

$\displaystyle\zeta(r_{1},r_{2};\hat{r}_{1}\cdot\hat{r}_{2})=\sum_{\ell=0}^{\infty}\zeta_{\ell}(r_{1},r_{2})\mathcal{P}_{\ell}(\hat{r}_{1}\cdot\hat{r}_{2}).$
(1)

We can form the full 3PCF as a volume-average ($V$ denotes volume) of “local”
estimates of the full 3PCF (indicated by a hat) about points $\vec{x}$:

$\displaystyle\zeta(r_{1},r_{2};\hat{r}_{1}\cdot\hat{r}_{2})=\int\frac{d^{3}\vec{x}}{V}\;\hat{\zeta}(r_{1},r_{2};\hat{r}_{1}\cdot\hat{r}_{2};\vec{x}).$
(2)

In turn, the local full 3PCF is constructed, for a density field $\delta$, as

$\displaystyle\hat{\zeta}(r_{1},r_{2};\hat{r}_{1}\cdot\hat{r}_{2};\vec{x})=\left<\delta(\vec{x})\delta(\vec{x}+\vec{r}_{1})\delta(\vec{x}+\vec{r}_{2})\right>_{\mathcal{R}}$
(3)

where the angle brackets with subscript $\mathcal{R}$ indicate an average over
joint rotations $\mathcal{R}$ of $\hat{r}_{1}$ and $\hat{r}_{2}$ about
$\vec{x}$. This rotation-averaging is lossless under the assumption of
isotropy about $\vec{x}$. Since projection onto the isotropic basis functions
is a linear operation, we then have that

$\displaystyle\zeta_{\ell}(r_{1},r_{2})=\int\frac{d^{3}\vec{x}}{V}\;\hat{\zeta}_{\ell}(r_{1},r_{2};\vec{x}).$
(4)

The relation above means that the overall radial coefficients of the full 3PCF
can be computed as the volume-average of the local estimates. We can write the
local estimate of the multipole coefficients as

$\hat{\zeta}_{\ell}(r_{1},r_{2};\vec{x})=\delta(\vec{x})\int
d\Omega_{1}d\Omega_{2}\;\mathcal{P}_{\ell}^{*}(\hat{r}_{1},\hat{r}_{2})\\\
\times\delta(\vec{x}+\vec{r_{1}})\delta(\vec{x}+\vec{r_{2}}).$ (5)

Using equations 1 and 3 of Cahn & Slepian (2020), the basis functions for the
full 3PCF can be written in terms of spherical harmonics:

$\displaystyle\mathcal{P}_{\ell}(\hat{r}_{1},\hat{r}_{2})=\frac{(-1)^{\ell}}{\sqrt{2\ell+1}}\sum_{m=-\ell}^{\ell}Y_{\ell
m}(\hat{r}_{1})Y_{\ell m}^{*}(\hat{r}_{2}),$ (6)

Combining these results, we can rewrite the local estimate of the multipole
coefficients as

$\hat{\zeta}_{\ell}(r_{1},r_{2};\vec{x})=\frac{(-1)^{\ell}}{\sqrt{2\ell+1}}\delta(\vec{x})\;\\\
\times\sum_{m=-\ell}^{\ell}\int d\Omega_{1}d\Omega_{2}\;Y_{\ell
m}^{*}(\hat{r}_{1})Y_{\ell
m}(\hat{r}_{2})\delta(\vec{x}+\vec{r_{1}})\delta(\vec{x}+\vec{r_{2}}).$ (7)

The angular integrals may then be split; defining our coefficients as

$\displaystyle a_{\ell m}(r_{i};\vec{x})=\int d\Omega_{i}\;Y_{\ell
m}^{*}(\hat{r}_{i})\delta(\vec{x}+\vec{r}_{i}),$ (8)

we find

$\displaystyle\hat{\zeta}_{\ell}(r_{1},r_{2};\vec{x})=\frac{(-1)^{\ell}}{\sqrt{2\ell+1}}\;\delta(\vec{x})\sum_{m=-\ell}^{\ell}a_{\ell
m}(r_{1};\vec{x})a_{\ell m}^{*}(r_{2};\vec{x}).$ (9)

We notice that the key quantity is $a_{\ell m}$ as given by equation 8, and
that it has the structure of a convolution. Furthermore, since $a_{\ell m}$
has this structure, it may be computed around all points $\vec{x}$ at once
using FFTs, leading to an algorithm scaling as $\mathcal{O}(N_{\rm g}\log
N_{\rm g})$ with $N_{\rm g}$ as the number of grid points used for the FFT.

Generally we compute the full 3PCF radial coefficients on bins in $r_{i}$, and
since binning is a linear operation, we can rewrite the $a_{\ell m}$ in
equation 8 as

$\displaystyle a_{\ell m}(r_{i};\vec{x})\to a^{\rm b}_{\ell m}(\vec{x})=\int
d\Omega\;Y_{\ell m}^{*}(\hat{r})\delta^{\rm b}(\hat{r};\vec{x})$ (10)

where superscript $\rm b$ denotes a binned quantity, and we have defined the
binned density field

$\displaystyle\delta^{\rm b}(\hat{r};\vec{x})\equiv\int
r^{2}dr\;\delta(\vec{x}+\vec{r})\Phi^{\rm b}(r).$ (11)

$\Phi^{\rm b}$ is a binning function demanding that $r$ is in the $\rm b^{\rm
th}$ bin. We define the binning function to be a Heaviside function normalized
by the volume of a given radial bin $V_{\mathrm{b}}=\int dr\;r^{2}\Theta^{\rm
b}(r)$ so that $\Phi^{\rm b}$ can be written as

$\displaystyle\Phi^{\rm b}(r)=\frac{\Theta^{\rm b}(r)}{V_{\rm b}}.$ (12)

With this binning scheme we can then write

$\displaystyle a^{\rm b}_{\ell m}(\vec{x})=\int d^{3}\vec{r}\;Y_{\ell
m}^{*}(\hat{r})\Phi^{\rm b}(r)\delta(\vec{x}+\vec{r}).$ (13)

The above still clearly has the structure of a convolution and consequently
can be obtained about all points $\vec{x}$ at once with an FFT.333equation 13
technically has the form of a correlation due to the addition of $\vec{x}$ and
$\vec{r}$ in $\delta$. In this work we also label it as a convolution to avoid
confusion with the correlation coefficients. Finally, we can rewrite equation
9 as

$\displaystyle\hat{\zeta}_{\ell}^{\;{\rm b}_{1}{\rm
b}_{2}}(\vec{x})=\frac{(-1)^{\ell}}{\sqrt{2\ell+1}}\;\delta(\vec{x})\sum_{m=-\ell}^{\ell}a_{\ell
m}^{\rm b_{1}}(\vec{x})a_{\ell m}^{\rm b_{2}\;*}(\vec{x}),$ (14)

where ${\rm b}_{1}$ and ${\rm b}_{2}$ represent the indices of the radial bins
$r_{1}$ and $r_{2}$ fall into. The multipole coefficients $\zeta_{\ell}$ of
the full 3PCF are then constructed as the volume average of equation 14.

###  2.2 Projected 3PCF

We now turn our attention to the projected 3PCF. Here the density is projected
onto a 2D plane and therefore becomes a function of a 2D vector $\vec{x}$,
which we then average over all 2D rotations. To do this we start with the full
(projected) 3PCF As in Section 2.1, the triangle is parameterized by side
lengths $\rho_{i}$ and the cosine of their enclosed angle:

$\displaystyle\hat{\rho}_{1}\cdot\hat{\rho}_{2}\equiv\cos\phi_{12}=\cos(\phi_{2}-\phi_{1})$
(15)

which we can rewrite as

$\displaystyle\cos(\phi_{2}-\phi_{1})=\frac{1}{2}\left[e^{i(\phi_{2}-\phi_{1})}+e^{-i(\phi_{2}-\phi_{1})}\right].$
(16)

Since the enclosed angle between the side lengths $\rho_{1}$ and $\rho_{2}$
can be rewritten in terms of complex exponentials, it is well motivated for us
to expand the angle-dependence of the projected 3PCF in the basis of Fourier
modes $\exp(im\phi_{12}$) yielding an analogous expression to equation 1 in
2D:

$\displaystyle\zeta_{\rm
proj}(\rho_{1},\rho_{2},\hat{\rho}_{1}\cdot\hat{\rho}_{2})=\sum_{m=0}^{\infty}\zeta_{m}(\rho_{1},\rho_{2})e^{im\phi_{12}}.$
(17)

We can rewrite the basis functions as

$\displaystyle e^{im\phi_{12}}=e^{im\phi_{2}}e^{-im\phi_{1}}.$ (18)

From here, the local projected 3PCF can be written as

$\displaystyle\hat{\zeta}_{\rm
proj}(\rho_{1},\rho_{2};\hat{\rho}_{1}\cdot\hat{\rho}_{2};\vec{x})=\left<\delta(\vec{x})\delta(\vec{x}+\vec{\rho}_{1})\delta(\vec{x}+\vec{\rho}_{2})\right>_{\mathcal{R}},$
(19)

where the angle brackets with subscript $\mathcal{R}$ indicate an average over
joint 2D rotations $\mathcal{R}$ of $\hat{\rho}_{1}$ and $\hat{\rho}_{2}$
about $\vec{x}$ and $\delta$ is the projected density field. The projected
3PCF coefficients can be constructed by an area average of the "local"
estimates of the projected 3PCF. In 2D we have

$\displaystyle\zeta_{\rm
proj}(\rho_{1},\rho_{2})=\int\frac{d^{2}\vec{x}}{A}\;\hat{\zeta}_{\rm
proj}(\rho_{1},\rho_{2};\vec{x}),$ (20)

and since projecting onto the Fourier basis is a linear operation, we have
that

$\displaystyle\zeta_{m}(\rho_{1},\rho_{2})=\int\frac{d^{2}\vec{x}}{A}\;\hat{\zeta}_{m}(\rho_{1},\rho_{2};\vec{x}).$
(21)

The radial coefficients of the projected 3PCF can be computed as an area
average of the local estimates. The local estimate of the multipole
coefficients can be written as

$\hat{\zeta}_{m}(\rho_{1},\rho_{2};\vec{x})=\int
d\phi_{1}d\phi_{2}\;e^{-im\phi_{2}}e^{im\phi_{1}}\\\
\times\delta(\vec{x})\delta(\vec{x}+\vec{\rho}_{1})\delta(\vec{x}+\vec{\rho}_{2}).$
(22)

Using orthogonality of the Fourier modes to project onto $\exp(im\phi)$, we
must integrate against its conjugate. We define

$\displaystyle c_{m}(\rho_{i};\vec{x})\equiv\int
d\phi\;e^{-im\phi}\delta(\vec{x}+\vec{\rho}_{i}),$ (23)

which has the structure of a convolution. We exploit the fact that
$c_{-m}=c_{m}^{*}$ as can be seen from making this replacement in equation 23.
These convolution coefficients can be directly measured with sarabande. The
$c_{m}(\rho_{i};\vec{x})$ are the projected equivalent of the $a_{\ell
m}(r_{i},\vec{x})$. From this, we can write $\hat{\zeta}_{m}$ as their
product:

$\hat{\zeta}_{m}(\rho_{1},\rho_{2};\vec{x})=\frac{\delta(\vec{x})}{(2\pi)^{2}}\;c_{m}(\rho_{1};\vec{x})c_{m}^{*}(\rho_{2};\vec{x}).$
(24)

Now we may proceed analogously to our work in equation 10 to bin these
coefficients.

$c_{m}(\rho_{i};\vec{x})\to c_{m}^{\rm b}(\vec{x})=\int\delta^{\rm
b}(\hat{\rho}_{i};\vec{x})e^{-im\phi}d\phi,$ (25)

where the binned density field is

$\delta^{\rm
b}(\hat{\rho}_{i};\vec{x})\equiv\int\rho_{i}d\rho_{i}\;\delta(\vec{x}+\vec{\rho}_{i})\Phi^{\rm
b}(\rho_{i}),$ (26)

and $\Phi^{\rm b}$ is a binning function. We define this binning function to
be a Heaviside function normalized by the area of a given radial bin
$A_{\mathrm{b}}=\int dr\;\Theta^{\rm b}(r)$ so that $\Phi^{\rm b}$ can be
written as

$\displaystyle\Phi^{\rm b}(r)=\frac{\Theta^{\rm b}(r)}{A_{\rm b}}.$ (27)

This binning function demands that $\rho_{i}$ is in a normalized radial bin
having integral unity. With the binning function we can write the multipole
coefficients of the projected 3PCF as

$\displaystyle\hat{\zeta}_{m}^{\mathrm{b_{1}}\mathrm{b_{2}}}(\vec{x})=\frac{\delta(\vec{x})}{(2\pi)^{2}}\;c_{m}^{\rm
b_{1}}(\vec{x})c_{m}^{\rm b_{2}\;*}(\vec{x}),$ (28)

where $\rm b_{1}$ and $\rm b_{2}$ represent the indices of the radial bins
$\rho_{1}$ and $\rho_{2}$ fall into. If we sum over all values of $m$ we
construct the projected 3PCF

$\displaystyle\zeta_{\rm
proj}^{\mathrm{b_{1}}\mathrm{b_{2}}}(\hat{\rho}_{1}\cdot\hat{\rho}_{2};\vec{x})=\frac{\delta(\vec{x})}{(2\pi)^{2}}\sum_{m=0}^{\infty}\;c_{m}^{\rm
b_{1}}(\vec{x})c_{m}^{\rm b_{2}\;*}(\vec{x})e^{im\phi_{12}}.$ (29)

These $c_{m}^{\rm b}$ coefficients can be interpreted as the convolution of
the density field with the kernel

$K_{m}^{\rm b}(\vec{\rho}_{i})\equiv\Phi^{\rm b}(\rho_{i})e^{-im\phi}.$ (30)

We can thus write the $c_{m}^{\rm b}$ coefficients as

$c_{m}^{\rm b}(\vec{x})=\left[\delta(\vec{\rho}_{i})\star K_{m}^{\rm
b}(\vec{\rho}_{i})\right](\vec{x}).$ (31)

We define the $\star$ symbol to denote a convolution.

###  2.3 Full 4PCF

Here we briefly discuss the derivation of the 4PCF. Further details are given
on the derivation in Philcox et al. (2022) and further information on the
isotropic basis functions can be found in Cahn & Slepian (2020). The formalism
of the full 4PCF is similar to that of the full 3PCF; now we just have four
points in the estimator as opposed to three. This translates to a tetrahedral
geometry instead of a triangular one. Instead of just two vectors: we now have
three $\vec{r}_{1}$, $\vec{r}_{2}$, and $\vec{r}_{3}$. We can write the full
4PCF in its most basic form as

$\zeta=\langle\delta(\vec{x})\delta(\vec{x}+\vec{r}_{1})\delta(\vec{x}+\vec{r}_{2})\delta(\vec{x}+\vec{r}_{3})\rangle_{\mathcal{R}}.$
(32)

We have left unstated the arguments on the left-hand side intentionally as
there are different ways to parameterize the full 4PCF, as we will shortly
discuss. We can characterize the tetrahedron by the lengths $r_{1}$, $r_{2}$,
and $r_{3}$ and the angles between $\vec{r}_{1}$, $\vec{r}_{2}$, and
$\vec{r}_{3}$. The dependence on these angles can be decomposed in the basis
of isotropic functions developed by Cahn & Slepian (2020). We then have the
full 4PCF as a sum of radial coefficients $\zeta_{\Lambda}$ times these
isotropic basis functions $\mathcal{P}_{\Lambda}$ for the angular dependence:

$\zeta=\sum_{\Lambda}\hat{\zeta}_{\Lambda}^{R}\mathcal{P}_{\Lambda}(\hat{R}).$
(33)

We define $\Lambda\equiv\\{\ell_{1},\ell_{2},\ell_{3}\\}$ to be the set of
angular momenta describing the angles of the quadruplets used for the full
4PCF. The isotropic orthonormal basis used in equation 33 is defined as

$\mathcal{P}_{\Lambda}(\hat{R})=\sum_{M}C^{\Lambda}_{M}Y_{\ell_{1}m_{1}}(\hat{r}_{1})Y_{\ell_{2}m_{2}}(\hat{r}_{2})Y_{\ell_{3}m_{3}}(\hat{r}_{3}),$
(34)

and $C^{\Lambda}_{M}$ is given by the Wigner 3-$j$ symbol with an additional
phase:

$C^{\Lambda}_{M}=(-1)^{\ell_{1}+\ell_{2}+\ell_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}.$ (35)

For compactness we have defined $M\equiv\\{m_{1},m_{2}.m_{3}\\}$,
$\hat{R}\equiv\\{\hat{r}_{1},\hat{r}_{2},\hat{r}_{3}\\}$, and
$R\equiv\\{r_{1},r_{2},r_{3}\\}$. Equation 33 shows the radial coefficients
$\hat{\zeta}_{\Lambda}^{R}$ must be computed to compute the full 4PCF. This
may be done in terms of binned convolution coefficients $a_{\ell m}^{\rm
b}(\vec{x})$ defined in equation 13, leading to the radial coefficients about
each point $\vec{x}$ as

$\hat{\zeta}_{\ell_{1}\ell_{2}\ell_{3}}^{\;\rm b_{1}\rm b_{2}\rm
b_{3}}(\vec{x})=\delta(\vec{x})\sum_{M}\mathcal{E}(\Lambda)\;C^{\Lambda}_{M}\\\
\times\;a_{\ell_{1}m_{1}}^{\rm b_{1}}(\vec{x})\;a_{\ell_{2}m_{2}}^{\rm
b_{2}}(\vec{x})\;a_{\ell_{3}m_{3}}^{\rm b_{3}}(\vec{x}),$ (36)

where $\mathcal{E}(\Lambda)$ is +1 for even parity and -1 for odd parity. Even
parity occurs when $\ell_{1}+\ell_{2}+\ell_{3}$ is an even number, and odd
parity occurs when the sum of the angular momenta is odd. For further details
on this derivation, see Section 3.2 of Philcox et al. (2022).

###  2.4 Projected 4PCF

For the projected 4PCF, we take a slightly different approach from that of the
previous subsections. We begin by setting up a coordinate system whose origin
is defined by one particle, which is always allowed by translational
invariance. We then expand the dependence on $\vec{\rho}_{i}$ for the three
remaining points in polar coordinates:

$\displaystyle\delta(\vec{x})\delta(\vec{x}+\vec{\rho}_{1})\delta(\vec{x}+\vec{\rho}_{2})\delta(\vec{x}+\vec{\rho}_{3})=$
(37)
$\displaystyle\delta(\vec{x})\sum_{M}c_{m_{1}}(\rho_{1};\vec{x})e^{im_{1}\phi_{1}}c_{m_{2}}(\rho_{2};\vec{x})e^{im_{2}\phi_{2}}c_{m_{3}}(\rho_{3};\vec{x})e^{im_{3}\phi_{3}}.$

This can be done because the Fourier modes are a basis for any function of
$\phi$. We capture the radial dependence on $\rho_{i}$ with the coefficients
$c_{m_{i}}$. These coefficients were previously defined in Section 2.2 (see
equation 23). We let $M\equiv\\{m_{1},m_{2},m_{3}\\}$ for compactness as done
in Section 2.2.

We now consider the constraints required for the projected 4PCF to be
rotation-invariant. The "local" estimates of the projected 4PCF are found by
averaging over joint rotations as done in the previous sections:

$\hat{\zeta}_{\mathrm{proj}}(\rho_{1},\rho_{2},\rho_{3};\hat{\rho}_{1},\hat{\rho}_{2},\hat{\rho}_{3};\vec{x})=\\\
\langle\delta(\vec{x})\delta(\vec{x}+\vec{\rho}_{1})\delta(\vec{x}+\vec{\rho}_{2})\delta(\vec{x}+\vec{\rho}_{3})\rangle_{\mathcal{R}}.$
(38)

Returning to equation 37, we can perform the average over rotations by adding
an arbitrary displacement in angle, $\phi$, to each $\phi_{i}$ and integrating
over it. We find

$\hat{\zeta}_{\mathrm{proj}}(\rho_{1},\rho_{2},\rho_{3};\hat{\rho}_{1},\hat{\rho}_{2},\hat{\rho}_{3};\vec{x})=\\\
\qquad\delta(\vec{x})\sum_{M}c_{m_{1}}(\rho_{1};\vec{x})c_{m_{2}}(\rho_{2};\vec{x})c_{m_{3}}(\rho_{3};\vec{x})\;\\\
\times\frac{1}{2\pi}\int
d\phi\;e^{im_{1}(\phi_{1}+\phi)}e^{im_{2}(\phi_{2}+\phi)}e^{im_{3}(\phi_{3}+\phi)}.$
(39)

From here we may rewrite the angular integral by factoring the exponentials

$\frac{1}{2\pi}\int
d\phi\;e^{im_{1}(\phi_{1}+\phi)}e^{im_{2}(\phi_{2}+\phi)}e^{im_{3}(\phi_{3}+\phi)}=\\\
\left(e^{im_{1}\phi_{1}}e^{im_{2}\phi_{2}}e^{im_{3}\phi_{3}}\right)\;\frac{1}{2\pi}\int
d\phi\;e^{i(m_{1}+m_{2}+m_{3})\phi}.$ (40)

The integral is only non-zero if $m_{1}+m_{2}+m_{3}=0$. Thus we may replace
$m_{3}=-(m_{1}+m_{2})$. We then have the local estimate of the coefficient of
the rotation-averaged projected 4PCF as

$\hat{\zeta}_{\mathrm{proj}}(\rho_{1},\rho_{2},\rho_{3};\vec{x})=\\\
\qquad\delta(\vec{x})\sum_{m_{1},\;m_{2}}c_{m_{1}}(\rho_{1};\vec{x})c_{m_{2}}(\rho_{2};\vec{x})c_{-(m_{1}+m_{2})}(\rho_{3};\vec{x})\\\
\times\left(e^{im_{1}\phi_{1}}e^{im_{2}\phi_{2}}e^{-i(m_{1}+m_{2})\phi_{3}}\right).$
(41)

From here we reuse equation 25 to bin the convolution coefficients
$c_{m_{i}}$. For compactness we define $\phi_{13}\equiv\phi_{1}-\phi_{3}$ and
$\phi_{23}\equiv\phi_{2}-\phi_{3}$. Finally, our local coefficient estimate of
the projected 4PCF becomes

$\displaystyle\hat{\zeta}_{\mathrm{proj}}^{\mathrm{b_{1}}\mathrm{b_{2}}\mathrm{b_{3}}}(\hat{\rho}_{1},\hat{\rho}_{2},\hat{\rho}_{3};\vec{x})=$
(42) $\displaystyle\qquad\delta(\vec{x})\sum_{m_{1},\;m_{2}}c_{m_{1}}^{\rm
b_{1}}(\vec{x})c_{m_{2}}^{\rm b_{2}}(\vec{x})c_{-(m_{1}+m_{2})}^{\rm
b_{3}}(\vec{x})\left(e^{im_{1}\phi_{13}}e^{im_{2}\phi_{23}}\right).$

The projected 4PCF coefficients are then

$\hat{\zeta}_{m_{1}m_{2}}^{\;\mathrm{b_{1}}\mathrm{b_{2}}\mathrm{b_{3}}}(\vec{x})=\delta(\vec{x})c_{m_{1}}^{\rm
b_{1}}(\vec{x})c_{m_{2}}^{\rm b_{2}}(\vec{x})c_{-(m_{1}+m_{2})}^{\rm
b_{3}}(\vec{x}).$ (43)

###  2.5 Normalization

Here we discuss the normalization of the FFT-based NPCF. We need to normalize
because by definition the NPCF is a measure of the excess probability over
random of finding a particle N-tuplet (e.g. the 2PCF is the probability excess
of finding a galaxy pair over a random distribution). In order to normalize
the measured NPCF coefficients to have this meaning we must follow an
analogous normalization scheme to that of codes such as encore.444encore
(Philcox et al., 2022) uses a generalization of the Landy-Szalay estimator
wherein random particles enable subtracting the mean and dividing by it. This
is instead of working with the density contrast $\delta$ directly. The
particle-based normalization of encore includes the number density to
represent the number of neighbors in a bin around a primary galaxy:
$\mathcal{N}=N_{\rm gal}\prod_{i=1}^{N-1}(\bar{n}V_{\rm bin,i})$, with $V_{\rm
bin}=(4\pi/3)(r_{+}^{3}-r_{-}^{3})$ with $r_{+}$ and $r_{-}$ respectively the
upper and lower bounds of a given spherical shell. Typical particle-based NPCF
codes measure the correlations on a set of $N_{\rm gal}$ galaxies within a
cosmological volume using an estimator. Our grid-based code does not have a
distribution of particles; instead, sarabande takes a density contrast field
$\delta$ as its input. We define $\delta=n/\bar{n}-1$ where $n$ is the number
density and $\bar{n}$ is the average number density.

####  2.5.1 Normalizing the Spherical-Shell Bins

In sarabande, the full 3/4 PCFs are measured on 3D bins in the side lengths of
the configuration about the primary. Each bin is normalized to have volume
unity as stated in equation 12. This normalization must be adjusted to reflect
the fact that the underlying space is on a grid. Here we show how to perform
this normalization. Since sarabande is computed using a grid mesh, we count
the number of cells within a radial bin and multiply that by the volume of
each cell. The volume of each cell in 3D is $V_{\rm cell}=(L/N_{\rm
mesh})^{3}$ where $L$ is the physical box size of our data. Therefore the bin
volume is defined as

$V_{\mathrm{bin}}=V_{\rm cell}\times N_{\textrm{cells in
bin}}=\left(\frac{L}{N_{\rm mesh}}\right)^{3}N_{\textrm{cells in bin}}.$ (44)

The text below equation 11 has additional information regarding our binning
scheme. We will divide the full 3/4 PCF coefficients by the product of bin
volumes to normalize as discussed below.

####  2.5.2 Total Number of Objects

Particle-based codes such as encore typically normalize by the total number of
galaxies $N_{\rm gal}$ in a cosmological volume. For particle-based codes this
is simply the number of objects (or their weighted sum, if weights are used
such as might be done to correct for survey systematics or survey geometry),
but for density fields the number of particles used is not as straightforward.
If we assume the input density field to have a mean of zero, then the
equivalent value is

$N_{\rm gal}=\sum_{i}\sum_{j}\sum_{k}[\delta(\vec{x})+1].$ (45)

In essence, we sum the density field over all the cells in the mesh, which are
indexed by ($i,j,k$), to get the equivalent of the number of particles.

To properly normalize the full 3/4 PCF coefficients in equation 14 and
equation 32 we divide the coefficients by a normalization coefficient
$\mathcal{N}$ which normalizes our bins to have volume unity and allows us to
compare our coefficients to particle-based codes such as encore as done in
Figure 2. The final normalization coefficient $\mathcal{N}$ reads

$\mathcal{N}=N_{\rm gal}\prod_{n=1}^{N-1}\bar{n}\;V_{\mathrm{bin},n}.$ (46)

We define $N$ is the order of the correlation function (e.g. $N=3$ for the
3PCF). We include the average number density $\bar{n}\equiv N_{\rm gal}/L^{3}$
in the product to represent the number of neighbors in a bin around primary
point.

####  2.5.3 Projected Case

If we follow the same normalization procedure as the full 3/4 PCF for the
projected 3/4 PCF, we find that the main difference is that instead of volumes
in 3D we are working with areas in 2D. In sarabande, the projected 3/4 PCFs
are measured on 2D bins in the side lengths of the configuration about the
primary. Each bin is normalized to have area unity as stated in equation 27.
This normalization must be adjusted to reflect the fact that the underlying
space is on a grid. Here we show how to perform this normalization. We still
have a grid mesh in the projected case so to get the total area of a bin we
have multiply the area of a cell by the number of cells within a bin. The area
of a cell in 2D is $A_{\rm cell}=(L/N_{\rm mesh})^{2}$. Therefore the bin area
is defined as

$A_{\mathrm{bin}}=A_{\rm cell}\times N_{\textrm{cells in
bin}}=\left(\frac{L}{N_{\rm mesh}}\right)^{2}N_{\textrm{cells in bin}}.$ (47)

The text below equation 26 has additional information regarding our binning
scheme in 2D.

We follow the same procedure for creating an equivalent expression for the
number of particles as in equation 45. In 2D we find

$N_{\rm gal}=\sum_{i}\sum_{j}[\delta(\vec{x})+1],$ (48)

where we sum the density field over all the cells in the mesh indexed by
($i,j$). If we assume the density field has a mean of zero.

To properly normalize the projected 3/4 PCF coefficients in equation 29 and
equation 43 we divide the coefficients by a normalization coefficient
$\mathcal{N_{\rm proj}}$ which normalizes our bins to have area unity. This
allows us to compare our coefficients to particle-based codes as done in
Figure 1. In Section 3.3, we describe the code including an optional flag to
turn on this normalization procedure. The final normalization coefficient
$\mathcal{N_{\rm proj}}$ reads

$\mathcal{N_{\rm proj}}=N_{\rm
gal}\prod_{n=1}^{N-1}\bar{n}\;A_{\mathrm{bin},n}.$ (49)

We define $N$ to be the order of the correlation function as done at the end
of Section 2.5.2. We include the average number density $\bar{n}\equiv N_{\rm
gal}/L^{2}$ in the product to represent the number of neighbors in a bin
around primary point.

Without normalization sarabande agrees with the particle-based direct counting
code to a precision of $\sim 10^{-13}$. When comparing sarabande to encore
without normalization we found the full 3/4 PCFs agreed with a precision of
$\sim 10^{-6}$ as seen in Figure 2. After our normalization scheme, the
agreement precision of sarabande with particle-based codes increases with
resolution. For a grid resolution of $N_{\rm g}=256^{2}$ sarabande agrees with
the particle-based direct counting code to a precision of $\sim 10^{-3}$ as
seen in Figure 1. As the resolution of the grid is increased the precision
increases as one would expect because the width of each cell decreases
allowing for a better approximation of the analytical bin volume/area
(equation 44 and equation 47). The same normalization behavior applies to both
the full and projected 3/4 PCF in sarabande. The particle-based direct
counting code used for testing is publicly available on the sarabande github
described in Section 3.3.

Figure 1: This figure shows the agreement between sarabande and a particle-
based direct counting (DC) code before and after normalization discussed in
Section 2.5 (equations 47-49). The left column displays the projected 3PCF
with $m=2$ as an example of how normalization impacts the agreement between
sarabande and particle-based codes. The right column shows the projected 4PCF
with $m_{1}=2$ and $m_{2}=2$ as an example. For all plots of projected 3/4 PCF
coefficients $\hat{\zeta}$ we add an arbitrary offset to the real parts and
subtract the offset from the imaginary parts for clarity. For all residual
plots the red dashed curve shows the residual of the real part of the
projected 3/4 PCF coefficients and the blue solid curve shows the residual of
the imaginary part. Before normalization sarabande has an agreement precision
of $\sim 10^{-13}$ with the direct-counting code, after normalization
sarabande has an agreement precision of $\sim 10^{-3}$ for this choice of
resolution. This precision after normalization improves as we increase the
resolution. The bottom row shows how we map the 2 or 3 radial bins to a 1D bin
index. For 2 radial bins we have a nested for loop where $\textrm{b}_{1}$
iterates from 0 to $N_{\rm bins}$ then $\textrm{b}_{2}$ iterates from
$\textrm{b}_{1}+1$ to $N_{\rm bins}$ (we avoid bin overlap where
$\textrm{b}_{1}=\textrm{b}_{2}$). Each iteration step is assigned to the
combination ($\textrm{b}_{1}$, $\textrm{b}_{2}$). For 3 radial bins we add an
additional nested for loop for $\textrm{b}_{3}$. With 3 radial bins we iterate
through the bin combinations such that $\textrm{b}_{1}$ iterates from 0 to
$N_{\rm bins}$ and $\textrm{b}_{1}+1<\textrm{b}_{2}+1<\textrm{b}_{3}<N_{\rm
bins}$ to avoid bin overlap. Each combination of ($\textrm{b}_{1}$,
$\textrm{b}_{2}$, $\textrm{b}_{3}$) is assigned to an iteration index. We
iterate with these conditions so that we only display the unique correlation
coefficients. We plot the coefficients as a function of 1D bin index to show
all the 3/4 PCF coefficients measured. The input data used for this figure was
a set of 41 randomly placed particles on a grid with a resolution of $N_{\rm
g}=256^{2}$ grid cells. To create the particles on the grid we assign cells
within the grid a weight equal to unity. We use a small sample of particles
because of the slow nature of the direct counting method. For this measurement
we use 10 linearly-spaced radial bins. Figure 2: This figure shows the
accuracy of sarabande relative to a naive direct counting code and encore. The
input data used was created by randomly placing 104 particles on a grid of
resolution $N_{\rm g}=128^{3}$. To create the particles on the grid we assign
cells within the grid a weight equal to unity. We compare non-normalized
coefficients for both the full 3PCF and the full 4PCF to show that the
measured coefficients agree. For the full 3/4 PCF we achieve a precision $\sim
10^{-6}$ without normalization as seen in the bottom of each plot. In both
plots we show the residual between sarabande and encore in red and the
residual between the direct counting code and encore in blue. The precision of
agreement between sarabande and the direct counting code is $\sim 10^{-14}$.
We use 10 linearly-spaced radial bins for all of the measurements made in this
figure. We use the same mapping of bins to a 1D index as discussed in the
bottom row of Figure 1.

##  3 Code structure and Usage

sarabande offers two pathways to calculate correlation functions: one pathway
is used to calculate full correlation functions, while the other is used to
calculate projected correlation functions. Both pathways resemble each other
and are outlined in Figure 3. All computations start with creating a grid and
radial bins for the input data (2D for projected and 3D for full). From there,
both pathways generate kernels based on their basis functions and radial bins
as seen in Figure 4. Once the basis kernels are created, then sarabande
convolves them with the data to yield the convolution coefficients ($a_{\ell
m}^{\mathrm{b}}(\vec{x}$) for 3D and $c_{m}^{\mathrm{b}}(\vec{x})$ for 2D)
necessary for computing the desired correlation function coefficients
$\hat{\zeta}$. The final step is to assemble the convolution coefficients into
products and sum those products over all cells. We next describe each pathway
in more detail in Section 3.1 and Section 3.2; we then give a guide to code
use in Section 3.3 and suggest how to explore the code to gain familiarity in
Section 3.4 using an online google colab notebook we have created.

Figure 3: Here we show how sarabande proceeds through the calculation of the
3/4 PCFs. The flowchart reads left to right starting with importing the
sarabande package. The code splits into two pathways depending on if the user
chooses to calculate the projected or full 3/4 PCF. The two pathways then
diverge again when deciding to compute the 3 or the 4 PCF (Algorithms 5-8).
The code follows a general structure such that we first create kernels out of
radially-binned basis functions ($Y_{\ell m}$ or $e^{im\phi}$) which we then
convolve with the density field to calculate our convolution coefficients
($a_{\ell m}^{\mathrm{b}}(\vec{x})$ or $c_{m}^{\mathrm{b}}(\vec{x})$). These
convolution coefficients are then combined in the final step of each pathway
to give the desired full or projected 3/4 PCF radial coefficients. Figure 4:
Visual example of the binning of spherical harmonics to create kernels for the
full 3/4 PCFs. In this figure we have a slice of $Y_{10}$ being binned into 4
separate radial bins numbered 0-3. We convolve these kernels with the original
data to create the convolution coefficients.

###  3.1 Full NPCF Pathway

The full 3/4 PCF calculation differs from the projected mostly in the handling
of the spherical harmonic kernel calculation. 3D gridded data has a much
larger memory footprint than 2D gridded data; thus to hold the computation of
$\zeta$ in memory we must offload intermediate calculations to disk and delete
them when they are no longer needed. The full pathway is divided into several
parts. The pathway starts with the generation of spherical harmonic kernels by
radially binning the spherical harmonics as described in Algorithm 1 and
saving their Fourier Transforms. The second part of the pathway convolves a
given dataset with the kernels to calculate the convolution coefficients
$a_{\ell m}^{\mathrm{b}}(\vec{x})$ as outlined in Algorithm 3. These
convolution coefficients are then assembled to yield the full 3/4 PCF as
desired. The only difference between calculating the full 3PCF and the full
4PCF is the final algorithm chosen for assembling convolution coefficients.
For the full 3PCF we use algorithm 5 which is adapted from Portillo et al.
(2018) and for the full 4PCF we use algorithm 6 adapted from Philcox et al.
(2022).

During the entire computation of the full 3PCF we only keep at most three
arrays, each of size $N_{\mathrm{g}}^{3}$ in memory at a time, where
$N_{\mathrm{g}}$ is the number of grid cells on each side of the data grid.
This approach requires saving intermediate calculations to disk one array at a
time. Only when we need an array for a subsequent calculation is it read back
in from disk. When computing the convolution coefficients via Algorithm 3, the
three arrays needed for the calculation are the Fourier Transform of the data,
the Fourier Transform for a single kernel, and an array to hold the product.
Then for the final combination of convolution coefficients in Algorithm 5 the
three needed arrays are the density field and two different convolution
coefficients.

For the full 4PCF we only keep up to four arrays of size $N_{\mathrm{g}}^{3}$
in memory at a time. The number of arrays in memory is the same as for the
3PCF except for the final step of combining convolution coefficients. The four
arrays in memory during this last step, given by Algorithm 6, are the density
field and three different convolution coefficients. Overall, this limited-
memory footprint approach allows sarabande to measure the 3/4 PCFs of high-
resolution gridded datasets and simulations.

###  3.2 Projected NPCF Pathway

For the projected 3/4 PCF calculation we can afford to keep intermediate
calculations in memory due to the much smaller memory footprint of 2D data.
This difference inherently means the projected calculations will be much
faster than the full 3/4 PCF calculations due to no time being spent on File
I/O operations. The projected pathway of sarabande follows the same structure
as the full pathway just in 2D instead of 3D.

The projected pathway also starts with the generation of kernels by radially
binning basis functions as described in Algorithm 2. These kernels use a
different set of basis functions; we call them the Fourier basis.
Specifically, they are complex exponentials, which are essentially a
projection of the spherical harmonics onto a 2D plane.555We recall that
$Y_{\ell m}\propto P_{\ell}^{m}(\cos\theta)\exp[im\phi]$, with $P_{\ell}^{m}$
an associated Legendre polynomial. Setting $\cos\theta=0$ as is the case in
the $xy$-plane renders this purely dependent on $\exp[im\phi]$. Sections 2.2
and 2.4 describe the use of this basis in more detail. Once the kernels have
been created the projected pathway convolves them with the density field to
calculate the convolution coefficients $c_{m}^{\mathrm{b}}(\vec{x})$ as
outlined in Algorithm 4. With these 2D convolution coefficients, we assemble
them into the projected 3/4 PCF as desired. We use Algorithm 7 for the
projected 3PCF (developed in Section 2.1 of Saydjari et al. 2021) and
Algorithm 8 for the projected 4PCF (which is entirely new).

###  3.3 Guide to Code Use

sarabande was designed to be simple to use while also allowing the user to
execute each step as desired. There are various functions listed in the
documentation on github: sarabande github. 666sarabande github:
https://github.com/James11222/sarabande sarabande makes use of object-oriented
programming in python to create a measure object. This object can be passed
into the zeta calculation function to calculate the desired correlation
function as follows:

⬇

import sarabande

NPCF_obj = sarabande.measure(**kwargs)

sarabande.calc_zeta(NPCF_obj)

zeta = NPCF_obj.zeta

where **kwargs are the possible arguments in the constructor function. There
are many arguments and comprehensive definitions for usage that can be found
in the documentation at the github repository linked above. The two most
important parameters are projected [boolean] and NPCF [integer] as they govern
which algorithms will be used to calculate $\zeta$. The variable type is
denoted by square brackets. We provide an optional normalized [boolean] flag
argument to activate the normalization scheme discussed in Section 2.5. This
calculation is then stored by adding the .zeta attribute to the measure
object.

###  3.4 Explore the Code

We have also provided tutorial colab notebooks to show the process of
calculating the 3/4 PCFs (projected and full). These notebooks walk through
each step of the code interactively without requiring one to examine the
source code of sarabande. The notebooks are located here: sarabande drive.
777https://drive.google.com/drive/folders/1oEum7DThj9kkAbyacx0eX-x9qY-8cZv0?usp=sharing

##  4 Performance and Validation

The attractive features of sarabande are that it is written in python, and
that it formally scales with resolution as $\mathcal{O}(N_{\rm g}\log N_{\rm
g})$, as shown in Figure 5. This scaling is made possible by the use of the
FFT algorithm originally discovered by Cooley & Tukey (1965) which has a
complexity scaling of $\mathcal{O}(N_{\rm g}\log N_{\rm g})$. The history of
the FFT algorithm is discussed in Heideman et al. (1984). The scaling of
sarabande enables computing higher-order correlation functions of high-
resolution gridded datasets and simulations. Using a particle-based code such
as encore (Philcox et al., 2022) the formal scaling would be
$\mathcal{O}(N_{\mathrm{p}}^{2})$, with $N_{\mathrm{p}}$ the number of
particles in the dataset.888In actual practice for test datasets, encore
scales linearly in number of particles for the full 4PCF, because it is not
the $a_{\ell m}$ computation that actually dominates the runtime, but
assembling the products of coefficients. encore, and the original particle-
based 3PCF algorithm Slepian & Eisenstein (2015a) on which it is based, does
scale as $N_{\mathrm{p}}^{2}$ for 3PCF. encore does not include projected
statistics. Previously the higher-order statistics provided by the N-point
correlation function were primarily applied in a cosmological setting in the
context of galaxy redshift surveys (Dawson et al., 2013; Eisenstein et al.,
2011). There has been some work in the application of the 3PCF and Bispectrum
on Magnetohydrodynamic (MHD) simulations of the ISM (Portillo et al., 2018;
O’Brien et al., 2022). sarabande allows us to measure higher-order statistics
of any gridded data set, opening the door for N-point statistical analysis in
other subfields of astronomy than just cosmology. As an example of this, we
show sample measurements made on a simulation of the ISM produced by the CATS
collaboration (Bialy & Burkhart, 2020; Burkhart et al., 2009; Cho & Lazarian,
2003; Portillo et al., 2018; Burkhart et al., 2020).

sarabande can take either 2D data slices or 3D data cubes as input to measure
either projected or full 3/4 PCFs. Figure 7 serves as a typical example output
one can expect when using sarabande in practice. The top half of the figure
shows example coefficients of a full 3/4 PCF and the bottom half shows the
process of measuring a projected 3/4 PCF. The density field used in Figure 7
is from a simulation of the ISM produced using a third-order-accurate hybrid
Essentially Non-Oscillatory (ENO) scheme developed in Cho & Lazarian (2003) to
solve the ideal magnetohydrodynamic equations in a periodic box with driven
turbulence. The simulation chosen for Figure 7 is characterized by an average
sonic Mach number of $\mathcal{M}_{s}$$\sim$$1.2$ and an average Alfvénic Mach
number of $\mathcal{M}_{A}$$\sim$$2.0$.

The sonic Mach number is defined as $\mathcal{M}_{s}\equiv|\vec{v}|/c_{s}$ and
the Alfvénic Mach number to be
$\mathcal{M}_{A}\equiv|\vec{v}|/\left<v_{A}\right>$. $\vec{v}$ is the velocity
field, $c_{s}$ the isothermal sound speed, $v_{A}$ the Alfvén speed, and angle
brackets denote a spatial average over the entire simulation box. The full
simulation chosen has a resolution of $256^{3}$ grid cells but for the
projected 3/4 PCF we only take a single slice of this density cube such that
there are a total of $256^{2}$ grid cells in the density field. This
simulation is one of many that has been described and used in previous work by
the CATS collaboration (Bialy & Burkhart, 2020; Burkhart et al., 2009; Cho &
Lazarian, 2003; Portillo et al., 2018; Burkhart et al., 2020) The CATS
Database webpage gives more details on the simulations. 999CATS Database:
www.mhdturbulence.com

The bottom half of Figure 7 shows example results for measuring the projected
3/4 PCFs using sarabande. The input for the bottom half of this figure is
simply a slice taken from the full density cube used in the top half of the
figure. For all measurements made we subtract out the mean from the data so
that all 3/4 PCF coefficients represent excesses and deficits relative to
random, depicted by red and blue in the coefficients respectively. A
noticeable difference in the outputs between the full and projected pathways
is that the projected 3/4 PCF coefficients are complex-valued while the full
3/4 PCF coefficients are entirely real-valued. Since the projected 3/4 PCFs
are far less expensive to measure than their full counterparts we can probe
significantly higher resolutions in a fraction of the time. It is apparent
that there is structure in the 3/4 PCF coefficients both in the full and
projected versions. This structure suggests further study which we are
pursuing in a subsequent paper (Sunseri et al. in prep).

Figure 5: Here we show the scaling of the creation of the convolution
coefficients in sarabande (Algorithms 1-4) as outlined in Figure 3. This plot
shows the scaling in clock time vs. $N_{\mathrm{g}}$ per side of a given data
set. $N_{\mathrm{g}}$ is defined as the total number of grid cells used for
the FFT, so the full 3/4 PCFs are using
$(N_{\mathrm{g}}\;\mathrm{per\;side})^{3}$ cells and the projected 3/4 PCFs
are using $(N_{\mathrm{g}}\;\mathrm{per\;side})^{2}$ cells. Both fitted lines
scale as $N_{\mathrm{g}}\log(N_{\mathrm{g}})$ which agrees with the predicted
scaling of our algorithms for the full and projected 3/4 PCFs. The full 3/4
PCFs are calculated with 4 radial bins and an $\ell_{\textrm{max}}=1$ and the
projected 3/4 PCFs are calculated with 10 radial bins and $m_{\textrm{max}}=5$
in this example. We omit the computational time spent combining the
convolution coefficients (Algorithms 5-8) because these algorithms are
strongly dependent on the number of File I/O operations instead of the
resolution of the input data. The amount of File I/O operations needed for
combining the convolution coefficients is determined by the choice of $N_{\rm
bins}$ and $\ell_{\rm max}$. The storage of convolution coefficients on disk
is necessary to keep the calculations in memory as discussed in Section 3.1.
Figure 6: Here we show the impact of parallelizing Algorithm 6 for the full
4PCF. The left pie chart shows the distribution of time spent computing the
full 4PCF without any parallelization. The right pie chart shows the same
calculation where we include parallelization of Algorithm 6. A substantial
amount of time is spent on the combination of convolution coefficients due to
the inefficiency of having a single processor reading in convolution
coefficient files from disk within 8 nested for loops. By unraveling the loops
with a lookup table we were able to distribute the loops amongst multiple
processors. With parallelization we see that the FFTs used to create the
convolution coefficients are the rate-limiting factor. The timings displayed
above are made using 5 radial bins and an $\ell_{\textrm{max}}=2$. Figure 7:
Here we show an example application of sarabande. The top half of the figure
shows the results of the full 3/4 PCF algorithm measurements on the 3D density
field shown on the left. The bottom half of the figure shows the results of
the projected 3/4 PCF algorithm measurement on the projected density field
displayed on the left. The full 3/4 PCF coefficients in the top half are
purely real while the projected coefficients in the bottom half are complex.
The projected 4PCF coefficients shown correspond to $\rm b_{3}=0$. For
measurements we subtract off the mean of the input density field so that
positive and negative coefficients correspond to excess and deficit relative
to random respectively. We used 10 radial bins for the full 3/4 PCFs and 20
radial bins for the projected 3/4 PCFs. For all cases we ignore bin overlap so
that the diagonals, where $\mathrm{b}_{1}=\mathrm{b}_{2}$ or
$\mathrm{b}_{1}=\mathrm{b}_{2}=\mathrm{b}_{3}$, have been set to zero.

###  4.1 Performance

sarabande is a valuable tool due to the improved computational speed gained by
using FFTs as opposed to direct counting methods for measuring 3/4 PCFs. In
Figure 5 we show how the creation of convolution coefficients scale in clock
time as a function of resolution of the input data (number of cells
$N_{\mathrm{g}}$). In Figure 5 we fit the timing data to scalings $\propto
N_{\mathrm{g}}\log N_{\mathrm{g}}$; this is the expected scaling that occurs
for both the projected and the full 3/4 PCFs. We do not include the
coefficient combination algorithms (Algorithms 5-8) in this plot because they
are primarily dependent on the number of radial bins chosen and the maximum
order of angular multipole rather than the grid resolution of the density
field.

As seen in Figure 6, it is clear that without multiprocessing the full 4PCF is
dominated by the combination of convolution coefficients rather than their
creation. Because of this, we parallelize this task of combining coefficients
in the 8 nested for loops of Algorithm 6 by unraveling the loops via multi-
process mapping using a lookup table approach. This approach allows for a
drastic improvement in overall time spent computing the 4PCF depending on how
many processors are available. As an example, for a 2.9 GHz Quad-Core Intel
Core i7 Macbook Pro we saw a roughly factor of 4 speed-up for Algorithm 6 and
expect performance to increase with the number of processors available for the
full 4PCF.

In Figure 8 we provide a breakdown of the computational time spent for each of
the four types of measurements for a given set of parameters. It is clear that
the rate limiting step for the projected 3/4 PCFs is not the calculation of
the convolution coefficients, but instead the combination of the convolution
coefficients to create the correlation function coefficients. The convolution
coefficients are calculated quickly because intermediate results are held in
memory instead of on disk, hence the slowest part of the pathway for the
projected 3/4 PCFs is the nested for loops in Algorithms 7 and 8.

As discussed in Sections 3.1 and 3.2, sarabande offloads intermediate
calculations to disk only in the full pathway due to the large memory
footprint of 3D arrays. The projected pathway outlined in Section 3.2 only
requires 2D arrays, which demand a far smaller amount of memory, allowing us
to keep the entire computation in memory without offloading calculations to
disk. The full 3PCF only requires three $N_{\rm g}^{3}$ arrays in memory at
once: the kernel FFT, the data FFT, and their product. If we take $N_{\rm g}$
to be 256 then we would expect each array to occupy $256^{3}\times 2\times 8$
bytes = 268 MB if we use complex double-precision arrays. Three of these
arrays in memory at any given time leads to a maximum memory usage of $3\times
268$ MB = 0.81 GB total. Similarly, for $N_{\rm g}=512$ we expect a required
total memory of 6.44 GB, and $N_{\rm g}=1024$ we would need 51.54 GB of
memory. For the full 4PCF the only difference is that we have at most four
$N_{\rm g}^{3}$ arrays in memory at once instead of three: three convolution
coefficients and the density field. In this case then if we have $N_{\rm
g}=256$ we would expect $4\times 268$ MB = 1.07 GB. For $N_{\rm g}=512$ we
expect a total memory usage of 8.59 GB and for $N_{\rm g}=1024$ we expect a
total memory usage of 68.72 GB.

The memory footprint is much smaller for the projected 3/4 PCF. If we use an
$N_{\rm g}=256$ in this case then we would expect the $N_{\rm g}^{2}$ complex
double-precision arrays to require $256^{2}\times 2\times 8$ bytes = 1.04 MB
per array. Even for a resolution of $N_{\rm g}=1024$, we expect each array to
use 16.78 MB of memory. This small memory footprint allows us to store all
convolution coefficients in memory at any given time without offloading to
disk. The number of these $N_{\rm g}^{2}$ arrays in memory is determined by
the number of bins used and the maximum value of $m$. The total memory
required is equal to $(\textrm{Memory per Array})\times(N_{\rm bins}\times
m_{\rm max}+1)$. For a resolution of $N_{\rm g}=1024$ we expect to use 1.69 GB
of memory total to store all the necessary arrays for $N_{\rm bins}=20$ and
$m_{\rm max}=5$.

###  4.2 Validation

To validate our code we compare sarabande to other codes that measure the 3/4
PCFs. For the full 3/4 PCFs we compare our code to both the current standard
particle-based code encore and a naive direct counting code which simply
counts triplets or quadruplets respectively. These codes should all be
equivalent for the same input data since they are measuring the same
correlation functions. Since encore does not compute projected 3/4 PCFs, we
can only compare our code to the naive direct counting code in this case.

To generate an equivalent input for both the particle-based codes and our grid
based code we start with an empty grid and randomly populate cells of this
grid with a weight of unity. We then treat the center of each of these
populated cells as the coordinate of the "particles" within the box. The box
is defined to have side length equal to the resolution of the grid (i.e a grid
of resolution $128^{3}$ will correspond to a box of side length $128$).
sarabande takes the original grid as its input while encore and the direct
counting code both take in a list of coordinates for the particles. These
inputs are inherently not precisely the same so we do expect minor
discrepancies that increase in magnitude with lower resolution of the mesh. As
the resolution increases our cells are better approximated as point particles.
We give examples of how the codes compare in Figure 1 and Figure 2. It is
evident from these figures that the codes are measuring the same correlation
coefficients except with minor disagreement. It is expected that with a higher
resolution of the grid these discrepancies will approach machine error
precision.

##  5 Future Applications

There have been previous efforts to use both the full and projected 3PCF to
better understand grid-based simulations of the turbulent ISM (Portillo et
al., 2018; Saydjari et al., 2021). The present work serves as a release of
open software available for making measurements of both full and projected 3/4
PCFs. In a future study (Sunseri et al. 2022 in prep.) we plan to report the
first application of the full and projected 4PCF on the ISM and other possible
applications using sarabande.

In addition to applications at galactic scales, sarabande can also be applied
to datasets on cosmological scales. As discussed elsewhere, there have been
several previous applications of full 3PCF and 4PCF to galaxies; here we
outline future possibilities. sarabande can be applied to galaxy redshift
surveys by obtaining either spectroscopic redshifts (BOSS collaboration et
al., 2017; eBOSS collaboration et al., 2021; DESI Collaboration et al., 2016)
or photometric redshifts (LSST Dark Energy Science Collaboration, 2012). In
addition to surveys with individually-resolved galaxies, sarabande can also be
applied to surveys using unresolved galaxies or the diffuse intergalactic
medium (FAST collaboration et al., 2011; LOFAR collaboration et al., 2013;
SPHEREX collaboration et al., 2014; CHIME collaboration et al., 2014; BINGO
collaboration et al., 2016) by measuring the intensity of a chosen emission
line. These surveys usually cover a large sky fraction but with relatively low
resolution compared to surveys that resolve individual galaxies. Within each
redshift shell, the measured intensity is a projected quantity; therefore,
sarabande with its projected 3/4 PCF implementation is useful for such
surveys.

Finally, in many cosmology settings, one computes the desired clustering
statistic on many mock catalogs to estimate covariance. Often these mocks are
grid-based and here using sarabande to evaluate the covariance would be
lossless, and the speed enabling. This would aid, for instance, a DESI 3PCF or
4PCF search for BAO (already detected in the 3PCF and bispectrum of DESI’s
predecessor Sloan Digital Sky Survey Baryon Oscillation Spectroscopic Survey
(BOSS), Slepian et al. 2017a, Slepian et al. 2017b, Pearson & Samushia 2018).
Especially given the development of fast methods for other parts of the BAO
fitting pipeline, e.g. Hansen et al. (2021), speeding up the covariance
calculation, which is often rate-limiting, is worthwhile. Given that higher-
order statistics such as 3PCF also have been shown to carry extra information
on the neutrino mass (Hahn et al. 2020, Hahn & Villaescusa-Navarro 2021,
Kamalinejad & Slepian 2020, Aviles et al. 2021) and on modified gravity (Alam
et al., 2021), two other major goals of DESI, new fast algorithms such as
sarabande are desirable.

Figure 8: Here we show an example of the fraction of time spent for each of
the different computations in the code with a given set of code parameters.
For the full 3/4 PCFs, the computational time is dominated by the computation
of the convolution coefficients (Algorithms 1 and 3) when our choice of
$\ell_{\rm max}$ and $N_{\rm bins}$ are small. As $\ell_{\rm max}$ and $N_{\rm
bins}$ increase, the number of File I/O operations increases causing the total
computational time of the full 3/4 PCFs to become dominated by combining the
coefficients (Algorithms 5 and 6). The computation of the convolution
coefficients in sarabande scales in complexity with the number of grid cells
used in the input data as $\mathcal{O}(N_{\rm g}\log N_{\rm g})$ as shown in
Figure 5. This scaling is important because this part of the total computation
is dominant for high resolution data measurements. The total computation time
of the projected 3/4 PCFs is dominated by the combination of convolution
coefficients (Algorithms 7 and 8). It is expected that the computation of the
convolution coefficients (Algorithms 2 and 4) does not dominate the total
computation time for the projected 3/4 PCFs because the number of grid cells
$N_{\rm g}$ is significantly smaller in 2D compared to 3D. This reduction in
the number of grid cells of the input data allows us to avoid File I/O
operations entirely. The combination of convolution coefficients is more
dominant in the projected 4PCF than the projected 3PCF due to the additional
nested for loops. The creation of kernels is nearly negligible for the
projected pathway, while it is more significant in the full pathway. For these
timings we used we used $N_{\rm bins}=5$, $\ell_{\textrm{max}}=2$, and a grid
resolution of $N_{\rm g}=128^{3}$ for the full 3/4 PCFs. For the timing of the
projected computation we used $N_{\rm bins}=5$, $m_{\textrm{max}}=5$, and a
grid resolution of $N_{\rm g}=128^{2}$.

##  6 Conclusion

sarabande is a new python package for measuring 3 and 4 Point Correlation
Functions on gridded data using Fast Fourier Transforms. The use of FFTs gives
sarabande a complexity scaling of $\mathcal{O}(N_{\mathrm{g}}\log
N_{\mathrm{g}})$ which offers the fastest method for measuring the 3/4 PCF to
date. The package allows for measuring the standard full 3/4 PCFs on a 3D data
cube while also giving the option to measure their projected counterparts on a
2D data slice. Since this package operates on gridded data, it allows users to
measure the 3/4 PCF in new contexts yet to be explored. The standard
algorithms use a particle-based approach by directly counting triplets and
quadruplets, which is slower than our FFT based approach which computes the
contributions of triplets or quadruplets simultaneously at every point
instead. This older direct-counting method was developed in the context of
cosmological surveys where galaxies are treated as points in a survey volume.
Now with sarabande we demonstrate how one can measure the 3/4 PCF and their
projected counterparts on datasets outside of cosmology. We give an example
measurement on a simulation of the the turbulent ISM provided by the CATS
collaboration (Bialy & Burkhart, 2020; Burkhart et al., 2009; Cho & Lazarian,
2003; Portillo et al., 2018).

In the future we plan to explore how the 3 and 4 Point Correlation Functions
can be applied to MHD simulations to better understand turbulence, shocks, and
other physical processes. We also plan to continue optimizing the code by
implementing more parallelization schemes and potentially harnessing GPUs for
further acceleration (e.g. Artiles & Saeed (2019)).

## Acknowledgements

We acknowledge the University of Florida Research Computing for providing
computational resources and support that have contributed to the research
results reported in this publication. We would like to thank Oliver Philcox,
Andrew Saydjari, and Matt Hansen for insightful discussions pertaining to the
development and optimization of sarabande. Support from Isabella Dougas and
Ronan Hix was crucial during the development of sarabande. Lastly, we thank
the rest of the Slepian Research Group at University of Florida for their
feedback and support on the project. Jiamin Hou 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 research would not be
possible without the funding provided by the National Science Foundation for
the Research Experience for Undergraduates at University of Florida.

## Data Availability

All data used for this work are publicly available and can be found on the
sarabande github.101010sarabande github:
https://github.com/James11222/sarabande

## References

  * Alam et al. (2021) Alam, S., Arnold, C., Aviles, A., Bean, R., Cai, Y.-C., Cautun, M., Cervantes-Cota, J. L., Cuesta-Lazaro, C., Devi, N. C., Eggemeier, A., Fromenteau, S., Gonzalez-Morales, A. X., Halenka, V., He, J.-h., Hellwing, W. A., Hernández-Aguayo, C., Ishak, M., Koyama, K., Li, B., de la Macorra, A., Meneses Rizo, J., Miller, C., Mueller, E.-M., Niz, G., Ntelis, P., Rodríguez Otero, M., Sabiu, C. G., Slepian, Z., Stark, A., Valenzuela, O., Valogiannis, G., Vargas-Magaña, M., Winther, H. A., Zarrouk, P., Zhao, G.-B., & Zheng, Y., 2021.  Towards testing the theory of gravity with DESI: summary statistics, model predictions and future simulation requirements, J. Cosmology Astropart. Phys, 2021(11), 050.
  * Artiles & Saeed (2019) Artiles, O. & Saeed, F., 2019.  Gpu-sfft: A gpu based parallel algorithm for computing the sparse fast fourier transform (sfft) of k-sparse signals, in 2019 IEEE International Conference on Big Data (Big Data), pp. 3303–3311.
  * Aviles et al. (2021) Aviles, A., Banerjee, A., Niz, G., & Slepian, Z., 2021.  Clustering in massive neutrino cosmologies via Eulerian Perturbation Theory, J. Cosmology Astropart. Phys, 2021(11), 028.
  * Bialy & Burkhart (2020) Bialy, S. & Burkhart, B., 2020.  The Driving Scale-Density Decorrelation Scale Relation in a Turbulent Medium, ApJ, 894(1), L2.
  * BINGO collaboration et al. (2016) BINGO collaboration, Battye, R., Browne, I., Chen, T., Dickinson, C., Harper, S., Olivari, L., Peel, M., Remazeilles, M., Roychowdhury, S., Wilkinson, P., Abdalla, E., Abramo, R., Ferreira, E., Wuensche, A., Vilella, T., Caldas, M., Tancredi, G., Refregier, A., Monstein, C., Abdalla, F., Pourtsidou, A., Maffei, B., Pisano, G., & Ma, Y.-Z., 2016.  Update on the BINGO 21cm intensity mapping experiment, arXiv e-prints, p. arXiv:1610.06826.
  * BOSS collaboration et al. (2017) BOSS collaboration, Alam, S., Ata, M., Bailey, S., Beutler, F., Bizyaev, D., Blazek, J. A., Bolton, A. S., Brownstein, J. R., Burden, A., Chuang, C.-H., Comparat, J., Cuesta, A. J., Dawson, K. S., Eisenstein, D. J., Escoffier, S., Gil-Marín, H., Grieb, J. N., Hand, N., Ho, S., Kinemuchi, K., Kirkby, D., Kitaura, F., Malanushenko, E., Malanushenko, V., Maraston, C., McBride, C. K., Nichol, R. C., Olmstead, M. D., Oravetz, D., Padmanabhan, N., Palanque-Delabrouille, N., Pan, K., Pellejero-Ibanez, M., Percival, W. J., Petitjean, P., Prada, F., Price-Whelan, A. M., Reid, B. A., Rodríguez-Torres, S. A., Roe, N. A., Ross, A. J., Ross, N. P., Rossi, G., Rubiño-Martín, J. A., Saito, S., Salazar-Albornoz, S., Samushia, L., Sánchez, A. G., Satpathy, S., Schlegel, D. J., Schneider, D. P., Scóccola, C. G., Seo, H.-J., Sheldon, E. S., Simmons, A., Slosar, A., Strauss, M. A., Swanson, M. E. C., Thomas, D., Tinker, J. L., Tojeiro, R., Magaña, M. V., Vazquez, J. A., Verde, L., Wake, D. A., Wang, Y., Weinberg, D. H., White, M., Wood-Vasey, W. M., Yèche, C., Zehavi, I., Zhai, Z., & Zhao, G.-B., 2017.  The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, MNRAS, 470(3), 2617–2652.
  * Burkhart & Lazarian (2010) Burkhart, B. & Lazarian, A., 2010.  Statistical tools of interstellar turbulence: connecting observations with theory, Proceedings of the International Astronomical Union, 6(S274), 365–368.
  * Burkhart et al. (2009) Burkhart, B., Falceta-Gonçalves, D., Kowal, G., & Lazarian, A., 2009\.  Density Studies of MHD Interstellar Turbulence: Statistical Moments, Correlations and Bispectrum, ApJ, 693(1), 250–266.
  * Burkhart et al. (2020) Burkhart, B., Appel, S. M., Bialy, S., Cho, J., Christensen, A. J., Collins, D., Federrath, C., Fielding, D. B., Finkbeiner, D., Hill, A. S., Ibáñez-Mejía, J. C., Krumholz, M. R., Lazarian, A., Li, M., Mocz, P., Mac Low, M. M., Naiman, J., Portillo, S. K. N., Shane, B., Slepian, Z., & Yuan, Y., 2020.  The Catalogue for Astrophysical Turbulence Simulations (CATS), ApJ, 905(1), 14.
  * Cahn & Slepian (2020) Cahn, R. N. & Slepian, Z., 2020.  Isotropic N-Point Basis Functions and Their Properties, arXiv e-prints, p. arXiv:2010.14418.
  * Cahn et al. (2021) Cahn, R. N., Slepian, Z., & Hou, J., 2021.  A Test for Cosmological Parity Violation Using the 3D Distribution of Galaxies, arXiv e-prints, p. arXiv:2110.12004.
  * CHIME collaboration et al. (2014) CHIME collaboration, Shaw, J. R., Sigurdson, K., Pen, U.-L., Stebbins, A., & Sitwell, M., 2014.  All-sky Interferometry with Spherical Harmonic Transit Telescopes, ApJ, 781(2), 57.
  * Cho & Lazarian (2003) Cho, J. & Lazarian, A., 2003.  Compressible magnetohydrodynamic turbulence: mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications, MNRAS, 345(12), 325–339.
  * Collis et al. (1998) Collis, W., White, P., & Hammond, J., 1998.  Higher-order spectra: the bispectrum and trispectrum, Mechanical systems and signal processing, 12(3), 375–394.
  * Cooley & Tukey (1965) Cooley, J. W. & Tukey, J. W., 1965.  An algorithm for the machine calculation of complex fourier series, Mathematics of computation, 19(90), 297–301.
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., Anderson, S. F., Aubourg, É., Bailey, S., Barkhouser, R. H., Bautista, J. E., Beifiori, A., Berlind, A. A., Bhardwaj, V., Bizyaev, D., Blake, C. H., Blanton, M. R., Blomqvist, M., Bolton, A. S., Borde, A., Bovy, J., Brandt, W. N., Brewington, H., Brinkmann, J., Brown, P. J., Brownstein, J. R., Bundy, K., Busca, N. G., Carithers, W., Carnero, A. R., Carr, M. A., Chen, Y., Comparat, J., Connolly, N., Cope, F., Croft, R. A. C., Cuesta, A. J., da Costa, L. N., Davenport, J. R. A., Delubac, T., de Putter, R., Dhital, S., Ealet, A., Ebelke, G. L., Eisenstein, D. J., Escoffier, S., Fan, X., Filiz Ak, N., Finley, H., Font-Ribera, A., Génova-Santos, R., Gunn, J. E., Guo, H., Haggard, D., Hall, P. B., Hamilton, J.-C., Harris, B., Harris, D. W., Ho, S., Hogg, D. W., Holder, D., Honscheid, K., Huehnerhoff, J., Jordan, B., Jordan, W. P., Kauffmann, G., Kazin, E. A., Kirkby, D., Klaene, M. A., Kneib, J.-P., Le Goff, J.-M., Lee, K.-G., Long, D. C., Loomis, C. P., Lundgren, B., Lupton, R. H., Maia, M. A. G., Makler, M., Malanushenko, E., Malanushenko, V., Mandelbaum, R., Manera, M., Maraston, C., Margala, D., Masters, K. L., McBride, C. K., McDonald, P., McGreer, I. D., McMahon, R. G., Mena, O., Miralda-Escudé, J., Montero-Dorta, A. D., Montesano, F., Muna, D., Myers, A. D., Naugle, T., Nichol, R. C., Noterdaeme, P., Nuza, S. E., Olmstead, M. D., Oravetz, A., Oravetz, D. J., Owen, R., Padmanabhan, N., Palanque-Delabrouille, N., Pan, K., Parejko, J. K., Pâris, I., Percival, W. J., Pérez-Fournon, I., Pérez-Ràfols, I., Petitjean, P., Pfaffenberger, R., Pforr, J., Pieri, M. M., Prada, F., Price-Whelan, A. M., Raddick, M. J., Rebolo, R., Rich, J., Richards, G. T., Rockosi, C. M., Roe, N. A., Ross, A. J., Ross, N. P., Rossi, G., Rubiño-Martin, J. A., Samushia, L., Sánchez, A. G., Sayres, C., Schmidt, S. J., Schneider, D. P., Scóccola, C. G., Seo, H.-J., Shelden, A., Sheldon, E., Shen, Y., Shu, Y., Slosar, A., Smee, S. A., Snedden, S. A., Stauffer, F., Steele, O., Strauss, M. A., Streblyanska, A., Suzuki, N., Swanson, M. E. C., Tal, T., Tanaka, M., Thomas, D., Tinker, J. L., Tojeiro, R., Tremonti, C. A., Vargas Magaña, M., Verde, L., Viel, M., Wake, D. A., Watson, M., Weaver, B. A., Weinberg, D. H., Weiner, B. J., West, A. A., White, M., Wood-Vasey, W. M., Yeche, C., Zehavi, I., Zhao, G.-B., & Zheng, Z., 2013.  The Baryon Oscillation Spectroscopic Survey of SDSS-III, AJ, 145(1), 10.
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., Ahlen, S., Alam, S., Allen, L. E., Allende Prieto, C., Annis, J., Bailey, S., Balland, C., Ballester, O., Baltay, C., Beaufore, L., Bebek, C., Beers, T. C., Bell, E. F., Bernal, J. L., Besuner, R., Beutler, F., Blake, C., Bleuler, H., Blomqvist, M., Blum, R., Bolton, A. S., Briceno, C., Brooks, D., Brownstein, J. R., Buckley-Geer, E., Burden, A., Burtin, E., Busca, N. G., Cahn, R. N., Cai, Y.-C., Cardiel-Sas, L., Carlberg, R. G., Carton, P.-H., Casas, R., Castander, F. J., Cervantes-Cota, J. L., Claybaugh, T. M., Close, M., Coker, C. T., Cole, S., Comparat, J., Cooper, A. P., Cousinou, M. C., Crocce, M., Cuby, J.-G., Cunningham, D. P., Davis, T. M., Dawson, K. S., de la Macorra, A., De Vicente, J., Delubac, T., Derwent, M., Dey, A., Dhungana, G., Ding, Z., Doel, P., Duan, Y. T., Ealet, A., Edelstein, J., Eftekharzadeh, S., Eisenstein, D. J., Elliott, A., Escoffier, S., Evatt, M., Fagrelius, P., Fan, X., Fanning, K., Farahi, A., Farihi, J., Favole, G., Feng, Y., Fernandez, E., Findlay, J. R., Finkbeiner, D. P., Fitzpatrick, M. J., Flaugher, B., Flender, S., Font-Ribera, A., Forero-Romero, J. E., Fosalba, P., Frenk, C. S., Fumagalli, M., Gaensicke, B. T., Gallo, G., Garcia-Bellido, J., Gaztanaga, E., Pietro Gentile Fusillo, N., Gerard, T., Gershkovich, I., Giannantonio, T., Gillet, D., Gonzalez-de-Rivera, G., Gonzalez-Perez, V., Gott, S., Graur, O., Gutierrez, G., Guy, J., Habib, S., Heetderks, H., Heetderks, I., Heitmann, K., Hellwing, W. A., Herrera, D. A., Ho, S., Holland, S., Honscheid, K., Huff, E., Hutchinson, T. A., Huterer, D., Hwang, H. S., Illa Laguna, J. M., Ishikawa, Y., Jacobs, D., Jeffrey, N., Jelinsky, P., Jennings, E., Jiang, L., Jimenez, J., Johnson, J., Joyce, R., Jullo, E., Juneau, S., Kama, S., Karcher, A., Karkar, S., Kehoe, R., Kennamer, N., Kent, S., Kilbinger, M., Kim, A. G., Kirkby, D., Kisner, T., Kitanidis, E., Kneib, J.-P., Koposov, S., Kovacs, E., Koyama, K., Kremin, A., Kron, R., Kronig, L., Kueter-Young, A., Lacey, C. G., Lafever, R., Lahav, O., Lambert, A., Lampton, M., Landriau, M., Lang, D., Lauer, T. R., Le Goff, J.-M., Le Guillou, L., Le Van Suu, A., Lee, J. H., Lee, S.-J., Leitner, D., Lesser, M., Levi, M. E., L’Huillier, B., Li, B., Liang, M., Lin, H., Linder, E., Loebman, S. R., Lukić, Z., Ma, J., MacCrann, N., Magneville, C., Makarem, L., Manera, M., Manser, C. J., Marshall, R., Martini, P., Massey, R., Matheson, T., McCauley, J., McDonald, P., McGreer, I. D., Meisner, A., Metcalfe, N., Miller, T. N., Miquel, R., Moustakas, J., Myers, A., Naik, M., Newman, J. A., Nichol, R. C., Nicola, A., Nicolati da Costa, L., Nie, J., Niz, G., Norberg, P., Nord, B., Norman, D., Nugent, P., O’Brien, T., Oh, M., Olsen, K. A. G., Padilla, C., Padmanabhan, H., Padmanabhan, N., Palanque-Delabrouille, N., Palmese, A., Pappalardo, D., Pâris, I., Park, C., Patej, A., Peacock, J. A., Peiris, H. V., Peng, X., Percival, W. J., Perruchot, S., Pieri, M. M., Pogge, R., Pollack, J. E., Poppett, C., Prada, F., Prakash, A., Probst, R. G., Rabinowitz, D., Raichoor, A., Ree, C. H., Refregier, A., Regal, X., Reid, B., Reil, K., Rezaie, M., Rockosi, C. M., Roe, N., Ronayette, S., Roodman, A., Ross, A. J., Ross, N. P., Rossi, G., Rozo, E., Ruhlmann-Kleider, V., Rykoff, E. S., Sabiu, C., Samushia, L., Sanchez, E., Sanchez, J., Schlegel, D. J., Schneider, M., Schubnell, M., Secroun, A., Seljak, U., Seo, H.-J., Serrano, S., Shafieloo, A., Shan, H., Sharples, R., Sholl, M. J., Shourt, W. V., Silber, J. H., Silva, D. R., Sirk, M. M., Slosar, A., Smith, A., Smoot, G. F., Som, D., Song, Y.-S., Sprayberry, D., Staten, R., Stefanik, A., Tarle, G., Sien Tie, S., Tinker, J. L., Tojeiro, R., Valdes, F., Valenzuela, O., Valluri, M., Vargas-Magana, M., Verde, L., Walker, A. R., Wang, J., Wang, Y., Weaver, B. A., Weaverdyck, C., Wechsler, R. H., Weinberg, D. H., White, M., Yang, Q., Yeche, C., Zhang, T., Zhao, G.-B., Zheng, Y., Zhou, X., Zhou, Z., Zhu, Y., Zou, H., & Zu, Y., 2016.  The DESI Experiment Part I: Science,Targeting, and Survey Design, arXiv e-prints, p. arXiv:1611.00036.
  * eBOSS collaboration et al. (2021) eBOSS collaboration, Alam, S., Aubert, M., Avila, S., Balland, C., Bautista, J. E., Bershady, M. A., Bizyaev, D., Blanton, M. R., Bolton, A. S., Bovy, J., Brinkmann, J., Brownstein, J. R., Burtin, E., Chabanier, S., Chapman, M. J., Choi, P. D., Chuang, C.-H., Comparat, J., Cousinou, M.-C., Cuceu, A., Dawson, K. S., de la Torre, S., de Mattia, A., Agathe, V. d. S., des Bourboux, H. d. M., Escoffier, S., Etourneau, T., Farr, J., Font-Ribera, A., Frinchaboy, P. M., Fromenteau, S., Gil-Marín, H., Le Goff, J.-M., Gonzalez-Morales, A. X., Gonzalez-Perez, V., Grabowski, K., Guy, J., Hawken, A. J., Hou, J., Kong, H., Parker, J., Klaene, M., Kneib, J.-P., Lin, S., Long, D., Lyke, B. W., de la Macorra, A., Martini, P., Masters, K., Mohammad, F. G., Moon, J., Mueller, E.-M., Muñoz-Gutiérrez, A., Myers, A. D., Nadathur, S., Neveux, R., Newman, J. A., Noterdaeme, P., Oravetz, A., Oravetz, D., Palanque-Delabrouille, N., Pan, K., Paviot, R., Percival, W. J., Pérez-Ràfols, I., Petitjean, P., Pieri, M. M., Prakash, A., Raichoor, A., Ravoux, C., Rezaie, M., Rich, J., Ross, A. J., Rossi, G., Ruggeri, R., Ruhlmann-Kleider, V., Sánchez, A. G., Sánchez, F. J., Sánchez-Gallego, J. R., Sayres, C., Schneider, D. P., Seo, H.-J., Shafieloo, A., Slosar, A., Smith, A., Stermer, J., Tamone, A., Tinker, J. L., Tojeiro, R., Vargas-Magaña, M., Variu, A., Wang, Y., Weaver, B. A., Weijmans, A.-M., Yèche, C., Zarrouk, P., Zhao, C., Zhao, G.-B., & Zheng, Z., 2021.  Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory, Phys. Rev. D, 103(8), 083533.
  * Eisenstein et al. (2011) Eisenstein, D. J., Weinberg, D. H., Agol, E., Aihara, H., Allende Prieto, C., Anderson, S. F., Arns, J. A., Aubourg, É., Bailey, S., Balbinot, E., Barkhouser, R., Beers, T. C., Berlind, A. A., Bickerton, S. J., Bizyaev, D., Blanton, M. R., Bochanski, J. J., Bolton, A. S., Bosman, C. T., Bovy, J., Brandt, W. N., Breslauer, B., Brewington, H. J., Brinkmann, J., Brown, P. J., Brownstein, J. R., Burger, D., Busca, N. G., Campbell, H., Cargile, P. A., Carithers, W. C., Carlberg, J. K., Carr, M. A., Chang, L., Chen, Y., Chiappini, C., Comparat, J., Connolly, N., Cortes, M., Croft, R. A. C., Cunha, K., da Costa, L. N., Davenport, J. R. A., Dawson, K., De Lee, N., Porto de Mello, G. F., de Simoni, F., Dean, J., Dhital, S., Ealet, A., Ebelke, G. L., Edmondson, E. M., Eiting, J. M., Escoffier, S., Esposito, M., Evans, M. L., Fan, X., Femenía Castellá, B., Dutra Ferreira, L., Fitzgerald, G., Fleming, S. W., Font-Ribera, A., Ford, E. B., Frinchaboy, P. M., García Pérez, A. E., Gaudi, B. S., Ge, J., Ghezzi, L., Gillespie, B. A., Gilmore, G., Girardi, L., Gott, J. R., Gould, A., Grebel, E. K., Gunn, J. E., Hamilton, J.-C., Harding, P., Harris, D. W., Hawley, S. L., Hearty, F. R., Hennawi, J. F., González Hernández, J. I., Ho, S., Hogg, D. W., Holtzman, J. A., Honscheid, K., Inada, N., Ivans, I. I., Jiang, L., Jiang, P., Johnson, J. A., Jordan, C., Jordan, W. P., Kauffmann, G., Kazin, E., Kirkby, D., Klaene, M. A., Knapp, G. R., Kneib, J.-P., Kochanek, C. S., Koesterke, L., Kollmeier, J. A., Kron, R. G., Lampeitl, H., Lang, D., Lawler, J. E., Le Goff, J.-M., Lee, B. L., Lee, Y. S., Leisenring, J. M., Lin, Y.-T., Liu, J., Long, D. C., Loomis, C. P., Lucatello, S., Lundgren, B., Lupton, R. H., Ma, B., Ma, Z., MacDonald, N., Mack, C., Mahadevan, S., Maia, M. A. G., Majewski, S. R., Makler, M., Malanushenko, E., Malanushenko, V., Mandelbaum, R., Maraston, C., Margala, D., Maseman, P., Masters, K. L., McBride, C. K., McDonald, P., McGreer, I. D., McMahon, R. G., Mena Requejo, O., Ménard, B., Miralda-Escudé, J., Morrison, H. L., Mullally, F., Muna, D., Murayama, H., Myers, A. D., Naugle, T., Neto, A. F., Nguyen, D. C., Nichol, R. C., Nidever, D. L., O’Connell, R. W., Ogando, R. L. C., Olmstead, M. D., Oravetz, D. J., Padmanabhan, N., Paegert, M., Palanque-Delabrouille, N., Pan, K., Pandey, P., Parejko, J. K., Pâris, I., Pellegrini, P., Pepper, J., Percival, W. J., Petitjean, P., Pfaffenberger, R., Pforr, J., Phleps, S., Pichon, C., Pieri, M. M., Prada, F., Price-Whelan, A. M., Raddick, M. J., Ramos, B. H. F., Reid, I. N., Reyle, C., Rich, J., Richards, G. T., Rieke, G. H., Rieke, M. J., Rix, H.-W., Robin, A. C., Rocha-Pinto, H. J., Rockosi, C. M., Roe, N. A., Rollinde, E., Ross, A. J., Ross, N. P., Rossetto, B., Sánchez, A. G., Santiago, B., Sayres, C., Schiavon, R., Schlegel, D. J., Schlesinger, K. J., Schmidt, S. J., Schneider, D. P., Sellgren, K., Shelden, A., Sheldon, E., Shetrone, M., Shu, Y., Silverman, J. D., Simmerer, J., Simmons, A. E., Sivarani, T., Skrutskie, M. F., Slosar, A., Smee, S., Smith, V. V., Snedden, S. A., Stassun, K. G., Steele, O., Steinmetz, M., Stockett, M. H., Stollberg, T., Strauss, M. A., Szalay, A. S., Tanaka, M., Thakar, A. R., Thomas, D., Tinker, J. L., Tofflemire, B. M., Tojeiro, R., Tremonti, C. A., Vargas Magaña, M., Verde, L., Vogt, N. P., Wake, D. A., Wan, X., Wang, J., Weaver, B. A., White, M., White, S. D. M., Wilson, J. C., Wisniewski, J. P., Wood-Vasey, W. M., Yanny, B., Yasuda, N., Yèche, C., York, D. G., Young, E., Zasowski, G., Zehavi, I., & Zhao, B., 2011.  SDSS-III: Massive Spectroscopic Surveys of the Distant Universe, the Milky Way, and Extra-Solar Planetary Systems, AJ, 142(3), 72.
  * FAST collaboration et al. (2011) FAST collaboration, Nan, R., Li, D., Jin, C., Wang, Q., Zhu, L., Zhu, W., Zhang, H., Yue, Y., & Qian, L., 2011.  The Five-Hundred Aperture Spherical Radio Telescope (fast) Project, International Journal of Modern Physics D, 20(6), 989–1024.
  * Fry & Peebles (1978) Fry, J. N. & Peebles, P. J. E., 1978.  Statistical analysis of catalogs of extragalactic objects. IX. The four-point galaxy correlation function., ApJ, 221, 19–33.
  * Gualdi et al. (2021) Gualdi, D., Novell, S., Gil-Marín, H., & Verde, L., 2021.  Matter trispectrum: theoretical modelling and comparison to n-body simulations, Journal of Cosmology and Astroparticle Physics, 2021(01), 015–015.
  * Hahn & Villaescusa-Navarro (2021) Hahn, C. & Villaescusa-Navarro, F., 2021.  Constraining Mν with the bispectrum. Part II. The information content of the galaxy bispectrum monopole, J. Cosmology Astropart. Phys, 2021(4), 029.
  * Hahn et al. (2020) Hahn, C., Villaescusa-Navarro, F., Castorina, E., & Scoccimarro, R., 2020\.  Constraining Mν with the bispectrum. Part I. Breaking parameter degeneracies, J. Cosmology Astropart. Phys, 2020(3), 040.
  * Hand et al. (2017) Hand, N., Li, Y., Slepian, Z., & Seljak, U., 2017.  An optimal FFT-based anisotropic power spectrum estimator, J. Cosmology Astropart. Phys, 2017(7), 002.
  * Hand et al. (2018) Hand, N., Feng, Y., Beutler, F., Li, Y., Modi, C., Seljak, U., & Slepian, Z., 2018.  nbodykit: An Open-source, Massively Parallel Toolkit for Large-scale Structure, AJ, 156(4), 160.
  * Hansen et al. (2021) Hansen, M., Krolewski, A., & Slepian, Z., 2021.  Accelerating BAO Scale Fitting Using Taylor Series, arXiv e-prints, p. arXiv:2112.06438.
  * Heideman et al. (1984) Heideman, M., Johnson, D., & Burrus, C., 1984.  Gauss and the history of the fast fourier transform, IEEE ASSP Magazine, 1(4), 14–21.
  * Hou et al. (2022a) Hou, J., Cahn, R. N., Philcox, O. H. E., & Slepian, Z., 2022a.  Analytic Gaussian covariance matrices for galaxy N -point correlation functions, Phys. Rev. D, 106(4), 043515.
  * Hou et al. (2022b) Hou, J., Slepian, Z., & Cahn, R. N., 2022b.  Measurement of Parity-Odd Modes in the Large-Scale 4-Point Correlation Function of SDSS BOSS DR12 CMASS and LOWZ Galaxies, arXiv e-prints, p. arXiv:2206.03625.
  * Jing (2005) Jing, Y. P., 2005.  Correcting for the Alias Effect When Measuring the Power Spectrum Using a Fast Fourier Transform, ApJ, 620(2), 559–563.
  * Kamalinejad & Slepian (2020) Kamalinejad, F. & Slepian, Z., 2020.  A Non-Degenerate Neutrino Mass Signature in the Galaxy Bispectrum, arXiv e-prints, p. arXiv:2011.00899.
  * Koch et al. (2019) Koch, E. W., Rosolowsky, E. W., Boyden, R. D., Burkhart, B., Ginsburg, A., Loeppky, J. L., & Offner, S. S. R., 2019.  TURBUSTAT: Turbulence Statistics in Python, AJ, 158(1), 1.
  * LOFAR collaboration et al. (2013) LOFAR collaboration, van Haarlem, M. P., Wise, M. W., Gunst, A. W., Heald, G., McKean, J. P., Hessels, J. W. T., de Bruyn, A. G., Nijboer, R., Swinbank, J., Fallows, R., Brentjens, M., Nelles, A., Beck, R., Falcke, H., Fender, R., Hörandel, J., Koopmans, L. V. E., Mann, G., Miley, G., Röttgering, H., Stappers, B. W., Wijers, R. A. M. J., Zaroubi, S., van den Akker, M., Alexov, A., Anderson, J., Anderson, K., van Ardenne, A., Arts, M., Asgekar, A., Avruch, I. M., Batejat, F., Bähren, L., Bell, M. E., Bell, M. R., van Bemmel, I., Bennema, P., Bentum, M. J., Bernardi, G., Best, P., Bîrzan, L., Bonafede, A., Boonstra, A. J., Braun, R., Bregman, J., Breitling, F., van de Brink, R. H., Broderick, J., Broekema, P. C., Brouw, W. N., Brüggen, M., Butcher, H. R., van Cappellen, W., Ciardi, B., Coenen, T., Conway, J., Coolen, A., Corstanje, A., Damstra, S., Davies, O., Deller, A. T., Dettmar, R. J., van Diepen, G., Dijkstra, K., Donker, P., Doorduin, A., Dromer, J., Drost, M., van Duin, A., Eislöffel, J., van Enst, J., Ferrari, C., Frieswijk, W., Gankema, H., Garrett, M. A., de Gasperin, F., Gerbers, M., de Geus, E., Grießmeier, J. M., Grit, T., Gruppen, P., Hamaker, J. P., Hassall, T., Hoeft, M., Holties, H. A., Horneffer, A., van der Horst, A., van Houwelingen, A., Huijgen, A., Iacobelli, M., Intema, H., Jackson, N., Jelic, V., de Jong, A., Juette, E., Kant, D., Karastergiou, A., Koers, A., Kollen, H., Kondratiev, V. I., Kooistra, E., Koopman, Y., Koster, A., Kuniyoshi, M., Kramer, M., Kuper, G., Lambropoulos, P., Law, C., van Leeuwen, J., Lemaitre, J., Loose, M., Maat, P., Macario, G., Markoff, S., Masters, J., McFadden, R. A., McKay-Bukowski, D., Meijering, H., Meulman, H., Mevius, M., Middelberg, E., Millenaar, R., Miller-Jones, J. C. A., Mohan, R. N., Mol, J. D., Morawietz, J., Morganti, R., Mulcahy, D. D., Mulder, E., Munk, H., Nieuwenhuis, L., van Nieuwpoort, R., Noordam, J. E., Norden, M., Noutsos, A., Offringa, A. R., Olofsson, H., Omar, A., Orrú, E., Overeem, R., Paas, H., Pandey-Pommier, M., Pandey, V. N., Pizzo, R., Polatidis, A., Rafferty, D., Rawlings, S., Reich, W., de Reijer, J. P., Reitsma, J., Renting, G. A., Riemers, P., Rol, E., Romein, J. W., Roosjen, J., Ruiter, M., Scaife, A., van der Schaaf, K., Scheers, B., Schellart, P., Schoenmakers, A., Schoonderbeek, G., Serylak, M., Shulevski, A., Sluman, J., Smirnov, O., Sobey, C., Spreeuw, H., Steinmetz, M., Sterks, C. G. M., Stiepel, H. J., Stuurwold, K., Tagger, M., Tang, Y., Tasse, C., Thomas, I., Thoudam, S., Toribio, M. C., van der Tol, B., Usov, O., van Veelen, M., van der Veen, A. J., ter Veen, S., Verbiest, J. P. W., Vermeulen, R., Vermaas, N., Vocks, C., Vogt, C., de Vos, M., van der Wal, E., van Weeren, R., Weggemans, H., Weltevrede, P., White, S., Wijnholds, S. J., Wilhelmsson, T., Wucknitz, O., Yatawatta, S., Zarka, P., Zensus, A., & van Zwieten, J., 2013.  LOFAR: The LOw-Frequency ARray, A&A, 556, A2.
  * LSST Dark Energy Science Collaboration (2012) LSST Dark Energy Science Collaboration, 2012.  Large Synoptic Survey Telescope: Dark Energy Science Collaboration,, ArXiv e-prints.
  * O’Brien et al. (2022) O’Brien, M. J., Burkhart, B., & Shelley, M. J., 2022.  Studying Interstellar Turbulence Driving Scales Using the Bispectrum, ApJ, 930(2), 149.
  * Pearson & Samushia (2018) Pearson, D. W. & Samushia, L., 2018.  A Detection of the Baryon Acoustic Oscillation features in the SDSS BOSS DR12 Galaxy Bispectrum, MNRAS, 478(4), 4500–4512.
  * Peebles (2001) Peebles, P., 2001.  The galaxy and mass n-point correlation functions: a blast from the past, arXiv preprint astro-ph/0103040.
  * Peebles & Ratra (2003) Peebles, P. J. E. & Ratra, B., 2003.  The cosmological constant and dark energy, Reviews of Modern Physics, 75(2), 559–606.
  * Philcox & Slepian (2021) Philcox, O. H. E. & Slepian, Z., 2021.  Beyond the Yamamoto approximation: Anisotropic power spectra and correlation functions with pairwise lines of sight, Phys. Rev. D, 103(12), 123509.
  * Philcox & Slepian (2022) Philcox, O. H. E. & Slepian, Z., 2022.  Efficient computation of N-point correlation functions in D dimensions, Proceedings of the National Academy of Science, 119(33), e2111366119.
  * Philcox et al. (2022) Philcox, O. H. E., Slepian, Z., Hou, J., Warner, C., Cahn, R. N., & Eisenstein, D. J., 2022.  ENCORE: an O (Ng2) estimator for galaxy N-point correlation functions, MNRAS, 509(2), 2457–2481.
  * Portillo et al. (2018) Portillo, S. K. N., Slepian, Z., Burkhart, B., Kahraman, S., & Finkbeiner, D. P., 2018.  Developing the 3-point Correlation Function for the Turbulent Interstellar Medium, ApJ, 862(2), 119.
  * Portillo et al. (2018) Portillo, S. K. N., Slepian, Z., Burkhart, B., Kahraman, S., & Finkbeiner, D. P., 2018.  Developing the 3-point correlation function for the turbulent interstellar medium, The Astrophysical Journal, 862(2), 119.
  * Sabiu et al. (2019) Sabiu, C. G., Hoyle, B., Kim, J., & Li, X.-D., 2019.  Graph Database Solution for Higher-order Spatial Statistics in the Era of Big Data, ApJS, 242(2), 29.
  * Samushia et al. (2021) Samushia, L., Slepian, Z., & Villaescusa-Navarro, F., 2021.  Information content of higher order galaxy correlation functions, MNRAS, 505(1), 628–641.
  * Saydjari et al. (2021) Saydjari, A. K., Portillo, S. K. N., Slepian, Z., Kahraman, S., Burkhart, B., & Finkbeiner, D. P., 2021.  Classification of magnetohydrodynamic simulations using wavelet scattering transforms, The Astrophysical Journal, 910(2), 122.
  * Slepian & Eisenstein (2015a) Slepian, Z. & Eisenstein, D. J., 2015a.  Computing the three-point correlation function of galaxies in O(N^2) time, MNRAS, 454, 4142–4158.
  * Slepian & Eisenstein (2015b) Slepian, Z. & Eisenstein, D. J., 2015b.  A new look at lines of sight: using Fourier methods for the wide-angle anisotropic 2-point correlation function, ArXiv e-prints.
  * Slepian & Eisenstein (2016) Slepian, Z. & Eisenstein, D. J., 2016.  Accelerating the two-point and three-point galaxy correlation functions using Fourier transforms, MNRAS, 455, L31–L35.
  * Slepian & Eisenstein (2018) Slepian, Z. & Eisenstein, D. J., 2018.  A practical computational method for the anisotropic redshift-space three-point correlation function, MNRAS, 478(2), 1468–1483.
  * Slepian et al. (2017a) Slepian, Z., Eisenstein, D. J., Beutler, F., Chuang, C.-H., Cuesta, A. J., Ge, J., Gil-Marín, H., Ho, S., Kitaura, F.-S., McBride, C. K., Nichol, R. C., Percival, W. J., Rodríguez-Torres, S., Ross, A. J., Scoccimarro, R., Seo, H.-J., Tinker, J., Tojeiro, R., & Vargas-Magaña, M., 2017a.  The large-scale three-point correlation function of the SDSS BOSS DR12 CMASS galaxies, MNRAS, 468(1), 1070–1083.
  * Slepian et al. (2017b) Slepian, Z., Eisenstein, D. J., Brownstein, J. R., Chuang, C.-H., Gil-Marín, H., Ho, S., Kitaura, F.-S., Percival, W. J., Ross, A. J., Rossi, G., Seo, H.-J., Slosar, A., & Vargas-Magaña, M., 2017b.  Detection of baryon acoustic oscillation features in the large-scale three-point correlation function of SDSS BOSS DR12 CMASS galaxies, MNRAS, 469(2), 1738–1751.
  * SPHEREX collaboration et al. (2014) SPHEREX collaboration, Doré, O., Bock, J., Ashby, M., Capak, P., Cooray, A., de Putter, R., Eifler, T., Flagey, N., Gong, Y., Habib, S., Heitmann, K., Hirata, C., Jeong, W.-S., Katti, R., Korngut, P., Krause, E., Lee, D.-H., Masters, D., Mauskopf, P., Melnick, G., Mennesson, B., Nguyen, H., Öberg, K., Pullen, A., Raccanelli, A., Smith, R., Song, Y.-S., Tolls, V., Unwin, S., Venumadhav, T., Viero, M., Werner, M., & Zemcov, M., 2014.  Cosmology with the SPHEREX All-Sky Spectral Survey, arXiv e-prints, p. arXiv:1412.4872.

## Appendix

Here we display each algorithm used in sarabande. Algorithms 1-4 are the
algorithms used to prepare kernels from basis functions and generate
convolution coefficients. Algorithms 5-8 are the algorithms used to compile
convolution coefficients into their corresponding correlation functions. All
algorithms are written heuristically for clarity. For details of the python
version, visit the sarabande github. 111111sarabande github:
https://github.com/James11222/sarabande

Result: The creation of $N_{\mathrm{g}}\times N_{\mathrm{g}}\times
N_{\mathrm{g}}$ arrays indexed by degree $\ell$, order $m$, and bin index
$\mathrm{b}$. The arrays that represent the Fourier Transform of the binned
spherical harmonic kernels are then saved to disk.

Notes: We initialize three $N_{\mathrm{g}}\times N_{\mathrm{g}}\times
N_{\mathrm{g}}$ arrays $x,y,z$ containing those coordinates centered on zero
then use these arrays to compute spherical harmonics from Cartesian
representations.

for  _$\ell$ = 0 to $\ell_{\rm max}$_ do

for  _m = 0 to $\ell$_ do

$Y_{\ell m}\leftarrow Y_{\ell m}(x,y,z)$;

for  _$\mathrm{b}$ = 1 to $N_{\rm bins}$_ do

$\Phi^{\mathrm{b}}\leftarrow r_{\rm
min}(\mathrm{b})<\sqrt{x^{2}+y^{2}+z^{2}}\leq r_{\rm max}(\mathrm{b})$;

$\tilde{Y}_{\ell m}^{\mathrm{b}}\leftarrow FFT(Y_{\ell
m}\times\Phi^{\mathrm{b}})$;

save $\tilde{Y}_{\ell m}^{\mathrm{b}}$ to disk;

end for

end for

end for

Algorithm 1 Generating Full 3/4 PCF Kernels

Result: The creation of $N_{\mathrm{g}}\times N_{\mathrm{g}}$ arrays that are
indexed by the angular index $m$ and the radial bin index $\mathrm{b}$. These
arrays represent the Fourier Transform of the binned Fourier basis kernels.

Notes: we initialize two $N_{\mathrm{g}}\times N_{\mathrm{g}}$ arrays $x,y$
containing those coordinates centered on zero, which we use in the following
identity for our basis: $e^{im\phi}=\left(\frac{x+iy}{r}\right)^{m}$

for  _$m=0$ to $m_{\mathrm{max}}$_ do

for  _$\mathrm{b}=1$ to $N_{\mathrm{bins}}$_ do

$\Phi^{\mathrm{b}}\leftarrow r_{\rm min}(\mathrm{b})<\sqrt{x^{2}+y^{2}}\leq
r_{\rm max}(\mathrm{b})$;

$K_{m}^{\mathrm{b}}\leftarrow FFT(e^{im\phi}\times\Phi^{\mathrm{b}}$) ;

end for

end for

Algorithm 2 Generate Projected 3/4 PCF Kernels

Result: Convolve binned spherical harmonic kernels to construct the $a_{\ell
m}^{\mathrm{b}}$ coefficients and save them to disk.

Notes: The $a_{\ell m}^{\mathrm{b}}$ coefficient arrays have the shape
$N_{\mathrm{g}}\times N_{\mathrm{g}}\times N_{\mathrm{g}}$ once indexed by
$\ell,m,\mathrm{b}$.

$\delta\leftarrow$ read in density field;

$\tilde{\delta}\leftarrow FFT(\delta)$

for  _$\ell$ = 0 to $\ell_{\rm max}$_ do

for  _m = 0 to $\ell$_ do

for  _$\mathrm{b}$ = 1 to $N_{\rm bins}$_ do

read $\tilde{Y}_{\ell m}^{\mathrm{b}}$ from disk;

$a_{\ell m}^{\mathrm{b}}\leftarrow
FFT^{-1}(\tilde{\delta}\times\tilde{Y}_{\ell m}^{\mathrm{b}})$;

save $a_{\ell m}^{\mathrm{b}}$ to disk;

end for

end for

end for

Algorithm 3 Convolve Full 3/4 PCF Kernels

Result: Convolve the data with the projected kernels to create
$c_{m}^{\mathrm{b}}$ coefficients which can be indexed by $m$ (angular index)
and $\mathrm{b}$ (radial bin index).

Notes: $c_{m}^{\mathrm{b}}$ coefficient arrays have the shape
$N_{\mathrm{g}}\times N_{\mathrm{g}}$ once indexed by $\mathrm{b}$ and $m$.

$\delta\leftarrow$ read in density field;

$\tilde{\delta}\leftarrow FFT(\delta)$;

for  _$m=0$ to $m_{\mathrm{max}}$_ do

for  _$\mathrm{b}=1$ to $N_{\mathrm{bins}}$_ do

$c_{m}^{\mathrm{b}}=FFT^{-1}(\tilde{\delta}\times K_{m}^{\mathrm{b}})$

end for

end for

Algorithm 4 Convolve Projected 3/4 PCF Kernels

Result: Computation of the unique full 3PCF coefficients from the $a_{\ell
m}^{\mathrm{b}}$ coefficients computed in algorithm 3. The end result is a 3D
array indexed by $\ell,\mathrm{b}_{1},\mathrm{b}_{2}$ where $\ell$ is the
order of the angular component of the triangles and $\mathrm{b}_{1}$ and
$\mathrm{b}_{2}$ denote the radial bin length of side 1 and 2 of the triangles
respectively.

Notes: We only calculate the unique coefficients since the structure of the
full 3PCF is symmetric under permutations of the indices $\mathrm{b}_{1}$ and
$\mathrm{b}_{2}$. To recover the full 3PCF one can just permute indices and
preserve corresponding values.

$\delta\leftarrow$ read in density cube

initialize $\zeta$ as an $(\ell_{\rm max}+1)\times N_{\rm bin}\times N_{\rm
bin}$ array of zeros;

for  _$\ell$ = 0 to $\ell_{\rm max}$_ do

for  _$\mathrm{b}_{1}$ = 1 to $N_{\rm bins}$_ do

for  _m = 0 to $\ell$_ do

load $a_{\ell m}^{\mathrm{b}_{1}}$ from disk;

for  _$\mathrm{b}_{2}$ = 1 to $\mathrm{b}_{1}$_ do

load $a_{\ell m}^{\mathrm{b}_{2}}$ from disk;

S $\leftarrow$ sum$\left(\delta\times a_{\ell m}^{\mathrm{b}_{1}}\times
a_{\ell m}^{\mathrm{b}_{2}*}\right)$;

if  _m > 0_ then

$S$ += $S^{*}$;

end if

$\zeta$[$\ell$, $\mathrm{b}_{1}$, $\mathrm{b}_{2}$] += S;

end for

end for

end for

end for

Algorithm 5 Combine Full 3PCF Coefficients

Result: Computation of the unique full 4PCF coefficients from the $a_{\ell
m}^{\mathrm{b}}$ coefficients computed in algorithm 3. The end result is a 6D
array indexed by
$\ell_{1},\ell_{2},\ell_{3},\mathrm{b}_{1},\mathrm{b}_{2},\mathrm{b}_{3}$
where $\ell_{i}$ is the order of the $i$th angular component of the
tetrahedrons and $\mathrm{b}_{1},\mathrm{b}_{2}$ and $\mathrm{b}_{3}$ denote
the radial bin length of each side of the tetrahedrons respectively.

Notes: We only calculate the unique coefficients since the structure of the
full 4PCF is symmetric under permutations of the indices
($\ell_{1},\ell_{2},\ell_{3})$ and
($\mathrm{b}_{1},\mathrm{b}_{2},\mathrm{b}_{3}$). To recover the full 4PCF one
can just permute indices and preserve corresponding values. We also define a
small function to reduce clutter with if statements below.

$\textrm{define }S(m)=\begin{cases}1/2&\text{if }m=0\\\
1&\text{else}\end{cases}$

$\delta\leftarrow$ read in density cube

for  _$\ell_{1}$ = 0 to $\ell_{\rm max}$_ do

for  _$\ell_{2}$ = 0 to $\ell_{\rm max}$_ do

for  _$\ell_{3}=|\ell_{1}-\ell_{2}|$ to $\rm min(\ell_{1}+\ell_{2},\ell_{\rm
max})$_ do

if  _$\ell_{1}+\ell_{2}+\ell_{3}\mathrm{\;is\;odd}$_ then

continue

end if

for  _$m_{1}=-\ell_{1}$ to $\ell_{1}$_ do

for  _$m_{2}=-\ell_{2}$ to $\ell_{2}$_ do

$m_{3}=-(m_{1}+m_{2})$

if  _$m_{3} >\ell_{3}$ or $m_{3}<0$_ then

continue

end if

$\mathcal{W}=\delta(\vec{x})C^{\Lambda}_{M}$

for  _$\mathrm{b}_{1}=1$ to $N_{\rm{bins}}$_ do

if  _$m_{1} >0$_ then

load $a_{\ell_{1}m_{1}}^{\mathrm{b}_{1}}$

else

load $(-1)^{m_{1}}a_{\ell_{1}m_{1}}^{\mathrm{b}_{1}*}$

end if

for  _$\mathrm{b}_{2}=\mathrm{b}_{1}+1$ to $N_{\rm{bins}}$_ do

if  _$m_{1} >0$_ then

load $a_{\ell_{2}m_{2}}^{\mathrm{b}_{2}}$

else

load $(-1)^{m_{2}}a_{\ell_{2}m_{2}}^{\mathrm{b}_{2}*}$

end if

for  _$\mathrm{b}_{3}=\mathrm{b}_{2}+1$ to $N_{\rm{bins}}$_ do

load $a_{\ell_{3}m_{3}}^{\mathrm{b}_{3}}$

$\zeta[\ell_{1},\ell_{2},\ell_{3},\mathrm{b}_{1},\mathrm{b}_{2},\mathrm{b}_{3}]\mathrel{+}=\textrm{sum}(2S(m_{3})\mathcal{W}\times\textrm{Re}[a_{\ell_{1}m_{1}}^{\mathrm{b}_{1}}a_{\ell_{2}m_{2}}^{\mathrm{b}_{2}}a_{\ell_{3}m_{3}}^{\mathrm{b}_{3}}])$

end for

end for

end for

end for

end for

end for

end for

end for

Algorithm 6 Combine Full 4PCF Coefficients

Result: Compute the unique final 3PCF coefficients by combining all the
$c_{m}^{\mathrm{b}}$ convolution coefficients. The final result is a 3D array
which can be indexed by $\mathrm{b}_{1},\mathrm{b}_{2}$ and $m$ where
$\mathrm{b}_{1}$ and $\mathrm{b}_{2}$ are the radial bin indices for each side
of the projected triangles, and $m$ is the index capturing the angular
multipole.

Notes: We only calculate the unique coefficients since the structure of the
projected 3PCF is symmetric under permutations of the indices $\mathrm{b}_{1}$
and $\mathrm{b}_{2}$. To recover the full projected 3PCF one can just permute
indices and preserve corresponding values.

for  _$\mathrm{b}_{1}=1$ to $N_{\mathrm{bins}}$_ do

for  _$\mathrm{b}_{2}=1$ to $\mathrm{b}_{1}$_ do

for  _$m=0$ to $m_{\mathrm{max}}$_ do

$\zeta_{\mathrm{proj}}[\mathrm{b}_{1},\mathrm{b}_{2},m]$ +=
sum$\left(\delta(\vec{x})\times
c_{m}^{\mathrm{b}_{1}}c_{m}^{\mathrm{b}_{2}*}\right)$

end for

end for

end for

Algorithm 7 Combine Projected 3PCF Coefficients

Result: Compute the unique final 4PCF coefficients by combining all the
$c_{m}^{\mathrm{b}}$ convolution coefficients. The final result is a 5D array
which can be indexed by $\mathrm{b}_{1},\mathrm{b}_{2},\mathrm{b}_{3},m_{1}$
and $m_{2}$ where $\mathrm{b}_{1},\mathrm{b}_{2}$ and $\mathrm{b}_{3}$ are the
radial bin indices for each side of the projected tetrahedrons, and
$m_{1},m_{2}$ are the indices capturing the angular multipoles of each free
angle.

Notes: We only calculate the unique coefficients since the structure of the
projected 4PCF is symmetric under permutations of the indices
($\mathrm{b}_{1},\mathrm{b}_{2},\mathrm{b}_{3}$) and ($m_{1},m_{2}$). To
recover the full projected 4PCF one can just permute indices and preserve the
corresponding values. We also do not index by $m_{3}$ due to the constraint
$m_{3}=-(m_{1}+m_{2})$.

for  _$\mathrm{b}_{1}=1$ to $N_{\mathrm{bins}}$_ do

for  _$\mathrm{b}_{2}=1$ to $\mathrm{b}_{1}$_ do

for  _$\mathrm{b}_{3}=1$ to $\mathrm{b}_{2}$_ do

for  _$m_{1}=-m_{\mathrm{max}}$ to $m_{\mathrm{max}}$_ do

if  _$m_{1} <0$_ then

$c_{m_{1}}^{\mathrm{b}_{1}}=(-1)^{m_{1}}c_{-m_{1}}^{\mathrm{b}_{1}}$

end if

for  _$m_{2}=-m_{\mathrm{max}}$ to $m_{\mathrm{max}}$_ do

if  _$m_{2} <0$_ then

$c_{m_{2}}^{\mathrm{b}_{2}}=(-1)^{m_{2}}c_{-m_{2}}^{\mathrm{b}_{2}}$

end if

$m_{3}=-(m_{1}+m_{2})$

if  _ $\mathrm{abs}(m_{3})\geq m_{\mathrm{max}}$_ then

continue

else

if  _$m_{3} <0$_ then

$c_{m_{3}}^{\mathrm{b}_{3}}=(-1)^{m_{3}}c_{-m_{3}}^{\mathrm{b}_{3}}$

end if

$\zeta_{\mathrm{proj}}[\mathrm{b}_{1},\mathrm{b}_{2},\mathrm{b}_{3},m_{1},m_{2}]$
+= sum$\left(\delta(\vec{x})\times
c_{m_{1}}^{\mathrm{b}_{1}}c_{m_{2}}^{\mathrm{b}_{2}}c_{m_{3}}^{\mathrm{b}_{3}}\right)$

end if

end for

end for

end for

end for

end for

Algorithm 8 Combine Projected 4PCF Coefficients