Algorithm to Produce a Density Field with Given Two, Three, and Four-Point Correlation Functions

Here we show how to produce a 3D density field with a given set of higher-order correlation functions. Our algorithm enables producing any desired two-point, three-point, and four-point functions, including odd-parity for the latter. We note that this algorithm produces the desired correlations about a set of ``primary'' points, matched to how the spherical-harmonic-based algorithms ENCORE and CADENZA measure them. These ``primary points'' must be used as those around which the correlation functions are measured. We also generalize the algorithm to i) $N$-point correlations with $N>4$, ii) dimensions other than 3, and iii) beyond scalar quantities. This algorithm should find use in verifying analysis pipelines for higher-order statistics in upcoming galaxy redshift surveys such as DESI, Euclid, Roman, and Spherex, as well as intensity mapping. In particular it may be helpful in searches for parity violation in the 4PCF of these samples, for which producing initial conditions for N-body simulations is both costly and highly model-dependent at present, and so alternative methods such as that developed here are desirable