{'fraction_non_alphanumeric': 0.03535771585068007, 'fraction_numerical': 0.02251244162400237, 'mean_word_length': 4.293913984311604, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strong Clustering of Lyman Break Galaxies around Luminous Quasars at z~4

In the standard picture of structure formation, the first massive galaxies are expected to form at the highest peaks of the density field, which constitute the cores of massive proto-clusters. Luminous quasars (QSOs) at z~4 are the most strongly clustered population known, and should thus reside in massive dark matter halos surrounded by large overdensities of galaxies, implying a strong QSO-galaxy cross-correlation function. We observed six z~4 QSO fields with VLT/FORS exploiting a novel set of narrow band filters custom designed to select Lyman Break Galaxies (LBGs) in a thin redshift slice of Delta_z~0.3, mitigating the projection effects that have limited the sensitivity of previous searches for galaxies around z>~4 QSOs. We find that LBGs are strongly clustered around QSOs, and present the first measurement of the QSO-LBG cross-correlation function at z~4, on scales of 0.1<~R<~9 Mpc/h (comoving). Assuming a power law form for the cross-correlation function xi=(r/r0_QG)^gamma, we measure r0_QG=8.83^{+1.39}_{-1.51} Mpc/h for a fixed slope of gamma=2.0. This result is in agreement with the expected cross-correlation length deduced from measurements of the QSO and LBG auto-correlation function, and assuming a linear bias model. We also measure a strong auto-correlation of LBGs in our QSO fields finding r0_GG=21.59^{+1.72}_{-1.69} Mpc/h for a fixed slope of gamma=1.5, which is ~4 times larger than the LBG auto-correlation length in random fields, providing further evidence that QSOs reside in overdensities of LBGs. Our results qualitatively support a picture where luminous QSOs inhabit exceptionally massive (M_halo>10^12 M_sun) dark matter halos at z~4.


INTRODUCTION
Our understanding of structure formation suggests that small inhomogeneities in the density field shortly after the Big Bang grew over cosmic time via gravitational instability (e.g. Dodelson 2003;Padmanabhan 2006;Schneider 2015) into massive dark matter halos at z = 0. As clusters of galaxies are the most massive, gravitationally bound structures in the Universe, they must have formed from the highest density peaks at early times. This make them ideal laboratories for studying the formation and evolution of cosmic structure.
Because of the small areas of sky surveyed at highredshift, and the low comoving number density ∼ 10 −7 Mpc −3 of local clusters (Gioia et al. 2001;Vikhlinin et al. 2009), the evolutionary link between these lowredshift clusters and high-redshift galaxies has been challenging to make. The progenitors of clusters are extremely difficult to identify when the density contrast between the forming cluster and its surroundings is small. For this reason, a commonly adopted approach is to cjgarci1@uc.cl 1 Based on observations collected at the European Organization for Astronomical Research in the Southern Hemisphere, Chile. Data obtained from the ESO Archive, Normal program, visitor mode. Program ID: 079.A-0644. 2 We dedicate this work to the memory of Josef Fried, who originally obtained and analyzed the data on which this work is based.
5 Department of Physics, University of California, Santa Barbara, CA 93106, USA search for these so-called proto-clusters around known massive galaxies at high redshift.
One very fruitful technique to find high-redshift protoclusters has been to use the presence of an active super massive black hole (BH) as a signpost for a massive galaxy and hence massive dark matter halo in the distant Universe (e.g. Venemans et al. 2007;Kashikawa et al. 2007;Overzier et al. 2008;Morselli et al. 2014). This technique is motivated by several considerations. First, the masses of supermassive BHs (M BH ) are known to tightly correlate with the bulge mass of their host galaxy (Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000), and possibly with the masses of their host dark halos (M halo ) (Ferrarese 2002, but see Kormendy & Bender 2011). Intriguingly, the most luminous quasars (QSOs) at z > 3 have M BH ∼ 1 − 6 × 10 9 M (Shen et al. 2011), comparable to the most massive known local BHs. If the present day M BH −M halo relation holds at early times, such BHs should reside in exceptionally massive halos. Second, some studies have suggested that the nuclear activity in active galactic nuclei (AGN) is triggered by processes related to the environment where they reside. For example, galaxy mergers could trigger the AGN activity (Bahcall et al. 1997;Wyithe & Loeb 2002;Hennawi et al. 2015), and galaxy mergers occur preferentially in dense environments (Lacey & Cole 1993). This would imply that the existence of an AGN requires a dense environment around it. Finally, another line of evidence that QSOs trace the rarest environments at high redshift arises from their extremely strong clustering. Indeed, Shen et al. (2007) determined that QSOs at z > 3.5 have a comoving auto-correlation length of r 0 = 24.3 h −1 Mpc (for a fixed correlation function slope of γ = 2.0), making them the most strongly clustered population in the universe, and demanding that they reside in the most massive M halo > 10 12 M dark matter halos at this epoch. Additionally, the Shen et al. (2007) correlation function, agrees with that required to explain the abundance of binary QSOs at z > 3.5 Shen et al. 2010), indicating that overdense structures around QSOs extend down to scales as small as 100 h −1 kpc. Since in hierarchical clustering models, QSOs and galaxies trace the same underlying dark matter density distribution, the generic prediction is that galaxies should be very strongly clustered around QSOs at z 3.5. Observationally this should be reflected as a strong QSO-galaxy cross-correlation function.
The QSO-galaxy cross-correlation function has been measured at z < 4 in the past. At z 1 it is found to be in good agreement with the auto-correlation of galaxies and QSOs, and it has been shown to be independent of the QSO luminosity, and weakly dependent on redshift (e.g. Padmanabhan et al. 2009;Coil et al. 2007). Adelberger & Steidel (2005) measured the AGN-galaxy cross-correlation function at higher redshifts (2 z 3) finding a cross-correlation length of r 0 ∼ 5 h −1 Mpc for a slope of γ = 1.6 which is similar to the auto-correlation of Lyman Break Galaxies (LBGs) at z ∼ 3 . They also claim an independence of the cross-correlation length with the AGN luminosity, implying that both faint and bright AGNs should be found in halos with similar masses. The highest redshift measurement of QSO environments is the work of Trainor & Steidel (2012), who quantified the clustering of LBGs around 15 hyper-luminous QSOs at z = 2.7. They find a QSO-LBG cross-correlation length of r 0 = 7.3 ± 1.3 h −1 Mpc for a fixed slope of γ = 1.5 and claim that this measurement is in agreement with the Adelberger & Steidel (2005) results. Additionally, they compute a halo mass for those QSOs of log(M halo /M ) = 12.3 ± 0.5, which is in agreement with the halos masses inferred for fainter QSOs at the same redshift .
Theoretical considerations suggest that high-redshift QSOs live in massive dark matter halos, but not necessarily the most massive ones (Fanidakis et al. 2013). However, a high signal to noise clustering analysis is necessary to confirm this hypothesis.
In addition to these statistical clustering analyses, many studies of individual AGN environments have been conducted. The population of AGNs whose environments have been studied most intensively are the highredshift radio galaxies (HzRGs) at z ∼ 2 − 4, which have been shown to often reside in proto-cluster environments (e.g. Venemans et al. 2007;Intema et al. 2006;Overzier et al. 2008;Hennawi et al. 2015). At higher redshifts the environments of other classes of AGN, such as opticallyselected QSOs, are currently less well constrained. Most previous work focuses on searching for galaxies around the most distant z 5 QSOs, and these results paint a diverse and rather confusing picture: Stiavelli et al. (2005), Zheng et al. (2006), Kashikawa et al. (2007), Utsumi et al. (2010), and Morselli et al. (2014) find a quite strong enhancement of galaxies compared to control fields around z ∼ 5 − 6 QSOs, whereas Willott et al. (2005) Bañados et al. (2013), Simpson et al. (2014) and Mazzucchelli et al. (2016), find no significant excess of galaxies around QSOs at z ∼ 6 − 7. Kim et al. (2009) studied five QSO fields at z ∼ 6 and reported a mix of overdensities and underdensities, and Husband et al. (2013) find galaxy overdensities in z ∼ 5 QSOs environments, but they note that even some randomly chosen patches of sky without AGN signposts ('blank fields') at the same redshift contain similar galaxy overdensities. Indeed, surveys of a few deg 2 for z ∼ 6 LBGs or LAEs have identified comparable overdensities in blank field pointings (e.g. Ouchi et al. 2005;Ota et al. 2008;Toshikawa et al. 2012). These mixed results at z 5 do not yet provide compelling evidence QSOs inhabit massive dark matter halos at the highest redshifts, and more work is clearly required.
One complication of these studies is that the majority of them are focused on dropout selection, which selects galaxies over a broad redshift range of ∆z ∼ 1 (e.g. Ouchi et al. 2004a), corresponding to ∼ 520 h −1 cMpc at z = 4. A large part of such a volume is unassociated with the QSO, which introduces projection effects that dilute the overdensity around the QSO making it much more difficult to detect. Furthermore, most work at the highest redshifts have focused their searches around a handful of individual QSOs, and given the poor statistics and large cosmic variance (which is typically not taken into account), this could preclude the detection of an overdensity.
In this paper we study the environs of QSOs at z ∼ 4. There are several advantages to working at this redshift. First, it is the highest redshift at which auto-correlation measurements exist for QSOs (Shen et al. 2007), establishing that they reside in massive dark matter halos. Second, the luminosity function and clustering properties of z ∼ 4 galaxies are also well known (e.g. Shen et al. 2007;Ouchi et al. 2004aOuchi et al. , 2008. The well-measured luminosity function allows us to accurately determine the background number density, essential for a robust clustering analysis. Furthermore, the fact that the autocorrelation of QSOs and galaxies are both known, gives us an idea of what the cross-correlation should be on large scales where linear bias models apply. In practical terms, redshift z ∼ 4 also represents a compromise since the dark matter halos hosting QSOs are still expected to be massive (Shen et al. 2007), while at the same time the characteristic galaxy luminosity L * can be imaged with much shorter exposure times than galaxies at z 5, allowing us to observe a larger statistical sample of QSO fields. Note that at z ∼ 4 the universe was only ∼ 1.5 Gyr old, and only 0.5 Gyr has elapsed since the end of reionization. Thus, our QSO targets are definitely young objects residing in large scale structures that are still forming.
Here we present VLT/FORS imaging of six z ∼ 4 luminous QSOs fields. Using a novel narrow band (NB) filter technique designed to select LBGs in a narrow redshift range (∆z ∼ 0.3) around QSOs. This minimizes the lineof-sight contamination, dramatically reducing the projection effects which are inherent in broad-band selection. We measure the QSO-LBG cross-correlation function at z ∼ 4 for the first time, to determine whether luminous QSOs at z ∼ 4 are surrounded by overdensities of galaxies. The large sample of QSOs studied allows us to beat down the noise from limited numbers of galaxies and cos-mic variance.
The outline of this paper is as follows. In section § 2 we describe the QSO target selection, we explain the novel NB imaging technique used to select LBGs, and we give details of the imaging observations, data reduction, and photometry. We present the color criteria used to select LBGs and compute the redshift selection function of the sample in section § 3. The measurement of the QSO-LBG cross-correlation function and LBG auto-correlation in QSO fields are presented in section § 4, where we also estimate the power law correlation function ξ(r) = (r/ r 0 ) −γ parameters r 0 and γ. We test the robustness of our results in section § 5, and summarize and conclude in section § 6.
Throughout this paper magnitudes are given in the AB system (Oke 1974;Fukugita et al. 1995) and we adopt a cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.26 and Ω Λ = 0.74 which is consistent with the nine-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations (Hinshaw et al. 2013). Comoving and proper Mpc are denoted as "cMpc" and "pMpc", respectively.

OBSERVATIONS AND DATA REDUCTION
The dataset presented in this section was obtained from the ESO Archive (Program ID: 079.A-0644, P.I: Rix). This program was designed to search for LBGs in z ∼ 4 QSOs environments using a novel NB filter technique. The aim was to test whether QSOs with the most massive BHs at z ∼ 4 live in the most massive dark matter halos. In this section we summarize the strategy used to select the targeted QSOs, we explain the NB technique used to select LBGs, and we provide details of the imaging observations, reduction process, and photometry.

QSO Target Selection
The PI of this program designed a custom set of filters (see § 2.2 for details) to search for LBGs in QSO environments. Using experiments with mock-catalogs, they showed that this filter set allowed one to select galaxies with z = 3.78 ± 0.08. Given this small redshift interval, and with the goal of stacking the galaxy number counts from several QSO fields, the QSO targets were selected to span a narrow redshift range of ∆z = 0.04, centered at z = 3.78.
Taking advantage of the large sample of QSOs from the Sloan Digital Sky Survey (SDSS; York et al. 2000), they first selected all QSOs in this redshift range. Given the goal of studying the most massive dark matter halos at z ∼ 4, believed to be correlated with the most massive BHs, only QSOs with M BH 10 9 M were selected. As is typical, the M BH was estimated from the emission line widths and continuum luminosities using the so-called single-epoch reverberation mapping technique (Vestergaard 2002). One of the targeted QSOs was not selected from SDSS, but it was added to the sample because it belongs to the redshift and M BH range of interest. The final sample is comprised of six bright QSOs with i < 20.2 mag.
We verified that none of the QSOs had a detected radio emission counterpart at 20cm by checking the Faint Images of the Radio Sky at Twenty-centimeters (FIRST Becker et al. 1995) catalog, since it is known that radio emission could strongly affect the galaxy clustering properties in AGN environments (e.g. Venemans et al. 2007;Shen et al. 2009). A summary of the QSO properties are listed in Table 1, where we show more recent M BH estimates taken from Shen et al. (2011).

A Novel Method to Select LBGs
The traditional Lyman break technique used to select high-redshift galaxies relies on the detection of the 912Å flux break (the so-called Lyman limit break) observed in galaxies due to the absorption of photons with λ < 912Å by neutral hydrogen in their interstellar and circumgalactic media. For this selection method, two bands are typically used bracketing the break, one located at λ < 912(1+z)Å, and the other at λ > 912(1+z)Å, where z is the redshift of the galaxies in question. Given this configuration, a non-detection is expected in the band blueward of the break, whereas a clear detection is expected redward of it, such that a very red color will be measured. Additionally, a third band is added at longer wavelengths in order to eliminate possible contaminants. This method was originally explored using the UGR filter system to detect galaxies at z ∼ 3 (Steidel et al. 1995(Steidel et al. , 1996, however, it was subsequently generalized to higher redshift (z ∼ 4 − 5) by using a filter set shifted to longer wavelengths (Steidel et al. 1999;Ouchi et al. 2004a).
At higher redshifts (z 4), a second break in galaxy spectra becomes important. The Lyα opacity of the intergalactic medium (IGM) rapidly increases with redshift, such that a large fraction of photons emitted by galaxies with λ < 1216Å are absorbed by neutral hydrogen. This implies a break at λ = 1216Å (the so-called Lyman alpha break), which can be used to select galaxies analogous to the traditional Lyman break technique described above. This Lyα break detection technique has been used to successfully identify galaxies and QSOs at z 6 (Fan et al. 2000;Bouwens et al. 2007Bouwens et al. , 2010Oesch et al. 2010;Bañados et al. 2016).
In order to achieve our goal of selecting galaxies physically associated with high-redshift QSOs, we need to select LBGs within a narrow redshift range centered on the QSO. However, the Lyman break method (using either the Lyman limit or Lyα breaks) efficiently selects LBGs in a broad redshift slice of ∆z ∼ 1 (e.g. Ouchi et al. 2004a;Bouwens et al. 2007Bouwens et al. , 2010, corresponding to ∼ 520 h −1 cMpc at z = 4. For such a broad redshift range, the overdensity signal around the QSO will be significantly diluted by the projection of galaxies at much larger distances, hundreds of comoving Mpc away. In order to address this problem, the PI proposed a novel selection technique analogous to the Lyα break method, but with the difference that the selection of LBGs is performed using two NB filters located very close to each other, instead of using broad bands. These filter curves are compared to those used for traditional LBG selection in Fig. 1. The advantage of using NB filters is that they allow one to select LBGs in a much narrower redshift range of ∆z ∼ 0.3 (∼ 167 cMpc at z = 3.78) (see section § 3.4), which is ∼ 3.3 times smaller than the redshift range covered when broad bands are used, allowing one to minimize line-of-sight projections from physically unassociated galaxies. This method has never been used before to select LBGs, and the filters used to perform the observations  shown on a LBG simulated spectrum at z = 3.78 (see section § 3.1 for the simulated spectra details). The NBs were designed specially for this program, in order to identify LBGs at z ∼ 3.78 by detecting the Lyman alpha break. This filter configuration selects galaxies in a quiet narrow redshift slice of ∆z ∼ 0.3. Lower panel: Example of a filter set used to identify galaxies with the standard Lyman break technique which is based in the detection of the Lyman limit break. The filter curves shown are those used by Ouchi et al. (2004a) to find LBGs at z ∼ 4 over a redshift slice of ∆z ∼ 1.0 were custom designed to select LBGs at z ∼ 3.78 centered on the redshift of our six QSO targets. The two NB filters used in this study are NB 571 (λ eff = 5657Å, FWHM = 187Å), and NB 596 (λ eff = 5947Å, FWHM = 116Å), which were designed to have a gap between them to exclude the Lyα emission line at z = 3.78. Then the galaxy selection is not influenced by the Lyα linestrength , but rather is sensitive to the Lyα break. Additionally data was also collected in the broad band filter r GUNN (λ eff = 6490Å) to help remove low-redshift interlopers.
Each QSO field was observed in the three filters shown in Fig. 1. The total exposure time for the filters was 8000s, 4000s, and 1800s for NB 571 , NB 596 , and r GUNN respectively. Observations were acquired in shorter individual dithered exposures, in order to fill the gap between the CCDs and to facilitate the data reduction process (cosmic ray and bad pixel rejection, building a superflat, etc). A spectrophotometric standard star was observed only in the second and third night. The typical seeing during the three nights was 0.6 -0.8 arcsec.
Science images were reduced using standard IRAF 6 tasks and our own custom codes written in the Interactive Data Language (IDL). The reduction process included bias subtraction and flat fielding. As our images exhibited illumination patterns, we performed the flat fielding with superflat images, created using the unregistered science frames. For that, we first masked all the objects out and then combined the science frames with an average sigma-clipping algorithm.
SExtractor (Bertin & Arnouts 1996) was used to create a source catalog for each individual image and then SCAMP (Bertin 2006) was used to compute an astrometric solution, using the SDSS-DR7 r-band star catalogs as the astrometric reference. Finally, the individual images were sky-subtracted, re-sampled, and median-combined using SWarp (Bertin et al. 2002), and then the noisy edges of the combined images were trimmed.
For the flux calibration, we only had observations of the spectrophotometric standard star SA109-949 at the beginning of the last two nights. The tabulated spectrum of this star has a coarse sampling of 25Å (Stone 1996) which is not suitable when NB filters are used. For the first night, spectrophotometric standard stars were not observed, but we took advantage of two existing SDSS star spectra in one of the fields taken during that night. The coordinates of the stars with available SDSS spectra are RA star1 = 21.014, DEC star1 = 0.740872 and RA star2 = 21.057, DEC star2 = 0.686577 and the median signal-to-noise ratio per angstrom of their spectra at the wavelengths of interest is 13.3 and 8.5 respectively.
The flux calibration process was as follows. For the first night calibration we convolved the SDSS star spectra with the three filters curves in order to obtain standard magnitudes. These magnitudes were compared with the stars instrumental magnitude (obtained using the MAG AUTO of SExtractor on the combined science images) to obtain the zero-points (ZPs) for each filter. A mean final ZP was computed from the two stars and the typical error for this ZP measurement was ∼ 0.08 mag. For the second and third night calibration, we used the spectrum of the observed spectrophotometric star to convolve it only with the broad-band filter curve to obtain the r GUNN ZP. The error in this computation was ∼ 0.02 mag. After that, the differential ZPs from the first night were used to determine the NB zero-points for the second and third nights for which we obtained a typical error of ∼ 0.11 mag.

Photometric Catalogs
Object detection and photometry were performed using SExtractor in dual mode, with the r GUNN image as the detection image. We set the parameters BACK SIZE and BACKPHOTO THICK such that the background was calculated in regions of 64 pixels in size and then recomputed locally in an annulus area of 24 pixels of thickness centered around the object. The parameters DETECT MINAREA and DETECT THRESH were set such that every group of at least five contiguous pixels having a value above 1.5σ (with σ the background RMS) was considered as an object.
In order to ensure an adequate color measurement we need to carry out photometry in the same object area for the three different filters. Therefore, we convolved our images with a Gaussian kernel to degrade its PSF to match it with the worst seeing image for each field. Then, the object magnitudes were estimated by the MAG APER parameter of SExtractor using a fixed aperture of 2 diameter. This magnitude is not necessarily the total magnitude of the object, but is used to compute the colors of galaxies. With this choice, if galaxies at z ∼ 4 are unresolved by the PSF, we are including the flux out to ∼ 3σ of the object's PSF (for a seeing of 0.8 ). This ensures that we measure the majority of the object's flux, as well as avoid contamination from other close sources. Magnitudes of objects not detected or detected with a signal-to-noise ratio (S/N) < 2 either in NB 571 or NB 596 were assigned the value of the corresponding 2σ limiting magnitude.
Here, the S/N of each object is defined as the ratio of counts in the 2 aperture, given by SExtractor, to the rms sky noise in the aperture. This rms sky noise is calculated using an IDL procedure which performs 2 aperture photometry in ∼ 5000 different random positions in the image (avoiding the locations of objects) to compute a robust measurement of the mean sky noise. The rms sky noise is calculated as the standard deviation of the distribution of mean values.
Magnitudes were corrected for extinction due to airmass using the atmospheric extinction curve for Cerro Paranal (Patat et al. 2011), and by galactic extinction calculated using the Schlegel et al. (1998) dust maps and extinction laws of Cardelli et al. (1989) with R V = 3.1. The error in the measured magnitude was computed by error propagation, with the object flux error given by the rms noise N in the aperture computed as we described above.
The mean 4σ limiting magnitude of the reduced images was of 26.06 for NB 571 , 25.53 for NB 596 and 25.82 for r GUNN for 2 diameter apertures. These limiting magnitudes are listed in Table 2 for each field.
For each field, we computed the completeness of the photometric catalogs for the image detection r GUNN . For that, we linearly fitted the logarithmic magnitude distribution in the magnitude range 21.0 < r GUNN < 24.5 where the photometric catalogs are assumed to be 100% complete. We extrapolated the linear fit to fainter magnitudes and we measured the completeness as a function of magnitude as the ratio of the histogram relative to that linear fit. We find that at our 4σ limiting magnitude the completeness is on average ∼ 12%.

LBG SELECTION AT Z=3.78
LBG candidates at z = 3.78 were selected using the Lyα break technique adapted to our custom filters, which target the Lyα break at λ rest−frame = (1 + z)1216Å. Our two NB filters were chosen to bracket this break, and thus we expect that LBGs at z = 3.78 will have red colors in NB 571 −NB 596 . But if we used only this color criteria, we could be including some low-redshift galaxy interlopers in the sample. In order to remove them, a third filter is used to give a measurement of the LBG continuum slope using the NB 596 − r GUNN color.
Since the filters used in this study are not standard, the color criteria to select LBGs is unknown. We also do not know what colors low-redshift galaxy contaminants have in this filter system. For this reason, we must explore how galaxies populate the color space in order to select a complete LBGs sample while avoiding low-redshift interlopers. Furthermore, in order to perform a LBGs clustering analysis in QSO fields we need to know the number density of LBGs expected at random locations in the universe. When a standard filter set is used (e.g. LBG selection using broad band filters), this number density can be computed directly from the LBG luminosity function measured from work using similar filters. However, in our case if we compute the number density from this LBG luminosity function, we have to correct this quantity to take into account the fact that our filter system is mapping a different survey volume and does not necessarily identify all of the LBGs selected by broad-band selection. Specifically, we need to a) determine what fraction of LBGs we are detecting at any redshift (i.e. the completeness) and b) determine the redshift range over which we are selecting LBGs (∆z). Both of these goals can be achieved by performing an accurate computation of the redshift selection function φ z (z), defined as the LBG completeness as a function of redshift.
In order to perform the optimal LBG selection and compute φ z (z), we conducted detailed simulations to model the distribution of LBG colors in the color-space. In this section we detail how the color modeling was performed, we study what contaminants could be affecting our LBGs selection, and we define a color criteria to select LBGs at z = 3.78. Finally, we present the redshift selection function providing the completeness as a func-tion of redshift for the sample.

LBG Color Modeling
We performed a Monte Carlo simulation of 1000 LBG spectra at each redshift, that were created to have different UV continuum slopes and Lyα equivalent widths (EW Lyα ), such that they reproduce the space of possible LBG spectra informed by our knowledge of LBG properties.
Each simulated rest-frame spectrum was created in the following way. As a starting point, we considered a template galaxy spectrum generated from Bruzual & Charlot (2003) population synthesis models 7 , corresponding to an instantaneous burst model with an age of 70Myr, a Chabrier (2003) IMF, and a metallicity of 0.4Z , as expected for LBGs at z ∼ 4 (Jones et al. 2012). We assumed a power law UV continuum for this template with amplitude A and a slope α BC , such that we modeled its flux as F BC (λ) = Aλ αBC . We fit this model to the template spectrum over the UV continuum range (here defined as 1300Å< λ <2000Å) by least-squares minimization to obtain the best fit A and α BC parameters.
First, we modified the UV slope of this template by multiplying its flux by λ α−αBC in order to obtain a spectrum with a power law UV continuum given by Aλ α . The new slope α was chosen as a value taken randomly from a Gaussian distribution with mean µ = −1.676 and σ = 0.39. These values are motivated by Bouwens et al. (2009), who presented the UV continuum slope distribution of LBGs at z ∼ 4 for samples selected in different magnitude ranges.
Second, we added a Gaussian Lyα line with rest-frame central wavelength λ Lyα = 1215.7Å, standard deviation σ Lyα and amplitude B which adjusts the intensity of the line. For all the simulated spectrum we used a fixed σ Lyα = 1Å which agrees with the σ Lyα of the composite spectrum of LBGs at z ∼ 4 (Jones et al. 2012). The B value was adjusted in order to model a Lyα line with a EW Lyα value drawn randomly from a distribution chosen to agree with observations of LBGs. The EW Lyα distribution was given by a Gaussian core plus a tail to large negative equivalent widths to represent strong line emitters. For the Gaussian core we adopted a mean µ = −25Å and standard deviation σ = 40Å (rest-frame), based on the measurements of Shapley et al. (2003), who studied the spectra of 811 LBGs at z ∼ 3. We thus assume that the Gaussian core of the LBG EW Lyα distribution does not evolve significantly from z ∼ 3 to z ∼ 4. For the tail representing strong line-emitters, we modified the Gaussian by adding an exponential function with rest-frame EW Lyα scale length of W 0 = −64Å, as presented in Ciardullo et al. (2012). In this way our model of line emission encompasses both LBG and LAE spectra. Fig. 2 shows the EW Lyα probability distribution function used to simulate our spectral models. The EW Lyα are defined as: where F Lyα is the flux of the Lyα line (with the continuum subtracted), which is given by a Gaussian with 7 Obtained from http://bruzual.org/ Normalized probability distribution function of EW Lyα used for the simulated spectra, where negative values correspond to emission lines. EW Lyα was chosen from a Gaussian distribution with rest-frame mean µ = −25Å and σ = 40Å  plus an exponential tail of high EW Lyα values with scale length of W 0 = −64Å (Ciardullo et al. 2012).
amplitude B, as we described above, and F cont is the flux of the continuum given by Aλ α . Note that we defined negative values of EW Lyα for emission lines and positive for absorption lines.
Once α and EW Lyα were chosen for a given simulated spectrum, we dust-attenuated it using the starburst reddening curve from Calzetti et al. (2000) and adopted a color excess value of E(B − V) = 0.16 according to the values estimated for LBGs at z ∼ 3 .
After the dust-attenuation is applied, we model the fact that only a small fraction of Lyman limit photons escape LBGs with an escape fraction parameter f λ<912 esc . Although this value is observationally poorly constrained, studies suggest it is in the range 0.04-0.14 (Fernández-Soto et al. 2003;Shapley et al. 2006;Ouchi et al. 2004a). We assumed a fixed value of f λ<912 esc = 0.05, and multiplied the spectrum at λ ≤ 912Å by this value. We also tested our results using different values of f λ<912 esc , finding that the colors of simulated galaxies are relatively insensitive to the exact value of f λ<912 esc used, because these wavelengths are subsequently significantly attenuated by the IGM transmission function (see below).
Finally, we redshifted each model spectrum to different redshifts on a grid with a grid spacing of 0.02 and ranging from z = 3.2 to z = 4.4. In the redshifting process we used the IGM transmission model T z (λ) for the corresponding redshift z from Worseck & Prochaska (2011) to attenuate the flux blueward of the Lyα line 8 . Note that in principle we should attenuate both the continuum blueward of the Lyα line and the line itself, however, the EW Lyα values used in this simulation are taken from the literature, which are observed values that are not corrected for IGM attenuation, such that this line emission is effectively already attenuated. In Fig. 3 we show some examples of our rest-frame simulated spectra, which have been normalized to have the same flux at λ = 1245Å.
At each redshift, we integrated the spectra against our three filter transmission curves to obtain the fluxes and then the LBG colors. In order to model the impact of noise, we added photometric errors to the simulated LBG photometry. To this end we first assigned an r GUNN magnitude to each simulated object by randomly drawing a value from the z ∼ 4 LBG luminosity function, integrated over the same magnitude range as our LBG sample (24.0 ≤ r GUNN ≤ 25.6 or 0.76 ≤ L/L * ≤ 3.5; see § 3.3) 9 . We also weighted the luminosity function by the completeness of the source detection at each apparent magnitude and for each field (computed in § 2.4), which takes into account the fact that the fraction of sources detected depends on their magnitude, such that the photometric catalog is complete for bright sources but less complete at the faint end. In this way the incompleteness of our photometry is also factored into our color modeling.
Based on the simulated LBG colors the chosen r GUNN value, we then determined the magnitude in the other two filters NB 571 and NB 596 for each spectrum in each redshift bin. In order to construct a noise model , we selected a galaxy sample from our photometric catalogs, and we computed the median magnitude error as a function of the magnitude for each filter (with the magnitude error computed as we explained in § 2.4). Finally, we assigned random Gaussian distributed magnitude errors using our median relations, and then added this noise to the model photometry which defined the final photometry of the simulated spectra. The colors for the 1000 simulated spectra at each redshift are shown in Fig. 4. We also computed the median of our 1000 rest-frame Monte Carlo spectra, redshifted it, and obtained the colors at each redshift to compute the median evolutionary track of LBG colors, shown as the black solid line in Fig. 4. Fig. 4 indicates that the median colors of LBGs at z = 3.78 are NB 571 −NB 596 = 1.05 and NB 596 −r GUNN = 9 Given that for each field we reached slightly different limiting magnitudes, we simulated the LBG photometry field by field according to their respective r GUNN limiting magnitude. This results in a slightly different redshift selection function for each field. 0.16. However, if we consider the intrinsic scatter in LBG properties (continuum slope and EW Lyα ) and photometric uncertainties, the z 3.78 LBGs (indicated by green points) span a wider color range with NB 571 − NB 596 0.5 and −0.6 NB 596 − r GUNN 0.8. In principle, we should select LBGs in this broad selection region to obtain a highly complete sample, however, we also need to take into account the colors of low-redshift galaxies in our filter system to define a final selection criteria. We perform this analysis in § 3.2, where we also test our LBG color modeling by reproducing the LBG evolutionary track presented in previous work using broad band LBG selection.

Low-Redshift Galaxy Colors
We use template galaxy spectra to develop a basic understanding of how low-redshift galaxies populate the color-color diagram in our new filters. We used a set of five commonly used templates for estimating photometric redshifts, such that they span the range of galaxy spectral energy distributions (SEDs). The templates are from the photo-z code EASY (Brammer et al. 2008), which are distilled from the PEGASE spectral synthesis models.
We redshifted these template spectra from z = 0 to z = 3, and integrated them over our filter transmission curves to generate their evolutionary track. Note that we need not attenuate these spectra by the IGM transmission function T z (λ), since our NB filters never cover rest-frame wavelengths lower than 1216Å for the low redshifts considered. In Fig. 5 we show the evolutionary tracks for different galaxy types together with the median LBG evolutionary track that we computed in § 3.1.
In order to test our Monte Carlo simulation as well as the evolutionary tracks for low-redshift galaxies, we have used our 1000 simulated spectra at each redshift to compute the median LBG evolutionary track in the standard BRi filter set used to select LBGs at z ∼ 4 (see Fig. 1) by Ouchi et al. (2004a). We also computed the evolutionary track of these low-redshift galaxies in the standard LBG filters in the same way as described above. These results are shown in Fig. 7, where we also overplot the selection region used by Ouchi et al. (2004a) to select z ∼ 4 LBGs. We find that the median LBG evolutionary track from our Monte Carlo model lies within the Ouchi et al. (2004a) selection region, and selects LBGs at z 3.5 as claimed. Note also that our LBG evolutionary track agrees well with the Ouchi et al. (2004a) evolutionary track (see Fig. 4 of their paper) determined from a much simpler model of LBG spectra and IGM transmission. In addition we see that the evolutionary tracks of lowredshift galaxies lie comfortably outside the BRi LBG selection region as claimed by Ouchi et al. (2004a).
However Fig. 5 shows that in our NB filter set, some of the low-redshift galaxies have similar colors as z = 3.78 LBGs, which suggests that our new filter configuration could make it challenging to select a sample of LBGs at z = 3.78 with high completeness and at the same time high purity. When we use NB filters the low-redshift galaxy colors are located in a wider region in the colorcolor plot in comparison with the location of the color locus of contaminants when broad band filters are used. We attribute this to sensitivity of the NB filters to features in the galaxy spectra such as emission or absorption lines. In the case of broad bands these features are diluted by averaging over large regions of spectra, but for NB the features result in large excursions in color with changing redshift, making the low-redshift galaxy locus extend over a larger region of color space that overlaps with the colors of z = 3.78 LBGs.
Given that LBG colors at z = 3.78 span the range NB 571 − NB 596 0.5 and −0.6 NB 596 − r GUNN 0.8 (see Fig. 4), there are several types of contaminants that could be affecting our LBG selection. Their colors are indicated by points on the respective low-redshift galaxy evolutionary tracks are labeled by with letters in Fig. 5, and some examples are shown in Fig. 6. The first type are red galaxies at z ∼ 0.45 having a large λ RF ∼ 4000Å Balmer break and strong Calcium H & K absorption. This break is located just between our two NBs, so they present red colors (point C in brown curve in Fig. 5). The second type of interlopers are star-forming galaxies at z ∼ 0.60 with strong [OII] 3727Å emission lines. If the NB 596 is located just over this line, and NB 571 over the continuum, we again detect red colors (point A on green curve in Fig. 5). The third type of interlopers are galaxies at z ∼ 1.04 with strong MgI and MgII absorption lines at λ RF = 2852Å, and λ RF = 2799Å respectively in com-bination with the λ RF ∼ 2900Å break. When the NB 571 is located over this absorption the filter NB 596 falls on the continuum, then, red colors are detected (points E and G on the magenta and brown curves, respectively, in Fig. 5). Other interlopers are galaxies with strong flux breaks redshifted just between our NB filters. One example are galaxies at z ∼ 1.23 with a large λ RF ∼ 2640Å break (points D and F on the magenta and brown curves, respectively, in Fig. 5) and galaxies at z ∼ 1.83 with a strong break at λ RF ∼ 2085Å (point B on the brown curve in Fig. 5).

Selection Region and LBG Sample
As we are interested in measuring the clustering properties of LBGs at z = 3.78, we need to select a sample with high completeness and purity. In order to avoid lowredshift contaminants, we were forced to choose a smaller selection region in the color-color diagram, which results in relatively low completeness, but it ensures that the sample is not highly contaminated.
First, we defined two vertical color cuts in Fig. 4, one to the left of the median LBG colors at z = 3.78 and one to the right. The first cut is meant to exclude LBGs located in the upper left region of the diagram, which  Figure 5. Evolutionary tracks of low-redshift galaxies redshifted from z = 0 to z = 3. We plot as brown, magenta, orange, blue, and red curves the evolutionary track of elliptical, Sa, Sb, Sc, and Irregular galaxies, respectively. We overplotted the track of LBGs computed as was explained in section § 3.1 as a black curve. Filled circles over the black curve indicate colors of LBGs from redshift 3.6 to 4.0, and the largest black point indicates the exact position of the color of LBGs at z = 3.78. Filled circles labeled with letters over the low-redshift galaxies evolutionary tracks are indicating the colors of some contaminants that could be affecting our selection: galaxies at z = 0.60 (A), z = 1.83 (B), z = 0.45 (C), z = 1.23 (D and F) and z = 1.04 (E and G). mostly corresponds to LBGs at z ∼ 3.9 with strong Lyα line emission. The second cut avoids LBGs at z > 3.9. A third color cut defines a lower limit for NB 571 − NB 596 , which ensures we are detecting the Lyα break, while at the same time avoiding LBGs at z 3.7. We used a diagonal color cut, to most effectively avoid the contamination of low-redshift galaxies (see Fig. 5), while at the same time including most of the LBGs at z = 3.78, thus  Figure 7. Same as Fig. 5 but using the filter system used by Ouchi et al. (2004a) to select LBGs at z ∼ 4 (broad band filters B, R, i). Filled circles over the black curve indicate colors of LBGs from redshift 3.2 to 3.9. Filled circles over the curve of the low-redshift galaxies indicate colors from redshift 0.0 to 3.0. The dashed line is indicating the region used Ouchi et al. (2004a) to select LBGs in their work.
maintaining the highest completeness possible.
We also tested several different color criteria to select LBGs. In section § 5.2 we will further discuss our color selection, contamination by low-redshift galaxies, and the impact that contamination can have on our clustering measurements. There we argue that the choice of color selection that we present here selects a reasonably complete LBG sample with high purity. Our final set of color cuts are shown in 4, and defined by the following relations: We selected LBGs based on our galaxy photometry, but required sources to have a S/N ≥ 4.0 in both the NB 596 and r GUNN filters, to ensure a solid detection of the LBG continuum. In order to reduce contamination by false detections, we only considered objects that have FLAGS= 0 in SExtractor, which excluded objects that were blended, saturated, truncated (too close to an image boundary), or affected by very bright neighboring objects. Bright stars in our images were masked in order to avoid spurious object detection due to contamination from their stellar flux. This procedure results in a set of masks indicating where we were able to detect galaxies, which we use later in our clustering analysis to compute the effective area of our survey.
We also imposed a lower limit on the magnitude in order to exclude bright low-redshift interloper galaxies from our selection. Thus we only considered objects with magnitudes fainter than r GUNN = 23.97, corresponding to LBGs with L ∼ 3.5 L * . We chose this value by computing the LBG luminosity function at z ∼ 4, and finding the bright end cut at which we would lose no more than 1% of the galaxies. In other words, 99% of the total number of LBGs have magnitudes between our bright end cut of r GUNN = 24.0 and the limiting magnitude r GUNN = 25.82 (mean limit magnitude at 4σ for a 2 Here the evolutionary track showed in Fig. 5 is plotted as redshift color-coded track according to the color bar. We have highlighted the selected LBGs as red points.
The magenta points indicate the color of each QSO in the filters. Arrows indicate lower limits for NB 571 − NB 596 color. These are cases in which the object was not detected in NB 571 filter at 2σ level and magnitude was replaced by the corresponding limit magnitude.
diameter aperture) of our images, which corresponds to L = 0.76 L * . In this way we can safely assume we are excluding only extremely rare bright LBGs. For the LBG luminosity function we used the Schechter parameters from Ouchi et al. (2004a) who studied the photometric properties based on a large sample of ∼ 2200 LBGs at z ∼ 4. The values used are φ * = 2.8 × 10 −3 h 3 70 Mpc −3 , M * 1700 = −20.6 mag, and α = −1.6. Given all of these selection criteria and the color cuts defined in eqn.
(2), we selected LBGs in each of our fields. We compute the total area of our survey by adding the effective area of each individual field, which is defined by subtracting the masked area from the total area of the image. The the total area of our survey is 232.7 arcmin 2 corresponding to an average area per image of 38.79 arcmin 2 (recall the FOV of FORS1 is 6.8 × 6.8 arcmin 2 or 46.24 arcmin 2 ). We show color-color diagrams of objects detected in all six of our fields in Fig. 8. We found a total of 44 LBGs (see Table 3) corresponding to a mean number density of 0.19 LBGs arcmin −2 . Image cutouts in our three filters for several of our selected LBGs are shown in Fig. 9. In Fig. 10 we show the spatial distribution of the LBGs relative to the QSO (red dot at zero) for our six fields. We also show the individual colorcolor diagrams and indicate the number of LBGs found in each individual field in Fig. 11. Note that the number of LBGs in the fields cannot be directly compared because each image has different limiting magnitude and different effective area (different reduced image size, masked region, etc). In Fig. 12 we show a false color image of the field around QSO SDSS J2301+0112 with the LBG candidate positions indicated.

Redshift Selection Function
We used the LBG color modeling machinery described in § 3.1 to compute the redshift selection function φ z (z) of our LBG color-selection criteria. At each redshift step, we redshift the 1000 rest-frame simulated LBG spectra into the observed frame, draw luminosities from the lumi- nosity function, compute magnitudes and colors, and add photometric errors. We then compute the completeness at each redshift by calculating the fraction of simulated LBGs that satisfy the selection criteria defined in § 3.3, namely: fulfill the color criteria in eqn. (2), and fulfill the magnitude constraints (given by the 4σ limiting magni-  tudes for NB 596 and r GUNN and by the bright end cut imposed for our selection, r GUNN > 24.0). Note that as the limiting magnitude of our fields are slightly different, we computed φ z (z) for each individual field, using their corresponding NB 596 and r GUNN limiting magnitudes. The final φ z (z) varied from field to field by a small amount, then we computed the median of φ z (z) over the six fields, which is shown in Fig. 13.
From Fig. 13 we see that our selection criteria recover a peak value of 26% of LBGs at z ∼ 3.78 over a small ∆z ∼ 0.04 interval. Our criteria also select 10% of LBGs over a wider redshift interval, ranging from z ∼ 3.65 to z ∼ 3.95 (∆z ∼ 0.3) which corresponds to ∼ 18, 800 km s −1 (∼ 167 h −1 cMpc) at z = 3.78. The NB technique selects LBGs over a much narrower redshift range compared to broad band LBG selection, which typically selects galaxies over a range ∆z ∼ 1.0 (e.g. Ouchi et al. 2004a;Bouwens et al. 2007Bouwens et al. , 2010, or ∼ 3.3 times larger than our selection. However, as low-redshift galaxies have similar colors as the LBGs in our filters, we have  The completeness per redshift bin was determined from 1000 simulated LBGs spectra with different EW Lyα , UV continuum slope, and r GUNN magnitudes. This is calculated by computing the fraction of the LBG simulated colors, per redshift bin, which was selected by the selection region. This defines the redshift selection function which is used for the clustering measurements in section § 4.
to adopt relatively conservative color cuts making our completeness relatively low.

CLUSTERING ANALYSIS
In this section we analyze the clustering of LBGs around our QSOs at z = 3.78. First, we present the measurement of the QSO-LBG cross-correlation function in § 4.1. We then estimate the correlation function parameters (r 0 , γ) assuming a power law form ξ(r) = (r/r 0 ) −γ . Our results are compared with theoretical expectations based on the auto-correlation of both LBGs and QSOs at z ∼ 4 assuming linear bias. In § 4.2 we present the LBG auto-correlation computed from our QSO fields, and finally in § 4.3 we compare our results with previous measurements.

The QSO-LBG Cross-Correlation Measurement
Following convention, we study clustering using the two-point correlation function ξ(r), which measures the excess probability over a random distribution, of finding an object at separation r from another random object, in a volume element dV (Peebles 1980). For the case of galaxies around QSOs this correlation function is defined by where ξ QG (r) is the QSO-galaxy cross-correlation function and n G is the mean number density of galaxies in the universe. Here r is real space comoving distance, which is however not the observable even when redshift information is available, as peculiar velocities induce redshift space distortions along the line-of-sight (Sargent & Turner 1977). Typically LBG clustering studies that lack redshifts measure the angular correlation function (e.g. Giavalisco et al. 1998;Ouchi et al. 2004b;Lee et al. 2006). Although we do not have redshifts of our LBGs, our NB selection technique selects LBGs over a narrow redshift interval ∆z 0.3 (see Fig. 13 in § 3.4) which allows us to measure clustering as a function of transverse comoving distance instead of angular distance (at z = 3.78, the angular diameter distance changes by just 3% over this redshift interval). Thus we write the real space separation r as r 2 = R 2 + Z 2 , where R is the transverse comoving distance between the QSO and the galaxy, and Z is the radial comoving distance between them approximately given by: where H(z) the Hubble constant evaluated at redshift z = 3.78, which we take to be a constant over the redshift interval considered (an approximation valid to 5%). We measure the volume-averaged projected crosscorrelation function between QSOs and LBGs χ(R min , R max ), which is a dimensionless quantity defined as the real space QSO-LBG cross-correlation function ξ QG (R, Z) integrated over a volume and then normalized by it (e.g. Hennawi et al. 2006Hennawi et al. , 2015: where V eff is a cylindrical volume defined by the radial bin [R min , R max ], the height ∆Z probed by our filter configuration, and modulated by the selection function of our survey. We measure χ(R min , R max ) in logarithmically spaced radial bins centered on the QSO for all fields using the estimator: where QG and QR are the number of QSO-LBG and QSO-random pairs in this cylindrical volume. The quantity QG is directly measured by counting the QSO-LBG pairs found in our images. The quantity QR is the expected random number of QSO-LBG pairs, which is given by: where V eff is the effective volume of the radial bin in question and n G (z, r lower GUNN < r GUNN < r limit GUNN ) is the mean number density of LBGs at redshift z in the magnitude range of our survey, which will be henceforth referred to as n G to simplify notation. Given that galaxy clustering measurements are normally performed in random locations of the universe, the mean number density measured from the survey is typically a good proxy for the mean number density of the universe, provided the survey volume is large enough. In such cases, the galaxy number counts for the random sample can be computed from the data itself, and one typically constructs random catalogs with a number density determined from the survey to estimate QR . However, in our case we are pointed towards a QSO situated in what is likely to be an overdensity. Therefore the mean number density of galaxies in our survey is not representative of the mean in random locations and we cannot follow the standard procedure for computing QR .
If we had observations of control fields (i.e not centered on QSOs) with our same filter configuration, then it would be possible to measure the background number density of LBGs directly and determine QR . Another alternative would be to measure this quantity directly from the outer parts of images, where the clustering becomes negligible, given a sufficiently large FOV instrument. Unfortunately, we do not have images of control fields, and the FOV of FORS1 is too small to provide a reliable measurement of the background. Thus our only alternative is to estimate QR from eqn. (7), where n G is calculated from the z ∼ 4 LBG luminosity function, and the effective volume V eff is determined from our Monte Carlo simulations of our selection function (see § 3.4) and the effective area covered by our survey. We provide further details of these computations in what follows.
To calculate n G we used the Schechter parameters from Ouchi et al. (2004a) z ∼ 4 LBG luminosity function. We integrated this luminosity function over the magnitude limits given by our LBG selection, and this magnitude integral was weighted by the photometric completeness fraction of our source detection following the same procedure described in § 3.1. Given that our fields all have slightly different limiting magnitudes and different source completeness, we compute the expected n G for each field. We assumed that n G is constant over the redshift ranged considered, which is a good approximation given the narrow redshift range ∆z 0.3 that we probe. The expected mean number density of LBGs in random fields, n G , computed for the magnitude range of each field is given in Table 4.
We define the effective volume of a radial bin as: where φ(R, Z) encodes the geometry of the survey, which can be separated into the radial R and the redshift (lineof-sight) Z selection function as φ(R, Z) = φ Z (Z)φ R (R). The redshift selection function of our survey φ z (z) was modeled in § 3.4 and we convert it to a redshift selection function in comoving units φ Z (Z) using eqn. (4). Then, we integrate it over the redshift range covered by our Monte Carlo modeling (corresponding to 3.2 ≤ z ≤ 4.4). The radial selection function φ R (R) is easily calculated using the detection masks for our images. We created catalogs with randomly distributed galaxies with number density n ran such that we had ∼ 100, 000 sources in the entire image. Then we computed φ R (R) in radial bins as the ratio between the number of randomly distributed galaxies and the expected number without masking n ran π(R 2 max − R 2 min ). The resulting φ R (R) then quantifies the fraction of the bin area where we could have detected LBGs. We computed the value of V eff for each radial bin in each field using eqn. (8). Summing the V eff over the radial bins, we obtain the total volume covered by each of our six fields V field , given in Table 4. We obtained that the total volume of our survey is 14,782 h −3 cMpc 3 .
To obtain a rough estimate of the LBG overdensity in our QSO fields, we calculated the expected number of random QSO-LBGs pairs QR field for each of our fields and compare to the number we find QG field . These results are tabulated in Table 4, where we also show the overdensity per field QG field / QR field . We see that five of our six fields exhibit an LBG overdensity of LBGs, while one appears to be underdense. Adding up the results for all six fields, we find that the random expectation is QR = 28.6 LBGs, whereas we detected a total of QG = 44 LBGs, giving an overall overdensity of 1.5, and indicating that our fields are on average overdense.
To explore the profile of this overdensity around QSOs, we computed QG and QR in bins of transverse distance for each of our six fields, and then summed them to determine the binned volume averaged cross-correlation function χ(R min , R max ) according to eqn. (6). These results are tabulated in Table 5 and plotted in Fig. 14. We estimate errors on χ(R min , R max ) assuming that shotnoise dominates the error budget, and use the one-sided Poisson confidence intervals for small number statistics from Gehrels (1986).
Given that the auto-correlation functions of both LBGs and QSOs at z ∼ 4 have been previously measured, we can compute the expected volume averaged QSO-LBG cross-correlation function χ(R min , R max ) assuming linear bias and compare it to our measurements. Since we are probing non-linear scales in our measurement where the linear bias assumption surely breaks down, the expected cross-correlation obtained in this manner is approximate, but nevertheless a useful reference. If we assume that both LBGs and QSOs trace the same underlying dark matter, and assume linear bias such that δ G = b G δ DM , and δ Q = b Q δ DM , then we can write ξ QG = ξ QQ ξ GG . Assuming a power law form ξ = (r/r 0 ) γ for the respective auto-correlations of QSOs and LBGs, and that  Shen et al. (2007) for z > 3.5 QSOs assuming a fixed γ = 1.8. Combining these implies r QG 0 = 9.6 h −1 cMpc for γ = 1.8. Plugging this power law LBG-QSO crosscorrelation function into eqn. (5) and integrating over the effective survey volume gives us the expected value of χ(R min , R max ), which is shown as a dashed line in Fig. 14. One sees that our QSO-LBG cross-correlation measurement is in reasonable agreement with the expected value of χ(R min , R max ) combining auto-correlation measurements and assuming linear bias. In § 4.1.1 we quantify this agreement by fitting our cross-correlation function.

Fitting the Cross-Correlation Function
Given the projected cross-correlation function measurement, we now determine the real-space crosscorrelation parameters r QG 0 and γ that best fit our data. To this end we use maximum likelihood estimator (MLE), and fit for the parameters which maximize the probability of the data we observe. Since we are dealing with a counting process with small number counts in each bin (see Table 5), we can assume that Poisson error dominates the error budget. Adopting the Poisson distribution for the counts in our cross-correlation function bins, we can write the likelihood of our data as Note. -(1) Field ID, (2) The mean number density of z ∼ 4 LBGs in units of (10 −3 h 3 cMpc −3 ), in the magnitude range of the survey r lower GU N N < rGUNN < r limit GUNN . Given that r limit GUNN and completeness in the source detection are different for each field, we obtain a number density slightly different for each one, (3) Total volume of the field in units of (h −3 cMpc 3 ), computed as V field = N bins i=1 V eff,i , (4) Total number of expected LBGs on the whole field computed as QR field = nGV field , (5) Total number of observed QSO-LBG pairs on the whole field, (6) Total overdensity per field, computed as QG field / QR field .

6018.75
Note. -We present the data for the volume-averaged projected cross-correlation function between QSOs and LBGs χ(Rmin, Rmax) shown in Fig. 14. This is measured in radial bins defined by Rmin and Rmax. QG is the observed number of QSO-LBG pairs per bin, and QR is the expected number of QSO-random pairs per bin, computed from eqn. (7). We also show the total volume of the bin added over the fields, computed as V eff,total = N fields i=0 where the product is over the N bins radial crosscorrelation function bins, x i is the number counts measured in the ith bin and λ i is the expected number counts in the ith bin for a given set of model parameters. In our case we have defined x = QG and λ = QG exp , where and is determined by the model parameters r QG 0 and γ. Taking the natural logarithm of both sides of eqn. (9), we obtain: where model independent terms have been dropped. We calculated the log-likelihood for a grid of (r QG 0 , γ) values which defines an uniform prior, ranging from 1.0 ≤ γ ≤ 5.0 and 1.0 ≤ r QG 0 ≤ 15.0 and maximized the likelihood to obtain r QG 0 = 6.93 h −1 cMpc and γ = 2.4. These values were used in eqn. (5) to calculate the expected χ(R min , R max ) value shown as the red line in Fig. 14. We also computed the 1σ and 2σ 2D confidence regions for these parameters, shown in the r QG 0 − γ plane in Fig. 15. We determined errors on the parameters by marginalization. Given that our grid of values is uniform, the normalized likelihood is the joint posterior distribution of the parameters P (r QG 0 , γ). Therefore, we marginalized out r QG 0 and γ to obtain the probability distributions P (γ) and P (r QG 0 ), respectively. From those probability distributions we computed 68% confidence regions about our MLE to define the error on the parameters. We find r QG 0 = 6.93 +2.13 −1.89 h −1 cMpc and γ = 2.4 +0.3 −0.5 . As shown in Fig. 15, our measurements are relatively noisy when we fit r QG 0 and γ simultaneously, and there is a clear degeneracy between these parameters. For that reason, following common practice, we also fit the correlation function with γ fixed. Independent measurements of QSO auto-correlation suggest a slope of γ = 2.0 (Shen et al. 2007), which lies within the 1σ confidence region of our measurement (see Fig. 15). Thus if we choose to fix the slope to this value, the maximum likelihood and the 1σ confidence interval for the cross-correlation length is r QG 0 = 8.83 +1.39 −1.51 h −1 cMpc. Note that in the analysis described above we have assumed that the error bars on our cross-correlation function are dominated by Poisson counting errors. This ignores cosmic variance fluctuations, and also assumes that the positions of the LBGs around the QSO are uncorrelated. However, as galaxies are not randomly dis- and γ parameters (in blue and red respectively), determined using a maximum likelihood estimator. The best estimation is shown as a white cross. tributed in the universe, but rather have significant autocorrelations, our binned measurements of QG are not truly independent. Given these correlations, our results could be somewhat sensitive to our choice of binning, and our error bars could also be somewhat underestimated. In principle, we should include the correlations and cosmic variance in our likelihood, analogous to computing the non-diagonal elements of the covariance matrix for a multivariate Gaussian likelihood. However, there is no simple analytical expression for the likelihood of a correlated Poisson process, and furthermore correctly modeling the cosmic variance would require the use of N-body simulations of massive QSO halos at z ∼ 4. Note however that while the positions of LBGs in the same field will be correlated, our QSO fields are separated by Gpc distances, and hence the positions of LBGs in different fields are completely independent. Given that our correlation function comes from six distinct fields, and the relatively large Poisson error bars, we believe that ignoring correlations and cosmic variance is a reasonable approximation.
Our measurement indicates a strong cross-correlation between QSOs and LBGs at z ∼ 4, implying that QSOs trace massive dark matter halos in the early universe, with detectable enhancements of LBGs. We expect that those halos evolve to the most massive cluster of galaxies at z = 0. Our results are in agreement with the expected cross-correlation function (r QG 0 = 9.6 h −1 cMpc for γ = 1.8) computed from the individual QSO and LBG autocorrelation functions assuming linear bias, as shown by the dashed line in Fig. 14.

Auto-Correlation of LBGs in QSO Fields
Another measure of the clustering of LBGs in QSO environments is the LBG auto-correlation function in our fields. If QSOs trace highly biased locations of the universe, then we expect the LBGs around them to be more highly clustered than LBGs in random fields, resulting in an enhancement of the LBG auto-correlation function. The auto-correlation function of z ∼ 4 LBGs in random fields was measured by Ouchi et al. (2004b), which we compare to our results.
To measure the LBG auto-correlation function we adopt the estimator: where GG is the number of observed LBG-LBG pairs, and RR is expected random number of LBG-LBG pairs, in a cylindrical volume defined by the radial bin [R min , R max ] and the height ∆Z. We measured GG directly from the images by counting the LBG pairs in each radial bin. Following the same argument in § 4.1, we used the LBG luminosity function to compute the background number density n G , rather than estimating it from our survey images.
We computed the expected random number of LBG pairs as (see e.g. Padmanabhan et al. 2007) where n G is the same quantity defined in § 4.1 and V eff is given by eqn. (8), but in this case using a different radial selection function φ R (R), because of the different binning used. Here N G is the expected number of LBGs for the entire volume in question in a random region of the universe, which is computed for each of our six fields as N G = n G V field . The radial selection function φ R (R) in this case is computed in an analogous way as for the cross-correlation: we created catalogs with N ran ∼ 100, 000 randomly distributed galaxies on our masked images, and then we computed φ R (R) as the ratio between the observed number of random galaxy pairs over the expected number of random galaxy pairs per radial bin. Here, the expected number of galaxy pairs per bin is computed by N ran n ran π(R 2 max − R 2 min ). Note according to eqn. (13) RR is proportional to the square of the LBG number density n G and to the square of the redshift selection function φ Z (Z) such that: We computed GG and RR for each individual field and then we stacked the counts to measure the binned χ(R min , R max ) value as in eqn. (12). We show the results in Fig. 16 and the numerical values are given in Table 6. We estimate errors on χ(R min , R max ) using the one-sided Poisson confidence intervals for small number statistics in the same way as in § 4.1.
Analogous to our approach for the cross-correlation, we used a MLE to fit our auto-correlation function. In this case the expected number of LBG-LBG pairs GG exp is modeled as: where ξ GG is the LBG auto-correlation function assumed to have a power law form with correlation length r GG 0 . For the fitting we used an uniform prior defined by 1.0 ≤ γ ≤ 2.5 and 5.0 ≤ r GG 0 ≤ 60.0. We show the 1σ and 2σ 2D confidence regions for the parameters in Fig. 17. We obtained that the maximum

4335.68
Note. -We present the data for the LBG auto-correlation function in QSOs fields χ(Rmin, Rmax) shown in Fig. 16. This is measured in radial bins defined by Rmin and Rmax. GG is the observed number of LBG-LBG pairs per bin, and RR is the expected number of random-random pairs per bin, computed from eqn. (13). We also show the total volume of the bin added over the fields, computed as V eff,total = Observed LBG Auto−correlation LBG Auto−correlation in random fields (Ouchi et al. 2004) Best fit Figure 16.
The data points are showing the LBGs autocorrelation measurement in QSO fields as we describe in section § 4.2. The solid red curve shows the best fit for our measurements given by r GG 0 = 21.59 h −1 cMpc and γ = 1.5. The dotted black curve shows the LBGs auto-correlation in blank fields at z ∼ 4 measured by Ouchi et al. (2004b). We find a stronger clustering in our fields in comparison with random fields, which suggests that QSOs are located in overdense regions.
likelihood and the 1σ confidence intervals are r GG 0 = 21.59 +3.73 −2.96 h −1 cMpc and γ = 1.5 +0.1 −0.2 , which is plotted as the red line in Fig. 16. Following the same arguments as in § 4.1, we also fit the auto-correlation function with γ fixed. LBG auto-correlation function measured in random locations at z ∼ 4 suggest a slope of γ = 1.8 (Ouchi et al. 2004b), which lies outside the 1σ confidence region of our measurement (see Fig. 17). We then prefer to fix γ to its maximum likelihood value γ = 1.5, which agrees with the LBG auto-correlation function slope measured in QSO fields at z = 2.7 (Trainor & Steidel 2012). After fixing γ, we obtain r GG 0 = 21.59 +1.72 −1.69 h −1 cMpc. In order to compare this clustering signal with the one computed in random fields, we use the LBG autocorrelation at z ∼ 4 measured by Ouchi et al. (2004b). Plugging their best fit values (r GG 0 = 4.1 h −1 cMpc and γ = 1.8) into eqn. (5) using a power law form for ξ GG (R, Z), gives the dotted line plotted in Fig. 16. To better compare our auto-correlation measurement with the Ouchi et al. (2004b) random field values we performed a fit with fixed γ = 1.8, obtaining a maxi- The same as Fig. 15 but for the r GG 0 and γ parameters corresponding to the LBG auto-correlation function. mum likelihood value and 1σ confidence interval given by r GG 0 = 16.86 +1.17 −1.14 h −1 cMpc, which is ∼ 4 times higher than the correlation length in random fields. The fact that our LBG auto-correlation measurement is higher, suggests that the LBGs in our fields are more clustered than LBGs in random fields, which provides another indication that QSOs fields trace regions of the universe that are denser than the cosmic average, confirming our findings from the cross-correlation measurement in section § 4.1.
While it may at first seem counter-intuitive that the LBG auto-correlation is enhanced by such a large factor (∼ 4) in QSO environs, this is actually exactly the expected behavior as we clarify here. The LBG autocorrelation function measures the radially binned profile of galaxy pairs, and it is not trivial to relate the autocorrelation to the cross-correlation. In order to build intuition, we will think in terms of the total number of galaxies detected in our survey (see Table 4). On average, we found 1.5 times more galaxies in QSO fields compared with the number of galaxies expected in random locations (which is estimated from our selection function and the number density of LBGs n G ), and that means that we should detect at least 1.5 2 more galaxy pairs in our fields compared with the expectation in random fields, simply because we are overdense by that factor. This implies that the auto-correlation function will never be less than 1.5 2 − 1 = 1.25.
To better illustrate what happens to the autocorrelation in an overdensity, we will consider a hypothetical scenario where galaxies are randomly distributed in the universe with number density n G , and are clustered only around QSOs, which are however rare objects in the universe. Now imagine that the number density of galaxies around QSOs is enhanced within a sphere of radius R QSO , but that galaxies are otherwise randomly distributed within the sphere. In other words, we imagine that in QSO fields the number density of galaxies is simply increased by a factor X, but that within the sphere they are unclustered. If QSOs are rare structures, then when averaging over large volumes of the universe, we expect that the observed number of galaxy pairs GG will be very close to the random expectation RR (computed using n G ), and then the galaxy auto-correlation function measured from random locations will be flat and very close to zero on small scales, r < R QSO . On larger scales, r ∼ R QSO , it would deviate more from zero, reflecting the clustering due to the top-hat overdensities around QSOs (however, if QSOs are very rare objects this positive correlation function would be diluted and could still be quite small). Note however that if we measure the galaxy auto-correlation around QSOs at radii r < R QSO , then GG will be X 2 times larger than RR (which is again computed using n G ), and then we would measure an auto-correlation function of roughly X 2 − 1. This could be much larger than the value measured in random locations which would be close to zero. This simple example illustrates that because RR is computed from the number density of galaxies in random regions n G , overdense fields will always result in an enhanced autocorrelation relative to that in random locations, and that these enhancements can be quite large. The situation clearly becomes more complicated if galaxies are intrinsically clustered with a power law profile, and there is no simple analytical relationship between the cross and auto correlation functions. To fully quantitatively understand the relationship between the cross and auto correlation functions in QSO environments, one would need to analyze cosmological simulations (see e.g. White et al. 2012). But the generic expectation is an enhancement of the auto-correlation function in QSO environs compared with blank field pointings, which exactly what we see in Fig. 16.

Comparison with Previous Measurements
The highest redshift for which the QSO-LBG crosscorrelation has been measured before is at z ∼ 3 by Trainor & Steidel (2012), who reported an overdensity of galaxies in QSO fields, and found a cross-correlation length of r QG 0,z∼3 = 7.3 ± 1.3 h −1 cMpc for a fixed γ = 1.5. At z ∼ 4 we find a steeper slope (we fixed γ = 2.0) than Trainor & Steidel (2012), but in order to facilitate a comparison with their results, we fit our crosscorrelation measurement for their same fixed γ = 1.5. This γ value, is near the border of our 2σ confidence region (see Fig. 15) and thus disfavored by our measurements, but we nevertheless proceed with this for comparison purposes. We obtain a cross-correlation length of r QG 0 = 10.73 +2.20 −2.41 h −1 cMpc, which is ∼ 1.5 times higher than their cross-correlation length at z ∼ 3, indicating that halos hosting QSOs are considerably more biased and highly clustered at z ∼ 4. This agrees with the result reported by Shen et al. (2007) who find that the QSO auto-correlation increase significantly from z = 3 to z = 4, and as such we expect to have found a larger cross-correlation.
At z > 4 only individual QSO fields have been studied so far. Some studies of QSO environments at z ∼ 6 − 7 find no enhancements of galaxies compared with the background (e.g Willott et al. 2005;Bañados et al. 2013;Simpson et al. 2014;Mazzucchelli et al. 2016), which could be suggesting that the strong QSO-galaxy crosscorrelation breaks down at those redshifts. The lack of QSO auto-correlation measurements at these high redshifts makes impossible to know masses of dark matter halos hosting z ∼ 6 QSOs, but if their masses are comparable to those hosting QSOs at z ∼ 4 (i.e M halo 10 12 M as suggested by the Shen et al. 2007 auto-correlation), then one would generically expect a strong QSO-galaxy clustering signal as we have detected here at z ∼ 4.

TESTING THE ROBUSTNESS OF OUR RESULTS
Two requirements must be fulfilled to ensure a robust clustering measurement: we need a low contamination level in the LBG sample and an accurate knowledge of the background number density of LBGs. Given that we used a novel NB technique to select LBGs, we need to carefully consider those requirements. In this section, we first discuss caveats related to the the use of this selection technique. Then we consider the effects of using a contaminated sample for clustering measurements, and finally we explore the impact of using different LBGs selection criteria on our results.

The Use of a NB Technique For LBG Selection:
Caveats A first complication of using our novel method for color-selecting LBGs is that the level of contamination of our sample is unknown. In principle, the purity of the sample can only be determined with follow-up spectroscopy, or detailed modeling of the population of contaminant galaxies. Both alternatives would be challenging to implement and beyond the scope of this paper, but as a compromise we qualitatively discuss the impact of contamination on the correlation function (see § 5.2) and we demonstrate the robustness of our results against contamination by exploring their sensitivity to the colorselection criteria (see § 5.3). Note however that we could excise contamination in our LBG sample if we had additional imaging on our fields using traditional broad band filters. This would allow us to confirm the presence of the Lyα break in our LBG candidates.
Another complication of using our novel color-selection is that we did not have an independent measurement of the background number density of LBG required to compute the clustering. This implied that we had to rely on Monte Carlo simulations to determine the LBG selection function φ Z (Z), and then our clustering results are sensitive to errors in this modeling. If the completeness of the sample were close to 100%, then 10% errors on φ Z (Z) would impact our measurement at the 10% level, whereas if the completeness were ∼ 20% (as is the case), then there could be 100% error on φ Z (Z), which could strongly impact the amplitude of the measured clustering. Note that the auto-correlation is even more sensitive to this quantity compared to the cross-correlation, because while QR is proportional to φ Z (Z), RR is proportional to the square of this quantity. In § 5.3 we test our redshift selection function to demonstrate that it is accurate and correctly modeled.

Impact of Contamination on the Clustering
Measurements One method to qualitatively check the contamination level in the sample is by studying the shape of the measured correlation function. For example, if we measure the cross-correlation function via eqn. (6) using a highly contaminated sample, the numerator in that equation would be overestimated, because of the inclusion of lowredshift contaminants which are taken to be real LBGs. However, since the denominator QR is simply computed from the LBG luminosity function and on our redshift selection function, this value does not include the extra number counts due to contamination. This implies that the measured cross-correlation will not behave like a power-law, but rather it will flatten toward larger scales. Quantitatively, for a contaminated sample, what we would actually measure is where QG and QR are given by eqns. (10) and (7) respectively, and N cont is the number of contaminants in the bin. Given that the contaminants are galaxies at different redshifts, the cross-correlation between them and the z ∼ 4 QSO is zero, then the number of contaminants will be given by N cont = n cont V eff,cont , where n cont is the number density of contaminants and V eff,cont is the effective volume of the bin, which is given by eqn. (8), but with the redshift selection function of the contaminants φ Z,cont (Z). Then the eqn. (16) reduces to: where χ true (R min , R max ) is the correlation function that we would measure from a non-contaminated sample (i. e, here χ true (R min , R max ) = QG / QR − 1), and D C (z) is the transverse comoving distance at redshift z. In the absence of contaminants, the second term in this equation would be zero, and we recover the correlation function defined in eqn. (5), which has a power law shape. However, if a large number of contaminants are included which span a large range in redshift, the second term becomes important, and given that it does not depend on radius, this same constant is added everywhere to the cross-correlation function flattening its shape, with the degree of flattening dependent on the level of contamination. This flattening effect will be even stronger for the autocorrelation function since it is proportional to the square of both the number density of contaminants and redshift range they cover. Then for a contaminated sample one obtains where ignored the clustering of the contaminants, which should be greatly diluted in projection if the contaminants span a large range of redshifts. To take clustering of contaminants into account, an additional term should be added to this equation to account for their auto-correlation. Therefore, the smoking gun of high contamination in the LBG sample would be a flat crosscorrelation and auto-correlation function in Figs. 14 and 16, respectively. Given that we measured a power law shape for both correlations, we believe that our LBG sample is not strongly affected by contamination.
We have explored a third way to check contamination which is also independent of our estimate of n G and φ Z (Z). For a highly contaminated sample that includes galaxies over a wide range of redshifts, it would be more appropriate to measure angular distances instead of transverse comoving distances. We thus compute the angular correlation function ω(θ) using the standard procedure, where RR is determined from the angular number density of the data itself, and we do not assume anything about the number density or selection function. In this case we only measure how clustered is our sample in comparison to a random distribution with the same number density as our sample. This angular correlation function calculation thus differs from our LBG auto-correlation function in § 4.2, where RR was computed from n G and our selection function φ Z (Z). For a highly contaminated sample we expect the angular correlation function to be close to zero on all angular scales, because the inclusion of uncorrelated galaxies over a broad redshift range, would dilute any real clustering signal. On the other hand, for a relatively pure sample composed primarily of LBGs at z = 3.78 ± 0.3, we expect to measure a power law angular auto-correlation because we would be selected only highly biased galaxies in a small volume. Note however that even for a pure LBG sample, the ω(θ) computed in this way is not the true angular correlation function of LBGs, because we are pointing towards overdense regions around QSOs.
We estimated the angular auto-correlation function of the LBGs as: where GG(θ) is the number of LBG-LBG pairs per angular bin, which is directly measured from our images, and RR(θ) is the number of random-random pairs per angular bin. The RR(θ) quantity was estimated using a random catalog of sources created as follows. First, we computed the total number of LBG candidates in all the fields, then we divided that by the total unmasked area to get the average angular number density of LBGs. Second, we multiplied the unmasked area per image by this average number density to determine the number of galaxies expected in each field. Finally, we increased the Angular auto-correlation function measurement for the LBGs sample. This measurement is used to test if the sample is contaminated. If the sample were highly contaminated ω(θ) = 0 at every scale and a power law fit with slope β would be consistent with zero. We find that the LBGs sample is not highly contaminated, given their power law shape in this plot, which is well fitted by ω(θ) = 21.56 θ −1.07 (red line). number of galaxies by a large factor F in order to decrease the noise in the measurement, and we randomly distributed those sources on the image and then we measured RR(θ) by counting the pairs of simulated galaxies per angular bin. We then re-scaled RR(θ) down by F 2 .
Our measurement of the angular correlation function is shown in Fig. 18. We see a non-flat correlation function which suggest that our LBG sample is not highly contaminated. Assuming a power law form given by ω(θ) = Aθ −β we performed a Levenberg-Marquardt least-squares fit to these data to quantify how consistent the measurement is with a flat shape (where β = 0). We obtain best fit parameters of A = 21.56 ± 39.54 and β = 1.07 ± 0.49. Given the large error bars in the measurement we are not able to discard a correlation function consistent with zero, however, as we show in the next subsection, if the LBG sample were highly contaminated then the angular correlation function would be much flatter. The fact that we measure a signal in Fig. 18 suggest that we are measuring real LBG clustering.

Robustness of Clustering Measurements Against
Changes in Color-Selection Here, we study the impact of using different colorselections on our clustering measurements to demonstrate that our results are not significantly impacted by contamination, and to show that our Monte Carlo simulation of the completeness is robust if we change the color-selection criteria. To this end we have defined different LBGs selection criteria, and for each one we compute the cross-correlation, auto-correlation, and angular auto-correlation function. The cross-correlation and auto-correlation functions were fitted in each case using a MLE following the same procedure described in § 4. We consider a progression of seven different selections, from the most permissive Case 1, which selects the majority of z ∼ 4 LBGs, but also likely incurs a large fraction of lowredshift contaminants, to the most conservative Case 7, which results in a low completeness for z ∼ 4 LBGs, but ensures low contamination. These results are shown in Fig. 19, and we tabulate the best fit values for each case in Table 7. There we also tabulate the best fit correlation lengths for a fixed γ = 2.0 for the QSO-LBG crosscorrelation, and γ = 1.5 for the LBG auto-correlation, in order to study how r 0 varies for the different cases. Note that for the three most permissive Cases 1-3, we measure a flat correlation function and hence do not quote fits for fixed γ. Additionally for those cases we had to use a different prior for the MLE fit, since flat correlation functions result in small values for the slope and large values for the correlation length, which are not covered by the prior used for the other cases. We only quote the best fit parameters for Cases 1-3 because the 1σ confidence region extends beyond the prior, precluding reliable error estimates.
We find that the cross-correlation function flattens and its amplitude increases for more permissive selections that increase the level of contamination, and the autocorrelation function shows a similar but even stronger tendency. This is the behavior that we expected as we describe in § 5.2 and according to eqns. (17) and (18). As for the angular correlation function, we find that the more conservative the selection, the steeper the slope of ω(θ) and the more significant its departure from zero. These are again the trends we expect because reduced contamination results in a more strongly clustered sample of z ∼ 4 galaxies, selected from a narrow redshift slice reducing the amount that the clustering is diluted by projection. Note however, that for the less conservative cases (i.e, Case 1 and Case 2), where the sample is dominated by contaminants, the angular correlation function is close to zero, but not perfectly consistent with β = 0. We believe that the measurement of a weak clustering signal in these cases results from the actual clustering of foreground contaminants, which is diluted by the line-ofsight projection, but nevertheless remains strong enough to not be perfectly consistent with zero.
The takeaway message from Fig. 19 is that we observe convergence of both the cross-correlation and autocorrelation functions for the more conservative selections. Specifically, we find stable results for Cases 5-7, with the only significant difference being the signal-to-noise ratio of the clustering measurements, resulting from the smaller sample of LBGs selected in the more conservative cases. In Fig. 20 and 21 we plot the values of the cross-correlation length r QG 0 (for fixed γ = 2.0) and autocorrelation length r GG 0 (for fixed γ = 1.5), respectively, for the four most conservative selections. The convergence of the correlation lengths demonstrates that: 1) we do not suffer large contamination and hence our results are robust against contamination, 2) that our Monte Carlo simulation of the selection function is reliable, since it results in consistent measurements as the colorselection and selection function are varied, 3) our results are largely independent of the exact color-selection region adopted. For these reasons we simply adopt Case 5 to present the final results in this paper.
Finally, we performed one last test to establish that the redshift selection function modeled from our Monte Carlo is essentially correct. We compared the total observed QSO-LBGs pairs in all the fields QG obs for each selection with the expected value QG exp based on our Table 7 Best fitted parameters for the cross-correlation and auto-correlation functions for the seven cases showed in Fig. 19.

Cross-Correlation
Auto clustering measurements, and our Monte Carlo determination of the redshift selection function. Specifically, for each selection the QG obs was measured by summing the observed QSO-LBGs pairs over the fields up to scales of R ∼ 9 h −1 cMpc, and Poisson errors were computed for this measurement. The expected value QG exp for each field was computed using eqn. (10), where we computed the corresponding φ Z (Z) using our Monte Carlo simulation method described in § 3.1 for each selection criteria. For all the cases, we used r QG 0 = 6.93 h −1 cMpc and γ = 2.4, which are the best fit parameters for our fiducial color-selection (Case 5; see Fig. 19) in the computation of QG exp . The total expected number of QSO-LBG pairs in the whole survey QG exp , was computed by summing QG exp over the bins and over the fields.
If the contamination is low, and the redshift selection function φ Z (Z) is correctly computed for each case, we expect that QG exp should equal to QG obs . As the sample becomes more contaminated we expect that QG obs will exceed QG exp and increasingly deviate from it for more permissive selections. The results of this test are shown in Fig. 22, where we plot QG obs versus QG exp for the seven color-selections we considered, and compare to the line QG exp = QG obs (solid line). We find that the total number of observed QSO-LBGs pairs is consistent with our expectations for the three more conservative selections Cases 5-7, but that QG obs exceeds QG exp for more permissive selections, with the deviations progressively increasing as more contaminants are included. Note that by construction we will have QG obs = QG exp for Case 5, since the clustering mea-surements ( QG obs ) were fit to determine the correlation function parameters, which go into the computation of QG exp . But the fact that expected QG exp matches the observed QG obs for the more conservative Cases 6 and 7 demonstrates that 1) the modeling of the redshift selection function φ Z (Z) is correct, 2) the contamination is insignificant, and 3) our clustering measurements are robust.

