© ESO, 2015

1. Introduction

Since the very first galaxy surveys have revealed the complex structure of the Universe, the mapping of its evolution has become an important part of cosmology. Large surveys revealed that the Universe is composed of dense regions like clusters and filaments, but also of almost empty voids. Additionally, the underlying dark matter structure and its evolution is shown to follow the visible baryonic matter, though the luminous matter is biased in relation to the dark matter distribution. It is still unclear what exactly this baryonic – dark matter relation looks like, especially in the early stages of galaxy formation, and how it changed through time.

The correlation function is a commonly used tool to describe how galaxies are clustered as a function of scale, and allows constraints to be put on the evolution of galaxies (Kauffmann et al. 1999; Zehavi et al. 2011). It is based on the very simple idea of measuring the probability of finding two galaxies at a given redshift and separation (Peebles 1980). The galaxy correlation function can be interpreted by introducing two different approximations. The first and most extensively used formalism is based on a simple power-law approximation of the correlation function of the form ξ(r) = (r/r₀)^{− γ} with two free parameters: the correlation length r₀, which describes how strongly galaxies are clustered, and the slope γ (Davis & Peebles 1983). The second, more recent and detailed approximation is based on halo occupation models (Seljak 2000; Peacock & Smith 2000; Magliocchetti & Porciani 2003; Zehavi et al. 2004; Zheng et al. 2005). In this framework the correlation function is built from two components that have their influence on different scales. The one-halo term dominates on small scales (≲ 1.5 h^-1 Mpc) and describes clustering of galaxies that reside within dark matter haloes. The two-halo term describes large-scale (≳ 3 h^-1 Mpc) galaxy clustering between different haloes.

Both formulations have been extensively used in the past years. Correlation function measurements have been produced for most large galaxy surveys, such as the Sloan Digital Sky Survey (SDSS, Connolly et al. 2002; Zehavi et al. 2004), 2dF Galaxy Redshift Survey (2dFGRS, Magliocchetti & Porciani 2003), VIMOS-VLT Sky Survey (VVDS, Le Fèvre et al. 2005; Pollo et al. 2006), VIMOS Public Extragalactic Survey (VIPERS, Marulli et al. 2013), and DEEP2 (Coil et al. 2006).

Thanks to this remarkable effort it is now very well established that the strength of the galaxy clustering for the general population of galaxies only mildly evolves from intermediate redshifts z ~ 1 to z ~ 0. It has been found that galaxy clustering depends on a variety of galaxy properties like the luminosity, morphology, colour, and spectral type of galaxies: luminous galaxies tend to be more clustered than faint ones and red galaxies with old stellar population are found to be more clustered than young blue ones (Norberg et al. 2002; Zehavi et al. 2004, 2011, 2012; Pollo et al. 2006; de la Torre et al. 2007; Coil et al. 2008; Quadri et al. 2008; Meneux et al. 2006, 2008, 2009; Hartley et al. 2008, 2010, 2013, 2015; Skibba et al. 2009; Abbas et al. 2010; Tinker & Wetzel 2010; Wake et al. 2011; Coupon et al. 2012).

At redshifts higher than z ~ 2 the situation is less clear. Various difficulties, mainly in collecting statistically significant and representative samples, need to be overcome. While a number of attempts to measure galaxy clustering have been made at early epochs, interpreting results is not straightforward, mainly because of uncertainties connected to redshift determination, small volumes covered, and different galaxy populations selected from a range of methods. Most of the measurements at high redshifts (z> 1) are produced using photometric surveys targeting specific classes of galaxies or applying specific observation techniques, like the BzK method (Kong et al. 2006; Lin et al. 2012) or Lyman-break galaxy (LBG) selection (Foucaud et al. 2003; Ouchi et al. 2004; Adelberger et al. 2005; Kashikawa et al. 2006; Savoy et al. 2011; Bielby et al. 2013). In general, within a given sample and based on angular correlation function measurements, some evidence has been presented at z> 1 showing that clustering also seems to depend on luminosity, stellar mass, or colour as observed at lower redshift, with bright galaxies clustering more strongly than faint ones (Savoy et al. 2011) and passive galaxies clustering more strongly than star-forming galaxies at a given stellar mass (Lin et al. 2012). However the relation of galaxy samples used in these analyses to the general population of galaxies is not well established, which makes the study of evolution of the galaxy clustering at z> 2 difficult. Moreover these galaxies cannot be easily connected with the galaxy populations at lower redshifts, which makes it difficult to conduct a consistent study of galaxy clustering evolution from high redshift to z = 0.

From the theoretical point of view, the correlation function should evolve with time – ξ = ξ(r,t) because the density field of the Universe evolves over time. In the framework of the Newtonian linear perturbation theory (and in the matter-dominated era of the history of the Universe) in the co-moving coordinates the density contrast can be decomposed into time-dependent and space-dependent factors: δ(r,t) = D(t)δ₀(r), where δ₀(r) is the present day value of the density contrast at a given location, and D(t) is referred to as the growth factor, which depends on the parameters of the assumed cosmological models. In practice, it means that the spatial shape of the density fluctuations in co-moving coordinates does not change, and only their amplitude increases. Consequently, the amplitude of the correlation function should increase over time: ξ(r,t) = D²(t)ξ(r,t₀), with a dependence on the cosmological parameters (Peebles 1980; Schneider 2006).

The consistent measurement of the evolution of a correlation function amplitude for a given galaxy population is therefore enough – in principle – to test the paradigm of the gravitational perturbations as the origin of the large-scale structure of the Universe (and a way to estimate cosmological parameters). Extending these measurements to as high redshifts as possible is necessary both for the theoretical framework of the large-scale structure evolution and for the galaxy formation and evolution models.

Fig. 1 Distribution of galaxies with spectroscopic redshifts 2 <z< 5 in two independent VUDS fields: COSMOS on the left panel and VVDS-02h on right panel. The blue squares indicates VIMOS pointing centres.

In this work we attempt to improve our current view of the evolution of the clustering of galaxies at z> 2, investigating how galaxy clustering evolved from the early phases of galaxy assembly to present times. We present a clustering analysis of the VIMOS Ultra Deep Survey (VUDS), the largest spectroscopic survey covering the redshift range 2 <z< 6 in a continuous way to date (Le Fèvre et al. 2014), using a sample of ~ 3000 galaxies with confirmed spectroscopic redshifts in the range 2 <z< 5. The sample covers a total area of 0.81deg² observed in two independent fields, COSMOS and VVDS-02h, with a mean redshift value of z ~ 3. We use a halo occupation distribution (HOD) model to put constrains on the properties of dark matter haloes hosting star-forming galaxies at redshifts z ~ 3.

The paper is organized as follows. In Sect. 2 we briefly describe the properties of the VUDS survey and our selected samples. The method used to measure the correlation function and derive power-law and HOD fits is presented in Sect. 3. Results are described in Sect. 4. We discuss our findings and compare them to other work, both at low and high redshifts, in Sect. 5, and we summarize our work in Sect. 6.

We adopt a flat ΛCDM cosmological model with Ω_m = 0.25, Ω_Λ = 0.75 (we note that using the latest cosmology parameters from e.g. Planck Collaboration XVI 2014 does not change conclusions drawn from clustering analysis, see Zehavi et al. 2011). The Hubble constant is normally parametrized via h = H₀/ 100 to ease comparison with previous works, while a value H₀ = 70 km s^-1 Mpc is used when computing absolute magnitudes and stellar masses. We report correlation length measurements in comoving coordinates and express magnitudes in the AB system.

2. Data

2.1. VUDS Survey summary

The VIMOS Ultra Deep Survey (VUDS, Le Fèvre et al. 2015) is a spectroscopic survey targetting ~ 10 000 galaxies performed with the VIMOS multi-object spectograph (Le Fèvre et al. 2003) at the European Southern Observatory Very Large Telescope. The main aim of the survey is to study early phases of galaxy formation and evolution at 2 <z< 6. The survey covers a total area of 1deg² in three independent fields, reducing the effect of cosmic variance, an important property for galaxy clustering measurements. The majority (~ 88%) of targets are selected based on photometric redshifts (z_phot> 2.4 ± 1σ) derived from deep multi-band photometry available in the VUDS fields. Some additional targets that not satisfy this criterion are selected by colour (mainly Lyman Break Galaxies, LBGs). Finally, in order to fill the remaining available space on the slit mask, targets which do not meet any of the two previous conditions but with the imposed limit of i_AB = 25 were selected. Details about the survey strategy, target selection, as well as data processing and redshift measurements are presented in Le Fèvre et al. (2015). Here we briefly describe VUDS features which are relevant for our work.

Redshift measurements were carried out in a way similar to that developed for the VVDS survey (Le Fèvre et al. 2005, 2013), the zCOSMOS survey (Lilly et al. 2007), and VIPERS survey (Guzzo et al. 2014). The core engine for redshift measurement is cross-correlation of the observed spectrum with reference templates using the EZ redshift measurement code (Garilli et al. 2010). At the end of the process, each redshift measurement is assigned a reliability flag which expresses the reliability of the measurement:

• Flag 0: No redshift could be assigned;

• Flag 1: 50−75% probability of being correct;

• Flag 2: 75−85% probability of being correct;

• Flag 3: 95−100% probability of being correct;

• Flag 4: 100% probability of being correct;

• Flag 9: ~ 80% probability of being correct; spectrum has a single emission line

The VUDS sample benefits from an extended multi-wavelength data set (see Le Fèvre et al. 2015). The multi-wavelength photometry is used to compute absolute magnitudes and stellar masses from SED fitting using the Le Phare code (Arnouts et al. 1999; Ilbert et al. 2006), as described in detail by Ilbert et al. (2013) and references therein.

2.2. The sample selection for clustering analysis

In this work we only use the objects with assigned reliability flag zflag = 2,3,4,9 in the redshift range 2 <z< 5. We use two independent fields, COSMOS and VVDS-02h, covering a total area 0.81deg², as data for the third VUDS field in the ECDFS are still too sparse to be used for clustering analysis. Covered area corresponds to a volume ~ 3 × 10⁷ Mpc³ sampled in the redshift range 2 <z< 5. Our sample consists of 3022 galaxies with reliable spectroscopic redshifts in the range 2 <z< 5. The sample is summarized in Table 1.

The spatial distribution of the spectroscopic galaxy sample in each of these fields is presented in Fig. 1, while Fig. 2 presents their redshift distributions. Mean redshift value for the whole sample is z ~ 2.96. In the following analysis we use three (not independent) redshift bins. The first sample includes the whole star-forming galaxy population observed by VUDS in the wide range 2 <z< 5 with a mean redshift z ~ 2.96. Then, a sample is split into two subsamples: a low redshift sample with a mean redshift in the range z = [2, 2.9], and a high redshift sample with in z = [2.9, 5]. These two bins cover approximately similar cosmic time periods. The general properties of these samples including median luminosity ${\mathit{M}}_{\mathit{U}}^{\mathrm{median}}$ and stellar mass $\log {\mathit{M}}_{\mathrm{\ast }}^{\mathrm{median}}$ are listed in Table 2.

