Integral constraints in spectroscopic surveys

TL;DR

This paper derives a comprehensive formalism for integral constraints in spectroscopic surveys, addressing their impact on clustering analyses, and proposes methods to mitigate systematic effects, especially on large scales.

Contribution

The paper introduces a complete derivation of integral constraints, including normalisation, shot noise, and wide-angle effects, and applies it to improve large-scale clustering analyses.

Findings

Radial integral constraints can significantly affect large-scale measurements.

Angular systematics can be mitigated by nulling density fluctuations on specific scales.

The formalism is validated with mock catalogues in RSD analysis.

Abstract

Clustering analyses of spectroscopic surveys are based upon density fluctuations, which are estimated by comparing the observed tracer density field to a selection function accounting for the survey density and geometry. However, this survey selection function is commonly partly inferred from the observed data itself, leading to so-called integral constraints, for which we propose a complete derivation. We discuss the normalisation of the introduced window functions, the shot noise contribution to the integral constraint corrections and wide-angle effects. Using this formalism, we review the well-known global integral constraint, arising when the expected mean galaxy density is taken to be the measured one. Another, stronger, constraint is imposed when the radial selection function itself is estimated from the data redshift distribution, as is often the case in the literature. We find that the impact of such a radial integral constraint can be as significant as the window function effect at large scales, depending on the survey geometry. Equations for this radial integral constraint are derived within our general formalism. We assess the validity of our approach by performing a Redshift Space Distortions (RSD) analysis on mock catalogues and emphasise that our results may be even more useful for analyses focusing on larger scales. Finally, as a further application, we show that unknown angular systematics can be mitigated by nulling the density fluctuations on a chosen angular scale. The induced loss of clustering is modelled by an angular integral constraint which can be combined with the radial one.

Tables

Table 1. Table 1 : Mean and standard deviation of the cosmological parameters fitted on the 84 N-series mocks, corresponding to figure  8 . Error bars should be divided by 84 ∼ 10 similar-to 84 10 \sqrt{84}\sim 10 to obtain errors on the mean of the mocks.

	Baseline global IC	Binned global IC	Binned radial IC	Expected
${\alpha }_{\parallel }$	$0.988\pm 0.031$	$0.979\pm 0.032$	$0.987\pm 0.031$	$0.989$
${\alpha }_{⟂}$	$0.980\pm 0.018$	$0.986\pm 0.019$	$0.980\pm 0.018$	$0.979$
$f{\sigma }_8$	$0.472\pm 0.037$	$0.459\pm 0.037$	$0.472\pm 0.037$	$0.470$
${\chi }^2$	$36.6\pm 9.7$	$42\pm 11$	$36.4\pm 9.6$	$\mathrm{ndof}=35$

Table 2. Table 2 : Mean and standard deviation of the cosmological parameters fitted on the 84 N-series mocks, corresponding to figure  11 . Error bars should be divided by 84 ∼ 10 similar-to 84 10 \sqrt{84}\sim 10 to obtain errors on the mean of the mocks.

	Baseline global IC	Pixelated angular IC	Systematics, pixelated angular IC	Expected
${\alpha }_{\parallel }$	$0.988\pm 0.031$	$0.991\pm 0.032$	$0.992\pm 0.032$	$0.989$
${\alpha }_{⟂}$	$0.980\pm 0.018$	$0.979\pm 0.020$	$0.979\pm 0.020$	$0.979$
$f{\sigma }_8$	$0.472\pm 0.037$	$0.471\pm 0.045$	$0.463\pm 0.044$	$0.470$
${\chi }^2$	$36.6\pm 9.7$	$36.9\pm 9.3$	$37.5\pm 9.7$	$\mathrm{ndof}=35$

Table 3. Table 3 : Mean and standard deviation of the cosmological parameters fitted on the 84 N-series mocks, corresponding to figure  12 . Error bars should be divided by 84 ∼ 10 similar-to 84 10 \sqrt{84}\sim 10 to obtain errors on the mean of the mocks.

	Baseline global IC	Binned, pixelated radial x angular IC	Systematics, binned, pixelated radial x angular IC	Expected
${\alpha }_{\parallel }$	$0.988\pm 0.031$	$0.988\pm 0.032$	$0.990\pm 0.032$	$0.989$
${\alpha }_{⟂}$	$0.980\pm 0.018$	$0.981\pm 0.020$	$0.980\pm 0.020$	$0.979$
$f{\sigma }_8$	$0.472\pm 0.037$	$0.475\pm 0.044$	$0.467\pm 0.043$	$0.470$
${\chi }^2$	$36.6\pm 9.7$	$36.5\pm 9.3$	$37.0\pm 9.6$	$\mathrm{ndof}=35$

Equations

$$F(\mathbf{r})=n_g(\mathbf{r})-\alpha n_s(\mathbf{r})$$

$$\alpha =\frac{\int d^{3}x{n}_g(\mathbf{x})}{\int d^{3}x{n}_s(\mathbf{x})}=\frac{{\sum }_{i=1}^{N_g}{w}_{g,i}}{{\sum }_{i=1}^{N_s}{w}_{s,i}}.$$

$$\begin{array}{rl}{\delta }^{\mathrm{cic}}(\mathbf{r})& =n_g(\mathbf{r})-\alpha n_s(\mathbf{r})\\ & =W(\mathbf{r})\left\{1+\delta (\mathbf{r})\right\}-W(\mathbf{r})\frac{\int d^{3}xW(\mathbf{x})\left\{1+\delta (\mathbf{x})\right\}}{\int d^{3}xW(\mathbf{x})}\\ & =W(\mathbf{r})\left\{\delta (\mathbf{r})-\int d^{3}x{W}_{\mathrm{glo}}(\mathbf{x})\delta (\mathbf{x})\right\},\end{array}$$

$$\alpha (r)=\frac{\int d^{3}x{n}_g(\mathbf{x}){ϵ}_{\mathrm{rad}}(r,x)}{\int d^{3}x{n}_s(\mathbf{x}){ϵ}_{\mathrm{rad}}(r,x)},$$

$$\frac{n_s(\mathbf{r})}{\int d^{3}x{n}_s(\mathbf{x}){ϵ}_{\mathrm{rad}}(r,x)}=\frac{W(\mathbf{r})}{\int d^{3}xW(\mathbf{x}){ϵ}_{\mathrm{rad}}(r,x)}.$$

$$\begin{array}{rl}{\delta }^{\mathrm{cic}}(\mathbf{r})& =n_g(\mathbf{r})-\alpha (r)n_s(\mathbf{r})\\ & =W(\mathbf{r})\left\{1+\delta (\mathbf{r})\right\}-W(\mathbf{r})\frac{\int d^{3}xW(\mathbf{x})\left\{1+\delta (\mathbf{x})\right\}{ϵ}_{\mathrm{rad}}(r,x)}{\int d^{3}xW(\mathbf{x}){ϵ}_{\mathrm{rad}}(r,x)}\\ & =W(\mathbf{r})\left\{\delta (\mathbf{r})-\int d^{3}x{W}_{\mathrm{rad}}(\mathbf{x})\delta (\mathbf{x}){ϵ}_{\mathrm{rad}}(r,x)\right\},\end{array}$$

$${\delta }^{\mathrm{cic}}(\mathbf{r})=W(\mathbf{r})\left\{\delta (\mathbf{r})-\int d^{3}x{W}_{\mathrm{ic}}(\mathbf{x})\delta (\mathbf{x}){ϵ}_{\mathrm{ic}}(\mathbf{r},\mathbf{x})\right\},$$

$$W_{\mathrm{ic}}(\mathbf{r})=\frac{W(\mathbf{r})}{\int d^{3}xW(\mathbf{x}){ϵ}_{\mathrm{ic}}(\mathbf{r},\mathbf{x})}.$$

$${\xi }^{\mathrm{cic}}(\mathbf{s})$$

$$-\int d^3\Delta \xi (\mathbf{\Delta })\int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})W_{\mathrm{ic}}(\mathbf{x}+\mathbf{\Delta }){ϵ}_{\mathrm{ic}}(\mathbf{x}+\mathbf{s},\mathbf{x}+\mathbf{\Delta })$$

$$-\int d^3\Delta \xi (\mathbf{\Delta })\int d^{3}xW(\mathbf{x})W(\mathbf{x}-\mathbf{s})W_{\mathrm{ic}}(\mathbf{x}+\mathbf{\Delta }){ϵ}_{\mathrm{ic}}(\mathbf{x}-\mathbf{s},\mathbf{x}+\mathbf{\Delta })$$

$$+\int d^3\Delta \xi (\mathbf{\Delta })\int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})$$

$$\int d^{3}y{W}_{\mathrm{ic}}(\mathbf{y})W_{\mathrm{ic}}(\mathbf{y}+\mathbf{\Delta }){ϵ}_{\mathrm{ic}}(\mathbf{x},\mathbf{y}){ϵ}_{\mathrm{ic}}(\mathbf{x}+\mathbf{s},\mathbf{y}+\mathbf{\Delta }).$$

