ENCORE: Estimating Galaxy $N$-point Correlation Functions in $\mathcal{O}(N_{\rm g}^2)$ Time

We present a new algorithm for efficiently computing the $N$-point correlation functions (NPCFs) of a 3D density field for arbitrary $N$. This can be applied both to a discrete galaxy survey and a continuous field. By expanding the statistics in a separable basis of isotropic functions built from spherical harmonics, the NPCFs can be estimated by counting pairs of particles in space, leading to an algorithm with complexity $\mathcal{O}(N_{\rm g}^2)$ for $N_{\rm g}$ particles, or $\mathcal{O}(N_{\rm FFT}\log N_{\rm FFT})$ when using a Fast Fourier Transform with $N_{\rm FFT}$ grid-points. In practice, the rate-limiting step for $N>3$ will often be the summation of the histogrammed spherical harmonic coefficients, particularly if the number of bins is large. In this case, the algorithm scales linearly with $N_{\rm g}$. The approach is implemented in the \textsc{encore} code, which can compute the 3PCF, 4PCF, 5PCF, and 6PCF of a BOSS-like galaxy survey in $\sim$ 100 CPU-hours, including the corrections necessary for non-uniform survey geometries. We discuss the implementation in depth, along with its GPU acceleration, and provide practical demonstration on realistic galaxy catalogs. Our approach can be straightforwardly applied to current and future datasets to unlock the potential of constraining cosmology from the higher-point functions.