Quasar clustering at redshift 6

Large-scale surveys over the last years have revealed about 300 QSOs at redshift above 6. Follow-up observations identified surprising properties, such as the very high black hole (BH) masses, spatial correlations with surrounding cold gas of the host galaxy, or high CIV-MgII velocity shifts. In particular, the discovery of luminous high-redshift quasars suggests that at least some black holes likely have large masses at birth and grow efficiently. We aim at quantifying quasar pairs at high redshift for a large sample of objects. This provides a new key constraint on a combination of parameters related to the origin and assembly for the most massive black holes: BH formation efficiency and clustering, growth efficiency and relative contribution of BH mergers. We observed 116 spectroscopically confirmed QSOs around redshift 6 with the simultaneous 7-channel imager GROND in order to search for companions. Applying identical colour-colour cuts as for those which led to the spectroscopically confirmed QSO, we perform LePHARE fits to the 26 best QSO pair candidates, and obtained spectroscopic observations for 11 of those. e do not find any QSO pair with a companion brighter than M1450(AB)<-26 mag within our 0.1-3.3 h^-1 cMpc search radius, in contrast to the serendipitous findings in the redshift range 4--5. However, a low fraction of such pairs at this luminosity and redshift is consistent with indications from present-day cosmological-scale galaxy evolution models. In turn, the incidence of L- and T-type brown dwarfs which occupy a similar colour space as z ~ 6 QSOs, is higher than expected, by a factor of 5 and 20, respectively.


Introduction
Clustering of quasars provides one of the most important observational constraints to probe their physical properties, formation and evolution (Haiman & Hui 2001) as well as the mass of their dark matter (DM) haloes (Sheth et al. 2001). As easily detectable objects throughout the Universe, quasars trace the underlying dark matter distribution (e.g. Haehnelt & Nusser 1999), and thus provide a powerful test of hierarchical structure formation theory (e.g. Fang 1989). How this structure formation in the early Universe proceeds in detail is far from understood. How many massive galaxies form in a single massive dark matter halo? Do all of these massive galaxies have a massive black hole (BH)? Do all these black holes accrete from the rich gas reservoir, or what are the duty cycles of QSO activity? Do massive galaxies form sub-structures and/or multiple massive black holes?
Particularly pressing questions over the last years are that of the formation channels and efficiencies of massive BHs, and of their growth rates with cosmic time (e.g. Johnson et al. 2013;Valiante et al. 2017Valiante et al. , 2018Inayoshi et al. 2020). Quantifying the number of QSO pairs at high-z is one observational avenue to constrain a combination of parameters of BH formation and BH growth rate. High-z QSOs can be as massive as the most massive BHs found in the local Universe (e.g. Bañados et al. 2018;Yang et al. 2020;Wang et al. 2021), and so must have grown Send offprint requests to: J. Greiner, jcg@mpe.mpg.de several orders of magnitude in less than 800 Myrs. Indeed, BH formation mechanisms predict the formation of BH seeds with at most masses of a few ∼10 5−6 M (Volonteri 2010;Inayoshi et al. 2020). If the growth of quasars was mostly through gas accretion, then they must have grown at the Eddington rate a large fraction of their lives. Finding QSO pairs, or the absence of, can help us understand if BH mergers play a role in their growth, and if dense environments can sustain efficient growth for several QSOs in the same region of the sky, or if baryonic processes (e.g., AGN feedback) prevent them from doing so (Wyithe et al. 2005;Kashikawa et al. 2007;Utsumi et al. 2010;Costa et al. 2014;Habouzit et al. 2019, for QSO feedback preventing galaxy formation/growth in their surroundings). Radiation from a QSO can impact a region from sub Mpc to tens of Mpc, depending on the QSO mass, accretion rate, and properties of its surrounding environment, and could therefore also inhibit the accretion activity of nearby BHs. In addition, quantifying the number of QSO pairs also provides us with constraints on BH formation efficiency, i.e. whether BH formation takes place in a large number of halos, but also on whether the formation of massive seeds could be clustered. The existence of many pairs of QSOs could indeed imply that massive seeds can form close to each other. Constraints from QSO pairs on QSO origins and growth with time are likely degenerated, but are key to improve our understanding of these extreme objects. Quasar clustering has been investigated at different scales and depths, among others using the Sloan Digital Sky Survey (Shen et al. 2007;Hennawi et al. 2010), the 2dF QSO Redshift and Dark Energy Survey (Porciani et al. 2007;Chehade et al. 2016), or the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) ). Studies at different depth for a given redshift range are important since a large luminosity dependence is theoretically predicted at higher redshifts (Hopkins et al. 2007). Clustering at large scales, sampling Mpc or larger distances, is stronger at high redshifts (z > 3.5) than at low redshifts (2.9 < z < 3.5) (Shen et al. 2007). Clustering at small scales, sampling the 10 kpc -1 Mpc scale environment, also increases from smaller (z ∼ 3) to larger (z ∼ 4) redshift, but shows a much stronger amplitude on kpc scales than on Mpc scales, e.g. a factor 4 stronger at z ∼ 2 (Eftekharzadeh et al. 2017).
Clustering measurements at high redshift (z > 3) are confronted with a problem: close quasar pairs are extremely rare. For comparison, the mean quasar separation at z > 3 is ∼ 150h −1 Mpc (given the standard luminosity function of Φ ∼ 10 −6 Mpc −1 mag −1 for M(AB) < −24; e.g. Willott et al. 2010;Ross et al. 2013;Matsuoka et al. 2018;Shen et al. 2020) and at small scales ( ∼ < 1 Mpc), the correlation function does not increase as fast as the volume decreases . However, direct observations of pairs might suggest differently -the three hitherto highest-redshift QSO pairs were all found serendipitously: (i) Schneider et al. (2000) reported a pair at z=4.26 with a separation of 33 (corresponding to 160h −1 kpc proper on the sky) which was found in the slit while spectroscopically confirming another QSO candidate; (ii) Djorgovski et al. (2003) discovered a z=4.96 quasar at a separation of 196 (or proper transverse separation of 0.9 h −1 Mpc) from a z=5.02 quasar discovered by Fan et al. (1999), and (iii) McGreer et al. (2016) reported a pair at z=5.02 with a separation of 21 (90h −1 kpc), discovered in a faint z ∼ 5 quasar (but not pair) search down to i = 23 mag, and selected simultaneously with identical colour criteria. Given the relatively small sample sizes which these discoveries are based upon, they suggest that quasars are even stronger clustered at redshifts > 4 on small scales. This in turn has been used to argue that feedback is very inefficient at high redshift (Willott et al. 2010;McGreer et al. 2016), since inefficient feedback would correspond to the correlation length to flatten out at high redshift, while concordant QSO and host halo mass growth would imply a sharply rising correlation length. In addition, merger models of galaxy formation predict a large fraction of binary QSOs, at least over some time span (e.g. Volonteri et al. 2003a), with the fraction of binary (off-centre and dual) QSOs increasing with redshift (Volonteri et al. 2016). This state of knowledge, and the importance of the conclusions which can be drawn, certainly warrants further searches for clustering properties at high redshift.
Here, we search a sample of 116 spectroscopically confirmed QSOs at redshift 6 (collected from the compilation of Bañados et al. 2016) for a companion. Throughout this paper we use the best-fit Planck cosmological model with Ω m = 0.286, Ω Λ = 0.714, and h = 0.69 (Planck Collaboration 2016). At a mean sample redshift of 6, an angular separation of 1 corresponds to a proper (comoving) transverse separation of 4 h −1 kpc (28 h −1 kpc). All magnitudes are in the AB system. (arcmin) 280 arcsec 1.1 h −1 Mpc proper transverse separation GROND g'r'i'z' GROND JHK Fig. 1. A schematic of the GROND optical (g r i z ) and near-infrared (JHK) channel spatial coverage at a redshift of 6. The dashed boxes are the respective g r i z FOV for each individual telescope dither position; a similar pattern applies to the JHK channels (not shown). With most of our targets being within 1 arcmin from the centre of the FOV, and since we require the r i z bands, our conservative search radius is 80 , corresponding to about 471 h −1 kpc proper, or 3.3 h −1 Mpc comoving.