Because of the i_AB = 25 limit of the VUDS selection, the galaxies in the high redshift sample are ~ 0.5 mag brighter and ~ 0.2 M_⊙ more massive than those in the low redshift subsample, which is the expected effect taking into account the VUDS selection strategy. In the analysis presented in this paper we do not operate on volume-limited subsamples. However our data is approximately volume-complete above an absolute luminosity corresponding to galaxies above the characteristic luminosity M_⋆ of galaxies at these redshifts with M_NUV< −20 for 2 <z< 3 and M_NUV< −20.5 to −21 over 3 <z< 5 (Fig. 18 of Le Fèvre et al. 2015). The impact of the sample properties on the clustering measurements of our three samples is fully discussed in Sect. 5.

Fig. 2 Redshift distribution of the used VUDS galaxy sample for redshift range 2.0 <z< 5.0. The filled grey histogram represents total sample of reliable galaxies (with zflag = 2−4,9, ~ 80% reliability), while the red and blue histograms represent the contribution from COSMOS and VVDS2h fields respectively.

2.3. VUDS mock catalogues

To estimate uncertainties of the correlation function (see Sect. 3.2) and to test robustness (see Appendix A) of our clustering measurements we make use of a large number of mock galaxy samples, which are designed to mirror the VUDS sample in the range 2 <z< 5 in the most realistic way.

The 66 independent mock samples have been created based on a method which includes the HOD and the stellar-to-halo mass relation (SHMR). This method is the same as developed and introduced by de la Torre et al. (2013), and we refer the reader to the details presented in this paper. We provide below a brief summary of its most important features.

In the first step the mock catalogues containing stellar masses have been created from the MultiDark N-body simulation (Prada et al. 2012) and Pinocchio halo lightcones (Monaco et al. 2002). We followed the stellar mass to halo mass ratio (SMHR) approach based on the assumption of a monotonic relation between halo/subhalo masses and the stellar masses of the galaxies associated with them. We first populated the haloes in lightcones with subhaloes. For this we randomly distributed subhaloes around each distinct halo following a NFW profile (Navarro et al. 1997), so that their number density satisfies the subhalo mass function proposed by Giocoli et al. (2010). Then we assigned a galaxy to each halo and subhalo, with a stellar mass given by the SHMR of Moster et al. (2013).

The procedure which was followed in the next step is similar to the one used in the VVDS, zCOSMOS and VIPERS surveys also based on VIMOS observations as the VUDS survey (Meneux et al. 2006; Iovino et al. 2010; de la Torre et al. 2011b, 2013). In order to obtain fully realistic VUDS mock catalogues we add the detailed survey selection criteria (see Sect. 2.1 and Le Fèvre et al. 2015). This procedure in the end produces mock parent photometric galaxy catalogues. We then apply the slit-positioning algorithm (SSPOC) with the same settings as used for the VUDS survey (Bottini et al. 2005). This allowed us to reproduce the VUDS footprint on the sky for a proper simulation of the small-scale angular incompleteness in the mocks.

As a result of this procedure, 66 realistic mock galaxy catalogues have been produced. They contain the detailed survey completeness function and observational biases imparted when taking spectra in the VUDS survey.

3. Clustering measurements

3.1. Two-point correlation function

The two point correlation function is one of the most frequently used statistical tool used to measure and investigate galaxy clustering. It measures the probability, above Poisson, of finding two galaxies separated by a given distance r (Peebles 1980), $\begin{array}{ccc}\mathit{P}\mathrm{=}\mathrm{[}\mathrm{1}\mathrm{+}\mathit{\xi }\mathrm{(}\mathit{r}\mathrm{\left)}\mathrm{\right]}\mathit{\rho }\mathrm{d}{\mathit{V}}_{\mathrm{1}}\mathrm{d}{\mathit{V}}_{\mathrm{2}}\mathit{,}& & \end{array}$(1)where P is the probability of finding these two galaxies in two infinitely small volumes dV₁ and dV₂ separated by the distance r, while the average density of galaxies is given by ρ.