$${\xi }_{\ell }^{\mathrm{c}}(s)=\sum _{p,q}{A}_{\ell p}^q\frac{2\ell +1}{2q+1}{\xi }_p(s){\mathcal{W}}_q^{\delta ,\delta }(s),$$

$${\mathcal{W}}_{\ell }^{\delta ,\delta }(s)=\frac{2\ell +1}{4\pi }\int d{\Omega }_s\int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s}){\mathcal{L}}_{\ell }({\hat{\eta }}_{\mathbf{s}}\cdot \hat{\mathbf{s}})$$

$${\mathcal{L}}_{\ell }(\mu ){\mathcal{L}}_p(\mu )=\sum _{q=0}^{\ell +p}{A}_{\ell p}^q{\mathcal{L}}_q(\mu ).$$

$$I{C}_{\ell }^{\delta ,\mathrm{ic}}(s)=\frac{2\ell +1}{4\pi }\int d{\Omega }_s\int d^3\Delta \xi (\mathbf{\Delta })\int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})W_{\mathrm{ic}}(\mathbf{x}+\mathbf{\Delta }){\mathcal{L}}_{\ell }({\hat{\eta }}_{\mathbf{s}}\cdot \hat{\mathbf{s}}){ϵ}_{\mathrm{ic}}(\mathbf{x}+\mathbf{s},\mathbf{x}+\mathbf{\Delta }).$$

$$\xi (\mathbf{\Delta })=\sum _{\ell }{\xi }_{\ell }(\Delta ){\mathcal{L}}_{\ell }({\hat{\eta }}_{\mathbf{\Delta }}\cdot \hat{\mathbf{\Delta }}),$$

$$\begin{array}{rl}I{C}_{\ell }^{\delta ,\mathrm{ic}}(s)& =\frac{2\ell +1}{4\pi }\int d{\Omega }_s\int {\Delta }^{2}d\Delta \int d{\Omega }_{\Delta }\sum _p{\xi }_p(\Delta )\\ & \int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})W_{\mathrm{ic}}(\mathbf{x}+\mathbf{\Delta }){\mathcal{L}}_{\ell }({\hat{\eta }}_{\mathbf{s}}\cdot \hat{\mathbf{s}}){\mathcal{L}}_p({\hat{\eta }}_{\mathbf{\Delta }}\cdot \hat{\mathbf{\Delta }}){ϵ}_{\mathrm{ic}}(\mathbf{x}+\mathbf{s},\mathbf{x}+\mathbf{\Delta }).\end{array}$$

$$\begin{array}{rl}{\mathcal{W}}_{\ell p}^{\delta ,\mathrm{ic}}(s,\Delta )& =\frac{(2\ell +1)(2p+1)}{(4\pi {)}^2}\int d{\Omega }_s\int d{\Omega }_{\Delta }\\ & \int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})W_{\mathrm{ic}}(\mathbf{x}+\mathbf{\Delta }){\mathcal{L}}_{\ell }({\hat{\eta }}_{\mathbf{s}}\cdot \hat{\mathbf{s}}){\mathcal{L}}_p({\hat{\eta }}_{\mathbf{\Delta }}\cdot \hat{\mathbf{\Delta }}){ϵ}_{\mathrm{ic}}(\mathbf{x}+\mathbf{s},\mathbf{x}+\mathbf{\Delta }),\end{array}$$

$$I{C}_{\ell }^{\delta ,\mathrm{ic}}(s)=\int {\Delta }^{2}d\Delta \sum _p\frac{4\pi }{2p+1}{\xi }_p(\Delta ){\mathcal{W}}_{\ell p}^{\delta ,\mathrm{ic}}(s,\Delta ).$$

$$I{C}_{\ell }^{\mathrm{ic},\mathrm{ic}}(s)=\int {\Delta }^{2}d\Delta \sum _p\frac{4\pi }{2p+1}{\xi }_p(\Delta ){\mathcal{W}}_{\ell p}^{\mathrm{ic},\mathrm{ic}}(s,\Delta ),$$

$$\begin{array}{rl}{\mathcal{W}}_{\ell p}^{\mathrm{ic},\mathrm{ic}}(s,\Delta )& =\frac{(2\ell +1)(2p+1)}{(4\pi {)}^2}\int d{\Omega }_s\int d{\Omega }_{\Delta }\int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})\\ & \int d^{3}y{W}_{\mathrm{ic}}(\mathbf{y})W_{\mathrm{ic}}(\mathbf{y}+\mathbf{\Delta }){\mathcal{L}}_{\ell }({\hat{\eta }}_{\mathbf{s}}\cdot \hat{\mathbf{s}}){\mathcal{L}}_p({\hat{\eta }}_{\mathbf{\Delta }}\cdot \hat{\mathbf{\Delta }}){ϵ}_{\mathrm{ic}}(\mathbf{x},\mathbf{y}){ϵ}_{\mathrm{ic}}(\mathbf{x}+\mathbf{s},\mathbf{y}+\mathbf{\Delta }).\end{array}$$

$${\xi }_{\ell }^{\mathrm{cic}}(s)={\xi }_{\ell }^{\mathrm{c}}(s)-I{C}_{\ell }^{\delta ,\mathrm{ic}}(s)-I{C}_{\ell }^{\mathrm{ic},\delta }(s)+I{C}_{\ell }^{\mathrm{ic},\mathrm{ic}}(s).$$

$$P_{\ell }^{\mathrm{cic}}(k)=4\pi (-i{)}^{\ell }\int s^{2}ds{j}_{\ell }(ks){\xi }_{\ell }^{\mathrm{cic}}(s),$$

$$I{C}_{\ell }^{\mathrm{glo},\mathrm{glo}}(s)=\frac{\int {\Delta }^{2}d\Delta {\xi }_0^{\mathrm{c}}(\Delta )}{\int {\Delta }^{2}d\Delta {\mathcal{W}}_0^{\delta ,\delta }(\Delta )}{\mathcal{W}}_{\ell }^{\delta ,\delta }(s).$$

$$P^{\mathrm{cic}}(\mathbf{k})=P^{\mathrm{c}}(\mathbf{k})-P^{\mathrm{c}}(0)|\tilde{W}(\mathbf{k}){|}^2$$

$$n_s(\mathbf{r})=W(\mathbf{r})\left\{1+\int d^{3}x{W}_{\mathrm{ic}}(\mathbf{x})\delta (\mathbf{x}){ϵ}_{\mathrm{ic}}(\mathbf{r},\mathbf{x})\right\},$$

$${\mathcal{S}}_{\ell }^{\delta ,\delta }(s)={\mathcal{W}}_{\ell }^{\delta ,\delta }(s)+I{C}_{\ell }^{\mathrm{ic},\mathrm{ic}}(s).$$

$$A={\alpha }^2\int d^{3}r{n}_s^2(\mathbf{r})={\alpha }^2\sum _{i=1}^{N_s}{w}_{s,i}{n}_{s,i},$$

$${\xi }_0^{\mathrm{c}}(s)\ni \int d^{3}xW(\mathbf{x})w(\mathbf{x}){\delta }^{3D}(s)$$

$$I{C}_{\ell }^{i,j}(s)\ni S{N}_{\ell }^{i,j}(s)(i,j)\in \{(\delta ,\mathrm{ic}),(\mathrm{ic},\delta ),(\mathrm{ic},\mathrm{ic})\},$$

Full text

1 Introduction

Spectroscopic surveys measure redshift-space positions of galaxies (or any other tracer of the matter density field). These surveys are not exhaustive and are characterised by a selection function $W(\mathbf{r})$, giving the expected density of observed galaxies at any redshift-space position $\mathbf{r}$ in the absence of clustering. In case one has a full knowledge of this selection function, density fluctuations can be fairly estimated by the difference of the observed density of galaxies to $W(\mathbf{r})$. In particular, density fluctuations averaged on the whole survey footprint can be non-zero, due to large-scale clustering modes.

However, the true survey selection function is difficult to determine in practice. Its norm, i.e. the expected mean density of galaxies, may be unknown. It is commonly taken to be the actually observed mean data density. Thus, the integral of the inferred density fluctuations over the whole survey footprint is set to zero, leading to a so-called integral constraint (IC) [1, 2], referred to as global IC in the following. Similarly, the radial survey selection function may be estimated from the data itself, as is currently done in the clustering analyses of the Baryon Oscillation Spectroscopic Survey (BOSS) (e.g. [3, 4, 5]) and its extended program eBOSS (e.g. [6, 7]). We propose to treat this effect as another integral constraint, dubbed radial IC.

Section 2 introduces the general formalism we develop to calculate these integral constraints. The computation of window functions required by our derivation is discussed in section 3.