Selection
From the sample of 173 spectroscopically confirmed QSOs at z > 5.6 compiled in Bañados et al. (2016) we have selected all 144 sources south of Decl. < +30 • , thus being visible from La Silla (Chile). We have added 6 objects from Mazzucchelli et al. (2017b), and 2 object from Jiang et al. (2016). For the observing campaigns, we have prioritised brighter sources, so that it is easier (less time consuming) to reach limiting magnitudes fainter than the spectroscopically confirmed QSOs. After weather and technical losses, we obtained useful new GROND observations for 77 sources, in addition to those 42 objects from the Bañados et al. (2016) sample, which had already been observed with GROND in 2013-2015 during the selection process towards their spectroscopic sample. Excluding in the following two QSOs at z > 6.6, our sample of 116 QSOs spans a redshift range from 5.6-6.6 and a z -band brightness of 18.6-23.1 mag.

Strategy, Data Analysis and Source Detection
Simultaneous imaging in g r i z JHK s with GROND   Since the field of view of the 4 visual channels g r i z is smaller than that of the NIR channels ( Fig. 1), the former define the search distance from the spectroscopically identified QSO to 80 distance. Thus, our maximum search radius for a QSO Colour-colour diagram in the GROND-intrinsic i z JHK s bands of the prime quasars (blue; triangles denote limits) compared to normal foreground stars and galaxies (grey) of the same fields, with selected brown dwarfs from earlier GROND observations indicated in green and purple, and labeled with their spectral type. Our 26 candidates are shown in red, with light red those 9 objects which have ALLWISE detections. The yellow-shaded region indicates the search region for the QSO pair. The arrow at the top right of each panel shows the amount by which a correction for A V = 1 would shift an object. Except for the grey objects, all are foreground-A V corrected. companion is 0.471 h −1 Mpc proper or 3.3 h −1 Mpc comoving (see Fig. 1). This accounts for the fact, that most of the primary QSOs are not exactly centred in the GROND field of view.
The question to answer before the start of our observing campaign was: how much deeper to look for a companion, relative to the prime QSO? We have picked a 2 mag difference as threshold, based on three reasons: (i) given the QSO luminosity function (Willott et al. 2010;Ross et al. 2013;Shen et al. 2020), this corresponds to a factor more than 10 in the number of objects, and should imply a large enough sample depth, (ii) optical variability is rare above ∼0.5 mag amplitude, ensuing that the effective difference is not diminished substantially, and (iii) the pair at z = 5, found by McGreer et al. (2016) also has such an i -band magnitude difference (19.4 vs. 21.4 mag). Thus, exposure lengths of typically 20 minutes was decided upon, ensuring that the 3σ limiting magnitude in the z -band was at least 2 mag deeper than the brightness of the spectroscopically identified QSOs (see below). While this 2 mag difference might suggest that the search could have been done just on Pan-STARRS data with its typical 5σ depth of z <22.3 mag, only 11 PS1 QSOs are bright enough to also survive the colour cut of iz > 1.5 mag at that faint z brightness (see below).
GROND data have been reduced in the standard manner ) using pyraf/IRAF (Tody 1993;Küpcü Yoldaş et al. 2008). The optical/NIR imaging was calibrated against the Sloan Digital Sky Survey (SDSS) 1 (Eisenstein et al. 2011), Pan-STARRS1 2 (Chambers et al. 2016) or the SkyMapper Survey 3 (Wolf et al. 2018) catalogs for g r i z , and the 2MASS catalog (Skrutskie et al. 2006) for the JHK s bands. This results in typical absolute accuracies of ±0.03 mag in g r i z and ±0.05 mag in JHK s . Since the GROND dichroics were built after the Sloan filter system , the colour terms are very small, below 0.01 mag, except for the i -band which is substantially narrower than the SDSS i'-band. All our analysis is done in the SDSS system, with the following transformations used: g S DS S − g GROND = (−0.006 ± 0.014) + (0.015 ± 0.025) * (g S DS S −r S DS S ), r S DS S −r GROND = (−0.004±0.004) For the fields calibrated against Skymapper and PS1, we use their conversion to SDSS (Scolnic et al. 2015). The PS1 z -band has a substantial colour term (z GROND = z PS 1 -0.214×(z PS 1 -y PS 1 ) due to the missing tail beyond 920 nm, but this becomes important only for redshift above 6.6 when Ly-α moves beyond the z -band limit.) The calibrated z -and J-band images for each field were stacked after re-sampling due to their different pixel sizes (0. 15 for the visual, and 0. 59 for the near-infrared); for the four objects at z > 6.5 we just used the J-band images. A source detection procedure on each stacked image then provides a list of all objects per field which then is used as input for 'forced' photometry on the individual 7-band images. This provides upper limits for the bands where the sources are not detected; otherwise it performs PSF photometry for g r i z , and aperture photometry in JHK (with a radius of 1×FWHM) at the position of the source in the master catalog which is allowed to re-centre by ±0.3 of the PSF width. Table 2 contains the details of the observation for each QSO, including the 3σ limiting magnitude in the z -band. Table 3 contains the foreground-A V corrected AB magnitudes of all 116 QSOs in the individual GROND filter bands.

