# The VIMOS Public Extragalactic Redshift Survey (VIPERS): Reconstruction of the redshift-space galaxy density field

B. R. Granett, E. Branchini, L. Guzzo, U. Abbas, C. Adami, S. Arnouts,

Julien Bel, M. Bolzonella, D. Bottini, A. Cappi, et al.

► To cite this version:

B. R. Granett, E. Branchini, L. Guzzo, U. Abbas, C. Adami, et al., The VIMOS Public Extragalactic Redshift Survey (VIPERS): Reconstruction of the redshift-space galaxy density field. Astronomy & Astrophysics - A&A, 2015, 583, pp.A61. 10.1051/0004-6361/201526330. hal-01260939

HAL Id: hal-01260939

https://hal.science/hal-01260939v1

Submitted on 4 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The VIMOS Public Extragalactic Redshift Survey
## Reconstruction of the redshift-space galaxy density field*

B. R. Granett¹,**, E. Branchini²,³,⁴, L. Guzzo¹,⁵, U. Abbas⁶, C. Adami⁷, S. Arnouts⁸, J. Bel¹, M. Bolzonella⁹,
D. Bottini¹⁰, A. Cappi⁹,¹¹, J. Coupon¹², O. Cucciati¹³,⁹, I. Davidzon⁹,¹³, G. De Lucia¹⁴, S. de la Torre⁷, A. Fritz¹⁰,
P. Franzetti¹⁰, M. Fumana¹⁰, B. Garilli¹⁰,⁷, O. Ilbert⁷, A. Iovino¹, J. Krywult¹⁵, V. Le Brun⁷, O. Le Fèvre⁷,
D. Maccagni¹⁰, K. Małek¹⁶, F. Marulli¹³,¹⁷,⁹, H. J. McCracken¹⁸, M. Polletta¹⁰, A. Pollo¹⁹,¹⁶, M. Scodeggio¹⁰,
L. A. M. Tasca⁷, R. Tojeiro²⁰, D. Vergani²¹,⁹, A. Zanichelli²², A. Burden²³, C. Di Porto⁹, A. Marchetti¹⁰, C. Marinoni²⁴,
Y. Mellier¹⁸, T. Moutard⁴, L. Moscardini¹³,¹⁷,⁹, R. C. Nichol²³, J. A. Peacock²⁵, W. J. Percival²³, and G. Zamorani⁹

*(Affiliations can be found after the references)*

Received 16 April 2015 / Accepted 19 June 2015

## ABSTRACT

**Aims.** Using the VIMOS Public Extragalactic Redshift Survey (VIPERS) we aim to jointly estimate the key parameters that describe the galaxy density field and its spatial correlations in redshift space.

**Methods.** We use the Bayesian formalism to jointly reconstruct the redshift-space galaxy density field, power spectrum, galaxy bias and galaxy luminosity function given the observations and survey selection function. The high-dimensional posterior distribution is explored using the Wiener filter within a Gibbs sampler. We validate the analysis using simulated catalogues and apply it to VIPERS data taking into consideration the inhomogeneous selection function.

**Results.** We present joint constraints on the anisotropic power spectrum, and the bias and number density of red and blue galaxy classes in luminosity and redshift bins as well as the measurement covariances of these quantities. We find that the inferred galaxy bias and number density parameters are strongly correlated although they are only weakly correlated with the galaxy power spectrum. The power spectrum and redshift-space distortion parameters are in agreement with previous VIPERS results with the value of the growth rate $f\sigma_8 = 0.38$ with 18% uncertainty at redshift 0.7.

**Key words.** large-scale structure of Universe – cosmology: observations – galaxies: statistics – cosmological parameters – methods: statistical – methods: data analysis

## 1. Introduction

The distribution of galaxies on large scales provides a fundamental test of the cosmological model. In the standard picture, galaxies trace an underlying matter density field and the statistical properties of this field such as its power spectrum and higher order moments are given by the theory (Peebles 1980). This clear view is confounded, however, by the sparse distributions of luminous galaxies mapped by surveys (Lahav & Suto 2004). Galaxies are complex systems; they are biased tracers of the non-linear matter field and their clustering strength depends on their properties and formation histories (Davis & Geller 1976; Blanton et al. 2005; Kaiser 1984; Bardeen et al. 1986; Mo & White 1996). Furthermore, their redshift gives a distorted view of distance, which is affected by coherent and random velocities (Kaiser 1987).

Upcoming galaxy surveys will be sensitive to subtleties in these trends requiring increasingly sophisticated modelling and numerical simulations to interpret the galaxy distribution in detail (The Dark Energy Survey Collaboration 2005; LSST Dark Energy Science Collaboration 2012; Levi et al. 2013; Laureijs et al. 2011). At the same time, these surveys will be sufficiently

large to be limited by minute selection effects that systematically and significantly alter the observed distribution of galaxies. Instrumental and observational artefacts can masquerade as genuine astrophysical effects and vice-versa. Thus the analyses will need to track a large number of instrumental and astrophysical parameters and be able to characterise the covariances between them. Reliable error estimation will require incorporating the set of both systematic and random uncertainties. The stakes are high as experiments promise highly precise constraints on the nature of gravity, dark energy and dark matter (Amendola et al. 2013).

Together, the physical and instrumental models compose the total data model. Given the large number of parameters, the Bayesian approach is often preferred over the frequentist one to set joint constraints on the relevant physical quantities (Trotta 2008). At the heart of this approach is the Bayes theorem which dictates a recipe for translating a set of observations into constraints on model parameters. Of fundamental importance is the incorporation of any prior knowledge of these parameters. This framework provides a natural means to jointly constrain physical parameters of interest while marginalising over a set of nuisance parameters. A paradigmatic example is the analysis of cosmic microwave background data (Jewell et al. 2004; Eriksen et al. 2004; Wandelt et al. 2004) as demonstrated through the ESA Planck mission results (Planck Collaboration I. 2015).

The Wiener filter is the first example of the application of Bayesian reconstruction techniques to galaxy surveys. The

<footnote>* Appendices are available in electronic form at http://www.aanda.org</footnote>

<footnote>** Corresponding author: Benjamin R. Granett, e-mail: ben.granett@brera.inaf.it</footnote>
Wiener solution corresponds to the maximum a posteriori solution given a Gaussian likelihood and prior. In general, for a signal contaminated by noise, the Wiener filter gives a reconstruction of the true signal with the minimum residual variance (Rybicki & Press 1992). This is also true for non-Gaussian signal and noise sources, and for this reason, since the galaxy field is not Gaussian (it is thought to tend toward Gaussianity on very large scales), the Wiener filter has seen significant use in reconstructing the density field from galaxy surveys (e.g. Lahav et al. 1994) and in particular to predict large-scale structures behind the Galactic plane (Zaroubi et al. 1995).

Wiener filtering is comparable to other adaptive density reconstruction techniques such as Delaunay tessellations (Schaap & van de Weygaert 2000), although Wiener filtering offers the advantage of naturally accounting for a complex survey selection function with inhomogeneous sampling. Examples of applications of Wiener filtering to galaxy surveys include the Two-degree Field Galaxy Redshift Survey (2dFGRS) in which the Wiener filter was used to identify galaxy clusters and voids (Erdogdu et al. 2004). Kitaura et al. (2009) present a Wiener density field reconstruction of the Sloan Digital Sky Survey (SDSS) main sample. Applied to the VIMOS Extragalactic Redshift Survey (VIPERS; Guzzo et al. 2014), the Wiener filter can naturally account for inhomogeneous sampling and survey gaps. In a comparison study of different density field estimators for VIPERS Cucciati et al. (2014) find that the Wiener filter performs well although it over-smooths the field in low-density environments affecting cell-count statistics.

Physically motivated probability distribution functions have been developed to improve on the Wiener filter and obtain unbiased density field reconstructions. Kitaura et al. (2010) demonstrate in a comparison study that the use of a Poisson sampling model for the galaxy counts with a log-normal prior on the density field allows better estimation of the lowest and highest density extremes on small scales. The generalisation of the model calls for a fully non-linear solver (Jasche & Kitaura 2010). The Poisson log-normal model was used to reconstruct the density field probed by the SDSS sample (Jasche et al. 2010a).

The Gaussian likelihood has also been used to construct maximum a posteriori estimators for the galaxy power spectrum (Efstathiou & Moody 2001; Tegmark et al. 2002; Pope et al. 2004; Granett et al. 2012). For the galaxy luminosity function estimates, maximum likelihood techniques have also enjoyed significant use (Ilbert et al. 2005; Blanton et al. 2003; Efstathiou et al. 1988).

Gaussian likelihood methods have only recently been developed to jointly infer the density field, power spectrum and luminosity function from galaxy surveys (Kitaura & Enßlin 2008; Enßlin et al. 2009). The first application to the Sloan Digital Sky Survey was demonstrated by Jasche et al. (2010b) who utilise a Gaussian likelihood and prior to jointly estimate the underlying galaxy field and power spectrum. This work was generalised to simultaneously estimate the linear galaxy bias and luminosity function (Jasche & Wandelt 2013b). Ata et al. (2015) further model a scale-dependent and stochastic galaxy bias using the log-normal Poisson model. The methodology has also been developed for photometric redshift surveys (Jasche & Wandelt 2012) as first proposed by Kitaura & Enßlin (2008).

The peculiar velocities of galaxies distort the density field inferred from redshift surveys (Kaiser 1987). The average effect may be accounted for by a convolution operation (Landy & Szalay 2002). This serves on large scales where the density field and velocity field may be inferred in a self-consistent manner (Kitaura et al. 2012b; Nusser & Davis 1994).

The full description of the galaxy field requires consideration of the higher order moments and depends on the physics of structure formation. Thus reconstruction methods have been developed that incorporate physical models based on second-order perturbation theory (Kitaura et al. 2012a; Kitaura 2013; Jasche & Wandelt 2013a; Jasche et al. 2015) or approximate n-body methods such as the particle-mesh code (Wang et al. 2014). Reconstructions of the local Universe have been used in novel ways, including to estimate the bias in the Hubble constant due to cosmic flows (Hess & Kitaura 2014).

In this work we carry out a Bayesian analysis of the VIMOS Extragalactic Redshift Survey (VIPERS; Guzzo et al. 2014). Our goal is to jointly estimate the key statistics including the matter power spectrum, galaxy biasing function and galaxy luminosity function. Our strategy is, given the observed number density of galaxies in the survey as a function of position N(RA, Dec, z), to compute the conditional probability distribution for the parameters, written schematically as: the matter over-density field δ, galaxy mean number density ¯N, galaxy bias b and the two-point correlation function S. The conditional probability distribution or posterior may be decomposed using Bayes theorem:

$$p(\delta, \bar{N}, b, S | N) \propto p(N | \delta, \bar{N}, b, S) p(\delta, \bar{N}, b, S). \quad (1)$$

