\section{Introduction}
Baryon Acoustic Oscillations (BAO) are a signal encoded in the two-point statistics of the CMB temperature and matter distributions and provide a standard ruler \citep{Blake2003ProbingRuler, Seo2005BaryonicSurveys} which can be used to trace the evolution of our Universe back to the recombination epoch, about $\num{380000}~\mathrm{yr}$ after Big-Bang. This BAO signal was first observed on galaxy distributions by \citet{Eisenstein2005} on Sloan Digital Sky Survey (SDSS) data by estimating the two-point correlation function (2PCF) of the galaxy sample and by analyzing the power spectrum of the 2df survey galaxy sample \citep{Cole2005TheImplications}. Since then, more observations have improved or complemented the early measurement such as the Baryon Oscillation Sky Survey (BOSS) \citep{Alam2017TheSample}, the extended BOSS (eBOSS) \citep{2020arXiv200708991E}. The ongoing Dark Energy Spectroscopic Instrument (DESI) survey \citep{Levi2019TheDESI},  will soon provide larger amounts of data than ever before. 

A variety of tracers can be used to estimate the matter density observationally. Their particular physical properties induce a bias in their clustering with respect to the underlying dark matter distribution. \citet{Gil-Marin2010ReducingStochasticity} have shown that a multi tracer approach can yield more information than a single tracer analysis. In particular, cosmic voids can be regarded as a secondary tracer that can be extracted directly from an existing tracer sample without further observations. 

In general we define cosmic voids as underdensities in the matter field. In practice, however, the formal definition depends strongly on how voids are extracted from galaxy catalogs. Some methods rely on a continuous density field and extract voids of arbitrary shape by merging underdense regions of the smoothed field (e.g. \citet{Platen2007}) or imposing spherical voids (e.g. \citet{Colberg2005}). Other methods use purely geometrical arguments to define voids from a discrete point distribution. The \textsc{Zobov} \citep{Neyrinck2008} void finder, defines voids using the Voronoi tessellation of the point space. In the present work, we use the Delaunay Triangulation Void findEr (\textsc{Dive}) \citep{Zhao2016} which defines voids as a subset of the circumspheres of the simplices resulting from the Delaunay Triangulation [DT]\citep{Delaunay1934} of the point space. This definition in particular, ensures that there are no galaxies within the spheres and poses no restrictions on the overlapping of spheres. 

\citet{Zhao2016} identify two distinct sphere populations that can be distinguished by their radius $R$. Spheres with $R<R_c$ behave as voids-in-clouds \citep{Sheth2003} since they are located in matter-dense regions. Hence, they do not fit the ``underdensities in the field'' definition of a cosmic void. Large spheres ($R>R_c$) are located in predominantly underdense regions and behave as voids-in voids. These are considered cosmic voids in this work. As an example, \citet{Kitaura2016, Liang2016} have already observed the BAO signal in the sample of overlapping large spheres.

The possibility to have overlapping spheres greatly increases the ``void'' sample size and opens the possibility of BAO measurements that had been impossible on voids due to the limited number of objects \citep{Zhao2016}.  Moreover, it has been shown that the void number density depends on the galaxy density \citep{Zhao2016} which would imply that the two sphere populations behave differently under such change, that is, the boundary between them changes too.

Some studies have shown that the void distribution encodes higher-order information \citep{Kitaura2016, Zhao2016} of the galaxy distribution. Furthermore, a multi-tracer analysis combining galaxies and voids \citep{Zhao2020} has shown that the void statistics can contribute to a decrease in the BAO scale measurement error in most mock realizations. In addition, there is a $\sim 10\%$ improvement in the distance measurement for observational data obtained by the Sloan Digital Sky Survey (SDSS) corresponding to the BOSS DR12 LOW-Z sample.

Redshift surveys face various observational challenges which result in an incomplete observed sample with respect to the number of known targets. Examples of this include the difficulties in observing galaxies separated by an angular distance smaller than the diameter of the fiber ferule in the spectrograph (e.g. close pairs), the discarded objects due to errors in the redshift measurement pipeline (redshift failures) or due to angular systematics. The latter are related to factors such as the seeing and stellar density which mask a fraction of the known targets. These examples affect the sample in the angular plane and result in a relatively small incompleteness. Nonetheless, even small angular incompleteness may result in biased clustering measurements if not carefully corrected (see \citet{Ross2011, Ross2017a}). Besides these angular effects, there is further incompleteness that comes from the radial selection function (see for example fig. 11 in \citet{Reid2016}) which is accounted for in the analysis by means of the Feldman, Kaiser \& Peacock (FKP) weights \citep{Feldman1994}. For a complete discussion on the systematics of the BOSS DR12 sample see \citet{Reid2016}.

In this paper we aim to study in more detail how the incompleteness in the tracer sample induced by observational systematics and radial selection affects the DT void sample and its clustering measurements. Moreover, we estimate how these may affect cosmological measurements involving voids.


In the next section, we describe the data used in this work. In \refsec{sec:gen-methods} we show the general methods we used for the analysis such as the definition of the SNR, 2PCF estimation, definitions for the combined analysis, BAO reconstruction and BAO fitting. \refsec{sec:radius-cut-dependence} explores the dependence of the $\mathrm{SNR}$-maximizing void radius cut on the tracer density. In \refsec{sec:bao-constraints} we analyze the selection effects on the different void samples defined in terms of BAO measurements. Finally, we conclude in \refsec{sec:conclusions} and outline future work.
\section{Data}
\label{sec:data}
Throughout this study we use two different data sets for cubic simulations and lightcone mocks. In addition we apply incompleteness masks to the cubic simulations in order to generate data sets that simulate angular systematics and radial selection on boxes. For all data, the DT sphere sample is extracted using the \textsc{DIVE} code or equivalently, its Python version \textsc{PyDIVE}. This section introduces the galaxy catalogs used.

\subsection{Simulation boxes}
We use 500 realizations of mocks produced using the PerturbAtion Theory Catalog generator of Halo and galaxY distributions (\textsc{Patchy}{}) code \citep{Kitaura2014}. The catalogs are produced using Augmented Lagrangian Perturbation Theory and are calibrated using the BigMultiDark (MD) n-body simulations \citet{Klypin2016}. These are the snapshots used for generating BOSS DR12 MD-\textsc{Patchy}{} mocks. We will refer to this data as \textsc{Patchy} boxes. The boxes have a volume of $L^3=\num{2500}^3~~h^{-1}\mathrm{Mpc}^3$ and have a galaxy density of  $\bar{n}_{\mathrm{gal}}\approx \num{4e-4}~h^3\mathrm{Mpc}^{-3}$.

The fiducial cosmology used for the analysis is \citep{Kitaura2016}
\begin{align}
    \label{eq:fiducial-om}
    &\Omega_{\rm m} = \num{0.307115},\\
    &\Omega_{\rm b} = \num{0.048206},\\
    &\sigma_8 = \num{0.8288},\\
    &n_{\rm s} = 0.9611,\\
    \label{eq:fiducial-h}
    &h = 0.677.
\end{align}
For \textsc{Patchy} boxes, the radial direction corresponds to the cartesian $Z$ coordinate. We also generate \textsc{Patchy}{} boxes in redshift space by applying Redshift Space Distortions (RSD) in this direction. The angular plane corresponds then to $(X,\ Y)$.
\subsection{Simulating incompleteness in boxes}
Based on the \textsc{Patchy}{} boxes we create incomplete data via random downsampling of the complete point distribution. There are three downsampling schemes: uniform, angular and radial. The respective data sets are dubbed \textsc{Patchy}{}-U, \textsc{Patchy}{}-NUA and \textsc{Patchy}{}-NUR. To avoid confusion, we will refer to the full \textsc{Patchy}{} boxes (without downsampling) as \textsc{Patchy}{}-F.
\subsubsection{Uniformly downsampled boxes: \textsc{Patchy}{}-U}
This sample is created by downsampling the \textsc{Patchy}{} boxes with a constant completeness $C$ in intervals of $0.05$, thus producing boxes with galaxy densities of $\bar{n}'_\mathrm{gal} = C \bar{n}_\mathrm{gal}$ and with the same clustering amplitude. Under these conditions, we are able to study the dependence of the void sample purely on the tracer density.