Colour-colour selection of Candidates
With these GROND data in hand, the goal was to search the field around each spectroscopically confirmed QSO for a nearby source with a similar colour, taking advantage of the simultaneous multi-filter optical/NIR spectral energy distributions (SEDs) of all objects inside the field of view (FOV).
We implicitly assume that the searches leading to the discovery of the QSOs have picked the brighter of the pairs; this is corroborated by modelling which suggests that current flux-limited surveys of QSOs at high redshift preferentially detect objects at their peak luminosity, and thus miss a substantial population of similarly massive black holes accreting at lower accretion rates (Costa et al. 2014). Thus, here we are looking for objects fainter than the original QSO.
In all 116 fields we apply the following two colour selections (all in AB): iz > 1.5 mag, and z -J > -0.7 mag. We checked that there are no detections in the g and r bands. For the z > 6.5 objects, we adjust the z -J colour to -0.5, add J − H < -0.5, and require a i non-detection. This returns 74 new highz quasar candidates. We excluded 48 candidates through visual inspection of the 7-channel image cut-outs of these candidates, based on either source confusion, i.e. either a nearby source with overlapping PSF, or an extended source which could either be an interloper galaxy, or a non-resolved equally bright pair of objects. The remaining 26 candidate sources are shown in Fig. 2 as red (light and dark) dots, overplotted over the distribution of the spectroscopically confirmed quasars (shown in blue), which by itself show that our colour cut criteria recover all but 13 objects (with most of these failing due to non-detections, i.e. upper limits, in the i -band, i.e. deeper i -band exposure would likely recover those as well).
For 9 of the remaining 26 candidate sources we found catalog entries in the ALLWISE all-sky catalog (Wright et al. 2010;Mainzer et al. 2011) within 1 arcsec distance, and after visual checking of the images accepted these as matching objects. Tables 5 and 6 contain the photometry of these 26 candidates, separately for those with and without ALLWISE matches.
We then performed Le PHARE (Arnouts et al. 1999;Ilbert et al. 2006) fits to the photometry of the 26 candidates (including the ALLWISE magnitudes when available). We included an alternative IGM opacity model (Songaila 2004) in addition to the default one (Madau et al. 1999), but otherwise used the default spectral templates except for adding about 200 cloned highz QSO templates based on stacked observed SDSS spectra of low-z QSOs, similar to . The goal of these Le PHARE fits was to distinguish late-type dwarfs from QSO candidates, where the distinction is based on two observational features: the steepness of the SED in the i -z -J range, and the SED shape in the HK range which is much redder for dwarfs. Since dwarfs in our sensitivity range would be at distances of order 100 pc, their photometry would not be affected by interstellar dust. We thus created two SEDs for each candidate, one without A V -correction (for Le PHARE fits with dwarf templates) and one with full galactic A V -correction (for Le PHARE fits with QSO templates). For 14 (out of the 26) candidates the best-fit template is that of a dwarf (see last column of Tabs. 5, 6), for 10 candidates that of a QSO at high redshift (labeled in boldface in Tabs. 6), for one a galaxy at intermediate redshift, and for 1 candidate (SDSSJ2054-0005_4) the reduced χ 2 does not provide a clear favorite. The Le PHARE fits of the 10 high-z quasar candidates (i.e. with A V -correction applied) are shown in Fig. 6, which also contains the best-fit photo-z (6th field in the red line of labels), all within 13% of the spectroscopic redshift except for PSO J055.4244-00.8035 (the source we did not get a spectrum of).

Spectroscopic Observations
We obtained spectroscopy of 11 objects (9 of the 10 QSO candidates, the galaxy, and the unidentified object), mostly with the GMOS instruments of the Gemini Observatory (under proposal IDs GN-2019B-FT-201 and GS-2019B-FT-202). One object was observed with the FIRE (Folded-port InfraRed Echellette) spectrograph in Prism mode at the Magellan Baade telescope. A log of the observations is given in Tab. 4. At Gemini-North (South), a Hamamatsu detector was used with the R150 grating and filter G5308 (G5326). The spectra cover the range from 615 (750) to 1080 nm. Wavelength calibration was provided by CuAr comparison lamps. Flux calibration was done against the standard stars G191B2B (LTT7987). For the FIRE spectrum, we used archival A0V stars for telluric and flux calibration. The data were reduced in the standard manner, using GMOS specific routines provided within the IRAF (Tody 1993) package, and a custommade python code for the FIRE observation.
For object PSOJ071.0322-04.5591 no trace is visible in the 2D spectrum; for the other sources the optimally-extracted 1D spectra are shown in Fig. 7. None shows a sign of the Ly-α 'step' at the wavelengths suggested by the photometric redshifts (vertical lines in Fig. 7).

