Measurement of Parity-Odd Modes in the Large-Scale 4-Point Correlation Function of SDSS BOSS DR12 CMASS and LOWZ Galaxies

A tetrahedron is the simplest shape that cannot be rotated into its mirror image in 3D. The 4-Point Correlation Function (4PCF), which quantifies excess clustering of quartets of galaxies over random, is the lowest-order statistic sensitive to parity violation. Each galaxy defines one vertex of the tetrahedron. Parity-odd modes of the 4PCF probe an imbalance between tetrahedra and their mirror images. We measure these modes from the largest currently available spectroscopic samples, the 280,067 Luminous Red Galaxies (LRGs) of the Baryon Oscillation Spectroscopic Survey (BOSS) DR12 LOWZ ($\bar{z} = 0.32$) and the 803,112 LRGS of BOSS DR12 CMASS ($\bar{z} = 0.57$). In LOWZ we find $3.1\sigma$ evidence for a non-zero parity-odd 4PCF, and in CMASS we detect a parity-odd 4PCF at $7.1\sigma$. Gravitational evolution alone does not produce this effect; parity-breaking in LSS, if cosmological in origin, must stem from the epoch of inflation. We have explored many sources of systematic error and found none that can produce a spurious parity-odd \textit{signal} sufficient to explain our result. Underestimation of the \textit{noise} could also lead to a spurious detection. Our reported significances presume that the mock catalogs used to calculate the covariance sufficiently capture the covariance of the true data. We have performed numerous tests to explore this issue. The odd-parity 4PCF opens a new avenue for probing new forces during the epoch of inflation with 3D LSS; such exploration is timely given large upcoming spectroscopic samples such as DESI and Euclid.