\reffig{fig:void-density-vs-radius} shows the void number density distributions as a function of the void radius (top panel) for the different galaxy densities in \textsc{Patchy}{}-U. The radius distribution functions shift towards larger voids as the galaxy density decreases, which is consistent with results shown in \citet{Zhao2016}. This confirms that the radius cut depends strongly on the galaxy density as well. This dependence could be parametrized as a power law, inspired by the definition of a dimensionless scale $\bar{n}_{\mathrm{gal}}^{1/3}R$. The bottom panel of the same figure shows that the distributions span a similar value range in this scale and the median void radii (black dotted lines) align at $R = \qty(0.8775\pm0.004)~\bar{n}_{\mathrm{gal}}^{-1/3}$. Previous work by \citet{zhao2017} also finds this scale relevant by showing that the DT sphere radius, $R_1$, at which spheres have a bias of 1, is precisely $R_1 = \bar{n}_{\mathrm{gal}}^{-1/3}$.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/redshift_both_scales_post.png}
    \caption{Void number distribution functions as a function of void radius for mocks of different galaxy densities. The top panel shows the distributions as a function of the radius cut $R$. The middle panel shows that in the dimensionless units $\bar{n}_{\mathrm{gal}}^{1/3}R$ the distributions span a similar range of dimensionless radii. The vertical dotted lines show the median void radius, the near-perfect alignment of the lines makes it difficult to distinguish them individually. Bottom panel shows the analogous plot to the mid panel for reconstructed boxes. Notice that in this case the mean values do not align as well as they do before reconstruction.}
    \label{fig:void-density-vs-radius}
\end{figure}
\subsubsection{Non-Uniform Angular downsampling: \textsc{Patchy}{}-NUA}
\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{plots/density_maps.png}
    \caption{Tracer densities for boxes with simulated incompleteness. Color maps show the density field $\delta+1 = \rho/\bar{\rho}$ for galaxies (top left) with angular downsampling and different radius bins of the DT spheres extracted from it. The last panel shows the radial density ($Z$ direction) for galaxies downsampled in this direction (dashed line) and the DT spheres extracted from it in different radius bins (continuous lines). The colorbar shows the radius  bin corresponding to each curve.}
    \label{fig:box-incompleteness}
\end{figure*}
Angular systematics on observations can be broadly described as sample incompleteness in the angular plane. Recent cosmological clustering analysis focus on regions of the sky with very high completeness, such that the sample incompleteness is very low. In this work we simulate this effect by downsampling the \textsc{Patchy}{} boxes in the $X,~Y-$plane with a 2D paraboloid mask whose maximum is placed in the center of the box. The minimum completeness is set to 80\% at the corners of the box. \reffig{fig:box-incompleteness} shows the $X,~Y$ map of the measured galaxy density (top-left panel) and DT sphere densities for different radii. It should be noted that smaller spheres show distributions closely resembling that of the galaxies of \textsc{Patchy}{}-NUA, confirming that these are correlated tracers. As the sphere radius increases the distribution becomes very uniform (bottom-left panel) and for even larger radii, the distribution becomes anti correlated to that of the galaxies (bottom middle-left panel).
\subsubsection{Non-Uniform Radial downsampling: \textsc{Patchy}{}-NUR}
The \textsc{Patchy}{}-NUR sample is obtained by downsampling the \textsc{Patchy}{}  boxes following a Gaussian completeness
\begin{equation}
\label{eq:radialincompleteness}
C(Z)=\frac{n'_{\mathrm{gal}}}{\bar{n}_{\mathrm{gal}}}(Z) = \exp{-\frac{(Z-c_ZL)^2}{2(\sigma_ZL)^2}};
\end{equation}
The parameters $\quad c_Z=0.5,~\sigma_Z=0.235$ define respectively the center and width of the distribution and $L$ is the box size. The value of $\sigma_Z$ is chosen in such a way that at the box boundaries, $10\%$ of the galaxies is kept. This distribution is meant to simulate a realistic radial selection function. The bottom-right panel of \reffig{fig:box-incompleteness} shows the measured radial density from \textsc{Patchy}{}-NUR boxes (dashed line) and the corresponding void density for different radius bins (solid lines) which correspond to the bins selected for the color maps. Analogous to the angular maps, the smaller DT spheres follow the galaxy distribution and flatten as the radius increases, then the distribution is inverted for the largest spheres.

\subsection{Light-cones}
In order to test our observations in the simulation boxes, we also analyze lightcone data, composed by the BOSS DR12 Luminous Red Galaxy (LRG) sample and the corresponding mock catalogs, the \textsc{MD-Patchy DR12} data. The BOSS DR12 LRG sample is composed by 1.2 million objects spanning $\num{9300}~\rm deg^2$ and redshifts $z\in(0.2,0.75)$ \citep{Ross2017a} and is the results of observations done from the Apache Point Observatory up to phase III of SDSS. Moreover, the corresponding mock catalogs are created by calibrating the Halo Abundance Matching (HAM) to match BigMD with BOSS DR12 LRG data \citep{Kitaura2016}. 

For our analysis, we use the combined (LOWZ and CMASS) sample and divide it in two redshift bins \citep{Alam2017TheSample, Zhao2020}. A low-$z$ bin with $z\in(0.2, 0.5)$ and a high-$z$ bin with $z\in(0.5, 0.75)$. To correct for known systematical effects in the galaxy sample, we also use weighting schemes as thoroughly described \citet{Reid2016, Ross2017a}, including FKP \citep{Feldman1994} weights to account for the radial dependence of the galaxy density. \reffig{fig:nz-caps} shows the radial distribution of galaxies on lightcones.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/nz_caps.png}
    \caption{Radial number densities for BOSS DR12 data in the North (NGC) and South (SGC) galactic caps. Measured from \textsc{Patchy}{} mocks.}
    \label{fig:nz-caps}
\end{figure} 

We use 2048 pre-reconstruction \textsc{MD-Patchy DR12} mocks and 1000 reconstructed lightcones to estimate the error on the observations.

To perform BAO fits on lightcones we use the fiducial cosmology $h=0.676$ and $\Omega_{\rm m}=0.31,\ \Omega_{\rm b} = \num{0.048},\ \sigma_8 = \num{0.8},\ n_{\rm s} = 0.97$.


\section{Methods}
\label{sec:gen-methods}
\subsection{BAO Reconstruction}

We use the BAO reconstruction method \citep{Eisenstein2007a} to analyze both cubic simulations and lightcone data. We employ the iterative Fourier-space method described in \citet{Burden2015} that we briefly describe here. 

First, using the Zeldovich approximation \citep{Zeldovich1970}, one can obtain the displacement field $\vb*{\Psi}(\vb*{k})$ from the linear real-space overdensity field $\delta_r(\vb*{x})$ by solving $\nabla\cdot \vb*{\Psi}(\vb*{q}) = -\delta_r(\vb*{x})/b$ where $\vb*{x} = \vb*{q} + \vb*{\Psi}(\vb*{q})$. By Fourier transforming the equation, it is evident that up to a first order
\begin{equation}
    \vb*{\Psi}(\vb*{k}) = i\frac{\vb*{k}}{k^2}\frac{\delta_r(\vb*{k})}{b},
\end{equation}
where $b$ is the galaxy bias, and $\vb*{\Psi(\vb*{k}})$ and $\delta_r(\vb*{k})$ are the Fourier-transformed displacement and real-space overdensity fields. However, to obtain the displacement field from a redshift-space overdensity, the process must account for RSD. The equation to solve is \citep{Nusser1993OnSamples}
\begin{equation} 
    \label{eq:displacement-rsd}
    \nabla\cdot\vb*{\Psi} + f\nabla\cdot(\vb*{\Psi}\cdot\vb*{\hat{r}})\vb*{\hat{r}} = -\frac{\delta_s}{b},
\end{equation}
with $f$ the growth rate, $\vb*{r}$ the radial coordinate vector, $\vb*{\hat{r}}$ the unitary vector in the same direction and $\delta_s$
is the redshift-space overdensity. This equation is solved for a potential $\phi$ such that $\vb*{\Psi} = \nabla\phi$. A first approximation can be obtained by replacing $\delta_r\rightarrow\delta_s$ in the Zeldovich approximation to estimate 
\begin{equation}
    \phi^{(1)}(\vb*{k}) = \frac{\delta_s}{k^2b} \equiv \frac{\delta_r^{(0)}}{k^2b},
\end{equation}
which can now be used to compute an estimate of $\delta_r^{(1)}$ using \refeq{eq:displacement-rsd}
\begin{equation}
    -\nabla\cdot\vb*{\Psi} = \frac{\delta_r^{(1)}}{b} =  f\nabla\cdot(\nabla\phi^{(1)}\cdot\vb*{\hat{r}})\vb*{\hat{r}}  + \frac{\delta_s}{b}.
\end{equation}
Given the estimate of $\delta_r^{(1)}$, it can be used in the Zeldovich approximation once again to compute $\phi^{(2)}(\vb*{k})$ which in turn is used to compute $\delta_r^{(2)}$. The process can be iteratively applied to obtain a good approximation to the real $\delta_r$. Once this is obtained, the Zeldovich approximation gives the displacement field necessary to perform the reconstruction.