Pair candidates
We do not find a single QSO pair candidate brighter than M 1450 (AB)< −26 mag in our GROND data of 116 spectroscopically confirmed redshift 6 quasars. Except for two candidate object (near PSOJ055.4244-00.8035 and PSOJ071.0322-04.5591) we have obtained optical spectroscopy for all other Le PHAREderived QSO pair candidates, but could not confirm the quasar nature as suggested by the colour selection and SED fitting. This is somewhat surprising, since the (post-facto) photometric redshifts we derived with Le PHARE for the prime QSOs where in good agreement with the spectroscopic redshifts reported by Bañados et al. (2016), with only a 5% fraction mismatch. We consider it very unlikely that there is an unrecognised quasar among our identified candidates in Tabs. 5 and 6. This is in contrast to our empirical expectation (at the start of the project) of 3 pairs in our sample, based on the serendipitously discovered pair at z=4.26 (Schneider et al. 2000) and the two pairs at z=5.02 (Djorgovski et al. 2003;McGreer et al. 2016) (see details in the introduction), which when averaged together (1/100, 1/14, 1/47) suggest one pair every 40 QSOs, or 2.5%. This line of reasoning includes a few simplifications. First, it uses post-facto arguments and guesses on the search volume of the above serendipitous discoveries. Next, it assumed that the quasar clustering fraction does not evolve with redshift -the question to be answered by the project. Also, it ignored the different luminosity ranges sampled at different redshifts -the integral over the luminosity function is certainly the most sensitive component contributing to the w p statistic (see below).
In order to put this in context, we follow Hennawi et al. (2006) and Shen et al. (2010, in particular their eq. 4) in estimating the projected correlation function, or w p statistics, as defined by e.g., Davis & Peebles (1983). Given the mean seeing in our images and correspondingly the ability to separate two sources, we take an inner radius of 100 h −1 kpc (corresponding to 3 ) for the projected comoving area of the cylindrical annulus (the outer one being our search radius of 3.3 h −1 Mpc). We use the latest luminosity function description of Shen et al. (2020), as well as their publicly provided tools (https://github.com/gkulkarni/QLF) to convert our observed z (AB) magnitude limit and mean redshift to M 1450 , and applying the corresponding bolometric correction (using their 'global fit A'; we verified our procedure against their Fig. 5). This results in a number density of 9×10 −10 Mpc −3 for quasars with M 1450 < −26 mag, for our redshift range 5.6 < z < 6.6. We obtain an upper limit of w p < 110 000 h −1 Mpc as shown in Fig. 3, where we use a formal upper limit of our search of 1 pair out of the 116 objects searched for, given the missing spectroscopy for PSOJ055.4244-00.8035. The figure also includes measurements at lower redshifts from Schneider et al. (2000), Shen et al. (2010) and McGreer et al. (2016) (modified by the updated luminosity function and corrected for our search scale according to the γ = 2 powerlaw dependence of w p (R) ∝ R −γ from Shen et al. (2010)). Overall, our upper limit corresponds to a similar w p as previous findings at lower redshift, suggesting that there is no strong rise of clustering towards high redshift.
Using eq. 5 in Shen et al. (2007) and γ = 2, our upper limit on w p corresponds to a limit of the auto-correlation length r o < 1000h −1 Mpc.  2016), also extrapolated to our search scale. Horizontal 'error' bars visualise the corresponding redshift range of the studies. The labels at the data points report the limiting i/z-band AB magnitude of the search and the corresponding absolute M 1450 (AB) limiting magnitude. The two black dashed lines depict w p (R) for the same rate of (7 and 8, respectively) pairs as in Shen et al. (2010) extrapolated to higher redshift with their same M 1450 (AB), and unchanged luminosity function, suggesting higher w p (R) at larger redshift. Similarly, the dotted line shows the extrapolation for the case of constant limiting z-band magnitude. The red line depicts the backwards extrapolation of our upper limit to lower redshifts, at our M 1450 (AB)< −26 mag sensitivity, demonstrating the mismatch with respect to the findings at lower redshift. All the lines are computed via integration over the luminosity function of Shen et al. (2020) in redshift-steps of 0.25, and proper redshift-dependent transformation between observed z and absolute magnitudes.