We detail a Redshift Space Distortions (RSD) analysis of realistic BOSS-like mocks in section 4, showing that inferring the radial selection function from data can significantly alter clustering measurements, and that the induced bias is completely removed when using our modelling of the radial integral constraint.

As another application of our formalism, we suggest in section 5 to mitigate potential unknown angular systematics by imposing an angular integral constraint, which can be combined with the radial one.

We conclude in section 6.

2 Modelling integral constraints

We detail the impact of global and radial integral constraints on the observed density fluctuations in sections 2.1 and 2.2, showing that both can be accounted for by a general formalism, which we use to derive the effect of integral constraints on 2-point statistics in section 2.3. We discuss and compare our findings to previous results in section 2.4.

2.1 Impact of the global integral constraint on the observed density fluctuations

Predicting the mean galaxy density of a survey is complex in practice. For example, following [8] the Yamamoto power spectrum estimator [9] makes use of the FKP field:

[TABLE]

where $n_{g}(\mathbf{r})$ and $n_{s}(\mathbf{r})$ denote the density of observed and random galaxies, respectively. Random galaxies come from a Poisson-sampled synthetic catalogue accounting for the survey selection function111We assume throughout the paper that the density of this synthetic catalogue is sufficiently high to avoid systematic bias from undersampling the survey selection function.. Observed and random densities may include weights, e.g. corrections for systematics effects or a redshift weighting scheme (such as FKP weights, see appendix A). The observed galaxy density is $n_{g}(\mathbf{r})=W(\mathbf{r})\left\{1+\delta(\mathbf{r})\right\}$, with $W(\mathbf{r})$ the survey selection function and $\delta(\mathbf{r})$ the density contrast. We recall that $W(\mathbf{r})$ describes the expected density of observed, possibly weighted, galaxies in the absence of clustering. We assume the shape of the survey selection function is known, and sampled by the (weighted) synthetic catalogue: $n_{s}(\mathbf{r})\varpropto W(\mathbf{r})$. The scaling $\alpha$ is defined by:

[TABLE]

with $N_{g}$, $N_{s}$ and $w_{g}$, $w_{s}$ the number and weights of observed and random galaxies, respectively. Then, the observed window-convolved222In configuration space, we actually mean window-multiplied., integral-constraint-corrected (subscript $\mathrm{cic}$) density fluctuations are [1, 10]:

[TABLE]

with $W_{\mathrm{glo}}(\mathbf{r})=\frac{W(\mathbf{r})}{\int d^{3}xW(\mathbf{x})}$.

We find two terms, $\delta^{\mathrm{c}}(\mathbf{r})=W(\mathbf{r})\delta(\mathbf{r})$ corresponding to the density contrast weighted by the selection function, and the integral constraint term $\int d^{3}rW_{\mathrm{glo}}(\mathbf{r})\delta(\mathbf{r})$. The normalisation of $W_{\mathrm{glo}}(\mathbf{r})$ ensures that $\int d^{3}r\delta_{\mathrm{cic}}(\mathbf{r})=0$ over the entire footprint: modes larger than the survey size are suppressed.

2.2 Extension to the radial integral constraint

The true radial selection function of a spectroscopic survey is often a complex function of the luminosity function, sky lines, spectrograph efficiency and redshift determination algorithm [11]. As an example, in BOSS, the radial distribution of the synthetic catalogue is directly inferred from the data. Various techniques exist: random redshifts can be picked from the whole data redshift distribution (the so-called shuffled scheme [12, 13, 14]), or, assuming the true radial selection function should be somewhat smooth, it can be fitted by a spline from which random redshifts are drawn [12, 13]. A third possibility (binned scheme) is to weight an arbitrary initial random redshift distribution by $\alpha(r)$ (with $r=|\mathbf{r}|$) to match the data radial density in redshift or comoving distance bins of size $\delta r$, below which variations of $W(r)$ are neglected. Here, $\alpha(r)$ is given by:

[TABLE]

where $\epsilon_{\mathrm{rad}}(r,x)$ is $1$ if $r$ and $x$ belong to the same radial bin, [math] otherwise. By construction, the radially-normalised random density and true selection function match:

[TABLE]

Then, the following density fluctuations are observed:

[TABLE]

with $W_{\mathrm{rad}}(\mathbf{r})=\frac{W(\mathbf{r})}{\int d^{3}xW(\mathbf{x})\epsilon_{\mathrm{rad}}(r,x)}$. Note that the global integral constraint is also automatically imposed.

Equations (2.3) for the global IC and (2.6) for the radial IC can be recast in the following form:

[TABLE]

with $\epsilon_{\mathrm{ic}}(\mathbf{r},\mathbf{x})$ some generic window function ($\epsilon_{\mathrm{glo}}(\mathbf{r},\mathbf{x})=1$ for the global integral constraint, $\epsilon_{\mathrm{rad}}(r,x)$ for the radial one), and:

[TABLE]

The general expression for the observed density fluctuations (2.7) will be used in the following section to calculate the effect of integral constraints on 2-point statistics.

2.3 Correlation function in the local plane-parallel approximation

Pursuing our calculation in configuration space, we want to express the even multipoles $\xi^{\mathrm{cic}}_{\ell}(s)$ of the correlation function $\xi^{\mathrm{cic}}(\mathbf{s})$ of the observed density fluctuations (2.7). $\xi^{\mathrm{cic}}(\mathbf{s})$ would correspond to the quantity estimated with the Landy-Szalay estimator [15] if the division by the window function represented by the random pair counts $RR(\mathbf{s})$ was omitted in the estimator.

Here we assume the local plane-parallel approximation [10], namely that the pair separation $s$ between two galaxies at positions $\mathbf{x}$ and $\mathbf{x}+\mathbf{s}$ is small compared to the distances $|\mathbf{x}|$, $|\mathbf{x}+\mathbf{s}|$ of the galaxies to the observer. Then, the cosine angle betwen $\mathbf{s}$ and $\mathbf{x}$ or $\mathbf{x}+\mathbf{s}$ can be well approximated by $\hat{\mathbf{\eta}}_{\mathbf{s}}\cdot\hat{\mathbf{s}}$, with $\mathbf{\eta}_{\mathbf{s}}=\mathbf{x}+\mathbf{s}/2$ the mid-point line-of-sight and $\hat{\mathbf{s}}$ the unit $\mathbf{s}$ vector.

As $W(\mathbf{r})$ is uncorrelated with $\delta(\mathbf{r})$, the correlation function of the observed density fluctuations (2.7) can be expressed as:

[TABLE]

Term (2.9a) is the true correlation function $\xi(\mathbf{s})$ multiplied by the window function $\mathcal{W}^{\delta,\delta}(\mathbf{s})=\int d^{3}xW(\mathbf{x})W(\mathbf{x}+\mathbf{s})$. We account for the multipoles of (2.9a) using the formalism of [2]:

[TABLE]

where we use $\xi_{p}(s)$ the multipoles of $\xi(\mathbf{s})$. The window function multipoles are obtained by integration over the Legendre polynomials $\mathcal{L}_{\ell}$:

[TABLE]

and coefficient $A_{\ell p}^{q}$ is defined by:

[TABLE]

Cross-terms (2.9b) and (2.9c) account for the correlation between the density field and the integral constraint term of eq. (2.7). We are only interested in even multipoles with respect to the mid-point line-of-sight, so term $IC_{\ell}^{\mathrm{ic},\delta}$ (2.9c) is equal to $IC_{\ell}^{\delta,\mathrm{ic}}$ (2.9b). We are thus left with evaluating the multipoles of term (2.9b):

[TABLE]

The redshift-space correlation function is fully described by its multipoles. We thus develop:

[TABLE]

leading to:

[TABLE]

Writing the 3-point correlation function of the selection function:

[TABLE]

we obtain the simple formula:

[TABLE]

Finally, term (2.9d) accounts for the auto-correlation of the integral term in eq. (2.7). Similarly, its multipoles are equal to:

[TABLE]

with:

[TABLE]

Putting all terms together, the observed correlation function multipoles are:

[TABLE]

The observed window-convolved, integral-constraint-corrected power spectrum multipoles can be inferred from eq. (2.20) by a simple Hankel transform [2]:

[TABLE]

with $j_{\ell}$ the spherical Bessel function of order $\ell$. $P^{\mathrm{cic}}_{\ell}(k)$ corresponds to the quantity measured by the Yamamoto estimator [9] in the local plane-parallel approximation, to some global normalisation factor discussed in section 3.2. We report the reader to appendix B for a full derivation in the case of the end-point line-of-sight for the Yamamoto estimator, taking into account wide-angle corrections to the plane-parallel approximation. In the specific case of the global integral constraint ($\epsilon_{\mathrm{ic}}(\mathbf{r},\mathbf{x})=1$), the window function (2.19) is separable and:

[TABLE]

We will see how this compares to previous results in the next section.

