We expect luminous M SUB1450SUB  26.5 highredshift quasars to trace the highestdensity peaks in the early Universe. Here we present observations of four z  6 quasar fields using JWSTNIRCam in the imaging and widefield slitless spectroscopy mode and report a wide range in the number of detected O IIIemitting galaxies in the quasars environments ranging between a density enhancement of   65 within a 2 cMpc radiusone of the largest protoclusters during the Epoch of Reionization discovered to dateto a density contrast consistent with zero indicating the presence of a UVluminous quasar in a region comparable to the average density of the Universe. By measuring the twopoint crosscorrelation function of quasars and their surrounding galaxies as well as the galaxy autocorrelation function we infer a correlation length of quasars at ltzgt  6.25 of inlineformula mmlmath overflowscrollmmlmsubsupmmlmrowmmlmirmmlmimmlmrowmmlmrowmmlmn0mmlmnmmlmrowmmlmrowmmlmiQQmmlmimmlmrowmmlmsubsupmmlmommlmommlmsubsupmmlmrowmmlmn22.0mmlmnmmlmrowmmlmrowmmlmommlmommlmn2.9mmlmnmmlmrowmmlmrowmmlmommlmommlmn3.0mmlmnmmlmrowmmlmsubsupmmlmspace width0.25emmmlmspacemmlmicMpcmmlmimmlmspace width0.25emmmlmspacemmlmsupmmlmrowmmlmihmmlmimmlmrowmmlmrowmmlmommlmommlmn1mmlmnmmlmrowmmlmsupmmlmath inlineformula while we obtain a correlation length of the O IIIemitting galaxies of inlineformula mmlmath overflowscrollmmlmsubsupmmlmrowmmlmirmmlmimmlmrowmmlmrowmmlmn0mmlmnmmlmrowmmlmrowmmlmiGGmmlmimmlmrowmmlmsubsupmmlmommlmommlmn4.1mmlmnmmlmspace width0.25emmmlmspacemmlmommlmommlmn0.3mmlmnmmlmspace width0.25emmmlmspacemmlmicMpcmmlmimmlmspace width0.25emmmlmspacemmlmsupmmlmrowmmlmihmmlmimmlmrowmmlmrowmmlmommlmommlmn1mmlmnmmlmrowmmlmsupmmlmath inlineformula. By comparing the correlation functions to darkmatteronly simulations we estimate the minimum mass of the quasars host dark matter halos to be inlineformula mmlmath overflowscrollmmlmsubmmlmrowmmlmilogmmlmimmlmrowmmlmrowmmlmn10mmlmnmmlmrowmmlmsubmmlmo stretchyfalsemmlmommlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmihalommlmimmlmommlmommlmiminmmlmimmlmrowmmlmsubmmlmrowmmlmo stretchytruemmlmommlmrowmmlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmsubmmlmo stretchyfalsemmlmommlmommlmommlmsubsupmmlmrowmmlmn12.43mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.15mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.13mmlmnmmlmrowmmlmsubsupmmlmath inlineformula and inlineformula mmlmath overflowscrollmmlmsubmmlmrowmmlmilogmmlmimmlmrowmmlmrowmmlmn10mmlmnmmlmrowmmlmsubmmlmo stretchyfalsemmlmommlmsubsupmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmihalommlmimmlmommlmommlmiminmmlmimmlmrowmmlmrowmmlmo stretchyfalsemmlmommlmiOIIImmlmimmlmo stretchyfalsemmlmommlmrowmmlmsubsupmmlmrowmmlmo stretchytruemmlmommlmrowmmlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmsubmmlmo stretchyfalsemmlmommlmspace width0.25emmmlmspacemmlmommlmommlmsubsupmmlmrowmmlmn10.56mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.03mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.05mmlmnmmlmrowmmlmsubsupmmlmath inlineformula for the O III emitters indicating that a luminous quasars do not necessarily reside within the most overdense regions in the early Universe and that b the UVluminous duty cycle of quasar activity at these redshifts is f SUBdutySUB  1. Such short quasar activity timescales challenge our understanding of early supermassive black hole growth and provide evidence for highly dustobscured growth phases or episodic radiatively inefficient accretion rates.
['Quasars', 'Supermassive black holes', 'Clustering', 'Astrostatistics techniques', 'Galaxy dark matter halos', 'High-redshift galaxy clusters', 'Large-scale structure of the universe', '1319', '1663', '1908', '1886', '1880', '2007', '902', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
EIGER. VI. The Correlation Function Host Halo Mass and Duty Cycle of Luminous Quasars at z  6
{'EIGER VI. The Correlation Function, Host Halo Mass and Duty Cycle of Luminous Quasars at z ≳ 6': 'Anna-Christina Eilers, 1, 2 Ruari Mackenzie, 3 Elia Pizzati, 4 Jorryt Matthee, 3, 5 Joseph F. Hennawi, 4, 6 Haowen Zhang, 7 Rongmon Bordoloi, 8 Daichi Kashino, 9 Simon J. Lilly, 3 Rohan P. Naidu, 2, ∗ Robert A. Simcoe, 1, 2 Minghao Yue, 2 Carlos S. Frenk, 10 John C. Helly, 10 Matthieu Schaller, 4, 11 and Joop Schaye 4  \n1  \nDepartment of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 2 MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 3 Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, Zurich, 8093, Switzerland 4 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands 5 Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria 6 Department of Physics, University of California, Santa Barbara, CA 93106-9530, USA 7 Department of Astronomy, University of Arizona, Tucson, AZ 85721, USA 8 Department of Physics, North Carolina State University, Raleigh, 27695, North Carolina, USA 9 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 10 Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham, DH1 3LE, UK 11 Lorentz Institute for Theoretical Physics, Leiden University, PO Box 9506, NL-2300 RA Leiden, The Netherlands  \n(Revised September 5, 2024)  \nSubmitted to ApJ', 'ABSTRACT': "We expect luminous ( M 1450 ≲ -26 . 5) high-redshift quasars to trace the highest density peaks in the early universe. Here, we present observations of four z ≳ 6 quasar fields using JWST/NIRCam in imaging and widefield slitless spectroscopy mode and report a wide range in the number of detected [O III ]-emitting galaxies in the quasars' environments, ranging between a density enhancement of δ ≈ 65 within a 2 cMpc radius - one of the largest proto-clusters during the Epoch of Reionization discovered to date - to a density contrast consistent with zero, indicating the presence of a UV-luminous quasar in a region comparable to the average density of the universe. By measuring the two-point cross-correlation function of quasars and their surrounding galaxies, as well as the galaxy auto-correlation function, we infer a correlation length of quasars at ⟨ z ⟩ = 6 . 25 of r QQ 0 = 22 . 0 +3 . 0 -2 . 9 cMpc h -1 , while we obtain a correlation length of the [O III ]-emitting galaxies of r GG 0 = 4 . 1 ± 0 . 3 cMpc h -1 . By comparing the correlation functions to dark-matter-only simulations we estimate the minimum mass of the quasars' host dark matter halos to be log 10 ( M halo , min /M ⊙ ) = 12 . 43 +0 . 13 -0 . 15 (and log 10 ( M [OIII] halo , min /M ⊙ ) = 10 . 56 +0 . 05 -0 . 03 for the [O III ]-emitters), indicating that (a) luminous quasars do not necessarily reside within the most overdense regions in the early universe, and that (b) the UV-luminous duty cycle of quasar activity at these redshifts is f duty ≪ 1. Such short quasar activity timescales challenge our understanding of early supermassive black hole growth and provide evidence for highly dust-obscured growth phases or episodic, radiatively inefficient accretion rates.  \nKeywords: dark ages, early universe - quasars: supermassive black holes - methods: data analysis", '1. INTRODUCTION': "Corresponding author: Anna-Christina Eilers eilers@mit.edu  \n∗ NHFP Hubble Fellow  \nThe existence of supermassive black holes (SMBHs) with ∼ 10 9 -10 10 M ⊙ harbored in the center of luminous quasars at redshifts z ≳ 6, when the universe was less than a billion years old, poses significant challenges  \nto our understanding of black hole growth (e.g. Mazzucchelli et al. 2017b; Ba˜nados et al. 2018; Yang et al. 2020; Wang et al. 2021; Yue et al. 2023; Bigwood et al. 2024). In the standard picture of SMBH growth, the rest mass energy of the accreted matter is divided between radiation and black hole growth via the so-called radiative efficiency ϵ , implying that the emission of quasar light is concurrent with the growth of the black hole. Every massive galaxy is thought to have undergone a luminous quasar phase, and thus high-redshift quasars represent the progenitors of the dormant SMBH population found in the centers of nearly all local bulge-dominated galaxies (Small & Blandford 1992; Yu & Tremaine 2002). This fundamental link between high-redshift quasars and local galaxies is supported by the tight correlations between the masses of SMBHs and various galaxy properties, such as the stellar velocity dispersion or bulge mass (e.g. Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese & Merritt 2000; Haring & Rix 2004).  \nIn the current paradigm, black holes grow exponentially starting from an initial black hole seed M seed during the lifetime t Q of the quasar, which is the time that galaxies spent in the luminous quasar state, i.e.  \nM BH ( t Q ) = M seed exp( t Q /t S ) . (1)  \nThe growth occurs on a characteristic timescale t S known as the e -folding time, or 'Salpeter' time (Salpeter 1964), which is approximately  \nt S ≈ 0 . 45 ( ϵ 1 -ϵ )( L L Edd ) -1 Gyr , (2)  \nwhere L and L Edd denote the bolometric and Eddington luminosity of the quasar, respectively. Assuming a thin accretion disk around the black hole, the radiative efficiency ϵ ranges between ∼ 6% -34%, the exact value depending on the spin of the black hole (Shakura & Sunyaev 1973; Thorne 1974).  \nThus, long timescales of quasar activity - comparable to the Hubble time t H - are required with nearly continuous Eddington-limited accretion, in order to explain the growth of the observed SMBHs at high redshift (e.g. Tanaka & Haiman 2009; Begelman & Volonteri 2017). It follows that the fraction of cosmic time that galaxies shine as luminous quasars, known as the duty cycle f duty = t Q /t H ( z ), is expected to be approximately unity for quasars in the early universe.  \nThis theoretical picture of SMBH formation predicts that early SMBHs form and grow in the centers of rare and massive dark matter halos, which constitute the highest density peaks of the cosmological density field (Efstathiou & Rees 1988; Cole & Kaiser 1989; Nusser & Silk 1993; Djorgovski et al. 2003; Volonteri & Rees  \n2006; Sijacki et al. 2009; Costa et al. 2014; Trebitsch et al. 2019; Costa 2023). Thus, luminous quasars are expected to reside in massive overdensities at early cosmic times, which implies that one should find a high number of companion galaxies in their vicinity according to our ΛCDM cosmological model (e.g. Mu˜noz & Loeb 2008; Tinker et al. 2010). The regions around some high-redshift quasars are predicted to constitute protoclusters, progenitors of the galaxy clusters that populate the local universe, albeit the growth of these structures across cosmic time sensitively depends on a variety of environmental factors (e.g. Angulo et al. 2012).  \nHowever, to date the search for the expected abundance of companion galaxies around z ≳ 6 quasars has led to inconclusive results. While some studies find strong galaxy enhancements in quasar fields using mostly narrow-band filters to capture the Ly α emission of galaxies in the quasars' environments (e.g. Kim et al. 2009; Morselli et al. 2014), others do not (e.g. Ba˜nados et al. 2013; Simpson et al. 2014; Mazzucchelli et al. 2017a). The first, recently published observations with JWST/NIRCam WFSS of two high-redshift quasar fields revealed an abundance of [O III ]-emitting galaxies in the quasars' environment, indicating that these particular quasars lie indeed in highly overdense regions in the universe (Wang et al. 2023; Kashino et al. 2023). It is unclear whether the conflicting results arise from the limited sensitivity, velocity offsets between the narrowband filter and the quasars' redshifts, or a too small field of view of some of the observations (e.g. Overzier 2016), or alternatively, at least some quasars may inhabit less massive hosts at z ≳ 6 than expected.  \nThe mass of the dark matter halos hosting these luminous quasars can be inferred by means of their clustering properties. The two-point correlation function is a powerful tool to probe the connection between the observed light and host dark matter halos (e.g. Osmer 1981; Shanks & Boyle 1994; Croom et al. 2001, 2005; Shen et al. 2007). At low and intermediate redshifts of 1 ≲ z ≲ 4 the two-point auto-correlation function has been studied using several thousands of quasars distributed over large areas on the sky. Across this wide redshift range most quasar clustering studies find similar masses of the dark matter halos that host luminous quasars of a few times 10 12 M ⊙ (e.g. Porciani et al. 2004; Croom et al. 2005; Shen et al. 2007; Myers et al. 2007; Coil et al. 2007; Ross et al. 2009; Padmanabhan et al. 2009; Trainor & Steidel 2013; Eftekharzadeh et al. 2015; He et al. 2018). This implies that quasars at intermediate redshifts are highly biased tracers of the underlying density field and reside within some of the most massive dark matter halos at these cosmic times. At z ∼ 1, how-  \ne dark matter halo mass is close to the characteristic mass in the Press-Schechter mass function (Press & Schechter 1974), indicating that quasars in the local universe are largely unbiased (Croom et al. 2005). Thus, we observe an overall increase in the quasars' bias factor with redshift, i.e. the ratio of the host halo mass for a typical quasar to the mean halo mass at the same epoch rises with increasing redshift.  \nIn turn, this implies an increase in the quasars' duty cycle from ≲ 1% at low redshifts, i.e. z ≲ 1 (Croom et al. 2005; Padmanabhan et al. 2009; Laurent et al. 2017), to a few tens of percent at z ∼ 4 (Shen et al. (2007); Pizzati et al. (2024a), but see Eftekharzadeh et al. (2015) who report a lower duty cycle at similar redshifts). Thus, the quasar activity timescales at these redshifts are expected to be relatively long, i.e. t Q > 10 8 years, therefore providing sufficient time to grow the observed SMBHs.  \nAt redshifts z > 6, however, this approach to determine the host dark matter halo masses by measuring the quasars' auto-correlation function becomes unfeasible as the space density of luminous ( M 1450 < -27) quasars drops dramatically to less than 1 Gpc -3 (Wang et al. 2019; Schindler et al. 2023). Using a much fainter ( M 1450 ≳ -25), but more numerous population of quasars a first attempt of the auto-correlation measurement at z ∼ 6 has recently been reported by the Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs) Collaboration (Arita et al. 2023), indicating host dark matter halo masses similar or slightly higher than those found at lower redshifts, but uncertainties are large. However, while auto-correlation measurements of quasars at these high redshifts remain challenging, the complementary quasar-galaxy cross-correlation measurement is still viable. Assuming that both quasars and galaxies trace the same underlying dark matter density distribution, their respective auto-correlation functions determine the cross-correlation function between these two classes of objects.  \nDuring the last year, JWST has revolutionized our ability to study the early universe, enabling us to reveal and spectroscopically confirm faint galaxies down to m AB ∼ 29 mag within the environment of high-redshift quasars. Using NIRCam in imaging and WFSS mode, JWST recently obtained large mosaic observations in five quasar fields at 6 ≲ z ≲ 7 as part of the Emission line galaxies and Intergalactic Gas in the Epoch of Reionization (EIGER) Collaboration (Kashino et al. 2023). The filter for the grism spectra was tuned to spectroscopically identify galaxies in the quasar fields by means of their [O III ] emission line doublet λλ 4960 , 5008. Using this data set, we present in this paper for the first time the quasar-galaxy cross-correlation measurement  \nat z ≳ 6. Combined with the auto-correlation of the [O III ]-emitting galaxies, we infer the clustering properties of these high-redshift luminous quasars, and determine their host dark matter halo mass as well as their duty cycle.  \nIn § 2 we will first describe the JWST/NIRCam observations of the high-redshift quasar fields analyzed in this study. We will briefly comment in § 3 on the large diversity in the observed abundance of galaxies that we observe in the environments of the quasars, before measuring the two-point cross-correlation function of the quasars and their surrounding galaxies, as well as the galaxy-galaxy auto-correlation function in § 4. In § 5 we infer the characteristic mass of the quasars' host dark matter halos based on their clustering properties, and estimate their duty cycle in § 6. We discuss the implications of our results in § 7 before summarizing in § 8.  \nThroughout this work we adapt a flat ΛCDM cosmology with H 0 = 67 . 66 kms -1 Mpc -1 , Ω M = 0 . 3097, and Ω Λ = 0 . 6889 (Planck Collaboration et al. 2018).", '2. JWST/NIRCAM OBSERVATIONS OF HIGH-REDSHIFT QUASAR FIELDS': "The EIGER Collaboration is observing a total of six z ≳ 6 quasar fields utilizing NIRCam in WFSS and imaging mode (ID: #1243). Of these six fields, the observations of five quasar fields were completed at the time of writing, and their analysis is presented in this work. Four out of five fields are used to study the quasargalaxy clustering (as in one field our setup does not cover the [O III ]-emission at the quasar's redshift), while all five quasar fields are used when measuring the galaxy auto-correlation (see § 4).", '2.1. Data Reduction': "An overview of the observing program as well as a detailed description of the data reduction can be found in Kashino et al. (2023) and Matthee et al. (2023a). We will briefly summarize the setup of the program and the main steps of the data reduction here.  \nThe NIRCam imaging and WFSS observations of the five quasar fields were conducted between August 22, 2022 and June 2, 2023. The observations consist of four visits per quasar field forming a 2 × 2 overlapping mosaic pattern, which provides an area of approximately 3 ' × 6 ' per field centered on the quasar. The central 40 '' × 40 '' region around the quasar is covered by all four mosaic tiles.  \nThe WFSS observations use the grism R in the long wavelengths (LW) channel combined with the F356W filter. This yields spectra covering the wavelengths between 3 . 1-4 . 0 µ mdispersed along the detector rows with  \n<!-- image -->  \nJ1148+5251  \n<!-- image -->  \nJ1030+0524  \n<!-- image -->  \n<!-- image -->  \nredshift 2  \nFigure 1. Observations of four z ≳ 6 quasar fields. Top panels show the RGB images constructed from NIRCam imaging data (F115W, F200W and F356W) of each quasar field. Bottom panels show the optical and NIR spectra of the quasars observed with X-Shooter/VLT and FIRE/Magellan (J0100+2802, J0148+0600, J1030+0524), or MOSFIRE+ESI/Keck (J1148+5251). The colored data points show the redshifts and projected distances of all [O III ]-emitting galaxies detected in the NIRCam grism spectra above a luminosity of L [O III] , 5008 = 10 42 erg s -1 ; their color corresponds to their velocity offset with respect to the quasars' redshifts and their size correlates with the logarithm of their luminosity, log 10 L [O III] , 5008 . The grey shaded area denotes the velocity interval of | ∆ v | ≤ v max around the redshift of the quasar. All galaxies that fall within this velocity range are also circled in the respective RGB images to show the spatial distribution around the quasar (indicated by a red arrow).  \n<!-- image -->  \na spectral resolution of R ≈ 1500. The two NIRCam modules A and B disperse the spectra in opposite directions. These reversed grism spectra facilitate significantly the identification of the source objects in the  \ndirect images, and enable us to remove contaminating lines from other sources.  \nOur observations employ a three-point primary dither using INTRAMODULEX and a four-point subpixel dither,  \n<!-- image -->  \n<!-- image -->  \nwhich yields a total of 12 exposures per visit and module for the direct images in the two short wavelength (SW) filters (F115W and F200W), and 24 WFSS images. By co-adding the multiple exposures, the total exposure time in a single visit amounts to 4380 s per filter for the SW imaging and 8760 s for WFSS. We obtained additional direct imaging in the LW channel F356W filter after the WFSS observations, with F200W in the SW channel. Given the small-scale dithering and larger-scale mosaic pattern, the total exposure time varies across the field, ranging from 1 . 5 ks at the edges sampled by only a single primary dither up to 13 ks in the center of each field for SW imaging, and from 2 . 9 ks up to 35 ks for WFSS (Kashino et al. 2023).  \nThe WFSS data were reduced with a combination of the jwst pipeline (version 1.9.4) and additional custom processing steps. We first process the raw exposure files with the Detector1 step from the pipeline and use the Spec2 step to assign the astrometric solution for each image. The data is flat-fielded by applying the Image2 step and sky background variations as well as 1 /f noise are removed using the median counts in each column. We construct an emission line image by applying a running median filter along each row to remove any continuum emission (see Kashino et al. 2023, for details).", '2.2. Identification of [O III ]-Emitting Galaxies': 'The spectral range covered by the filter F356W of 3 . 1 ≤ λ ≤ 4 . 0 µ m enables us to identify [O III ]-emitters within the redshift range 5 . 33 ≲ z ≲ 6 . 96. This implies that we can study the [O III ]-emitting galaxies in the environment of four of the targeted quasars, while for one quasar, J1120+0641 at z = 7 . 0848, the [O III ] emission lines lie outside of the observed redshift range. The galaxies observed in this field, however, will be included in the measurement of the galaxy-galaxy autocorrelation function ( § 4.2).  \nWe employ two complementary methods to identify [O III ]-emitting galaxies from the emission line images. The first approach starts with the co-added emission line image and searches for combinations of emission lines that could plausibly arise from doublets ([O III ] λλ 4960 , 5008) or triplets ([O III ] λλ 4960 , 5008 + H β λ 4861) using SExtractor . To evaluate the validity of these detections, we then try to identify a source galaxy in the F356W direct image using the grism model, which enables us to estimate the location on the sky for a given position in the grism image.  \nAs a second method to identify [O III ]-emitters, we employ an inverse approach and use the catalog of sources identified in the F356W images as a starting point. For each source we extract a 2D spectrum from  \nthe median-filtered grism images and run SExtractor to identify emission lines. [O III ]-emitters are identified if emission line doublets or triplets are detected at ≥ 3 σ close to the expected locus of the spectral trace of the object. The list of detected [O III ]-emitters from both approaches are afterwards combined and cross-validated. We estimate the systematic uncertainty on the derived redshifts to be about 100 km s -1 due to the inaccuracy of the grism model (Kashino et al. 2023).  \nNote that our identification algorithm for [O III ]-emitters has been fine-tuned and improved since the publication of the first results from the EIGER Collaboration (Kashino et al. 2023; Matthee et al. 2023a). This results in an increased number of discovered [O III ]-emitting galaxies. Within the field around quasar J0100+2802 we now discover a total of 180 [O III ]-emitting galaxy systems (213 individual galaxies, see § 2.3 for details on how we define a group) compared to the previously reported 117. Our total sample of individual [O III ]-emitters in all five quasar fields amounts to 861 (Kashino et al. in prep.).', '2.3. Groups of [O III ]-Emitting Galaxies': "When analyzing the distribution of [O III ]-emitting galaxy pairs as a function of their angular separation, we find a significant excess in the number of 'clumps', i.e. pairs or multiples, on very small scales at separations ≲ 2 '' , which corresponds to approximately 12 kpc at z = 6 (see Fig. 4 in Matthee et al. 2023a). All of these [O III ]-emitting systems within 2 '' separation are similarly closely separated in redshift ( | ∆ v | ≲ 600 kms -1 ), suggesting that these components are physically associated with each other. Consequently, we merge all [O III ]-emitters within an angular separation < 2 '' and | ∆ v | < 1000 kms -1 together into a single system. The motivation for this is to mitigate uncertainties in the identification of individual objects that can arise due to the choice of deblending parameters used in SExtractor .  \nThis merges about ∼ 13% of the detected indivdiual [O III ]-emitters and results in a total of 750 [O III ]-emitting systems discovered in the five quasar fields. Since the completeness of our detections drops significantly below a luminosity limit of L [OIII] , 5008 < 10 42 erg s -1 (see § 2.4), we restrict ourselves in this study to the 533 systems more luminous than this threshold (see Tab. 1). Furthermore, we restrict our analysis to systems that have been found using the second method to identify [O III ]-emitters starting from the source catalog identified in the F356W image (see § 2.2), since the completeness of our observations is assessed using this method (see § 2.4). This step only removes 11 objects, with a total of 522 systems remaining for our analysis.  \nTable 1. Overview of the number of [O III ]-emitting galaxies in the quasar fields.  \nNote -The columns denote the quasar field, the quasars' coordinates RA and DEC, the quasars' redshifts (determined by means of a : [C II ] or b : Mg II emission lines) and absolute magnitude at 1450 ˚ A in the rest-frame, the masses of the quasars' SMBHs determined via the H β emission line (Yue et al. 2023) without any systematic uncertainties of ∼ 0 . 5 dex, the number of [O III ]-emitting galaxy systems above a luminosity limit of L [O III] , 5008 = 10 42 erg s -1 detected in the whole quasar field, the number of these galaxies within | ∆ v | ≤ 1000 km s -1 of the quasar's redshift, the total observed survey volume V around the quasar within ± 1000 km s -1 , and the observed overdensity δ within a radius of 2 cMpc around the quasar.  \nFig. 1 shows our observations for four of the quasar fields, for which the quasars' environment can be studied with our setup. For each field, we show in the top panels the NIRCam images (F115W, F200W and F356W) indicating the location of the quasar in the center as well as all [O III ]-emitting galaxies within the quasar's environment, i.e. within | ∆ v | ≤ 1000 kms -1 from the quasar's systemic redshift. The bottom panels show the quasar spectra observed with the ground-based spectrographs FIRE (Simcoe et al. 2008) on the Magellan Telescopes, X-Shooter (Vernet et al. 2011) on the Very Large Telescope and MOSFIRE (McLean et al. 2012) and ESI (Sheinis et al. 2002) on the Keck Telescopes (see Eilers et al. 2023; ˇ Durovˇc'ıkov'a et al. 2024, for details on the ground-based quasar spectra), as well as all detected [O III ]-emitting systems in the quasar field as a function of redshift and projected distance from the quasar.", '2.4. Modeling the completeness of our observations and creating random galaxy catalogs': "In order to measure the clustering in the quasar fields as well as the excess in galaxies in the quasars' environments compared to a blank field, we need to construct a set of 'random' galaxy catalogs. These catalogs contain a large number of sources with similar properties as the observed galaxies, but randomly distributed across the survey volume.  \nOur selection function is not trivial due to the 4point mosaic, the different sensitivities of NIRCam's module A and B, as well as the wavelength dependence of the noise over our field of view. To this end, we inject random sources in our survey volume and forward model their spectra through our survey coverage and completeness. Their [O III ]-luminosities are randomly  \ndrawn from the [O III ]-luminosity distribution function (Matthee et al. 2023a, with updated model parameters from Mackenzie et al. in prep.) and a minimum luminosity of L [OIII] , 5008 = 10 42 erg s -1 in the [O III ]-emission line at 5008 ˚ A. Note that the [O III ]-luminosity function is estimated over the entire redshift range, but masking all galaxies at the quasars' redshifts, such that no biases due to potential overdensities in the quasars' environments arise. The spatial positions and redshifts are randomly chosen within the observed survey volume. We then check the spatial and spectral coverage of this source in our survey, and ensure that both emission lines of the [O III ]-doublet, i.e. at 4960 ˚ A and 5008 ˚ A, would have been observed in our setup. Since both our methods to identify [O III ]-emitting galaxies are based on the observability of the line doublet (see § 2.2), we would miss all sources in our catalog for which only a single line would be observed.  \nTo model the completeness of our line-selected sample into the random catalogs, we apply the completeness models that will be described in a forthcoming companion paper (Mackenzie et al. in prep.). In short, each field has an individual model, which consists of a noise cube and a completeness interpolating function. The noise cube represents the typical propagated error in a 2D stacked spectrum as a function of position and wavelength. The function is constructed by injecting millions of mock sources into the 2D spectra and re-running the [O III ]-emitter identification step. Thus, the measured completeness is not a step function at some limiting flux, but instead is a smooth function which we interpolate from the injection experiments. This function then translates the flux of the weaker of the two  \nFigure 2. Density contrast of [O III ]-emitters around the four z ≳ 6 quasars analyzed in this work. The overdensity is measured in cylinders around the quasars with radius r and a depth of | ∆ v | = 1000kms -1 around the quasars' redshifts.  \n<!-- image -->  \n[O III ] emission lines, i.e. 4960 ˚ A, divided by the local error, to a completeness fraction. Dividing by the error essentially re-scales the flux such that the completeness function is independent of the local sensitivity or redshift. To create the random galaxy catalogs we calculate these completeness maps of each EIGER field, and then draw random samples using the completeness as a probability threshold. Since our observations become highly incomplete for sources below a luminosity of L [OIII] , 5008 = 10 42 erg s -1 , we only consider galaxies above this luminosity limit for our analysis, which corresponds to a line flux of for the weaker [O III ]-emission line of f [OIII] , 4960 = 0 . 7 × 10 -18 erg s -1 cm -2 at z ≈ 6 . 3.  \nWe estimate an average completeness of ∼ 41% for galaxies with a luminosity of L [OIII] , 5008 = 10 42 erg s -1 at 6 . 0 < z < 6 . 4, which increases to ∼ 98% for L [OIII] , 5008 = 10 43 erg s -1 . Note that these completeness estimates are averaged over the five quasar fields and the full spatial coverage of each mosaic, for regions closer to the quasar the completeness is higher. For more details on the completeness evaluation of our observations we refer the reader to (Matthee et al. 2023a, Mackenzie et al. in prep.).", '3. LARGE DIVERSITY IN QUASAR ENVIRONMENTS': "Current SMBH formation models predict that luminous quasars in the high-redshift universe trace the highest-density peaks and thus reside in a significantly overdense environment (e.g. Sijacki et al. 2009). To this end we estimate the density contrast, i.e. δ = N gal /N exp -1, by comparing the number of galaxies detected in close vicinity to the quasars N gal , to  \nthe number of galaxies one would expect within the same observed volume N exp in a blank 'random' part of the universe. In this case, N gal denotes simply the number counts of [O III ]-emitting galaxy systems with | ∆ v | ≤ 1000 kms -1 from the quasars' systemic redshift and within different radii. The expected number of objects N exp is estimated by multiplying the galaxy density n gal inferred from the [O III ]-luminosity function (Matthee et al. 2023a, Mackenzie et al. in prep.) by the completeness map derived in § 2.4 and by the observed volume V around the quasar, i.e. a cylinder with a given radius and depth of ± 1000 kms -1 . We estimate the uncertainties on the density contrast from Poisson noise on the galaxy counts.  \nOur results are shown as a function of radius r for the four quasar fields in Fig. 2 and suggest a strong diversity among the observed quasar fields. We find significant overdensities around three of the four observed quasars. In close vicinity of the quasars within r < 2 cMpc we find density enhancements of δ ( r < 2 cMpc) = 29 ± 10 and δ ( r < 2 cMpc) = 16 ± 8 around the quasars J0100+2802 and J1148+5251, respectively. With a total of N gal = 47 galaxies observed in the complete observed volume around the quasar J0148+0600, this quasar resides in a spectacular proto-cluster within the Epoch of Reionization with an overdensity of δ ( r < 2 cMpc) = 65 ± 15.  \nInterestingly, one of the targeted quasar fields around J1030+0542, shows only 2 [O III ]-emitting galaxies within the | ∆ v | ≤ 1000 kms -1 velocity interval from the quasar and therefore no galaxy enhancement compared to a blank field, i.e. δ ( r < 2 cMpc) = 3 ± 4. We will discuss the possible implications of the highly diverse quasar environments in § 7. All measurements and galaxy counts across the whole observed volume, as well as the density enhancements at r < 2 cMpc, are listed in Tab. 1.", '4. THE CLUSTERING PROPERTIES OF HIGH-REDSHIFT QUASARS': "In this section we present measurements of the crosscorrelation function between the quasars and the [O III ]-emitting galaxies in their environment ( § 4.1), as well as of the auto-correlation function of [O III ]-emitters outside of the quasars' environments ( § 4.2). By jointly fitting these two correlation functions we will infer the clustering properties of the luminous quasars themselves ( § 4.3 and § 4.4).", '4.1. The Quasar-Galaxy Cross-Correlation Function': "We first measure the volume-averaged projected cross-correlation function between the quasars and the surrounding [O III ]-emitting galaxies, i.e.  \nTable 2. Quasar-galaxy cross-correlation and galaxy-galaxy auto-correlation measurements at ⟨ z ⟩ = 6 . 25.  \nχ QG ( R min , R max ), which is a dimensionless quantity defined as the real-space quasar-galaxy cross-correlation function ξ QG ( R,Z ) integrated over and normalized by a comoving volume V (e.g. Hennawi et al. 2006; Garc'ıa-Vergara et al. 2017), i.e.  \nχ QG ( R min , R max ) = 2 V ∫ R max R min ∫ Z max 0 ξ QG ( R,Z ) 2 πR d R d Z. (3)  \nWe consider galaxies with a maximum velocity difference from the central quasar of | ∆ v | ≤ v max with v max = 1000kms -1 , and thus Z max = v max (1+ z ) /H ( z ), where H ( z ) denotes the Hubble parameter at redshift z . The integrated volume V is a cylindrical shell with inner and outer radii R min and R max , respectively, as well as a depth of Z max . We integrate over the redshift range, in order to avoid the need to model redshift-space distortions and redshift uncertainties.  \nFollowing previous work (e.g. Garc'ıa-Vergara et al. 2017) we choose the estimator  \nχ QG ( R min , R max ) = ⟨ QG ⟩ ⟨ QR ⟩ -1 , (4)  \nwhere ⟨ QG ⟩ denotes the number of detected galaxies in radial bins around the quasar, while ⟨ QR ⟩ describes the number of expected galaxies within the same volume in a random 'blank' field. Thus, in order to infer ⟨ QR ⟩ , we calculate the number of pairs between the quasar and galaxies from our random catalog within the same radial bins and | ∆ v | ≤ v max , and normalize them to the number of galaxies we would have expected to find within the same observed cylindrical volume in a 'blank' field. As before, we estimate the expected quasar-galaxy pair counts by means of the number density of galaxies inferred from the [O III ]-luminosity function multiplied  \nby the observed volume in each shell and the completeness of our survey volume. Uncertainties on the crosscorrelation measurements are estimated using Poisson noise on the detected galaxy counts 1 . Our measurements for χ QG are listed in Tab. 2 and shown in Fig. 3.", '4.2. Galaxy-Galaxy Auto-Correlation Function': "Analogous to Eqn. 3, we calculate the volumeaveraged projected auto-correlation function of galaxies χ GG ( R min , R max ) by integrating the real space galaxygalaxy auto-correlation function ξ GG ( R,Z ) in shells of volume V . We consider all galaxies within the redshift interval 5 . 95 ≤ z ≤ 6 . 55, which is approximately centered around the mean redshift of the quasars, ⟨ z ⟩ ≈ 6 . 25. Again, galaxies within the quasars' environments are excluded from the auto-correlation measurement. We then calculate pairs of galaxies separated by | ∆ v | ≤ v max .  \nWe choose the Landy-Szalay estimator (Landy & Szalay 1993), i.e.  \nχ GG ( R min , R max ) = ⟨ GG ⟩ -2 ⟨ GR ⟩ + ⟨ RR ⟩ ⟨ RR ⟩ -1 , (5)  \nwhere ⟨ GG ⟩ , ⟨ GR ⟩ and ⟨ RR ⟩ denote the number of galaxy pairs, pairs of detected galaxies and galaxies in the random catalog, as well as pairs of galaxies within the random catalog, respectively. The terms ⟨ GR ⟩ and ⟨ RR ⟩ are normalized such that the number of random galaxies would match the number of detected galaxies in  \n1 Note that in reality the uncertainties on the correlation measurements arise from a combination of Poisson counting errors as well as cosmic variance. The latter errors are, however, challenging to estimate and require extensive cosmological simulations, and thus we approximate the uncertainties by the Poisson errors alone.  \nFigure 3. The volume-averaged projected quasar-galaxy cross-correlation function is shown as black data points, while the blue data points show the galaxy-galaxy autocorrelation function at an average redshift of ⟨ z ⟩ = 6 . 25. The dashed lines and shaded areas show the best power-law fits with 2 σ uncertainties when jointly fitting both sets of measurements while keeping the slopes γ GG = 1 . 8 and γ QG = 2 . 0 fixed (see § 4.4).  \n<!-- image -->  \neach quasar field within the considered redshift interval. Uncertainties on the measurements are approximated by the Poisson noise on the number of detected galaxies. The measurements of the volume-averaged projected galaxy-galaxy auto-correlation function are shown in Fig. 3 and listed in Tab. 2.", '4.3. Fitting the Cross- and Auto-Correlation Functions': "The volume-averaged projected correlation functions χ QG and χ GG represent integrals of the real-space correlation functions ξ QG and ξ GG as shown in Eqn. 3. We parameterize the real-space cross- and auto-correlation functions by a power-law, i.e.  \nξ ( r ) = ( r/r 0 ) -γ (6)  \nwhere r = √ R 2 + Z 2 . Thus, the real-space correlation functions are governed each by two parameters, i.e. their slopes γ , denoted as γ QG and γ GG for the cross- and auto-correlation, respectively, and their scale lengths r 0 , denoted as r QG 0 and r QQ 0 .  \nWe use a Markov Chain Monte Carlo algorithm (Foreman-Mackey et al. 2013) with a Gaussian likelihood function and uniform priors for the correlation  \nscale length of r 0 ∈ [1 , 30] cMpc h -1 and the slope γ ∈ [1 , 3] to fit our measurements. The median and 16th and 84th percentiles of the posterior distribution define our the best estimate.  \nWith both the slope and scale length as free model parameters, we obtain r QG 0 = 12 . 4 +3 . 8 -2 . 3 cMpc h -1 and γ QG = 1 . 6 +0 . 2 -0 . 3 as the best fits for the quasar-galaxy cross-correlation function. When fixing the slope to γ QG = 2 . 0 to enable a comparison to other studies (see § 7), we obtain r QG 0 = 9 . 1 +0 . 5 -0 . 6 cMpc h -1 .  \nSimilarly, for the galaxy-galaxy auto-correlation, we obtain r GG 0 = 3 . 9 ± 0 . 4 cMpc h -1 and γ GG = 1 . 9 ± 0 . 1 when fitting for both scale length and slope. When fixing the slope to γ GG = 1 . 8 as in previous studies (e.g. Garc'ıa-Vergara et al. 2017), we obtain r GG 0 = 4 . 1 ± 0 . 3 cMpc h -1 .", '4.4. Jointly Fitting both Correlation Functions': "We can now attempt to study the clustering properties of the quasars themselves and infer the real-space quasar-quasar auto-correlation function ξ QQ by means of ξ QG and ξ GG . To this end, we follow a deterministic bias model (see e.g. Croom et al. 2001; Garc'ıaVergara et al. 2017) and denote the density contrast of galaxies and quasars as δ Q and δ G . Thus, the crosscorrelation between quasars and galaxies is defined as ξ QG = ⟨ δ Q δ G ⟩ . Assuming that galaxies and quasars trace the same underlying dark matter distribution, the galaxy and quasar density contrast can be expressed as δ G = b G ( δ DM ) δ DM and δ Q = b Q ( δ DM ) δ DM , respectively, where b G ( δ DM ) and b Q ( δ DM ) denote the galaxy and quasar bias, which are (possibly non-linear) functions of the dark matter density contrast δ DM . If we describe δ Q and δ G as two stochastic processes, their cross-correlation coefficient ρ can be expressed as  \nρ = ⟨ δ Q δ G ⟩ √ ⟨ δ Q δ Q ⟩⟨ δ G δ G ⟩ . (7)  \nHowever, given the deterministic relations for δ Q and δ G as defined above, it follows that ρ = 1. Thus, the quasar-quasar auto-correlation ξ QQ can be expressed as a combination of the other two correlation functions, ξ QG and ξ GG , i.e.  \nξ QQ = ξ 2 QG /ξ GG . (8)  \nAnalogously to Eqn. 6, we parameterize all autoand cross-correlation functions with a power-law, but now perform a joint fit to both the auto- and crosscorrelation measurements. To this end, we fit ξ GG as a function of γ GG and r GG 0 , as well as ξ QG parameterized by γ QQ , γ GG , r QQ 0 and r GG 0 using Eqn. 8.  \nTo enable a comparison to previous work (see § 7.1), we keep the slope of the two auto-correlation functions fixed, i.e. γ GG = 1 . 8 (Garc'ıa-Vergara et al. 2017), and γ QQ = 2 . 0 (Shen et al. 2007; Eftekharzadeh et al. 2015), and only fit the two remaining free parameters, i.e. the scale lengths r GG 0 and r QQ 0 . Our best estimates determined by the median and 16th and 84th percentile of the posterior probability distributions for the scale lengths of the galaxy auto-correlation function at ⟨ z ⟩ ≈ 6 . 25 is r GG 0 = 4 . 1 ± 0 . 3 cMpc h -1 , while we obtain r QQ 0 = 22 . 0 +3 . 0 -2 . 9 cMpc h -1 for the scale length of the quasars' auto-correlation function. These estimates for the power-laws and our measurements are shown in Fig. 3.  \nWe also attempt to jointly fit for all four model parameters, i.e. two slopes and two correlation scale lengths. We caution that the parameters are strongly degenerate and approach the prior boundaries, but the marginalized posteriors suggest r GG 0 = 4 . 0 +0 . 4 -0 . 3 cMpc h -1 , r QQ 0 = 23 . 7 +4 . 3 -5 . 6 cMpc h -1 , γ GG = 1 . 9 ± 0 . 1 and γ QQ = 1 . 9 ± 0 . 2.", "5. MASS ESTIMATES OF THE QUASARS' HOST DARK MATTER HALOS": "In this section we will estimate the mass of the dark matter halos that harbor luminous quasars at z ≳ 6 based on the quasars' clustering properties. Our cosmological model ΛCDM provides a direct link between the clustering strength of a population of objects and the characteristic mass of their host dark matter halos, since more massive halos have a higher clustering bias (e.g. Kaiser 1984; Mo & White 1996; Tinker et al. 2010).  \nWe will present three different approaches to constrain the host dark matter halo masses, all based on darkmatter-only simulations. We use N-body cosmological simulations to leverage the accurate representation of the halo clustering properties that they provide. Accurately modeling the distribution of halos in the universe is particularly important at the small scales probed by our observations (i.e. at distances of r ≲ 10 cMpc h -1 ), as these scales are highly non-linear in overdense environments. Alternative approaches based on non-linear analytical extensions of linear theory to small scales (such as the 'halo model' framework, e.g., Cooray & Sheth 2002) would likely return results that are not accurate enough for our purposes.  \nWe first use a new N-body simulation Flamingo10k , which is part of the Flamingo project (Schaye et al. 2023, see § 5.1), and assume simple 'step-function' models for the halo occupation distributions (HODs; e.g. White et al. 2008) of quasars and galaxies. We then also compare to two empirical models for galaxies and SMBHs based on the Uchuu N-body simulation  \n(Ishiyama et al. 2021), namely the UniverseMachine (Behroozi et al. 2019, see § 5.2), as well as the Trinity model (Zhang et al. 2023a, see § 5.3).", '5.1. Host dark matter halo masses based on the Flamingo-10k cosmological simulation': "As a first approach, we use the N-body cosmological simulation Flamingo-10k to directly infer the host halo masses of quasars and [O III ]-emitters by matching the observed quasar-galaxy cross-correlation and galaxy-galaxy auto-correlation measurements with the clustering properties of dark matter halos. Flamingo10k is a new dark-matter-only simulation part of the Flamingo suite (Schaye et al. 2023), which makes use of the open-source code Swift (Schaller et al. 2023), a highly-parallel gravity and smoothed particle hydrodynamics solver. The simulation was run with 10080 3 cold dark matter particles and 5600 3 neutrino particles, and has the same volume as the Flamingo flagship hydrodynamical run (2 . 8 cGpc on a side) but eight times more dark matter and neutrino particles, i.e. with a dark matter particle mass of 8 . 4 × 10 8 M ⊙ . The run adopts the '3x2pt + all' cosmology from Abbott et al. (2022) (Ω m = 0 . 306, Ω b = 0 . 0486, σ 8 = 0 . 807, H 0 = 68 . 1 kms -1 Mpc -1 , n s = 0 . 967), with a summed neutrino mass of 0 . 06 eV. Based on this simulation, we can model the clustering properties of any subset of z ≈ 6 dark matter halos in the mass range M = 10 10 . 5 -10 13 M ⊙ . Our simulation covers a total of ≈ 10 8 dark matter halos over this mass range, 60 , 857 (947, 4) of which above a halo mass of 10 12 M ⊙ (10 12 . 5 M ⊙ , 10 13 M ⊙ ). Thus, the number of lower-mass halos hosting [O III ]-emitters is very large, but we need to ensure that our simulations include sufficient high-mass halos to limit the effects of cosmic variance. Assuming that ∼ 10 -20 of such halos are sufficient, the crosscorrelation function estimated from the simulations will be independent of sampling effects due to cosmic variance for halo masses up to ∼ 10 12 . 8 -12 . 9 M ⊙ . For details on the modeling procedure and the simulation, we refer the reader to (Pizzati et al. 2024a,b).  \nAssuming a prescription for how quasars and [O III ]-emitters populate dark matter halos, we model the quasar-galaxy cross-correlation, χ QG , as well as the galaxy-galaxy auto-correlation, χ GG . Here, we assume a 'step-function' HOD model for both quasars and [O III ]-emitters, where halos are populated by quasars ([O III ]-emitters) only above a minimum threshold mass M halo , min ( M [OIII] halo , min , ). The two minimum host halo masses, M halo , min and M [OIII] halo , min , , are free parameters in our model, which we infer by jointly fitting the observed correlation functions, χ QG and χ GG . We assume that  \nthese two correlation functions are independent, and write down a 2-d Gaussian likelihood for the parameters M halo , min and M [OIII] halo , min , . By computing the median and 16th-84th percentiles of the marginalized likelihood distributions, we infer the minimum host dark matter halo mass of quasars to be log 10 ( M halo , min /M ⊙ ) = 12 . 43 +0 . 13 -0 . 15 , while we find that [O III ]-emitters are hosted in halos of mass log 10 ( M [OIII] halo , min /M ⊙ ) = 10 . 56 +0 . 05 -0 . 03 . In a companion paper (Pizzati et al. 2024b), we model the full mass distribution of the dark matter halos hosting quasars and galaxies by jointly fitting the clustering properties as well as the luminosity functions of quasars and [O III ]-emitters, and find results that are in excellent agreement with the model presented here.  \nIn the left panel of Fig. 4, we compare the predictions from the Flamingo-10k simulation with the observed quasar-galaxy cross-correlation measurements, for different values of the minimum halo mass for the quasar hosts, M halo , min . The cross-correlation terms are obtained by assuming a minimum host halo mass for [O III ]-emitters equal to the median value determined above, i.e., log 10 ( M [OIII] halo , min /M ⊙ ) = 10 . 56. The best fit between the simulated and observed cross- and autocorrelation functions is obtained with a minimum quasar host halo mass of log 10 ( M halo , min /M ⊙ ) ≈ 12 . 43.  \nNote that the slope of the cross-correlation functions at very small scales, i.e. r ≲ 200 ckpc h -1 , appears to differ between the observations and the simulations. This could either be due to small number statistics given our current data set, as we only observe a total of 5 galaxies for the inner most two data points, or alternatively, uncertainties in the modeling procedure could cause the discrepancy. For instance, the identification of subhalos becomes increasingly challenging at such small separations to a massive structure, and could thus alternate the shape of the correlation functions at small scales. However, our conclusions are not affected by the two inner most data points, as the uncertainties on these measurements are high and thus do not influence the best fit significantly.  \nAdditional differences at the smallest scales arise since we do not apply the same merging procedure to the subhalos in the simulations that we apply to the detected [O III ]-emitters as described in § 2.3. Merging the subhalos would lead to a slightly shallower slope in the simulated cross-correlation function at the smallest scales of ≲ 3 '' , but since the uncertainties on our measurements on the smallest scales are large due to the low number of detected galaxies, this has a negligible effect on our estimates. When applying the same merging procedure to subhalos in the Uchuu simulations (see § 5.2) as a check, we found that  \n∼ 2% of the subhalos were affected compared to the ∼ 13% of [O III ]-emitters affected by merging, which indicates that the [O III ]-emitters show mostly clumpy substructure within an individual subhalo.  \nIn the following, we will briefly discuss two other methods to determine the host dark matter halo masses from the quasars' clustering properties to ensure consistency between different simulations and models, but moving forward, we will consider the measurement obtained here as our best estimate for the minimum host dark matter halo mass of z ≳ 6 quasars.", '5.2. Host dark matter halo masses based on the UniverseMachine': "A second method to determine the dark matter halo masses of the observed quasars is based on finding analogous systems to the observed quasar environments using the UniverseMachine (Behroozi et al. 2019) applied to the Uchuu N-body simulations. The details of this approach will be presented in detail in a forthcoming companion paper (Mackenzie et al. in prep.). In short, by comparing the number of simulated [O III ]-emitters around halos of different masses to the number of companion galaxies observed in the quasar fields, a distribution of possible host dark matter halo masses can be determined for each observed quasar.  \nWe first construct a catalog of simulated environments around each halo, and select all halos above 10 10 h -1 M ⊙ (at z = 6 . 35). In total this simulation contains ∼ 10 8 . 5 such halos, approximately 1 . 3 × 10 6 (24 , 525, 1 , 860, 85, 8) of which are above a mass threshold of 10 11 . 5 M ⊙ (10 12 . 1 M ⊙ , 10 12 . 4 M ⊙ , 10 12 . 7 M ⊙ , 10 12 . 9 M ⊙ ). Each halo has its own realisation where the halo position becomes the center of a mock observation (i.e. at the quasar position). The surrounding galaxies are then forward modeled through our completeness function and survey geometry. By calculating the joint distribution of the number of detected neighbouring galaxies vs. halo mass, and applying Bayes equation with a flat prior on halo mass, one can estimate the mass probability distribution for an observed number of [O III ]-emitting neighbours. The resulting most probable dark matter halo masses range between 10 10 . 85 M ⊙ ≲ M halo ≲ 10 12 . 85 M ⊙ with a median of ⟨ M halo ⟩ = 10 12 . 4 ± 0 . 5 M ⊙ , where the wide range of host halo masses is largely driven by the cosmic variance observed between the different quasar fields. The estimate for the minimum host dark matter halo masses using constraints from the UniverseMachine are well consistent with the estimate based on the Flamingo10k simulations.  \n5.3. Host dark matter halo masses based on Trinity  \n<!-- image -->  \nFigure 4. Left: Comparison of the quasar-galaxy cross-correlation measurement to models from the Flamingo-10k simulation for different minimum quasar host dark matter halo masses, M halo , min . By comparing the correlation functions of dark matter halos from the simulations to the measured quasar-galaxy cross-correlation and galaxy-galaxy auto-correlation, we obtain a best estimate for the quasars' minimum host halo mass of log 10 ( M halo , min /M ⊙ ) = 12 . 43 +0 . 13 -0 . 15 . Right: Comparison between the measurements of the quasar-galaxy cross-correlation with predictions (median: dashed yellow curve; 16th-84th percentile: shaded region) for the observed quasar fields from the Trinity model, which imply an average host dark matter halo mass of log 10 ( M halo /M ⊙ ) = 12 . 14 +0 . 24 -0 . 26 .  \n<!-- image -->  \nFinally, we compare the observed χ QG with the prediction from the Trinity model (Zhang et al. 2023a,b) of the halo-galaxy-SMBH connection, to infer the typical host dark matter halo mass of the quasars. Briefly, Trinity reconstructs empirical galaxy and black hole mass growth rate distributions as functions of halo mass from redshifts z = 0 -10, and forward models observables such as galaxy mass functions and quasar luminosity functions. SMBH properties such as radiative efficiencies, Eddington ratio distributions, and duty cycles are parameterized and constrained by observational data. For details on the Trinity implementation, we refer the reader to Zhang et al. (2023a).  \nWe utilize the Trinity mock catalogs of galaxies and SMBHs, with position, velocity, halo mass, galaxy mass, star formation rate, rest-frame UV luminosity, as well as black hole mass and bolometric luminosity (if hosting a SMBH) for each galaxy. These mock catalogs are generated also from the Uchuu N-body simulation (box size: 2 cGpc h -1 , mass resolution: 3 . 27 × 10 8 M ⊙ h -1 , Ishiyama et al. 2021), and will be made publicly available in a forthcoming Trinity publication (Zhang et al. in prep.). In order to calculate χ QG and host halo mass distributions, we use all quasars in our z = 6 . 35 mock catalog, and assign each mock quasar a weight to quan-  \nits degree of similarity to the here analyzed quasars. Specifically, the weight w for a simulated quasar with black hole mass and bolometric luminosity ( M BH , L bol ) is a sum of the log-normal probabilities given the observed mass and luminosities, { ( M BH ,i , L bol ,i ) } N =4 i =1 , of the EIGER quasars, i.e.  \nw = N =4 ∑ i =1 P [( M BH , L bol ) | ( M BH ,i , L bol ,i )] = N =4 ∑ i =1 1 2 πσ L bol ,i σ M BH ,i × exp [ -(log M BH -log M BH ,i ) 2 2 σ 2 M BH ,i ] × exp [ -(log L bol -log L bol ,i ) 2 2 σ 2 L bol ,i ] , (9)  \nwhere σ M BH ,i and σ L bol ,i are the measurement uncertainties in SMBH mass and luminosity of the i th quasar, respectively. For σ M BH ,i , we also add a scatter of 0.5 dex in quadrature to account for the random scatter in virial estimates of M BH (McLure & Jarvis 2002; Vestergaard & Peterson 2006).  \nWhen calculating the predicted χ QG we convert the galaxy UV luminosity, L UV , into L [OIII] , 5008 by applying the scaling relation from Matthee et al. (2023a), and include all mock galaxies with L [OIII] , 5008 ≥ 10 42 erg s -1 . To quantify χ QG uncertainties from Trinity , we simulate the observation 3000 times, i.e. we choose four mock quasars for every observation, where each quasar is randomly selected within the 1 . 5 σ range of each observed quasar in the (log M BH , log L bol ) space. We then calculate the 16 th -84 th percentile range of the measured χ QG 's as the uncertainties from Trinity .  \nAs shown in right panel of Fig. 4, the resulting χ QG from Trinity are consistent with the observed values within the measurement uncertainties. This consistency allows us to infer the host dark matter halo masses of the quasars by examining the host halo mass distribution of the mock quasars, weighted by the weights defined in Eqn. 9. From this exercise, we obtain a median (with 16 th -84 th percentile) host halo mass of log 10 ( M halo /M ⊙ ) = 12 . 14 +0 . 24 -0 . 26 .", "6. ESTIMATING THE QUASARS' DUTY CYCLE AND LIFETIME": "Our cosmological model ΛCDM allows us to determine the number density of the quasars' host dark matter halos n hosts from their characteristic mass determined in the previous section (e.g. Mo & White 1996; Tinker et al. 2010), which in turn enables constraints on the quasars' duty cycle and lifetime. As quasars temporarily subsample their hosts, their number density n Q can be expressed as  \nn Q ≃ t Q t H ( z ) n hosts , (10)  \nwhere t Q is the quasars' UV-luminous lifetime, and t H ( z ) denotes the Hubble time at redshift z , and hence f duty = t Q /t H ( z ) is the quasars' duty cycle at the given redshift (Martini & Weinberg 2001; Haiman & Hui 2001). We will now determine these timescales of quasar activity, which are essential for understanding the concomitant growth phases of early SMBHs.  \nUsing the Flamingo-10k simulation we find that the number density of halos with a mass larger than log 10 ( M halo , min /M ⊙ ) = 12 . 43 +0 . 13 -0 . 15 is n host = 5 . 3 +16 . 7 -3 . 9 × 10 -8 cMpc -3 . For estimating the number density of luminous quasars we use the quasar luminosity function (QLF) at ⟨ z ⟩ ≈ 6 . 25 by Schindler et al. (2023) and integrate it over the magnitude range -26 . 5 > M 1450 > -29 . 5 bracketing the observed magnitudes of the analyzed quasars (see Tab. 1), which results in a quasar number density of n Q ≈ 0 . 23 cGpc -3 . Using Eqn. 10 the two number densities determine a duty cycle of  \nf duty = 0 . 4 +1 . 2 -0 . 3 %, and hence a UV-luminous quasar lifetime of approximately t Q ∼ 10 6 . 5 yr at ⟨ z ⟩ = 6 . 25.  \nHowever, we caution that these timescale estimates depend sensitively on the number density of the host dark matter halos given the characteristic minimum halo mass derived in § 5. Thus even relatively small differences in the host halo mass estimates of approximately ± 0 . 25 dex can lead to variations in the derived timescales of ± 1 dex, i.e. a host dark matter halo mass of log 10 ( M halo , min /M ⊙ ) ≈ 12 . 65 (log 10 ( M halo , min /M ⊙ ) ≈ 12 . 18) with a number density of a factor of 10 lower (higher) than our best estimate would correspond to a UV-luminous quasar lifetime of t Q ∼ 10 7 . 5 yr ( t Q ∼ 10 5 . 5 yr). And yet, even with this uncertainty, these inferred UV-luminous lifetime estimates represent a major departure from the expected long timescales of t Q ∼ 10 9 yr and f duty ∼ 1, required to grow the quasars' SMBHs. In order to match the abundance of host dark matter halos to the abundance of luminous quasars to ensure a duty cycle of unity, the required host dark matter halo mass would have to be ∼ 10 13 M ⊙ .", '7.1. Comparison to previous studies': "There is only one study to date measuring the autocorrelation function of z ≳ 6 quasars, which is the recent result by the SHELLQs Collaboration (Arita et al. 2023). The authors analyze the auto-correlation function of approximately one hundred faint ( M 1450 ≳ -25) quasars at z ∼ 6 distributed over 891 deg 2 . They probe scales between r ∼ 10 -1000 cMpc h -1 , i.e. much larger scales than the scales of r < 10 cMpc h -1 we probe with our NIRCam WFSS data, and measure an autocorrelation scale length of r QQ 0 = 23 . 7 ± 11 cMpc h -1 , which is in excellent agreement to our inferred estimate of r QQ 0 = 22 . 0 +3 . 0 -2 . 9 cMpc h -1 . They infer a host dark matter halo mass using the linear halo model software HaloMod (Murray et al. 2013, 2021), and obtain a significantly larger host dark matter halo mass estimate of log 10 ( M halo , min /M ⊙ ) = 12 . 9 +0 . 4 -0 . 7 .  \nAn independent approach to infer the host dark matter halos of a sample of luminous ( M 1450 < -26 . 5) quasars at z ∼ 6 with similar properties to the ones analyzed in this study is presented in Chen et al. (2022). The authors use the flux transmission profile observed in the rest-frame UV quasar spectra to infer the density field around the quasars and obtain an estimate of the typical host dark matter halo masses of log 10 ( M halo /M ⊙ ) = 12 . 5 +0 . 4 -0 . 7 , which is consistent with our results.  \nIn Fig. 5 we compare our results for the autocorrelation length, host dark matter halo mass and duty cycle to these previous analyses at z ∼ 6. We also add  \nFigure 5. Redshift evolution of the quasars' auto-correlation scale length ( top ), the minimum host dark matter halo mass ( middle ), and the quasars' duty cycle ( bottom ). Our inferred scale length is well consistent with the recent auto-correlation measurement by Arita et al. (2023), but we infer a smaller host dark matter halo mass using dark matter only simulations (see § 5). Our results indicate a very small duty cycle of quasars at z ≳ 6, which agrees with independent duty cycle measurements inferred from Ly α damping wings (diamond shaped data points; Davies et al. 2019; ˇ Durovˇc'ıkov'a et al. 2024). The host dark matter halo mass estimate by Chen et al. (2022) is also inferred from the transmitted flux along quasar sightlines, while all other measurements are derived from quasar clustering studies (circular data points; Shen et al. 2007; White et al. 2012; Eftekharzadeh et al. 2015; Laurent et al. 2017; Arita et al. 2023; Pizzati et al. 2024a). The grey dotted lines in the bottom panel indicate different quasar lifetimes.  \n<!-- image -->  \nmeasurements at lower redshifts to study the redshift evolution of the parameters. We confirm previous findings that the dark matter halo mass of quasars appears to be approximately constant throughout cosmic time, suggesting that there is a characteristic halo mass of a few times M halo ∼ 10 12 M ⊙ where quasars can be activated (e.g. Croom et al. 2005; Arita et al. 2023).  \nInterestingly, we know that galaxies typically quench at stellar masses M ⋆ ∼ 10 11 M ⊙ across a wide range of redshifts out to z ∼ 5 (e.g. Peng et al. 2010; Behroozi  \net al. 2013; Caplar et al. 2015), corresponding to a halo mass M halo ∼ 10 12 . 0 -12 . 5 M ⊙ , assuming an approximately constant stellar mass to halo mass ratio. Since there is considerable evidence that the quenching of galaxies is associated with, or possibly caused by, quasar activity (e.g. Cano-D'ıaz et al. 2012; Caplar et al. 2018), the identification of this same halo mass as the mass of the EIGER quasars' host halos at z ∼ 6 may not be unexpected.", '7.2. Overdensities around high-redshift quasars?': 'The numerous previous inconclusive results aiming to determine the existence of an overdense environment around high-redshift quasars have cast doubt on whether these quasars indeed trace the highest-density peaks in our universe. Our clustering measurement suggests that quasars reside in dark matter halos with a minimum mass of log 10 ( M halo , min /M ⊙ ) = 12 . 43 +0 . 13 -0 . 15 , which implies that they do on average trace large overdensities but not necessarily the rarest and highest density peaks.  \nOur sample of four quasar fields has revealed a wide range of overdensities providing evidence for large cosmic variance. The observed density contrasts in the quasar fields range from δ ( r < 2 cMpc) = 65 ± 15 to δ ( r < 2 cMpc) = 3 ± 4. While J0148+0600 with 47 [O III ]-emitting galaxies in its surroundings likely resides in one of the highest-density peaks, and traces one of the largest high-redshift proto-cluster reported to date, J1030+0542 on the other hand appears to reside in an average density region with no enhancement compared to a random field (which is in contrast to the protocluster reported in this field by Mignoli et al. 2020).  \nIt is interesting to note that the [C II ] emission arising from the molecular gas of the host galaxy of J1030+0542 was not detected with ALMA, indicating that the FIR luminosity of this quasar is significantly fainter than for most other z ∼ 6 quasars with similar rest-frame UV magnitudes (Decarli et al. 2018). This result, combined with the recent discoveries of large overdensities around sub-mm bright galaxies at z > 5 (Herard-Demanche et al. 2023; Sun et al. 2024), provides circumstantial evidence that FIR luminosities might correlate more strongly with dark matter halo mass than the UV luminosities of high-redshift objects do.', '7.3. Implications for the early growth of SMBHs': "As the emission of quasar light is concomitant with black hole accretion, the timescales of quasar activity have implications for the growth of SMBHs. Previous constraints on the quasars' duty cycle based on clustering studies at redshifts of 1 ≲ z ≲ 4 suggest broadly  \nan increase towards higher redshifts 2 as shown in the bottom panel of Fig. 5, implying relatively long UVluminous quasar lifetimes of t Q ∼ 10 7 -10 8 yr (e.g. Martini 2004).  \nOur standard black hole growth model requires such long quasar lifetimes in order to explain the growth of the observed billion solar mass black holes of highredshift quasars, and consequently high quasar duty cycles of f duty ∼ 1 have been postulated for quasars in the early universe (e.g. Martini 2004; Volonteri 2012). However, the clustering analysis of z ≳ 6 quasars in this study surprisingly provides evidence for the contrary: we find a very small duty cycle, f duty ≪ 1, indicating that galaxies do not shine as UV-luminous quasars throughout the whole Hubble time at these redshifts, thus posing additional challenges to our understanding of early SMBH growth.  \nInterestingly, such low values are consistent with completely independent constraints on the UV-luminous duty cycle of quasars at similar redshifts derived from the extents of their proximity zones (e.g. Eilers et al. 2017, 2018, 2020; Andika et al. 2020; Morey et al. 2021; Satyavolu et al. 2023b) and damping wing features (e.g. Davies et al. 2019; ˇ Durovˇc'ıkov'a et al. 2024) that are observed in the quasars' rest-frame UV spectra. These studies estimate the number of ionizing photons emitted by a quasar that would results in the observed flux transmission profile along the line-of-sight, and thus enable constraints on their emission timescales for individual objects 3 . These observations have previously suggested UV-luminous quasar lifetimes of only t Q ∼ 10 6 years, which imply similarly low duty cycles of f duty ≪ 1 for high-redshift quasars (see also Khrykin et al. 2019, for a similar anlysis of He II proximity zones at z ∼ 4).  \nThus, two independent arguments - one based on the quasars' clustering properties and the other based on the observed flux transmission profiles in their spectra provide evidence for very short UV-luminous quasar ac-  \n2 We note that the increase in the quasars' duty cycle with redshift is not found by Eftekharzadeh et al. (2015), who finds significantly lower duty cycles than Shen et al. (2007) at z ∼ 3 -4. However, their quasar sample using BOSS includes much fainter objects, which makes the selection of quasars at these redshifts as well as the completeness evaluation of their sample more challenging.  \n3 Note that for quasars that are embedded in an already highly ionized IGM, with a neutral gas fraction of x HI ∼ 10 -4 at z ∼ 6, the proximity zone size reflects the lifetime of the last quasar episode and hence the duty cycle presents an upper limit on these episodic lifetime estimates. For quasars at z ≳ 7 that are embedded in a still fairly neutral IGM, x HI ≳ 0 . 1, the observed damping wing signature in the proximity zone allows for an estimate of the integrated quasar lifetime, i.e. the duty cycle, due to the long recombination time (see e.g. Davies et al. 2019).  \ntivity timescales in the early universe. This suggests that accretion onto the SMBHs occurs in highly radiatively inefficient phases with a radiative efficiency ϵ of less than the fiducial ∼ 10% suggested by thin accretion disk models (e.g. Davies et al. 2019), or, alternatively, a significant fraction of the black hole growth happens in UVobscured, dust-enshrouded environments (e.g. Satyavolu et al. 2023a).  \nThe latter is a particularly intriguing scenario in light of the recent discovery of an abundance of faint and highly dust-reddened active galactic nuclei (AGN) present at early cosmic times (Harikane et al. 2023; Matthee et al. 2023b; Kocevski et al. 2023; Maiolino et al. 2023; Greene et al. 2023), which might provide evidence for a highly obscured population of growing SMBHs. However, studies in the mid-infrared such as the recently published Systematic Mid-infrared Instrument Legacy Extragalactic Survey (SMILES; Lyu et al. 2023) are necessary to identify completely dust-obscured growing SMBHs and securely determine the fraction of obscured quasars in the early universe. If black hole growth in completely dust-enshrouded, UV-obscured environments is indeed responsible for the discrepancy between the theoretically expected long black hole growth timescales and the observed short UV-luminous growth phases, we would expect ∼ 100 -1000 dust-obscured quasars for every unobscured one, indicating the presence of a very large population of dust-obscured quasars in the early universe that we yet have to detect. While this would imply a very high obscuration fraction of > 99%, cosmological simulations (Trebitsch et al. 2019; Ni et al. 2020; Vito et al. 2022; Bennett et al. 2024) as well as a recent study using ALMA observations (Gilli et al. 2022) have shown that obscuration fraction of at least ∼ 80 -90% are expected for high-redshift quasars due to the presence of high-column density gas within the innermost regions of their host galaxies, resulting in low observable UV-luminous duty cycles.", '8. SUMMARY': "In this paper we present the first clustering measurement of luminous quasars at z ≳ 6 and [O III ]-emitting galaxies in their environment in four quasar fields observed with JWST/NIRCam in imaging and WFSS mode as part of the EIGER project. We identify 522 [O III ]-emitting galaxy systems above a luminosity of L [OIII] , 5008 = 10 42 erg s -1 in all five observed quasar fields, 92 of which reside in close vicinity, i.e. | ∆ v | ≤ 1000 kms -1 , to the luminous quasars.  \nWe find a large diversity in the observed density contrast between the different quasar fields, ranging from density enhancements of δ ( r < 2 cMpc) = 65 ± 15 - one  \nof the most spectacular proto-clusters within the Epoch of Reionization discovered to date encompassing 47 [O III ]-emitting galaxies around the quasar J0148+0600 - to δ ( r < 2 cMpc) = 3 ± 4, indicating no significant galaxy enhancement around the quasar J1030+0564. The other two fields around the super-luminous quasar J0100+2802 and the quasar J1148+5251, show density contrasts of δ ( r < 2 cMpc) = 29 ± 10 and δ ( r < 2 cMpc) = 16 ± 8, respectively.  \nWe present measurements of the volume-averaged projected galaxy-galaxy auto-correlation function and the quasar-galaxy cross-correlation function at an average redshift of ⟨ z ⟩ = 6 . 25. We parameterize the real-space correlation functions with a power-law and obtain estimates for the auto- and cross-correlation lengths of r GG 0 = 4 . 1 ± 0 . 3 cMpc h -1 and r QG 0 = 9 . 1 +0 . 5 -0 . 6 cMpc h -1 at ⟨ z ⟩ = 6 . 25, respectively (keeping the slopes γ GG = 1 . 8 and γ QG = 2 . 0 fixed). When jointly fitting the two correlation functions to infer the scale length of the quasar-quasar auto-correlation, we measure r QQ 0 = 22 . 0 +3 . 0 -2 . 9 cMpc h -1 (with fixed slope γ QQ = 2 . 0).  \nBy comparing our measurements to dark matter only simulations and different empirical models, we find consistent estimates for the host dark matter halo mass: We determine a minimum host halo mass of log 10 ( M halo , min /M ⊙ ) = 12 . 43 +0 . 13 -0 . 15 as the best fit between our measurements and the cross-correlation function of dark matter halos in the Flamingo-10k simulation. When searching for analog systems to the quasar fields in UniverseMachine we find that such quasars are hosted in dark matter halos with average masses of log 10 ( M halo /M ⊙ ) = 12 . 4 ± 0 . 5. Lastly, when comparing the cross-correlation function to predictions from the Trinity model, we find host dark matter halo masses with a median of log 10 ( M halo /M ⊙ ) = 12 . 14 +0 . 24 -0 . 26 . All three estimates agree with each other within the uncertainties despite their very different approaches.  \nThese results for the mass of the quasars' host dark matter halos suggest that quasars do not necessarily reside in the rarest, highest-density peaks, and could explain why some previous studies did not find large overdensities of galaxies around high-redshift quasars. It also indicates a smaller than expected duty cycle of quasars in the early universe, i.e. f duty ≪ 1. While the short UV-luminous duty cycles are consistent with independent results analyzing the flux transmission observed in rest-frame UV quasar spectra, the findings pose significant challenges for the growth of the early SMBHs which power the quasars' emission. Such short UV-luminous quasar lifetimes imply that (a) either the enormous ≳ 10 9 M ⊙ black holes must have build up their mass extremely rapidly, e.g. via radiatively inef-  \nt, 'super-Eddington' accretion phases; or (b) a significant fraction ( ≳ 99%) of the black hole growth happens in dust-enshrouded and UV-obscured environments, indicating the existence of a vast population of obscured quasars in the early universe.  \nThe authors would like to thank the anonymous referee for thoughtful comments, which significantly improved our manuscript, and Jan-Torge Schindler, Jiamu Huang and Feige Wang for helpful discussions.  \nJFH and EP acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 885301). JM acknowledges support from the European Union (ERC, AGENTS, 101076224).  \nThis work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The JWST data presented in this article were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. The specific observations analyzed are associated with program #1243, and can be accessed via DOI: 10.17909/m5mp5v90.  \nThis work used the DiRAC Memory Intensive service (Cosma8) at the University of Durham, which is part of the STFC DiRAC HPC Facility (www. dirac.ac.uk). Access to DiRAC resources was granted through a Director's Discretionary Time allocation in 2023/24, under the auspices of the UKRI-funded DiRAC Federation Project. The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.  \nWe thank Instituto de Astrofisica de Andalucia (IAA-CSIC), Centro de Supercomputacion de Galicia (CESGA) and the Spanish academic and research network (RedIRIS) in Spain for hosting Uchuu DR1, DR2 and DR3 in the Skies & Universes site for cosmological simulations. The Uchuu simulations were carried out on Aterui II supercomputer at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan, and the K computer at the RIKEN Advanced Institute for Computational Science. The Uchuu Data Releases efforts have made use of the skunIAA RedIRIS and skun6IAA computer facilities managed by the IAA-CSIC in Spain (MICINN EUFeder grant EQC2018-004366-P).  \nSoftware: numpy(Harrisetal.2020),scipy(Jonesetal. 2001),matplotlib(Hunter2007),astropy(AstropyCollab-  \noration et al. 2013, 2018, 2022), emcee (Foreman-Mackey et al. 2013)"}