Completeness and biases
In the past, clustering studies have typically relied on selecting pairs from large catalogs (e.g. SDSS), and the most important question is about the completeness of the underlying catalog (e.g. Hennawi et al. 2006;Eftekharzadeh et al. 2017). Here, we have obtained follow-up observations of known quasars at known redshift. Thus, in our case the question of completeness reduces to the questions of (i) whether or not the selection criteria which led to the sample of these known redshift quasars had a bias against selecting binary quasars? and (ii) whether our search strategy misses QSOs? As to item (i), one potential bias is that the seeing limit in the original survey data affects the discovery of QSO on 1-2 scales due to nearby contaminating sources. As this is a statistical problem, this is not expected to affect the fraction of binary QSO in the spectroscopically confirmed QSO sample.
As to item (ii), this splits into two sub-questions: First, what is the completeness of our colour-colour cuts? Since we ran our colour-colour selection blindly on all sources in each observa-tion, we can answer this with the recovery rate of the spectroscopically confirmed QSO: out of the 116 QSO, we miss 9, with 3 of those not detected in z , and the other 6 due to insufficiently deep i -band limits. Since the 3 non-detections are due to variability, they are not counted here (see below on the variability bias), so our incompleteness of the colour-colour selection is about 5%. Second, what is the completeness of the Le PHARE fits in recognizing the QSO nature? Again, this can be answered with the sample of the spectroscopically confirmed QSOs: ignoring sources with detections in less than 3 photometric bands, we ran Le PHARE fits on the remaining 97 spectroscopically confirmed QSOs in the same manner as with the candidates (i.e. two runs, one with and one without A V correction).
Only 5 sources are classified as dwarfs, all other 92 as QSOs. We double-checked that the colours and photometric errors of those 5 are not outliers with respect to the sample distribution of the 92 correctly classified QSOs. So this corresponds to only <5% incompleteness.
Further, we note that it is unlikely that even substantial optical variability of the QSOs would affect the fraction of pairs, since statistically one expects the same number of "risers" vs. "faders". The only exception would be effects due to Malmquist bias, but since only a handful of objects is detected with a z -band uncertainty larger than 0.2 mag, this imbalance is unimportant.
Lastly, a potential bias could be due to dust obscuration of AGN, in our case obscuring the companion of the spectroscopically confirmed QSO. Again, this is not expected to be a major issue. First, previous suggestions were that obscuration could primarily happen in the highest mass accretiion phases, one would first expect the brighter object of the pair to be obscured, rather than the lower-luminosity companions of the unobscured spectroscopicall confirmed QSOs. Second, while ALMA observations revealed substantial dust at high redshift, there is no evidence that these dusty objects are obscured AGN (see sect. 3.4 below). Finally, theory of dust production fails to explain dust in large amounts at z > 5 (Lesniewska & Michalowski 2019), preventing quantitative estimates or simulations.
Another potential bias could come in if black holes of similar mass have different Eddington ratios at different redshift, such that only the optically brightest QSOs are found at high redshift. However, Mazzucchelli et al. (2017b) have shown that (at least for the biased high-z sample they considered) BH masses and Eddington ratios of the z > 6.5 and SDSS 0.35 < z < 2.25 samples are consistent in a luminosity-matched sample. Thus, our assumption was that the clustering fraction at redshift 4-5 (Schneider et al. 2000;McGreer et al. 2016) is not much different than at redshift 6, leading to our expectation of more pairs as compared to our findings.
Further, our search comes with the biases, that the luminosity difference between the two members of the pair can of course be larger than our selected 2 mag, due to at least two different reasons: firstly, the BH masses and/or accretion rates can differ substantially, and secondly the Lyα emission which contributes to the z-band magnitude, can be drastically different between the two members of the pair, replicating the differences in the spectroscopically confirmed QSO sample (see e.g. Bañados et al. 2016). As to the first reason, simulations suggest differently. Based on the large-scale cosmological hydrodynamical Horizon-AGN simulations, Volonteri et al. (2016) find that the luminosity ratio between the central BH (the more massive one) and an off-centre BH is within a factor of ten (rms = 5x) for the cases of L bol > 10 43 erg/s and BH masses up to 10 8 M . While this corresponds to the ratio we employed, it should be noted that this ratio depends on selected input parameters for the simula-Article number, page 5 of 20 A&A proofs: manuscript no. QSOPairs_ed tion, among others a mass ratio <1:6 for the merging galaxies at which AGN activity is triggered in both black holes. Similarly, Bhowmick et al. (2019) find that in the MassiveBlack II simulations the satellite quasar luminosity is similar to that of central quasars. Finally, we note that for 41 of our targets the 3σ limiting z -magnitude reached in our observations is actually 3 mag fainter than that of the spectroscopically confirmed QSO, which still would lead to an expectation of one pair in our sample.