The first and second factors on the right hand side are the data likelihood and the parameter prior. We will account for observational systematics such as the survey selection function in the data model, but we will not propagate their uncertainties. For VIPERS, the uncertainties in the selection function are subdominant compared with the statistical errors and so the inclusion of the uncertainties will be reserved for future work. We adopt the Gibbs sampling algorithm to sample from the posterior distribution (Marin & Robert 2007). With this approach the complex joint probability distribution is broken up in a number of simpler, individual conditional distributions. Sampling these distributions allows us to build up a Markov chain that rapidly converges to the joint distribution.

VIPERS has mapped the galaxy field to redshift 1 with unprecedented fidelity (Guzzo et al. 2014; Garilli et al. 2014). So far, VIPERS data have been used to constrain the growth rate of structure through the shape of the redshift-space galaxy correlation function (de la Torre et al. 2013). The cosmological interpretation of the galaxy power spectrum monopole has been presented by Rota et al. (in prep.). The biasing function that links the galaxy and dark matter density has been estimated using the luminosity-dependent correlation function (Marulli et al. 2013) and the shape of the one-point probability distribution function of galaxy counts in cells (Di Porto et al. 2014). Moreover VIPERS has tightened the constraints on the galaxy luminosity and stellar mass functions (Davidzon et al. 2013; Fritz et al. 2014). These measurements firmly anchor models of galaxy formation at redshift 1.

Beyond the one- and two-point statistics of the galaxy field, galaxies are organised into a cosmic web of knots, filaments and walls that surround large empty voids. In VIPERS the higher order moments of the galaxy counts in cells distribution function have been measured (Cappi et al. 2015). While ongoing efforts are being made to measure the morphologies of cosmic structures. A catalogue of voids has been constructed with VIPERS (Micheletti et al. 2014; Hawken et al., in prep.). More generally, the Minkowski functionals can be used to characterise the topology of large-scale structure as a function of scale (Schimd, in prep.). These measurements typically call for precise reconstructions of the density field corrected for observational systematics such as survey gaps and inhomogeneous sampling.
This work extends previous analyses by considering the joint
distribution of galaxy luminosity, colour and clustering bias with
the spatial power spectrum and density field. We begin in Sect. 2
with an overview of VIPERS and the parameterisation of the
selection function. The data model is described in Sect. 3, and
the method is outlined in Sect. 4. In Sect. 5 we present the con-
straints from the VIPERS data.

We assume the following fiducial cosmology $\Omega_m = 0.27$, $\Omega_b = 0.0469$, $\Omega_\Lambda = 0.73$, $n_s = 0.95$, $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}$, $\sigma_8 = 0.80$. This coincides with the MultiDark simulation run (Prada et al. 2012) that was used to construct the mock VIPERS catalogues. Magnitudes are in the AB system unless noted. The absolute magnitudes used were computed under a flat cosmology with $\Omega_m = 0.30$; however we transform all magnitudes to the $\Omega_m = 0.27$ cosmology.

## 2. VIPERS

The VIMOS Public Extragalactic Redshift Survey (VIPERS) is an ESO programme on VLT (European Southern Observatory – Very Large Telescope; Guzzo et al. 2014; Garilli et al. 2014). The survey targets galaxies for medium resolution spectroscopy using VIMOS (VIsible Multi-Object Spectrograph; Le Fèvre et al. 2003) within two regions of the W1 and W4 fields of the CFHTLS-Wide Survey (Canada-France-Hawaii Telescope Legacy Wide; Cuillandre et al. 2012). Targets are chosen based upon colour selection to be in the redshift range 0.5–1.2. The final expected sky coverage of VIPERS is 24 deg².

For each galaxy, the *B*-band rest-frame magnitude was estimated following the Spectral Energy Distribution (SED) fitting method described in Davidzon et al. (2013) and adopted to define volume limited samples. The choice of *B*-band rest-frame is natural, corresponding to the observed *I*-band at redshift ~0.8. We derived *K*-corrections from the best-fitting SED templates using all available photometry including near-UV, optical, and near-infrared.

### 2.1. Sample selection

This analysis is based on the VIPERS v5 internal data release which represents 77% of the final survey. We select sources from the VIPERS catalogue in the redshift range 0.6–1.0 with redshift confidence >95% (redshift flags 2,3,4,9). The redshift distribution is shown in Fig. 1. The sources have estimated rest frame Buser *B*-band magnitudes and *U* − *V* colours defined in the Johnson-Cousins-Kron system as described by Fritz et al. (2014). The total number of sources used in the analysis is 36928. We construct subsamples of galaxies in bins of redshift (Δz = 0.1), luminosity (ΔMB = 0.5), and colour as illustrated in Fig. 2. We separate red and blue galaxy classes using the cut defined by Fritz et al. (2014) at (M*U* − M*V*)Vega + 0.25z = 1.1 in the Vega system. In total there are 37 bins, 19 for blue and 18 for red galaxy subsamples, including those that are not complete, see the discussion in Sect. 3.3.

### 2.2. Survey completeness

The VIPERS survey coverage is characterised by an angular mask (Guzzo et al. 2014). The mask is made up of a mosaic of VIMOS pointings, each consisting of four quadrants. Regions around bright stars and of poor photometric quality in the CFHTLS photometric catalogue have been removed.

<figure><img src="image_2.png" /></figure>

**Fig. 1.** VIPERS galaxy number density of the v5 internal release sample. The curves show the effect of the completeness corrections including the spectroscopic success rate (SSR) and target sampling rate (TSR). Our analysis uses the redshift range 0.6–1.0 (vertical lines).

<figure><img src="image_3.png" /></figure>

**Fig. 2.** Subsamples of galaxies in colour, absolute magnitude and redshift bins. Left: the absolute magnitude-colour plane. The histogram at top shows the distribution of colour. The sample is divided into blue and red classes following Fritz et al. (2014) at $M_U - M_V + 0.25z = 1.1$. Right: the sample is further binned by redshift and absolute magnitude. The luminosity bins account for the mean evolutionary trend. The faintest luminosity bins are not volume limited and the thick blue and dashed red curves show the limiting magnitudes for blue and red galaxy classes.

Within an observed quadrant there are many factors including the intrinsic source properties, instrument response and observing conditions that determine the final selection function (Garilli et al. 2014). The fraction of sources out of the parent photometric sample that are targeted for spectroscopy is referred to as the target sampling rate (TSR). Among the targeted sources, not all will give a reliable redshift measurement. We refer to this fraction as the spectroscopic success rate (SSR). The sampling rate is the product of the TSR and SSR: $r = r_{\text{TSR}}r_{\text{SSR}}$.

The arrangement of slits in VIMOS is strongly constrained since the spectra cannot overlap on the imaging plane (Bottini et al. 2005). In VIPERS, the result is that the number of targeted sources in a pointing is approximately constant, damping the galaxy clustering signal on both small and large scales (Pollo et al. 2005; de la Torre et al. 2013).

As described in Sect. 3, in our data model we bin the galaxies onto a cubic grid, so it is only necessary to estimate the TSR correction on the scale of the cubic cell. For 5 $h^{-1}$ Mpc cells, this corresponds to 10 arcmin at $z = 0.7$, which is larger than the VIMOS quadrants. The TSR is estimated on a fine grid as the fraction of targets out of the parent sample within a 3 arcmin circular aperture: $N_{\text{target}}/N_{\text{parent}}$. The fine grid is then down-sampled to determine the average TSR in each grid cell. The colours in Fig. 3 indicate the TSR measurements as a function of angular position (at the positions of observed galaxies).

The spectroscopic success rate is primarily correlated to the conditions at the time of observation and so varies with pointing.
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Distribution of VIPERS targets on the sky is plotted for the two fields W1 and W4. The points are coloured according to the target sampling rate (TSR) which is defined as the ratio of the number of targeted galaxies in a patch of sky over the total number available in the parent photometric catalogue. In this work we estimate the TSR within a circular aperture with radius 3 arcmin. The TSR depends on the projected density of targets on the sky. TSR is higher in low-density fields with few potential targets while in high-density fields slit positioning constraints severely limit the number of sources that may be targeted. The inset histogram shows the cumulative fraction of targets with TSR below a given value. The median TSR over the survey is 48%.

<figure><img src="image_5.png" /></figure>

**Fig. 4.** The spectroscopic success rate (SSR) quantifies the fraction of targets for which the redshift could be measured with >95% confidence. On the left we show the mean SSR of each pointing divided by quadrant (7×8 arcmin). The inset histogram gives the cumulative number of quadrants with SSR below the given value. On the right, we show the SSR as a function of $i_{AB}$ magnitude (solid curves) with an analytic fit (dashed curves). The sample is divided based on the overall quality of the quadrants quantified by SSR. The quadrants are ranked by mean SSR and the curves are computed for each decile. The range of SSR and the number of sources in the bin are given in the figure legend.

For a particular source the SSR depends on the apparent flux as well as the spectral features that are available to make the redshift measurement. We find that the primary contribution comes from the apparent flux, and we quantify the mean SSR, defined as $N_{measured}/N_{target}$ in each quadrant, as a function of the *i* band magnitude (Guzzo et al. 2014; Garilli et al. 2014). The degradation is most severe in poor observing conditions, so we compute the SSR separately according to the quality of the quadrant. We rank the quadrants based on mean SSR and compute separately $r_{SSR}(m_i)$ in each decile. The SSR is fit with an analytic form: $r_{SSR}(m) = a(1 - e^{c(m-b)})$. In the right panel of Fig. 4 we see that the SSR depends strongly on the quadrant quality. For the top 10% of quadrants (shown by the red curve in Fig. 4), the SSR remains >90%, but it drops quickly as quadrant quality falls. What is important to note here is that the shape of the $r_{SSR}(i_{AB})$ curves changes as a function of quadrant quality.

The effect of the weights on the redshift distribution is shown in Fig. 1.

## 2.3. Mock catalogues

We use a set of simulated (mock) galaxy catalogues constructed to match the VIPERS observing strategy. The catalogues are built on the MultiDark *N*-body simulation (Prada et al. 2012) using the Halo Occupation Distribution (HOD) technique. Details of the construction may be found in (de la Torre et al. 2013; de la Torre & Peacock 2013).

Galaxies were added to the dark matter halos in the simulation according to a luminosity dependent HOD model. The correlation function and number counts in luminosity bins were set to match measurements made at $0.5 < z < 1.2$ in CFHTLS, VVDS and earlier releases of VIPERS. Each mock galaxy is
characterised by its angular coordinate, comoving distance, observed redshift including its errors and an absolute magnitude in the B band.

We partition the mock catalogues into bins of redshift and luminosity, but not in colour, as we do for the VPERS data. The step size in redshift and luminosity are $\Delta z = 0.1$ and $\Delta M_B = 0.2$ over the redshift range $0.6 < z < 1.0$. We have 46 bins including those bins that are not complete owing to the apparent flux limit.

The mock galaxies match the number density of the VIPERS observations, but do not include the slit placement constraints that we correct in the data with the TSR weights. The catalogues do not simulate the spectroscopic sampling rate.

## 3. Data model

### 3.1. Galaxy number counts

We overlay a three-dimensional cartesian grid on the survey. The number of galaxies in a given sample observed within a cell indexed by $i$ is related to an underlying continuous galaxy density field $\delta_G$ by

$$N_i = \bar{N}w_i(1 + \delta_{G,i}) + \epsilon_i, \quad (2)$$