INTRODUCTION
The laws of nature respect certain symmetries; the physical processes governed by them are invariant under the corresponding transformations.Parity transformation (P), which reverses the sign of each coordinate axis, had been thought to be such a symmetry.Indeed, the electromagnetic and strong interactions are invariant under P.However, this symmetry is broken in the weak interaction (Lee & Yang 1956;Wu et al. 1957).Sakharov (1967) showed that the matter-antimatter asymmetry of the Universe requires that the combination CP of P and charge-conjugation (C) symmetry be broken.The currently-known CP violation is inadequate to explain the observed matter-antimatter asymmetry.Whatever additional CP violation is responsible may involve pure P violation as well.
A number of mechanisms producing parity violation at cosmological scales have been presented in the literature.For instance, one can add a Chern-Simons coupling to the standard cosmological paradigm at early or at late times.This term typically describes an interaction between a pseudo-scalar field and a spin-1 field (Sorbo 2011; Barnaby et al. 2011;Özsoy 2021) or a r2 r1 r3 < l a t e x i t s h a 1 _ b a s e 6 4 = " S M z 9 / k 4 u w / G A Q U c H u t V O K S q 1 G D Y = " > A A A B 8 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 i U A u L o I 1 l B P M B l + P Y 2 + w l S / Z 2 j 9 0 9 I Y T 8 D B s L R W z 9 N X b + G z f J F Z r 4 Y O D x 3 g w z 8 6 K U M 2 1 c 9 9 s p r K 1 v b G 4 V t 0 s 7 u 3 v 7 B + X D o 7 a W m S K 0 R S S X q h t h T T k T t G W Y 4 b S b K o q T i N N O N L q b + Z 0 n q j S T 4 t G M U x o k e C B Y z A g 2 V v J V 6 N 2 o s G b r I i x X 3 K o 7 B 1 o l X k 4 q k K M Z l r 9 6 f U m y h A p D O N b a 9 9 z U B B O s D C O c T k u 9 T N M U k x E e U N 9 S g R O q g 8 n 8 5 C k 6 s 0 o f x V L Z E g b N 1 d 8 T E 5 x o P U 4 i 2 5 l g M 9 T L 3 k z 8 z / M z E 1 8 H E y b S z F B B F o v i j C M j 0 e x / 1 G e K E s P H l m C i m L 0 V k S F W m B i b U s m G 4 C 2 / v E r a t a p 3 W a 0 / 1 C u N 2 z y O I p z A K Z y D B 1 f Q g H t o Q g s I S H i G V 3 h z j P P i v D s f i 9 a C k 8 8 c w x 8 4 n z 8 G s 5 B y < / l a t e x i t > r 1 < r 2 < r 3 r2 r1 r3 Figure 1.Parity transformation applied to a tetrahedron formed by a quartet of galaxies.Each vertex represents a galaxy.Choosing one galaxy (red dot) as our primary, the quartet is defined by the three vectors to the remaining vertices, r 1 , r 2 , and r 3 .For a quartet, the subscripts are fixed by requiring that r 1 ≤ r 2 ≤ r 3 .When viewing the tetrahedron from the primary (red) looking down along each vector r i , the direction in which one reads going from smallest to largest side (r 1 to r 3 ) defines a handedness, either clockwise or counterclockwise.
Here, the tetrahedron on the left, as viewed from the primary, is clockwise.Parity transformation in 3D is a reflection about a plane and then a 180 • rotation about the vector perpendicular to that plane, and converts the clockwise tetrahedron at left to the counterclockwise one at right.When one averages over rotations, as in this work, only the mirroring matters.
spin-2 field (Jackiw & Pi 2003;Alexander & Yunes 2009;Soda et al. 2011;Dyda et al. 2012).The pseudo-scalar field can be axion-like and if present at late times, can play the role of dark matter or dark energy.In this case, the Chern-Simons coupling can rotate the polarizations of initially linearly-polarized CMB photons. 1 In contrast, if the axion-like field plays the role of the inflaton, the coupling can give rise to non-vanishing parity-odd polyspectra of the primordial curvature perturbations (Bartolo et al. 2015;Shiraishi 2016).Since the curvature perturbations seed the subsequent formation of large-scale structure, the primordial parity-odd polyspectra would produce the same in the late-time distribution of galaxies.Recently, Cahn et al. (2023) made the novel proposal of using the galaxy 4-Point Correlation Function (4PCF) to probe parity violation in 3D large-scale structure.Four galaxies can be taken as the vertices of a tetrahedron, the lowest-order 3D shape that cannot be rotated into its mirror image, rendering the 4PCF sensitive to parity violation.An illustration of a galaxy quartet and how we define parity on it, is shown in Fig. 1.Following Cahn et al. (2023), we expand the 4PCF in the isotropic (i.e.rotation-averaged) basis functions of Cahn & Slepian (2020).In the standard inflationary paradigm (Albrecht & Steinhardt 1982;Linde 1982Linde , 1983)), we would not expect a parity-odd 4PCF.The initial density fluctuations are a Gaussian Random Field (Bardeen 1980;Starobinsky 1982), which then evolves under gravity and forms galaxies at late times.Gravity, and even the baryonic physics of galaxy formation, is parity-conserving.Hence, the detection of a parity-odd 4PCF of cosmological origin would be evidence that parity violation was present before the known forces dominated the evolution of the matter distribution.
We present here a measurement of the parity-odd modes of the 4PCF measured using the Baryon Oscillation Spectroscopic Survey (BOSS; BOSS collaboration et al. 2017) of Sloan Digital Sky Survey (SDSS)-III (Dawson et al. 2013;Eisenstein et al. 2011).Philcox et al. 2021b presented the parity-even 4PCF measurement on the same dataset and found an 8.1σ detection of a non-Gaussian 4PCF (expected in the standard picture of gravitationally evolved structure formation, e.g.Bernardeau et al. 2002).A progenitor of the algorithm and approach here used was employed to measure the 3-Point Correlation Function (3PCF) of BOSS DR12 CMASS (Slepian et al. 2017a,b,c, Sugiyama et al. 2019, 2021), and extended to use Fourier transforms in Slepian & Eisenstein (2015c); Portillo et al. (2018) and to the anisotropic 3PCF in Slepian & Eisenstein (2018); Friesen et al. (2017); Garcia & Slepian (2020).The 4PCF has been measured before in just a few works (in 2D;Fry & Peebles 1978; averaged over internal angles (Sabiu et al. 2019), as well as in Fourier space and with a degree of compression (integrated trispectrum; Gualdi & Verde 2022), but never separated into parity-odd modes; the history of N-Point Correlation Functions (NPCFs) is reviewed in Peebles (2001).
Given that as yet no detection of parity-odd physics in large-scale structure has been made and that a number of proposed theoretical models can produce it, in this work we pursue a model-independent analysis.Hence we lack an a priori expectation for the shape of the signal.While using such a model could strengthen any detection significance by correlating data in different modes, it would inevitably tie our detection significance to a particular model, which we wish to avoid.In contrast to typical analyses e.g. of the 2PCF for Baryon Acoustic Oscillations (BAO) or the 3PCF for BAO or galaxy biasing, in this work we cannot identify systematic errors simply by observing departures from an expected template for the cosmological signal.We must thus pay especial attention to systematics.We make extensive use of both mock catalogs and analytics to assess whether any systematics can produce spurious parity-odd 4PCF modes.
Our fiducial cosmology here matches that adopted by BOSS (BOSS collaboration et al. 2017).In particular, we take a geometrically flat ΛCDM model with redshift-zero matter density (in units of the critical density) Ωm = 0.31, baryon density (in units of the critical density) Ω b = 0.048, Hubble constant h ≡ H0/(100 km/s/Mpc) = 0.676, root-mean-square density fluctuations (within 8 h −1 Mpc spheres) of σ8 = 0.8, scalar spectral tilt ns = 0.96 and sum of the neutrino masses i mν,i = 0.06 eV.
The present work is organized as follows.§3 outlines the multiple analyses of varying complexity used to increase confidence in our measurements.§2 reviews the basis used to decompose the 4PCF.We also present a toy model to illustrate the relation between parity and tetrahedra.§4 describes the data, simulations, and covariance matrix.§5 then presents two different paths to obtaining a detection significance and their outcomes.We also present an analysis of cross-correlations between spatially separated regions.§6 outlines our systematics tests on mocks as well as analytic work.§7 concludes.A number of Appendices present details of the work.

METHOD FOR 4PCF MEASUREMENT
The 4PCF estimator, indicated by a hat, is ζ(r1, r2, r3) ≡ δ(s)δ(s + r1)δ(s + r2)δ(s + r3) where angle brackets denotes an ensemble average of the density fluctuations field δ(s) ≡ ρ(s)/ρ − 1, with ρ(s) the density field and ρ the average density.2Invoking ergodicity, the ensemble average may be replaced by an integral of spatial position s over the volume V .This integration results in a function that depends only on the relative separation vectors r1, r2, and r3.
We also average over joint rotations of these vectors.In practice, the density fluctuation field is computed from discrete galaxy data, appropriately weighted.Since the 3D distribution of galaxies is assumed to be isotropic on cosmological scales (ignoring redshift-space distortions), the isotropic basis (Cahn & Slepian 2020) is an efficient means of systematically extracting cosmological information.The isotropic basis functions required to measure an NPCF are given by products of (N − 1) spherical harmonics Y m (r) combined according to angular momentum addition.In particular, for the 4PCF (N = 4) we require the three-argument basis functions, which are (2) The factorizability of these functions is important to the speed-up of the 4PCF algorithm (Philcox et al. 2021a); in practice, it enables us to compute the 4PCF as a sum over the spherical harmonic coefficients a i m i of the density field about a given primary galaxy at s.Each unit vector ri is associated with one total angular momentum i, with z-component mi. 3 The key point is simply that spherical harmonics are two-index tensors, and conventionally the total angular momentum and its z-component are chosen to represent them.This point is further discussed in Cahn & Slepian (2020) around their equation 2. The weight is The 3-j symbol enforces the triangular inequality | 1 − 2| ≤ 3 ≤ 1 + 2, because the total angular momentum must be zero for an isotropic function. 4Under the parity operator, denoted P, a spherical harmonic Y m transforms as (−1) , so the three-argument isotropic functions transform as P [P 1 2 3 (r1, r2, r3)] ≡ P 1 2 3 (−r1, −r2, −r3) = (−1) 1 + 2 + 3 P 1 2 3 (r1, r2, r3) = P * 1 2 3 (r1, r2, r3).( 4) Thus the basis functions are real if 1 + 2 + 3 is even and imaginary if the sum is odd.The isotropic functions also satisfy an orthonormality relation, which follows from that of the spherical harmonics.We have dr1dr2dr3 P 1 2 3 (r1, r2, r3) P * 1 2 3 (r1, r2, r3) = δ K 1 1 δ K 2 2 δ K 3 3 . (5) The Kronecker delta δ K i i is unity when its subscripts are equal and zero otherwise.The 4PCF estimator defined in Eq. ( 1) can be expanded into the basis of isotropic functions (see Eq. 2), where the expansion coefficients depend only on the ri and are given by orthogonality as ζ 1 2 3 (r1, r2, r3) = ds V δ(s) dr1dr2dr3 δ(s + r1)δ(s + r2)δ(s + r3) P * 1 2 3 (r1, r2, r3).( 6) To avoid an over-complete basis, the radial arguments ri are ordered as r1 ≤ r2 ≤ r3, as further discussed in Cahn & Slepian (2020).
To construct a density fluctuation field from the discrete galaxy counts, and also to account for the survey geometry, we use a generalized Landy-Szalay estimator (Landy & Szalay 1993;Szapudi & Szalay 1998, see also Kerscher et al. 2000) as first outlined for the angular momentum basis in Slepian & Eisenstein (2015b) and further developed in Philcox et al. (2021a,b).It is where N ≡ (D − R) 4 and R ≡ R 4 , and these powers are shorthand for expanding by the binomial theorem and letting each D and R be evaluated at a different spatial position.D means a particle drawn from the "data" and R means a particle drawn from the "random" catalog (a spatially uniform catalog cut by the survey geometry).As outlined in Slepian & Eisenstein (2015b), we may estimate the numerator and denominator separately (i.e.compute each separately averaging over the whole survey).Doing so gives optimally weighted estimates of each in the shot-noise limit, as discussed in Slepian & Eisenstein 2015b §4, equations 24-26 and surrounding text.Multiplying Eq. ( 7) through by R, expanding each side of the resulting relation in the isotropic basis, reducing a product of two isotropic basis functions to a sum over single ones, and finally taking an inverse to solve the linear system so obtained (Slepian & Eisenstein 2015b;Philcox et al. 2021a), we find the edge-corrected 4PCF estimator as We note that there is no mixing in the radial variables; survey geometry does not change lengths.Our notation indicates a given element of M −1 .This latter is the inverse of the coupling matrix M describing how survey geometry breaks the orthogonality of our basis functions, much as Fourier modes are orthogonal only on an infinite domain.N denotes a measurement from the "data-minus-random" catalog, and R000 is the 1 = 2 = 3 = 0 expansion coefficient of R ≡ R 4 , i.e. the randoms' 4PCF.N and R are evaluated by replacing δ in Eq. ( 6) with their definitions, given below Eq. ( 7).
The coupling matrix has elements with the coefficient which depends on the product of the primary (hence the superscript "P") angular momenta. 5The matrix in curly brackets in Eq. ( 10) is a Wigner 9-j symbol.The factor C i L i i

000
(defined in Eq. 3 preceding it guarantees that i, i , and Li can be combined to make a zero total angular momentum state.It also requires that i + i + Li is even. Regarding the edge correction, we note that formally M is infinite, but we have found in practice truncating it at one angular momentum beyond that used for the physical analysis is suitable (Slepian & Eisenstein 2015b, Philcox et al. 2021a).In this work we use max = 4 for our analysis but work to = 5 on all i for the edge correction.6Further details regarding the suitability of, when performing the edge correction, truncating at an one above that used for the analysis, are in Slepian & Eisenstein (2015b).Ultimately this suitability stems from the rough tri-diagonality of the edge-correction matrix (see their §4.2. and our Fig.28).

Illustration With Toy Tetrahedra
To understand the parity-odd measurement more intuitively, we study cubic boxes of side length L box = 1000 h −1 Mpc with tetrahedra tuned to produce particular parity signals.To fill the boxes as fully as possible, yet at the same time have tetrahedra with a minimum side length of order 10 h −1 Mpc, which is similar to the situation in our BOSS dataset, we choose the three sides extending from the primary to be roughly r1 ∼ 10 h −1 Mpc, r2 ∼ 20 h −1 Mpc, and r3 ∼ 30 h −1 Mpc.We require that the minimum separation between primaries be twice the longest side of the tetrahedron (i.e. 60 h −1 Mpc) in order to minimize any overlap between tetrahedra.Finally, we have Ntets ∼ 1, 500 tetrahedra within each cubic box.Fig. 2 shows an example box and also an example tetrahedron and its partner under parity transformation.In 3D, parity transformation is equivalent to a mirror reflection across a 2D plane (which flips the sign of the coordinate axis perpendicular to that plane) plus a 180 • rotation around this latter axis.For simplicity, in Fig. 2 we depict the parity transformation simply as a mirror reflection, since our basis, being isotropic, always averages over 3D rotations and is thus insensitive to a 180 • rotation.Put another way, only the mirroring fundamentally alters the shape of a tetrahedron; the 180 • rotation only changes its orientation in absolute space. 7In particular, one can imagine the mirroring as taking one side and "pulling it through" the tetrahedron from being on one side of the plane formed by the other two sides, to being on the other side of this plane.
In practice, we must allow all four vertices of each tetrahedron a chance to be the primary (Philcox et al. 2021a).However, for this toy tetrahedron illustration we restricted to bins in side length (radial bins) such that 9 h −1 Mpc < ri < 30 h −1 Mpc for all ri so that nearly always only one of the four vertices will satisfy the radial bins required for each tetrahedron.This means that the contribution to the signal will stem only from the isotropic function evaluated on unit vectors extending from a single primary; this renders it easier to understand the measured signal in this toy model.
We set r1 to be along x, r2 to be along ŷ, and r3 to be along ẑ.The tetrahedron is therefore clockwise at the only primary allowed in this toy model; our convention is presented in Fig. 1.A parity transformation simply means that we interchange the sides, so that r2 aligns with the x-axis and r1 aligns with the y-axis, while r3 is unchanged.Again, characterization as "clockwise" or "counterclockwise" depends on at which vertex one sits, as discussed in Fig. 1, but here is unambiguous.By restricting the radial bins we use, we force only one galaxy to be the primary; the other choices of primary will not lead to sides extending from them that fall within our chosen bins.Thus, we can ensure that the tetrahedra for this test are perfectly counterclockwise as viewed from the single primary allowed.
We produce toy boxes in three configurations.The first is filled with only counterclockwise tetrahedra (enforced by the radial bin restriction); the second has only clockwise tetrahedra, and the third has an equal mix.To render each toy box somewhat more realistic, we randomly rotate each vertex by an angle θ ∈ [−180 • , 180 • ] around each of the three Cartesian coordinates.We also add random numbers ∆r1 ∈ [0, 1], ∆r2 ∈ [−2, 2], and ∆r3 ∈ [−1, 0] (in h −1 Mpc) to, respectively, r1, r2, and r3.We choose these ranges such that we always have r1 < r2 < r3.Hence the parity of a given tetrahedron will not be flipped by these additions.
Fig. 3 shows the 4PCF of these illustrative boxes.We may approximate each tetrahedron as a sphere about the primary of radius roughly 20h −1 Mpc, with volume Vtet, to estimate the expected 4PCF amplitude.We define n as the local number density due to a given tetrahedron, n = 1/Vtet, and n as the average number density in the box, n = Ntets/L 3 box .We find (Cahn & Slepian 2020).As expected, the counterclockwise-only box has a positive projection onto this function, while the clockwise-only box has a negative projection.The mixed box is consistent with zero projection onto P111 on average.We can analytically predict the ratios among the 4PCF coefficients for different channels 1, 2, 3. We compare the mean ratios of the measured 4PCF coefficients for several combinations to these predictions and find good agreement.8This also serves as an additional test of our code (the code is further discussed in §4.1).

Internal Cancellation
For a given tetrahedron, in practice (but not in our illustrative boxes above), each of the four vertices gets a chance to serve as the primary about which the isotropic basis function expansion is computed.Some of these vertices will be "clockwise", and some "counterclockwise".Hence if co-added into the same channel and triple-bin there will be "internal cancellation" and consequent reduction of any parity-odd signal.However, if the radial bins are made fine enough then each vertex, in virtue of the presumably unique lengths of the sides extending from it, will be accumulated to a different triple-bin.9Thus finer binning can reduce the internal cancellation and increase the signal.This is much the same as in a configuration-space BAO search, Each of them has a unique primary (black) from which the three sides with respective lengths r 1 ∼ 10h −1 Mpc (red), r 2 ∼ 20h −1 Mpc (yellow), and r 3 ∼30h −1 Mpc (blue) extend.The larger panel on the left is a zoom-in on the full box to display the tetrahedra more clearly.Right: A sketch of a "clockwise" tetrahedron and its "counterclockwise" mirror image.The primary is in red.Our convention on clockwise and counterclockwise is detailed in Fig. 1.On the left, the red point is closer to us than all the others.Thus the tetrahedron on the left is clockwise, as looking down from the primary, we go clockwise as we move from the smallest side to the largest.On the right, the primary (in red) is behind the other galaxies, so looking down from it towards them will reverse the handedness.Thus the rightmost tetrahedron is counterclockwise as viewed from the primary.
. By construction the scalar triple product r1 • (r 2 × r3 ) is unity for all three of our illustrative boxes.The amplitudes on the righthand side are not strictly zero because there are still a small number of tetrahedra formed by connecting secondaries around one primary with secondaries around another.Additionally, the amplitude fluctuations divide into two envelopes at bin index 36.This bin index corresponds to raising the index for r 1 (the smallest side and the slowest-varying one).
where fine enough bins must be chosen that the BAO feature in the 2PCF is not averaged out by all being added into a single bin.

GUIDING PRINCIPLES FOR THE ANALYSIS
(This section used to be directly after the introduction; now it has been moved here.)To isolate the potentially parity-violating component of the 4PCF, we expand the correlation function in two distinct sets of isotropic basis functions, one that is parity-even and one that is parity-odd.These are constructed from products of three spherical harmonics with angular-momentum indices 1, 2, 3. Isotropy requires that the i satisfy the triangular inequality.If the sum of the i is even, the product is parity-even and if the sum is odd, the product is parity-odd.In this analysis, only the parity-odd elements are used.Each basis element is a function of three radial distances, the length of the sides from a chosen vertex among the four defining a tetrahedron.In practice, the radial distances are binned.
Ideally, one would like to capture as much information as possible, using many narrow radial bins, and working to some high value of the i.This would increase the difficulty of evaluating all the independent amplitudes, but in fact, the technique of Slepian & Eisenstein (2015b); Philcox et al. (2021a) makes this manageable, as we discuss in §2.The real challenge is in determining the covariance matrix.This is especially critical when looking for parity violation, where establishing a statistically significant non-zero signal is the crux and hence understanding the inevitable statistical fluctuations is essential.
A fairly modest choice of separating the radial variable into ten bins and including i only up to max = 4 results in 23 × 120 = 2, 760 independent amplitudes; we also consider eighteen radial bins, which produces 816 × 120 = 18, 768 4PCF amplitudes.Determining the covariance matrix is thus a formidable challenge.In order to invert a sampling covariance matrix, as required for calculating χ 2 , we need at least many mocks as the dimensionality of the data vector, which is in excess of the 2,000 available to us for BOSS.
We have chosen three ways to obtain the covariance matrix.First, the NPCF covariance matrix can be calculated analytically under the assumption of a Gaussian Random Field (GRF) as shown in Hou et al. (2021a).The GRF does not on average have any parity-violation at the signal level, but it can still have non-zero fluctuations in parity-odd modes.Therefore the GRF is still the leading-order contribution to the covariance of the parity-odd modes.This is simply the statement that a signal may be zero but its root-mean-square may not.We fit the analytic template covariance matrix by varying the number density and volume with respect to the covariance matrix derived from the mocks.
With this analytic template in hand, we may (i) directly compute the χ 2 of the data using the adjusted analytic covariance matrix.An alternative (ii) is to compress the data vector to reduce its dimensionality.The eigenvectors of the analytical covariance matrix with the smallest eigenvalues represent the linear combinations of basis functions that have the smallest statistical uncertainties.We then expand the measured 4PCF using just the Neig best expansion functions.We may also determine the Neig × Neig covariance matrix directly from the mocks.Since Neig is much less than the number of mocks this covariance matrix is invertible.Finally, we may (iii) use the empirical covariance matrix from the mocks directly, with no involvement of the analytic template at all, by considering many fewer channels than in (i) and (ii) (we lower max to be max = 2).A substantial reduction of the number of channels is required to enable inverting the empirical covariance matrix employed in this approach, and so we lose statistical power.We thus treat (iii) as a test rather than as giving us the main result of our analysis.The reliability of all three approaches above can be assessed using the mocks themselves by verifying that their χ 2 (or T 2 ) values match the expected distribution.

DATASET AND COVARIANCE
We use the final galaxy catalog of the Baryon Oscillation Spectroscopic Survey (BOSS), from the twelfth data release (DR12) of the Sloan Digital Sky Survey-III (SDSS-III).The catalog is split into the North Galactic Cap (NGC) and the South Galactic Cap (SGC).The catalog contains two samples, CMASS and LOWZ, which were selected via the SDSS multicolor photometry and cover a redshift range of 0.15 < z < 0.7.CMASS and LOWZ use similar target selection algorithms (Eisenstein et al. 2001;Cannon et al. 2006).The target selection algorithm provides samples that are mainly composed of Luminous Red Galaxies (LRGs).For CMASS the selection algorithm is further tuned to select massive objects uniformly in redshift (Reid et al. 2016), which results in an approximately mass-limited sample down to a stellar mass M ∼ 10 11.3 h −1 M (Maraston et al. 2013;Leauthaud et al. 2016;Saito et al. 2016;Bundy et al. 2017).The majority of CMASS is LRGs (∼ 74 percent), while the rest is late-type spirals (Masters et al. 2011).LOWZ consists primarily of LRGs (Parejko et al. 2013).Despite the difference in target selection, the LRGs in the two samples have similar stellar mass distribution (Maraston et al. 2013).We apply a redshift cut of 0.43 < z < 0.7 to CMASS, which results in a redshift tail from LOWZ at redshift 0.43 < z < 0.5.This tail (∼ 15% of the entire "CMASS" sample) slightly raises the purity of the CMASS sample by adding more LRGs.To ensure that the LOWZ sample as used here is independent of the CMASS one, we apply a redshift cut of 0.2 < z < 0.4 to the former.This cut produces a separation of about 70h −1 Mpc between the lower edge of CMASS and the upper edge of LOWZ.Fig. 4 shows the number density for each sample as a function of redshift.Finally, we note that the early LOWZ target selection was not uniform due to the use of different iterations of the galaxy-star separation algorithm (Reid et al. 2016).Therefore we do not include those early chunks in this analysis.

4PCF from Data and Mocks
For our main analyses, we considered tetrahedra with side lengths r1, r2 and r3 from the primary ranging from 20 h −1 Mpc to 160 h −1 Mpc inclusive, split into ten linearly-spaced bins, and also from 20 h −1 Mpc to 164 h −1 Mpc inclusive, split into eighteen  The LOWZ (0.2 < z < 0.4) North Galactic Cap (NGC) is brown and South Galactic Cap (SGC) is orange.For CMASS (0.43 < z < 0.7), the NGC is in purple and the SGC is in lavender.We intentionally do not allow redshift overlap between the two samples, and the redshift gap ∆z = 0.03 corresponds to a comoving radial separation of 73 h −1 Mpc.This separation means that the samples are fairly independent; the 2PCF ξ between a point in LOWZ and CMASS would be of order 1% at this scale.The covariance between two "worst-case" tetrahedra, one in each sample (where each is very close to the respective edge), is of order ξ 4 .Few tetrahedra are near enough on either edge to be significantly correlated with those in the other slice.The plot shows that LOWZ has both a more uniform selection function and a somewhat higher average number density than CMASS.It also lacks the strongly decaying tail with an increasing redshift that CMASS displays.These points are important when assessing the possible impact of systematics on each sample and when addressing any differences between the detection significances in the two samples.
linearly-spaced bins.We also explored a coarser binning (six bins), presented as a test in Appendix E. As a result, the sides of the tetrahedron that do not include the primary can range from zero to 320 h −1 Mpc.We expand the parity-odd 4PCF in the 23 angular channels with i ≤ 4 given in Cartesian form in Appendix A. For the edge corrections, we include all functions with i ≤ 5, as further discussed in §2.We also compute the even-parity modes (for which we do not reproduce the basis functions here), as they are needed within the edge-correction, despite that our actual analysis focuses solely on the parity-odd ones.
To each galaxy, we apply a total weight w given by w = w fkp wsys(wnoz + wcp − 1), (11) with the systematic weight wsys (a combined weight for stellar density and seeing), the redshift failure weight wnoz, and the fiber collision weight wcp (subscript "cp" for "close pairs").The FKP weight (Feldman et al. 1994) is is the weighted number density at the given redshift.We use the public random catalog provided by BOSS10 , which is fifty times the size of the data catalog.We split it into 32 chunks such that each chunk is about 1.5 to 2 times the size of the data; the reason is discussed in Slepian & Eisenstein (2015b) and Philcox et al. (2021b) though the detailed splitting does not affect the 4PCF algorithm speed notably.Each chunk is first randomly shuffled and then normalized so that its weighted sum matches the sum of the completeness weights wc of the data.We use 2,000 public MultiDark Patchy lightcone mocks (Patchy hereafter; Kitaura et al. 2016;Rodríguez-Torres et al. 2017).These mocks use second-order Lagrangian Perturbation Theory (2LPT) plus a spherical collapse prescription calibrated on full N-body simulations, and the mocks are additionally calibrated to match the 2-point and some 3-point statistics of the observed BOSS CMASS and LOWZ samples.The Patchy mocks include realistic survey geometry as well as many observational effects, further detailed in Rodríguez-Torres et al. (2017).
The 4PCF measurement for the BOSS data catalog and the 32 (D − R) catalogs alone takes 0.33 GPU-hours for the measurement where we use ten radial bins, and 0.65 GPU-hours for the measurement where we use eighteen radial bins; both timings are using 69 NVIDIA A100 GPUs (further discussed below) and include edge-correction, which is done on the CPU linked to the GPU and takes negligible extra time.
All calculations were done with the GPU-accelerated NPCF code CADENZA (Slepian et al. 2021), built on the CPUbased NPCF code ENCORE (Philcox et al. 2021a).CADENZA uses the representation of the basis functions in spherical harmonics, for reasons outlined in §2.Tests of ENCORE are described in (Philcox et al. 2021a), and CADENZA was verified to obtain exactly the same outputs on a fiducial data set used for testing.As noted already, we used 69 NVIDIA A100 GPUs simultaneously on the HiPerGator cluster at the University of Florida.Overall, all of the computations performed for this work would have taken roughly 7.5 million CPU hours, or about 20 months of continuously running on a cluster of 500 CPU cores.
Parity-Odd 4PCF Detection 9 Radial Bin Index  Here for legibility we focus on the ten-bin results; the eighteen-bin results are in Appendix G, as are the remaining channels for the ten-bin analysis.We have mapped the three radial bins to a single index, with r 3 varying the fastest and r 1 the slowest.This is done in all similar plots that follow.
Lower: The mean of 2, 000 Patchy mocks (solid curves) for both NGC (red) and SGC (blue).The shaded region is the rms expected for a single mock; this set of panels is intended to display the region around zero that the 4PCF of a dataset with no true parity-odd signal might inhabit.
In practice, since CADENZA is about 140× faster than ENCORE for the 4PCF, we required about 54,000 GPU-hours, or about month if one ran continuously on 69 GPUs.

Analytic Covariance
The covariance matrix for the 4PCF may be computed analytically if we assume that the density is a Gaussian Random Field (Hou et al. 2021a).It is expressed in terms of f -functions with j (kr) the spherical Bessel function of order and P (k) the galaxy power spectrum.The power spectrum is taken to be isotropic; the effect of RSD (after averaging over the line of sight) is simply to boost the amplitude, which is embedded in the input power spectrum and also partially absorbed in our later fitting for the best amplitude.Fast approaches do exist to evaluate these f -integrals (e.g.Slepian et al. 2019) but here we use direct computation.
The analytic covariance matrix couples each vertex of an (unprimed) tetrahedron to a vertex in another (primed) tetrahedron.Since in our approach, the primary vertex is distinguished from the other three, in the covariance expression there are two pieces: in one the two primary vertices, unprimed and primed, are paired with each other; in the other, they are not.The  1 2 3 (r1, r2, r3) First and second rows: measurement of the parity-odd 4PCF using ten radial bins for the BOSS LOWZ data in the NGC (brown) and SGC (blue).We show the same six channels as in Fig. 5, and the error bars are again the rms over the mocks.Again for legibility, we focus on the ten-bin results; the eighteen-bin results are in Appendix G, as are the remaining channels for the ten-bin analysis.Third and fourth rows: The mean of 2, 000 Patchy mocks (solid curves) for both NGC (brown) and SGC (blue).The shaded region is the rms expected for a single mock.
final result will be a sum given by products involving Wigner-nj symbols and the f -functions of Eq. ( 12).Below we simply reproduce the final result derived in Hou et al. (2021a), for brevity here we define Λ ≡ { 1, 2, 3} and Λ ≡ { 1 , 2 , 3 }.
In Case I, the sum over G includes all six permutations, while in Case II the sum over G includes only cyclic permutations.The sum over H includes all six permutations.EG is the Levi-Civita symbol, giving a negative sign if the permutation G is odd.Σ(Λ ) ≡ 1 + 2 + 3 .The coefficients C 1 2 3 m 1 m 2 m 3 and D P 1 2 3 are given in Eq. ( 3) and Eq. ( 10), respectively.For the parity-even modes of the 4PCF, the phases in Eqs. ( 14) and ( 15) are identically unity because the sum of the three angular momenta is always even.However, these phases do matter for the parity-odd 4PCF.
We also note that, apart from the leading Gaussian contribution, the covariance should also contain non-Gaussian contributions as well as the coupling of the redshift-space power spectrum multiple moments to different angular channels of the covariance.Any combination of multipoles that after multiplying out is rotation-invariant (i.e. has total angular momentum zero) can contribute.An expression for this RSD-added analytic covariance is in Hou et al. (2021a) Appendix E., where it is also shown that these redshift-space induced coupling can also increase the off-diagonals of the covariance.Hou et al. (2021a) showed that the rms of the off-diagonals after comparison of the analytic covariance to the mock-based covariance is exactly the same in real-space and redshift-space, both for lognormal mocks and for quijote simulations (cf.figure 6, 7 and figure 10, 11).Therefore, off-diagonal enhancement by redshift-space-induced coupling is in principle not significant enough to greatly affect our analysis.Furthermore, some of the redshift-space-induced effects can be absorbed in the overall normalization, as much of their effect after isotropic averaging is to rescale the overall power spectrum.Nevertheless, we caution that missing these other terms could in principle reduce the accuracy of our analytic covariance relative to that from the mocks.
Fitting the volume and number density to match the mocks' covariance can be done without inverting this latter, important to avoid as that inverse will be biased unless the number of degrees of freedom is much smaller than the number of mocks.For the fitting, we maximize a likelihood based on the Kullback-Leibler (KL) divergence (Kullback & Leibler 1951).One can show that with the likelihood Eq. ( 55) of Hou et al. (2021a), the optimal volume at a given number density is uniquely determined, so the optimization is only a 1D problem, making it efficient.There is no guarantee that the best-fitting number density and volume for the covariance of the even-parity 4PCF (e.g. as found in Philcox et al. 2021b) and for the covariance of the parity-odd 4PCF should be the same.We summarize the values obtained for each sample and each radial binning in Table 1.For comparison, we also give the effective volume estimated from BOSS using equation ( 50) of Reid et al. (2016); the survey area used in that equation is taken from Table 2 of Reid et al. (2016).As discussed above, the template covariance employs several approximations, including the assumption that the density field is purely Gaussian Random, and the omission of higher power spectrum multipole moments (w/rt line of sight).Such simplifying assumptions may cause an underestimate of the true covariance.To address this issue, an overall pre-factor is introduced and is interpreted as the effective volume V eff , where this fitted effective volume could be lower than the one reported in Reid et al. (2016).However, we want to point out that this volume-like prefactor is only a highly approximated estimation.To fully account for the effective volume requires including geometry effects in the analytic covariance template.
In Fig. 7 we compare for CMASS NGC the parity-odd 4PCF analytic covariance matrix to that from Patchy for ten radial bins.We have mapped the channel and triple-bin to a 1D index, so the covariance may be plotted as a 2D matrix (one 1D index for the unprimed channels and triple-bins, another for the primed).In the left panel, we show a direct comparison of the correlation matrix, which is the covariance matrix normalized by its diagonal.The similarity between the analytic correlation matrix (upper triangle) and the mock correlation matrix (lower triangle) shows that the analytic approach can well capture the covariance's features.In the right panel, we compare the diagonal elements of the covariance matrices and also show the ratio between the mock and analytic covariance diagonals.Again this uses our mapping of the angular momenta and radial bins into a 1D index.Despite the similarity between the two diagonals in their behaviour with increasing index, there is non-negligible variation.On average, the mean of the ratio between the mock and analytic covariance diagonal elements is roughly 1.05 for CMASS NGC.Fig. 8 shows the same comparisons for the LOWZ NGC.Again, we see that the analytic covariance can well describe that of the mocks.The mean of the ratio is 1.01 for the LOWZ NGC.
To further quantify the similarity between the analytic and mock covariance, we define a half-inverse test matrix as where analytic is the square root of the inverse of the analytic covariance matrix, and C patchy is the covariance from the mocks.We use half-inverses so that the test is symmetric.Were the two covariance matrices identical, Φ would be the zero matrix (we note that the identity matrix is subtracted).One can show that the optimal volume at a given number density (as discussed above) will imply that Φ is traceless.
Fig. 9 shows the lower half of the matrix Φ for the half-inverse test of CMASS NGC (left) and LOWZ NGC (right), with the standard deviation shown on the upper triangle.The standard deviation for all elements is σ all = 0.06 for both CMASS NGC and LOWZ NGC, which is expected as it scales as 1/ √ N mocks (we used N mocks = 2, 000 for this test).Each element of the covariance matrix follows a Wishart (1928) distribution, and the variance of the diagonal elements for this distribution is two times that of the off-diagonal elements.We do find that the standard deviation σ diag is approximately two times that of the off-diagonal, σ non−diag .However, we also find residuals in the diagonal elements in Φ for both the left and right panels, and larger residuals for LOWZ NGC than for CMASS NGC.These indicate that fitting the volume and number density to best match the mock covariance is imperfect.As mentioned in §3, one of our analysis approaches ("compressed") will only require a smooth estimate of the covariance and turns out to be insensitive to its amplitude.Furthermore, even when employing the   13), and the volume and number density here used are given in Table 1.Left: analytic correlation matrix (upper triangle) and correlation matrix from the Patchy mocks (lower triangle).The sub-blocks represent the covariance between a pair of channels; each is 120 × 120 as for ten radial bins, there are 120 radial bin combinations in each channel.These sub-blocks show how the covariance changes as the side lengths of each tetrahedron vary at fixed channel pairs.The similarity between the upper and lower triangles shows that the analytic covariance captures the cross-covariance between different radial bins as well as between different channels.Right: diagonals of the Patchy covariance matrix (solid black) and the analytic covariance (dashed red) for CMASS NGC.The average ratio is 1.05.analytical covariance matrix directly, we may compare the χ 2 from data only to the distribution of χ 2 values for the mocks, both computed using the same covariance matrix.This comparison should minimize any bias in the detection significance computed using the analytic covariance.

Mock Covariance Matrix
In this work, we assume that the Patchy mocks match well the covariance structure of the observed data.This is a reasonable assumption as the Patchy mocks are calibrated to match the 2-point and some 3-point statistics measured from BOSS.Here we outline in detail aspects of Patchy that could cause a mismatch between the mock covariance and the data's true covariance.In the end, we do not find any strong reason to believe that such a mismatch is present, but that is after performing a number of tests motivated by the below points.
(i) Patchy uses approximate methods to simulate nonlinear structure formation, and underestimates small-scale power (s 20 h −1 Mpc), especially for galaxies with high stellar mass (Kitaura et al. 2016; see their Fig. 4).Lacking power at smaller scales could have two impacts on the covariance.First, it reduces the "squeezed" tetrahedra (where the side lengths from the primary are close to each other in size and hence galaxies can become very close to each other as well as the angles at the primary vertex become small).This could produce an underestimate of the covariance on small scales.Second, this lack of small-scale power also may reduce the derived covariance at all scales, as the covariance comes from averaging the product of unprimed and primed tetrahedra over all possible separations between them, including small separations wherein points on the primed one may be close to those on the unprimed one.Indeed, this is the origin of s in the integral Eq. ( 12).
To assess the impact of the above, we do an analysis where we force each side to be different in length from the others by at least one full bin width; this at least addresses the first point in (i).We still find a significant detection of parity-odd 4PCF in all samples.The second point is more challenging to deal with, although we did not have a direct test against it; its impact can be inferred from the parity-even connected 4PCF measurement (which tracks only the contribution due to non-linear evolution).The level of agreement between Patchy and BOSS CMASS suggests that Patchy actually does well-capture the non-linearity on the scales relevant to our analysis.
(ii) Patchy may not reproduce all systematics that may be present in data.We computed the 2PCF and found that the LOWZ SGC has an amplitude that is higher on average than the mean of the Patchy mocks by 1σ.However, it is quite consistent with the spread we see in the measured 2PCF for a set of 1,000 mocks.In addition, as discussed in Ross et al. (2016) and Kitaura et al. (2016) there are deviations between Patchy and the BOSS 2PCF at scales greater than roughly 100 h −1 Mpc.These 2PCF bins are highly correlated though and the actual deviation may be less severe than the visual impression.Again we found that the spread of the 2PCF of 1,000 mocks covered the BOSS data well.Nonetheless, to make our analysis robust, we also performed our parity-odd analysis with a cut on the maximum side lengths from the primary so that no side of the tetrahedron could exceed 160 h −1 Mpc.We still found evidence for a parity-odd 4PCF consistent with what one would expect given the significant reduction in the constraining power in this analysis due to the much smaller number of triple-bins it permits.
We also briefly consider that the number density of the survey is not equal to the true number density of the Universe; this would be a failure of the total integral constraint.Such a failure would produce a correction to the observed power spectrum, such that a term Pic ≥ 0 is subtracted (see equation ( 29) of Beutler et al. 2014).The Patchy mocks would not have this correction, and hence have a larger power spectrum than the observed data.Thus, the covariance estimated from them would be larger than it should be.Since this error, if present, would cause an over-estimation of the covariance (and thus a spuriously decreased detection significance), and also since given BOSS's large volume it is expected to be a small effect, we do not explore it further.

Impact of possible systematics
As will be discussed further in §6, distortions in the radial or the angular directions that are not captured in the mocks can mean that a mock-based covariance underestimates the true covariance.Here we study the sensitivity of the parity-even sector to such contamination.For a distortion of the radial selection function, we apply the contamination directly to the mocks and Table 1.Summary of the datasets and binning schemes used.We give the effective redshift z eff and the fitted effective volume V eff,fit and number density n derived by fitting our analytic covariance to that from the mocks.The fitted effective volume and number density are different for each number of radial bins because we re-fit the covariance for each radial binning scheme.In the last row, we also list the effective volume using equation ( 50) from Reid et al. (2016) 2.50 0.95 0.66 0.29 find their mean detection significance of a connected-even parity 4PCF (when computed with an undistorted covariance, as would be the case if the data had an unaccounted-for change) is shifted by ∼ 0.6σ for ten radial bins and ∼ 2.4σ for eighteen radial bins.We also explore angular contamination as described in §6.1.5;this increases the even-parity detection significance difference between data and mean of the mocks to 0.9σ for the ten bin and 1.8σ for the eighteen bin.

Self-calibration of covariance
Our detection significance relies on correctly estimating the covariance.Since the Patchy mocks have no intrinsic parityviolating mechanism, they cannot be used to check our pipeline regarding signal, but only to estimate the covariance.Yet it is possible that the mocks do not contain all the systematics present in the real data, and this lack could result in mis-estimation of the covariance.However, any unknown data-only systematic might well impact the parity-even modes as well as the parity-odd modes.Hence, we can see if the mocks and data are consistent in the parity-even sector at the signal level, where both do have a signal due to non-linear evolution.If we see that the even-parity signal from the data is consistent with that from the mocks, this suggests that there are unlikely to be unaccounted-for systematics.It also would argue that we have reasonably estimated the even-sector covariance, as an underestimate of it would show up as a severe tension between the signal in the mocks and in the data.Since our procedures for obtaining the even and odd-sector covariances are exactly the same, such a finding would in turn build confidence that our odd-sector covariance is correct.
Furthermore, as a highly conservative approach, we can ask what factor the covariance in the even-sector would need to be rescaled by to force the data to agree within e.g.3σ or 1σ with the mocks.We can then apply this rescaling to our odd-parity covariance and ask how it impacts our detection significance in the odd sector.We emphasize that this is a highly conservative check, not a correction procedure; it is perfectly possible that the even-parity data signal is e.g.3σ different from that in the mocks just by chance, without implying that there is any error that should be corrected in either sector.Thus, while we summarize the results of this idea in Table 4, one should not take these rescalings too literally.
For CMASS with ten radial bins, we found excellent agreement between the data and mock distribution both using the analytic covariance and the compressed method (also see Philcox et al. (2021b), where a slightly different sample selection was applied).Therefore, we do not present any rescaling factor for the detection significance for the parity-odd modes.
For CMASS with eighteen radial bins, we found a discrepancy in the even-parity sector between the data and the mock distribution of 4.9σ when using the analytic covariance.Rescaling the even covariance such that the data and mocks agree at the 3σ level, we then propagate the same factor to the parity-odd measurement.We also repeat this procedure to enforce the agreement at the 1σ level in the even sector.We summarize these results in Table 4.
Regarding the discrepancy found in the eighteen-bin test, we ascribe it to the possible breakdown of the GRF assumption of the analytic covariance as one goes to smaller scales, especially for the even sector, as discussed in §4.2.In particular, the parity-even covariance requires including the leading higher-order covariance and connected estimator contribution.These contributions are naturally absent when estimating the parity-odd covariance.Further, when using the compressed method with Neig = 800, the discrepancy is reduced to 1.3σ.
Table 2.The detection significances using the analytic covariance for the split sky (into NGC and SGC) and joint sky (NGC + SGC).The latter treats a given channel and radial bin combination in N and S as two separate degrees of freedom, i.e. we assume vanishing covariance between N and S. All of the quoted detection significances are from comparing the χ 2 of the data to a Gaussian fit to the histogram of the mocks' χ 2 ; see lower rightmost panels of Figs.10-13 (ten and eighteen bin results) and of Fig. E1  Here we summarize the rescaled detection significance in the parity-odd modes by forcing agreement in the parity-even sector at the levels indicated in the leftmost column.These results are for CMASS and eighteen radial bins; see §4.2.3.We explore this rescaling both fr our analytic covariance analysis and our compressed analysis with 800 eigenvalues; in the last line, no rescaling is required to attain 3σ consistency in the even sector.Generally, it is notable that an odd detection more or less remains in all cases even after rescaling.

RESULTS ON CMASS AND LOWZ DATA
Here, we first compute the significance of our observed parity-odd 4PCF using the methods briefly outlined in §3, and then explore the cross-correlation of NGC and SGC in each sample.

Statistical Significance
We quantify the statistical significance of the parity-violating amplitudes using three approaches: a direct analysis harnessing the analytic covariance calibrated against Patchy ( §5.1.1),a compressed analysis using the analytic covariance just to select a reduced basis with few enough degrees of freedom that the mocks can then provide the covariance ( §5.1.2),and finally, by lowering max enough so that we may use the direct approach but with covariance from the mocks ( §5.1.5).

Direct Analysis
With the analytic covariance, adjusted in density and volume to best fit the covariance given by the mocks, we evaluate where ζ and ζ model are respectively the various parity-violating amplitudes, ζ 1 2 3 (r1, r2, r3) and ζ model is a model for them.To investigate the null hypothesis that there is no parity violation, we set all elements of ζ model to zero.C −1 is the inverse covariance matrix.If the underlying data have Gaussian behavior with N d degrees of freedom, the resulting distribution is the The adequacy of the covariance obtained by fitting the analytic covariance to the data sets can be assessed by examining the distribution of χ 2 values from the mocks when run through our analysis, and seeing if it matches the predicted χ 2 distribution.11

Compressed Analysis
The alternative means of obtaining an invertible covariance matrix by reducing the dimensionality of the problem, the "compressed" analysis, is described in §3.This scheme was introduced by Scoccimarro et al. (1999) to quantify detection significance, and the same method was applied to the even-parity connected 4PCF analysis of Philcox et al. (2021b).We diagonalize the analytic covariance, writing C = ODO −1 where D is diagonal and O is an orthogonal matrix (we note that The columns of O then provide an eigenbasis for the analytic covariance.We would like to select the subset of these eigenvectors that will be most useful in detecting a signal.If one has a model for the expected signal, one can rank the eigenvectors by their signal-to-noise, and choose Neig of them with the best S/N.In the absence of an anticipated signal, we choose the Neig eigenvectors associated with the smallest eigenvalues.These eigenvectors have the lowest noise.The analysis is then restricted to the space spanned by these eigenvectors.With Neig N mocks , the covariance matrix for this restricted analysis can be determined entirely from the mocks, while retaining the domain with the greatest information.
While this procedure provides an invertible covariance matrix, the inverse of the mock covariance matrix is unbiased only when Neig N mocks , where here Neig is the number of degrees of freedom, i.e. plays the role of N d in §5.1.1.Consequently, the χ 2 -distribution we used in §5.1.1 is modified.To distinguish this distribution from the χ 2 distribution, we refer to the variable T 2 , defined analogously to χ 2 , as whose distribution is (Sellentin & Heavens 2016) This reduces to the χ 2 distribution for N mocks → ∞.
Since the data from CMASS and LOWZ in the NGC and SGC are statistically independent given the large physical separations of both NGC and SGC and of CMASS and LOWZ, we can simply add their values of χ 2 or T 2 , while increasing the number of degrees of freedom correspondingly.

Detection Significances from Direct & Compressed Methods
Fig. 10 displays the T 2 or χ 2 distributions for CMASS (both NGC and SGC) for the two approaches described above (direct and compressed).The detection significance changes as we vary the number of eigenvalues, but this is expected: as Neig rises, more information is added (unless the S/N is pathological), but at the same time the potential bias of the inverse covariance is increasing.This latter can be assessed by looking for when the distribution of the mocks' T 2 values begins to deviate from the T 2 distribution expected (Eq.20).Fig. 11 shows the same information but for LOWZ, which generally shows comparable detection significances compared to CMASS both at each number of eigenvalues in the compressed analysis, and in the direct analysis.

Using Finer Binning
Since two canceling contributions can occur from a single tetrahedron when they fall into the same radial bins ( §2.2), increasing the number of bins could increase the power of the analysis, at the cost of further enlarging the covariance matrix.So motivated, we now use eighteen linearly-spaced radial bins from rmin = 20 h −1 Mpc to rmax = 164 h −1 Mpc, leading to a bin width of 8 h −1 Mpc, roughly double the radial resolution of our previous analyses.We refit the volumes and number densities for the covariance matrices of CMASS NGC and SGC and LOWZ NGC and SGC as in §4.2, and these numbers are in Table 1.Other than these new inputs to the covariance matrix, all else is the same in our analysis, and we display the results in Fig. 12 50 CMASS, 10 bins, max = 4 Figure 10.T 2 and χ 2 distributions for CMASS data using ten radial bins.The first five plots are obtained using the data compression scheme of §5.1.2with different numbers of eigenvalues N eig as indicated in the plot titles.If there is no parity-breaking, we expect a T 2 distribution (Eq.20) with N eig degrees of freedom (solid black).The vertical lines in the first five plots are the T 2 values from CMASS for the NGC (red), SGC (blue), and joint NGC and SGC (black).For comparison, we show histograms of the T 2 values from the Patchy mocks for each case.We interpret the increasing deviations of the mocks' histograms from the T 2 distribution as N eig rises as evidence that the covariance becomes less well-determined for larger N eig .This is expected as for N eig = 800, there are only about three mocks per degree of freedom (we use 2,000 mocks).In the lower rightmost plot, we use the direct approach of §5.1.1,which employs the analytic covariance.Here, the expected distribution if there is no parity-breaking is a χ 2 distribution (Eq.18; solid black).We show a histogram of the mocks' χ 2 values for comparison (orange) as well as a Gaussian fit to it (dot-dashed).This histogram is wider than the predicted χ 2 distribution (which the mocks should match as they have no parity-breaking).We interpret this as an indication that the analytic covariance may not be unbiased.To be conservative, we compute a detection significance σ G by comparing the data χ 2 values to the width of the Gaussian fit to the mocks.In fact, we quote σ G in all panels above, but it generally agrees with that computed from the T 2 distribution, σ T 2 .We do not quote a detection significance computed with respect to the χ 2 distribution (as would naively be appropriate in the lower rightmost panel) because the distribution of mocks is clearly mismatched with a χ 2 distribution.
(CMASS) and Fig. 13 (LOWZ).Generally, for each sample, the detection significance rises at each Neig relative to the ten-bin analysis, and the same is true for the "direct" significances (where the analytic covariance was used).This increase is expected given our arguments about "internal cancellation".We summarize the significances in Table 2.

Reducedmax Analysis
A pure-mock covariance is ideal for deriving the detection significance, yet the high dimensional data vector leads to a noninvertible sampling matrix due to the limited number of mocks.Reducing the number of channels considered makes it possible to use a covariance matrix taken directly from the mocks, with no need for the analytic one.
For max = 1, there is just one parity-odd channel, and for max = 2 there are four.Fig. 14 shows the detection significances for CMASS and LOWZ with max = 1 in the first and second rows and max = 2 in the third and fourth rows.For max = 1 the histograms of the mocks' T 2 values agree well with the expected T 2 distribution, both for the compressed and the pure mock-based covariance approaches.However, the histogram of the mocks' χ 2 values from the direct approach with analytic covariance deviates from the expected χ 2 distribution, which indicates that the analytic covariance is imperfect.Nonetheless, all three methods show consistent detection significance.For max = 2, all three methods are still consistent.However, we notice the changes in detection significance of CMASS and LOWZ are different when going from max = 1 to max = 2.We suspect this is due to the difference in the two samples, given their redshifts and number density, which we will leave for future exploration.After all, this reduced analysis serves as a robustness check of the pure analytic covariance matrix approach.
LOWZ, 10 bins, max = 4  10 but for LOWZ; again using ten radial bins.As for CMASS, in the compressed approach, the detection significance first increases with an increasing number of eigenvalues and then drops.In general, the detection significance is comparable to CMASS for both compressed and direct analyses and for split and joint sky.Again in the direct method, we see a noticeable deviation of the mocks' histogram from a χ 2 distribution; this again suggests that the analytic covariance may not be unbiased.We also see noticeable deviation at N eig = 800 (and even at N eig = 500), to a much greater extent than for the same N eig values in CMASS (cf.Fig. 10).Using finer binning we find even higher detection significance, both in our compressed and direct analyses (cf.Fig. 10).Unlike the ten-bin CMASS analysis, which peaked at N eig = 200, here the detection significance monotonically increases with N eig .This behavior is likely due to how the covariance matrix structure and signal-to-noise in each channel and bin combination balance with the penalty paid for adding more degrees of freedom if they do not contain much additional S/N.Given that this balance can change depending on the radial binning, we would not have expected that the behavior of the ten-bin and eighteen-bin compressed analyses with N eig would be uniform.11 but for LOWZ with eighteen bins.As is the case for CMASS, the detection significance increases when using the finer binning (cf.Fig. 11).However, the increase is more moderate than the increase for CMASS.The detection significance rises monotonically with N eig , in agreement with the trend for the ten-bin LOWZ analysis of Fig. 11.

Cross-Correlation Between NGC and SGC in Each Sample
A cosmological parity-odd correlation is expected to appear isotropically when decomposed using our choice of basis function.This would manifest itself as a correlation between the signals in the SGC and NGC.Noise in the observed 4PCF, would however not be correlated between hemispheres or between samples, and would reduce any underlying cross-correlation.We compute the Pearson correlation coefficient rp between the signals in the SGC and NGC, with some modifications which we now describe.The naive Pearson coefficient does not incorporate covariance between different points within each data stream being cross-correlated, nor does it account for varying errors from point to point.For instance, two neighboring triple-bins within a given channel in NGC are likely to be much more correlated with each other than two very different triple-bins.Thus, if it happens that one such triple-bin is highly correlated with its analog in SGC, then the neighboring combination in NGC is also likely correlated with the neighboring combination in SGC.This is not unphysical, but it should be taken as less strong evidence for correlation than if two highly independent triple-bins in NGC each showed a strong correlation with their analogs in SGC.The Pearson coefficient would not distinguish.This motivates us to rotate our measured data streams from NGC and from SGC to a basis where each data vector (within the NGC or SGC) is independent; we thus work in the eigenbasis of the analytic covariance matrix.This procedure will be affected by any issues with the covariance matrix estimation, which is present for both analytic and mock covariance.
Another issue with naive Pearson correlation is that some channels and triple-bins have larger statistical errors than others.They may therefore more easily be outliers and drive a spurious Pearson correlation.To correct this, once we are in the "independent" basis as above, taking the square root of the associated covariance matrix eigenvalue as its error bar, we divide each value by its error bar.
For the above procedure, if we had the true covariance, rotating into the true eigenbasis would preserve any correlations perfectly.However, as discussed in §4.2 the analytic covariance is likely imperfect.We believe it to be close enough to the true covariance to preserve correlations within a given sample (CMASS or LOWZ) after rotation, but that trying to rotate into the different eigenbases implied by the covariances of the two different samples and then cross-comparing would not be robust.In particular, the different number densities in CMASS and LOWZ (see Table 1) will non-trivially impact the eigenvectors for each, likely randomizing any underlying correlations.Given this issue (explored more fully in Appendix C), we do not report the results of a CMASS cross LOWZ analysis in this work.
Having outlined how we render our measured 4PF coefficients suitable for a cross-correlation analysis, we now define the .Reduced max analysis for ten radial bins.We use three different methods of obtaining the covariance matrix but always with max restricted to 1 or 2. The left column uses the compression scheme.The central column uses the covariance matrix obtained directly from the mocks, which is possible for this reduced max analysis.The third column uses the covariance obtained from adjusting the analytical covariance to fit, as well as possible, the purely mock-based covariance.The solid curves show the expected distributions for T 2 or χ 2 while the histograms show the distributions obtained from the mocks.The solid lines indicate results from the data.There is overall good agreement in detection significance between the purely mock-based covariance (central column) and the analytic covariance (right column). .Left: Scatter plot for the 4PCF with ten bins measured from CMASS in the decorrelated basis (described in §5.2) and scaled by the covariance matrix eigenvalues, with both Pearson (rp) and Spearman (rs) correlation coefficients given in the legend.The extended panes in the left panel display the histogram overplotted with a Gaussian fit.Right: PDF f (r) (black curve; Eq. 22) for the Pearson correlation coefficient and the Pearson and Spearman correlation coefficients obtained from Patchy (orange and blue histograms).The 95% confidence level (dashed vertical lines) corresponds to an rp = 0.037.Our Pearson rp of -0.014 corresponds to a p-value of 0.460, so there is no statistically significant correlation detected.statistic we use, the Pearson correlation coefficient rp.
where di is the i th element of the 4PCF data vector di with i = 1 for the NGC and i = 2 for the SGC.Both are in the decorrelated, normalized basis for their respective hemisphere. . . .all denotes the average over the elements of the vector, and N d is the number of degrees of freedom.
The PDF for the Pearson correlation coefficient if there is no correlation is (Kenney 1947;Hotelling 1953) where B is the Beta function.Fig. 15 shows the correlation analysis results for CMASS using ten bins.The left panel is a scatter plot of the (decorrelated and scaled) data for CMASS NGC and SGC.The histograms in the extended panels can be well-fit by a Gaussian, showing this assumption of the Pearson coefficient is satisfied.We find no statistically significant correlation at the 95% confidence level.We also report Spearman's rank coefficient, rs, in the Figure (see footnote 10).The right panel shows the PDF for rp (black curve; see Eq. 21) if the two datasets are both normally-distributed.As noted above, a correlation at some level might be expected if the signal were of cosmological origin.Meanwhile, there are some systematics where such a correlation would not be expected, and the same goes for some instrumental artifacts.On the other hand, artifacts at the telescope could correlate NGC and SGC.
Fig. 16 is similar to Fig. 15 but for the LOWZ sample.Fig. 17 shows the correlation analysis results using eighteen bins, where the data here also follow a Gaussian distribution (although we do not show them explicitly).In the right panel, for LOWZ, the correlation coefficient is a negative number at more than 95% confidence.However, we also find that the correlation coefficients are sensitive to the covariance matrix.Changing the number density used in the analytic covariance matrix by 10% can flip the LOWZ correlation coefficient to a positive number (but one consistent with zero).Overall, we thus conclude that there is no evidence for a positive correlation between the N and S for LOWZ or for CMASS when using eighteen bins.Our correlation study results for ten and eighteen bins and LOWZ and CMASS are summarized in Table 3. .PDF of f (rp) for eighteen bins; the left panel displays results using CMASS, the right using LOWZ.The correlation coefficient for CMASS is consistent with zero, but we see a negative correlation for LOWZ; this is further discussed in the main text.

Linear Toy Model in the Eigenbases to Explore χ 2 and rp
To gain an intuition for the high detection significance and the low correlation between the NGC and SGC in the LOWZ, we construct a toy model.Henceforth we use the Pearson correlation coefficients and χ 2 obtained from ten bins results.The purpose is not to constrain real physics, but simply to gain an intuition for whether the values of χ 2 and rp found in our data are plausibly consistent with each other.
We assume the signal is given by a toy model that is linear in the eigenbases of the analytic covariance matrices parametrized by a slope m and an intercept b as s = mx + b, where x is the index of each eigenvalue.We produce two sets of fake data.Each is drawn randomly from a Gaussian distribution, the mean of which is the linear model, and the width of which is the square root of the eigenvalue of the analytic covariance matrix for either CMASS or LOWZ as appropriate.
We separately maximize the log-likelihood for the data vectors dcmass = {rp,cmass, χ 2 cmass } and d lowz = {r p,lowz , χ 2 p, lowz }, respectively.If there is a true primordial signal, we might expect that the parity-odd 4PCF in LOWZ has an overall amplitude 2× that in CMASS, folding in the difference in growth rates and linear biases.Therefore we search over different ranges in the slope and the intercept for each: mcmass ∈ [1 × 10 −8 , 2 × 10 −7 ], the intercept bcmass ∈ [1 × 10 −5 , 5 × 10 −4 ] for CMASS and In Fig. 19 is the log-likelihood plotted for slope as a function of intercept.We find the best fit slope and intercept (the bluest region in the log-likelihood): mcmass = 5.3 × 10 −8 and bcmass = 1.1 × 10 −4 for CMASS; m lowz = 1.5 × 10 −7 and b lowz = 2.7 × 10 −4 for LOWZ.In Fig. 20 we show a comparison of the correlation coefficients and χ 2 values obtained from the real data to the distribution of those from the fake data.The fake data log 10 ( ln ) Figure 19.Log-likelihood for our two-parameter linear toy model s = mx + b ("fake signal").We search for m and b so as to make the observed χ 2 and rp consistent with each other in each sample.The fitting is independent for CMASS and LOWZ.In the eigenbasis of the 4PCF analytic covariance for either CMASS or LOWZ, we produce our linear model and then add Gaussian noise drawn using these eigenvalues to produce 200 realizations of our "fake signal", each with an associated χ 2 and rp.This then lets us obtain the covariance matrix of χ 2 and rp in this toy model.We may then compute, for a given m, and b, their likelihood given the actual observed χ 2 and rp from the data.The blue region is the highest likelihood.We display these results for both CMASS NGC cross SGC and LOWZ NGC cross SGC.It is notable that the best-fit parameters for LOWZ are roughly two times those for CMASS; this did not have to be the case, as the fitting is independent for each sample, but is in fact the scaling we would naively expect were the parity-odd signal truly of primordial origin.Our interpretation of this linear toy model is further discussed in §5.2.
distribution in different colors is generated with (i) a zero signal (ii) a signal with the best-fitting parameters and (iii) a signal with both a steeper slope and a higher intercept.This plot shows that it is possible to explain the low correlation and noticeable detection significance with a simple linear model.13This model's purpose is not to constrain any non-standard inflationary physics.Rather we wish to demonstrate that it is generically possible to have high χ 2 and low correlation coefficients, even with a simple functional form for the model.19), and from a larger "fake" signal in purple to indicate the possible range of outcomes.The fake data at each index (see Fig. 18) is drawn from a Gaussian distribution with mean given by the linear model and width given by the square-root of the eigenvalue of the analytic covariance matrices.CMASS results are shown in the upper row and LOWZ in the lower.We see that there exists a "fake signal" that does explain the observed rp and χ 2 from the real data.While our "fake signal" should not be taken too literally, this does indicate that, generically, our observed χ 2 and rp are not necessarily in tension with each other.

POTENTIAL SYSTEMATICS
Parity-breaking in large-scale structure is not expected in the standard cosmological paradigm.Hence it is important to assess whether systematic errors could produce an apparent parity-breaking signal.We must consider systematics both at the signal level and at the covariance level.These latter must be considered because they could cause an underestimation of the covariance, thus leading to an overestimation of the detection significance.
We divide these systematics into three categories: survey-related effects, observer-related effects, and procedure-related effects.We discuss these in respectively §6.1, §6.2, and §6.3.We have already considered the possible impact of systematics on the covariance matrix in §4.2.2.

Survey-Related Effects
Spectroscopic surveys provide galaxies' 3D positions, but estimating the correlation functions requires the excess probability over random of finding an N-tuplet of galaxies in a given configuration.This in turn demands precise knowledge of the selection function.

A Toy Model for Errors in the Selection Function
Despite that great care has been taken in BOSS to model the impact of astrophysical foregrounds (e.g.stellar density), or observational conditions on the survey selection function, we still seek an intuition for the impact of an imperfectly-modeled selection function on the parity-odd modes.In particular, our goal is to see if an even-parity 4PCF is converted into an odd parity one by an improperly corrected selection function.For this goal, it is sufficient to begin with a spatially unclustered density field (which would produce only an even-parity 4PCF).
The selection function encapsulates how our selection criteria (e.g.color and magnitude cuts) samples the underlying distribution of galaxies in the universe.The underlying density field can be inferred from the observed number density n(s) and the selection function f obs (s).Since we have no access to the underlying distribution, we cannot disentangle how our selection criteria shape the sample we end up observing from any underlying variation in the universe (usually specifically along the line of sight).In this discussion we allow the deviation of the selection function to be general and depend on 3D position s.
Our estimate of the true number density of objects is where n is the observed number density and f obs our estimate of the selection function.The actual true number density is where f obs is the true selection function.
Hence we may relate our estimated true number density to the actual true number density as Our density fluctuation field is then We have assumed that we correctly estimated n and also that n true (s) ≡ n, i.e. in this simple toy model, there is no underlying clustering.We require that the integral of g over the survey will be unity; we use this condition to normalize g shortly.
We now make a toy model for g to enable assessment of its possible impact.We take it that it is factorizable as and that along each direction it is a power law as gx = I −1/3 (x0 + ∆x) px , gy = I −1/3 (y0 + ∆y) py , gz = I −1/3 (z0 + ∆z) pz , (28) where we have written ∆x, ∆y, and ∆z as displacements around a central point (x0, y0, z0) and have defined the normalization I as This normalization ensures that the integral of g over the survey is unity and the Li are the lengths of the box sides in each direction.The primary galaxy for our computation of the spherical harmonic coefficients a m (see §2, as well as Philcox et al. 2021a for a detailed discussion of how these coefficients enter our algorithm) is at s0 and each secondary galaxy is displaced by ri from the primary.If we now assume that the true density field is uniform, its change due to our modulation will be (for the primary) For the secondaries (in this case, the three field positions around the primary, which, with it, build up that primary's contribution to a given 4PCF channel) we have, by taking a Taylor series for g to leading order: with ri ≡ si − s0.We now assess whether this systematic can contribute to the parity-odd 4PCF.When we multiply out the three secondaries, we can produce terms involving one, two, or three factors of ∇g.All of the parity-odd basis functions are proportional to r1 • (r2 × r3), so only the term involving three factors of ∇g may possibly contribute.We now compute the projection of this term onto our basis.Denoting it ζg, we have where we used a vector integral identity to obtain the second line.
Normalized number count distribution along the x-, y-, and z-axes of mocks with a spatially-random underlying distribution but modulated by a power-law in each direction ( §6.1.1,Eq. 28), as indicated in the legend; the function f shows the modulation along each axis, and the variable after the colon which axis is plotted.We show two cases along each of the three axes; powers px = py = pz = 1 and then px = 1, py = 2, pz = 3.
Several parity-odd 4PCF coefficients as measured from our spatially-random power-law-modulated mocks ( §6.1.1).We show results from 100 mocks with number density n = 3 × 10 −4 [h −1 Mpc] −3 for no modulation (tan), x 1 y 2 z 3 (red), and xyz (grey).We find that the χ 2 for the power-law-modulated mocks is smaller than that for the unmodulated ones, implying that this modulation does not produce a significant parity-odd 4PCF at the signal level.
Though the above shows that on average (in an infinite volume), a power-law modulation would produce no spurious parityodd signal, such modulation could increase the fluctuations observed.If such a systematic were present in the data but not the mocks, then the covariance used (from mocks or calibrated by mocks) would be an underestimate, and this could produce an apparent detection.Hence, we explore this power-law toy model numerically to assess if such an effect might occur.We construct 100 mocks whose points have a spatially random distribution, each in a cubic box with side length L box = 500 h −1 Mpc and number density n = 3 × 10 −4 h −1 Mpc −3 .The distribution of the mock number density along the x-, y-, and z-axes is shown in Fig. 21.We investigate two cases, px = 1, py = 1, pz = 1 and px = 1, py = 2, pz = 3.We also create a set of 100 mocks with the same box size and number density but with a uniform distribution along the three axes.For all cases, we use the same random catalog with a uniform distribution at all spatial positions to convert the density field into a density fluctuation field (the only relevant part of edge correction for a periodic box).Fig. 22 shows the 4PCF coefficients measured from the no-clustering mocks.From left to right, we plot the channels In each panel, we compare a parity-odd 4PCF coefficient for an unmodulated sample and two differently-modulated samples.As expected, the 4PCFs measured from the uniform sample have very small amplitudes.The χ 2 for the uniform sample is comparable to those modulated samples.Overall, this means that the power-law modulation would not produce a spurious detection of parity-odd modes at the signal level.

Density Distortions, Selection Function, and Fiber Collisions
In theory, one can expect parity-odd signals to be induced solely by 3D modulations of the density field because such modulation could convert a parity-even mode into a parity-odd mode.A 2D screen uncorrelated with the underlying density field can be separately averaged over rotations, and the rotation will force any parity-odd piece of the screen to vanish.An example of an uncorrelated screen is that given by the veto mask, which was designed to remove angular regions contaminated by foreground bright stars, plate holes, etc. and is not correlated with the underlying 3D galaxy density.
Another 2D effect is "fiber collisions" (FC).BOSS uses fibers to guide the light of the observed objects from the focal plane to the spectrograph.Each fiber has a limited angular diameter (62 ) and multiple galaxies that fall within this radius cannot be resolved within one exposure.In the data, fiber collisions are corrected by upweighting the nearest angular neighbor (Ross et al. 2012;Reid et al. 2016).This approach can potentially introduce long wavelength mode-coupling because it identifies neighbors only on an angular basis.Although one can model the fiber-collision effect as a top-hat function (Hahn et al. 2017), the exact physical scale could still be redshift-dependent (see Hou et al. 2021b for quasars, which span a wide redshift range).
To demonstrate the impact of these angular effects on the parity-odd modes, we turn off the veto mask weights implemented in the Patchy mocks as well as the fiber collision weights.By switching on and off the fiber collision weight we can see the impact of the nearest-neighbor upweighting on the parity-odd detection significance.The left panel of Fig. 24 shows the distribution of the resulting T 2 values compared to that expected.When the covariance is inferred from the contaminated mocks themselves, there is good agreement between the contaminated mocks' T 2 distribution and that for the standard case.We, therefore, conclude that these effects do not induce parity-odd modes at the signal level.
We also explored the effect of the fiber collision weights on the signal as found in the real CMASS data.We set these weights to be identically unity on the data, thus undoing the upweighting of the nearest angular neighbour of a galaxy lost to fiber collision.We also make this same change to the fiber collision weights on the mocks, and then remeasure their 4PCF and refit our analytic covariance matrix template using it.
We found that with this setup, the ten-bin parity-odd detection significance increased by 6σ.We also computed the parityeven connected 4PCF, to assess if an unaccounted-for error of this type in fiber collision weights on the data would create tension between data and mocks.Indeed, we found that with all fiber-collision weights on the data set to unity (as if one had made a gross error in the data weights), but with standard weights on the mocks (and also in fitting the covariance template), there would be a 5σ disagreement in the parity-even sector between mocks and data.This argues that any serious error in the fiber-collision weights that impacts the parity-odd sector can be controlled by enforcing consistency in the parity-even sector.
We now briefly discuss what we regard as the most likely reason the parity-odd significance in the data increased in the test above.We know from the tests on the mocks alone that an error in the fiber-collision weights cannot induce a true parity-odd signal where there was initially none (left panel of Fig. 24).Thus, the increased detection significance in the data is presumed to result from the change in the fiber collision weights' increasing the variance in the data more than it does in the mocks used to estimate the covariance.
We suspect this disproportionate increase is due to the following.In the real survey, each observational tile is allowed to be overlapped with others to maximize the fraction of targets that can be assigned fibers.Thus, uncorrected fiber collisions likely have a lower impact on the data than on the mocks.In particular, even with the fiber collision weights set to unity, there is likely better sampling of the densest regions on the sky in the data than in the mocks.This would tend to increase the data's variance relative to that of the mocks in the situation where neither are corrected for fiber collision.
In addition to the angular effects, we also distort the radial selection function by introducing a 10% scatter in the n(z) of the randoms.At each redshift in Fig. 23, we draw a value from a Gaussian with mean 10% and width 5%, and then decide randomly and with equal probability whether to give it a plus or a minus sign.We then use this value to set a new probability of a random particle at that redshift being selected.This procedure results in the red curve in Fig. 23 for the shuffled randoms' n(z).The FKP weights for the randoms are also recomputed after this.
The distribution of T 2 values for the Patchy mocks with shuffled n(z) when we use the compression scheme with the covariance inferred from the contaminated mocks agrees well with that obtained using the standard weights (black histogram) and also with the predicted T 2 distribution (black curve) as shown in the left panel of Fig. 24.However, the radially-distorted randoms lead to an increase in the covariance, which is interestingly in contrast to mocks with angular contamination.As a result, the mock distribution shifts towards higher χ 2 when we use an analytic covariance matrix calibrated with respect to the standard mocks (right panel of Fig. 24).In other words, had this systematic been present in our BOSS data, the covariance matrix inferred from the standard Patchy mocks would have been underestimated.In table 5 we estimate how much this effect could have been affecting our detection significance for ten and eighteen bin tests, respectively.

Magnification Bias
We consider magnification bias qualitatively.If the underlying density field produces non-zero amplitudes only in even-parity 4PCF channels, it is sufficient to ask whether magnification bias will preferentially affect tetrahedrons of one handedness.Much like the method of image charges, we may consider a notional pair: a tetrahedron and its mirror image.As long as we can prove that for any such pair, there is no effect, no parity-breaking will be introduced.
Consider for simplicity a tetrahedron for which the three galaxies defining a triangular base are at the same redshift, and the fourth point is at a higher redshift.If the three "base" points are closely enough concentrated, they may magnify the fourth galaxy behind them and on average make it more likely to be selected than otherwise.However, this will be the case also for the The selection function applied to the Patchy mocks is designed to match that from the BOSS data for both the mock data catalog (black histogram) and the random catalog (grey histogram).To test whether an imperfect inference of the radial selection function can lead to a spurious parity-odd signal, we deliberately distorted the n(z) for the random catalog (red histogram), by shifting it by 10% on average, where the percentage shift is drawn from a Gaussian with that mean and width of 5%.We show the distribution of T 2 values for Patchy mocks analyzed using these "shuffled" randoms in Fig. 24; it is consistent with no detection of a parity-odd 4PCF.T 2 values calculated with 50 eigenvalues using our data compression scheme ( §5.1.2) with the covariance matrix inferred from the systematics-integrated mocks themselves.The black curve is the expected T 2 distribution.The histogram for the Patchy mocks with the standard weighting scheme is in black, with the close pair weight turned off to assess the impact of fiber collisions in blue, with the veto mask turned off to assess the impact of the angular mask in orange, and finally, in red, the results of shuffling the randoms' redshift distribution by 10% as described in Fig. 23.None of these alterations to the 4PCF analysis produce a distribution of T 2 values that are inconsistent with that expected under the null hypothesis of no-parity-odd 4PCF (black curve).Right: χ 2 values calculated using the analytic covariance calibrated to mocks analyzed with the standard weights.Here we can see that the χ 2 distribution has a higher mean for the shuffled mocks than for the standard mocks, suggesting that although such a systematic does not induce a parity-odd signal, it might lead to underestimation of the covariance.
mirror image of this tetrahedron.Indeed, magnification will be invariant to the sign-flip of the two plane-of-sky coordinates; at most, it can be parity-odd in z, but any such contribution from the z behavior alone will vanish after rotation-averaging.

Splitting the Sky into Angular Patches
To assess if our detection is coming from one part of the sky preferentially, which might indicate a systematic, we split the sky into eight angular patches, each of roughly 40 × 35 deg 2 (right ascension, RA, times declination, DEC), and run our analysis on each independently.Each angular patch implies a subvolume (hence SV) as we extend it out in the redshift direction.The For six parity-odd channels, we display the difference between the 4PCF coefficients of the eight sub-volumes created by extending these angular patches in the redshift direction, and those from the full NGC sky.We divide these differences by the standard deviation of the Patchy results.The colors for the eight sub-volumes match in the left and right panels.Overall, we see that no one angular region stands out as producing a great deal more parity-odd signal than another.We offer a more rigorous demonstration of this claim in Fig. 26, which shows the χ 2 for each sub-volume.
right panel shows some examples of the 4PCF measurement and compares them to the standard deviation of each subvolume.
We estimate the covariance of each subvolume by splitting the Patchy mocks in the same way.Since the subvolumes receive modulations from Fourier modes larger than them, "beat-coupling" leads to additional corrections to the covariance matrix (Li et al. 2014;Putter et al. 2012).Thus a naive rescaling of the covariance by its volume relative to the total survey volume would likely give an underestimate and thence produce artificially high T 2 values.Fig. 26 shows the T 2 computed from the entire BOSS CMASS sky (thick dashed vertical black line) and the eight subvolumes (thin colored vertical lines).The expected T 2 distribution is the black curve.We combine the eight subvolumes by taking an inverse-variance weighted sum of the T 2 values.Using the diagonals of the covariance matrix (see the right panel of Fig. 26) we find the combined T 2 = 128, which differs from the one from the entire BOSS CMASS sky T 2 = 138.We repeat this using the full covariance of the subvolumes and find T 2 = 129, quite similar to the result using just the diagonal.We see that subvolume 2 and 3 have the highest detection significance compared to the average; given the errorbar of the χ 2 shown in the right panel of Fig. 26, this corresponds to a 1.8σ difference.We also note that combining the signal at the estimator level is not a simple linear operation due to the edge correction (see also Eq. 8).Moreover, since we use the same analytic covariance for different subvolumes, the amount of information picked out by the ranked eigenvalues could differ for each subvolume.We, therefore, expect that there is a difference between analyzing the data on the whole sky or splitting it and then combining it.
We also repeat this same test on the SGC and on the LOWZ sample.Overall, we do not find any subvolume that has a particularly high or low detection significance compared to the average.

Redshift Failure
The BOSS fibers are not completely randomly positioned in an observing tile.The resolving power of the fibers is impacted by their positions on the CCD cameras.Those fibers near the edge of the spectrograph slit are less likely to obtain a high-quality spectrum than are more central fibers (see §4.2.2 of Smee et al. 2013).Consequently, redshift failure occurs in a way that depends on 2D position on the telescope plate.This dependence is not wholly symmetric under a parity transformation of the x and y axes on the plate.Now, when a redshift failure occurs, the nearest angular neighbor of the failed galaxy is upweighted.The neighbor's redshift should be a fair draw from the redshift distribution of the survey.However, redshift failure is statistically slightly more likely for a galaxy at higher redshift, as it will on average be fainter, and hence more affected by the reduction in spectrum quality from being at the edge of a CCD.Thus, replacing it by upweighting a galaxy fairly drawn from the survey's redshift distribution can actually produce a bias.Effectively, we will systematically be pulling galaxies towards lower redshift as a function of 2D position on the plate.This 2D function, as already noted, is not fully symmetric under (x, y) → (−x, −y) and hence, in combination with the coupling to z produced by the above, can in principle imprint a parity-odd signal on the density.
We now seek to explore this issue.CMASS has 1.8% of its objects undergo redshift failures, while LOWZ has only 0.3% do so (Ross et al. 2012).Thus we mainly focus on CMASS.Since the redshift failure galaxies are already removed from the catalog, we explore this effect by creating more redshift failures and seeing whether it appreciably alters our parity-odd signal.25.The T 2 value for the entire CMASS NGC is the thick black line, the eight subvolumes' T 2 values are in thin, colored lines, and the combined T 2 from the eight subvolumes are in dashed black.The predicted T 2 distribution (under the null hypothesis of no parity-odd signal) with 100 degrees of freedom is the black curve; we have used the compressed analysis described in §5.1.2.Right: Covariance for the T 2 values of the eight sub-volumes, estimated using the Patchy mocks.This covariance is needed to properly co-add the T 2 values for the eight sub-volumes to report the combined T 2 given in the lefthand panel (dashed line there).The slight non-vanishing off-diagonal indicates correlations between the eight sub-volumes, but this turns out not to affect substantially the combined T 2 value.We may observe that the combined T 2 value differs from that for the whole NGC by about 10%; we believe this is due to missing the produced by quartets that are split across more than one sub-volume, which go uncounted in our sub-volume analysis.We would expect that adding them in (as the whole NGC analysis does) would indeed raise the χ 2 .
We first infer the redshift efficiency according to the inverse of redshift failure weight, then we apply a fit and assign each fiber a redshift efficiency rate accordingly.We apply a magnitude cut with mag cut = 19.20 and from which a galaxy subsample is selected weighted by the fiber efficiency rate such that the subsample is 1.8% of the total galaxies.We then identify the faint and bad fiber galaxies' neighbors and select the ones that are closest in angular distance.Finally, we upweight those galaxies' redshift failure weights by one and renormalize the random catalog's weights such that it matches the total weights of the new data catalog.We repeat the above steps for both NGC and SGC.
Here we do not transfer the total systematic weights, as they cannot be simply propagated to the nearest neighbours.Assuming the total systematic weights to be unity for those galaxies has an impact on the total weights of only 0.1%, so we do not expect an issue to arise due to this.We also retain the current estimation of the sector completeness, as this cannot be easily redone at the catalog level.
We found that doubling the redshift failures increased the total parity-odd detection on the joint sky by 2.5σ when using eighteen radial bins and 1.5σ when using ten radial bins, both with the analytic covariance; these results are summarized in Table 5.This increase could be due to the 3D modulation of the density field producing a parity-odd mode as discussed above, but it could also simply be due to increased variance in the data (i.e. that any parity-odd signal so produced would average to zero in the infinite limit).This increased variance could stem from our inability to correct the sector completeness to account for the removal of more objects.
Given that the angular separation between galaxies and their nearest neighbours is typically much smaller than the sky coverage of a sector, and most of the upweighted galaxies and the removed galaxies are still within the same sector and this, we might expect that the impact of uncorrected sector completeness in our test setup is small.However, we need to quantify its impact.In order to do so, we applied a uniform selection in redshift and selected 1.8% of those galaxies that are close to bad fibers.We proceed by upweighting their neighbours.There we observed 2.2σ and 0.5σ increase in detection significance for the joint sky using eighteen bins and ten bins, respectively.Hence the increase in the detection significance is actually dominated by the imperfect compensation between the random catalog and the real data produced by the uncorrected sector completeness in our test setup, which leads to a spurious variation in the density field.This is qualitatively similar to what we find with the power-law toy model in §6.1.1.)Furthermore, redshift failures also impact the parity-even modes; this impact (for our "doubling" test) is given in Table 5.The good agreement between this mode in the true BOSS data and in the distribution of the mocks (Philcox et al. 2021b) demonstrates that redshift failure at the level it is actually occurring in BOSS should not have a strong impact on the detection significance.6.2 Observer-Related Effects 6.2.1 Redshift-Space Distortions (RSD) Galaxies' radial distances are inferred from their redshifts, however, their peculiar velocities caused by the local gravitational environment add a component in addition to their cosmological redshifts.This distortion in galaxies' positions is known as redshift-space distortions (RSD).In order to understand the impact of the RSD on parity-odd modes, we first consider the "global" line of sight approximation and take the ẑ as the line of sight.RSD modulates the density by the usual Kaiser (1987) factor, which is even-parity with respect to the line of sight.The primary galaxy itself is at the origin and so does not carry any angular momentum; multiplying it with three density fields, with each of them being parity-even, will be parity-even.Given that rotation averaging does not alter the parity, RSD considered in the "global" line of sight approximation cannot produce parity-odd modes.We can easily extend this logic to a "local" line of sight by taking the line of sight to a quartet to be that to the primary.The modulation in the overdensity field due to the RSD will be proportional to (n • k) 2 , but again results only in parity-even combinations.Given that rotation averaging does not alter the parity, certainly, RSD considered in the "local" line of sight approximation cannot produce parity-odd modes.14 .We also note Jeong & Schmidt (2019) have considered parity-odd bispectrum induced by galaxies' radial velocity by introducing explicit coupling to the line of sight degree of freedom, by doing so the total rotational invariant system includes the additional line of sight vector, which by itself has no constraint on parity as long as the combined system has total zero angular momentum.In our work, instead, the line of sight dependencies are already averaged over and we are not sensitive to the parity-odd mode due to the RSD as in Jeong & Schmidt (2019). 15o investigate the impact of RSD, we also use the Patchy mocks in real and redshift space.For this test, the Patchy mocks used are in a cubic box of side length L box = 2, 500h −1 Mpc and with number density n = 6.6 [h −1 Mpc] −3 at redshift z = 0.57.The galaxies are shifted using the exact 3D velocities of each galaxy in the mock.A comparison of the 4PCF in real-and in redshift-space is shown in Fig. 27.We do not observe an overall enhancement in the amplitude for the 4PCF.

Tidal Alignment
The large-scale tidal field impacts galaxies' orientations, which might lead to an anisotropic galaxy selection.For instance, in the case of spirals, an imaging survey might be more likely to detect a face-on rather than edge-on disk.Even LRGs are stretched ellipsoids, with aspect ratios of order 2:1, so orientation could impact the selection, though the effect would likely be smaller.For LOWZ, the sample is quite pure (LRGs), but for CMASS, the sample is roughly 26% disky blue galaxies.We might thus expect any anisotropic selection would manifest more strongly in CMASS.
However, given that the observer's orientation is random with respect to the tidal field (and thus the orientation of the galaxies being selected), we suspect that such an effect could not produce parity-odd modes.If it did, we might expect a stronger detection of them in CMASS for the reason noted above.
However, for completeness, we here work out what a tidal model might do if the selection depended upon it.Following Hirata (2009), the observed galaxy density fluctuation can be expressed as where is the selection function at position x, and it depends on the viewing direction n.The average of the selection function over the unit sphere is zero.By treating the large-scale tidal field as a perturbation, the lowest-order Fourier-space contribution (after Taylor-expanding) is (Hirata 2009 equation (30)): where A is an amplitude and δm is the matter overdensity.We see that the selection function is parity-even, and so multiplying it onto a single density does not alter the parity of any spherical harmonic expansion coefficient associated with that density field.Hence, if there is on average signal in only parity-even channels in the underlying galaxy 4PCF, modulating each density point by (n|k) above will not convert them into a signal in parity-odd channels.
6.3 Procedure-Related Effects

Fiducial Cosmology
In order to convert the sky coordinates to Cartesian coordinates, we need to assume a fiducial cosmology.Using a fiducial cosmology that differs from the true cosmology may lead to distortions in the clustering.At the signal level, if any galaxy quadruplets have their parity flipped by such an error, it is equally likely to happen for a given tetrahedron and for its mirror image.Therefore, on average, even if a "parity-flip" occurs, it will be canceled out after averaging over all galaxies in the survey.The radial bin separation ∆r we choose is 14 h −1 Mpc when we use ten radial bins and 8h −1 Mpc when we use eighteen radial bins.Beginning with the cosmology of Planck Collaboration et al. ( 2020) and moving ±3σ, the change in the angular diameter distance is less than 1%, which will not be resolved by our radial bin choice.At the covariance level, a different cosmology could lead to a different estimate of the best-fit number density and volume (for the analytic covariance) and would also change the power spectrum shape.Even though we use the same cosmology to convert the galaxies in the real data and the Patchy mocks, the statistical fluctuations do not necessarily scale the same way as we assume a fiducial cosmology that differs from the true one.We can see a 6% change in the volume when distorting our fiducial cosmology by 3× the constraints on Ωm obtained in Planck Collaboration et al. (2020).Such a difference corresponds to an overestimation of the detection significance by ∼ 2.5σ (see Table 5).However, if there is an obvious difference between our choice of fiducial cosmology and the underlying one, it will also be visible in the comparison of the parity-even 4PCF measurement.The fact that we saw good agreement between the observed parity-even 4PCF and the mock distribution suggests wrong fiducial cosmology is not relevant here.

Survey Geometry and Edge Correction
As discussed in §2, the 4PCF measurement requires correction for the survey geometry as given in Eq. ( 8).In principle, there can be a mixing between parity-even and parity-odd channels due to the survey geometry.The left panel of Fig. 28 shows the coupling matrix M 1 2 3 , 1 2 3 as in Eq. ( 9) from Patchy mocks with CMASS NGC geometry at a radial bin combination (r1, r2, r3) = (27, 55, 97) [h −1 Mpc].This coupling matrix has off-diagonal elements, which means that the parity-even mode can couple to the parity-odd mode.To assess the correction of the parity-odd modes due to the survey geometry, we plot the ratio of the random coefficients to the lowest-order one, R 1 2 3 /R000 (right panel of Fig. 28).These ratios are of order 10 −6 ; the dominant edge correction is thus simply R000, which does not mix parity.In addition, the fact that the Patchy mocks (which also get edge-corrected) are consistent with the predicted distribution indicates that survey geometry correction is not inducing spurious parity-odd modes.
We also may develop a toy model for the edge-correction factors based on Slepian & Eisenstein (2015b) §4.3, where the impact of a planar edge on a uniform density sphere was considered.We consider a sphere of radius R some distance away from a planar edge of the survey, and orient the z-axis perpendicular to that plane.We may define a critical angle cosine µc; for µ < µc, a point will fall outside the survey.We may compute the a m as a function of µc, then form the 4PCF estimate, and finally average over all possible µc, which corresponds to averaging over how far away the sphere's center is from the edge.The choice of z-axis means that only m = 0 a m are non-zero.The setup is described more extensively in Slepian & Eisenstein (2015b) §4.3.
Proceeding in this way, we may find a general form for a 0 , which is given in Slepian & Eisenstein (2015b).Here, we will only compute the first few leading edge correction factors, so we find it more useful to simply quote the first few a 0 explicitly.We have a00  R L 1 L 2 L 3 (r 1 , r 2 , r 3 )/R 000 (r 1 , r 2 , r 3 ) 1e 6 The odd coefficients of the randoms normalized by R 000 .We have mapped the odd angular momentum triplets given by Λ to a 1D index (vertical axis), and show the radial bin also so mapped on the horizontal axis.Here we see the dimensionless parity-even-to-parity-odd mixing is of order 10 −6 , showing that the dominant change to the odd-parity modes from edge-correction is due to R 000 , i.e. overall renormalization from dividing by the mean number density.
Table 5.Here we summarize the shifts in detection significance, in units of σ, due to various systematics for both the parity-odd (top two rows) and parity-even (bottom two rows) modes, including redshift failure ( §6.1.5),the impact of the fiducial cosmology ( §6.3.1), and a distorted selection function ( §6.1.2).We used the N+S of CMASS and the analytic covariance for these calculations.These lead to, respectively, f011 = 8%, f112 = −1%, and f022 = 2%, where these are the ratios of the R with the indicated angular momenta to R000.Furthermore, we note that only of order 20% of spheres in the survey impinge on a boundary at all for the BOSS footprint; hence these factors would be expected to scale down by roughly a factor of five.These estimates show that the edge correction is not a substantial alteration to the measured 4PCF, and indeed that the dominant effect is that of R000, with the higher-order edge-correction factors contributing of order a few percent at most (after division by the factor of five).

CONCLUSION
In this work, we have used a novel, recently-suggested method (Cahn et al. 2023) to search for parity-violation in 3D large-scale structure.This method exploits the fact that the lowest-order shape that cannot be rotated into its mirror image in 3D is a tetrahedron, and thus counting tetrahedra and looking for excesses of one type over another can reveal parity violation.We have detected strong evidence for it in both the CMASS and LOWZ samples of BOSS, and in both the NGC and SGC of each, summarized in Table 2.We have explored numerous systematics to assess whether any could produce either a spurious parity-odd signal, or a substantial change to the covariance matrix used such that we would underestimate the statistical error bars so as to make a spurious detection.We have not found any systematics that seems to do either.We have also performed a number of tests of our covariance matrix and pursued a number of analysis variations both as regards covariance matrix and radial binning, scales, and angular momenta used.Overall, the compact takeaway from this work is that CMASS with eighteen radial bins has a 7σ detection of parity-odd modes when using the analytic covariance with all degrees of freedom, and a roughly 4σ detection when using an order of 500 eigenvalues plus a covariance estimated directly from the mocks.Given a large number of degrees of freedom in our analyticcovariance analysis, underestimating the noise by 20% in every degree of freedom would be enough to produce this detection spuriously.In the "compressed" analysis with 500 eigenvalues, we would have had to underestimate the covariance by 30% in every degree of freedom to spuriously produce our 4σ result.Such a substantial underestimate seems unlikely but given that the result, if cosmological in origin, would be significant, caution remains warranted. 16e now briefly revisit the key technical aspects of our analysis and of systematics search.

Radial Binning
In this work, we tried three different radial binning schemes and found different detection significances.We hypothesize that the increasing detection significance is due to the reduction of internal cancellation.The finest binning we use in this work has ∆s = 8 h −1 Mpc, which is roughly half of the average inter-particle separation for BOSS.In this work, our main goal is to explore a possible parity-odd signal; we leave the exploration of the optimal radial binning to future work.
CMASS vs. LOWZ Detections CMASS and LOWZ are both mainly composed of LRGs, but the two samples behave differently with regard to i) sensitivity of the detection significance to varying radial binning and ii) sensitivity of the detection significance to varying max.These different behaviours may be due to differences in the number density, mean redshift, galaxy population, and redshift failure rates of the two samples.