Comparison to simulations
The upper limit to the projected correlation function for bright QSO pairs at z ∼ 6 presented in Section 3.1 serves as an important constraint for galaxy evolution models. Currently, cosmological-scale simulations typically either focus on quasar properties below z ∼ 5 due to improved statistics ( Bhowmick et al. (2019) find no pairs of QSOs with g < 20.85 within a search radius of 4 Mpc at z = 2, and only three such systems by z = 0.6.
At higher redshift, work on the BlueTides simulation has shown that QSOs bright enough to match observations are also very rare, even when accounting for large black hole seeds (Di Matteo et al. 2017). There is only one black hole of mass above 10 8 M present in the BlueTides (400 h −1 Mpc) 3 box at z = 8, which has a luminosity of L bol ∼ 2.2 × 10 46 erg/s (Di Matteo et al. 2017;Tenneti et al. 2019).
However, despite the lack of pairs, galaxy evolution models still find that brighter QSOs cluster more strongly than their fainter counterparts, particularly on smaller scales. For example, on scales of ∼ 10 kpc at z = 4, DeGraf & Sijacki (2017) find the two-point auto-correlation function, ξ(r), for AGN with L bol ≥ 10 44 erg/s in the Illustris simulation is higher by a fac-tor of three compared to that for all AGN. On the larger scales probed by our observational measurement (∼ 0.1 − 3.3 h −1 Mpc), this enhancement is reduced but still remains present for the brightest objects, as also seen in the L-Galaxies semi-analytic model (Bonoli et al. 2009) and the Horizon-AGN simulation (Volonteri et al. 2016). All three of these models report a value of ξ(r = 1.2 h −1 Mpc) ∼ 10 at z = 4 for AGN with L bol ≥ 10 44 erg/s, with the expectation that this value will be higher at z = 6, as also seen in the observations discussed in Section 1.
In order to better probe the exact epoch and search volume of our observational QSO sample, we have carried-out a census of bright AGN pairs at z = 6 in the Illustris-TNG300 (hereafter, TNG300) and L-Galaxies 2020 galaxy evolution models. TNG300 Springel et al. 2018;Naiman et al. 2018;Marinacci et al. 2018;Nelson et al. 2018) is a magnetohydrodynamical simulation run in a (205 h −1 Mpc) 3 box, using a Bondi-Hoyle-Lyttleton BH accretion model with no boost factor but an assumed BH seed mass of 8 × 10 5 M (see Habouzit et al. 2021). L-Galaxies 2020 (Henriques et al. 2020) is a semianalytic model run on the larger (480 h −1 Mpc) 3 Millennium-I box (Springel et al. 2005), but with a more simplified phenomenological model of quasar-and radio-mode BH accretion and no BH seed mass (see Croton et al. 2006). In both cases, we assume that all AGN are radiatively efficient when calculating their luminosities, which provides an upper limit on the expected number of bright pairs. Therefore, the AGN bolometric luminosity for both models is calculated here as where r is the assumed radiative efficiency of the accretion disc,Ṁ BH is the BH accretion rate, and c is the speed of light. Given that the detection limit for the fainter of the two QSOs in our observational search is M B ∼ − 26 mag (i.e. L B ∼ 8.25 × 10 45 erg/s), we assume a minimum bolometric luminosity of L bol,min = 4 × 10 46 erg/s when searching for BH pairs in the models, given the bolometric conversion factor of 5.13 recommended by Duras et al. (2020).
For the same search area as our observational sample (i.e. an inner radius of 125 h −1 kpc and outer radius of 3.3 h −1 Mpc comoving, see Section 3.1), we find no BH pairs with L bol ≥ 4×10 46 erg/s in either TNG300 or L-Galaxies 2020 at z = 6. This holds true even when doubling the assumed radiative efficiency from r = 0.1 to 0.2. Indeed, there are no BHs at that luminosity at all within the TNG300 box above z = 5 when assuming r = 0.2. This theoretical result is nicely consistent with the null result returned by our observational search at z = 6, especially when considering the relatively large box sizes of L-Galaxies 2020 and TNG300.
We do not find any BH pairs at z = 5 in TNG300 or L-Galaxies 2020 that match the serendipitous discovery by Mc-Greer et al. (2016) of a bright pair with M 1450 ≤ − 25.4 (i.e. L bol ≥ 2.4 × 10 46 erg/s) separated by 21 (i.e. 816 ckpc) at z = 5.02. These simulations are therefore somewhat in tension with this discovery, although the exact search area to use is not well constrained. However, we do find pairs in TNG300 that match the serendipitous discovery on a BH pair with M B ≤ −23.7 mag (i.e. L bol ≥ 5.1×10 45 erg/s) by Schneider et al. (2000) at z = 4.26. When assuming r = 0.1, and a standard search volume around the mean deprojected separation of r = 1.92 cMpc with inner radius r − (r/4) = 1.44 cMpc and outer radius r+(r/4) = 2.40 cMpc, we find 14 BH pairs with L bol ≥ 5.1×10 45 erg/s at z = 4 in TNG300. This equates to an upper limit on the pair fraction (i.e. the number of BH pairs divided by the total number of BHs at that luminosity) of ∼ 3.1%. We are also able to reproduce a fairly redshift invariant pair fraction of ∼ 1% between z ∼ 4 and 1 (see Fig. 4) for BH pairs with a primary BH luminosity of L bol,pri ≥ 1×10 45.3 erg/s, when selecting systems as Silverman et al. (2020) do for the Horizon-AGN simulation (see their section 5.1). No such pairs were found in TNG300 above z = 4.
In conclusion, the observed number of bright QSO pairs seen at high redshift appears largely consistent with that predicted by current cosmological-sized galaxy evolution models. Although, further modelling efforts at the redshift of our observational sample are still required to confirm this.
3.4. The number of QSO pairs as a constraint on massive BH origins and assembly The origin of M BH ∼ 10 8−9 M black holes in massive galaxies at z ∼ 6 − 7 is largely unconstrained (e.g. Johnson et al. 2013;Valiante et al. 2017Valiante et al. , 2018Inayoshi et al. 2020). The cosmic time since the Big Bang is likely too short for light seeds with masses of ∼100 M (e.g., PopIII remnant model, Madau & Rees 2001;Bromm & Loeb 2003;Volonteri et al. 2003b) to grow to such masses even if the seeds accrete at the Eddington rate most of their lifetime. Instead, more massive seeds with masses of order ∼ 10 4 M would have a somewhat less challenging growth history to reach the masses of the high-z QSOs.
Several theoretical models have been put forward to explain the formation of such massive BH seeds. By analogy with the nuclear clusters that we commonly observe in low-redshift galaxies, compact nuclear clusters in metal-poor environments are predicted to form very massive stars by runaway stellar collisions (Omukai et al. 2008;Devecchi & Volonteri 2009;Regan & Haehnelt 2009). These very massive stars would collapse onto BH seeds of M BH ∼ 10 2 − 10 3 M . Seeds of about the same masses can be formed by runaway black hole mergers in metalpoor star clusters in the centre of high-redshift galaxies. The direct collapse model predicts the most massive seeds, called direct-collapse BHs (DCBHs): atomic cooling halos form single supermassive stars, which can collapse onto massive seeds of M BH ∼ 10 4 −10 6 M (Loeb & Rasio 1994;Bromm & Loeb 2003;Spaans & Silk 2006;Begelman et al. 2006;Lodato & Narayan 2006;Dijkstra et al. 2008;Visbal et al. 2014;Latif & Volonteri 2015;Habouzit et al. 2016b). In order to avoid the fragmentation of the gas into multiple and less massive stars, the presence or formation of efficient gas coolants (i.e., molecular hydrogen and metals) must be prevented. Therefore, the formation of the first DCBHs depends on the number of surrounding haloes as the Lyman-Werner radiation (photons with energy in the 11.2-13.6 eV range) needed to prevent H2 formation comes from surrounding star-forming galaxies (Habouzit et al. 2016a).
These mechanisms yield different galaxy occupation fractions, i.e. the probability for a galaxy to host a BH, and therefore could a priori lead to different predictions on the number of QSO pairs at high redshift. While the compact nuclear stellar cluster model and the model of runaway black hole mergers in clusters predict potentially large occupation fractions, i.e. the presence of BHs in many galaxies, the number density of direct collapse BHs is predicted to be small (e.g., Dijkstra et al. 2008;Habouzit et al. 2016b). However, DCBHs could appear clustered. Indeed, whether the large Lyman-Werner radiation is produced by several or a single star-forming halo, the high intensity required could irradiate several nearby halos, potentially forming DCBHs in some of these halos if metal-poor conditions are met. The need for clustered halo environments is not required for all the variations of the direct collapse model (e.g., the synchronised pair model, Visbal et al. 2014). It has also been proposed that DCBHs themselves could ignite a runaway process of further DCBH formation (Yue et al. 2014). Indeed, DCBHs are expected to be Compton thick, and their reprocessed radiation outshines small high-redshift galaxies by more than a factor 10. Thus, once the first DCBHs form, their Lyman-Werner radiation could help providing surrounding halos with the needed conditions to form new DCBHs (Yue et al. 2014). This process comes to an end when the host atomic-cooling halos are rapidly photo-evaporated by ionizing photons. An analytical model of these processes (Yue et al. 2014) shows that a highly clustered spatial distribution of haloes is the pre-requisite to triggering the formation of several DCBHs in the same region.
The different mechanisms of BH formation are not a priori mutually exclusive in the Universe, which therefore makes the interpretation of the QSO clustering difficult. Moreover, successive mergers of halos with time alter the initial halo occupation fraction and clustering of BHs at birth.
Observing several QSO pairs in our sample would likely suggest that massive BH seeds can form in the same regions (on scales of up to a Mpc), and that they can grow efficiently in the dense environments in which we often find high-redshift QSOs. The fact that we find no QSO pairs, i.e. no bright companion to any of our 116 spectroscopically confirmed QSOs could mean that (i) most QSOs are isolated objects with no massive galaxy in their surroundings, (ii) QSOs are often located in dense environments but nearby galaxies are devoid of BHs or devoid of massive BHs, or (iii) QSOs are often located in dense environments and nearby galaxies also host massive BHs, however these BHs have high variability/low duty cycle and no companion BHs were active when observed.
While the question of whether QSOs are embedded in the densest environments is still subject to debate in observations (Kim et al. 2009;Husband et al. 2013;Mazzucchelli et al. 2017a;McGreer et al. 2016;Habouzit et al. 2019, and references therein), the hypothesis (i) of a large fraction of QSOs being isolated objects is unlikely. Theoretical models suggest that the first massive galaxies formed through mergers of gasrich galaxies at very high redshift. Recently, Decarli et al. (2017) and Neeleman et al. (2019) identified six z > 6 QSOs with close, gas-rich companions through ALMA observations of [CII] and dust emission; 5 of these are in our sample: SDSSJ0842+1218, PSOJ167.6415-13.4960, SDSSJ1306+0356, PSOJ308.0416-21.2339 and CFHQSJ2100-1715. These QSOgalaxy pairs prove that z > 6 QSOs are not living isolated from the surrounding matter in their halos. This QSO-galaxy clustering is consistent with the quasar/LBG clustering at z ∼ 4, see Fig. 3 in Decarli et al. (2017). Recent deep X-ray observations of some of the above sources did not reveal potentially obscured quasars. So far, there is no strong evidence of these obscured companions of z ∼ 6 quasars to be obscured AGN (Connor et al. 2020;Vito et al. 2021).
Therefore, our non-detection of bright QSO-QSO pairs could hint towards an evolutionary scenario which prevents the formation of a massive accreting BH in the potentially gas-rich QSO companion. This could be due to tidal interaction of the primary QSO being more rapid than the accretion on the central BH of the companion. Feedback from the QSO could also reduce the gas reservoir of the companion galaxy through its lifetime, therefore diminishing the ability of a BH to efficiently accrete and become as massive as the BH powering the QSO. Additionally, efficient stripping of gas from companion galaxies of massive systems is seen out to several virial radii at low redshift in cosmological simulations (Ayromlou et al. 2021a,b). If this effect is present in the neighbourhood of the brightest QSOs at high redshift, it could also contribute to a reduction in the gas reservoir of companion BHs. The absence of QSO pairs could also be explained in case of low efficiency of BH formation, so that early massive BH formation happens only in a small fraction of galaxies, or alternatively if the formation of light seeds is dominant in the Universe. Due to their low masses, light seeds would have a hard time growing and even sinking efficiently to the potential well of their host galaxies. Indeed, searching for companions of the ten most massive BHs in TNG300 reveals none with masses M BH ≥ 10 7.7 M at z=5 and above. There is just one companion BH with a mass of M BH ∼ 10 7.3 M at z=6, but with a luminosity of only L bol = 10 43.3 erg/s, which is not even the brightest luminosity in the companion sample (which is 10 45.1 erg/s), and below our search sensitivity. With its limited volume of (205 h −1 cMpc) 3 there is no environment in TNG300 that was able to build up two extremely massive and luminous BHs in the same region by z ∼ 6.
Finally, this leads us to our hypothesis (iii) of very short coincident activity periods of both QSOs (see, e.g., Haiman & Hui (2001) for an early account of the lifetime-dependence of QSO clustering). In that case, QSO companion galaxies would also host massive BHs, but these BHs would not be efficiently accreting at the time of observation. In that case, the companion massive BHs would have grown by several orders of magnitude, potentially releasing a significant amount of energy in their host galaxies through cosmic time, which could explain an accretion history with a succession of "on and off" phases. In the observations presented here, we can only detect very bright objects. Further investigations at larger telescopes (e.g. with the GRONDdescendant SCORPIO instrument at Gemini Observatory) will shed new light onto the degenerate scenarios to explain the absence of QSO pairs with M 1450 (AB) < −26 mag.