where $w_i$ is the spatial selection function and $\bar{N}_i$ is the mean density providing the normalisation. The stochastic nature of galaxy counts is captured by the random variable $\epsilon_i$ and is dominated by Poisson noise except in the highest density peaks (Di Porto et al. 2014). The cells are defined in comoving redshift-space coordinates and we adopt a fiducial cosmology to define the relationship between redshift and comoving distance.

The selection function $w$ in Eq. (2) gives the likelihood of observing a galaxy at a given grid point. It accounts for the angular geometry of the survey, sampling rate and redshift distribution. In this analysis we separate the angular and line-of-sight components. As described in Sect. 5.1, the angular dependence is determined from the survey mask and TSR while the redshift distribution is computed assuming the luminosity function and apparent flux limit for the given subsample of galaxies.

The expected number of galaxies in a cell is given by the product of $\bar{N}$ and the selection function, giving

$$\langle N_i \rangle = \bar{N}w_i. \quad (3)$$

The selection function described can account for spatial variations but cannot describe sampling dependencies on galaxy type or apparent flux. For VIPERS data we will up-weight galaxies based on the inverse SSR depending on quadrant and apparent $i$ band magnitude. These weights are only indirectly correlated with the density field so they result in an amplification of the shot noise level. The SSR weight of a galaxy is $w_{SSR}$ and the noise amplification factor is $\alpha = \langle w_{SSR} \rangle$ averaged over all galaxies in the subsample. Therefore, in Eq. (2) $N_i$ represents the weighted count of galaxies and, consequently, the variance of the stochastic term, $\sigma_{\epsilon}^2$, is boosted to $\alpha\sigma_{\epsilon}^2$.

In this analysis we discretise the galaxy field onto a coarse spatial grid as well as onto a finite grid of Fourier modes. This process introduces an error in the density field arising from the aliasing of structures: small-scale structures with spatial frequencies higher than the Nyquist frequency become imprinted on larger scales (Hockney & Eastwood 1988). The effect may be corrected for in the power spectrum by assuming the spectral shape above the Nyquist frequency (Jing 2005). However, to accurately reconstruct the density field without making such assumptions, we may use a mass-assignment scheme or anti-aliasing filter discussed in Appendix A. The smoothing effect

of mass-assignment schemes introduces a convolution in Eq. (2) which invalidates the simple count model. An alternative is to use the super-sampling method of Jasche et al. (2009) that approximates the ideal anti-aliasing filter and does not damp small-scale power. Since we desire a compact window in both configuration and Fourier space, we adapt this technique with a soft cut as described in Appendix A. This approach reduces the aliased signal to the level of the triangle-shaped-cell scheme while preserving small-scale power to $k \sim 0.7k_{\text{Nyquist}}$.

The convolution introduced by the anti-aliasing filter modifies the noise properties such that the Poisson expectation, $\sigma_i^2 = \bar{N}_i w_i$, cannot be assumed. Instead we use a re-scaled Poisson variance characterised by the factor $v_i \equiv \sigma_i^2 / \bar{N}$. As described in Sect. 5, the factor may be estimated in a Monte Carlo fashion given the mask and anti-aliasing filter.

A cell partially cut by the mask will become strongly coupled with its neighbour through the anti-aliasing filter. In practice we neglect the additional off-diagonal cell-cell contributions in the noise covariance matrix. However, we found that it is necessary to regularise the noise matrix by increasing the noise level to achieve stable results. We set a lower limit on the cell variance through a parameter $v_{\text{threshold}}$ such that the scaled shot noise has a floor set by $v_i \leftarrow \max(v_{\text{threshold}}, v_i)$. The procedure is identically applied both to mock and real data.

### 3.2. Galaxy bias

We assume a constant linear biasing model such that

$$\delta_{G,i} = bD(z_i)\delta_i, \quad (4)$$

where $z_i$ is the redshift of the cell indexed by $i$. The bias factor $b$ depends on the luminosity and colour of the galaxy subsample. We give explicitly the growth of matter fluctuations with time according to the linear growth factor $D(z) \equiv D(z)/D(z_{\text{ref}})$ with $\delta(z) = \delta(z_{\text{ref}})D(z)$. For the VIPERS sample the reference redshift is set to $z_{\text{ref}} = 0.7$. Since VIPERS covers an extended redshift range this factor brings large-scale density modes to a common epoch.

### 3.3. Number density

The number density of galaxies in a luminosity bin is given by the integral of the luminosity function:

$$n(z) = \int_{M_{\text{bright}}}^{M_{\text{faint}}(z)} n(M, z) dM. \quad (5)$$

We parameterise the luminosity function using the Schechter function (Schechter 1976) in terms of magnitudes and the parameters $(\phi_{\star}, M_{\star}, \alpha)$:

$$n(M) = 0.4 \ln 10 \phi_{\star} (10^{0.4(M_{\star} - M)})^{\alpha+1} \exp(10^{0.4(M_{\star} - M)}). \quad (6)$$

The characteristic magnitude evolves as $M_{\star}(z) = M_{\star}(0) + \mathcal{E}z$ with $\mathcal{E} \approx -1$ for red galaxies in VIPERS (Fritz et al. 2014) confirming the findings of previous studies at moderate redshift (Ilbert et al. 2005; Zucca et al. 2009).

The number density observed is further reduced by the survey completeness. In VIPERS, galaxies are targeted to an apparent magnitude limit of $m_{\text{lim}} = 22.5$ in the $i_{AB}$ photometric band. This sets an absolute magnitude limit for a given class of galaxy

$$M_{\text{limit}}(z) = m_{\text{lim}} - D_{\text{m}}(z) - \mathcal{K}(z), \quad (7)$$
where $D_m(z)$ is the distance modulus which depends only on the background cosmology and $\mathcal{K}$ is the $K$-correction term which depends on the particular type of galaxy targeted. The absolute magnitudes were computed for each galaxy by fitting spectral energy distribution templates to broadband photometry as described by Davidzon et al. (2013) and from the absolute magnitudes we infer the $K$-corrections. We parameterise the trend of $K$-correction with redshift as $\mathcal{K}(z) = \mathcal{K}_0 + \Delta_K(z)$. The value of $\mathcal{K}_0$ is estimated from the median value of galaxies within a given subsample while $\Delta_K(z)$ is a polynomial fit with the following coefficients, different for red and blue galaxy types:

$$ \Delta_{K,\text{blue}}(z) = 1.784(z - 0.7)^2 + 0.440(z - 0.7) - 0.678 \quad (8) $$

$$ \Delta_{K,\text{red}}(z) = 2.144(z - 0.7)^2 + 1.745(z - 0.7) - 0.720. \quad (9) $$

Using this parameterisation the inferred magnitude limits corresponding to 50% completeness are indicated in Fig. 2 for the blue and red samples by the solid and dashed lines. For mock samples, the $K$-correction term and its evolution are fixed to $\mathcal{K}_{\text{mock}} = z - 1.3$.

With these ingredients we model the redshift distribution of each galaxy subsample by integrating the luminosity function with Eqs. (7) and (5). We leave the mean density of each subsample free to set the normalisation of the redshift distribution. We then take the shape given by the Schechter function to interpolate the luminosity function across the bin. The parameters $M_\star$ and $\alpha$ in each colour and redshift bin are fixed to the values measured in VIPERS (Fritz et al. 2014). Since the precise luminosity evolution is not known, the evolution term, $\mathcal{E}$, is allowed to vary as a function of colour and redshift. This gives a characteristic magnitude $M_\star(z) = M_\star(z_{\text{ref}}) + \mathcal{E}(z - z_{\text{ref}})$, where $z_{\text{ref}}$ is taken to be the midpoint of the redshift bin. Changing $\mathcal{E}$ modifies the shape of the redshift distribution.

### 3.4. Power spectrum

The matter power spectrum in real space $P(k) = \langle |\delta_k|^2 \rangle$ is assumed to be isotropic. Seen in redshift-space, it is distorted along the line-of-sight direction (Hamilton 1998). We model the signal on the cartesian Fourier grid as

$$ S(k, \mu; \beta, \sigma_v, \sigma_{\text{obs}}) = A \frac{(1 + \beta \mu^2)^2}{1 + k_{\text{los}}^2 \sigma_v^2} e^{-\frac{k_{\text{los}}^2 \sigma_{\text{obs}}^2}{2}} \mathcal{B}^2(k_x, k_y, k_z) P(k), \quad (10) $$

where $\mu \equiv k_{\text{los}}/k$ and $k = \sqrt{k_x^2 + k_y^2 + k_z^2}$. The line-of-sight direction is aligned with the grid such that $k_{\text{los}} = k_z$ taking the plane-parallel approximation.

The coherent motions of galaxies on large scales are described by the Kaiser (1987) factor with $\beta = f/b_g$, where the growth rate in $\Lambda$CDM is $f(z) = d \log D/d \log a$. On small scales, velocities randomise and may be modelled by an exponential pairwise velocity dispersion giving a Lorentzian profile in Fourier space which we refer to as the dispersion model (Ballinger et al. 1996). The velocity dispersion term, $\sigma_v$ in Eq. (10), has units of $h^{-1}$ Mpc. The conversion to velocity units is $H(z)/(1+z)/\sqrt{2} \approx 60.0 h \text{ Mpc}^{-1} \text{ km s}^{-1}$, which is nearly constant over the redshift range of interest. We add a Gaussian term along the line of sight to characterise redshift measurement errors where $\sigma_{\text{obs}} = \sigma_{cz}/H(z)$ and $\sigma_{cz}$ is the redshift error. For VIPERS the estimated redshift error is $\sigma_{cz} = 141(1+z) \text{ km s}^{-1}$ (Guzzo et al. 2014) and $\sigma_{\text{obs}} = 1.67 h^{-1} \text{ Mpc}$ and is nearly constant over the redshift range 0.6–1.0.

<table><caption>**Table 1.** Accounting of the free parameters in the data model.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Dimension</th></tr></thead><tbody><tr><td>Overdensity field</td><td>δ</td><td>2 × 72 × 16 × 172<br>(5 h<sup>-1</sup> Mpc cubic cells)</td></tr><tr><td>Power spectrum</td><td>P</td><td>109</td></tr><tr><td>Distortion factor</td><td>β</td><td>1</td></tr><tr><td>Velocity dispersion</td><td>σ<sub>v</sub></td><td>1</td></tr><tr><td>Galaxy bias</td><td>b</td><td>37 (19 blue, 18 red)</td></tr><tr><td>Mean number density</td><td>N̄</td><td>37</td></tr><tr><td>Luminosity evolution</td><td>ε</td><td>8</td></tr></tbody></table>

The factor $\mathcal{B}(k_x, k_y, k_z)$ accounts for the cell window function arising from the anti-aliasing filter and is given by Eq. (A.2). In this analysis the absolute amplitude of the power spectrum is not constrained. So we set the amplitude $A$ in Eq. (10) to fix $\sigma_8 = 0.8$, the variance computed on a scale of $R = 8 h^{-1}$ Mpc integrated to the Nyquist frequency.