Correlation Across the Sky
We do not observe a statistically significant positive cross-correlation between the NGC and SGC in each dataset, and we find that this test is sensitive to radial binning.We explored a simple phenomenological model to understand if high detection significance but low cross-correlation is generically possible, and it seems so.These two measures are algebraically independent so this finding was not outside the realm of expectation.

Systematics
While we have shown that nearly all the systematics considered do not produce parity-odd modes at the signal level (i.e. they would contribute vanishingly in the infinite-volume limit), some can enhance the variance of the observed data relative to the covariance from the mocks if they are present in the data but not the mocks.If substantial enough, this in principle could produce a spurious detection.In practice, we have not found any systematic that, even if grossly larger than what we believe to be present in the data, could increase the variance enough to explain our results.We have also looked for consistency in the even-parity sector as a way to guard against any systematics' unbeknownst to us increasing the variance.

Implications
Given the unexpected nature of the detected signal, it is premature to speculate on the implications.Following the advice of Newton: hypotheses non fingo.17Nonetheless, it is worth placing our results in the context of possible theoretical models that could explain them, should the signal here be of genuinely cosmological origin.Besides Chern-Simons-like interactions (Deser et al. 1982b,a;Witten 1989) with couplings of different forms and various extensions (Freidel et al. 2005;Alexander et al. 2006;Liu et al. 2020;Li & Zhao 2022), string-sourced perturbations (Pogosian & Wyman 2008), primordial vortices (Vilenkin 1978;Vilenkin & Leahy 1982;Brizard et al. 2000), broken symmetry during phase transition ('t Hooft 1974;Quashnock et al. 1989;Baym et al. 1996), anomalies in gravity (Alvarez-Gaume & Witten 1984;Contaldi et al. 2008;Kuntz 2019) can all produce parity-violation.Among these mechanisms, the simplest relevant for us is where there is a non-vanishing parity-odd trispectrum (Fourier-space analog of the 4PCF) in the scalar sector (Shiraishi et al. 2014;Shiraishi 2016), as curvature perturbations directly seed the galaxy density perturbations.Extending such models to include the full wave-vector dependence as well as "late-time" effects such as galaxy biasing and redshift-space distortions will be an important step on the road to comparing models to data and determining if the signal here observed is consistent with any such models.Especially given that the detection significance we report here is a cumulative result from hundreds to thousands of degrees of freedom, a model comparison will provide much more constraining power as it will enable the use of shape information.

Outlook
There are two main directions we suggest on the data side to further pursue the parity-odd sector.First, it is worth validating the measurement on different surveys with different designs.Upcoming surveys such as Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016), Euclid (Laureijs et al. 2011), andRubin (LSST Science Collaboration et al. 2009), which will provide much-larger samples and with different systematics, will offer expanded insight into the signal's true origin.In particular, if our signal is of truly cosmological origin, it would likely be detectable in DESI at much higher statistical significance due to the larger volume.Thus, while in the present paper, making an error in our covariance matrix at the 20-30% level would produce a spurious detection, in DESI, we expect one would have to make a much larger error in the covariance.Since this is unlikely, finding evidence for this signal in DESI would be a strong indicator that either the signal is real, or that it is due to some as yet-undiscovered parity-breaking systematic shared between BOSS and DESI.Yet, given that DESI has a different imaging survey, a different instrument, and a method for fibering objects in taking spectra, and is on a different telescope than BOSS, this latter possibility also seems unlikely.Second, it is worth pursuing further approaches to search out systematics.For example, assuming a linear relationship between the systematics and the observed galaxy density field, one can study cross-correlations between them to reveal residual systematics (Ho et al. 2012).Alternatively, likelihood-based forward modelling can also be applied to search for unknown systematics (Lavaux et al. 2019).The first approach would require a non-trivial extension of our current NPCF algorithm in order to compute cross-correlations between fields, and the second requires injecting inferred systematics maps into mock catalogs.We leave this for future work.
Whether the detection of the current work ultimately turns out to be cosmological in origin or not, we hope that the careful and extensive treatment of systematics here offers a well-delineated avenue forward for future studies of parity-breaking early-Universe physics using 3D large-scale structure.Figure C1.Comparison of the CMASS and LOWZ analytic covariance matrices, but with their volumes set to be equal.We do this because an overall scaling will not adversely impact a cross-correlation analysis, but non-trivial differences in the matrices do and it is these we want to display.We have mapped each set of channels and radial bins into 1D indices to make a 2D plot.The lefthand panel shows a symmetrized test; notably, the diagonal is non-vanishing after subtracting the identity matrix.The non-vanishing off-diagonal elements present in the lefthand panel indicate that the eigenbases of the two matrices also differ, as we might expect given that the shot noise impacts the covariance non-trivially and the two matrices are computed with different shot-noises (see Table 1).The righthand panel shows just the diagonal of our test matrix; the residual is likely due to the slightly different number densities.We have split the sample into three bins in imaging depth.Right: Same as the left but for r-band.We see that the three bins in i-and r-band depths have very similar n(z), implying that the imaging depth is unlikely to impact our analysis.