Incidence of Brown dwarfs
The exposure for the spectroscopic observations was determined to provide >5σ detections of flux redwards of Ly-α (if it existed) binned at 200 Å. Thus, the signal-to-noise ratio in these spectra is good enough to distinguish brown dwarfs and lowredshift galaxies from QSOs, but way too small to perform a spectral classification of the dwarfs or galaxies. We therefore revert to the results of the Le PHARE fits, and provide the spectral types of the best-fitting template in the last columns of Tabs. 5, 6. For the conversion from the Le PHARE best-fit temperature to spectral type we used the "mean" value from the compilation of Kellogg (2017). The associated error is ±2 sub-types over most of the range.
Since after the selection of our targets for spectroscopic follow-up the unWISE (Schlafly et al. 2019) and CatWISE (Marocco et al. 2021) catalogs became available, we run a new cross-correlation of our targets against these catalogs. From our list of 17 candidates without ALLWISE counterparts (Tab. 6) we find 6 objects with W1 and W2 entries. Re-running Le PHARE fits with this extended photometry yields consistent results for three dwarf identifications (within the above mentioned ±2 subtype error), but deviates for the other three objects: (i) two SEDs, previously best fit with a QSO or M8 dwarf template are now best fit with a low-redshift (z ∼ 1.2 in both cases) galaxy, similar to ULASJ0148+0600 (Tab. 5), (ii) for another source, SDSSJ2054-0005_3, the best-fit template is that of a galaxy at redshift 6, which is inconsistent with our absolute photometry. The second-best template (nearly identical χ 2 ) is that of a QSO at redshift 6, which is unlikely given our GMOS spectrum, and the third-best template is that of a L2 dwarf, but with substantially worse χ 2 -we therefore do not assign an identification in Tab. 7.
For the candidate of PSOJ055.4244-00.8035, which also had a QSO template as best fit in Le PHARE but for which we did not obtain a spectrum, we refrain from changing the ID, though Occam's razor would argue against a QSO.
The combination of our spectroscopic results and the Le PHARE fits results in identifications as follows: 20 dwarfs, 3 low-redshift galaxies, 1 QSO, and 2 unsettled cases (Tabs. 5, 6, 7). Their spectral range is rather wide, from M5-T8, though accurate spectral typing remains to be done. The large incidence rate of such dwarfs is a surprise, both, because our search was not tuned towards late dwarfs, and also given our small search size of just 0.52 square degrees (summing over all 116 candidates). For instance, Kakazu et al. (2010) found 6 faint (19 ∼ < J ∼ < 20) ultracool T dwarfs over a 9.3 deg 2 area with a limiting magnitude of z (AB) ∼ < 23.3 mag. This is even slightly deeper than our mean depth of z (AB) ∼ < 23.1 mag, yet we find 2 T-dwarfs per 0.52 deg 2 , as compared to their 0.6 deg −2 , a factor 6 higher. Given that we pointed at known high-z QSOs, our search area similarly avoids the galactic plane. In a more recent compilation, Ryan & Reid (2016) used thick/thin exponential disk models and the luminosity function of ultracool dwarfs to predict their surface density for per spectral type. While this depends strongly on the galactic coordinates, the prediction for T-dwarfs has only a small dependence on galactic latitude, suggesting an even higher excess in our data, about a factor of 20, with even the Kakazu et al. (2010) surface density being higher by a factor of ∼3. For early L-type dwarfs, the excess is smaller, and for the M8-9 dwarfs we find consistency. Fig. 5 summarises the findings.
Since we observe a substantial excess of T dwarfs, search completeness cannot be an issue. Also foreground extinction correction should not be an issue, since most of these dwarfs are within of order 100 pc distance. Even if we adopt a twice as large (systematic?) error in the temperature estimates from the Le PHARE fits, this would not change our numbers. One potential reason could be the appropriateness of the templates in Le PHARE, i.e. that dwarfs are more numerous than e.g. lowredshift red galaxies. This is contrary to the case of high-z QSOs, for which late dwarfs form a numerous class of false positives. This finding of a larger than expected T dwarf number density is a surprising puzzle, which warrants further study. Existing GROND data could provide a first step in this direction, though spectroscopic typing of the candidate dwarfs is the ultimate step. QSOs, compared to the cumulative surface density as modeled by Ryan & Reid (2016), plotted here for their COSMOS sample field (central line) as well as the fields with the smallest and largest numbers (corresponding to different galactic latitudes) providing the covered range (shaded area). Our data points are plotted at the faintest magnitude per spectral type range, and the horizontal bar covers the magnitude range up to the brightest. While the observed number of M dwarfs is compatible with the known surface density, that of T dwarfs is substantially larger.   et al. (2010). The 4th column is the date at the start of the observing night, and the 5th column indicates the nth observing night of that target (first number) and which observation blocks (OB)s were used in the analysis (numbers after the underscore): more than one number indicates stacking of multiple OBs. Column six provides the 3σ limiting z magnitude of the (stacked) observation, not corrected for galactic foreground extinction (which is given in column 7, based on (Schlafly & Finkbeiner 2011)). The letter attached to the z limiting magnitude denotes the catalog used for photometric calibration of the optical channels g r i z : S=SDSS (DR12), P=PanSTARRS1, SM=SkyMapper (DR1.1).     Table 6. Pair candidates without ALLWISE photometry, not corrected for foreground-A V , all in the AB system. The ID column gives the bestfitting Le PHARE spectral template. The objects marked in boldface are those where the Le PHARE fit preferred a QSO high-z template, and we obtained optical spectroscopy. Except for PSOJ071.0322-04.5591 our spectra argue against a QSO nature, and we therefore provide the next best-fitting template, mostly a dwarf star (see text). For objects with a * in the last column see Tab