We ignore geometric distortions arising from the choice of the fiducial cosmology (Alcock & Paczynski 1979). The resulting bias is not significant when compared with the statistical uncertainties of the VIPERS redshift-space clustering measurements (de la Torre et al. 2013). However, when carrying out a model test, we may rescale the density field and two point statistics to transform from the fiducial to the test cosmology as carried out for the VIPERS power spectrum analysis by Rota et al. (in prep.), but this is not done here.

## 4. Gibbs sampler

We present a brief overview of the Gibbs sampler. Since our implementation differs from that of Jasche & Wandelt (2013b) we provide a detailed description in Appendix B. The full parameter set introduced in the previous section is summarised in Table 1. We use the Gibbs sampling method to sample from the joint posterior of the parameter set. This is performed by iteratively drawing samples from each conditional probability distribution in the following steps (where $\leftarrow$ indicates that a sample is drawn from the given distribution):

1. generate $\delta^{s+1} \leftarrow p(\delta | \tilde{N}^s, b^s, P^s, \beta^s, \sigma_v^s, N)$;

2. generate $P^{s+1} \leftarrow p(P | \delta^{s+1}, \tilde{N}^s, b^s, \beta^s, \sigma_v^s, N)$;

3. normalise power spectrum $P^{s+1}$;

4. generate $\beta^{s+1}, \sigma_v^{s+1} \leftarrow p(\beta, \sigma_v | P^{s+1}, \delta^{s+1}, \tilde{N}^s, b^s, N)$;

5. generate $\tilde{N}^{s+1} \leftarrow p(\tilde{N} | b^{s+1}, P^{s+1}, \beta^{s+1}, \sigma_v^{s+1}, \delta^{s+1}, N)$;

6. generate $\mathcal{E}^{s+1} \leftarrow p(\mathcal{E} | \tilde{N}^s, b^{s+1}, P^{s+1}, \beta^{s+1}, \sigma_v^{s+1}, \delta^{s+1}, N)$;

7. generate $b^{s+1} \leftarrow p(b | P^{s+1}, \beta^{s+1}, \sigma_v^{s+1}, \delta^{s+1}, \tilde{N}^{s+1}, N)$.

These steps are repeated forming a Markov chain and after an initial burn-in period we can expect that the samples are representative of the joint posterior distribution.

In the first step, we sample from the conditional probability distribution for the density field in a two-stage procedure. First, the Wiener filter is used to compute the maximum a priori field $\delta_{\text{WF}}$ (Kitaura et al. 2010). The Wiener filter solution is a smoothed field that gives an underestimate of the true power. To generate a realisation of the density field a random component that is uncorrelated with the observations $\delta_{\text{random}}$ is added (Jewell et al. 2004). The final field is thus the sum $\delta = \delta_{\text{WF}} + \delta_{\text{random}}$.

After constructing a realisation of the density field, the second step is to sample the power spectrum. We put a Gaussian prior on the first bin at $k < 0.01 h \text{ Mpc}^{-1}$ setting the mean and variance to the fiducial value and sample variance expectation.
This aids the stability of the chain. A uniform prior is used for the bins at $k > 0.01 h \text{ Mpc}^{-1}$. We use two approaches to sample the power spectrum detailed in Appendix B.2. First, we draw samples from the inverse-gamma distribution, see e.g. Jasche et al. (2010b); however, this produces very small steps in the low signal-to-noise regime and so can be inefficient at small scales. Therefore, on alternative steps we carry out a Metropolis-Hastings routine to draw samples of the power spectrum according to the likelihood, Eq. (B.2). We find consistent sampling of the power spectrum using the two methods.

Since we cannot constrain the absolute normalisation of the power spectrum, we normalise to the desired value of $\sigma_8$. We then draw the redshift-space distortion parameters $\beta, \sigma_v$ which are independent of the power spectrum amplitude.

Next, we sample from the mean density conditional probability distribution for each galaxy sample which includes the evolution factor $\mathcal{E}$. Here we use a Poisson distribution, as described in Appendix B.4.

Finally we sample from the bias conditional probability distribution for each galaxy sample. This distribution is Gaussian for the bias parameter (see Appendix B.3). In this method, the bias is computed on the redshift-space grid, which in our case has a resolution of $5 h^{-1} \text{ Mpc}$. For physical interpretation it is interesting to estimate the bias averaged on larger scales. So, in estimating the bias we first down-grade the grid resolution by a factor of two, such that the bias is averaged over a scale of $10 h^{-1} \text{ Mpc}$. We impose a uniform prior for the bias values of $0.5 < b < 4$.

# 5. Application to VIPERS

## 5.1. Set-up

The data and mock catalogues are processed similarly, although the construction of galaxy subsamples differs. The mock catalogues do not include the inhomogeneous incompleteness corrected for in the data by the SSR and TSR factors. The uncertainties introduced by these corrections are negligible compared with statistical uncertainties in VIPERS.

1. The two survey fields, W1 and W4 are separately embedded into rectangular boxes. The grids have dimensions $72 \times 16 \times 172$ cells and each cubic cell has comoving size $5 h^{-1} \text{ Mpc}$. We align the grid such that at the field centre the three axes correspond to the right ascension, declination and line-of-sight directions. The co-moving coordinates are computed using the fiducial cosmological model. In the real catalogue alone, galaxies are up-weighted by the inverse SSR depending on quadrant and apparent $i$-band magnitude.

2. We compute the density on the grid using the anti-aliasing filter based on the super-sampling method proposed by Jasche et al. (2009) with a soft $k$-space cut-off as described in Appendix A.

3. The angular ($\alpha, \delta$) and radial ($z$) components of the selection function are computed separately on the grid: $w(\alpha_i, \delta_i, z_i) = w(\alpha_i, \delta_i)w(z_i)$.

4. For the angular component, we generate a uniform grid of test points that over-sample the grid by a factor of 8 and reject points outside the survey angular mask. For VIPERS data the remaining are down-sampled by the TSR. The points are then assigned to the grid points using the anti-aliasing mass-assignment scheme. The selection function $w_i$ is then given by the normalised density of test points on the grid.

5. The radial component of the selection function is estimated in bins of redshift, luminosity, and colour. We estimate the

median $K$-correction term for galaxies within each bin and use Eq. (5) to compute the unnormalised $N(z)$.

6. We estimate the generalised shot noise variance $v_i \equiv \sigma_i^2 / \bar{N}_i$ which depends on the mask through the anti-aliasing filter in Monte Carlo fashion. We generate a set of 1000 shot noise maps by distributing random points over the survey volume. For VIPERS data the points are down-sampled by TSR. We then compute the variance for each cell of the map over the 1000 realisations. To regularise the noise covariance matrix, a threshold is set $v_i = \max(v_{\text{thresh}}, v_i)$, where $v_{\text{thresh}} = 0.3$ for mocks and 0.15 for data to account for TSR.

With the galaxy number density map, selection function and noise map, we have all the components of the data model required to estimate the posterior probability distribution in Eq. (1). To sample this probability distribution we run the Gibbs sampler Markov chain for 2000 steps and, allowing for a burn-in period, begin the analysis from step 1000. The convergence properties and justification for the burn-in period are shown in Appendix C. For the VIPERS data we ran seven independent chains for 2000 steps each providing 7000 post-burn-in samples for analysis.

Taking the variance of the Markov chain gives us an internal error estimate on the parameters. The runs on mock catalogues show that the chain variance corresponds to the expected sample variance for the power spectrum. However, this is not necessarily true for the other statistics. For instance the luminosity function quantifies the distribution of observed galaxies that remains fixed in the chain. It is only indirectly dependent on the underlying density.

## 5.2. Density field

The density field taken from a single step (1500) in the Markov chain is shown in Fig. 5. It represents the application of the Wiener filter on a bias-weighted combination of the galaxy subsamples. The reconstruction is based on the redshift-space distortion model and so the resulting field is anisotropic and is characterised by effective redshift-space distortion parameters averaged over the galaxy samples.

The result of the Wiener filter is an adaptively smoothed field that extrapolates structures over the correlation length of a few megaparsecs. To build a full realisation of the structures we add a Gaussian constrained realisation that fills in the gaps and gives the full variance. The structures outside the survey boundary are generated from a random Gaussian realisation although the phases are properly aligned at the boundary. The true galaxy density field on scales of $5 h \text{ Mpc}^{-1}$ is far from Gaussian and the difference is visible by eye.

In Fig. 5 we can recognise the cosmic web of structures including knots, filaments and void regions. The structures are richest where the sampling is highest at lower redshift. At redshift $z > 0.8$ we see fewer coherent structures and the contribution from the constrained Gaussian realisation is larger. Each step of the Markov chain gives a reconstruction of the field with different realisations of noise and large-scale modes. Once the chain has passed the burn-in period (see Appendix C), we can consider these realisations to represent Gaussian perturbations around the observed galaxy field.

## 5.3. Redshift-space power spectrum

The galaxy power spectrum in redshift space is parameterised in terms of the real-space matter power spectrum, bias, and
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** VIPERS cone diagrams for the fields W1 (top) and W4 (bottom). The left panels show the redshift-space positions of observed galaxies. The marker colour indicates the blue or red colour class and the marker size scales with B-band luminosity. The depth of the slice is $10 h^{-1}$ Mpc. The orange line traces the field boundaries cut in the redshift direction at $0.6 < z < 1.0$. At right we show a slice of the density field taken from one step in the Markov chain. It represents the anisotropic Wiener reconstruction from the weighted combination of galaxy tracers. The field is filled with a constrained Gaussian realisation. The field has been smoothed with a Gaussian kernel with a full width half maximum of $10 h^{-1}$ Mpc. The colour scale gives the over-density value.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Constraints on the real-space power spectrum. The lower panel shows the relative difference with the fiducial model. The black dots give our estimates of the binned real-space power spectrum taken from the median of the Markov chain. Overplotted is the fiducial model adopted in this study (black dashed curve). We find agreement with the best-fit model using VIPERS data by Rota et al. (in prep.) (purple dot-dashed curve). The pink dashed curve is the mean of the power spectrum estimates taken from the 27 mock catalogues. We present three error estimates: the internal chain variance determined from VIPERS data (grey steps), the chain variance determined from mock catalogues (blue steps) and the variance of the individual estimates from the 27 mock catalogues (red steps). The error corridors show 70% confidence intervals.</figcaption></figure>

redshift-space distortion factors (Eq. (10)). We bin the power spectrum linearly with bin size $\Delta k = 0.01$ giving 109 bins. The redshift-space distortion parameters are fit to $k < 0.4 h \text{ Mpc}^{-1}$. This limit corresponds to $k\sigma_v \approx 1$ where we can expect the dispersion model to break down.

The Markov chain provides joint samples of the parameters. In Fig. 6 we show the median over the power spectrum chain (black dots). The confidence corridor gives the $1\sigma$ confidence interval estimated from the chain variance (grey steps).

We find good agreement with the model computed with CLASS and Halofit (black dashed) with $\Omega_m = 0.27$ (Lesgourgues 2011; Smith et al. 2003; Takahashi et al. 2012). The best-fitting model determined by Rota et al. (in prep.) has $\Omega_m = 0.272 \pm .03$ (over-plotted with purple dot-dashed curve).