2.4 Discussion

The global integral constraint was first discussed in [1, 10, 2], which expressed the observed window-convolved, integral-constraint-corrected power spectrum as:

[TABLE]

with $P^{\mathrm{c}}(\mathbf{k})$ the window-convolved power spectrum. $\tilde{W}(\mathbf{k})$ is the (Fourier transform of) the survey selection function, rescaled by imposing $|\tilde{W}(\mathbf{k})|^{2}=1$ when $k\ll 1\,h/\mathrm{Mpc}$, so that $P^{\mathrm{cic}}(k)\rightarrow 0$ as $k\rightarrow 0$. Only one correction term is assumed, $P^{\mathrm{c}}(\mathbf{0})|\tilde{W}(\mathbf{k})|^{2}$, while our complete derivation leads to 3 terms (2.20).

The multipoles of the integral constraint correction $P^{\mathrm{c}}(\mathbf{0})|\tilde{W}(\mathbf{k})|^{2}$ of eq. (2.23) are actually exactly equal to our last correction term, $IC_{\ell}^{\mathrm{glo},\mathrm{glo}}(s)$ after a Hankel transform. However, eq. (2.23) introduces a minus sign in front of $P^{\mathrm{c}}(\mathbf{0})|\tilde{W}(\mathbf{k})|^{2}$ to obtain a damping of the power spectrum at large scales. Instead, in our eq. (2.20), $IC_{\ell}^{\mathrm{glo},\mathrm{glo}}(s)$, which accounts for the auto-correlation of the density contrast integral over the footprint, leads to an increase of power at large scales. Damping of power comes from the two other terms $-IC_{\ell}^{\delta,\mathrm{glo}}(s)-IC_{\ell}^{\mathrm{glo},\delta}(s)$ which correspond to the cross-correlations between the density contrast and its integral over the footprint, respectively. Finally, note that our full derivation of the integral constraint effect makes it clear that the correction $P^{\mathrm{c}}(\mathbf{0})|\tilde{W}(\mathbf{k})|^{2}$, with $W$ describing the full survey footprint, is no longer valid if this footprint is composed of different chunks with separate normalisation factors (2.2). Indeed, in this case, $\epsilon_{\mathrm{glo}}(\mathbf{r},\mathbf{x})$ is $1$ if $\mathbf{r}$ and $\mathbf{x}$ belong to the same chunk and integrals of term (2.19) cannot be split as in eq. (2.22). The resulting global IC effect would be larger.

We detail the different terms contributing to the global integral constraint on the left panel of figure 1. As expected, on large scales, the window-convolved, integral-constraint-corrected monopole reaches [math]. The effect of the complete global integral constraint is negligible at the scales involved in a RSD analysis ($k\gtrsim 0.01\,h/\mathrm{Mpc}$), but is significant at large scales. As also shown in this panel, taking only $-IC_{\ell}^{\mathrm{glo},\mathrm{glo}}$, as in eq. (2.23) (with chunk-splitting — see section 4), for the global integral constraint correction (2.20) appears to be a very legitimate approximation in the illustrated survey case333We checked that this remains a good approximation for the radial integral constraint, as can be inferred from the right panel of figure 1 (dot-dashed and dashed curves). However, we observed that this approximation fails to account for the angular integral constraint described in section 5..

The right panel of figure 1 displays the different contributions to the radial integral constraint. Compared to the global IC, the radial IC has a larger effect on a broader range of wavenumbers, especially in the quadrupole and hexadecapole. We may thus expect a non-negligible impact of the radial integral constraint on clustering measurements in the illustrated case.

3 Computing window functions

This section is devoted to the practical computation of window functions $\mathcal{W}_{\ell}^{\delta,\delta}(s)$ and $\mathcal{W}_{\ell p}^{i,j}(s,\Delta)$ ($(i,j)\in\{(\delta,\mathrm{ic}),(\mathrm{ic},\delta),(\mathrm{ic},\mathrm{ic})\}$) required in section 2.3. We argue that a reasonable approximation for the true survey selection function can be inferred from data in section 3.1. We then discuss the normalisation of window functions, and detail their response to a Poisson shot noise in sections 3.2 and 3.3. Algorithms required to compute these window functions from a synthetic catalogue sampling the survey selection function are given in section 3.4.

3.1 Approaching the true survey selection function

Our derivation of integral constraint corrections makes use of the true, underlying, survey selection function $W(\mathbf{r})$. However, as we mentioned previously, one has often only access to a synthetic catalogue whose density $n_{s}(\mathbf{r})$ is tuned to match the actual data in certain regions (chunks, or radial bins for the radial IC). Thus, clustering modes inprinted in the data distribution are propagated to this synthetic catalogue. Using $n_{s}(\mathbf{r})$ for $W(\mathbf{r})$ does not change the value of $W_{\mathrm{ic}}(\mathbf{r})$ as it is normalised in regions where the integral constraint is imposed, i.e. where the expected mean of $n_{s}(\mathbf{r})$ is unknown. However, if that tuned $n_{s}(\mathbf{r})$ is taken for $W(\mathbf{r})$, the calculated window effect and integral constraint corrections statistically differ from the truth.

We can provide some estimate for this bias in the dominant term (2.10). The density $n_{s}(\mathbf{r})$ tuned on the actual data can be expressed as (see equations (2.3) and (2.6)):

[TABLE]

where $W_{\mathrm{ic}}(\mathbf{r})$ can be described by $n_{s,\mathrm{ic}}(\mathbf{r})=\frac{n_{s}(\mathbf{r})}{\int d^{3}xn_{s}(\mathbf{x})\epsilon_{\mathrm{ic}}(\mathbf{r},\mathbf{x})}$. Then, the multipoles of the 2-point window function estimated from $n_{s}(\mathbf{r})$ can be written:

[TABLE]

Taking $n_{s}(\mathbf{r})$ for $W(\mathbf{r})$ in $IC_{\ell}^{\mathrm{ic},\mathrm{ic}}(s)$, eq. (3.2) provides an estimate for $\mathcal{W}_{\ell}^{\delta,\delta}(s)$ at first order in the integral constraint corrections.

Note that such a dependence of the estimation of window effect and integral constraint corrections in the observed clustering also leads to a change in the variance of cosmological measurements, whose study is beyond the scope of this paper.

3.2 Normalisation

In order to recover the true (unconvolved) power spectrum at small scales, [2] suggested to normalise the 2-point window function in the limit $s\rightarrow 0$. However, this limit may be ill-defined in practice, since its estimation requires a very fine sampling of the survey selection function.

Similarly, in current Fourier space analyses, power spectrum measurements are normalised by [4]:

[TABLE]

using the notations of eq. (2.2). This also corresponds to the ill-defined limit $\mathcal{W}_{0}^{\delta,\delta}(0)=\int d^{3}rW^{2}(\mathbf{r})$ if $W(\mathbf{r})$ is represented by the same synthetic catalogue as in eq. (3.3).

Besides, the estimation of the density $n_{s}(\mathbf{r})$ in eq. (3.3) is in practice non-trivial when accounting for various survey selection effects. For example, $n_{s}(\mathbf{r})$ is commonly taken to be the redshift density $n(z)$, computed by binning (weighted) data in redshift slices [14], while it can also be a function of the angular position on the sky.

We thus also use term (3.3) in the normalisation of window functions $\mathcal{W}_{\ell}^{\delta,\delta}(s)$ and $\mathcal{W}_{\ell p}^{i,j}(s,\Delta)$, so that $A$ terms divide both the power spectrum measurements and model, and compensate. Therefore, the estimation of $A$ does not impact the estimation of cosmological parameters. Integral constraint terms are properly normalised using eq. (2.8).

3.3 Shot noise

In principle, formulae given in section 2.3 (and eq. (3.2)) fully describe integral constraint corrections. However, they are derived in real space, where accounting for a Poisson shot noise term (a spike at $s=0$) is numerically challenging.

The shot noise contribution to the correlation function weighted by the window function only shows up in the monopole [8, 9]:

[TABLE]

where $w(\mathbf{x})$ is the total (data and random444The total shot noise estimate (equations (3.4) and (3.5)) is valid in the limit of infinite random density, but appears to be a legitimate approximation if the random density is high enough, as shown by our cosmological fits in sections 4.3, 5.2, and 5.3.) weight at position $\mathbf{x}$. In Fourier space, this translate in a simple offset in the monopole. However, the shot noise contribution to the integral constraint correction leaks on all scales of all multipoles:

[TABLE]

with:

[TABLE]

and

[TABLE]

where we used $w_{\mathrm{ic}}(\mathbf{r})=\frac{w(\mathbf{r})}{\int d^{3}xW(\mathbf{x})\epsilon_{\mathrm{ic}}(\mathbf{r},\mathbf{x})}$.