Although the idea is relatively simple, this definition cannot be straightforwardly applied to compute the correlation function from real data samples because of the limitations of galaxy surveys themselves. In principle, to retrieve ideal correlation functions one would need an unlimited survey which covers the whole sky, and includes all galaxies. Naturally, in practice such surveys are unreachable. A number of estimators of ξ(r), aimed at minimizing the effects related to the limited number of objects and limited areas covered by available surveys, have been proposed (see e.g. Davis & Peebles 1983; Hamilton 1993). The estimator the most commonly used because of its well established capability to minimize the above mentioned observational limitations is the Landy & Szalay (1993) estimator, which we also apply in this work: $\mathit{\xi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{R}}\mathrm{(}{\mathit{N}}_{\mathrm{R}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{GG}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{G}}\mathrm{(}{\mathit{N}}_{\mathrm{G}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{RR}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{-}\mathrm{2}\frac{\mathrm{(}{\mathit{N}}_{\mathrm{R}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{GR}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{G}}\mathit{RR}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{+}\mathrm{1.}$(2)Here, N_G and N_R represent, respectively, numbers of galaxies in the galaxy sample and randomly distributed objects in the same volume as observed in the survey; GG(r) is the number of distinct galaxy-galaxy pairs with separations lying in the interval (r,r + dr). Similarly, RR(r) and GR(r) are the numbers of random-random pairs and galaxy-random pairs, respectively, in the same intervals. The galaxy pairs GG, GR and RR are normalized by N_G(N_G−1), N_GN_R and N_R(N_R−1) respectively, where N_G and N_R are the number of galaxies in the data and random catalogues.

In this work we make galaxy clustering measurements using combined data from two independent VUDS survey fields - COSMOS and VVDS2h. The final ξ^final(r) is computed for all fields simultaneously by using Landy& Szalay estimator (Eq. (2)), and the differences in size and numbers between the fields are accounted for by an appropriate weighting scheme. In particular, each pair is multiplied by the number of galaxies per unit volume for the given field, ${\mathit{\xi }}^{\mathrm{final}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\sum _{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{{\mathit{n}}_{\mathrm{field}}}{\mathit{w}}_{\mathit{i}}\mathrm{·}(\mathit{G}{\mathit{G}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}\mathrm{2}\mathit{G}{\mathit{R}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\mathit{R}{\mathit{R}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)})}{\sum _{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{{\mathit{n}}_{\mathrm{field}}}{\mathit{w}}_{\mathit{i}}\mathrm{·}\mathit{R}{\mathit{R}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(3)where ${\mathit{w}}_{\mathit{i}}\mathrm{=}{{}^{\mathrm{(}}\mathit{N}_{\mathrm{gal}\mathit{,i}}\mathit{/}{\mathit{V}}_{\mathit{i}}{}^{\mathrm{)}}}^{\mathrm{2}}$.

The Landy & Szalay estimator requires the creation of a random catalogue, which follows the geometrical properties of the corresponding galaxy sample and is distributed in the same volume (or on the same area on the sky, in its angular version). In our method the random catalogue is created separately for each field in order to mirror the real redshift (not only angular) distribution of the data in a given field. The composite correlation function ξ_final (see Eq. (3)) is then estimated using random-random pairs from both random catalogues merged with the weighting scheme applied. The number of objects generated within this catalogue is also crucial – in practice it should be significantly more numerous than the real data sample. When there are too few random object, shot noise is introduced into the ξ(r) measurements, especially on small scales (r_p< 3 h^-1 Mpc). After a set of tests, for the measurements presented in this paper we decided to use N_R = 100 000, for which additional noise was found to be negligible

Another problem one has to face while measuring the real space correlation function is related to peculiar velocities of galaxies. In redshift space, these peculiar velocities distort distances computed from the redshifts z, and, as a result, they affect the shape of the correlation function ξ(r) itself. The corresponding distortions are known as the coherent infall and “fingers of God” (Kaiser 1987). These differences, however, apply only to the radial separations and do not have any influence on the measured position of a galaxy on the sky. A commonly applied method to eliminate them is to split the comoving redshift space separations into two components – parallel π and perpendicular r_p to the line of sight – thus re-defining the redshift-space correlation function as ξ(r_p,π). Integrating ξ(r_p,π) along the line of sight gives us a projected correlation function w_p(r_p), which is the two-dimensional counterpart of the real-space correlation function, free from the redshift-space distortions (Davis & Peebles 1983): ${\mathit{w}}_{\mathrm{p}}\mathrm{\left(}{\mathit{r}}_{\mathrm{p}}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\xi }\mathrm{\left(}{\mathit{r}}_{\mathrm{p}}\mathit{,\pi }\mathrm{\right)}\mathrm{d}\mathit{\pi }\mathrm{=}\mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\xi }\left(\mathrm{(}{\mathit{r}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{d}\mathit{y}\mathit{.}$(4)Here, y is the real-space separation along the line of sight, and ξ(r) stands for real space correlation function computed for r=$\sqrt{{\mathit{r}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$. In practice, the upper integral limit π_max has to be finite in order to avoid adding noise to the result. After performing a number of tests for different π_max values in the range 15 <π_max< 35h^-1 Mpc, we have chosen π_max = 20h^-1 Mpc, which is the maximum value for which the correlation function measurement is not significantly affected by noise.

3.2. Error estimates

Estimating errors on the two-point correlation function is the subject of numerous discussions in the literature since the time of the very first measurements (see Hamilton 1993; Fisher et al. 1994; Bernardeau et al. 2002). The usage of a properly constrained covariance matrix is needed to account for the fact that the values of w_p for different separations r_p are not independent – pair counts in different r_p bins include partially the same galaxies. For the VUDS data we decided to apply a combined method, which makes use of the so-called blockwise bootstrap re-sampling (Barrow et al. 1984) coupled to mock catalogues (see Sect. 2.3), similar to the method proposed by Pollo et al. (2005). This decision has been motivated by a number of tests performed on 66 mock catalogues. Errors were computed by using separately (1) bootstrap re-sampling; (2) mocks and (3) the combination of these two in order to check which method results in the least disturbed covariance matrix. It appeared that the last method, i.e. the combination of mocks and bootstraps, leads to the least noisy covariance matrix.

The blockwise bootstrap re-sampling is based on the idea of perturbing the data set by randomly creating a large number of comparable pseudo data sets, which differ only slightly from the original sample. A number of objects are randomly selected from the data sample (objects are allowed to be drawn multiple times), and the correlation function is computed for each of these subsamples. This procedure is repeated N_boot times, giving as a result the variation around the real result, which is used as an error estimate.

The second method makes use of sets of independent mock surveys – the simulated catalogues based on large N-body simulations coupled with physical definitions of galaxies, inserted e.g. by semi-analytic models. The correlation functions computed for each independent mock catalogue give the variation of the result at each scale, and are used as error estimates.

For both the bootstrap re-sampling and mock methods the associated covariance matrix C between the values of w_p on ith and kth scale can be computed as ${{C}}_{\mathit{ik}}\mathrm{=}⟨\mathrm{\left(}{\mathit{w}}_{\mathrm{p}}^{\mathit{j}}\mathrm{\right(}{\mathit{r}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathrm{⟨}{\mathit{w}}_{\mathrm{p}}^{\mathit{j}}\mathrm{(}{\mathit{r}}_{\mathit{i}}\mathrm{\left)}{\mathrm{⟩}}_{\mathit{j}}\mathrm{\right)}\mathrm{\left(}{\mathit{w}}_{\mathrm{p}}^{\mathit{j}}\mathrm{\right(}{\mathit{r}}_{\mathit{k}}\mathrm{)}\mathrm{-}\mathrm{⟨}{\mathit{w}}_{\mathrm{p}}^{\mathit{j}}\mathrm{(}{\mathit{r}}_{\mathit{k}}\mathrm{\left)}{\mathrm{⟩}}_{\mathit{j}}\mathrm{\right)}{⟩}_{\mathit{j}}\mathit{,}$(5)where “⟨⟩” indicates an average over all bootstrap or mock realizations, the ${\mathit{w}}_{\mathrm{p}}^{\mathit{j}}\mathrm{\left(}{\mathit{r}}_{\mathit{k}}\mathrm{\right)}$ is the value of w_p computed at r_p = r_i in the cone j, where 1 <j<N_mock for the VUDS mocks and 1 <j<N_boot for the bootstrap data.

In our case the covariance matrix reconstructed from N_mock = 66 VUDS mock catalogues could not be directly applied to the observed data because it caused the fit to be often unstable, and not able to properly converge. This instability was caused by the fact that the diagonal elements of the matrix have realistic values but the off-diagonal non-zero elements differ significantly from those pertaining to the data sample. In other words, the detailed statistical properties of the VUDS mock catalogues are not close enough to the real VUDS data to allow for such an operation directly. For this reason we computed error bars using the scatter between the VUDS mock catalogues (since we believe it is more realistic than the scatter between the bootstrap realizations), but the off-diagonal elements of the covariance matrix were computed using the classical bootstrap method. This means that for each redshift range we measured the correlation function w_p(r_p) from 1) the VUDS data; 2) N_mock = 66 VUDS mock catalogues; 3) N_boot = 100 bootstrap re-samplings of the VUDS data (to determine this number we increased N_boot until we reached the value above which the estimated errors do not change even if we increase a number of bootstrap representations even further) and in the fitting and error estimation we used the combination of these measurements.

3.3. Systematics in the CF measurement

Before introducing correlation function measurements from our data several tests have been performed on mock catalogues to investigate the influence of various VUDS survey properties on the correlation function. As the result of these tests, the correlation function correction scheme has been performed as described in the Appendix A.

In the lower panel of Fig. A.1 it is shown that the final measurement of the correlation function w_p(r_p), obtained after introducing the full correction scheme, is still slightly underestimated, by as much as ~ 10% on small scales, with respect to the true correlation function obtained for the mock parent samples (see details in Appendix A). This systematic underestimation may affect the value of the correlation length r₀, and all HOD parameters including an underestimation of the average halo mass ⟨ M_h ⟩ with respect to the true underlying values. To estimate the size of this effect, we performed a power-law function fit on the true correlation function (computed for the parent mock) as well as on the observed (the mock catalogue after SSPOC selection), with the full correcting scheme, both measured as the average from 66 VUDS mock catalogues. We found that the true correlation length is on average larger than the observed r₀ by Δr₀ = 0.42 ± 0.34, and that the observed slope was underestimated by Δγ = 0.06 ± 0.04. This, of course, would affect the galaxy bias estimation, since it is derived directly from these parameters, which we estimated to be underestimated by Δb = 0.35 ± 0.27.

To address this problem we compute the ratio ${\mathit{w}}_{\mathrm{p}}^{\mathrm{par}}\mathit{/}{\mathit{w}}_{\mathrm{p}}^{\mathrm{obs}}$ between the true ${\mathit{w}}_{\mathrm{p}}^{\mathrm{par}}\mathrm{\left(}{\mathit{r}}_{\mathrm{p}}\mathrm{\right)}$ and observed ${\mathit{w}}_{\mathrm{p}}^{\mathrm{obs}}\mathrm{\left(}{\mathit{r}}_{\mathrm{p}}\mathrm{\right)}$ correlation function (based on the VUDS mock, see Appendix A) on every scale r_p for each redshift range separately. The correlation function measurements are then multiplied by this ratio for each separation r_p, resulting in the final correlation function measurement w_p(r_p) on which the power-law and HOD fitting have been performed.

Another systematic underestimation of the correlation function may be related to possible wrong galaxy redshifts within the sample. It can by introduced mainly by the presence of the flag 2 galaxies in our sample (see Sect. 2.1 for details about redshift reliability flags). To estimate the possible effect of this uncertainty andcheck how a fraction of the sample with possibly wrongly measured redshifts could influence the final correlation function measurements, we performed a set of tests fully described in Appendix B.

Based on these tests, we conclude that the correlation function measurements performed in this study are robust against a possible contamination from incorrect redshift measurements. Assuming the worst possible case, i.e. the lower limit for confidence level of flag 2 galaxies being 70%, the resultant clustering strength r₀ will be underestimated only by Δr₀ = 0.13, corresponding to ~ 3.3% of the measured clustering amplitude. To take this effect into account, the upper error bar of r0 measurements is rescaled by the corresponding value.

All results presented in this paper include the correlation function corrections described above.

3.4. Power-law model

In most cases, especially in the local Universe, the correlation function ξ(r) is well described by a power law function ξ(r) = (r/r₀)^{− γ}, where r₀ and γ are correlation length and slope, respectively. With this parametrization, the integral in Eq. (4) can be computed analytically and w_p(r_p) can be expressed as ${\mathit{w}}_{\mathrm{p}}\mathrm{\left(}{\mathit{r}}_{\mathrm{p}}\mathrm{\right)}\mathrm{=}{\mathit{r}}_{\mathrm{p}}{\left(\frac{{\mathit{r}}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{p}}}\right)}^{\mathit{\gamma }}\frac{\mathrm{\Gamma }\mathrm{\left(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\right)}\mathrm{\Gamma }\mathrm{\left(}\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{\mathrm{2}}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}\frac{\mathit{\gamma }}{\mathrm{2}}\mathrm{\right)}}\mathit{,}$(6)where Γ is the Euler’s Gamma Function.

The values of w_p are not independent at different separations, hence it is not possible to use simple χ² minimization to find the best-fit parameters. However, the covariance matrix C is symmetric and real, and therefore it can be inverted (if it is not singular). We can fit w_p by minimizing a generalized χ², which is defined as ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i,j}}^{{\mathit{N}}_{\mathrm{D}}}\mathrm{\left(}{\mathit{w}}_{\mathrm{p}}^{\mathrm{obs}}\mathrm{\right(}{\mathit{r}}_{\mathit{i}}\mathrm{)}\mathrm{-}{\mathit{w}}_{\mathrm{p}}^{\mathrm{mod}}\mathrm{(}{\mathit{r}}_{\mathit{i}}{\mathrm{\left)}}^{\mathrm{\right)}}{{C}}_{\mathit{ij}}^{-1}\mathrm{\left(}{\mathit{w}}_{\mathrm{p}}^{\mathrm{obs}}\mathrm{\right(}{\mathit{r}}_{\mathit{j}}\mathrm{)}\mathrm{-}{\mathit{w}}_{\mathrm{p}}^{\mathrm{mod}}\mathrm{(}{\mathit{r}}_{\mathit{j}}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{.}$(7)Here the covariance matrix C is computed using the method described in Sect. 3.2 (Eq. (5)). The fitting procedure to estimate the power-law parameters r₀ and γ for the projected correlation function w_p(r_p) follows Fisher et al. (1994) (see also Guzzo et al. 1997 and Pollo et al. 2005).

A more detailed description of the shape of the real-space correlation function can be done e.g. in the framework of HOD models, which we discuss in the following Sect. 3.5. Nevertheless, the use of a power-law model to describe the correlation function remains an efficient and simple approximation of galaxy clustering properties.

3.5. The halo occupation distribution (HOD) model

As we mentioned above, in practice ξ(r) is often well fitted by a power law. However, there is no strict theoretical reason to force the galaxy correlation function to assume a power law shape. Indeed, recent measurements show deviations from a power law shape. The strongest deviations are observed on small scales, especially for the most luminous galaxies (Coil et al. 2006; Pollo et al. 2006; Zehavi et al. 2011), while probably the most famous ones are baryonic acoustic oscillations (BAOs) at large scales (Eisenstein et al. 2005; Percival et al. 2010; Anderson et al. 2014).

The small-scale behaviour of the correlation function of galaxies can be well interpreted in the framework of the HOD models based on the relations between the distribution of dark matter and galaxies. The HOD models describe bias in terms of the probability P(N | M) that a dark matter halo of a mass M contains N galaxies of a given type. Recently this approach has been very successfully used to model the shape of the two-point correlation function (e.g. Skibba et al. 2009; Abbas et al. 2010; Zehavi et al. 2011; Coupon et al. 2012; Kim et al. 2014).

In this work we apply the halo occupation model following the current commonly used analytical prescriptions, so that our results can be easily compared to results in the literature. We describe below the most important features and components of the HOD model we use, while for the general overview of the HOD philosophy we encourage readers to look into the review paper by Cooray & Sheth (2002). A very thorough description is also given by Coupon et al. (2012).

In the HOD framework the correlation function can be split into two components. The one-halo component ξ^1h(r) dominates on scales smaller than the size of dark matter haloes (usually < 1.5 h^-1 Mpc), while the two-halo component ξ^2h(r) dominates on larger scales. The one-halo term arises from pairs of galaxies located within the same halo, while the two-halo terms is built by the pairs of galaxies hosted by different haloes. Consequently, the correlation function can be written as $\mathit{\xi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\xi }}^{\mathrm{1}\mathrm{h}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{\xi }}^{\mathrm{2}\mathrm{h}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(8)The properties of the first term, i.e. ξ^1h(r) can be described with the use of the halo occupation model which tries to follow the distribution of galaxies located in one dark matter halo. The second term ξ^2h(r), depends on the statistical properties of the large-scale dark matter density field and the distribution of dark matter haloes with respect to this density field.

The most important properties of the dark matter halo from this point of view are the halo mass function n(M,z), for which we adopted the fitting formula proposed by Tinker et al. (2008), the halo density profile ρ(r | M), for which we assumed the form described by Navarro et al. (1997), and the halo bias b_h(M,z,r), for which we use the formula proposed by Tinker et al. (2010) with the scale dependency from Tinker et al. (2005).

We parametrized our halo occupation model in the way used by Zehavi et al. (2005b) and motivated by Kravtsov et al. (2004). In particular, we express the halo occupation function, i.e. the number of galaxies per halo, as a sum of a central galaxy and satellite galaxies. The mean occupation function for central galaxies N_c(M) is represented by a step function, while the mean halo occupation function for satellite galaxies N_s(M) is approximated by a Poisson distribution with the mean being a power-law, $\begin{array}{ccc}\mathrm{⟨}{\mathit{N}}_{\mathrm{g}}\mathrm{|}\mathit{M}\mathrm{⟩}& \mathrm{=}& \mathrm{1}\mathrm{+}{\left(\frac{\mathit{M}}{{\mathit{M}}_{\mathrm{1}}}\right)}^{\mathit{\alpha }} \mathrm{for} \mathit{M}\mathit{>}{\mathit{M}}_{\min }\\ & \mathrm{=}& \mathrm{0} \mathrm{otherwise}\mathit{,}\end{array}$(9)where M_min is the minimum mass needed for a halo to host one central galaxy, and M₁ is the mass of a halo having on average one satellite galaxy, while α is the power law slope of the satellite mean occupation function.

The HOD model used in this work is simplified in comparison with the ones widely used at lower redshift ranges, e.g. by Zehavi et al. (2011) for fitting the SDSS data. However, based on the available statistics of the VUDS sample, with lower numbers of galaxies than e.g. in the SDSS, the three parameter model appears to be the best solution to retrieve robust measurements of M_min and M₁. A similar approach was taken by Abbas et al. (2010) in the case of the measurements from the VVDS survey at z ~ 1 with a similar sample size (Abbas et al. 2010). Since the details of the interplay between galaxy evolution and large-scale structure growth at epochs as early as the ones probed by VUDS are still poorly understood, it is reasonable to keep a simplified description and not to dilute our analysis with models introducing too many free variables. The main results presented in this paper, however, are not expected to differ much from results where more complex HOD models are applied.

From the best-fit HOD parameters it is possible to obtain quantities describing halo and galaxy properties, like the average host halo mass M_h, $\mathrm{⟨}{\mathit{M}}_{\mathrm{h}}\mathrm{|}\mathit{g}\mathrm{⟩}\mathrm{(}\mathit{z}\mathrm{)}\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{M}\mathit{M}\mathit{n}\mathrm{(}\mathit{M,z}\mathrm{)}\frac{\mathrm{⟨}{\mathit{N}}_{\mathrm{g}}\mathrm{|}\mathit{M}\mathrm{⟩}}{{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathit{,}$(10)and the large-scale galaxy bias ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}$, ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{M}{\mathit{b}}_{\mathrm{h}}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathit{n}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\frac{\mathrm{⟨}{\mathit{N}}_{\mathrm{g}}\mathrm{|}\mathit{M}\mathrm{⟩}}{{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathit{,}$(11)where n(M,z) is the dark matter mass function, b_h(M,z) is the large-scale halo bias, and the n_g given by ${\mathit{n}}_{\mathrm{g}}\mathrm{=}\mathrm{\int }\mathit{n}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{⟨}{\mathit{N}}_{\mathrm{g}}\mathrm{|}\mathit{M}\mathrm{⟩}\mathrm{d}\mathit{M}$(12)represents the number density of galaxies.

3.6. Fitting the HOD model to the correlation function

To explore the HOD parameter space we implemented our HOD model with the publicly available CosmoPMC¹ code, which uses the Population Monte Carlo (PMC) technique to sample likelihood space (Wraith et al. 2009; Kilbinger et al. 2011). PMC is the adaptive importance-sampling technique (Cappé & Moulines 2007) which allows efficient sampling of the parameters space for the large number of samples. For each galaxy sample we fit the projected correlation function w_p(r_p) and the number density of galaxies n_g, by summing both contributions in log-likelihood obtained by $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{i,j}}^{{\mathit{N}}_{\mathrm{D}}}\mathrm{(}{\mathit{w}}_{\mathrm{p}}^{\mathrm{obs}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}\mathrm{-}{\mathit{w}}_{\mathrm{p}}^{\mathrm{mod}}\mathrm{(}{\mathit{r}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{)}}{{C}}_{\mathit{ij}}^{-1}\mathrm{(}{\mathit{w}}_{\mathrm{p}}^{\mathrm{obs}}\mathrm{\left(}{\mathit{r}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}{\mathit{w}}_{\mathrm{p}}^{\mathrm{mod}}\mathrm{(}{\mathit{r}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{)}}\\ & & \mathrm{+}\frac{{\mathrm{(}{\mathit{n}}_{\mathrm{g}}^{\mathrm{obs}}\mathrm{-}{{\mathit{n}}_{\mathrm{g}}^{\mathrm{mod}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathit{n}}_{\mathrm{g}}}^{\mathrm{2}}}\mathit{,}\end{array}$(13)where ${\mathit{n}}_{\mathrm{g}}^{\mathrm{mod}}$ is given by Eq. (12) at the mean redshift of the sample. The data covariance matrix is approximated taking Eq. (5). The error on the galaxy number density σ_ng contains Poisson noise and cosmic variance.

3.7. Large-scale galaxy bias

According to the current cosmological paradigm of structure formation, galaxies form and evolve inside dark matter haloes (White & Rees 1978). In other words, there exists a connection between the dark matter distribution and galaxies in the dense dark matter regions, where the halo clustering is important, galaxies should be more clustered. The galaxy spatial distribution, however, is biased with respect to the dark matter density field. The strength of this effect is referred to as galaxy bias.

As a first approximation we use a linear model for the large-scale galaxy bias. It assumes a linear relation between galaxy σ_R,g(z) and mass σ_R,m rms at a given redshift, ${\mathit{\sigma }}_{\mathit{R,}\mathrm{g}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{b}{\mathit{\sigma }}_{\mathit{R,}\mathrm{m}}\mathrm{\left(}\mathit{z}\mathrm{\right)}$(14)where b is the galaxy bias. We also assume that b is independent of scale R, which is true especially for large scales. Usually the adopted value for the scale on which σ is measured is R = 8 h^-1 Mpc. Locally, at z = 0, the mass fluctuations reach the value σ_8,m(z = 0) = 0.83 (Planck Collaboration XVI 2014). In our model the redshift evolution of this quantity is described as ${\mathit{\sigma }}_{\mathrm{8}\mathit{,}\mathrm{m}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathrm{8}\mathit{,}\mathrm{m}}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{D}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}$(15)where $\mathit{D}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathit{g}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{g}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}$(16)and g(z) is the normalized growth factor, which describes how fast the linear perturbations grow with the scale factor. Its value depends on the assumed cosmological model, i.e. cosmological parameters. After Carroll et al. (1992) we write $\mathit{g}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{2.5}{\mathrm{\Omega }}_{\mathrm{m}}{\left[{\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{4}\mathit{/}\mathrm{7}}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{\Lambda }}\mathrm{+}\left(\mathrm{1}\mathrm{+}\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{\mathrm{2}}\right)\left(\mathrm{1}\mathrm{+}\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{\mathrm{70}}\right)\right]}^{-1}\mathit{,}$(17)where Ω_m and Ω_Λ are the matter and dark energy density parameters at a given redshift.

The galaxy rms at a given scale R can be retrieved from the power-law fits to the correlation function as ${\mathit{\sigma }}_{\mathrm{8}\mathit{,g}}\mathrm{=}\sqrt{{\mathit{C}}_{\mathit{\gamma }}{\left(\frac{{\mathit{r}}_{\mathrm{0}}}{\mathrm{8}{\mathit{h}}^{-1}\hspace{0.17em}\mathrm{Mpc}}\right)}^{\mathit{\gamma }}}\mathit{,}$(18)with ${\mathit{C}}_{\mathit{\gamma }}\mathrm{=}\frac{\mathrm{72}}{{\mathrm{2}}^{\mathit{\gamma }}\mathrm{(}\mathrm{3}\mathrm{-}\mathit{\gamma }\mathrm{)}\mathrm{(}\mathrm{4}\mathrm{-}\mathit{\gamma }\mathrm{)}\mathrm{(}\mathrm{6}\mathrm{-}\mathit{\gamma }\mathrm{)}}\mathit{,}$(19)where r₀ and γ are the correlation length and the slope of power-law approximation of the correlation function.

In the following sections, we refer to the large-scale galaxy bias computed by this method as ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}$. This is also used to compute the values based on measurements of the correlation function from other surveys when the b_L estimations were not presented by the authors; in these cases we use their power-law best fit parameters r₀ and γ (e.g. from Le Fèvre et al. 2005). Another estimate of linear large-scale galaxy bias can be obtain from the HOD model. We compute the large-scale galaxy bias ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}$ using the fit of the halo occupation function (Eq. (11)).

The comparison of these two bias values and differences between these two procedures are discussed in Sect. 5.

4. Results

4.1. The VUDS correlation function – general population

The projected two-point correlation function w_p(r_p) is computed for the galaxy population with i_AB< 25 at 2 <z< 5 separated into three redshift ranges (see Sect. 2.2). We are using a total number of 3022 spectroscopically confirmed VUDS galaxies located in two fields. In Fig. 3 we present composite correlation function w_p(r_p) measurements for the full data sample in the redshift range 2 <z< 5, along with the measurements made for each field independently. It is worth stressing that it is the first time that the correlation function has been computed with such a high accuracy from spectroscopic data at redshift z ~ 3. As expected, the values of the correlation function measured for galaxies located in separate fields are slightly different due to cosmic variance. However, for both fields the measurements are consistent with respect to the error bars and the signal is retrieved on all scales 0.3 <r_p< 16 h^-1 Mpc, so that the composite function is not dominated by the signal from any of the individual fields at any separation r_p. The minimum and maximum scale r_p were chosen by testing different scale limits for which we were able to retrieve correlation function signal from our data. The lower r_p limit is set for the lowest scale for which the measurement is not underestimated due to the low number statistics. The maximum scale was chosen by testing different scale limits for which we were able to retrieve a correlation function signal and preserve a reasonable number of bins (in r_p) necessary to fit the HOD model at the same time.

The correlation functions for the full sample z = [2,5] and the two z = [2,2.9] and z = [2.9,5] redshift intervals are presented in Fig. 4.

Fig. 3 Projected two-point correlation function wp(rp) in each individual VUDS field for the redshift range 2 <z< 5. Blue and red lines correspond to the VVDS and COSMOS field, respectively. Black points and line indicate the combined correlation function from the measurements preformed on galaxies from both fields simultaneously.

Fig. 4 Left panel: projected two-point correlation function wp(rp). The symbols and error bars denote measurements of the composite correlation function from the VUDS survey. Right panel: associated best-fit power-law parameters r0 and γ along with error contours for the general galaxy population separated in three redshift ranges. Top panels present a measurement for a wide redshift range 2 < z < 5 sample, while lower panels present the correlation function measurements for two independent low and high redshift galaxy samples.

Fig. 5 Correlation length r0 as a function of redshift from various surveys of different types of objects. Stars indicate measurements obtained from this work in three redshift ranges (magenta for a wide redshift range sample, and black for a low and high redshift samples). Different colours indicate different galaxy populations targeted by using various techniques, as indicated in the upper right corner. Open symbols indicate measurements based on photometric data, while filled symbols are for measurements from spectroscopic data. Blue: LBG galaxies (open squares – Foucaud et al. 2003; open circles – Ouchi et al. 2004; open triangles – Adelberger et al. 2005; open reversed triangles – Kashikawa et al. 2006; open diamonds – Savoy et al. 2011; filled diamonds – Bielby et al. 2013; open pentagon – Barone-Nugent et al. 2014). Orange: BzK galaxies (open circles – Blanc et al. 2008; open triangles – Hartley et al. 2010; open reversed triangle – McCracken et al. 2010; open diamonds – Lin et al. 2012). Green: galaxy samples from surveys limited in luminosity (filled squares – Norberg et al. 2002; open circles – Coil et al. 2006; filled triangles – Le Fèvre et al. 2005; filled reversed triangles – Pollo et al. 2006, filled diamonds – Zehavi et al. 2011; filled pentagons – Marulli et al. 2013; crosses – Skibba et al. 2009). Red: EROs or massive red galaxies (open squares – Daddi et al. 2003; filled squares – Zehavi et al. 2005a; open circles – Brown et al. 2008). The five solid curves show the correlation length of dark haloes with different minimum masses, as labelled.

4.2. General population – Power law modelling

We fit a power-law to w_p(r_p), with two free parameters r₀ and γ, to quantify the clustering strength. All the points in the range 0.3 h^-1 Mpc <r_p< 16 h^-1 Mpc are used (see Sect. 3.4 for the fitting method details). The measured best-fit parameters are listed in Table 2 and the power-law fit of the correlation function for each redshift subsample is shown in Fig. 4 together with the (r₀,γ) error contours.

Within our sample we find ${\mathit{r}}_{\mathrm{0}}\mathrm{~}\mathrm{3.9}{\mathrm{5}}_{-0.54}^{\mathrm{+}\mathrm{0.48}} {\mathit{h}}^{-1}\hspace{0.17em}\mathrm{Mpc}$ in z = [2.0,2.9] and ${\mathit{r}}_{\mathrm{0}}\mathrm{~}\mathrm{4.3}{\mathrm{5}}_{-0.76}^{\mathrm{+}\mathrm{0.60}}{\mathit{h}}^{-1}\hspace{0.17em}\mathrm{Mpc}$ in z = [2.9,5]. The slight increase in r₀ with redshift is probably due to the small luminosity differences between samples (see Sect. 5 for discussion), although the difference is marginally significant given measurement errors. For these samples, the slope γ varies between $\mathit{\gamma }\mathrm{=}\mathrm{1.6}{\mathrm{0}}_{-0.13}^{\mathrm{+}\mathrm{0.12}}$ and $\mathit{\gamma }\mathrm{=}\mathrm{1.8}{\mathrm{1}}_{-0.06}^{\mathrm{+}\mathrm{0.02}}$, hence showing a tendency to decrease with redshift.

4.3. General population - HOD modelling

The best-fit HOD model defined in Sect. 3.5 is fit to each of three redshift subsamples as presented in the left panels of Fig. 6. The HOD model provides a good fit to the real-space correlation function w_p(r_p) in all the redshift ranges, particularly in z = [2,2.9] where it is an excellent representation of the data. For the highest redshift bin z = [2.9,5] the fit is not as good. This is mainly due to the less accurately constrained one-halo term of the correlation function in this redshift range in comparison to the other subsamples.

The three HOD parameters, the minimum halo mass M_min, the satellite occupation halo mass M₁, and the high halo mass slope α (see Sect. 3.5), are inferred from the full covariance matrix (see Sect. 3.6) and are listed in Table 3. The parameters are allowed to vary within a large range set from previous observations at high and low redshift: 10 < log M_min< 14, 10 < log M₁< 15, 0.6 <α< 2.0. The right panel of Fig. 6 presents the halo occupation function obtained from the HOD fit within three redshift ranges. We do not observe any significant difference in halo masses between the two redshift ranges probed for the general VUDS galaxy population. The minimum halo mass, for which a halo hosts at least one central galaxy M_min, has comparable value for the z ~ 2.5 and z ~ 3.5 samples.

The satellite occupation mass M₁ is noticeably higher in the higher redshift bin than in the lower redshift sample, although the errors are quite large. These uncertainties are related to the weak one-halo term signal in the correlation function for the higher redshift measurement.

The high halo mass slope α of the satellite mean occupation function is around ~ 1.3 for z = [2,2.9], significantly steeper than α = 0.73 found in z = [2.9,5].

The average host halo mass ⟨ M_h | g ⟩, large-scale galaxy bias ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}$ and associated 1σ errors, are reported in Table 3. The redshift evolution of these parameters is presented in Fig. 7 – for the average host halo mass, and in Fig. 8 – for the large-scale galaxy bias. These results are fully discussed in Sect. 5.

Fig. 6 Left panels: correlation function for three redshift subsamples as in Fig. 4. Top panel represents measurements for a wide redshift range 2 <z < 5 sample, while lower panels present correlation function measurements for two independent low and high redshift galaxy samples. The solid lines present the best-fit halo occupation model. Right panel: the evolution of the halo occupation for the three redshift ranges.

Fig. 7 Evolution of the number-weighted average host halo mass given by Eq. (10) for the three redshift ranges analysed in this study. The red filled square indicate mass estimations from VUDS for the wide range redshift sample, while red filled circles for low and high redshift ones. Black and grey symbols represent the results of previous work based on spectroscopic and photometric surveys respectively. The solid black lines indicate how a host halo of a given mass M0 at z = 0 evolves with redshift, according to the model given by van den Bosch (2002). The solid red line represents the halo mass evolution derived using Eq. (21), with the HOD parameters obtained from the best-fit HOD model at a redshift z ~ 3. The dashed red line is using the HOD best-fit parameters for z ~ 2.5. VUDS galaxies with a typical L⋆ luminosity are likely to evolve into galaxies with a luminosity >L⋆ today.

4.4. Linear large-scale bias measurements

The large-scale linear galaxy bias ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}$ has been computed using the method described in Sect. 3.7, with the best-fit correlation parameters obtained for each galaxy subsample. The bias values obtained and the associated 1σ errors are listed in Table 2. We find that the galaxy bias is increasing from 2.39 ± 0.32 at z ~ 2.5 to 3.26 ± 0.47 at z ~ 3.5, galaxies in the higher redshift sample appearing to be more biased than those in the lower redshift sample. Additionally, we note that the bias reaches a relatively high value at z ~ 3.5. This is further discussed in Sect. 5.

5. Discussion

5.1. Evolution of the clustering length r₀ from z ~ 3.5 until the present epoch

In this paragraph we discuss the evolution of clustering properties of the general galaxy population since a redshift z ~ 3.5 (see Sect. 4) in the view of other measurements presented in the literature.

Within our sample, we observe the clustering strength to increase slightly with redshift over the range z = [2.0−5.0], with r₀ being the highest at the highest redshift z ~ 3.5. However, this increase is marginally significant given the measurement errors. The observed increase in r₀ with z, surprising at the first glance, may be caused by the luminosity difference between two samples at z ~ 2.5 and z ~ 3.5. Because of the magnitude selection i_AB ≤ 25 of the VUDS sample we observe systematically more luminous galaxies at higher z. As a result, our high redshift sample contains on average more luminous galaxies, with median M_U = −21.81, while for the lower redshift sample the median M_U = −21.31. The luminosity dependence of galaxy clustering is a well established fact; it has been observed both at early and late cosmic epochs (e.g. Zehavi et al. 2011; Coupon et al. 2012; Abbas et al. 2010; Pollo et al. 2006) and luminous galaxies at all z investigated so far were found to be more strongly clustered than the fainter ones which is usually related to the fact that they are located in more massive dark matter haloes and, consequently, stronger fluctuations of the dark matter density field. If this trend is also true for the redshift range investigated in this work, the observed increase in r₀ with z may be caused by the fact that our higher redshift sample is on average more luminous that the lower one. The dependency of clustering on luminosity at these high redshifts will be investigated in a future paper (Durkalec et al., in prep.).

The correlation power-law parameters r₀ = 3.97 ± 0.48 and γ = 1.70 ± 0.09 at z ~ 3 are comparable with the clustering parameters for the star forming galaxies at intermediate (z ~ 1) redshift (e.g. Le Fèvre et al. 2005; Pollo et al. 2006), but the clustering length is lower than found for galaxies in the local Universe (e.g. Zehavi et al. 2011; Norberg et al. 2002). This has to be placed in the global perspective of the clustering expected to grow stronger along cosmic time.

A direct comparison of clustering strength at different epochs is a difficult task for several reasons. Different measurements are based on various galaxy samples, each with a unique selection function. Since galaxies cluster differently depending on their properties (like luminosity, colours or stellar masses), and in practice every catalogue contains a galaxy population somewhat different, the differences in clustering measurements may primarily be attributed to different galaxy sample properties as much as to evolutionary effects. The inhomogeneities between galaxy samples are especially pronounced for the high redshift samples. At z> 2 correlation functions have mostly been measured from surveys targeting specific classes of galaxies, e.g. extremely massive red objects or sources selected using a Lyman-break technique selecting star-forming galaxies. In addition, until recently most clustering measurements at z> 2 were produced from the angular correlation function computed on photometric samples with the knowledge of the redshift being only approximate and based on colour selection (BzK e.g. Daddi et al. 2003; McCracken et al. 2010; LBG e.g. Adelberger et al. 2005; LAE e.g. Ouchi et al. 2004). Compared to spectroscopic samples, the lower accuracy and precision in redshift measurements in these samples and contamination in a given redshift bin from galaxies at other redshifts introduces a blurring of large-scale structures which is difficult to overcome even when using large samples.

Fig. 8 Large-scale linear galaxy bias retrieved by using two methods. Red filled square indicate measurements of ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}$ for the wide range redshift sample, while the red filled circles for the low and high redshift samples. The open red square and circles represent ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}$ in three redshift ranges respectivley. Green and blue point represents the results from previous works, retrieved by using HOD (green circles – Blake et al. 2008; green diamonds – Wake et al. 2008; green reversed triangles – Ross & Brunner 2009; green polygons – Sawangwit et al. 2011; green squares – Kim et al. 2014; green triangles – Bielby et al. 2013) and power-law frameworks (blue squares – Marinoni et al. 2005; blue triangles – Adelberger et al. 2005; blue reversed triangles – Geach et al. 2012; blue diamonds – Barone-Nugent et al. 2014) We show 3 different models for the evolution of bias: a passive evolution model (red dashed line) a merging model (red dot-dashed line), a star-forming model (dotted line), as described in the Sect. 5.3.

Our r₀ estimations are in excellent agreement with other available results at similar high redshifts (e.g. Daddi et al. 2003; Ouchi et al. 2004; Adelberger et al. 2005; Kashikawa et al. 2006; Hartley et al. 2010; Savoy et al. 2011), as presented in Fig. 5. While r₀ in the VUDS sample is on the lower side of previous clustering measurements at these redshifts, most other samples have specific selection functions from which is it expected that they have a higher clustering bias than our sample, particularly for the redder or more luminous samples. Our measurements can be compared to the measurements by Bielby et al. (2013), based on a large sample of 2135 spectroscopically confirmed LBGs in ~ 10 000 arcmin² at z = 2.8, the largest sample with a spectroscopically-based measurement from a large volume at these redshift prior to our measurements. They find a clustering length r₀ = 3.46 ± 0.41 at z = 2.79 and a slope γ = 1.5−1.6, flatter than for local galaxies, in good agreement with our measurements given the associated errors. Savoy et al. (2011) measured a correlation length r₀ = 3.5 ± 0.5 from ~ 1400 UV-selected galaxies with photometric redshifts z ~ 3 and 23 <R< 27 in the Keck Deep Fields, and Geach et al. (2012) found r₀ = 3.7 ± 0.3 from the clustering of ~ 370Hα emitters (HAEs) at z = 2.23 from HiZELS, both results being very similar to ours at comparable redshifts. We conclude that our measurements are in general agreement with the most recent and least biased measurements available in the literature. Our results are among the most reliable today given the VUDS survey area and the number of galaxies used in our study.

5.2. The host halo mass

We draw in Fig. 5 the model predictions for the evolution of r₀ of dark matter haloes following the prescription given by Mo & White (2002) (which is based on the Press & Schechter (1974) formalism). According to this model, dark matter haloes of an average mass ⟨ M_h ⟩ ~ 10^11.5 M_⊙ have the same r₀ as observed VUDS galaxies at z ~ 3. If we take the rough assumption that at z ~ 3 haloes are occupied mostly by only one galaxy (central), this mass would be the average mass of a halo hosting the typical VUDS galaxy. However, this result is only a first approximation, since it is based on simplified assumptions of the Press & Schechter (1974) theory and needs to be refined.

The HOD model is expected to give more accurate predictions, since it accounts for the actual number and distribution of galaxies occupying haloes. We obtain the number weighted average mass ⟨ M_h ⟩ of haloes hosting VUDS galaxies using the best-fit HOD parameters (Eq. (10)), in three redshift ranges. We find that the halo mass decreases with redshift (see Fig. 7), evolving from $\log \mathrm{⟨}{\mathit{M}}_{\mathrm{h}}\mathrm{⟩}\mathrm{=}\mathrm{11.6}{\mathrm{1}}_{-0.40}^{\mathrm{+}\mathrm{0.46}}{\mathit{h}}^{-1}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$ at z ~ 3.5 to $\log \mathrm{⟨}{\mathit{M}}_{\mathrm{h}}\mathrm{⟩}\mathrm{=}\mathrm{12.0}{\mathrm{1}}_{-0.18}^{\mathrm{+}\mathrm{0.16}}{\mathit{h}}^{-1}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$ at z ~ 2.5. These estimations are in reasonable agreement with the recent observations by Bielby et al. (2013) who found an average mass of log ⟨ M_h ⟩ ~ 11.57 ± 0.15 M_⊙ for haloes hosting LBG galaxies at a redshift z ~ 2.79.

It can also be seen from Fig. 7 that our estimated host halo masses are much lower than the ones typically observed in the local Universe. At z ~ 3 the typical host halo masses are located in the range between ⟨ M_h ⟩ ~ 10¹¹−10¹² M_⊙ while locally these masses usually reach the values between 10¹³−10¹⁴ M_⊙. These results indicate that at some stage dark matter haloes must have experienced a rapid accretion phase, as expected in the framework of the hierarchical mass assembly.

A universal formula for the mass growth of cold dark matter haloes has been derived by van den Bosch (2002). The history of a halo with a given mass M₀ at z = 0 can be traced back in time using a simple formula, $\log \mathrm{⟨}\mathrm{\Psi }\mathrm{(}{\mathit{M}}_{\mathrm{0}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\mathrm{⟩}\mathrm{=}\mathrm{-}\mathrm{0.301}{\left[\frac{\log \mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}^{\mathrm{)}}}{\log \mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathit{f}}\mathrm{)}}\right]}^{\mathit{\nu }}\mathit{,}$(20)where z_f and ν are free fitting parameters that depend on halo mass and cosmological parameters. To obtain these fitting parameters we follow the analytical formula from van den Bosch (2002). We then trace the evolution of the galaxy population sampled by VUDS from z ~ 3 to the present epoch to predict the mass of haloes hosting the present day descendants of the VUDS galaxies. We find that in this model the typical VUDS halo with a mass ⟨ M_H ⟩ ~ 10^11.75 M_⊙ at z ~ 3 should evolve into a halo with a mass ⟨ M_H ⟩ ~ 10^13.5 M_⊙ at z = 0. In the local SDSS galaxy sample Zehavi et al. (2011) found that haloes with these masses are typically occupied by star forming galaxies with luminosity M_r< −20.5. According to the van den Bosch (2002) model the star forming galaxies at z ~ 3 in VUDS with a typical characteristic luminosity L_⋆ would likely evolve into galaxies equivalent to or brighter than L_⋆ at the present day.

The above comparison assumes that each halo is occupied by only one galaxy. This is not expected to generally be the case, and the above picture, while broadly correct might need to be adjusted. In order to trace the evolution of dark matter haloes and the hosted galaxy population in a more realistic way we use both the halo mass growth model Ψ(M₀,z) and the halo occupation function ⟨ N_g | M ⟩ at redshift z = 3. The average halo mass M_h as a function of redshift z is measured taking (see Eq. (10) for comparison) $\mathrm{⟨}{\mathit{M}}_{\mathrm{h}}\mathrm{⟩}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{M}{\mathrm{\Psi }}^{-1}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathit{n}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\frac{\mathrm{⟨}{\mathit{N}}_{\mathrm{g}}\mathrm{|}\mathit{M}\mathrm{⟩}}{{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathit{,}$(21)where the Ψ^-1(M,z) is the inverse mass growth function proposed by van den Bosch (2002) (Eq. (20)), n(M,z) is the dark matter mass function, and n_g the galaxy number density. This allows us to trace the history of the typical dark matter halo hosting the average VUDS galaxy from a redshift z = 3 to the present day, as presented in Fig. 7. We find that the typical halo of mass ⟨ M_h ⟩ = 10^11.75 at z ~ 3 should evolve into a halo of mass ⟨ M_h ⟩ = 10^13.9 at redshift z = 0. Comparing again our results from this improved model with Zehavi et al. (2011), haloes with these masses are typically occupied by star forming galaxies with M_r< −21.5, above the characteristic luminosity M_⋆ at z ~ 0. This means that the general galaxy population observed at redshift z = 3 evolves into the brightest and most massive galaxies, which occupy the most massive haloes observed in the local Universe, as is expected in the hierarchical mass growth paradigm.

Through this study it has been found that the VUDS star forming galaxies observed at redshift z ~ 3 typically occupy haloes of masses ~ 10¹² M_⊙, as presented in Table 3. Similar halo masses have been measured for star-forming galaxies across other redshift ranges: the halo mass 10¹² M_⊙ have been found for star forming galaxy populations at z ~ 0 and z ~ 1 (Zheng et al. 2007; Ross et al. 2010; Zehavi et al. 2011; Hartley et al. 2013). At the same time, the kinematic galaxy models (Shankar et al. 2006) predict the highest efficiency of galaxy star formation (more precisely, of conversion of primordial gas into stars) for similar halo masses ~ 10¹² M_⊙ – above this critical mass the feedback effects from stars and AGNs tend to suppress the star formation. Our measurement confirms that a halo mass of ~ 10¹² M_⊙ is also a characteristic halo mass scale for star forming galaxies up to redshifts z ~ 5.

5.3. Galaxy bias

The halo mass is directly related with the galaxy bias. The hierarchical model of galaxy formation predicts an increase in the galaxy bias with redshift Mo & White (1996). In this section we discuss the relation between the distribution of galaxies and the underlying dark matter density field at different redshifts, by comparing the galaxy bias values measured from VUDS to the bias of galaxy populations at different redshifts from the literature.

Figure 8 shows the linear large-scale galaxy bias b_L, computed using methods based on power-law ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}$ (see Sect. 3.7) and HOD ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}$ approximations (see Sect. 3.5). The results obtained from both methods are comparable within errors, and small differences are likely due to the less accurate power-law approximation. The general VUDS galaxy population shows a bias decreasing with decreasing redshift. According to our measurements, from z ~ 3.5 to z ~ 2.5 the linear bias decreased by Δb_L ~ 0.9. Our values at z = 2.5 and z = 3.5 are both significantly higher than measurements at later epochs reported in the literature. For comparison, we plot the linear large–scale bias inferred for the intermediate (z = [0.5−1.5]) redshift range by Marinoni et al. (2005), based on VVDS-Deep survey, for z> 2 by Geach et al. (2012) based on observations of Hα emitters, and by Barone-Nugent et al. (2014) for z ~ 3.8LBGs, along with the values we computed using the best-fit power-law parameters r₀ and γ from Le Fèvre et al. (2005); Adelberger et al. (2005); and Bielby et al. (2013, the method is described in Sect. 3.7). The values presented in Table 3 are similar to the bias measurements recently observed for LBGs at redshift z> 2. Bielby et al. (2013) found b_gal = 2.33 ± 0.17 from their LBG-selected spectroscopic sample, and at an even higher redshift z ~ 3.8Barone-Nugent et al. (2014) estimated b_L = 3.0 ± 0.2 for similarly selected galaxies.

The visible decrease in the galaxy bias with cosmic time can be explained in terms of the hierarchical scenario of structure formation. In this framework galaxies formed at early epochs in the most dense regions (highest peaks of the density field) which are the most biased with respect to the underlying average mass density field. As the mass density field evolved with time, galaxy formation is then expected to systematically move into less dense, hence less biased, regions. Additionally, it is likely that the gas in the most dense regions became too hot for collapse, and thus star formation could not take place in the strongest over-densities any more (e.g. Blanton et al. 1999). Therefore, the observed decrease in galaxy bias with cosmic time may be explained as galaxy formation is moving to less dense areas as cosmic time increases in a downsizing trend (De Lucia et al. 2006).

In Fig. 8 we also draw the predicted evolution of the large-scale galaxy bias using different prescriptions. We consider three different theoretical descriptions of the biasing functions, based on different ideas of how evolution may proceed: the passive, the merging, and the star forming biasing models (see e.g. Marinoni et al. 2005). In the first model the number of galaxies (halo occupation) is conserved as a function of time and the bias evolves assuming the form presented in Eq. (11). This model allows us to trace the expected bias of the galaxy population observed at z ~ 3 when evolved to lower redshifts. The two other models predict the galaxy bias which would be measured for the galaxy population of the same properties as the one observed at z ~ 3, but placed at different redshifts. The second model takes into account galaxy merging, as proposed by Mo & White (1996), who gave analytical prescriptions for computing the bias of haloes using the Press & Schechter formalism. If we assume that galaxies can be identified with dark matter haloes, an approximate expression for the biasing of all haloes of mass >M existing at redshift z (but which collapsed at redshift greater than the observation redshift, see discussion in Matarrese et al. 1997) is given by $\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{c}}}\left(\frac{{\mathit{\delta }}_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}}\mathrm{-}\mathrm{1}\right)\mathit{,}$(22)where δ_c ≃ 1.69 is the linear overdensity of a sphere which collapses in an Einstein-de Sitter Universe and σ(M,z) is the linear rms fluctuations on scales corresponding to mass M at the redshift of observation. The star-forming model assumes that the distribution of galaxies with luminosity >L is well traced by haloes with mass >M, and predicts the biasing of objects that just collapsed at the redshift of observation (e.g. Blanton et al. 2000). In this model, $\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{{\mathit{\delta }}_{\mathrm{c}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}}$(23)represents the biasing of galaxies that formed in a narrow time interval around redshift z (i.e. galaxies which experienced recent star formation at redshift z).

From the halo mass analysis presented in Sect. 5.2 we infer that VUDS galaxies would evolve into galaxies with a luminosity M_r< −21.5 in the local Universe (Sect. 5.2), for which the bias is b_L = 1.7. This is consistent with the predictions of the passive model in which VUDS galaxies with b_L ~ 2.8 at z = 3 would evolve to b_L = 1.7 at z ~ 0, which is indeed the value measured for a local galaxy population with M_r< −21.5. However, it is likely that more sophisticated models would need to be developed to reproduce the joint evolution of DM halo mass and galaxy bias. On the other hand, two remaining models predict that the galaxy population of properties similar to the VUDS galaxy population are now weakly biased class of galaxies.

6. Summary

We examined the galaxy clustering properties of high redshift galaxies with 2 <z< 5, by measuring the two-point real-space correlation function w_p(r_p) for 3022 galaxies from the VUDS survey. These measurements are among the first performed from spectroscopically confirmed galaxies at such high redshifts and have allowed clustering evolution and galaxy host haloes to be probed since very early times. We quantify our observations in terms of a power-law approximation and a three-parameter HOD model fitted to the correlation function computed in three redshift ranges, with median z ~ 2.5, z ~ 3.0 and z ~ 3.5. Our work therefore complements and extends to higher redshifts previous galaxy clustering analyses.

The results and conclusions can be summarized as follows:

• We observe a mild clustering evolution for the star-forming popu-lation in the VUDS survey. The correlation length increases from r₀ = 3.95 ± 0.48at z ~ 2.5 to reach r₀ = 4.35 ± 0.60 at z = 3.5. We attribute this slight increase to the differencein absolute magnitudes between these two samples, the higherredshift sample being more luminous. The slope of the corre-lation function is found to be $\mathit{\gamma }\mathrm{=}\mathrm{1.8}{\mathrm{1}}_{-0.06}^{\mathrm{+}\mathrm{0.02}}$ and $\mathrm{1.6}{\mathrm{0}}_{-0.13}^{\mathrm{+}\mathrm{0.12}}$ at z = 2.5 and 3.5 respectively, similar orslightly flatter than the canonical γ = 1.8 observed at low redshifts. Thesemeasurements are in general agreement with the most recent andleast biased measurements available in the literature at similarredshifts. We find that the power-law parameters we derive fromthe correlation function measurements are comparable to theclustering parameters for star-forming galaxies at intermediateredshifts, but the clustering length is somewhat lower than forgalaxies in the local Universe.

• The average halo mass obtained from our HOD model fit reaches ⟨ M_h ⟩ = 10^11.75 h^-1M_⊙ at redshift z = 3, which is smaller than observed locally. Assuming the mass growth model from van den Bosch (2002) these haloes should evolve into haloes with a mass ⟨ M_h ⟩ (z = 0) = 10^13.9 h^-1M_⊙ at z = 0. As haloes with such a mass are occupied by star forming galaxies with M_r< −21.5 in the local Universe. We infer that the star-forming galaxy population observed at a redshift z ~ 3 has evolved into the massive and bright galaxy population observed today, as predicted in the hierarchical mass growth paradigm.

• We found that the linear large-scale galaxy bias b_L, computed using both the HOD and power-law best-fit parameters, increases with redshift from ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}\mathrm{=}\mathrm{2.68}\mathrm{\pm }\mathrm{0.22}$ at z = 2.5 to ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}\mathrm{=}\mathrm{3.26}\mathrm{\pm }\mathrm{0.47}$ at z = 3.5. These bias values are comparable with measurements obtained at similar redshifts but are much higher than bias values observed at low and intermediate redshift ranges where b_L reaches values in the range 1−2. We discuss simple passive model to link the z ~ 3 galaxies in our sample to lower redshift populations, and find that this passive halo evolution model maintaining halo occupation parameters constant will lead to a bias value close to b_L = 1.7 at z ~ 0 we inferred from the halo mass evolution. However it is likely that more sophisticated models are needed to fully interpret this evolution.

Our data provide a robust reference for the clustering of galaxies at an epoch before the peak in star formation rate. We infer that the star-forming population of galaxies at z ~ 3 should evolve into the massive and bright (M_r< −21.5) galaxy population which typically occupies haloes of mass ⟨ M_h ⟩ = 10^13.9 at redshift z = 0. More complete modelling seems necessary to relate high redshift to low redshift galaxy populations, possibly taking into account non-linear effects, to better explain the galaxy clustering properties at early epochs in the life of the Universe.

Online material

Appendix A: Correlation function corrections scheme

Fig. A.1 Projected two point correlation function wp(rp) measured for 66 VUDS mock catalogs in three redshift ranges. The points correspond to the mean measurement of all 66 mock catalogues, while the errors are computed as their standard deviations. The true wp(rp) computed for the whole parent sample (open circles) is compared to that measured from the observed sample (filled circles). Upper panel: without corrections (2) and (3). Lower panel: after our full correction scheme was applied.

In the galaxy surveys, there exists a number of observational biases which may influence the clustering measurements. These biases are related to the construction of the survey. Both the parent photometric catalogue, from which target galaxies are selected, and the final spectroscopic catalogue may be the source of these biases. Below, we list the most important survey features that cause the biases. Next, we describe the corrections applied to minimize their effects. Since the VUDS survey is performed on the same multi-slit spectrograph VIMOS as the VVDS, zCOSMOS and VIPERS surveys, the correction scheme adopted here is based on the same methodology and the ones used for the other VIMOS-based surveys in the past. In particular, the main corrections used here were proposed and fully described by Pollo et al. (2005) and de la Torre et al. (2011a).

The main observational biases in the VUDS data can be listed as follows.

• 1.

  As shown in Fig. 1, some of the areas are excisedfrom the observations, due to the VIMOS layout (see A.2). Thefield of view of this instrument consists of four quadrantsseparated by 2′ gaps. No galaxies are observed in these gapsbetween quadrants, which influences the pair counts.

• 2.

  Obviously, not all galaxies from the photometric target candidate sample can be spectroscopically observed. Each galaxy spectrum occupies a certain area on the CCD detector, calculated from the spatial extent of the slit and the length of the spectrum due to spectral dispersion, which imposes geometrical constraints in target selection as galaxies too close to each other cannot be targeted simultaneously. It means that the galaxies from the parent photometric sample must be chosen is specific way to effectively allocate as many as possible spectral slits for each observation. In the case of the VIMOS observations, the slit allocation is performed automatically by the Super-SPOC code (Bottini et al. 2005). As a result, the spectroscopically observed sample is not a random representation of the parent photometric sample, and the introduced bias is relatively complex, especially on small scales.

These biases, as well as a few other, less important, effects are described in detail by Pollo et al. (2005).

However, these biases can be minimized by a combination of corrections, which were presented by Pollo et al. (2005) and de la Torre et al. (2011a), and which we found the best working for the VUDS data:

• 1.

  Random sample construction. The first, most basic part of thecorrection scheme is the appropriate construction of the randomcatalogue which needs to be geometrically identical (with theexception of the small-scale non-linear biases introduced by theSSPOC) with the spectroscopic sample. Generating randomobjects in this catalogue, we take into account the shape of thesingle pointings and quadrants, which are related to the shape ofthe VIMOS spectrograph. Additionally, we exclude the regionsremoved from the parent photometric sample (e.g. due to thepresence of a bright star), by applying the same photometric maskto the random sample. These first-order corrections reduce mostof the negative effects on the correlation function.

• 2.

  A global correction. In order to account for the missing pairs (due to the VIMOS limitations and the SSPOC strategy) we assign a global weight to each galaxy-galaxy pair. Assuming that the parent catalogue is free from angular incompleteness, we define a weighting function f(θ) as a ratio between the mean number of pairs in the parent photometric catalogue and the main number of pairs in the spectroscopic sample, as a function of angular separation (de la Torre et al. 2011a): $\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}{\mathit{w}}_{\mathrm{par}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{1}\mathrm{+}{\mathit{w}}_{\mathrm{spec}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\mathit{,}$(A.1)where w_parθ and w_spec(θ) are the angular correlation functions of the parent photometric and spectroscopic samples, respectively. Then, each pair from the spectroscopic sample separated by the angular distance θ is weighted by this ratio f(θ).

• 3.

  Small-scale corrections. As it comes out, even this strategy does not fully account for the small-scale angular effects of the SSPOC target selection strategy. In order to account for these local small-scale biases we additionally use the local weighting scheme which is also using the parent photometric catalogue as a reference. In this case, we count how many galaxies around the targeted galaxy we are missing due to the limited space for spectroscopic slits. Each targeted galaxy pair is then weighted proportionally to its representativeness of the surrounding density field. Similarly to Pollo et al. (2005) we found that the optimal size of the weighting area is the circle with a radius ~ 40′′.

The correcting scheme including all these ingredients was confirmed to be optimal after a series of tests performed using all 66 VUDS mock catalogues. During these tests we have been switching on and off different combinations of corrections listed in Sect. A1 (i.e. changing parameters controlling random catalogue properties, global correction scheme, weighting on small scales, etc.). All of these changes naturally influence the correlation function measured on different scales r_p. Hence, for each case we measured the deviations with respect to the real correlation function computed for the parent mock catalogues. As a result of this procedure we choose the set of corrections which – when combined together – minimized the differences between real and observed correlation functions at each separation r_p most efficiently. In Fig. A.1 we present the results of the correlation function measurements for the mock catalogues in our three redshift ranges. All the points correspond to the means of the measurements from 66 VUDS mock catalogues, while the errors are computed as the standard deviations between them. The three upper panels show the comparison between (1) the true (computed for whole mock parent sample) correlation function and (2) the observed correlation function computed only with the most obvious correction number (1), i.e. the appropriate geometrically cut random sample. The lower panels present the comparison of the same true correlation function and (3) the observed correlation function when the full correcting scheme described above was applied.

It is clearly seen that the introduced corrections significantly improve the retrieved signal from the observed sample at all the spatial scales, and in all three redshift bins considered. However, it can be seen that on small scales (r_p< 3 h^-1 Mpc), there is still a signal missing, even after introducing the full optimal correcting scheme (see lower panels of Fig. A.1). For these separations, our measurement is on average underestimated by 10% ± 3% with respect to the true value of w_p(r_p) computed from the parent sample. The possible influence of this effect on the final results is discussed in Sect. 3.3. At separations larger than r_p ~ 3 h^-1 Mpc the value of w_p(r_p) is also slightly underestimated, but only by about 1.6% ± 0.9%.

Our correction scheme provides a significant improvement in comparison to the average 28% ± 5% bias observed without applying it. However, the remaining systematic underestimation of the w_p(r_p) has to be taken into account during the correlation function measurements (see Sect. 3.3).

Fig. A.2 Layout of the VIMOS field of view. Picture from Bottini et al. (2005).

Appendix B: Impact of the reliability of redshift measurements on the correlation function accuracy

The precision and reliability of the galaxy correlation function measurement strongly depends on the reliability and accuracy of the measurement of their radial distance, i.e. on the quality of the redshift measurement in the survey. In this section we investigate the effect that a small fraction of the VUDS data with incorrect redshifts could have on galaxy clustering measurements. In this paper we use the VUDS galaxies with the best quality flags 2,3,4 and 9 (see Sect. 2.1 for the description of the quality flags), which in turn should result in a robust estimation of the correlation function. As the reliability of each flag has been independently estimated (Le Fèvre et al. 2015) the confidence level of the whole sample can be computed using the numbers of galaxies with different flags and their corresponding reliability estimates (95% at worst for flags 3+4, 75% for flag 2 and 80% for flag 9): $\mathit{P}\mathrm{=}\mathrm{(}{\mathit{N}}_{\mathit{f}\mathrm{4}}\mathrm{+}\mathrm{0.95}{\mathit{N}}_{\mathit{f}\mathrm{3}}\mathrm{+}\mathrm{0.75}{\mathit{N}}_{\mathit{f}\mathrm{2}}\mathrm{+}\mathrm{0.8}{\mathit{N}}_{\mathit{f}\mathrm{9}}\mathrm{)}\mathit{/}{\mathit{N}}_{\mathrm{g}}\mathit{,}$(B.1)where N_fi is the number of galaxies with a given i flag and N_g is the total number of galaxies within the sample. From this formula, for our selected VUDS sample the total confidence level is as high as ~ 86%, which is unprecedented at these redshifts and at this depth i_AB = 25.

How could the remaining ~ 14% fraction of galaxy with wrong redshifts possibly affect the correlation function measurement? To clarify this we performed a set of tests, as described below.

The presence of galaxies with wrongly measured redshift in the data sample can distort the correlation function w_p(r_p) at each r_p. In particular, we expect a decrease in the measured correlation length r₀ as the correlation amplitude will be diluted by a number of pairs of incorrect separations. In the mock catalogues, we can simulate these distortions by changing the spatial distribution of a corresponding fraction of mock galaxies. For our simulations we assume that galaxies with reliability flag 3 and 4 do not introduce any significant distortions to the correlation function measurements, thanks to the 95−100% confidence level. We also assume that the possible influence of flag 9 galaxies is negligible, due to their low number in the sample (see Table 1). Thus, the largest fraction of galaxies with possibly wrongly measured redshift is expected to come from flag 2 galaxies. Hence, we focused our attention on these galaxies.

Fig. B.1 Influence of the possible incorrect redshift measurements on the correlation function. Solid lines represent the projected two point correlation function wp(rp). Colours indicate the different levels of sample randomization (as labelled). The black solid line represents the true undistorted correlation function.

We measured how strongly the correlation function would be affected if 15%, 25%, 30% or 50% (highly unlikely case) of flag 2 galaxies would have wrongly measured redshifts. Within the VUDS galaxy sample selected for this study, galaxies with flag 2 represent 48% of the total sample (see Table 1). Therefore, in the mock catalogues we randomly selected a corresponding fraction of objects to play the role of flag 2 galaxies (i.e. ${\mathit{N}}_{\mathit{f}\mathrm{2}}^{\mathrm{mock}}\mathrm{=}\mathrm{0.48}{\mathit{N}}_{\mathrm{tot}}^{\mathrm{mock}}$). Then, in order to mimic their possibly wrong redshift measurements, we randomized the spatial coordinates (α and δ) of 15%, 25%, 30% and 50% of these objects, to obtain differently distorted samples. We then computed the projected correlation function w_p(r_p) and fitted it with a two parameter power-law.

The procedure was repeated for 10 independent VUDS mock catalogues (see Sect. 2.3) and the final estimation of w_p(r_p) was computed as the average. The results are presented in Fig. B.1. As expected, the amplitudes of the estimated projected correlation functions for the samples with wrongly measured redshifts are underestimated with respect to the unaffected one (black line), and the underestimation increases with the increasing fraction of wrongly estimated redshifts. However, the changes are not large except of the unrealistic 50% case: the correlation length for the fraction of 30% randomized galaxies differ form the real r₀ only by Δr₀ = 0.13 ± 0.04.

Based on these simple tests we conclude that the correlation function measurements performed in this study are robust, in spite of a possible contamination from wrong redshift measurements. Assuming the worst possible case, i.e. the lower limit for flag 2 galaxies confidence level being 70%, the resultant clustering strength r₀ will be underestimated only by Δr₀ = 0.13, which is ~ 3.3% of the measured amplitude. Therefore, the upper error bar of the measured r0 has been rescaled by a corresponding value, to fully account for the uncertainties introduced by possible wrongly measured redshifts in the data.

Appendix C: Projected correlation function measurements

For every galaxy subsample used in this work we present the tables of the projected correlation function measurements w_p(r_p) with associated 1σ errors at different separations r_p (in units h^-1 Mpc).

Acknowledgments

We thank Martin Kilbinger for assistance with implementing our HOD model into the CosmoPMC code. This work is supported by funding from the European Research Council Advanced Grant ERC-2010-AdG-268107-EARLY and by INAF Grants PRIN 2010, PRIN 2012 and PICS 2013. A.C., O.C., M.T. and V.S. acknowledge the grant MIUR PRIN 2010–2011. This work is based on data products made available at the CESAM data centre, Laboratoire d’Astrophysique de Marseille. This work partly uses 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. This work is supported by the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02). A.P. is supported by the Polish National Science Centre grant UMO-2012/07/B/ST9/04425 and the Polish-Swiss Astro Project. Research conducted within the scope of the HECOLS International Associated Laboratory, supported in part by the Polish NCN grant DEC-2013/08/M/ST9/00664.

References

All Tables

All Figures

	Fig. 1 Distribution of galaxies with spectroscopic redshifts 2 <z< 5 in two independent VUDS fields: COSMOS on the left panel and VVDS-02h on right panel. The blue squares indicates VIMOS pointing centres.
In the text

	Fig. 2 Redshift distribution of the used VUDS galaxy sample for redshift range 2.0 <z< 5.0. The filled grey histogram represents total sample of reliable galaxies (with zflag = 2−4,9, ~ 80% reliability), while the red and blue histograms represent the contribution from COSMOS and VVDS2h fields respectively.
In the text

	Fig. 3 Projected two-point correlation function wp(rp) in each individual VUDS field for the redshift range 2 <z< 5. Blue and red lines correspond to the VVDS and COSMOS field, respectively. Black points and line indicate the combined correlation function from the measurements preformed on galaxies from both fields simultaneously.
In the text

Fig. 4 Left panel: projected two-point correlation function wp(rp). The symbols and error bars denote measurements of the composite correlation function from the VUDS survey. Right panel: associated best-fit power-law parameters r0 and γ along with error contours for the general galaxy population separated in three redshift ranges. Top panels present a measurement for a wide redshift range 2 < z < 5 sample, while lower panels present the correlation function measurements for two independent low and high redshift galaxy samples.

In the text

Fig. 5 Correlation length r0 as a function of redshift from various surveys of different types of objects. Stars indicate measurements obtained from this work in three redshift ranges (magenta for a wide redshift range sample, and black for a low and high redshift samples). Different colours indicate different galaxy populations targeted by using various techniques, as indicated in the upper right corner. Open symbols indicate measurements based on photometric data, while filled symbols are for measurements from spectroscopic data. Blue: LBG galaxies (open squares – Foucaud et al. 2003; open circles – Ouchi et al. 2004; open triangles – Adelberger et al. 2005; open reversed triangles – Kashikawa et al. 2006; open diamonds – Savoy et al. 2011; filled diamonds – Bielby et al. 2013; open pentagon – Barone-Nugent et al. 2014). Orange: BzK galaxies (open circles – Blanc et al. 2008; open triangles – Hartley et al. 2010; open reversed triangle – McCracken et al. 2010; open diamonds – Lin et al. 2012). Green: galaxy samples from surveys limited in luminosity (filled squares – Norberg et al. 2002; open circles – Coil et al. 2006; filled triangles – Le Fèvre et al. 2005; filled reversed triangles – Pollo et al. 2006, filled diamonds – Zehavi et al. 2011; filled pentagons – Marulli et al. 2013; crosses – Skibba et al. 2009). Red: EROs or massive red galaxies (open squares – Daddi et al. 2003; filled squares – Zehavi et al. 2005a; open circles – Brown et al. 2008). The five solid curves show the correlation length of dark haloes with different minimum masses, as labelled.

In the text

Fig. 6 Left panels: correlation function for three redshift subsamples as in Fig. 4. Top panel represents measurements for a wide redshift range 2 <z < 5 sample, while lower panels present correlation function measurements for two independent low and high redshift galaxy samples. The solid lines present the best-fit halo occupation model. Right panel: the evolution of the halo occupation for the three redshift ranges.

In the text

Fig. 7 Evolution of the number-weighted average host halo mass given by Eq. (10) for the three redshift ranges analysed in this study. The red filled square indicate mass estimations from VUDS for the wide range redshift sample, while red filled circles for low and high redshift ones. Black and grey symbols represent the results of previous work based on spectroscopic and photometric surveys respectively. The solid black lines indicate how a host halo of a given mass M0 at z = 0 evolves with redshift, according to the model given by van den Bosch (2002). The solid red line represents the halo mass evolution derived using Eq. (21), with the HOD parameters obtained from the best-fit HOD model at a redshift z ~ 3. The dashed red line is using the HOD best-fit parameters for z ~ 2.5. VUDS galaxies with a typical L⋆ luminosity are likely to evolve into galaxies with a luminosity >L⋆ today.

In the text

Fig. 8 Large-scale linear galaxy bias retrieved by using two methods. Red filled square indicate measurements of ${\mathit{b}}_{\mathrm{L}}^{\mathrm{HOD}}$ for the wide range redshift sample, while the red filled circles for the low and high redshift samples. The open red square and circles represent ${\mathit{b}}_{\mathrm{L}}^{\mathrm{PL}}$ in three redshift ranges respectivley. Green and blue point represents the results from previous works, retrieved by using HOD (green circles – Blake et al. 2008; green diamonds – Wake et al. 2008; green reversed triangles – Ross & Brunner 2009; green polygons – Sawangwit et al. 2011; green squares – Kim et al. 2014; green triangles – Bielby et al. 2013) and power-law frameworks (blue squares – Marinoni et al. 2005; blue triangles – Adelberger et al. 2005; blue reversed triangles – Geach et al. 2012; blue diamonds – Barone-Nugent et al. 2014) We show 3 different models for the evolution of bias: a passive evolution model (red dashed line) a merging model (red dot-dashed line), a star-forming model (dotted line), as described in the Sect. 5.3.

In the text

Fig. A.1 Projected two point correlation function wp(rp) measured for 66 VUDS mock catalogs in three redshift ranges. The points correspond to the mean measurement of all 66 mock catalogues, while the errors are computed as their standard deviations. The true wp(rp) computed for the whole parent sample (open circles) is compared to that measured from the observed sample (filled circles). Upper panel: without corrections (2) and (3). Lower panel: after our full correction scheme was applied.

In the text

	Fig. A.2 Layout of the VIMOS field of view. Picture from Bottini et al. (2005).
In the text

	Fig. B.1 Influence of the possible incorrect redshift measurements on the correlation function. Solid lines represent the projected two point correlation function wp(rp). Colours indicate the different levels of sample randomization (as labelled). The black solid line represents the true undistorted correlation function.
In the text