On small scales $k > 0.45 h \text{ Mpc}^{-1}$ (0.75× Nyquist frequency) the power drops. This is due to neglecting correlations between cells that arise because of the anti-aliasing filter. On large scales, there is a dip in power at $k = 0.05 h \text{ Mpc}^{-1}$ seen in both mock
<table><caption>**Table 2.** Constraints on redshift-space distortion parameters.</caption><thead><tr><th></th><th>β</th><th>σ<sub>β,chain</sub></th><th>σ<sub>β,mock</sub></th><th>b<sub>eff</sub></th><th>σ<sub>b,chain</sub></th><th>σ<sub>b,mock</sub></th><th>fσ<sub>8</sub></th><th>σ<sub>f,chain</sub></th><th>σ<sub>f,mock</sub></th></tr></thead><tbody><tr><td>Mock</td><td>0.47</td><td>-0.12/+0.09</td><td>0.09</td><td>1.54</td><td>-0.04/+0.04</td><td>0.03</td><td>0.46</td><td>+0.09/-0.12</td><td>0.08</td></tr><tr><td>VIPERS</td><td>0.41</td><td>-0.08/+0.07</td><td></td><td>1.44</td><td>-0.03/+0.02</td><td></td><td>0.38</td><td>-0.07/+0.06</td><td></td></tr></tbody></table>

**Notes.** We give the 68% confidence intervals from the chains and the standard deviation among the 26 mock catalogues. The fiducial value is $f\sigma_8(z = 0.7) = 0.45$.

catalogues and data, although it biases the estimate only at the 1σ level. Scales at $k < 0.06 h Mpc^{-1}$ are only measured in the line-of-sight direction with VIPERS so the inability to reconstruct them properly without a prior constraint is not surprising.

The median values we find for the redshift-space distortion parameters are $\beta_{\text{VIPERS}} = 0.41$ and $\beta_{\text{mock}} = 0.47$. The within-chain variance is $\sigma_{\text{chain}} = (-0.12, +0.10)$ (68% confidence interval), while the scatter of the 26 mocks gives standard deviation $\sigma_{\text{mock}} = 0.09$.

We compute the growth rate through the relation

$$f\sigma_8(z) = \beta\sigma_{8,\text{galaxy}}(z), \qquad (11)$$

where $\sigma_{8,\text{galaxy}} = b_{\text{eff}}\sigma_8$. In this analysis we have fixed the amplitude of the matter power spectrum at redshift $z = 0.7$ with $\sigma_8(z = 0.7) = 0.643$ which corresponds to $\sigma_8 = 0.8$ at $z = 0$ in the fiducial cosmology.

We compute the effective galaxy bias as the number-weighted average over the galaxy samples,

$$b_{\text{eff}}^2 = \frac{\sum_{l,i} \bar{N}_l w_{l,i} b_l^2}{\sum_{l,i} \bar{N}_l w_{l,i}}, \qquad (12)$$

where the sums are over the galaxy subsamples and selection function grids.

We summarise the constraints on the growth rate in Table 2. We find $(f\sigma_8)_{\text{VIPERS}} = 0.38_{-0.07}^{+0.06}$ and $(f\sigma_8)_{\text{mock}} = 0.46_{-0.12}^{+0.09}$ at $z = 0.7$ where we quote the chain variance. The scatter between the mocks gives a standard deviation $\sigma_{f,\text{mock}} = 0.08$. Thus these constraints on the growth rate are in agreement with the VIPERS correlation function measurement by de la Torre et al. (2013). There, the error was 16% on the growth rate $f\sigma_8 = 0.48$ at $z = 0.8$. In this work we find an error of 18%. We attribute the higher error in this analysis to the fact that we marginalise over the real-space power spectrum, while in the previous analysis it was fixed to a fiducial cosmology.

The correlations between a subset of the power spectrum and redshift-space distortion parameters are shown in Fig. 7. The star symbols mark the median values of the parameters estimated from the VIPERS Markov chain while the filled contours give the 70% and 90% confidence intervals. The median value and marginalised 70% uncertainty on each parameter are labelled. We find that $f\sigma_8$ (18% relative error) is better constrained than $\beta$ (20% error). This is due to the anti-correlation between $b_{\text{eff}}$ and $\beta$ with correlation coefficient $\rho = -0.58$. These correlations arise from the specific parameterisation adopted and would be modified under a different data model.

The black dots in Fig. 7 represent the median values estimated from individual mock catalogues. We find that the value of the galaxy bias is different within the mocks ($b_{\text{eff}} = 1.55$) and VIPERS ($b_{\text{eff}} = 1.44$). The bias of the mock galaxies is determined by the luminosity-dependent HOD prescription (de la Torre et al. 2013) and so the minor difference from real data is not unexpected. Accounting for the difference in bias, we find excellent agreement between the distribution of mocks and

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Degeneracies between RSD parameters $\beta$, $\sigma_v$, effective bias and the power spectrum at $k = 0.4 h Mpc^{-1}$. The shaded regions mark the 68% and 95% confidence intervals from VIPERS chain and the star symbols mark the mean value. The points give the distribution of mean values derived from mock catalogues. The histograms along the diagonal give the marginalised distributions of each parameter chain. The filled histogram gives the distribution from the VIPERS chain, while the solid line is the distribution of mean values derived from the mock catalogues.</figcaption></figure>

the parameters estimated from real data. The similarity of the probability distribution function shapes also gives us confidence in the analysis method and error estimates. Furthermore, since the mocks do not include many of the selection effects present in the data the agreement suggests that these sources of systematic uncertainties do not influence our conclusions.

## 5.4. Colour and luminosity dependent galaxy bias

We compute the galaxy bias from the variance of the galaxy counts on the grid. However, we first down-sample the grid by a factor of two such that the bias is computed on a scale of $10 h^{-1}$ Mpc.

In Fig. 8 we show the median bias values of the Markov chain and the confidence intervals are given by the chain variance. The bias is computed in bins of redshift, luminosity, and colour. We find a colour bimodality with red galaxies more strongly biased than blue. This corresponds to the well-known galaxy morphology-density relationship that early type galaxies are predominantly found in high density environments (e.g. Cucciati et al. 2006; Dressler 1980; Davis & Geller 1976).
<figure><img src="image_9.png" /><figcaption>Fig. 8. VIPERS galaxy bias parameters in redshift, luminosity, and colour bins. A colour bimodality is seen in each redshift bin. The trend with luminosity is most striking in the lower redshift bins for both blue and red galaxies.</figcaption></figure>

Similarly, we expect to find that galaxy bias increases with galaxy luminosity since more massive and more luminous galaxies tend to form in more massive dark matter clumps (Coupon et al. 2012).

Previous studies with VIPERS data estimated the galaxy bias of the full galaxy sample as a function of luminosity and redshift. Marulli et al. (2013) measured the projected galaxy correlation function in luminosity and redshift bins to constrain the mean bias averaged over scales 5–20 h⁻¹ Mpc. Di Porto et al. (2014) modelled the counts-in-cells probability distribution function to estimate the linear bias. To compare our estimate of the galaxy bias with these previous results we construct luminosity threshold samples counting both red and blue galaxies. We cross-correlate these number density maps with the Wiener density field from the core analysis and estimate the measurement uncertainty from the chain variance. The resulting bias values are shown by the markers with error bars in Fig. 9. We find excellent agreement with the previous analyses, although our redshift bins differ. The bias values from Di Porto et al. (2014) have been taken on a scale R = 8 h⁻¹ Mpc while those of Marulli et al. (2013) are sensitive to smaller scales. The disagreement at z > 0.9 may indicate that the bias of luminous galaxies is scale dependent at high redshifts probed differently by the three studies.

### 5.5. Luminosity function

In Fig. 10 we show the derived luminosity function based on the mean galaxy number density (dots with error bars) for different galaxy types. We compare the result to the analysis from Fritz et al. (2014) based on the Sandage-Tammann-Yahil (STY, Ilbert et al. 2005; Sandage et al. 1979) estimator (dashed curves). We can expect to find a difference in the two estimates arising from how the galaxy bias is treated. The STY estimate is designed to be independent of the underlying density field under the approximation that the galaxy luminosity is uncorrelated with density. A strong luminosity dependence of the bias can systematically tilt the inferred luminosity function (Smith 2012; Cole 2011). The agreement between our analysis and the STY measurement

indicates that the VIPERS volume is large enough such that the sample variance does not significantly alter the amplitude and that any effects due to the luminosity dependence of bias are weak.

### 5.6. Parameter covariance

We estimate the covariance of the statistics with the Markov chain. Figure 11 shows the normalised correlation matrix determined in our analysis,

$$ \rho_{ij}^2 = \frac{C_{ij}}{\sigma_i \sigma_j} \qquad (13) $$

The bias and mean number density parameters are ordered first by luminosity and colour and then by redshift bin. The appearance of blocks in the matrix indicates that within redshift bins the statistics are strongly correlated. We also find that bias and mean density are anti-correlated, that is, increasing bias necessitates decreasing mean density to preserve the same fluctuation. The bias and mean density parameters are weakly correlated with the power spectrum measurement.

The bins of the power spectrum (spacing Δk = 0.01 h Mpc⁻¹) show independence on large scales, as is expected for the Gaussian data model, but they become correlated at k > 0.3 h Mpc⁻¹. The correlations arise from the redshift-space parametrisation that couples the amplitude of the power spectrum to β and σv. The upper right square in the figure represents the nearly 100% correlation between these two parameters.

## 6. Conclusions

Using VIPERS we have demonstrated a method of reconstructing the galaxy density field jointly with the redshift-space power spectrum, galaxy biasing function and galaxy luminosity function with minimal priors on these parameters. The Bayesian framework naturally accounts for the correlations between these observables. We adopt a likelihood function for the galaxy number counts that is given by a multivariate Gaussian and set a Gaussian prior on the density field. The solution that maximises the posterior distribution is given by the classical Wiener filter. To sample from the posterior distribution we add a Gaussian constrained realisation. Incorporating this density field estimator within a Gibbs sampler, we jointly sample the full posterior distribution including the power spectrum, bias and luminosity function parameters. We find encouraging results by using a multivariate Gaussian model for the likelihood and prior distributions, although more theoretically motivated descriptions may be used (Kitaura et al. 2012a; Jasche & Wandelt 2013a).

There are clear gains when jointly estimating correlated parameters. For instance the galaxy colour-density relation can be used to improve estimates of the density field. Furthermore it is well known that bias weighting galaxies when estimating the power spectrum leads to improved accuracy (Percival et al. 2004) and greater statistical power (Cai et al. 2011).