APPENDIX F: MAXIMUM SCALE CUT AND MINIMUM BIN SEPARATION
To explore whether the use of small scales makes a difference in our analysis (as these are the scales on which the mocks may be most likely to imperfectly mirror the data, due to approximate treatment of non-linear structure formation), we test the detection significances by applying further restrictions to the radial bins.First, we force the minimum radial bin separation to be ∆r ≥ 15 h −1 Mpc such that the results are less sensitive to small scales.Second, we vary the maximum radial bin rmax = 90 h −1 Mpc and rmax = 130 h −1 Mpc such that we are less affected by any mismatch between mocks and data around the BAO position and towards larger scales (as was seen in some of the 2PCF measurements).
The results are shown in Figs.F1 and F2.In all these cases the detection significance reduces compared to the fiducial case (without minimum bin separation or maximum radial bin) given that the number of degrees of freedom is reduced.However, we still observe non-negligible detection significances for all cases.It is worth pointing out that there are in total only 35×23 = 805 degrees of freedom when using rmax = 90 h −1 Mpc, which allows us to use the mock covariance directly.Despite the mismatch

Figure 2 .
Figure2.Left: A box with side length L box = 1000h −1 Mpc filled with tetrahedra (upper left, small panel).Each of them has a unique primary (black) from which the three sides with respective lengths r 1 ∼ 10h −1 Mpc (red), r 2 ∼ 20h −1 Mpc (yellow), and r 3 ∼30h −1 Mpc (blue) extend.The larger panel on the left is a zoom-in on the full box to display the tetrahedra more clearly.Right: A sketch of a "clockwise" tetrahedron and its "counterclockwise" mirror image.The primary is in red.Our convention on clockwise and counterclockwise is detailed in Fig.1.On the left, the red point is closer to us than all the others.Thus the tetrahedron on the left is clockwise, as looking down from the primary, we go clockwise as we move from the smallest side to the largest.On the right, the primary (in red) is behind the other galaxies, so looking down from it towards them will reverse the handedness.Thus the rightmost tetrahedron is counterclockwise as viewed from the primary.