Integrals over $\mathbf{s}$ of shot noise contributions to $\xi^{\mathrm{c}}_{0}(s)$ and $IC_{\ell}^{i,j}(s)$ are equal, such that total shot noise contributions vanish as $k\rightarrow 0$, as expected. In practice, we choose to normalise terms (3.4) and (3.5) by their integral over $\mathbf{s}$555In the case $W$ is sampled by a synthetic catalogue, the normalisation for $SN_{\ell}^{\mathrm{ic},\delta}(s)$ and $SN_{\ell}^{\mathrm{ic},\mathrm{ic}}(s)$ is the (weighted) number of correlated pairs., and multiply them by the shot noise measured by the Yamamoto estimator [9, 10].

Normalised shot noise contributions to the global and radial integral constraints are shown in figure 2. By definition, they reach $1$ as $k\rightarrow 0$ in the monopole. Their impact on the quadrupole and hexadecapole should not be ignored in the case of the radial IC.

3.4 Calculation from a synthetic catalogue

The survey selection function entering window functions $\mathcal{W}_{\ell}^{\delta,\delta}(s)$ and $\mathcal{W}_{\ell p}^{i,j}(s,\Delta)$ ($(i,j)\in\{(\delta,\mathrm{ic}),(\mathrm{ic},\delta),(\mathrm{ic},\mathrm{ic})\}$) required in section 2.3 can be randomly sampled by a synthetic catalogue. We use a classic pair-count algorithm to compute the anisotropic 2-point correlation and implement the algorithm from [16] to compute the anisotropic 3-point correlation of the synthetic catalogue. As detailed in appendix B, lines-of-sight are defined according to the FKP power spectrum estimator, which we use for clustering measurements.

Pairs (or triplets) of random objects are binned in comoving distance bins of typical size $1\,\mathrm{Mpc}/h$ (for $\mathcal{W}_{\ell}^{\delta,\delta}(s)$) to $4\,\mathrm{Mpc}/h$ (for $\mathcal{W}_{\ell p}^{i,j}(s,\Delta)$); the result is divided by the comoving volume of each distance bin to recover density, renormalised to data counts (by multiplying by $\alpha^{2}$), and divided by the normalisation factor $A$ as stated in section 3.2.

Window function calculations $\mathcal{W}_{\ell p}^{i,j}(s,\Delta)$ involving $\epsilon_{\mathrm{ic}}(\mathbf{r},\mathbf{x})$ are performed by correlating random objects at position $\mathbf{r}$ with random objects at position $\mathbf{x}$ satisfying $\epsilon_{\mathrm{ic}}(\mathbf{r},\mathbf{x})\neq 0$, i.e. residing in the same spatial bin.

We show examples of $\mathcal{W}_{\ell}^{\delta,\delta}(s)$ and $\mathcal{W}_{\ell p}^{\mathrm{rad},\mathrm{rad}}(s,\Delta)$ in figure 3. All window functions (for the global and radial integral constraints) could be accurately computed in $\sim 1500$ CPU hours.

4 Radial integral constraint in RSD analyses

In this section we discuss the impact of the radial integral constraint on Fourier-space RSD measurements, and show that eq. (2.20), along with practical details developed in section 3, helps us recover an unbiased estimation of cosmological parameters. We consider three sets of mock catalogues: the baseline mocks for measurements including only the global IC, and the binned and shuffled mocks for the radial one. These mocks are described in section 4.1, while details about power spectrum measurements and the RSD modelling are given in section 4.2. Cosmological fits are performed in section 4.3. We discuss the potential impact of the radial integral constraint on other clustering analyses in section 4.4.

4.1 Mock catalogues

We work on a set of 84 high fidelity N-series mocks used for the LRG clustering mock challenge [17] of the BOSS Data Release 12 (DR12) [5]. These mocks were built from the 3 projections of 7 independent, periodic box realisations of side $2600\,\mathrm{Mpc}/h$ at redshift $z=0.5$. The simulated cosmology is a flat $\Lambda$CDM model with $\Omega_{m}=0.286$, $\Omega_{\Lambda}=0.714$, $\sigma_{8}=0.82$, $n_{s}=0.96$ and $h=0.7$. Simulations were run with the GADGET2 code [18], and a HOD modelling was used to populate dark matter halos with galaxies, mimicking the observed clustering in data.

The precision matrix required for cosmological fits is built from 2048 Multidark Patchy mocks provided by the BOSS collaboration. These approximate mocks were calibrated on BigMultiDark simulations [19]; halo abundance matching was applied to reproduce the 2- and 3-point clustering measurements and mocks at different redshifts were combined into light cones [20]. We use V6C catalogues, which were adjusted to reproduce the data clustering measurements. These mocks are also used to compute the standard deviation (blue shaded area) shown in figures 5, 9 and 13 and the error bars plotted in figure 7.

Both N-series and Multidark Patchy mocks were trimmed following the DR12 CMASS NGC (hereafter abbreviated by CMASS) selection function using the make_survey software [21]. We impose integral constraints consistently on N-series and Multidark Patchy mocks in each series of cosmological fits. Baseline mocks are obtained by using a random catalogue sampling the true survey selection function, normalised according to eq. (2.2), so that those mocks are impacted by the global integral constraint only. The binned or shuffled schemes are used to impose the radial integral constraint. In the shuffled scheme, random redshifts are drawn from each mock data realisation. In the binned scheme, random galaxies are weighted to reproduce the mock data radial distribution in comoving distance bins of size $\delta r=2\,\mathrm{Mpc}/h$.

To enhance the radial integral constraint effect coming from the binned or shuffled schemes, we divide the CMASS footprint ($\simeq 7420\deg^{2}$) in $6$ smaller chunks of size ranging from $\simeq 980\deg^{2}$ to $\simeq 1570\deg^{2}$ (see figure 4), areas representative of the 4 chunks of the eBOSS ELG survey [22]. The binned or shuffled schemes are applied to the $6$ chunks separately, before these are recombined in a single catalogue used for the power spectrum measurement. This procedure, though not representative of the real BOSS DR12 analysis, allows us to test our modelling of the radial integral constraint in stringent conditions.

One can appreciate the difference between the true survey selection function (continuous black curve) and the ones measured from one mock realisation (dashed black curve for full CMASS and colored dashed curves for each chunk). Note that the radial selection function of future spectroscopic surveys like DESI [23] and Euclid [24] may vary on the sky; modelling this effect may indeed require to combine small patches with locally constant radial selection functions.

We apply FKP weights [8] $w_{\mathrm{FKP}}=\frac{1}{1+n(z)P_{0}}$ to both mock data and randoms using $P_{0}=20000\,(\mathrm{Mpc}/h)^{3}$. The redshift density $n(z)$, which also enters the normalisation of power spectrum measurements (see eq. (3.3)), is calculated in bins of $\Delta z=0.005$.

The redshift density is computed according to the type of mocks.

When the radial IC (binned or shuffled mocks) is imposed, the redshift density $n(z)$ is computed for each CMASS mock data realisation, as is done in an actual data analysis (see e.g. [14]). As we divide the CMASS footprint into $6$ chunks, one may want to measure $n(z)$ in the $6$ chunks separately (see figure 4, right-panel). However, we show in appendix A that doing so significantly biases clustering measurements as FKP weights, using local $n(z)$, smooth out clustering. We therefore choose to measure $n(z)$ from the whole CMASS sample, making the bias due to FKP weights negligible.

Similarly, when the global IC (baseline mocks) is imposed, we compute $n(z)$ on the random catalogue accounting for the true survey selection function (see the continuous black curve on the left panel of figure 4) but rescale it according to the weighted number of data for each mock realisation.

4.2 Analysis methods

We use the implementation of the Yamamoto estimator in the Python toolkit nbodykit [25] to compute mock data power spectra. We take a large box size of $4000\,\mathrm{Mpc}/h$ to reduce sampling effects in low $k$-bins. The FKP field is interpolated on a $512^{3}$ mesh following the Triangular Shaped Cloud (TSC) scheme. The Nyquist frequency is thus $k\simeq 0.4\,h/\mathrm{Mpc}$, more than twice larger than the maximum wavenumber used in our analysis. We employ the interlacing technique [26] to mitigate aliasing effects. Power spectrum multipoles are measured in bins of $\Delta k=0.01\,h/\mathrm{Mpc}$, from $k=0\,h/\mathrm{Mpc}$.

Power spectrum multipoles are calculated on a discrete $k$-space grid, making the angular modes distribution irregular at large scales. We correct for this effect using the technique employed in [4]. The correction is very small, given the large box size used to measure power spectra. Hankel transforms between power spectrum and correlation functions multipoles (e.g. eq. (2.21)) are performed using the FFTLog [27] software. As in [4] we only consider correlation function multipoles up to $\ell=4$ in our calculations. We checked that adding $\xi_{6}(s)$ has a completely neligible impact on the model prediction.