In this work we use the \textsc{Revolver}\footnote{\url{https://github.com/seshnadathur/Revolver}} \citep{Nadathur2019REVOLVER:Reconstruction} implementation of this algorithm with 3 iterations. For simulation boxes, we use $b=2.2$ for the bias of the LRG sample and $f=0.743$, corresponding to the growth rate at redshift $z=0.466$ (the redshift of the boxes). For lightcones we use the same bias value but $f(z=0.5)=0.757$. The overdensity field is put in a grid of size $512^3$ using the Cloud-In-Cell (CIC) method. A smoothing scale of $15~h^{-1}\mathrm{Mpc}$ is used to avoid non-linearities in the overdensity field.


\subsection{2PCF estimation}
\label{sec:2pcf-estimation}
When working with \textsc{Patchy}{} boxes we use the natural estimator \citep{Peebles1974}
\begin{equation}
    \xi^{\rm ab}(s) = \frac{D_{\rm a}D_{\rm b}}{RR} - 1,
\end{equation}
where $RR(s)$ is analytically computed and the $D_{\rm a}D_{\rm b}$ are the normalized pair counts of data catalogs of tracers a and b. Nonetheless, the \textsc{Patchy}{}-NUA, \textsc{Patchy}{}-NUR and \textsc{MD-Patchy} are incomplete, which makes it necessary to use random catalogs to more accurately sample the domain geometry. In this cases we use the Landy-Szalay estimator \citep{Landy1993, Szapudi1997ACorrelations, Padmanabhan2012}
\begin{equation}
    \xi^{\rm ab}(s) = \frac{D_{\rm a}D_{\rm b} - D_{\rm a}R_{\rm b} - D_{\rm b}R_{\rm a} + R_{\rm a}R_{\rm b}}{R_{\rm a}R_{\rm b}}.
\end{equation}
This estimator is also used for reconstructed data, in which case we replace the galaxy random catalogs $R_{\rm g}$ with the shifted random $S_{\rm g}$.

For lightcone data we will combine the void and galaxy samples as introduced in \citet{Zhao2020}. To do this we use a modified LS estimator
\begin{equation}
    \xi_{\mathrm{comb}} = \frac{n_{\rm g}^2\xi_{\rm gg}R_{\rm g}R_{\rm g} + 2n_{\rm g}n_{\rm v}w\xi_{\rm gv}R_{\rm g}R_{\rm v}+n_{\rm v}^2w^2\xi_{\rm vv}R_{\rm v}R_{\rm v}}{n_{\rm g}^2R_{\rm g}R_{\rm g} + 2n_{\rm g}n_{\rm v}wR_{\rm g}R_{\rm v}+n_{\rm v}^2w^2R_{\rm v}R_{\rm v}}.
\end{equation}
The quantities $\xi_{\rm gg}$, $\xi_{\rm vv}$ and $\xi_{\rm gv}$ correspond to the galaxy autocorrelation, void autocorrelation and galaxy-void cross correlation, and the $n_{\rm g,v}$ denote the number of galaxies and voids respectively.

To generate random catalogs for voids $R_{\rm v}$ we apply a slightly modified version of the algorithm presented by \citet{Liang2016}. More precisely, we change step (4) of the process: split each subsample in two parts and shuffle separately. Let the catalog have columns ($X,Y,Z,f(X,Y,Z)$) with $f$ a function of the coordinates. If the function is dependent on a single or a subset of the coordinates, then the column $f$ is to be shuffled alongside that coordinate(s) to avoid altering the dependence. For example, if $f = f(X,Y)$ then the columns should be split in $(X,Y,f(X,Y))$ and $(Z)$ subsets, the lines of each should be shuffled separately. This slight modification generalizes the algorithm in \citet{Liang2016} when we recognize that the void radius, $R$, can be regarded as a function of the tracer density (see \reffig{fig:void-density-vs-radius}), which is in turn most dependent on the redshift (see \ref{fig:nz-caps}), then in that case the sample is split in $(RA,\ DEC)$ and $(z, R(n_{\rm g}(z))$.

As a metric of the ``goodness'' of the BAO signal in the 2PCF, we use the BAO signal $S$ as defined in \citet{Liang2016}.  
\begin{equation}
S = \xi_0(s^{\mathrm{BAO}}) - \frac{\xi_0(s_1^{\mathrm{dl}}) + \xi_0(s_2^{\mathrm{dl}}) + \xi_0(s_1^{\mathrm{dr}}) +  \xi_0(s_2^{\mathrm{dr}})}{4},
\end{equation}
with $s^{\mathrm{BAO}} = 102.5$, $s_1^{\mathrm{dl}} = 82.5$, $s_2^{\mathrm{dl}}=87.5$, $s_1^{\mathrm{dr}} = 117.5$ and $s_2^{\mathrm{dr}}=122.5~h^{-1}\mathrm{Mpc}$. From this definition, we can extract the $\mathrm{SNR}$ as 
\begin{equation}
\mathrm{SNR}(R_c) = \frac{\langle S \rangle}{\sigma_S},
\end{equation}
where the statistics are computed over various mock realizations. The error reported in the $\mathrm{SNR}$ measurements is estimated using jackknife, i.e. the reported values are the average and standard deviation of 499 measurements of the $\mathrm{SNR}$, each done by leaving one of the realizations out.

Note that in the case of anticorrelated samples, the signal, $S$, will be negative so the $\mathrm{SNR}$ will also be. When optimizing this metric, we will consider its absolute value.
\subsection{BAO fitting}
\subsubsection{BAO model}
To fit the BAO feature we use the model proposed by \citet{Xu2012} for a given fiducial cosmology (Eq. (\ref{eq:fiducial-om}), (\ref{eq:fiducial-h})), $\xi_{\rm t}$. The model correlation function is
\begin{equation}
\xi_{\rm t}(s) = \int\frac{k^2\dd k}{2\pi^2}\frac{\sin ks}{ks}P_{\rm t}(k)\exp(-k^2a^2),
\end{equation}
where the parameter $a = 1~h^{-1}\mathrm{Mpc}$ \citep{Xu2012, Zhao2020} defines the damping. The function $P_{\rm t}(k)$ is the template power spectrum, obtained via the relation
\begin{equation}
P_{\rm t}(k) = \qty[P_{\mathrm{lin}}(k) - P_{\mathrm{nw}}(k)]\exp\qty(\frac{-\Sigma_{\mathrm{nl}}^2k^2}{2}) + P_{\mathrm{nw}}(k).
\label{eq:pkt}
\end{equation}
The linear power spectrum for lightcones is computed using CAMB \citep{Lewis2002CosmologicalApproach} while for boxes it corresponds to the input power spectrum used to create the simulations.

To fit our data we introduce the model 2PCF \citep{Zhao2020},
\begin{equation}
\xi_{\mathrm{model}}(s) \equiv B^2 \xi_{\rm t}(\alpha s) + A(s).
\label{eq:model2pcf}
\end{equation}
The \textit{nuisance parameters}, $a_i,\ i=1, 2, 3$ encoded in the polynomial
\begin{equation}
A(s) = \frac{a_1}{s^2} + \frac{a_2}{s} + a_3
\end{equation}
are meant to fit the broad shape of the 2PCF. The other free parameters in the model are then $B$, for normalization, $\Sigma_{\mathrm{nl}}$ to account for the dampening of the BAO peak; and $\alpha$, the dilation parameter.

Nonetheless, to fit voids it is necessary to modify this model. The ratio between the tracer and linear non-wiggled power spectra, $P_{\mathrm{tracer,nw}}(k)/P_{\mathrm{lin,nw}}(k)$, shows a nontrivial behavior in the case of voids, while for halos it is almost flat \citep[see Figure 7 in ][]{Zhao2020}.

The non-linear behavior of such quantity is modeled as an extra quadratic term in $k$ \citep[see also figure 10 in][]{Tegmark2006}. The modified template spectrum is, more generally \citep{Zhao2020}
\begin{equation}
\begin{split}
    P_{\rm t}^v(k) =& \qty{\qty[P_{\mathrm{lin}}(k) - P_{\mathrm{lin,nw}}(k)]\exp\qty(\frac{-\Sigma_{\mathrm{nl}}^2k^2}{2}) + P_{\mathrm{lin,nw}}(k)}\times\\&\frac{P_{\mathrm{t,nw}}(k)}{P_{\mathrm{lin,nw}}(k)},
\end{split}
\label{eq:pktmod}
\end{equation}
with
\begin{equation}
\label{eq:pktemplate}
\frac{P_{\mathrm{t,nw}}(k)}{P_{\mathrm{lin,nw}}(k)}=1+ck^2,
\end{equation}
where $c$ is a fit parameter.
\subsubsection{Parameter inference}
We infer the posterior distribution, $p(\Theta|X)$, for each parameter, $\Theta$:
\begin{equation}
p(\Theta|X) = \frac{p(X|\Theta)p(\Theta)}{p(X)},
\label{eq:bayes}
\end{equation}
where we denote $X$ the data. The distribution $p(X|\Theta)\equiv\mathcal{L}$ is the \textit{likelihood}
\begin{equation}
\mathcal{L}(\Theta) \equiv A\exp(-\frac{\chi(\Theta)^2}{2}),
\end{equation}
where $A$ is a normalization constant and $\chi^2(\Theta)$ is defined as
\begin{equation}
\chi^2(\theta) = \qty[\xi_{\mathrm{data}} - \xi_{\mathrm{model}}(\Theta)]^T\vb*{C}^{-1}\qty[\xi_{\mathrm{data}} - \xi_{\mathrm{model}}(\Theta)].
\end{equation}
Here, $\xi_{\mathrm{data}}$ is the measured correlation function and $\xi_{\mathrm{model}}$ is the model 2PCF in \refeq{eq:model2pcf}. The matrix $\vb*{C}^{-1}$ is the inverse unbiased covariance matrix estimated from the mocks \citep{Hartlap2007}
\begin{equation}
\vb*{C}^{-1} = \vb*{C}_s^{-1}\frac{N_m - N_{\mathrm{bins}} - 2}{N_m - 1}.
\label{eq:covmat}
\end{equation}
The matrix $\vb*{C}_s^{-1}$ is the sample covariance matrix, 
\begin{equation}
\vb*{C}_{s,\ ij} =  \frac{1}{N_m -1}\sum_{k=1}^{N_m}[\xi_k(s_i) - \bar{\xi}(s_i)][\xi_k(s_j) - \bar{\xi}(s_j)],
\end{equation}
$N_m$ is the number of mocks used, $N_{\mathrm{bins}}$ is the number of separation bins used for the fit, $\xi_k(s)$ is the 2PCF as obtained from the $k$-th catalog, $\bar{\xi}(s)$ is the average 2PCF over all mocks and the indices $i,\ j = 1,\ \dots,\ N_{\mathrm{bins}}$.

The evidence $p(X)\equiv\mathcal{Z}(X)$ can be computed as the normalization of the posterior
\begin{equation}
p(X)\equiv\mathcal{Z} = \int\dd\Theta p(\Theta)p(X|\Theta).
\end{equation}
It corresponds to the average likelihood over the prior distribution, $p(\Theta)$.

The \textsc{Multinest} algorithm \citep{Feroz2008} based on \textit{Nested Sampling} \citep{Skilling2004} is used to compute the evidence and the posteriors in this work. The reported values of $\alpha$ and $\sigma_\alpha$ are the median and 1$\sigma$ CI of the distribution as obtained from \texttt{pymultinest}'s Analyzer \citep{Buchner2014}.
From the 6 (7) free parameters in the galaxy (void) model, we only use \textsc{Multinest} to fit the non-nuisance parameters. The $a_i$ are fitted by a simple least squares procedure.

\subsubsection{Priors}
\label{sec:priors}
The prior distribution $p(\Theta)$ of parameter $\Theta$ depends entirely on our knowledge.

Following the results in \citet{Zhao2020}, we use the fitting range $s\in[60, 150]~~h^{-1}\mathrm{Mpc}$ with $\num{1000}$ live points and tolerance of $\num{0.1}$.


In the case of galaxies, we set the following priors
\begin{align}
&p(\Sigma_{\mathrm{nl}}) = \mathcal{U}(0, 25)\\
&p(B) = \mathcal{U}(0, 50)\\
&p(\alpha)=\mathcal{U}(0.8, 1.2)
\end{align}

For voids we choose the following priors
\begin{align}
&p(\Sigma_{\mathrm{nl}}) = \mathcal{U}(0, 25)\\
&p(B) = \mathcal{U}(0, 50)\\
&p(\alpha)=\mathcal{U}(0.8, 1.2)\\
&p(c) = \bar{\mathcal{N}}(-100, 900, 100)
\end{align}
Here, $\mathcal{U}(\Theta_{\mathrm{min}}, \Theta_{\mathrm{max}})$ denotes a uniform distribution between $\Theta_{\mathrm{min}}$ and $\Theta_{\mathrm{max}}$; $\mathcal{N}(\mu_\Theta, \sigma_\Theta)$ is a Gaussian distribution with mean $\mu_\Theta$ and width $\sigma_\Theta$ and $\bar{\mathcal{N}}(\Theta_{\mathrm{flat\ min}}, \Theta_{\mathrm{flat\ max}}, \sigma_\Theta)$ is a flat distribution in the range $[\Theta_{\mathrm{flat\ min}}, \Theta_{\mathrm{flat\ max}}]$ but with Gaussian tails of width $\sigma_\Theta$:
\begin{equation}
\begin{split}
\bar{\mathcal{N}}(\Theta_{\mathrm{flat\ min}}, &\Theta_{\mathrm{flat\ max}}, \sigma_\Theta) = \\
&\begin{cases}
A', \hspace{2.65cm}  \Theta\in[\Theta_{\mathrm{flat\ min}}, \Theta_{\mathrm{flat\ max}}]\\
A'\exp(-\frac{(\Theta-\Theta_{\mathrm{flat\ min}})^2}{2\sigma_\Theta^2}), \Theta<\Theta_{\mathrm{flat\ min}}\\
A'\exp(-\frac{(\Theta-\Theta_{\mathrm{flat\ max}})^2}{2\sigma_\Theta^2}), \Theta>\Theta_{\mathrm{flat\ max}},
\end{cases}
\end{split}
\end{equation}
where $A'$ is a normalization constant.
\section{Radius cut dependence on galaxy density}
\label{sec:radius-cut-dependence}
To infer the relation between the optimal radius cut, $R^*$, and the galaxy density, we use the redshift-space \textsc{Patchy}{}-U data. We often present radius values in dimensionless scale, in which case we also provide a reference dimensionful value for a BOSS-like density of $\bar{n}_{\mathrm{gal}}\approx\num{3.5e-4}~h^3\mathrm{Mpc} ^{-3}$.
\subsection{Pre-reconstruction}
In order to optimize the $\mathrm{SNR}$ with respect to the radius cut \citep{Liang2016}, we perform a search in the dimensionless void radius space. Values of $\bar{n}_{\mathrm{gal}}^{1/3}R_c \in [0.5, 1.33]$ ($R_c \in [7, 19] ~h^{-1}\mathrm{Mpc}$) are used to define the radius cut for large spheres ($R>R_c$) when computing the void 2PCF. For small spheres ($R<R_c$), we used cuts in $\bar{n}_{\mathrm{gal}}^{1/3}R_c \in [0.3, 0.9]$ ($R_c \in [4.2, 12] ~h^{-1}\mathrm{Mpc}$ ). We obtain the optimal radius cut
\begin{equation}
R_c^*(\bar{n}_{\mathrm{gal}}) = \argmax_{R_c}\left(\abs{\mathrm{SNR}(R_c, \bar{n}_{\mathrm{gal}})}\right)    
\end{equation}
for every curve of constant $\bar{n}_{\mathrm{gal}}$ with a parabolic fit to the points near the maximum measured $\abs{\mathrm{SNR}}$ value. We use the \textsc{iminuit} package\footnote{\url{https://github.com/scikit-hep/iminuit}} to estimate the error in the position of the maximum $\abs{\mathrm{SNR}}$.

To generalize our measurements, we fit a power-law model of the form 
\begin{equation}
    R^*_c = \gamma \bar{n}_{\mathrm{gal}}^\nu,
\end{equation}
where we have weighted the importance of our measurements with their measured $\mathrm{SNR}$ value, thus giving more importance to the higher-density measurements which are more reliable. We use again \textsc{iminuit} to estimate errors on these parameters.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{plots/snr_all_pre.png}
    \caption{Optimum radius cut for pre-reconstruction boxes of different densities. The coloured region corresponds to the BAO $\mathrm{SNR}(R_c, \bar{n}_{\mathrm{gal}})$ map. The stars show the measured $R_c^*$  as a function of the galaxy density and the dashed line shows the fit to the model $R^*_c = \gamma \bar{n}_{\mathrm{gal}}^\nu$. The top panels show results for large-spheres while the bottom panels show measurements for small-spheres. Left side panels correspond to autocorrelations and right side panels to cross-correlations. The dotted lines on the latter show the fit to the measurements on the left panels. }
    \label{fig:snr-all-r-n-pre}
\end{figure*}
\subsubsection{Large-sphere autocorrelation}
The results of the $\mathrm{SNR}$ map for void-void autocorrelations are shown in the top-left panel of \reffig{fig:snr-all-r-n-pre}. Our model fit yields
\begin{equation}
    R_c^* \approx (2.39\pm0.02)~\bar{n}_{\mathrm{gal}}^{-0.237\pm0.001} \quad(R_c^* = 15.76 \pm 0.18 ~h^{-1}\mathrm{Mpc}).
    \label{eq:pre-rdependence}
\end{equation}
The maximum $\mathrm{SNR}$ we obtain is $\sim 8$, smaller than the galaxy-galaxy correlation $\mathrm{SNR}$ for this case, which is closer to $\sim 9.5$. Notice that the $\mathrm{SNR}$ is most sensitive to radius cuts for the large sphere autocorrelation.


\subsubsection{Small-sphere autocorrelation}
In \reffig{fig:snr-all-r-n-pre} (bottom-left panel) we present the results of the $\mathrm{SNR}$ optimization for the small sphere sample. Our fit shows that for densities $\lesssim \num{e-4}~h^3~\mathrm{Mpc}^{-3}$ there is no BAO signal from small spheres. For larger densities, the measured $\mathrm{SNR} \lesssim 4$, meaning that the small sphere autocorrelation does not yield much more information. For small densities and large radius cuts we find a large negative $\mathrm{SNR}$, which shows that in this region there are negative correlations due to the inclusion of large-enough spheres which actually trace cosmic voids. Finally we observe that the $R^*_c$ for this sample has a very weak dependence on the galaxy density and is given by
\begin{equation}
    R^*_c \approx (8.60\pm0.18)~~\bar{n}_{\mathrm{gal}}^{-0.022\pm0.002}\quad(R_c^* = 10.25 \pm 0.27 ~h^{-1}\mathrm{Mpc}).
\end{equation}
\subsubsection{Small-sphere-galaxy cross correlation}
Even though the small sphere autocorrelation does not show a very significant signal, these objects do show a stronger cross correlation with the galaxy sample, with a maximum $\mathrm{SNR}_{\mathrm{max}} \sim 7$, close to that of the large spheres' autocorrelation. In the bottom-right panel of \reffig{fig:snr-all-r-n-pre} we observe that there is indeed a BAO signal for most of the density domain probed by our experiments. Our fit results show that the optimum small sphere sample to cross correlate with galaxies is defined differently than that of the autocorrelation fit. In this case we obtain that 
\begin{equation}
    R^*_c \approx (4.90\pm0.09)~\bar{n}_{\mathrm{gal}}^{-0.096\pm0.002}\quad(R_c^* = 10.52 \pm 0.26 ~h^{-1}\mathrm{Mpc}),
\end{equation}
with a much closer fit in all densities. Even though the dependence on the galaxy density is small with $\nu = -0.096\pm0.002$, it is more significant than the result of small-sphere autocorrelation.
\subsubsection{Large-sphere-galaxy cross correlation}
The top-right panel of \reffig{fig:snr-all-r-n-pre} shows the $\mathrm{SNR}$ map obtained from the cross correlations between galaxies and large spheres. In this case, the maximum $\abs{\mathrm{SNR}} \sim 12$ which is notably higher than that of the galaxy and large sphere autocorrelations. It is then the correlation that carries the most information, potentially outperforming the galaxy autocorrelation measurement.

Our model fitted to the whole density range yields
\begin{equation}
    R_c^* \approx (4.88\pm0.08)~\bar{n}_{\mathrm{gal}}^{-0.138\pm0.002}\quad(R_c^* = 14.63 \pm 0.33 ~h^{-1}\mathrm{Mpc}).
\end{equation}

Given that the galaxy-density dependence of the $\mathrm{SNR}$ is weak in the small-sphere correlations and the signal is not as strong as the large-sphere counterparts, we will hereon focus on the latter, specially the autocorrelation, which shows the strongest density dependence. Therefore, we will only be using the autocorrelation parameters in \refeq{eq:pre-rdependence} to define the void clustering sample to be used in auto and cross correlations. In what follows we will also refer to the sample of large DT spheres as voids.


\subsection{Post-reconstruction}
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{plots/snr_vv_large_post.png}
    \caption{Optimal radius cut as a function of $\bar{n}_{\mathrm{gal}}$ for reconstructed downsampled boxes. The line shows the power law fit to the points. The underlying map shows the SNR.}
    \label{fig:snr-r-n-post}
\end{figure}
We now use the algorithm described in \refsec{sec:gen-methods} to reconstruct the redshift-space \textsc{Patchy}{}-U data. In this case we obtain the map shown in \reffig{fig:snr-r-n-post}, which yields the relation
\begin{equation}
    R_c^* \approx (2.10\pm0.01) ~\bar{n}_{\mathrm{gal}}^{-0.242\pm0.001}\quad(R_c^* = 14.41 \pm 0.13 ~h^{-1}\mathrm{Mpc}).
    \label{eq:post-rdependence}
\end{equation}
In this case the $\gamma$ parameter is different, however the power parameter, $\nu\approx-0.24$ is similar to the one obtained from pre-reconstruction voids. \citep{zhao2017} has shown that the DT spheres with a radius $R_1 = \bar{n}_{\mathrm{gal}}^{-1/3}$ have a bias of 1 and that this relation does not evolve significantly with redshift. Moreover, the same work shows that spheres with $R=R_0\approx1.8\bar{n}_{\mathrm{gal}}^{-0.29}$ have a bias of 0. In this sense, we could interpret our power parameter $\nu\approx-0.24$ as defining a sphere population with negative bias, which is a desired property in the cosmic void sample \citep{Zhao2016}. The redshift evolution of our power parameter, $\nu$, is left for a future work.


\section{Clustering of Voids with non-uniform incompleteness}
In this section, we show how voids extracted from the \textsc{Patchy}{}-NUA, \textsc{Patchy}{}-NUR and lightcone data sets are clustered.
\subsection{Simulation boxes}
\subsubsection{The constant cut void sample}
We show the clustering results in \reffig{fig:tpcf-angular-avg-loc}. First, the void clustering from \textsc{Patchy}{} boxes (marked No syst.), defined with a radius cut of $R_c = 15.6~h^{-1}\mathrm{Mpc}$, are shown as the shaded area for reference. The error bars show the clustering of the void sample $\mathcal{V}_0 \equiv \qty{\mathtt{s}\in\mathcal{S} | R(\mathtt{s}) > R_c}$, that is, a unique and constant radius cut $R_c$, is defined for the whole box. This is the approach that has been used in previous studies using DT voids \citep{Zhao2016, Liang2016, Kitaura2016, Zhao2020}. In the present study we define $R_c = \gamma\bar{n}_{\mathrm{gal}}^{\nu}$ using the parameters shown in \refeq{eq:pre-rdependence}, which is equivalent to $R_c  = 15.5~h^{-1}\mathrm{Mpc}$ for \textsc{Patchy}{}-NUA and $R_c  = 17.5~h^{-1}\mathrm{Mpc}$ for \textsc{Patchy}{}-NUR. The use of different radius cuts is motivated by the strong dependence of the void radius with the galaxy density as shown in \reffig{fig:void-density-vs-radius}. While the difference between \textsc{Patchy}{} and \textsc{Patchy}{}-NUA is very small, \textsc{Patchy}{}-NUR has about half the density, which means that using the same cut as \textsc{Patchy}{}-F would likely include more small spheres, thus contaminating the void BAO signal.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/tpcf_void_gal_angular_avg-loc.png}
    \caption{Two-point correlation functions of large-spheres from artificially downsampled boxes simulating angular systematics. Top panel shows the galaxy clustering with (errorbars) and without (shaded area) downsampling. Mid panel shows void clustering. ``Loc. cut'' marks the use of a local radius cut and ``Cst. cut'', the use of a constant one. ``rand'' marks the use of random catalogs (random sample takes into account incompleteness), otherwise, analytically computed random pair counts (uniform random sample, no incompleteness) were used. Bottom panel shows the galaxy-void cross correlation. All curves but ``No. syst'' are correspond to incomplete boxes.}
    \label{fig:tpcf-angular-avg-loc}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/tpcf_void_gal_radial_avg-loc.png}
    \caption{Two-point correlation functions of large-spheres from artificially downsampled boxes simulating radial selection. Convention is the same as \reffig{fig:tpcf-angular-avg-loc}. $y$-axis is in symmetric log scale with linear thresholds of 50, 15 and 15 (top to bottom) to show the very large bias introduced to the 2PCF due to the high incompleteness.}
    \label{fig:tpcf-radial-avg-loc}
\end{figure}

The \textsc{Patchy}{}-NUA correlations (\reffig{fig:tpcf-angular-avg-loc}) show that a small incompleteness is enough to introduce a noticeable bias in the galaxy autocorrelation (top panel) when the incompleteness is not considered by the random sample. By contrast, the void autocorrelations remain unscathed when using a constant radius cut. Analogously, the significantly larger incompleteness introduced by the radial selection function also induces significant bias in the galaxy-galaxy correlation (top panel, \reffig{fig:tpcf-radial-avg-loc}), albeit on a much larger scale, owing to the significantly higher incompleteness of the \textsc{Patchy}{}-NUR sample. Contrary to to the angular case, the void-void and galaxy-void correlations also show an equally large bias when the random samples do not take into account the incompleteness. When they do, the correlations are much closer to \textsc{Patchy}{}-F's, as shown by the black error bars. Nonetheless, this correlation shows a significant difference from the \textsc{Patchy}{}-F boxes at low scales $\sim 50~h^{-1}\mathrm{Mpc}$. This is possibly due to the very low density in some regions of the box ($Z$-edges), which yields DT spheres with large radii, that present a void exclusion effect on larger scales. It should be noted that a similar behavior is expected of any sparse tracer density and does not seem to be linked to the fact that the incompleteness is in any particular direction. In fact, very large exclusion features are found in low-density samples such as eBOSS LRG \textcolor{blue}{(Zhao et al. in preparation)} and QSO.
\subsubsection{The local cut void sample}
We use the fitting formulae in \refeq{eq:pre-rdependence} to define a different void sample $\mathcal{V}_l \equiv \qty{\mathtt{s}\in\mathcal{S} | R(\mathtt{s}) > \gamma n_{\mathrm{gal}}(\vb*{x}_{\mathtt{s}})^\nu}$, where $\vb*{x}_{\mathtt{s}}$ denotes the sphere position. This sample will be often referred to as the local cut sample given that it depends on the DT sphere position. The clustering of these void samples is shown as a dashed-line-bounded region. For \textsc{Patchy}{}-NUA with analytical randoms, using the $\mathcal{V}_l$ sample suffers a similar bias to that seen in the galaxies. This is expected since we are introducing more precise information about the local density and the randoms should reflect this more precise selection too.  Nonetheless, no bias is observed when using the shuffled void randoms according to the algorithm mentioned in \refsec{sec:2pcf-estimation}.  

With the \textsc{Patchy}{}-NUR set, the local cut proposed dampens significantly the large exclusion feature observed with a constant cut and yields clustering more consistent with that of the complete \textsc{Patchy}{} data. This can be understood by considering the local cut as a more nuanced selection. For example, on the edges of the box, where the galaxy number density is the smallest, the average radius is  $~23~h^{-1}\mathrm{Mpc}$. The constant cut at $R_c = 17.5~h^{-1}\mathrm{Mpc}$ selects a large number of DT spheres in the region which yield lots of pairs with large scale exclusions (\reffig{fig:tpcf-radial-avg-loc}). However, including the local density information, the cut becomes $R_c\approx25~h^{-1}\mathrm{Mpc}$ in this region, successfully ignoring average-sized DT spheres that ``pollute'' the sample and making the exclusion effect less significant. On the other hand, this reasoning also allows us to understand why voids with a constant cut ($\mathcal{V}_0$) are less sensitive to incompleteness, thus showing a smaller large-scale bias in the 2PCF than the $\mathcal{V}_l$ sample. By using a constant cut we do not discard these average-sized spheres in poorly sampled regions, which appear to fill the gaps in the observed density. Indeed bottom left-most panel in \reffig{fig:box-incompleteness} shows how uniform the void field is for these radii. As a result, the overall void distribution might be well approximated by a uniform random field as seen for \textsc{Patchy}{}-NUA, this however is not enough to approximate the distribution of the void sample from \textsc{Patchy}{}-NUR. Performing a cut with local information will discard these objects so that the statistics are more vulnerable to their distribution thus inducing the bias when a uniform random field is used for normalization.

The bottom panels in \reffig{fig:tpcf-angular-avg-loc} and \ref{fig:tpcf-radial-avg-loc} show the cross correlations analogous to the void autocorrelations described above, which are all consistent with those results. It should be noted that no galaxy weights were used to compute these correlations. Curves marked with ``rand'' were obtained however by using both galaxy and void randoms, such that the tracer incompleteness is taken into account, while for the remaining curves the analytical $RR(s)$ counts were used, thus a uniform tracer distribution is assumed in this case. It is remarkable that the cross correlations are still unbiased in large scales when using the $\mathcal{V}_0$ sample. This, coupled to the realization that this correlation has stronger $\mathrm{SNR}$ than even the galaxy sample alone means that it could yield better BAO measurements than current analysis. This is explored in \textcolor{blue}{Zhao et al. (in preparation)}. 

Moreover, we have used \textsc{Patchy}{}-NUR data and the corresponding voids along with uniform randoms for galaxies and the shuffled randoms for voids (hence including incompleteness) to investigate how robust the local cut can be to vastly inadequate randoms for galaxies (see top panel in \reffig{fig:tpcf-radial-avg-loc}). It is evident then that despite the larger variance, the correlation function shows no bias, even using the local cut, which indicates that even if we are unable to correct for unknown systematics in the galaxy data, the cross correlation will be robust using the local cut and the void randoms obtained through the algorithm in \citet{Liang2016}.

We have observed that the constant cut is very robust to low sample incompleteness, like the one that can be found in the angular plane of a galaxy survey, and allows to recover the correlations that are expected from a complete sample even using suboptimal randoms for normalization. On the contrary, when the incompleteness is very high, the void sample suffers from the same systematics as galaxies if the geometry of the survey is not well accounted for in the randoms. Furthermore, even using appropriate randoms, the void-void correlation using the $\mathcal{V}_0$ sample shows a large exclusion peak due to the inclusion of very large spheres from undersampled regions. In this high incompleteness case, the local cut proposed here avoids these DT spheres and yields a correlation that is compatible with the complete-box one on all scales. We have also shown that the cross correlation is robust even to very large incompleteness in the galaxy sample as long as appropriate void randoms (shuffled) are used.

\subsection{Light cones}
\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{plots/void_tpcf_avg.png}
    \caption{Average void two-point correlation functions from mock lightcone data. The shaded area shows the clustering obtained using the void sample extracted with a constant cut ($\mathcal{V}_0$). In this case we have chosen $R_c=16~h^{-1}\mathrm{Mpc}$ to contrast with the results of \citet{Zhao2020}. The error bars show the clustering results using a local cut to define the void sample.}
    \label{fig:void-tpcf-avg}
\end{figure*}
 The situation in lightcones is very similar to what is observed in the \textsc{Patchy}{}-NUA(R) data. The radial selection function depicted in \reffig{fig:nz-caps} shows that the low-$z$ region has a more stable density than the high-$z$ bin, in which density drops significantly towards $z\rightarrow0.75$. This density variation is much larger than any angular incompleteness, so we will assume for lightcones that the latter is constant and approximate $n_{\mathrm{gal}}(\vb*{x}) \approx n_{\mathrm{gal}}(z)$.
 
 The 2PCFs in \reffig{fig:void-tpcf-avg} mirror the difference in the behavior of the density in both bins. The top-left panel shows void clustering results for the low-$z$ bin, where there are no significant differences between the local and constant cuts. Contrarily, the high-$z$ constant cut clustering in the bottom panel of the same figure shows the same low-scale feature observed in the \textsc{Patchy}{}-NUR 2PCF. Moreover, the use of a local cut (\refeq{eq:pre-rdependence}) succeeds in better selecting a void sample, discarding the large spheres in coming from low sampling regions that enhance the exclusion effect.
 
 When we instead use reconstructed lightcones and the model in \refeq{eq:post-rdependence} we obtain results consistent with pre-reconstruction clustering. It is also clear that the void exclusion in the high-$z$ 2PCF (bottom-right panel) using the constant cut is not as wide as without reconstruction, which makes it less intrusive at BAO scales. This already enhances the BAO peak SNR, as would be expected from BAO reconstruction. In fact, the void size distributions post-reconstruction are narrower (see \reffig{fig:void-density-vs-radius}, bottom panel), showing that the processing moves galaxies into previously large empty regions, thus explaining the narrower exclusion feature. Nonetheless, the use of the local cut is still able to decrease the void exclusion amplitude, signaling that the sample contains less of the very-large spheres.
 
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{plots/bossdr12_void_snr.png}
    \caption{$\mathrm{SNR}$ computed from \textsc{Patchy} lightcone voids using both the constant radius cut (squares) and the local cut (triangles, Eqs. \ref{eq:pre-rdependence}, \ref{eq:post-rdependence}). Solid markers denote the post-reconstruction results while open markers show pre-reconstruction. The error bars correspond to the jackknife 95\% confidence interval}
    \label{fig:snr-voids}
\end{figure}
 These effects are quantified in \reffig{fig:snr-voids}, where it is evident that the local cut yields higher $\mathrm{SNR}$ in both bins before and after BAO reconstruction. As expected, the largest improvement is found in the high-$z$ bin before reconstruction. After reconstruction, the improvement is not as noticeable, given that this process already supresses large-scale exclusion effects
 
 We observe then that even though the constant cut post reconstruction is already good enough at mitigating void-exclusion effects for this sample, thus the use of the local cut here is not necessary.
 
\section{BAO constraints with incomplete samples}
\label{sec:bao-constraints}
Large enough sample incompleteness has been seen to affect the void clustering in the void exclusion scales, which depend in turn on how the void sample is selected. In this section we will evaluate how this impacts the cosmological measurements on the \textsc{Patchy}{} samples.

We start by fitting the \textsc{Patchy}{}-NUA and \textsc{Patchy}{}-NUR datasets using the parabolic model introduced in \refsec{sec:gen-methods}. \reffig{fig:baostats-nu} shows the SNR measurements and BAO fit results of void autocorrelations. We observe in the top panel that the use of a local cut is advantageous in terms of $\mathrm{SNR}$ in the large incompleteness (radial selection) situation. The difference in SNR on the right side is ten times smaller.

On the same figure, the second row shows that the $\alpha$ parameter is recovered in all cases with negligible differences. There are some samples with inferred $\alpha\approx0.8, 1.2$, which corresponds to the edges of the prior. This is a known issue of this fitting method, and better methods are being developed \textcolor{blue}{(Variu et al.  in preparation)}. The issue is visible also on the third row, where we show the uncertainty on the estimation of $\alpha$, $\sigma_\alpha$. Nonetheless, the bulk of the distribution is very consistent among \textsc{Patchy}{} and \textsc{Patchy}{}-NUA(R). Similarly, the goodness-of-fit metrics $\chi^2$ and $\ln{\mathcal{Z}}$ do not show significant differences, with distributions compatible in all situations.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/baostats_void_cst-loc_samps.png}
    \caption{Void-only BAO fit results and SNR measurements. Top row shows the SNR measurements on \textsc{Patchy}{}-NUR (left) and \textsc{Patchy}{}-NUA (right) data. Palens below show BAO statistics such as the $\alpha$ parameter measurement, the error on it's estimation $\sigma_\alpha$, the goodness of fit $\chi^2/\mathrm{dof}$ with $\mathrm{dof}=11$ (number of degrees of freedom in the fit) and the Bayesian evidence expressed as $\ln{\mathcal{Z}}$. Distributions are obtained from fitting each individual realization in the data sets. The left half of the violin plots (hatched, dashed line) represents the distribution of the statistic obtained using a constant cut. The shaded distributions (solid black line) on the right corresponding to the analogous using a local cut. The violin with a continuous line (red) shows the BAO fit of complete boxes. Finally the horizontal lines are from top to bottom: maximum, median and minimum}
    \label{fig:baostats-nu}
\end{figure}

Nonetheless, the void autocorrelation is rarely fit on its own. The combined galaxy-void sample incorporates more information and helps mitigate the effect of the suboptimal BAO fitting method on the results. We now apply this method to lightcone data. First we optimize the clustering measurements by tuning the $w$ parameter of the modified LS estimator. To do this we explore various $w$ values and perform the BAO fit on the average of the mocks as well as evaluating the BAO $\mathrm{SNR}$ for each $w$.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{plots/bossdr12_comb_snr.png}
    \caption{$\mathrm{SNR}$ as a function of $w$ computed from lightcone $\xi_{\mathrm{comb}}$ for both high (continuous lines) and low (dashed lines) $z$ bins. Red curves show results from local radius cut and black curves from constant radius cut. The error bars correspond to the jackknife 95\% CI}
    \label{fig:snr-comb}
\end{figure}
\begin{figure}
   
    \includegraphics[width=\linewidth]{plots/bossdr12_comb_bao.png}
    \caption{$\sigma_\alpha$, $\chi^2$ and $\ln\mathcal{Z}$ obtained from the fit to the mean of 2048 (1000) pre (post) reconstruction 2PCF as a function of $w$ computed from lightcone $\xi_{\mathrm{comb}}$. Line and color convention as in \reffig{fig:snr-comb}. Top row displays pre-reconstruction and bottom row shows post-reconstruction results. The horizontal lines in the leftmost panels show the value of $\sigma_\alpha(w=0)$ in the low $z$ (dashed) and high $z$ (continuous) bins}
    \label{fig:bao-avg-comb}
\end{figure}
\iffalse
\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{plots/combine_tpcf_avg.png}
    \caption{Combined Galaxy and Void 2PCF showing the difference between the constant radius cut used before \citep{Zhao2020} (shaded region) and the local radius cut proposed in this paper. The $1\sigma$ deviation is shown.}
    \label{fig:combine-tpcf-avg}
\end{figure*}
\fi

The top panel in \reffig{fig:snr-comb} shows the $\mathrm{SNR}$ vs. $w$ relations with both the local and constant cuts for catalogs before reconstruction. It is evident that in the low $z$ bin, the relation is very similar when comparing constant and local cuts, as would be expected given the similarities between the void correlations in this case. For the high $z$ bin however it is possible to see a significant $~95\%$ confidence interval (CI) increase in the $\mathrm{SNR}$, again compatible with the results from the void $\mathrm{SNR}$ only. The bottom panel in the same figure shows the corresponding post-reconstruction results. Notice that contrary to what is seen pre-reconstruction, the use of the local cut makes the $\mathrm{SNR}$ peak move towards smaller values of $\abs{w}$, almost two orders of magnitude smaller than the optimal $\abs{w}$ when using a constant cut. This is probably related to the fact that the number of voids obtained in the clustering sample, $\abs{\mathcal{V}_l}$, is $\sim1.15$ times that of the sample obtained from the constant cut. 

\reffig{fig:bao-avg-comb} shows the results of the BAO fit performed with the different $w$ parameters. We show the fitted $\alpha$, $\sigma_\alpha$, as well as the \textsc{Multinest} best-fit $\chi^2$ and Bayesian evidence $\mathcal{Z}$. 

Before reconstruction $\sigma_\alpha$ is very similar regardless of the radius cut chosen and suggests an optimal value of $w\approx-0.05$. Notice that to the left, the error on $\alpha$ increases significantly as the contribution of voids ($\sim\abs{w}$) does. The best-fit $\chi^2$ however shows a minimum around $w=-0.25$ for most cases, except for the constant cut in the high $z$ bin. 

Post-reconstruction BAO fits show a similar trend to that seen in the $\mathrm{SNR}$. For the constant cut the minimum $\sigma_\alpha$ remains around $w\approx-0.05$ while for the local cut the minimum lies at $w\approx-0.005$. The best-fit $\chi^2$ in this case is very low when using the constant cut, even when $w\rightarrow0$.

For pre reconstruction then, we use the values of $w=-0.05$ for the constant and local cuts. For post-reconstruction, we use $w=-0.05$ for the constant cut and $w=-0.005$ for the local.

Using these $w$ values, we perform individual fits on each mock to analyze the distribution of the $\sigma_\alpha$. The comparison of the constant and local cuts on \textsc{Patchy}{} BOSS DR12 lightcones and BOSS DR12 observations for both redshift bins before and after reconstruction are shown in \reffig{fig:bothcut-comparison}. The results corresponding to the constant cut here shown are consistent with those in \citet{Zhao2020}. Moreover, our fits using the local cuts show no significant difference in terms of BAO measurements. The BOSS observations (red star) show a slight improvement pre reconstruction but this is compatible with statistical fluctuations, based on the fits to the individual mocks. 
\begin{figure}
    \centering
    \includegraphics[width=0.49\linewidth]{plots/s_alpha_highz_bothcut_wt-0.05_pre.png}
    \includegraphics[width=0.49\linewidth]{plots/s_alpha_lowz_bothcut_wt-0.05_pre.png}
    \includegraphics[width=0.49\linewidth]{plots/s_alpha_highz_bothcut_wt-0.005_post.png}
    \includegraphics[width=0.49\linewidth]{plots/s_alpha_lowz_bothcut_wt-0.005_post.png}
    \caption{Distribution of $\sigma_\alpha$ fit from \textsc{Patchy}{} lightcones and BOSS DR12 data. Top plots show results before and bottom after reconstruction. Axis have been cut to show the largest probability density. Dotted lines show the mean of the distributions. Dashed lines and star show the BOSS DR12 results. Magenta line corresponds to the diagonal and cyan line to $\sigma_\alpha^{\mathrm{const.}} = m \sigma_\alpha^{\mathrm{local}}$.}
    \label{fig:bothcut-comparison}
\end{figure}


\section{Conclusions}
\label{sec:conclusions}
Throughout this work we have analyzed the response of DT void sample and its clustering to systematic incompleteness in the tracer distribution. In order to do this we have prepared the Non-Uniformly Angular downsampled and Non-Uniformly Radial downsampled data sets (\textsc{Patchy}{}-NUR and \textsc{Patchy}{}-NUA respectively). These correspond to \textsc{Patchy}{} boxes downsampled with a Gaussian selection function in the $Z$ direction and with a paraboloid in the $X,Y$ plane respectively. By doing this we aim at modeling radial selection effects and angular systematics as seen in galaxy surveys. In addition, we have studied the dependence of the DT void radii distribution with the tracer density before and after the reconstruction process. To this end we created the Uniformly downsampled (\textsc{Patchy}{}-U) data set.

We find fitting formulae to describe the BAO $\mathrm{SNR}$ maximizing radii selection criteria $R_c(n_{\mathrm{gal}};\gamma,\nu)$ using a power-law fit to the optimal radius vs. the galaxy density, where $\gamma$ indicates the amplitude and $\nu$ is the power index. We do this for auto and cross correlations with either small or large spheres. The BAO significance of large sphere -- under-dense regions that are typically referred to as cosmic voids -- autocorrelation is found to depend the most on the tracer density. The fit parameter $\nu$ is apparently shared in pre- and post- reconstruction fits, while $\gamma$ varies. This would allow to fix the $\nu$ parameter and grid-search for an appropriate $\mathrm{SNR}$-maximizing $\gamma$ for other samples. This sample will take into account local observed density information, we define it as $\mathcal{V}_l \equiv \qty{\mathtt{s}\in\mathcal{S} | R(\mathtt{s}) > \gamma n_{\mathrm{gal}}(\vb*{x}_{\mathtt{s}})^\nu}$.
Moreover, we observe that the pre-reconstruction large void-galaxy correlation has a larger $\abs{\mathrm{SNR}}$ than even the galaxy-galaxy correlation, which is the current standard for BAO measurements. 

In terms of the broad-band clustering, we show that the large-void sample defined as $\mathcal{V}_0 \equiv \qty{\mathtt{s}\in\mathcal{S} | R(\mathtt{s}) > R_c}$ (constant cut), contrary to what is observed from galaxies, is not sensitive to systematical effects that manifest as a small sample incompleteness such as the one generated by angular observational systematics of redshift surveys. On the other hand, the effects of a very large incompleteness, such as the one generated by the radial selection function, do induce large biases in the galaxy-galaxy, galaxy-void and void-void correlations, which means that these need a more careful estimation of the survey geometry, e.g. generation of random catalogs. Taking this into account, we find that the $\mathcal{V}_0$ sample shows a large void exclusion peak at low scales that is not consistent with the expectation from the correlations of complete \textsc{Patchy}{} data.  An alternative definition of the void sample $\mathcal{V}_l$ that takes into account the local observed density, greatly reduces the amplitude of this feature and yields clustering that is much more consistent with the complete \textsc{Patchy}{} sample. We explain this excess probability by noting that the constant cut defining the $\mathcal{V}_0$ sample allows for the selection of large spheres in undersampled regions of space that do not necessarily correspond to matter-field underdensities.

We compare the clustering of the $\mathcal{V}_0$ and $\mathcal{V}_l$ samples from \textsc{Patchy}{} lightcones pre- and post-reconstruction. We find that similar to what is seen in boxes, the pre-reconstruction void-void $\mathcal{V}_l$ correlations show a significantly lower exclusion in the high-$z$ bin, which contributes to a larger $\mathrm{SNR}$. In fact, the local cut enhances the $\mathrm{SNR}$ both pre- and post-reconstruction, nonetheless the reconstruction itself seems to already mitigate the exclusion effect in the 2PCF in such a way that the constant cut is already good enough.

Furthermore,  we evaluate the impact of the different void selection mechanisms on cosmological measurements. We observe that on simulations as well as on lightcones, the differences in terms of BAO measurement precision are negligible and that the $\alpha$ parameter is always recovered. Moreover, we observe that for the reconstructed samples, the exclusion feature is less intrusive upon BAO scales, even when using the simpler constant cut selection scheme. Therefore, we conclude that for reconstructed samples with BOSS (LRG-like) density, such as \textsc{Patchy}{} and BOSS DR12, the constant cut already delivers a robust clustering sample that allows for an unbiased measurement of the BAO peak position. In particular, this sample is shown to be robust even to unknown systematics that represent a moderate incompleteness (<20\%), for which the galaxy autocorrelation cannot be corrected, producing not only an unbiased void autocorrelation, but also an unbiased void-galaxy correlation with a generally more significant BAO signal than either autocorrelation.

Nonetheless, more sparse observations may be susceptible to large exclusion patterns that can intrude upon BAO scales. In these cases the local radius cut may play a more significant role in the void BAO measurement. Alternatively, the presence of this feature could be modeled and taken into account more precisely in the BAO fits. We leave this analysis for a future work. 

We believe that the approach described here, can in principle provide tighter BAO constraints than just the galaxy autocorrelation measurement, and we plan to apply this methodology to the upcoming DESI data.

\section*{Acknowledgements}
DFS, CZ, AV and AT acknowledge support from the SNF grant 200020\_175751. FSK acknowledges the grants SEV-2015-0548, RYC2015-18693, and AYA2017-89891-P.

\section*{Data Availability}

The light cone mocks and BOSS DR12 data are publicly available. The \textsc{Patchy}{} boxes used can be provided upon request to CZ.






\bibliographystyle{mnras}

\subsection{Bias dependence}
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/bias.png}
    \caption{Top panel: Void bias measured from the ``optimal'' void sample for different galaxy densities shown the colorbar. Bottom panel: Galaxy-galaxy power spectrum. Dashed lines correspond to results from the 20 and 50\% most massive galaxies where the voids are obtained through the LRG model for large void autocorrelations (\refeq{eq:pre-rdependence}). The dotted line shows the 20\% most massive sample but the voids are selected according to a re-fitted analogous model.}
    \label{fig:ll-bias}
\end{figure}

We computed the void bias $b_{gv}(k) = P_{gv}(k) / P_{g}(k)$ for samples of boxes with different galaxy densities, results are shown in \reffig{fig:ll-bias}. Here, $P_{gv}$ stands for the void-galaxy cross-power spectrum and $P_{g}$ the galaxy auto power spectrum.

Given that our datasets are tuned to closely emulate the BOSS LRG distribution, they have a set bias. We illustrate this in the bottom panel of \reffig{fig:ll-bias} where the continuous colored lines show the bias measurements for boxes where the lower density was obtained through random downsampling as explained before in this section. By randomly downsampling, we do not alter the bias of the resulting galaxy population and the large-sphere radius cut shows that the resulting void population has a constant negative linear bias with respect to the tracer density. 

We therefore perform a 20\% and 50\% largest-mass cuts in the complete tracer distribution to generate a sample with different bias as illustrated by the dashed lines in the bottom panel of \reffig{fig:ll-bias}. The top panel of the same figure shows that the resulting void population extracted using the radius cut in \refeq{eq:pre-rdependence} (labeled hereon as LRG model) produces different void bias sample in both top 20\% and 50\% massive galaxy samples. For the former in particular, we obtained a positively biased void sample which is not expected from underdensities \citep{Zhao2016}. In this case, we re-optimized the SNR as described before. We find that the resulting radius cut $R^* = (2.54\pm0.03)~\bar{n}_{\mathrm{gal}}^{-0.237\pm0.001}$ differs from the expected one (\refeq{eq:pre-rdependence}) in $(0.063\pm0.017)\%$. Even though the difference is a small percentage, it is significant over the 1$\sigma$ level and the bias measurement shows that it yields a negatively-biased void sample as shown by the dotted line in the same figure.

This analysis shows that maximizing the void autocorrelation BAO SNR is an effective way to obtain an optimized, negatively biased Cosmic Void sample. However our model so far seems to be able to fit only samples with a fixed bias. A more complex model is left for a future study.