Figure 3 .
Figure3.Here we display 4PCF measurements from the illustrative toy boxes.The left and right panels each show a different channel, as indicated in their titles.In both panels, the "clockwise-only" box results are in blue, the "counterclockwise-only" box results are in brown, and the "mixed" box results are in grey.Left: Projection onto P 111 .As discussed in the main text, we expect a negative projection for the clockwise box and a positive projection for the counterclockwise box.The mixed will have on average zero projection.These expectations are borne out and indeed analytic calculation yields agreement with the measured signals.Right: Projection ontoP 122 (r 1 , r2 , r3 ) ∝ ir 1 • (r 2 × r3 )(r 2 • r3 ).By construction the scalar triple product r1 • (r 2 × r3 ) is unity for all three of our illustrative boxes.The amplitudes on the righthand side are not strictly zero because there are still a small number of tetrahedra formed by connecting secondaries around one primary with secondaries around another.Additionally, the amplitude fluctuations divide into two envelopes at bin index 36.This bin index corresponds to raising the index for r 1 (the smallest side and the slowest-varying one).

Figure 4 .
Figure4.The number density n as a function of redshift z for the two BOSS samples used in this work.The LOWZ (0.2 < z < 0.4) North Galactic Cap (NGC) is brown and South Galactic Cap (SGC) is orange.For CMASS (0.43 < z < 0.7), the NGC is in purple and the SGC is in lavender.We intentionally do not allow redshift overlap between the two samples, and the redshift gap ∆z = 0.03 corresponds to a comoving radial separation of 73 h −1 Mpc.This separation means that the samples are fairly independent; the 2PCF ξ between a point in LOWZ and CMASS would be of order 1% at this scale.The covariance between two "worst-case" tetrahedra, one in each sample (where each is very close to the respective edge), is of order ξ 4 .Few tetrahedra are near enough on either edge to be significantly correlated with those in the other slice.The plot shows that LOWZ has both a more uniform selection function and a somewhat higher average number density than CMASS.It also lacks the strongly decaying tail with an increasing redshift that CMASS displays.These points are important when assessing the possible impact of systematics on each sample and when addressing any differences between the detection significances in the two samples.