Mock data are fitted by the RSD TNS [28] model as implemented in [4]. However, contrary to [4], we adopt a Lorentzian form for the Finger-of-God effect [29]:

[TABLE]

with $\sigma_{v}$ the velocity dispersion. A common practice of Fourier space clustering analyses is to remove the shot noise contribution from the power spectrum monopole measurement. Instead, we fit the power spectrum monopole with its shot noise and add up in the model the shot noise term $N$ measured from the data, including contributions from integral constraints as shown in section 3.3. As in [4], we let the constant galaxy stochastic term (see [30]) free.

Power spectrum monopole, quadrupole and hexadecapole are fitted from $0.01\,h/\mathrm{Mpc}$ to $0.15\,h/\mathrm{Mpc}$. Fitted cosmological parameters are the logarithmic growth rate of structure $f$ and the Alcock-Paczynski [31] parameters $\alpha_{\parallel}$ and $\alpha_{\perp}$. As $f$ is almost completely degenerate with the power spectrum normalisation $\sigma_{8}$, we quote the combination $f\sigma_{8}$. We consider $4$ nuisance parameters: the linear and second order biases $b_{1}$ and $b_{2}$, the velocity dispersion $\sigma_{v}$ and $A_{g}=N_{g}/N$, with $N_{g}$ the constant galaxy stochastic term and $N$ the measured Poisson shot noise. Minimisations are performed using the algorithm Minuit [32, 33], taking large variation intervals for all parameters. We checked that the fitted parameters do not reach the input boundaries.

We use the fiducial BOSS DR12 cosmology in our analysis:

[TABLE]

In all figures showing the power spectrum model alone (1, 6), we use $f=0.75,\,b_{1}=2,\,b_{2}=1,\,\sigma_{v}=4\,\mathrm{Mpc}/h$.

The radial integral constraint correction is consistently calculated using in $\epsilon_{\mathrm{rad}}(\mathbf{r},\mathbf{x})$ the same comoving distance bins of size $\delta r=2\,\mathrm{Mpc}/h$ as used in the binned scheme. We account for the chunk-splitting of the CMASS footprint by adding the condition that $\mathbf{r}$ and $\mathbf{x}$ should belong to the same chunk for $\epsilon_{\mathrm{rad}}(\mathbf{r},\mathbf{x})$ to be non-zero.

Window function calculations are based on different estimations of the survey selection function, according to the type of mocks.

In the radial IC case (binned or shuffled mocks), to mimic an actual data analysis, the radial part of $W(\mathbf{r})$ should in principle be estimated on each realisation of the mocks, for which specific window function and integral constraint corrections should be derived. However, this would require a large computation time. Thus, in this work, window function estimations for the radial IC case are based on the selection function $n_{s}(\mathbf{r})$ calculated with a synthetic catalogue tuned to match the radial distribution (using the binned scheme) in each chunk of one realisation of mock data. All illustrations of window functions (figure 3) and integral constraint corrections (figures 1, 2, 6, 10) are provided for this estimation of the survey selection function.

Similarly, in the global IC case (baseline mocks), $W(\mathbf{r})$ entering window function calculations is taken to be the true selection function, which should be normalised using eq. (2.2) for each realisation of the mocks, as done in a real data analysis. However, as for the radial IC case, we normalise the selection function on one realisation of the mocks only (the same realisation as in the radial IC case).

We checked that using in the model the estimate provided by eq. (3.2) (with $\mathrm{ic}=\mathrm{glo}$ in the baseline case, $\mathrm{ic}=\mathrm{rad}$ in the binned or shuffled cases) for $\mathcal{W}_{\ell}^{\delta,\delta}(s)$ has a negligible impact on the measurement of cosmological parameters in our analysis. Hence, we do not include this correction in our cosmological fits.

4.3 Cosmological fits with the radial integral constraint

Figure 5 illustrates the effect of tuning the radial selection function on data. Both the binned scheme (using comoving distance bins of size $\delta r=2\,\mathrm{Mpc}/h$) and the shuffled scheme result in a loss of power parallel to the line-of-sight, thus mostly affecting the power spectrum quadrupole and hexadecapole, compared to the baseline relying on the true survey selection function. The effect of the radial integral constraint is significantly increased when it is imposed to the $6$ chunks separately (orange and red curves), which is the case considered in the following cosmological fits. We thus expect that neglecting the radial integral constraint would lead to a significant bias in the fitted cosmological parameters. Figure 5 also shows that power spectrum measurements obtained with the binned and shuffled scheme are very similar. It is thus fair to account for e.g. the shuffled scheme with the binned scheme, for which we derived integral constraint corrections in section 2.2. We also checked that differences between these two schemes in the measurement of cosmological parameters are negligible (i.e. small compared to the uncertainty on the mean of the mocks).

Figure 6 displays the model for the global and radial integral constraint corrections corresponding to the CMASS footprint divided in 6 chunks. By construction, the power spectrum monopoles converge to zero at large scales when the global or radial IC is applied. We stress again that window functions are normalised according to section 3.2, without using any low-$k$ nor low-$s$ limit. In the specific, illustrated case, the effect of the radial integral constraint is found to be large (comparable to the window function effect alone) in both Fourier and configuration space correlation functions. Note however that the latter prediction cannot be directly compared to measurements using e.g. the Landy-Szalay estimator [15] from which the window function effect is already removed.

As a preview, figure 7 illustrates the agreement between the mean of the 84 N-series mocks and the appropriate model in the baseline (left panel) and the binned scheme (right panel). In the baseline case, we fit the TNS model including the global integral constraint. The model performs well enough, given the high measurement precision (reduced $\chi^{2}=77.4/(42-7)=2.21$). Note that same cosmological and nuisance parameters are kept for the model in the right-hand plot, where the binned scheme is applied to the mocks, and the radial integral constraint is included in the model. Here again, the agreement betwen the model and the mean of the mocks is correct (reduced $\chi^{2}=82.2/(42-7)=2.35$). We thus anticipate that our radial integral constraint correction will very well account for the binned and shuffled schemes in the cosmological fits.

Figures 8 and table 1 present the cosmological fit results obtained in three different cases.

The baseline cosmological fits (column 1 in table 1) are obtained with the power spectrum measurements on the baseline mocks and with the global integral constraint applied in the model. The expected parameter values (column 4 in table 1), predicted from the mock fiducial cosmology, are recovered to the statistical uncertainty on the mean of the mocks.

Applying the binned scheme to the mocks (see the orange curve in figure 5) while modelling the global integral constraint only results in a bias (compared to the baseline case) on all cosmological parameters of roughly $30\%$ of the statistical error on one realisation (column 2 in table 1). The goodness-of-fit (probed by the $\chi^{2}$ distribution, figure 8) is significantly degraded ($\Delta\chi^{2}\simeq 5$).

The modelling of the radial integral constraint successfully removes the bias to better than the statistical uncertainty on the mean of the mocks (column 3 in table 1), and the goodness-of-fit is well recovered. No increase on cosmological parameter errors is detected. Similar results (not reported here) were obtained with the shuffled scheme.

4.4 Discussion

We showed that the tuning of the radial selection function on data may significantly bias RSD measurements. We recall that the radial IC effect has been purposedly enhanced to test the robustness of our predictions by applying the binned scheme separately in $6$ small chunks cut into the CMASS footprint. After modelling of the induced radial integral constraint, the bias on $f\sigma_{8}$ and on Alcock-Paczynski parameters is well below the uncertainty on the mean of the mocks. More statistics would be required to test our modelling further, but we note that the obtained bias is already below the expected $\simeq 1\%$ statistical uncertainty on $f\sigma_{8}$ provided by future galaxy surveys [23, 24].

The red curve in figure 5 shows that the shuffled (or binned) scheme applied to the CMASS NGC footprint as a whole (i.e. without the division in $6$ chunks) induces a small radial integral constraint on clustering measurements, as already noted on SDSS DR7 and DR9 by [12, 13]. We checked that not accounting for the radial IC in the model results in a negligible shift of cosmological parameters (compared to the statistical uncertainty) $\Delta\alpha_{\parallel}=-0.003$, $\Delta\alpha_{\perp}=0.001$ and $\Delta f\sigma_{8}=-0.001$. Note that, contrary to our test case, the BOSS DR12 analysis used the NGC and SGC CMASS and LOWZ samples. These were divided in three overlapping redshift slices [5], thus probably decreasing the radial integral constraint correction with respect to the pure window function effect. We thus infer that conclusions of the BOSS DR12 analysis are likely unaffected by the radial integral constraint effect.

Nonetheless, RSD measurements based on deep, pencil-like or small surveys such as the eBOSS ELG program [22] may suffer from the radial integral constraint, if not correctly accounted for. Also, analyses focusing on large scales, aiming at e.g. setting constraints on primordial non-Gaussianity [34, 35] may benefit from accounting for the radial IC.

5 Angular integral constraint to mitigate angular systematics