SUMMARY AND CONCLUSIONS
The strong observed auto-correlation of QSOs at z > 3.5 indicates that they inhabit massive dark matter halos with M halo > 10 12 M , which implies QSO environments should exhibit an enhancement of galaxies manifest as a strong QSO-galaxy cross-correlation function. We characterized the environments of six QSO fields at z = 3.78 that were chosen to host massive BHs ( 10 9 M ). The fields were imaged using VLT/FORS1 with two custom NB filters, and the broad band r GUNN , to identify LBGs using a novel technique which selects them in a redshift range ∼ 3.3 times smaller than the range typically probed when selecting LBGs with broad band filters. This significantly reduces the line-of-sight projection effects that have hampered previous searches for overdensities around z 5 QSOs.
Since we used a non-standard filter set to select LBGs, we performed detailed Monte Carlo simulations to model LBG colors, define our selection criteria, and compute the redshift selection function and volume probed by our survey. This new method effectively selects LBGs in a narrow redshift range, but the color loci of z 3.78 . Impact of contaminants on the clustering measurements. We show seven different selection criteria and their respective clustering measurements. From left to right we show the color-color plot showing the color cuts used and the photometry. The cross-correlation including the best fit according to our MLE estimator and the theoretical expectation of χ(R min , Rmax) calculated from the QSO and LBGs auto-correlation functions and assuming a linear bias model (dashed line). The auto-correlation measurement including the best fit according to our MLE estimator, and the LBGs auto-correlation in blank fields at z ∼ 4 measured by Ouchi et al. (2004b) (dotted line). Finally, the angular auto-correlation function, with a power law fit ω(θ) = Aθ −β , with the A, β values indicated in the top right corner. In the three correlation function plots we show the best fit in each case as a red curve and the best fit for our fiducial selection (i. e Case 5) as a blue curve. From top to bottom we show selections progressively more conservative, and then less contaminated. The selection used in this paper to measure the clustering properties of z ∼ 4 LBGs in QSO environments correspond to the Case 5. We detected a convergence in the clustering measurements for the three last cases.  values for a fixed γ = 1.5 for the four more conservative selections shown in Fig. 19. We detected a convergence of the correlation length. The horizontal dashed line indicate the best fitted r GG 0 value for the Case 5.
LBGs and low-redshift galaxies overlap more than with traditional LBG selection using broader filters. Defining a pure sample free low-redshift contaminants required adopting stricter color cuts, which decreased the completeness of the resulting LBG sample. We devised selection criteria which resulted in ∼ 26% completeness at z = 3.78, and detected 44 LBGs in our six fields, corresponding to a number density of 0.19 LBGs arcmin −2 . Our survey probes ∆z 0.3, and covers a volume equal to 14,782 h −3 cMpc 3 within R < 9 h −1 cMpc from the QSO, and we find on average 1.5 times more galaxies than expected in random locations of the universe.
Our work resulted in the first volume-averaged projected QSO-LBG cross-correlation function at z ∼ 4.
We fit our measurements with a (real-space) power-law cross-correlation function, and found r QG 0 = 6.93 +2.13 −1.89 h −1 cMpc and γ = 2.4 +0.3 −0.5 . When we fix the 10 100 <QG> exp 10 100 1000 <QG> obs C a s e 1 C a s e 2 C a s e 3 C a s e 4 C a s e 5 C a s e 6 C a s e 7 Figure 22. The total observed QSO-LBGs pairs in the whole survey for our selections with Poisson error bars compared with the expected number of QSO-LBGs pairs, assuming the cross-correlation parameters for the Case 5 and using the redshift selection function φ Z (Z) computed for each case as we described in § 3. The solid line indicates the case where QG obs = QG exp . Our results show that the three more conservative cases (noncontaminated samples) are consistent with the expectation, which means that the modeled redshift selection function φ Z (Z) is sensitive and that the level of contamination in those selections are negligible. slope at γ = 2.0 we find r QG 0 = 8.83 +1.39 −1.51 h −1 cMpc. This strong cross-correlation function is in agreement with the theoretical expectation for the cross-correlation assuming a linear bias model, which can be estimated using the auto-correlation of both LBGs and QSOs at z ∼ 4.
We also measured the auto-correlation function of LBGs near these QSOs and found an auto-correlation length of r GG 0 = 21.59 +1.72 −1.69 h −1 cMpc for a fixed slope of γ = 1.5 which is ∼ 4 times higher than the measured auto-correlation length of LBGs in random fields at the same redshift. Our measurement of an enhanced LBG auto-correlation in QSO environments, and the strong QSO-LBG cross-correlation both indicate that QSOs at z ∼ 4 trace massive dark matter halos in the early universe, which are the likely progenitors of massive cluster of galaxies at z = 0.
We demonstrated that our results are robust against contamination and that our selection function modeling is reliable, by varying our color-selection criteria and showing that the cross-correlation and auto-correlation functions are converged. Spectroscopic follow-up of our LBGs candidates would provide an additional and definitive test of the reliability of our novel color-selection technique. However, because the colors of LBGs in our filters do not separate as cleanly from contaminants as for broad-band LBG selection, we had to choose relatively conservative color cuts that recovered only 26% of LBGs. As such, we believe that the preferred approach to search for overdensities around z ∼ 4 − 6 QSOs using NB filters is to perform traditional LAE selection. Although LAE selection also only selects a fraction of the total population of high-redshift galaxies (Stark et al. 2010(Stark et al. , 2011Curtis-Lake et al. 2012), the primary advantages are: 1) contamination (from low-redshift line-emitters) is very low, 2) the background number density is known from wide-field observations of blank fields (Hu et al. 2004;