Figure 5 .
Figure5.Upper: The parity-odd 4PCF for the BOSS CMASS data with ten radial bins, with NGC in red and SGC in blue, showing only the six lowest-lying channels.The error bars are the root-mean-square of the Patchy mocks (over our set of 2,000).Here for legibility we focus on the ten-bin results; the eighteen-bin results are in Appendix G, as are the remaining channels for the ten-bin analysis.We have mapped the three radial bins to a single index, with r 3 varying the fastest and r 1 the slowest.This is done in all similar plots that follow.Lower: The mean of 2, 000 Patchy mocks (solid curves) for both NGC (red) and SGC (blue).The shaded region is the rms expected for a single mock; this set of panels is intended to display the region around zero that the 4PCF of a dataset with no true parity-odd signal might inhabit.

Figure 7 .
Figure7.Comparisons of the analytic parity-odd 4PCF covariance with that from the Patchy mocks, all for CMASS NGC and ten radial bins.For greater legibility, we show only 11 of the 23 channels used in our analysis, and we have mapped each unprimed combination of angular momenta and radial bins to a 1D index, and the same for each primed; the full set of variables and indices involved is given in Eq. (13), and the volume and number density here used are given in Table1.Left: analytic correlation matrix (upper triangle) and correlation matrix from the Patchy mocks (lower triangle).The sub-blocks represent the covariance between a pair of channels; each is 120 × 120 as for ten radial bins, there are 120 radial bin combinations in each channel.These sub-blocks show how the covariance changes as the side lengths of each tetrahedron vary at fixed channel pairs.The similarity between the upper and lower triangles shows that the analytic covariance captures the cross-covariance between different radial bins as well as between different channels.Right: diagonals of the Patchy covariance matrix (solid black) and the analytic covariance (dashed red) for CMASS NGC.The average ratio is 1.05.