In this section we show that potential unknown angular systematics (section 5.1) can be mitigated using an angular integral constraint (section 5.2) which can be combined with the radial one (section 5.3). Though we note in section 5.4 that only the additive part of systematics is removed in this way, we emphasise that the technique described hereafter may be simply used as a consistency check in clustering analyses.

5.1 Problem statement

As the radial selection function, the angular selection function of a spectroscopic survey can be difficult to evaluate because of residual photometric calibration errors or other potential photometric systematics. We consider these systematics to be completely unknown; that is, they are not and cannot be corrected by any photometric template in the analysis.

As a case study, we inject photometric systematics into our mocks as a function of right ascension (R.A.) and declination (Dec.), using a weight:

[TABLE]

The amplitude of the systematics ($\pm 20\%$) is very large compared to the typical requirements of a target selection (e.g. $\pm 7.5\%$ in [36]). $w_{\mathrm{sys}}$ has a drastic impact on power spectrum measurements as can be seen on the left-hand plots of figure 9 (orange and blue curves). A standard RSD analysis, using scales $0.01\,h/\mathrm{Mpc}<k<0.15\,h/\mathrm{Mpc}$ would be impossible. We will see that these systematics can be strongly reduced using a similar procedure as in section 4.

5.2 Angular integral constraint

We suggest to remove the contaminated modes by weighting randoms from the synthetic catalogue by $\sum_{i\in\mathrm{pixel}}w_{g,i}/\sum_{i\in\mathrm{pixel}}w_{s,i}$ in pixels (the pixelated scheme), thus nulling the density fluctuations in each pixel. We use a HEALPix666http://healpix.jpl.nasa.gov/ [37] map with $\mathrm{nside}=64$ (pixel area of $\simeq 0.84\deg^{2}$), for which the contaminated and uncontaminated mocks look similar (see red and green curves in figure 9, left). The pixelated scheme completely mitigates angular systematics in the quadrupole and hexadecapole. However, a difference close to a scale factor and statistically significant at small scales remains in the monopole.

We model the pixelated scheme by an angular integral constraint. A similar idea is developed in [38] to mitigate the impact of the DESI [23] fiber assignment.

Formulae are directly deduced from section 2 by replacing $\epsilon_{\mathrm{ic}}(\mathbf{x},\mathbf{y})$ by $\epsilon_{\mathrm{ang}}(\mathbf{x},\mathbf{y})$, being non-zero if $\mathbf{x}$ and $\mathbf{y}$ lie within the same pixel. Contrary to the radial IC case (binned or shuffled schemes), for the angular IC (pixelated scheme) one has access to an estimate of the survey selection function accounting for known systematics (weighted synthetic catalogue), independent of the observed data. Thus, in the following, window function calculations are based on the true survey selection function, as in the baseline case. The angular IC correction is calculated from the same window functions for both uncontaminated and contaminated pixelated mocks, though in the latter case the normalisation of the survey selection function would be inferred from contaminated data in a practical analysis. However, we do not expect any systematic bias from this simplification, as introduced systematics are independent of the actual clustering.

The shot noise contribution to the angular integral constraint is very large in the quadrupole and hexadecapole, as can be seen in figure 10.

In all this section, we use the same model and fitting $k$-range as in section 4. Applying the pixelated scheme to uncontaminated mocks, cosmological parameters are well recovered when modelling the angular integral constraint: as shown in figure 11 and table 2, differences with the baseline case (same as in section 4) are within the uncertainty on the mean of the mocks (column 2 in table 2). The error bar on $f\sigma_{8}$ increases by $22\%$.

The magnitude of angular systematics makes it impossible to perform any relevant standard cosmological fit: we measure biases of $\Delta\alpha_{\parallel}=0.043$, $\Delta\alpha_{\perp}=-0.037$ and $\Delta f\sigma_{8}=-0.17$. On the contrary, the cosmological analysis remains possible when applying the pixelated scheme to the contaminated mocks and the angular integral constraint in the model (column 3 in table 2). Alcock-Paczynski parameters are recovered within the statistical uncertainty on the mean of the mocks, and a bias of $\simeq 20\%$ of the error on a single realisation can be seen on $f\sigma_{8}$. A lower bias on the growth rate would be expected with smaller, more realistic angular systematics.

5.3 Combining radial and angular integral constraints

One would probably like to combine the radial and angular integral constraints, to account for both unknown radial and angular selection functions. A prerequisite to applying the radial integral constraint is that the redshift distribution does not depend on the angular position on the sky in each chunk of the survey. Thus, both radial and angular integral constraints can be imposed at the same time since radial and angular selection functions are independent.

Right-hand plots of figure 9 show the effects of the combined binned and pixelated schemes on N-series mocks. They add up in the hexadecapole, and partially cancel in the quadrupole. This is expected as $\mathcal{L}_{2}$ is negative around $\mu=0$ (where the angular integral constraint removes signal) and positive around $\mu=1$ (where the radial integral constraint plays up). As with the pixelated scheme alone (left-hand plots of figure 9), angular systematics are mitigated in the quadrupole and hexadecapole, but a multiplicative effect remains in the monopole.

Let us model the radial and angular integral constraints in a row:

[TABLE]

Since the radial and angular parts of $W(\mathbf{r})$ are independent, the last term is just the integral of the density contrast over the whole (chunk) footprint, i.e. the global integral constraint. Thus, the two integral constraints commute. Then, building up the correlation function multipoles of the observed density fluctuations (5.2), we find 16 terms:

[TABLE]

Terms (5.3a), (5.3b) and (5.3c) correspond to the radial, angular and global integral constraints, while terms (5.3) are the cross-integral constraints, given by a formula similar to equations (2.17) and (2.18), with:

[TABLE]

As in section 4.2, window function calculations are based on the survey selection function radially-tuned on one realisation of mock data using the binned scheme. As in section 5.2, the same window functions are used to calculate the angular IC correction for both uncontaminated and contaminated pixelated mocks.

As shown in figure 12 and table 3, the combined (radial x angular) integral constraint accounts well for the binned and pixelated schemes in the mocks: cosmological parameters are recovered within the uncertainty on the mean of the mocks (column 2 in table 3). Cosmological fits of contaminated mocks are not further degraded by adding the radial integral constraint on top of the angular one. As in section 5.2, a bias of $12\%$ to $20\%$ of the error on a single realisation can be seen on $f\sigma_{8}$ (comparing column 3 to columns 2 and 1).

5.4 Caveat: multiplicative systematics

We emphasise that the angular (pixel) integral constraint can only mitigate the additive part of angular systematics. Let us call $c(\mathbf{r})$ the contamination signal, assumed constant over a pixel. First, let us suppose $c$ to be purely additive. Applying the pixelated scheme, we would measure the power spectrum of the density fluctuations:

[TABLE]

Then, as $c(\mathbf{r})$ is constant over a pixel, $c(\mathbf{x})\epsilon_{\mathrm{ang}}(\mathbf{r},\mathbf{x})=c(\mathbf{r})\epsilon_{\mathrm{ang}}(\mathbf{r},\mathbf{x})$ and:

[TABLE]

As expected, $c(\mathbf{r})$ disappears from the analysis.

Now, let us consider $c(\mathbf{r})$ to be multiplicative (as implemented in our contamination model (5.1)):

[TABLE]

i.e. $c(\mathbf{r})$ multiplies the selection function $W(\mathbf{r})$. However, by definition, $c(\mathbf{r})$ is unknow and therefore cannot be taken into account in $W(\mathbf{r})$. The resulting multiplicative systematics can explain the observed bias on $f\sigma_{8}$ (compensating for a higher $b_{1}\sigma_{8}$). A possible way to alleviate this effect may consist in estimating the angular survey selection function (including $c(\hat{\mathbf{r}})$) with the pixelated scheme. This would induce a bias which can be estimated using the method presented in section 3.1. We report the reader to e.g. [39] for a fully coherent treatment of unknown multiplicative systematics.

We showed that one could mitigate potential, unknown, systematics at the price of a small increase of statistical errors by applying an angular integral constraint in the data. Such a technique may be helpful to reduce the impact of poorly understood angular target density variations with photometric conditions, or can be considered as a regular consistency test to check for remaining angular systematics.

6 Conclusions

This paper revisited the notion of integral constraints which we showed to be useful in clustering analyses to account for biases related to calibrating the survey selection function partly on data.

We presented a general formalism to account for integral constraints that completes the existing models for the global integral constraint [1, 10, 2]. We discussed the window function normalisation and the shot noise contribution to the integral constraints and included wide-angle effects following [40]. Our formalism was used to model the radial integral constraint effect induced by inferring the radial selection function from actual data.

We indeed noticed that the common practice of drawing random redshifts based on data may significantly bias large-scale clustering measurements. This effect can be particularly large compared to the window function effect if the survey is composed of several patches whose radial selection functions must be treated separately. In particular, its impact is likely to be large in analyses focusing on large scales, e.g. dedicated to primordial non-Gaussianity [34, 35].