Moreover, the Bayesian framework provides a recipe for propagating uncertainties in the measurement, incorporating prior knowledge and constraints from the data, and it guarantees reliable error estimates. In VIPERS we account for inhomogeneous sampling and detailed angular masks. We correct for the selection function of VIPERS by up-weighting galaxies according to the magnitude-dependent spectroscopic sampling rate, while including the target sampling rate in the angular dependence of the survey selection function. These corrections are
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Galaxy bias measured from the full (red and blue combined) VIPERS galaxy sample in luminosity threshold bins. Reference data are taken from the VIPERS projected correlation function analysis (Marulli et al. 2013) and counts-in-cells probability distribution function analysis (Di Porto et al. 2014). We note that the redshift ranges differ.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Galaxy luminosity function inferred from the mean density Markov chain for red, blue and combined samples in redshift bins. Markers are plotted at the median value of the chain and the height of the rectangles indicates the 68% confidence interval. The Schechter function fits from Fritz et al. (2014) are overplotted for comparison.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Normalised correlation matrix of the parameters computed from the VIPERS Markov chain. The blocks represent the mean density, galaxy bias, power spectrum and RSD parameters. The structure in the covariance arises from the data model parameterisation. The values of luminosity and colour dependencies of galaxy bias and mean density within a redshift bin are strongly correlated, while they are only weakly correlated across redshift. On large scales the power spectrum covariance is diagonal, but at $k > 0.3 \hbar \text{Mpc}^{-1}$ the bins become correlated owing to coupling of the small-scale power with the redshift-space distortion parameters.</figcaption></figure>

fixed in our analysis, although for upcoming surveys it will be important to propagate the uncertainties in the selection function to the data products.

Investigating the covariances between parameters, we find strong correlations between galaxy bias and number density parameters within a given redshift bin. This is not unexpected since both these parameters depend on the one-point probability distribution function of the density field. On the other hand the correlation with the power spectrum is weak.

Our estimate of the power spectrum is effectively deconvolved from the survey window function (see Rota et al., in prep.) and we find that the covariance of the power spectrum bins is diagonal on large scales as expected from an unmasked Gaussian random field. On small scales, $k > 0.3 \hbar \text{Mpc}^{-1}$ we find significant correlations between power spectrum bins. On these scales correlations are expected due to the physical processes of structure formation; however, in this case the correlations

arise from the parameterisation of the data model. There is a degeneracy between the redshift-space distortion factors $\beta$ and $\sigma_v$ and the amplitude of the power spectrum on small scales. Nevertheless, the error estimate given by the Gibbs sampler closely matches the expectation of cosmic variance estimated from mock catalogues.

Our results are in good agreement with previous VIPERS measurements. We find values of the redshift-space distortion factor $\beta$ that are consistent with the correlation function analysis (de la Torre et al. 2013). Our values of luminosity dependent bias follow the trends expected from Marulli et al. (2013) and Di Porto et al. (2014) at $z < 0.9$. We further estimate the galaxy bias for colour samples finding a more pronounced dependence on luminosity for red galaxies than blue. The luminosity function we infer from the mean number density matches well with those found by Fritz et al. (2014) using the STY estimator.
In our analysis we have left the power spectrum, galaxy bias and number density without parameterisation. Despite this freedom, we find that the resulting errors in the key quantities such as the distortion parameter are only marginally larger than those given by traditional methods which can be strongly dependent on parameterisation.

Our methodology can be extended to jointly analyse multiple datasets in a self-consistent manner. A particular challenge when considering multiple surveys is dealing with the differences in angular coverage, sampling rates and galaxy types. The Bayesian approach provides a method of homogenising datasets allowing for consistent measurements.

For future studies with VIPERS we can consider the joint analysis with the VVDS-Wide spectroscopic survey (Garilli et al. 2008). Although the two surveys partially overlap, the selection function and sampling rates differ prohibiting their direct combination. However, through the Bayesian framework, the joint analysis becomes natural. We may further add constraints given by the density reconstructions in the gaps by the ZADE algorithm (Cucciati et al. 2014) or photometric redshift samples from the full CFHTLS Wide fields (Granett et al. 2012; Coupon et al. 2012). Sheer measurements in these fields can provide additional constraints on the underlying matter density providing a powerful probe in combination (Coupon et al. 2015). For upcoming surveys, this strategy will guarantee a complete and self-consistent picture of the Universe.

**Acknowledgements.** This work is based on observations collected at the European Southern Observatory, Cerro Paranal, Chile, using the Very Large Telescope under programs 182.A-0886 and partly 070.A-9007. Also based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work is based in part on data products produced at TERAPIX and the Canadian Astronomy Data Centre as part of the Canada-France-Hawaii Telescope Legacy Survey, a collaborative project of NRC and CNRS. The VIPERS web site is http://www.vipers.inaf.it/. We acknowledge the crucial contribution of the ESO staff for the management of service observations. In particular, we are deeply grateful to M. Hilker for his constant help and support of this program. Italian participation to VIPERS has been funded by INAF through PRIN 2008 and 2010 programs. L.G., and B.R.G. acknowledge support of the European Research Council through the Darklight ERC Advanced Research Grant (# 291521). A.P., K.M., and J.K. have been supported by the National Science Centre (grants UMO-2012/07/B/ST9/04425 and UMO-2013/09/D/ST9/04030), the Polish-Swiss Astro Project (co-financed by a grant from Switzerland, through the Swiss Contribution to the enlarged European Union), and the European Associated Laboratory Astrophysics Poland-France HECOLS. K.M. was supported by the Strategic Young Researcher Overseas Visits Program for Accelerating Brain Circulation (# R2405). O.L.F. acknowledges support of the European Research Council through the EARLY ERC Advanced Research Grant (# 268107). G.D.L. acknowledges financial support from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement # 202781. W.J.P. and R.T. acknowledge financial support from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement # 202686. W.J.P. is also grateful for support from the UK Science and Technology Facilities Council through the grant ST/1001204/1. E.B., F.M. and L.M. acknowledge the support from grants ASI-INAF I/023/12/0 and PRIN MIUR 2010-2011. L.M. also acknowledges financial support from PRIN INAF 2012. Y.M. acknowledges support from CNRS/INSU (Institut National des Sciences de l'Univers) and the Programme National Galaxies et Cosmologie (PNCG). C.M. is grateful for support from specific project funding of the Institut Universitaire de France and the LABEX OCEVU.

## References

* Alcock, C., & Paczynski, B. 1979, Nature, 281, 358

* Amendola, L., Appleby, S., Bacon, D., et al. 2013, Living Rev. Relativity, 16, 6

* Ata, M., Kitaura, F.-S., & Müller, V. 2015, MNRAS, 446, 4250

* Ballinger, W. E., Peacock, J. A., & Heavens, A. F. 1996, MNRAS, 282, 877

* Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15

* Blanton, M. R., Eisenstein, D., Hogg, D. W., Schlegel, D. J., & Brinkmann, J. 2005, ApJ, 629, 143

* Blanton, M. R., Hogg, D. W., Bahcall, N. A., et al. 2003, ApJ, 592, 819

* Bottini, D., Garilli, B., Maccagni, D., et al. 2005, PASP, 117, 996

* Brooks, S. P., & Gelman, A. 1998, J. Comp. Graph. Statist., 7, 434

* Cai, Y.-C., Bernstein, G., & Sheth, R. K. 2011, MNRAS, 412, 995

* Cappi, A., Marulli, F., Bel, J., et al. 2015, A&A, 579, A70

* Cole, S. 2011, MNRAS, 416, 739

* Coupon, J., Kilbinger, M., McCracken, H. J., et al. 2012, A&A, 542, A5

* Coupon, J., Arnouts, S., van Waerbeke, L., et al. 2015, MNRAS, 449, 1352

* Cucciati, O., Iovino, A., Marinoni, C., et al. 2006, A&A, 458, 39

* Cucciati, O., Granett, B. R., Branchini, E., et al. 2014, A&A, 565, A67

* Cui, W., Liu, L., Yang, X., et al. 2008, ApJ, 687, 738

* Guillardre, J.-C. J., Withington, K., Hudelot, P., et al. 2012, in SPIE Conf. Ser., 8448

* Davidzon, I., Bolzonella, M., Coupon, J., et al. 2013, A&A, 558, A23

* Davis, M., & Geller, M. J. 1976, ApJ, 208, 13

* de la Torre, S., & Peacock, J. A. 2013, MNRAS, 435, 743

* de la Torre, S., Guzzo, L., Peacock, J. A., et al. 2013, A&A, 557, A54

* Di Porto, C., Branchini, E., Bel, J., et al. 2014, ArXiv e-prints [arXiv:1406.6692]

* Dressler, A. 1980, ApJ, 236, 351

* Efstathiou, G., & Moody, S. J. 2001, MNRAS, 325, 1603

* Efstathiou, G., Ellis, R. S., & Peterson, B. A. 1988, MNRAS, 232, 431

* Enßlin, T. A., Frommert, M., & Kitaura, F. S. 2009, Phys. Rev. D, 80, 105005

* Erdoğan, P., Lahav, O., Zaroubi, S., et al. 2004, MNRAS, 352, 939

* Eriksen, H. K., O'Dwyer, I. J., Jewell, J. B., et al. 2004, ApJS, 155, 227

* Fritz, A., Scodeggio, M., Ilbert, O., et al. 2014, A&A, 563, A92

* Garilli, B., Le Fèvre, O., Guzzo, L., et al. 2008, A&A, 486, 683

* Garilli, B., Guzzo, L., Scodeggio, M., et al. 2014, A&A, 562, A23

* Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457

* Granett, B. R., Guzzo, L., Coupon, J., et al. 2012, MNRAS, 421, 251

* Guzzo, L., Scodeggio, M., Garilli, B., et al. 2014, A&A, 566, A108

* Hamilton, A. J. S. 1998, in The evolving universe, ed. D. Hamilton, Astrophys. Space Sci. Libr., 231, 185

* Hess, S., & Kitaura, F.-S. 2014, ArXiv e-prints [arXiv:1412.7310]

* Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles (Bristol, PA, USA: Taylor & Francis, Inc.)

* Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863

* Jasche, J., & Kitaura, F. S. 2010, MNRAS, 407, 29

* Jasche, J., & Wandelt, B. D. 2012, MNRAS, 425, 1042

* Jasche, J., & Wandelt, B. D. 2013a, MNRAS, 432, 894

* Jasche, J., & Wandelt, B. D. 2013b, ApJ, 779, 15

* Jasche, J., Kitaura, F. S., & Enßlin, T. A. 2009, ArXiv e-prints [arXiv:0901.3043]

* Jasche, J., Kitaura, F. S., Li, C., & Enßlin, T. A. 2010a, MNRAS, 409, 355

* Jasche, J., Kitaura, F. S., Wandelt, B. D., & Enßlin, T. A. 2010b, MNRAS, 406, 60

* Jasche, J., Leclercq, F., & Wandelt, B. D. 2015, J. Cosmology Astropart. Phys., 1, 36

* Jewell, J., Levin, S., & Anderson, C. H. 2004, ApJ, 609, 1

* Jewell, J. B., Eriksen, H. K., Wandelt, B. D., et al. 2009, ApJ, 697, 258

* Jing, Y. P. 2005, ApJ, 620, 559

* Kaiser, N. 1984, ApJ, 284, L9

* Kaiser, N. 1987, MNRAS, 227, 1

* Kitaura, F.-S. 2013, MNRAS, 429, L84

* Kitaura, F. S., & Enßlin, T. A. 2008, MNRAS, 389, 497