Figure 8 .
Figure 8.Here we show the same comparisons as made in Fig. 7 but for the LOWZ NGC.Again, in the left panel, we see that the analytic and mock correlation matrices have very similar patterns.The average ratio between the mock and analytic covariance diagonals is 1.01 as indicated by the dotted black horizontal line in the lower right panel.

Figure 9 .
Figure 9. Half-inverse test Φ (Eq.16) of the analytic covariance matrix at {Λ, Λ } = {111, 111} for CMASS NGC (left) and LOWZ NGC(right).We show only half of the matrix in the lower triangle and indicate the standard deviation of the matrix elements in the upper triangle; the test is symmetric so there is no information lost by deleting the latter.Here we compare the analytic covariance matrix to that from 2, 000 Patchy mocks; the overall standard deviation is roughly 1/ √ N mocks ≈ 0.02.It is expected that the standard deviation of the diagonal is double, as further discussed in the main text; we see a slight deviation from this relation in both panels.

Figure 11 .
Figure11.Similar to Fig.10but for LOWZ; again using ten radial bins.As for CMASS, in the compressed approach, the detection significance first increases with an increasing number of eigenvalues and then drops.In general, the detection significance is comparable to CMASS for both compressed and direct analyses and for split and joint sky.Again in the direct method, we see a noticeable deviation of the mocks' histogram from a χ 2 distribution; this again suggests that the analytic covariance may not be unbiased.We also see noticeable deviation at N eig = 800 (and even at N eig = 500), to a much greater extent than for the same N eig values in CMASS (cf.Fig.10).