Applying our modelling of the radial integral constraint to a RSD analysis, the bias induced on cosmological parameters was shown to be successfully reduced. Though more statistics would be required to push our validation further, the potentially remaining bias is already below the expected statistical uncertainty of future galaxy surveys.

As a further application, we showed that we could similarly apply an angular integral constraint to help mitigating unknown angular systematics. This angular integral constraint can be combined to the radial one, as required if the radial selection function is also estimated from the data. Though ony additive systematics can be fully accounted for, our scheme performs well enough with multiplicative systematics on a BOSS CMASS-like survey.

Further developments may concern the use of a fully analytic selection function, as it can often be split into radial and angular components. This would enable a faster estimation of the window functions required in our analysis.

We noted that a potential bias can emerge when the selection function used in the window function calculations for integral constraint corrections is estimated from the data itself. We suggested a way to estimate this bias, which remained very subdominant in our analysis. A workaround for future analyses may consist in predicting the radial selection function from first principles (e.g. from the luminosity function), without relying on the observed data, while marginalising on possible unknowns in the cosmological fits.

Acknowledgments

We are very grateful to the anonymous referee for their useful comments. We thank Etienne Burtin, Will J. Percival and Ashley J. Ross for encouraging us in our study of the radial integral constraint. We are grateful to Florian Beutler and Emanuele Castorina for discussions about wide-angle corrections. Power spectrum multipoles were measured using nbodykit [25]. Minimisations were performed with the MINUIT algorithm [32] via the iminuit [33] Python interface. Some of our results required the use of the HEALPix package [37]. We acknowledge support from the P2IO LabEx (ANR-10-LABX-0038) in the framework ”Investissements d’Avenir” (ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR, France).

Appendix A FKP weights

We apply FKP weights [8] to both mock data and randoms. They are calculated according to:

[TABLE]

using $P_{0}=20000\,(\mathrm{Mpc}/h)^{3}$, which should be representative of the value of the power spectrum at the scales of interest for the clustering analysis. These weights require an estimation of the true redshift density $n(z)$ in absence of clustering. If $n(z)$ is computed from the data itself, using narrow bins in $z$, clustering overdensities may be smoothed out along the line-of-sight.

As mentioned in section 4.1, the redshift density $n(z)$ is computed in bins $\Delta z=0.005$. Figure 13 shows power spectrum measurements obtained with the binned scheme, with different $n(z)$ estimations. Measuring $n(z)$ from the full CMASS sample induces a negligible bias with respect to taking the true $n(z)$. However, measuring $n(z)$ separately in each of the $6$ chunks dividing the CMASS footprint, as would be natural to do, leads to an additional bias on all scales of all multipoles. Though the density $n(z)$ also enters the power spectrum normalisation (see eq. (3.3)), the main effect comes from the smoothing of the clustering along the line-of-sight, as can be seen from the loss of power in the quadrupole and hexadecapole. This effect is multiplicative and would be difficult to take into account.

A simple way to prevent FKP weights from biasing clustering measurements would be to fit a simple spline to the redshift distribution in wide redshift bins to reduce the correlation between $w_{\mathrm{FKP}}$ and the density field. As stated in section 4.1, we choose for simplicity to estimate $n(z)$ from each data realisation, using the full CMASS sample (instead of the $6$ different chunks), making the bias from $w_{\mathrm{FKP}}$ almost invisible.

Appendix B Modelling wide-angle effects

Wide-angle effects arise when a single line-of-sight must be chosen for a galaxy pair [41]. They can be accounted for by expanding both correlation and survey window functions in powers of $x=s/d$, with $s$ the separation of a pair of galaxies and $d$ its comoving distance to the observer. To allow a considerable reduction of computation time, the Yamamoto power spectrum estimator [9] uses one galaxy of a pair as line-of-sight: the so-called end-point line-of-sight (see figure 14). As shown in [40], this choice leads to detectable wide-angle effects, mainly in odd multipoles. The window-convolved power spectrum multipoles with respect to the end-point ($\mathrm{ep}$) line-of-sight can be written [40]:

[TABLE]

where:777Note that we now use the same triangle definition as in [41] (see figure 14) $\mathbf{s}=\mathbf{r}_{1}-\mathbf{r}_{2}$, which explains the minus signs in the selection functions.

[TABLE]

$\xi_{p}^{\mathrm{ep},(n)}$ can be predicted within linear pertubation theory [41, 40].

Similar expansions can be derived for $\tilde{IC}_{\ell}^{i,j}(k),\,(i,j)\in\{(\delta,\mathrm{ic}),(\mathrm{ic},\delta),(\mathrm{ic},\mathrm{ic})\}$. The first cross-term $(\delta,\mathrm{ic})$ reads:

[TABLE]

with:

[TABLE]

Using the Rayleigh plane wave expansion:

[TABLE]

and

[TABLE]

one gets:

[TABLE]

with $\mathbf{s}=\mathbf{r}_{1}-\mathbf{r}_{2}$. $\tilde{IC}_{\ell}^{\delta,\mathrm{ic}}(k)$ is simply equal to the Hankel transform of a term similar to eq. (2.17):

[TABLE]

with $(i,j)=(\delta,\mathrm{ic})$, if we define:

[TABLE]

The second cross-term $(\mathrm{ic},\delta)$ is not trivially equal to the first one for even multipoles. Indeed:

[TABLE]

gives:

[TABLE]

with $\mathbf{s}=\mathbf{r}_{1}-\mathbf{r}_{2}$, $\mathbf{\Delta}=\mathbf{r}_{2}-\mathbf{r}_{3}$. Taking the opposite $\mathbf{s}\rightarrow-\mathbf{s}$ and defining:

[TABLE]

we obtain the integral constraint correction eq. (B.8), with $(i,j)=(\mathrm{ic},\delta)$. Let us move to the last term $(\mathrm{ic},\mathrm{ic})$:

[TABLE]

gives:

[TABLE]

with $\mathbf{s}=\mathbf{r}_{1}-\mathbf{r}_{2}$, $\mathbf{\Delta}=\mathbf{r}_{3}-\mathbf{r}_{4}$. This is eq. (B.8) with $(i,j)=(\mathrm{ic},\mathrm{ic})$ if we define:

[TABLE]

In particular, results for the global integral constraint are obtained by replacing $\mathrm{ic}\rightarrow\mathrm{glo}$ and $\epsilon_{\mathrm{ic}}\rightarrow 1$. Then, eq. (2.22) is recovered at any order $n$ if one takes $\mathcal{W}_{\ell}^{\delta,\delta,(0)}(s)$ for $\mathcal{W}_{\ell}^{\delta,\delta}(s)$, as stressed out by [40].

Contrary to the line-of-sight definition in $\mathcal{W}_{\ell}^{\delta,\delta,\mathrm{ep},(n)}(s)$ (and the first line-of-sight in $\mathcal{W}_{\ell p}^{i,j,(n)}(s,\Delta)$), which should be the same as in the power spectrum (or correlation function) estimator, the line-of-sight connecting $\xi_{p}^{\mathrm{ep},(n)}(\Delta)$ to $\mathcal{W}_{\ell p}^{i,j,(n)}(s,\Delta)$ is a purely practical choice. Our calculations use the end-point line-of-sight, but the derivation with any line-of-sight $\mathbf{d}$ is straightforward by replacing $x^{-n}$ in eq. (B.9), (B.12) and (B.15) by $d^{-n}$, changing the arguments of Legendre polynomials accordingly and taking the corresponding $\xi_{p}^{\mathbf{d},(n)}(\Delta)$ in eq. (B.8). In particular, taking the mid-point line-of-sight makes first-order ($n=1$) wide-angle corrections vanish.

For clarity, we have used the end-point line-of-sight as second line-of-sight in $\mathcal{W}_{\ell p}^{i,j,(n)}(s,\Delta)$ throughout the paper. With this choice, by definition, the integral constraint corrections completely cancel with the convolved power spectrum monopole $P_{0}^{\mathrm{c}}(k)$ on large scales. This would only be asymptotically true (as the wide-angle correction order $n\rightarrow\infty$) if a different line-of-sight definition was used in $\mathcal{W}_{\ell}^{\delta,\delta,(n)}(s)$ and $\mathcal{W}_{\ell p}^{i,j,(n)}(s,\Delta)$. We show wide-angle contributions to the convolved power spectrum multipoles and the radial integral constraint in figure 15. They are significant for $k\lesssim 10^{-2}\,h/\mathrm{Mpc}$. However, they remain small compared to the radial integral constraint correction. Then, for simplicity, we do not include these wide-angle corrections ($n\geq 1$) in our analysis. However, wide-angle contributions may dominate over radial integral constraint corrections for a large survey with a constant radial selection function.