* Kitaura, F. S., Jasche, J., Li, C., et al. 2009, MNRAS, 400, 183

* Kitaura, F.-S., Jasche, J., & Metcalf, R. B. 2010, MNRAS, 403, 589

* Kitaura, F.-S., Erdoğan, P., Nuza, S. E., et al. 2012a, MNRAS, 427, L35

* Kitaura, F.-S., Gallerani, S., & Ferrara, A. 2012b, MNRAS, 420, 61

* Lahav, O., & Suto, Y. 2004, Living Rev. Relativity, 7, 8

* Lahav, O., Fisher, K. B., Hoffman, Y., Scharf, C. A., & Zaroubi, S. 1994, ApJ, 423, L93

* Landy, S. D., & Szalay, A. S. 2002, ApJ, 579, 76

* Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]

* Le Fèvre, O., Saisse, M., Mancini, D., et al. 2003, in SPIE Conf. Ser. 4841, eds. M. Iye, & A. F. M. Moorwood, 1670

* Lesgourgues, J. 2011, ArXiv e-prints [arXiv:1104.2932]

* Levi, M., Bebek, C., Beers, T., et al. 2013, ArXiv e-prints [arXiv:1308.0847]

* LSST Dark Energy Science Collaboration. 2012, ArXiv e-prints [arXiv:1211.0310]

* Marin, J.-M., & Robert, C. P. 2007, Bayesian Core: A Practical Approach to Computational Bayesian Statistics (Springer Texts in Statistics), 1st edn. (Springer)
Marulli, F., Bolzonella, M., Branchini, E., et al. 2013, A&A, 557, A17

Micheletti, D., Iovino, A., Hawken, A. J., et al. 2014, A&A, 570, A106

Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347

Nusser, A., & Davis, M. 1994, ApJ, 421, L1

Peebles, P. J. E. 1980, The large-scale structure of the universe (Princeton University Press), 435

Percival, W. J., Verde, L., & Peacock, J. A. 2004, MNRAS, 347, 645

Planck Collaboration I. 2015, ArXiv e-prints [arXiv:1502.01582]

Pollo, A., Meneux, B., Guzzo, L., et al. 2005, A&A, 439, 887

Pope, A. C., Matsubara, T., Szalay, A. S., et al. 2004, ApJ, 607, 655

Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018

Rybicki, G. B., & Press, W. H. 1992, ApJ, 398, 169

Sandage, A., Tammann, G. A., & Yahil, A. 1979, ApJ, 232, 352

Schaap, W. E., & van de Weygaert, R. 2000, A&A, 363, L29

Schechter, P. 1976, ApJ, 203, 297

Smith, R. E. 2012, MNRAS, 426, 531

Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Tegmark, M., Hamilton, A. J. S., & Xu, Y. 2002, MNRAS, 335, 887

The Dark Energy Survey Collaboration. 2005, ArXiv e-prints
[arXiv:astro-ph/0510346]

Trotta, R. 2008, Contemporary Phys., 49, 71

Wandelt, B. D., Larson, D. L., & Lakshminarayanan, A. 2004, Phys. Rev. D, 70, 083511

Wang, H., Mo, H. J., Yang, X., Jing, Y. P., & Lin, W. P. 2014, ApJ, 794, 94

Zaroubi, S., Hoffman, Y., Fisher, K. B., & Lahav, O. 1995, ApJ, 449, 446

Zucca, E., Bardelli, S., Bolzonella, M., et al. 2009, A&A, 508, 1217

7 Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France

8 Canada-France-Hawaii Telescope, 65-1238 Mamalahoa Highway, Kamuela, HI 96743, USA

9 INAF-Osservatorio Astronomico di Bologna, via Ranzani 1, 40127 Bologna, Italy

10 INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via Bassini 15, 20133 Milano, Italy

11 Laboratoire Lagrange, UMR 7293, Université de Nice Sophia Antipolis, CNRS, Observatoire de la Côte d'Azur, 06300 Nice, France

12 Astronomical Observatory of the University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

13 Dipartimento di Fisica e Astronomia-Alma Mater Studiorum Università di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy

14 INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, 34143 Trieste, Italy

15 Institute of Physics, Jan Kochanowski University, ul. Swietokrzyska 15, 25-406 Kielce, Poland

16 National Centre for Nuclear Research, ul. Hoza 69, 00-681 Warszawa, Poland

17 INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy

18 Institute d'Astrophysique de Paris, UMR 7095 CNRS, Université Pierre et Marie Curie, 98 bis boulevard Arago, 75014 Paris, France

19 Astronomical Observatory of the Jagiellonian University, Orla 171, 30-001 Cracow, Poland

20 School of Physics and Astronomy, University of St. Andrews, North Haugh, St Andrews KY16 9SS, UK

21 INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica Bologna, via Gobetti 101, 40129 Bologna, Italy

22 INAF-Istituto di Radioastronomia, via Gobetti 101, 40129 Bologna, Italy

23 Institute of Cosmology and Gravitation, Dennis Sciama Building, University of Portsmouth, Burnaby Road, Portsmouth, PO1 3FX, UK

24 Centre de Physique Théorique, UMR 6207 CNRS-Université de Provence, Case 907, 13288 Marseille, France

25 SUPA, Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

<footnote n="1">1 INAF-Osservatorio Astronomico di Brera, via Brera 28, 20122 Milano, via E. Bianchi 46, 23807 Merate, Italy</footnote>

<footnote n="2">2 Dipartimento di Matematica e Fisica, Università degli Studi Roma Tre, via della Vasca Navale 84, 00146 Roma, Italy</footnote>

<footnote n="3">3 INFN, Sezione di Roma Tre, via della Vasca Navale 84, 00146 Roma, Italy</footnote>

<footnote n="4">4 INAF-Osservatorio Astronomico di Roma, via Frascati 33, 00040 Monte Porzio Catone (RM), Italy</footnote>

<footnote n="5">5 Dipartimento di Fisica, Università di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy</footnote>

<footnote n="6">6 INAF-Osservatorio Astronomico di Torino, 10025 Pino Torinese, Italy</footnote>
Appendix A: Anti-aliasing filter

Binning the continuous galaxy density field onto a grid may be characterised as a smoothing operation followed by discrete sampling (Hockney & Eastwood 1988). In the simplest approach, if particles are assigned to the nearest grid point (NGP) the smoothing kernel W(x) has a top-hat shape with the dimension of the grid cell. More extended kernels can be chosen that distribute the weight of the particle over multiple cells. Common assignment schemes include cloud-in-cell (CIC) and triangle-shaped-cell (TSC) which correspond to iterative smoothing operations with the same top-hat filter. Smoothing the field damps small scale power, but better localises the signal in k-space which is beneficial in Fourier analyses.

A consequence of using the fast Fourier transform (FFT) algorithm is that the signal must be discretised onto a finite number of wave modes or equivalently that periodic boundary conditions are imposed. The observed power is thus a sum of all harmonics of the fundamental wavelength which are known as aliases (Jing 2005; Hockney & Eastwood 1988):

$$P'(k) = \sum_{n} |W(k + 2k_{N}n)|^{2} P(k + 2k_{N}n). \qquad (A.1)$$

As long as the signal is sufficiently well sampled and band-limited, meaning that there is no power above the Nyquist frequency, the signal may be recovered without error. However, most signals of interest are not band limited, and so after sampling onto the grid, the harmonics overlap and the signal cannot be exactly recovered.

The ideal anti-aliasing filter in k-space is the top-hat that cuts all power above the Nyquist frequency. However, the signal cannot be localised both in k-space and in position-space and a sharp cut in k-space will distribute the mass of a particle over the entire grid.

A practical implementation of the ideal filter was developed by Jasche et al. (2009). The method involves first super-sampling the field by assigning particles to a grid with resolution higher than the target grid. Transforming to Fourier space via FFT, the high resolution grid is filtered setting to 0 all modes above the target Nyquist frequency. Finally the grid is transformed back to position-space and down sampled to the target grid. The technique is very efficient with the additional trick that the down-sampling step can be carried out by the inverse FFT by cutting and reshaping the Fourier grid.

The algorithm of Jasche et al. (2009) is limited only by the memory needed to apply the FFT on the high resolution grid, so it is competitive with common position-space cell assignment schemes. However, the sharp cut in k-space leads to an extended, oscillatory distribution in position-space. Such a decentralised and unphysical cell assignment is undesirable for estimation of the density field.

An alternative approach was taken by Cui et al. (2008) who introduce the use of Daubechies wavelets to approximate the ideal filter while keeping compactness in position-space. We do not consider this technique here because the centre of mass of the wavelets adopted is offset from the particle position and so they are not well suited for estimation of the density field.

As a compromise, we test a smooth cutoff with a k-space filter with the form:

$$W(k_x, k_y, k_z) = \prod_{i=(x,y,z)} \frac{1}{1 + |k_i/k_N|^\alpha} \qquad (A.2)$$

The parameter $\alpha$ adjusts the sharpness of the cut: we test $\alpha = 100$ (Sup hard) and $\alpha = 8$ (Sup soft).

<figure><img src="image_13.png" /><figcaption>Fig. A.1. Top panel: an illustration of the mass-assignment kernels in two-dimensions: a hard k-space cut (Sup hard), a soft k-space cut (Sup soft), nearest-grid-point (NGP) and cloud-in-cell (CIC). The lower three panels give a comparison of the kernel functions (in one dimension) for various mass assignment schemes. Top: the kernels in Fourier space. The Nyquist frequency is 1 pixel$^{-1}$ on the scale. Middle: the kernels in position space. Bottom: the cumulative power of each kernel in position space.</figcaption></figure>

Figure A.1 shows the comparison of the mass assignment schemes including the super-sampling technique. We see that the hard cut (Sup hard) leads to ringing in position-space (the transform is the sinc function). The amplitude is significantly reduced when the parameter $\alpha = 8$ (Sup soft) and the first peak is more compact than CIC. In k-space, the filter remains nearly 1 to ~0.7kN while NGP, CIC and TSC functions are dropping. Beyond the Nyquist frequency, the soft cut-off kernel drops faster than TSC without oscillatory behaviour. The bottom panel of A.1 shows the cumulative power $S(R) = \int_{0}^{R} |W(r)|^{2}dr$. We see that soft-cut off scheme has compactness similar to CIC and NGP, while the latter schemes severely damp the total power.

Figure A.2 compares the aliased power in the power spectrum measurement. The top four panels show the shot noise power multiplied by the kernels. The principal contribution ($n = 0$) and the first harmonic ($n = -1$) are plotted along with the full sum of all harmonics. We see for instance, that the power measured with the NGP scheme mixes power from all harmonics. The hard cut-off filter performs nearly ideally with no aliasing effect. The soft cut-off also performs well matching the aliasing characteristics of the TSC scheme.
<figure><img src="image_14.png" /><figcaption>Fig. A.2. A comparison of the shot noise power as a function of wavenumber under various mass assignment schemes. The Nyquist frequency is 1 pixel⁻¹ on the scale. The top four plots show the principal component (thick line), the *n* = −1 harmonic and the sum of all harmonics (thin line). The bottom frame shows the fraction of aliased power for each assignment scheme.</figcaption></figure>