Figure 12 .
Figure12.Similar to Fig.10but for CMASS with eighteen bins.Using finer binning we find even higher detection significance, both in our compressed and direct analyses (cf.Fig.10).Unlike the ten-bin CMASS analysis, which peaked at N eig = 200, here the detection significance monotonically increases with N eig .This behavior is likely due to how the covariance matrix structure and signal-to-noise in each channel and bin combination balance with the penalty paid for adding more degrees of freedom if they do not contain much additional S/N.Given that this balance can change depending on the radial binning, we would not have expected that the behavior of the ten-bin and eighteen-bin compressed analyses with N eig would be uniform.

Figure 13 .
Figure13.Similar to Fig.11but for LOWZ with eighteen bins.As is the case for CMASS, the detection significance increases when using the finer binning (cf.Fig.11).However, the increase is more moderate than the increase for CMASS.The detection significance rises monotonically with N eig , in agreement with the trend for the ten-bin LOWZ analysis of Fig.11.
Figure14.Reduced max analysis for ten radial bins.We use three different methods of obtaining the covariance matrix but always with max restricted to 1 or 2. The left column uses the compression scheme.The central column uses the covariance matrix obtained directly from the mocks, which is possible for this reduced max analysis.The third column uses the covariance obtained from adjusting the analytical covariance to fit, as well as possible, the purely mock-based covariance.The solid curves show the expected distributions for T 2 or χ 2 while the histograms show the distributions obtained from the mocks.The solid lines indicate results from the data.There is overall good agreement in detection significance between the purely mock-based covariance (central column) and the analytic covariance (right column).
Figure15.Left: Scatter plot for the 4PCF with ten bins measured from CMASS in the decorrelated basis (described in §5.2) and scaled by the covariance matrix eigenvalues, with both Pearson (rp) and Spearman (rs) correlation coefficients given in the legend.The extended panes in the left panel display the histogram overplotted with a Gaussian fit.Right: PDF f (r) (black curve; Eq. 22) for the Pearson correlation coefficient and the Pearson and Spearman correlation coefficients obtained from Patchy (orange and blue histograms).The 95% confidence level (dashed vertical lines) corresponds to an rp = 0.037.Our Pearson rp of -0.014 corresponds to a p-value of 0.460, so there is no statistically significant correlation detected.

Figure 16 .
Figure16.Same as Fig.15but for LOWZ, again with ten bins.Again we detect no statistically significant correlation.

Figure 18 .
Figure18.Exploration of a linear model for the underlying signal (vs. the eigenvalue index) for both CMASS and LOWZ (with ten radial bins).Left: Toy models for CMASS with m = 5.3 × 10 −8 , b = 1.1 × 10 −4 (dashed red) and LOWZ (solid purple) with m = 1.5 × 10 −7 , b = 2.7 × 10 −4 (these are our best-fit values; see Fig.19for discussion of the fitting and likelihood).Middle: Toy model (black curve; "fake signal") transformed back to the observed basis plotted with the measured CMASS 4PCF for 1 = 2 = 3 ; upper row is NGC and lower is SGC.Right: Analog of the middle panel but for LOWZ.

Figure 20 .
Figure20.The Pearson correlation coefficients rp (left column) and χ 2 values (right column) observed in our data (black lines) shown over the distribution of the same from the "fake signal" generated by our "linear-in-eigenbasis" toy model.We show results from m = b = 0 in blue, from the best-fit m and b in red (Fig.18displays these models and gives the best-fit values in the caption, and the likelihood used is shown in Fig.19), and from a larger "fake" signal in purple to indicate the possible range of outcomes.The fake data at each index (see Fig.18) is drawn from a Gaussian distribution with mean given by the linear model and width given by the square-root of the eigenvalue of the analytic covariance matrices.CMASS results are shown in the upper row and LOWZ in the lower.We see that there exists a "fake signal" that does explain the observed rp and χ 2 from the real data.While our "fake signal" should not be taken too literally, this does indicate that, generically, our observed χ 2 and rp are not necessarily in tension with each other.

Figure 23 .
Figure23.PDF of the radial selection function for Patchy NGC.The selection function applied to the Patchy mocks is designed to match that from the BOSS data for both the mock data catalog (black histogram) and the random catalog (grey histogram).To test whether an imperfect inference of the radial selection function can lead to a spurious parity-odd signal, we deliberately distorted the n(z) for the random catalog (red histogram), by shifting it by 10% on average, where the percentage shift is drawn from a Gaussian with that mean and width of 5%.We show the distribution of T 2 values for Patchy mocks analyzed using these "shuffled" randoms in Fig.24; it is consistent with no detection of a parity-odd 4PCF.

Figure 24 .
Figure24.Histograms of statistics for Patchy mocks' 4PCFs analyzed with different weighting schemes.Left: T 2 values calculated with 50 eigenvalues using our data compression scheme ( §5.1.2) with the covariance matrix inferred from the systematics-integrated mocks themselves.The black curve is the expected T 2 distribution.The histogram for the Patchy mocks with the standard weighting scheme is in black, with the close pair weight turned off to assess the impact of fiber collisions in blue, with the veto mask turned off to assess the impact of the angular mask in orange, and finally, in red, the results of shuffling the randoms' redshift distribution by 10% as described in Fig.23.None of these alterations to the 4PCF analysis produce a distribution of T 2 values that are inconsistent with that expected under the null hypothesis of no-parity-odd 4PCF (black curve).Right: χ 2 values calculated using the analytic covariance calibrated to mocks analyzed with the standard weights.Here we can see that the χ 2 distribution has a higher mean for the shuffled mocks than for the standard mocks, suggesting that although such a systematic does not induce a parity-odd signal, it might lead to underestimation of the covariance.

Figure 25 .
Figure 25.Left: Footprint of CMASS NGC showing how we split the data into eight patches, each roughly 40 • in right ascension and 35• in declination."SV" stands for sub-volume, as each patch extends into the redshift direction.Right: For six parity-odd channels, we display the difference between the 4PCF coefficients of the eight sub-volumes created by extending these angular patches in the redshift direction, and those from the full NGC sky.We divide these differences by the standard deviation of the Patchy results.The colors for the eight sub-volumes match in the left and right panels.Overall, we see that no one angular region stands out as producing a great deal more parity-odd signal than another.We offer a more rigorous demonstration of this claim in Fig.26, which shows the χ 2 for each sub-volume.

Figure 26 .
Figure26.Left: Here we show the T 2 values from our angular patch test as described in Fig.25.The T 2 value for the entire CMASS NGC is the thick black line, the eight subvolumes' T 2 values are in thin, colored lines, and the combined T 2 from the eight subvolumes are in dashed black.The predicted T 2 distribution (under the null hypothesis of no parity-odd signal) with 100 degrees of freedom is the black curve; we have used the compressed analysis described in §5.1.2.Right: Covariance for the T 2 values of the eight sub-volumes, estimated using the Patchy mocks.This covariance is needed to properly co-add the T 2 values for the eight sub-volumes to report the combined T 2 given in the lefthand panel (dashed line there).The slight non-vanishing off-diagonal indicates correlations between the eight sub-volumes, but this turns out not to affect substantially the combined T 2 value.We may observe that the combined T 2 value differs from that for the whole NGC by about 10%; we believe this is due to missing the produced by quartets that are split across more than one sub-volume, which go uncounted in our sub-volume analysis.We would expect that adding them in (as the whole NGC analysis does) would indeed raise the χ 2 .

Figure 28 .
Figure 28.Left: Edge-correction matrix from Patchy mocks with CMASS NGC geometry at a given radial bin (r 1 , r2 , r3 ) = (27, 55, 97) [h −1 Mpc].Each Λ = { 1 , 2 , 3 } denotes a specific angular channel for both even and odd parity.The non-vanishing offdiagonal features indicate that the survey geometry couples even and odd modes.Right:The odd coefficients of the randoms normalized by R 000 .We have mapped the odd angular momentum triplets given by Λ to a 1D index (vertical axis), and show the radial bin also so mapped on the horizontal axis.Here we see the dimensionless parity-even-to-parity-odd mixing is of order 10 −6 , showing that the dominant change to the odd-parity modes from edge-correction is due to R 000 , i.e. overall renormalization from dividing by the mean number density.

Figure D1 .
Figure D1.Redshift-distribution dependence on imaging depth for CMASS NGC.Left: Normalized galaxy number counts as a function of redshift for i-band.We have split the sample into three bins in imaging depth.Right: Same as the left but for r-band.We see that the three bins in i-and r-band depths have very similar n(z), implying that the imaging depth is unlikely to impact our analysis.

Figure G2 .
Figure G2.The parity-odd 4PCF for the BOSS LOWZ data including all the angular channels for ten radial bins.NGC is in brown and SGC in blue; the error bars are the rms of the Patchy mocks.

Figure G4 .
Figure G4.The parity-odd 4PCF for the BOSS LOWZ data including all the angular channels for eighteen radial bins.NGC is in brown and SGC in blue; the error bars are the rms of the Patchy mocks.
for comparison.

Table 3 .
(six bin results, CMASS  only).Correlation coefficients between NGC and SGC for CMASS and LOWZ samples; see §5.2 and Figs.15-17.We do not see significant cross-correlation between NGC and SGC in either binning scheme or sample; the possible reason is further discussed in the main text, as is te reason we cannot perform a CMASS × LOWZ analysis.