Conclusions
No candidate QSO pair has been found in a search of 116 spectroscopically confirmed redshift 6 quasars within our 0.1-3.3 h −1 cMpc search radius with a companion brighter than M 1450 (AB)< −26 mag. Although this result is consistent with cosmologicalscale galaxy evolution simulations (given the redshift, search volume, and luminosity limits), the serendipitously discovered pairs at z=4.26 (Schneider et al. 2000) and z=5.02 ) suggested one pair every 40 QSOs, or 2.5%, implying an empirical expectation of 3 pairs in our sample of 116 search targets. At first glance, this result suggests little clustering at redshift 6. However, with updated knowledge of the faint end of the luminosity function, the w p statistics returns more than an order of magnitude higher value than inferred from the low-redshift clustered QSOs; thus, given the still low statistics in our search and the uncertainties (in Table 7. Pair candidates from Tab. 6 with unWISE or CatWISE photometry (AB system). The ID columns give the best-fitting Le PHARE spectral template without (ID1) and with un/Cal-WISE photometry (ID2 particular in the unknown search volumes) of the serendipitous detections, we are not able to claim that the clustering at redshift 6 is indeed smaller than up to redshift 5. As a consequence of our search, we do find a higher than expected rate of ultracool dwarfs. While the excess of L dwarfs is marginally consistent with statistics and could be due to the accidental lack of a single source about 1 mag fainter, this argument does not hold for the T dwarfs. The observed excess of a factor of ∼20 suggests either (i) strong surface density variations over the sky, (ii) a lack of appropriate templates for other source types (e.g. low-redshift red galaxies) such that the T templates just happen to be the closest in shape to the observed SEDs, or (iii) an under-estimate of the true rate of T dwarfs. A larger survey at our depth would be needed to settle this question. Fig. 6. Le PHARE fits for the 10 objects from Tabs. 5 & 6, for which a QSO template provides the better fit, plus the one with a galaxy template fitting best; from top left to right bottom: ULASJ0148+0600, CFHQSJ1509-1749, PSOJ009.7355-10.4316, PSOJ011.3899+09.0325_C2, PSOJ045.1840-22.5408, PSOJ071.0322-04.5591, PSOJ261.0364+19.0286_C1, PSOJ267.0021+22.7812_C2, SDSSJ2054-0005_C3, SDSSJ2054-0005_C4, IMSJ2204+0012. These are the SEDs after correction of the galactic foreground-A V , and most of the dwarf fits return an even worse reduced χ 2 for the no-A V SEDs than shown here.