To implement the soft cut-off filter we use the following practical super-sampling recipe:

1. Assign particles to a grid with resolution increased by a factor *f*. The assignment is done with the CIC scheme. We set *f* = 8.

2. Transform the field to Fourier space with the FFT.

3. Multiply Fourier modes by the *k*-space filter.

4. Transform back to position space with the inverse-FFT.

5. Take a sample of the grid every *f* cells to down-sample to the target resolution.

## Appendix B: Gibbs sampler

Here we give the algorithms used to sample from the conditional probability functions.

### B.1. Sampling the density field

Using a galaxy survey we count galaxies in a given sample *l* and construct the spatial field *N*<sub>*l*</sub> which is a vector of length *n*<sub>cells</sub>. We will write the set of *m* galaxy samples (e.g. different luminosity and redshift bins) as {*N*<sub>*l*</sub>} ≡ {*N*<sup>l</sup><sub>1</sub>, *N*<sup>l</sup><sub>2</sub>, ..., *N*<sup>l<sub>m</sub></sup>}. Each sample *l* has a corresponding mean density *N̄*<sub>*l*</sub> and bias *b*<sub>*l*</sub>.

We write the conditional probability for the underlying density field *δ* as

$$p(\delta|\{N_l\}, \{\bar{N}_l\}, \{b_l\}, \mathbf{S}) \propto p(\{N_l\}|\delta, \{\bar{N}_l\}, \{b_l\}, \mathbf{S})p(\delta|\mathbf{S}) \\ = \prod_{l=1}^{N_{\text{samples}}} p(N_l|\delta, \bar{N}_l, b_l, \mathbf{S})p(\delta|\mathbf{S}). \quad (\text{B.1})$$

To write the last line we use the property that the number counts of different samples *l* are conditionally independent but depend on the common underlying density field.

For clarity we now explicitly index the cells with subscript *i*. We adopt a Gaussian model for the number counts and write the log of the likelihood $p \propto \chi^2$ as

$$-2 \log p(N_l|\bar{N}_l, b_l, \delta) = \sum_{i=1}^{n_{\text{cells}}} \frac{(N_{l,i} - \bar{N}_l w_{l,i} (1 + b_l D(z_i) \delta_i))^2}{\sigma_{l,i}^2} + \log(2\pi \sigma_{l,i}^2). \quad (\text{B.2})$$

The variance of counts is generalised as $\sigma_{l,i}^2$ and may be different from the Poisson expectation $\bar{N}w$ due to the mask and aliasing filter. However, we neglect noise correlations between cells which would introduce off-diagonal terms. The sums are carried out over cells with variance $\sigma_{l,i}^2 > 0$.

Next we consider the density field prior. We set a Gaussian model giving

$$-2 \log p(\delta|\mathbf{S}) = \sum_{j=1}^{n_{\text{cells}}} \sum_{i=1}^{n_{\text{cells}}} (\mathbf{S}^{-1})_{ij} \delta_i \delta_j + \log(2\pi \det \mathbf{S}). \quad (\text{B.3})$$

The correlation function of *δ* is given by the anisotropic covariance matrix *S*. In Fourier space the matrix is diagonal and given by Eq. (10). When writing the posterior we now drop the terms that do not depend on *δ* giving

$$-2 \log p(\delta|\{N_l\}, \{\bar{N}_l\}, \{b_l\}, \mathbf{S}) = \sum_{l=1}^{m_{\text{samples}}} \sum_{i=1}^{n_{\text{cells}}} \frac{(N_{l,i} - \bar{N}_l w_{l,i} (1 + b_l D_i \delta_i))^2}{\sigma_{l,i}^2} + \sum_{j=1}^{n_{\text{cells}}} (\mathbf{S}^{-1})_{ij} \delta_i \delta_j. \quad (\text{B.4})$$

Differentiating the log posterior with respect to *δ* we find the equation for the maximum a posteriori estimator which is also called the Wiener filter. The estimate $\hat{\delta}$ is given by

$$\sum_{l=1}^{m_{\text{samples}}} \frac{(\bar{N}_l w_{l,i} b_l D_i)^2}{\sigma_{l,i}^2} \hat{\delta}_i + \sum_{j=1}^{n_{\text{cells}}} (\mathbf{S}^{-1})_{i,j} \hat{\delta}_j = \sum_{l=1}^{m_{\text{samples}}} \frac{\bar{N}_l w_{l,i} b_l D_i}{\sigma_{l,i}^2} (N_{l,i} - \bar{N}_l w_{l,i}). \quad (\text{B.5})$$

It is informative to point out how the Wiener filter operation combines the galaxy subsamples. The right hand side of Eq. (B.5) shows the weighted combination, given by

$$\sum_{l=1}^{m_{\text{samples}}} \frac{\bar{N}_l w_{l,i} b_l D_i}{\sigma_{l,i}^2} (N_{l,i} - \bar{N}_l w_{l,i}) \quad (\text{B.6})$$

where we consider cell *i* and the sum is over galaxy samples indexed by *l*. The subsamples are being weighted by their relative biases *b*<sub>*l*</sub>. This form of weighting for the density field was derived by Cai et al. (2011) and is optimal in the case of Poisson sampling but also matches the weights for the power spectrum (Percival et al. 2004).

To generate a residual field $\delta_r$ with the correct covariance we follow the method of Jewell et al. (2004). We draw two sets of Gaussian distributed random variables with zero mean and unit
variance: $w_{1,i}$ and $w_{2,i}$. The residual field is then found by solving the following linear equations.

$$ \zeta_i = \sum_{l=1}^{m_{\text{samples}}} \frac{\bar{N}_l w_{l,i} b_l D_i}{\sigma_{l,i}} w_{1,i} \quad (B.7) $$

$$ \eta_i = \sum_{j=1}^{n_{\text{cells}}} (\sqrt{S^{-1}})_{i,j} w_{2,j} \quad (B.8) $$

$$ \sum_{l=1}^{m_{\text{samples}}} \frac{(\bar{N}_l w_{l,i} b_l D_i)^2}{\sigma_{l,i}^2} \delta_{r,i} + \sum_{j=1}^{n_{\text{cells}}} (S^{-1})_{i,j} \delta_{r,j} = \zeta_i + \eta_i. \quad (B.9) $$

The final constrained realisation is given by the sum $\delta_{CR} = \hat{\delta} + \delta_r$.

The terms including the signal covariance matrix are computed in Fourier space where the matrix is diagonal. The transformation is done using the Fast fourier transform algorithm. We solve Eqs. (B.5) and (B.9) using the linear conjugate gradient solver bicgstab provided in the SciPy python library<fnref n="1" />.

## B.2. Sampling the signal

The posterior distribution for the signal is

$$ p(S|\delta) \propto p(\delta|S)p(S). \quad (B.10) $$

We take a flat prior leaving $p(S|\delta) \propto p(\delta|S)$ which is given by Eq. (B.3). This is recognised as an inverse-Gamma distribution for $S$ which can be sampled directly (Kitaura & Enßlin 2008; Jasche et al. 2010b).

We parametrisate the signal in Fourier space as the product of the real-space power spectrum and redshift-space distortion model (Eq. (10)). We sample the parameters in two steps. First we fix $\beta$ and $\sigma_v$ and draw a real-space power spectrum from the inverse-gamma distribution.

The normalisation of the power spectrum is fixed by

$$ \sigma_R^2 = \int d^3 k P(k) |W(k)|^2. \quad (B.11) $$

In the regime where shot noise is more important than cosmic variance, this method of sampling is inefficient. An alternative sampling scheme was proposed to improve the convergence (Jewell et al. 2009; Jasche et al. 2010b). The new power spectrum is drawn making a step in both $P$ and $\delta$ such that $\delta^{s+1} = \sqrt{P^{s+1}/P^s \delta^s}$. The consequence is that the conditional posterior $p(P|\delta)$ is unchanged but the data likelihood is modified through $\delta$. We perform a sequence of Metropolis-Hastings steps to jointly draw $P$ and $\delta$. We alternate between the two modes for sampling $P$ on each step of the Gibbs sampler.

The redshift-space distortion parameters $\beta$ and $\sigma_v$ are next sampled jointly. These parameters are strongly degenerate on the scales we consider. To efficiently sample these we compute the eigen decomposition of the Hessian matrix to rotate into a coordinate system with two orthogonal parameters. We evaluate the joint probability over a two-dimensional grid in the orthogonal space. Finally we use a rejection sampling algorithm to jointly draw values of $\beta$ and $\sigma_v$.

## B.3. Sampling galaxy bias

The bias is sampled in the manner described by Jasche & Wandelt (2013b). The conditional probability distribution for the

bias of a given galaxy sample is a Gaussian with mean given by:

$$ \mu_b = \frac{\sum \bar{N}_i w_i D_i \delta_i (N_i - \bar{N}_i) / \sigma_i^2}{\sum_i (\bar{N}_i w_i D_i \delta_i)^2 / \sigma_i^2} \quad (B.12) $$

and variance

$$ \sigma_b^2 = \frac{1}{\sum_i (\bar{N}_i w_i D_i \delta_i)^2 / \sigma_i^2} \quad (B.13) $$

After drawing a value from a Gaussian distribution with the given mean and variance, it is limited to the range $0.5 < b < 4.0$.

## B.4. Sampling mean density

We take a Poisson model for the mean density. Taking a uniform prior we have

$$ p(\bar{N}|\delta, N, b, S) \propto p(N|\delta, b, \bar{N}, S) p(\bar{N}) \quad (B.14) $$

$$ \propto \prod_i \frac{[\bar{N}w_i (1 + bD\delta)]^{N_i}}{N_i!} e^{-\bar{N}w_i(1+bD\delta_i)}. \quad (B.15) $$

Keeping only terms that depend on $\bar{N}$ we have

$$ \log p(\bar{N}) = \sum_i N_i \log \bar{N} - \bar{N}w_i (1 + bD\delta_i). \quad (B.16) $$

We draw a sample from this distribution by computing the likelihood over a grid of values and using a rejection sampling algorithm.

# Appendix C: Convergence analysis

The Gelman-Rubin convergence diagnostic (Gelman & Rubin 1992; Brooks & Gelman 1998) involves running multiple independent Markov chains and comparing the single-chain variance to the between-chain variance. The scale reduction factor $R$ is the ratio of the two variance estimates. In Fig. C.1 we show this convergence diagnostic for three parameters: the power spectrum at $k = 0.2 h \text{Mpc}^{-1}$, the distortion parameter $\beta$ and the velocity dispersion $\sigma_v$. These parameters are the slowest to converge. The $R$ factor is computed up to a given maximum step $i$ in the chain after discarding the first $i/2$ steps. We consider the chain to be properly “burned-in” when $R < 1.2$. Based on this diagnostic we set the maximum step in our analysis to be 2000 and discard the first 1000 steps.

<figure><img src="image_15.png" /><figcaption>Fig. C.1. Gelman-Rubin convergence diagnostic $R$ computed from a single mock catalogue with 10 chains.</figcaption></figure>

<footnote n="1">¹ http://www.scipy.org/</footnote>