The Cosmic Ultraviolet Baryon Survey (CUBS). VIII. Group Environment of the Most Luminous Quasars at z ≈ 1

We investigate the group-scale environment of 15 luminous quasars (luminosity L 3000 > 1046 erg s−1) from the Cosmic Ultraviolet Baryon Survey (CUBS) at redshift z ≈ 1. Using the Multi Unit Spectroscopic Explorer integral field spectrograph on the Very Large Telescope, we conduct a deep galaxy redshift survey in the CUBS quasar fields to identify group members and measure the physical properties of individual galaxies and galaxy groups. We find that the CUBS quasars reside in diverse environments. The majority (11 out of 15) of the CUBS quasars reside in overdense environments with typical halo masses exceeding 1013 M ⊙, while the remaining quasars reside in moderate-size galaxy groups. No correlation is observed between overdensity and redshift, black hole (BH) mass, or luminosity. Radio-loud quasars (5 out of 15 CUBS quasars) are more likely to be in overdense environments than their radio-quiet counterparts in the sample, consistent with the mean trends from previous statistical observations and clustering analyses. Nonetheless, we also observe radio-loud quasars in moderate groups and radio-quiet quasars in overdense environments, indicating a large scatter in the connection between radio properties and environment. We find that the most UV luminous quasars might be outliers in the stellar mass-to-halo mass relations or may represent departures from the standard single-epoch BH relations.


INTRODUCTION
In hierarchical structure formation models, small initial density fluctuations grow to form galaxies and galaxy clusters in the present Universe.Residing at the centers of massive galaxies, quasars can provide substantial feedback into the surrounding environment through radio jets and direct heating of the interstellar and circumgalactic media (ISM/CGM) (e.g., Silk & Rees 1998;Heckman & Best 2014).The observed tight correlations between SMBHs and host galaxy properties support the BH-galaxy co-evolution models (e.g., Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Kormendy & Ho 2013, and references therein).Studying the environment quasars reside in provides in-sights into galaxy/BH evolution, the interplay between quasars, their surrounding gas reservoirs, and hierarchical structure formation.
Previous studies of quasars and their environments often relied on wide-field optical/infrared imaging surveys and clustering analyses.The overdensity of quasar environments is estimated by comparing the galaxy number counts within certain projected distances around known quasars and background fields.Many studies found that the galaxy number density around quasars is higher than around inactive galaxies of similar masses (e.g., Serber et al. 2006;Wylezalek et al. 2013;Karhunen et al. 2014).This could indicate that quasars reside in relatively massive halos since environmental density is tightly corre-Li et al. lated with halo mass (Haas et al. 2012).Radio-loud quasars, in particular, are found to reside in denser environments compared to radio-quiet quasars and inactive galaxies of similar masses over a wide redshift range (e.g., Ramos Almeida et al. 2013;Wylezalek et al. 2013;Hatch et al. 2014).If radio-loudness is connected to the overdensity of quasar environments or quasar properties, then this could be evidence connecting merger and galaxy interaction to the triggering of nuclear activity and jet formation (e.g., Di Matteo et al. 2005;Hopkins et al. 2006).On the other hand, it is unclear if overdensity correlates with other quasar properties, e.g., redshift, quasar luminosity, and BH mass (e.g., Serber et al. 2006;Karhunen et al. 2014).These results are broadly consistent with clustering analyses at low redshift, where there is a weak luminosity dependency in quasar clustering and radio-loud quasars typically reside in more clustered environments than radio-quiet quasars at fixed optical luminosity (e.g., Shen et al. 2009;Zhang et al. 2013;Shen et al. 2013).However, these methods can only provide a statistical view of the quasar environments, since the foreground or background galaxies in the fields cannot be removed without measuring individual galaxy redshifts.Stott et al. (2020) found that 8 out of 12 quasar fields display galaxy overdensities at 1 < z < 2 using a deep HST grism spectroscopy survey.However, Wethers et al. (2022) found that quasars and galaxies reside in similar-sized galaxy groups when controlled for stellar mass and redshift, arguing against the scenario of merger/interaction triggered quasar activity.Instead, they find quasars are more likely to be central galaxies, which potentially links quasar activities to either gas accretion or rich group environments.In addition, Stone et al. (2023) found that the galaxy group members around quasars and their inactive counterparts have similar stellar properties, suggesting that quasar feedback does not have a strong influence on the galaxy group members.With these mixed results, more comprehensive studies of individual quasars, their group environment and neighbors, and surrounding gas flow are needed to understand the interplay between quasars and their environments.
Wide-field integral field spectrographs, like the Multi Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT, Bacon et al. 2010), are capable of capturing the galaxy group environment by providing redshifts for all galaxies in the field of view (above a certain flux threshold).In this paper, we characterize the environment of 15 quasars observed by MUSE in the Cosmic Ultraviolet Baryon Survey (CUBS, Chen et al. 2020), and how they depend on the central BH properties.The 15 CUBS quasars are selected from the brightest quasars in the GALEX near-UV bandpass (NUV; 1770-2730 Å) at z ≈ 1 in the Dark Energy Survey footprint (Sevilla-Noarbe et al. 2021).The quasar fields are selected without prior knowledge of the surrounding galactic environment or radio properties.
This paper is organized as follows.We describe the CUBS survey and relevant follow-up observations in Sec- tion 2 and our data analysis in Section 3. The main results are presented in Section 4. We discuss our results in Section 5 and conclude in Section 6.Throughout this paper, we adopt a flat ΛCDM cosmology with Ω M = 0.3, Ω λ = 0.7, and H 0 = 70 km s −1 Mpc −1 .

OBSERVATIONS
The main goal of CUBS is to map the cosmic baryonic reservoir at intermediate redshift (0.8 ≲ z ≲ 1.4) using high-quality UV absorption spectroscopy and deep ground-based optical and near-infrared observations.The UV absorption spectroscopy is obtained from the Cosmic Origins Spectrograph (COS, Green et al. 2012) in a large Hubble Space Telescope (HST) Cycle 25 General Observer Program (GO-CUBS; PID = 15163; PI: Chen).The ground-based spectroscopic observations consist of three main components using the Inamori-Magellan Areal Camera and Spectrograph (IMACS; Dressler et al. 2011) on the Magellan Baade Telescope, the Low Dispersion Survey Spectrograph 3 (LDSS3) on the Magellan Clay Telescope, and the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT) that cover different depths and projected angular radius from the quasar sight lines.All quasar fields are covered by the Dark Energy Survey footprint, and additional H-band photometry is obtained by the Four Star Infrared Camera (Persson et al. 2013) on the Magellan Telescopes.The detailed survey design is described in Chen et al. (2020), here we summarize the observations and data analysis relevant to this work.

MUSE Observation
The MUSE galaxy survey is the deepest and most compact component of our ground-based optical spectroscopic observations, aiming to target galaxies as faint as 0.01 L * at z = 1 within ≲ 30 ′′ from the central quasar.The MUSE survey is carried out on the VLT UT4 in service mode under program ID 0104.A-0147 (PI: Chen).The observations cover 1 ′ × 1 ′ fields of view with plate scales of 0.2 ′′ and spectral resolution of ≈ 120 km s −1 at 7000 Å.The observations are in the wide-field mode (WFM) with ground layer adaptive optics (AO) assistance to ensure uniform image quality of < 0.8 ′′ across all fields.The limiting magnitude (3σ) of the detected sources in the MUSE datacubes are AB(r) ≈ 26 in the psuedo-r band (Qu et al. 2023), sufficient to identify faint galaxies down to M ⋆ ≈ 10 8 − 10 9 M ⊙ .The MUSE data cubes are reduced using the standard ESO MUSE pipeline (Weilbacher et al. 2020) and the custom data reduction package CUBEXTRACTOR (Cantalupo et al. 2019).Additional details of the MUSE observations are described in Chen et al. (2020).

FourStar H-Band Imaging
Deep H-band images were obtained in October 2017 using the FourStar Infrared Camera (Persson et al. 2013) on the Magellan Telescopes.The data reduction was performed using the FourCLift custom package (details described in Kelson et al. 2014) and following the procedures in Kelson et al. (2014) and Rudie et al. (2017).The exposure in each field is roughly 2000−3000 seconds, and the limiting magnitude (3σ) of the detected sources are ≈ 24.5 mag (Qu et al. 2023).The median seeing of the final images is ≈0.5 ′′ .

Supplementary Spectra
While using the MUSE WFM AO observation setup ensures a uniform imaging resolution, the use of wavefront lasers produces a gap at 5800 to 5965 Å in the spectra.We obtained supplementary spectra from the Magellan Echellette (MagE) Spectrograph for quasars with the MgII line in/near the sodium gap (J0333, J2135, J2245, J2308) to ensure good spectral coverage around the broad MgII line for BH mass estimation.The observations were carried out on September 29 and 30 2021 on the Magellan Baade Telescope.Two to three exposures of 300-600 seconds, depending on the quasar luminosity, were taken for each source to mitigate cosmic rays.The final reduced spectral range is 4000 Å to 9900 Å with median spectral resolution of ≈1 Å and spectral sensitivity (1σ) of 1-4×10 −17 erg cm −2 s −1 Å −1 .We follow the standard data reduction pipeline with CarPy (Kelson et al. 2000;Kelson 2003) to reduce the data and match the flux scale to the MUSE spectra.

MUSE Redshift Survey
We identify all continuum sources in the MUSE whitelight image with Source Extractor (Bertin & Arnouts 1996) for redshift measurements.For the quasars, we extract the optical spectra within apertures that include >95% of the quasar light and measure the redshifts based on narrow [O II] emission lines, which results in a typical redshift uncertainty of 20-30 km s −1 (Hewett & Wild 2010).For galaxies, the spectra of each source are extracted within a ≈0.6 ′′ radius using the Python MUSE data analysis package mpdaf (Bacon et al. 2016).We then follow the procedures described in Johnson et al. (2018) and Helton et al. (2021) to determine the redshift of each source.To briefly summarize, each spectrum is fitted with a linear combination of SDSS galaxy eigenspectra over a wide redshift grid, and the best-fit redshift is identified at the global minimum in the resulting χ 2 grid.We then visually inspect the best-fit spectra to assign a quality ranking for each redshift.We define the robust "double-feature" redshifts to have good fits around at least two spectral features (including both emission and absorption, visually evaluated by trained team members) and the "single-feature" redshift to have good fits around one spectral feature, and the rest of the redshifts are discarded.It is possible that the single-feature redshifts are not unique solutions, in these cases, we only include them in the final redshift measurement after all other possible solutions (e.g., a different emission line is well-fitted to the spectral feature) are ruled out.
Figure 1 shows the redshift distribution.Across all 15 quasar fields, we identified 1542 continuum sources (excluding stars) and measured 712 "double-feature" redshifts and 109 "single-feature" redshifts.The redshift success rates are approximately 40-65% across all fields.The primary reason for not measuring a redshift is usually because of low signal-to-noise ratios in the spectra or the lack of spectral features.Our MUSE redshift survey is most sensitive to redshifts z < 1.5, where strong line features (e.g., Hα, Hβ, OIII, [OII] lines, and stellar absorption features) are within the MUSE spectral coverage.

BH Masses
BH masses are estimated through the single-epoch method (Shen et al. 2011).Assuming the broad-line region (BLR) is virialized, we can measure the BH mass from a single spectrum by estimating the BLR size from the quasar luminosity with the radius-luminosity relation (Bentz et al. 2013) and the BLR virial velocity from the width of a broad emission line (e.g., the 2800 Å MgII line).We extract quasar spectra from apertures of 10 pixels centered on the quasar positions in the MUSE data cube, which roughly encompass 95% of the total light.For the quasars where the MgII line falls into the sodium gap, we use the supplementary quasar spectra from the MagE observations.We fit the quasar spectra with a power-law continuum, the FeII pseudocontinuum, and a series of Gaussian profiles for each broad line using the code PyQSOFit (Guo et al. 2018).The host galaxy light is expected to be very faint com-pared to the bright core for these quasars, so we do not perform host decomposition in PyQSOFit.We include up to three broad components, one narrow component, and an additional wing component to model the MgII line.Finally, we follow the recipe from Shen et al. (2011) to estimate BH masses using the luminosity at 3000 Å and the combined full-width-half-maximum (FWHM) from the MgII line for each quasar.The BH mass uncertain- Note-The table columns are quasar name, number of group galaxy members (including the quasar host galaxy), overdensity, significance of overdensity (σ δ = δg/∆δg), mean velocity of the galaxy group with respect to the quasar, velocity dispersion of the galaxy group, the corrected velocity dispersion (see Section 5), and the estimated halo mass.a The group measurements of J2339 when assuming a bimodal velocity distribution, the group with mean velocity closer to the quasar velocity is reported here.
ties from the single-epoch BH mass estimator with the MgII line is roughly 0.4 dex (Shen et al. 2023), though it might be larger for the CUBS quasars since they are outliers in luminosity compared to the quasar population used in reverberation mapping studies.The virial BH masses and the luminosity at 3000 Å (as functions of redshift) are plotted in Figure 2 and tabulated in Table 1.

Radio Properties
We obtain the radio properties of the quasars by cross-matching with the Rapid Australian Square Kilometre Array Pathfinder (ASKAP) Continuum Survey (RACS) DR1 catalog (Hale et al. 2021).ASKAP surveys the sky south of +41 deg declination at a central frequency of 887.5 MHz with a spatial resolution of ∼15 ′′ (convolved to 25 ′′ for the source catalog).We first match the quasar positions in the RACS DR1 catalog for all sources within 30 ′′ .A total of seven quasars have radio counterparts; two (J0110, J2245) show extended/lobed detection (lobe-dominated), and the rest (J0114, J0248, J0357, J0454, J2339) show unresolved single detection (core-dominated or unresolved).For the lobed-dominated sources, we use the total source flux as the radio flux.We calculate the radio-loudness (R) as the ratio of flux density at 6 cm and 2500 Å in the rest frame, R = f ν,6cm /f ν,2500 (Stocke et al. 1992).The UV flux density f ν,2500 is calculated from the best-fit quasar spectra from PyQSOFit.The rest-frame 6 cm flux density is extrapolated from the RACS radio flux by assuming a spectral index α (F ν ∝ ν α ) of −0.5 .The RACS DR1 catalog requires a 5σ detection for sources to be included in the catalog and has an overall 95% point source completeness at an integrated flux density ≈ 3 mJy.Adopting 3 mJy as the upper limit in radio flux density, we find most of the non-detected sources have an upper limit of R ≲ 10.Typically, quasars with R ≳ 10 are defined as radio-loud quasars (Kellermann et al. 1989).We identify five radio-loud quasars, and the remaining quasars are undetected or radio-quiet.For the quasars above declination of −40 degrees, we also crosscheck with the FIRST and VLASS surveys and find consistent radio properties at 3 cm.The radio-loudness parameter R is tabulated in Table 1.

Group Environment around CUBS quasars
To identify member galaxies in the quasar galaxy groups, we search for galaxies close to the quasar redshift in the MUSE field of view.We calculate the group mean velocity (µ Vel ) and velocity dispersion (σ Vel ) by fitting the galaxy relative velocity within |∆V | < 1500 km s −1 of the quasar redshift with a normal distribution.The  .Velocity distribution of galaxies within each quasar field, with the quasar at ∆ V = 0.The black histogram shows all galaxies around the quasar redshift, and the orange (grey) histograms show the double(single)-line redshift within each group identified by the selection criteria described in Section 4.1.For galaxy groups with more than 5 galaxies, we fit the velocity distribution with a normal distribution to measure the mean velocity µ Vel and velocity dispersion σ Vel .The best-fit normal distributions (black solid lines) are scaled to the peak of the velocity distribution.The bimodal fit for the field J2339 is displayed in grey dashed lines.initial velocity range (|∆V | < 1500 km s −1 ) roughly equals the velocity dispersion of the most massive galaxy clusters.For more robust identification, we re-calculate the new µ Vel and σ Vel after sigma-clipping at 3σ and remove all sources beyond 3σ of the original µ Vel and σ Vel , this calculation is repeated for 2-3 times until the group member identification has converged.Figure 3 shows the MUSE white-light images and the galaxies identified within each quasar group.Figure 4 shows the velocity distribution of each group.All galaxies within |∆V | < 3000 km s −1 of the quasar are shown in open histograms, and the identified group members shown in orange (double-feature redshifts) and grey (singlefeature redshifts) histograms.
The groups include 3-23 galaxies around the quasars (including the quasar host galaxy), indicating the CUBS quasars reside in poor to massive groups.For groups with N gal > 5, we report the final µ Vel and σ Vel of the group.The velocity dispersion ranges from 150-600 km s −1 (except for J2339), which suggests our initial velocity search range is sufficient.While there could still be chance projections in the identified group members, the majority of the galaxies in N gal > 5 groups are projected within the virial radii, estimated with the derived halo masses (see Section 5.1).In addition, most of the identified group members are likely bound in N gal > 5 groups, their relative velocities to the quasars are less than the escape velocities at the corresponding galaxy distances.Removing the few galaxies near the virial radius or above to the escape velocity does not change the estimated velocity dispersion or halo mass significantly.The final group size, mean velocity, and velocity dispersion are tabulated in Table 2.
The velocities of galaxies in the field of J2339 show a bimodal distribution, suggesting there might be two individual galaxy groups perhaps starting or undergoing a merger event.When fitted with two normal distributions, we find two distinct groups of 7 (closer to the quasar redshift) and 6 galaxies with similar velocity dispersion of ≈ 350 km s −1 .While ∼1000 km s −1 (σ Vel of J2339 fitted with a single Gaussian) is still within a reasonable range for a galaxy group/cluster, ≈300 km s −1 is more consistent with other groups and better describes the bimodal velocity distribution.We report the group properties of the group closer to the quasar redshift in Table 2 and plot both velocities in the remaining figures of this paper.

Overdensities
To study the quasar environments, we search for overdensities around the central quasars using the relative velocity distribution.Following Stott et al. (2020) calculate the overdensity δ g using the equation: where N group is the number of objects in each quasar group and N bkg is the expected number of background objects without the presence of any structure.We calculate N bkg by stacking the relative velocity distribution for all quasar fields with a bin size of 3000 km s −1 (i.e., ±1500 km s −1 ), centered on the quasar systematic velocity.The velocity bin size is chosen to ensure all identified group galaxies are within the center bin of the stacked histogram (see Figure 4).The full ∆V range of ± 36, 000 km s −1 shown in Figure 5 roughly corresponds to 0.2-0.3 in ∆z depending on the quasar redshift.
Figure 5 shows the stacked distribution of relative velocity from the central quasars.The median (16th, 84th percentiles) number counts of the background galaxies are N bkg = 9.5 (6.0, 14.3) when excluding the central bin.The stacked δ g across all 15 quasar fields is 12.4 +1.2 −1.1 (i.e., ≈ 11σ), and the individual δ g of each quasar field ranges ≈4-35 (tabulated in Table 2) using the galaxy counts in each field and the average background galaxy count across 15 fields, N bkg = 0.63.The uncertainties of δ g are calculated through the Bayesian approach described in Kraft et al. (1991) for determining the Poisson confidence level for a number count given by a background.Following the convention in the literature, we find 11 quasars reside in >2σ overdense environments (i.e., galaxy groups with ≳ 5 members), and 5 quasars reside in >3σ overdensity (≳ 9 members).When excluding the "single-feature" redshifts, the results remain unchanged; the stacked δ g is 14.9 +1.4 −1.3 (i.e., ≈ 11σ), and the same 11(5) quasars are above the > 2σ(3σ) threshold.While N bkg and δ g depend on the exact choice of the velocity grid, the significance of overdensity (σ δ ) remains consistent with bin sizes of 2500 − 4500 km s −1 .We note that the decline at ∆V > 25, 000 km s −1 in the stacked relative velocity (Figure 5) could be due to the lower completeness in redshift measurement at z > 1.However, the overdensity measurements are not sensitive to the exact choice of the full ∆V range, provided that it is large enough for accurately estimating the background galaxy counts.For our sample, the median background galaxy count remains consistent at N bkg = 9 − 10 with full ∆V ranges of ± 24, 000 − 36, 000 km s −1 .
The MUSE FOV (1 ′ ≈ 500 kpc at z ≈ 1) might not be large enough to survey all galaxies in the galaxy groups/clusters.We find 1−3 additional galaxy group members within a projected distance of 500 kpc around half of the CUBS quasars in the LDSS3/IMACS redshift surveys.These additional galaxies are all found in the fields of the > 2σ overdense quasars and the measured velocity dispersion is consistent within 10% if these galaxies are included.The only two exceptions are the fields J0114 and J0119.The additional galaxy found in the J0114 field lies at the edge of the velocity distribution and thus increased the velocity dispersion by ≈30% when included.For the field of J0119, six additional galaxies were found in the LDSS3 FOV near the quasar redshift, which made J0119 one of the most overdense groups in the sample.The estimated velocity dispersion of the galaxy group around J0119 decreases by ≈20% when including the six additional galaxies.However, the source selection of the LDSS3/IMACS redshift surveys is based on color selection to prioritize galaxies in the quasar foreground and the redshift surveys are much shallower than the MUSE observation.Therefore, it is non-trivial to quantify the selection effect in overdensity estimation, and we choose to focus on the overdensity in the MUSE FOV only.

Stellar Masses of the Group Members
We estimate the stellar masses for the galaxies within the CUBS quasar groups by galaxy spectral modeling with the code BAGPIPES (Carnall et al. 2018(Carnall et al. , 2019) ) using both spectroscopic and photometric observations.We extract the galaxy spectra within a ≈0.6 ′′ radius using mpdaf, i.e., the same spectra as used for the redshift measurements.We extract deep psuedo-photometry in g, r, and i-band by convolving the MUSE datacubes with the g, r, i psuedo-bandpasses (i.e., boxcar functions around 4800-5800, 6000-7500, and 7500-9000 Å, following previous CUBS papers).For bright galaxies (i ≳ 20.5 mag), we supplement with the g, r, i, z, Y -band photometry from the Dark Energy Survey Y3 GOLD catalog (Sevilla-Noarbe et al. 2021).Finally, we include the H-band photometry from our Four Star observations.For the MUSE and H-band photometry, we first perform Source Extractor on the pseudo-r band image and use the pseudo-r band segmentation map to extract photometry from each MUSE and H-band images to ensure the photometry is obtained from the same region.The photometric uncertainties are set to be 10%.
BAGPIPES is a Bayesian fitting code that can model galaxy spectra with spectroscopic and photometric data simultaneously.We model the galaxy spectra with the Kroupa (2001) initial mass function, Bruzual & Charlot (2003) stellar population models, an exponential star formation history, and include nebular emission and dust extinction models from Calzetti et al. (2000).To fit spectroscopic and photometric data simultaneously, we allow zeroth to second order calibrations and a noise scaling parameter that accounts for relative flux calibration as suggested by Bagpipes, e.g., corrections for aperture or template mismatch, underestimation of uncertainties, and wavelength-dependent flux calibration in the spectroscopic data.As a robust check, we tested different star formation history (exponential and double power law) and dust extinction models (Calzetti et al. (2000) and Charlot & Fall (2000)) and found the estimated galaxy masses to be consistent regardless of the model and parameter setup.
Galaxies closest to the central quasar are often heavily contaminated, or even out-shined, by the quasar light.For galaxies very close to the quasars, we perform quasar light subtraction using the methods from Johnson et al. (2018) and Helton et al. (2021) to recover the quasarsubtracted spectra and photometry for the galaxies closest to the quasars.In short, we model the light around the quasars using a linear combination of galaxy and quasar eigenspectra to decompose the quasar and galaxy light.For the H-band image, we use GALFIT (Peng et al. 2002) to perform 2D surface brightness decomposition to remove the central quasar light.Finally, we follow the same procedures described previously to extract MUSE spectra, pseudo-photometry, and H-band photometry from the quasar-subtracted datacubes and images.Of the 11 galaxies that are strongly contaminated by the quasar light, we were able to recover the galaxy photometry and stellar mass for four galaxies.
The derived properties of the galaxies in each quasar group are tabulated in Table 3, along with their rel-Li et al.
ative velocities and projected distances from the central quasar.The estimated galaxy masses are 8.5 < log(M ⋆ /M ⊙ ) < 11.6 , with median log(M ⋆ /M ⊙ ) = 9.9, and the typical mass uncertainty from Bagpipes is ≈0.1 dex, excluding systematic uncertainties.

BH-galaxy-halo relation
In this section, we discuss the relations between the central SMBH, the host galaxy, and the surrounding galaxy cluster/group.Observations of local and intermediate (z < 2) redshift quasars have revealed tight correlations between the SMBH mass and host galaxy mass (e.g., Jahnke et al. 2009;Bennert et al. 2011), luminosity (e.g., Laor 1998;Peng et al. 2006;Decarli et al. 2010), velocity dispersion (e.g., Treu et al. 2004;Woo et al. 2006Woo et al. , 2010)), as well as halo mass (M h ) and temperature (e.g., Ferrarese 2002;Gaspari et al. 2019;Donahue & Voit 2022).The correlation between BH, host galaxy, and dark matter halo can be explained through selfregulated feedback processes, i.e., galaxy growth halts as the cool, star-forming ISM gas heats up or gets expelled to CGM scale when the central SMBH releases energy into the surrounding galaxy and halo.If the BH mass is controlled by the halo binding energy, then (Silk & Rees 1998;Booth & Schaye 2010).Though some works argue against a direct correlation between M BH and M h (e.g., Kormendy & Ho 2013).In addition, through abundance matching and direct measurements of galaxy clusters, we can probe the efficiency of turning baryonic matter into stars using the ratio of galaxy stellar and halo masses.Maximum star formation efficiency occurs around halo masses of typical L * galaxies (M h ∼ 10 12 M ⊙ ), and star formation efficiency decreases at both low and high mass ends due to strong feedback from star formation and AGN activities (e.g., Behroozi et al. 2013;Kravtsov et al. 2018;Wechsler & Tinker 2018).

Galaxy and halo mass estimation
We estimate the bulge (or galaxy) stellar masses of the quasar hosts to be around 10 11 -10 12.5 M ⊙ using the M BH − M ⋆,bulge relation from Kormendy & Ho (2013).For the purpose of this work, the difference between host galaxy mass and bulge mass is negligible.For all 15 quasar fields, the quasar host galaxies are more massive than other group members and encompass ≈ 50%−99% of the total stellar mass.Therefore, the quasar host galaxies are very likely the central galaxies in each of the galaxy groups.Central galaxies of galaxy groups, similar to brightest cluster galaxies (BCGs) in galaxy clusters, hold a special position in structure formation .Velocity dispersion of the galaxy groups as a function of overdensity.The velocity dispersion is only calculated for galaxy groups with more than five galaxies.The radioloud quasars are labeled with red squares around the data points.The grey point shows the overdensity and velocity dispersion of the galaxy group around J2339 when the velocity distribution is fitted with a single Gaussian distribution, and is connected to the two Gaussian measurements with a dotted line.The vertical dashed line shows the 2σ overdensity limit.The velocity dispersion is expected to correlate positively with overdensity, as velocity dispersion traces the dynamical mass of the dark matter halo around the central quasar.However, the correlation is not statistically significant for our sample.
as they are centered in the group/cluster halo (e.g., Lin & Mohr 2004).
Figure 6 shows the relation between the velocity dispersion and the overdensity parameter.The velocity dispersion traces the dynamical mass of the dark matter halo around the central quasar and the surrounding neighbors, therefore the velocity dispersion is expected to correlate positively with the overdensity and the number of group galaxies.Though the correlation is not statistically significant in our sample (the Pearson correlation coefficient is 0.39, with a p-value of 0.26, for the groups with N gal > 5).Assuming the galaxy groups are virialized, the halo mass can also be traced dynamically from the velocity dispersion in the galaxy groups.Following Equation 1 from Munari et al. (2013), we convert the line-of-sight velocity dispersion to halo mass: where A = 1177 and α = 0.364 for their simulations using galaxies as tracers and including AGN feedback.The halo mass here refers to M 200 , the dark matter halo mass within R 200 , where the average density is 200 times the critical density of the Universe.
In the small N gal regime, all velocity dispersion estimators are statistically biased by small number statistics.We calculate the halo mass for groups with more than 5 members using σ Vel and follow the parametric corrections for small N gal from Ferragamo et al. (2020) to derive the unbiased standard deviation σ Vel,corr .The intrinsic scatter of Equation 2 is relatively small (≈ 5%), so the halo mass uncertainty (≈ 0.1 − 0.3 dex) is dominated by the measurement uncertainty of σ Vel (assuming ≈ 10 − 30%).We also note that there are additional biases in the halo mass estimation from incomplete sampling of member galaxies and interloper contamination (Ferragamo et al. 2020), which are relatively small and thus not included in our calculation.The estimated halo mass ranges from 10 13 −10 14 M ⊙ , suggesting the most luminous quasars reside in more massive dark matter halos than the general quasar population.In contrast, observations show typical optically-selected quasars reside in M h ∼ 10 12 − 10 13 M ⊙ (e.g., Hickox et al. 2009).Results from halo occupation modeling and clustering analyses also showed that the median masses of quasar halos are M h = 4.1 × 10 12 M ⊙ for central quasars at z ≈ 1.4, and that halos with M h ∼ 10 13 M ⊙ are ≈ 1σ deviation from the full mass distribution of central quasar halos (Richardson et al. 2012).2023) and convert M 500 derived from the halo gas temperature to M h assuming a Navarro-Frenk-White mass profile with a concentration parameter c 200 ≈ 4. Despite the limited dynamical range in BH mass and large scatter, the BH mass and halo mass for the majority of the sample are roughly consistent with the local observed relations and the predicted M BH ∝ M 5/3 h relation.The J0248 galaxy group appears to be an outlier of the sample, the velocity dispersion is much smaller than other galaxy groups with a similar number of galaxies, and the halo mass ∼ 10 12.1 M ⊙ might be underestimated.The observed velocity dispersion may be much smaller by chance, or J0248 could be an outlier in the BH-galaxy-halo relations, but these hypotheses can not be tested without more data on individual quasars.

MBH − M
The right panel of Figure 7 shows the stellar masshalo mass relation.We note that the host galaxy stellar mass and halo mass are derived independently, the former from the M BH − M ⋆,bulge relation and the latter from group velocity dispersion.We compare our sample to the stellar mass-halo mass relation relations from direct measurements of low-redshift galaxy clusters (Kravtsov et al. 2018) and the abundance matching ansatz (Behroozi et al. 2013).Across all halo mass ranges, the stellar mass-halo mass relation has relatively small scatter in stellar mass (≈ 0.2 dex) at fixed halo masses (Tinker et al. 2013).While the Kravtsov et al. (2018) relation is derived for z < 0.1 only, the stellar mass-halo mass relation is not expected to evolve much with redshift up to at least z ≈ 4 as shown in Behroozi et al. (2013Behroozi et al. ( , 2019)).Despite large uncertainties in the stellar masses inferred from the M BH − M ⋆,bulge relation (e.g., ≈ 0.5 dex considering measurement uncertainty and intrinsic scatter), most of our sample lies above the mean stellar mass-halo mass relation from abundance matching, hinting that luminous quasar hosts may represent outliers in the M BH −M ⋆ and M BH −M h relations or have less reliable virial BH masses.

Overdensity and its dependency on quasar properties
At z < 0.5, Serber et al. (2006) and Karhunen et al. (2014) found that quasars reside in denser environments at hundreds of kpc scale, but the overdensity vanishes at ≈ 1 Mpc scale compared to galaxies with comparable masses.However, using spectroscopic redshifts, Wethers et al. (2022) found quasars tend to reside in moderate groups, and there is no statistical difference compared to the redshift-and stellar mass-matched control galaxy sample at the sub-Mpc scale.At 1 < z < 3, the focus of quasar environment studies shifts to identifying protoclusters using radio-loud galaxies and AGN.For example, Wylezalek et al. (2013) found that the majority of radio-loud AGN resides in richer environment than average, with 55% (10%) of their sample in overdensities at > 2σ (> 5σ) level on the Mpc-scale.In a redshift survey with HST WFC3 grism, Stott et al. (2020) found an overdensity in 8 out of 12 quasar fields in their sample, with a stacked overdensity of ≈ 6σ.Similarly, Trainor & Steidel (2012) revealed a significant stacked overdensity on the scale of 200 kpc around the most luminous quasars at z ≃ 2.7 in the Keck Baryonic Structure Survey (KBSS).In this work, we find that the majority (11 out of 15) of CUBS quasars reside in > 2σ overdensities within a projected distance of 250 kpc (or < 1 Mpc including the LDSS3/IMACS FOV) and has a stacked overdensity of ≈ 11σ, consistent with the overdensity at both lower and higher redshifts in the literature.The median group size is 7 galaxies, roughly ≈ 2σ overdensity above the background.In addition, we show that there is substantial diversity in the environments in which quasars reside.In this section, we explore the dependency of overdensity on various quasar and group properties.
Figure 8 shows the overdensity as a function of the SMBH mass (left panel) and quasar luminosity (middle panel).If the BH mass is a good tracer of the halo mass, as previously mentioned in Section 5.1, we expect a correlation between the BH mass and overdensity.However, possibly due to the limited dynamical range and large uncertainties associated with the virial BH mass and halo mass, we do not see a correlation between overdensity and BH mass in our sample.The Pearson correlation coefficients are 0.21 and −0.13, with p-values of 0.46 and 0.65, for the M BH − δ g and L 3000 − δ g relations, respectively.There is no consensus on whether quasar environments depend on luminosity and BH mass in the literature.Karhunen et al. (2014) find that galaxy number density within 1 Mpc has no dependency on redshift, quasar and host galaxy luminosity, or BH mass for lowredshift (z < 0.5) quasars, and Zhang et al. (2013) find that clustering amplitude does not depend on redshift, luminosity, or BH mass, in SDSS Stripe 82 quasars over a wide redshift range of 0.6 < z < 1.2.On the other hand, while there is no luminosity trend for the entire sample, Shen et al. (2009) found that the brightest 10% of quasars are more strongly clustered than the remaining 90% of quasars, suggesting only the brightest quasars reside in the rarest, most massive groups.Similarly, Shen et al. (2013) finds a weak luminosity dependency in quasar analysis at z ≈ 0.5.The CUBS quasars are no doubt some of the brightest quasars at z ≈ 1 (see Figure 2), and while the overall stacked overdensity is significant, some CUBS quasars remain in moderate-topoor environments (i.e., no overdensity compared to the background).
Many studies have shown that radio-loud AGN/quasars tend to reside in denser regions than their radio-quiet counterparts and radio galaxies (e.g., Wylezalek et al. 2013;Ramos Almeida et al. 2013;Hatch et al. 2014;Stott et al. 2020), but others have found no significant difference (e.g., McLure & Dunlop 2001;Karhunen et al. 2014).Figure 8 (right panel) shows the majority (4/5) of our radio-loud quasars reside in regions of high overdensity, while only roughly half (6/10) of the radio-quiet quasars show > 2σ overdensity.Similar to The orange triangle shows the 3σ upper limit of the radio-loudness.The vertical and horizontal dashed lines indicate the 2σ overdensity and the radio-loud limit.We do not see a correlation between overdensity and BH mass, quasar luminosity, or radio-loudness in our sample.
the literature finding radio-loud AGN in overdense environments, the connection between radio-loudness and overdensity is not straightforward and has significant scatter.In the simple AGN feedback evolution model from Hickox et al. (2009), the central SMBH/galaxy grows mass from an optical/IR bright phase until reaching a critical halo mass around 10 12 − 10 13 M ⊙ , then shifts into a more quiescent phase with intermittent radio activity, which could be triggered by mergers and interactions.In this scenario, denser environments could lead to more mergers and interactions, and thus radio-loud quasars preferentially reside in overdense regions.

CONCLUSIONS
Using a MUSE redshift survey and photometric follow-up observations, we study the group environments of the 15 UV-luminous quasars from the CUBS survey.We first identify galaxy group members around the CUBS quasars and measure the group properties and stellar properties of individual member galaxies.
We find that the CUBS quasars reside in groups of sizes ranging from 3-26 galaxies, with halo masses of 10 13 − 10 14 M ⊙ .Out of 15 CUBS quasars, 11 are located in environments with overdensities greater than 2σ relative to the background (i.e., without the presence of a quasar).The overdensity of the quasar environment does not correlate with BH properties (BH mass, bolometric luminosity, and radio-loudness), but there is a tendency for radio-loud quasars to reside in overdense regions.Despite large uncertainties and scatter, the CUBS quasars deviate from the average stellar mass-halo mass relation, suggesting they could be outliers in the BH-galaxy-halo relations or have less reliable BH masses.
This study shows that quasars reside in diverse environments and that the relations between quasars and their environments are complicated.More deep, spectroscopic surveys around individual quasars spanning wider luminosity and redshift ranges with wide-field IFS like MUSE, are crucial to quantify the relations, as well as the intrinsic scatter in these relations, between quasars and their environments.
We thank the anonymous referee for useful comments that helped improve this work

Figure 1 .
Figure 1.Redshift distribution for the CUBS quasars (black), and the double-feature (orange) and single-feature (grey) galaxy redshifts identified in 15 quasar fields.

Figure 2 .
Figure2.BH mass (top), luminosity at 3000 (bottom), and redshift distribution of the CUBS quasars.The shaded contours show the quasar properties from the DR7 quasar catalog, containing quasars brighter than Mi = −22 with a redshift range of 0 < z < 5(Shen et al. 2011).The error bars include measurement uncertainties from spectral fitting and systematic uncertainties (0.4 dex for BH mass and 0.05 dex for luminosity).The CUBS quasars are roughly ≈ 0.5 dex and ≈ 1 dex higher than the general quasar population in BH mass and luminosity, respectively.

Figure 3 .
Figure3.The white-light images of each quasar field constructed from the MUSE data.The CUBS quasars are located at the plus signs, and the identified quasar group galaxies are labeled (the circles show the galaxies with double-feature redshifts, and the squares show the galaxies with single-feature redshifts) and colored by the relative velocity from the quasar systematic redshift.The horizontal bar on the lower right corner of each panel indicates a physical scale of 50 kpc.

Figure 4
Figure4.Velocity distribution of galaxies within each quasar field, with the quasar at ∆ V = 0.The black histogram shows all galaxies around the quasar redshift, and the orange (grey) histograms show the double(single)-line redshift within each group identified by the selection criteria described in Section 4.1.For galaxy groups with more than 5 galaxies, we fit the velocity distribution with a normal distribution to measure the mean velocity µ Vel and velocity dispersion σ Vel .The best-fit normal distributions (black solid lines) are scaled to the peak of the velocity distribution.The bimodal fit for the field J2339 is displayed in grey dashed lines.

Figure 5 .
Figure 5. Stacked distribution of the relative velocity from the central quasars.The solid (and dashed) line shows the median (16th and 84th percentiles) number count of the background galaxies (i.e., excluding the central bin).
Figure6.Velocity dispersion of the galaxy groups as a function of overdensity.The velocity dispersion is only calculated for galaxy groups with more than five galaxies.The radioloud quasars are labeled with red squares around the data points.The grey point shows the overdensity and velocity dispersion of the galaxy group around J2339 when the velocity distribution is fitted with a single Gaussian distribution, and is connected to the two Gaussian measurements with a dotted line.The vertical dashed line shows the 2σ overdensity limit.The velocity dispersion is expected to correlate positively with overdensity, as velocity dispersion traces the dynamical mass of the dark matter halo around the central quasar.However, the correlation is not statistically significant for our sample.

Figure 7 .
Figure7.BH mass (left) and stellar mass-to-halo mass ratio (right) as functions of halo mass.In the left panel, the lines show the best-fit relations from local observations (blue dotted lines: Equation 4 and 6 from(Ferrarese 2002), blue dotted-dash line:Gaspari et al. (2019)) and the theoretical prediction from Booth & Schaye (2010) (black solid line).In the right panel, the blue solid line shows theKravtsov et al. (2018) stellar mass-to-halo mass relation from direct observations and abundance matching at z < 0.1, and the blue dashed and dotted lines show the relations at z ≈ 0 and z ≈ 1 fromBehroozi et al. (2013).J0248, which lies above the average stellar mass-halo mass relations and the observed cosmic baryon fraction (grey shaded area, Planck Collaboration et al. 2020), could be an outlier in the BH-galaxy-halo relations or have less accurate halo mass estimation (see Section 5.1.2).The typical uncertainties (dominated by systematic uncertainties) are shown in each panel.There is little redshift evolution in the stellar mass-to-halo mass relation.Most of our sample lies above the average stellar mass-halo mass relations from abundance matching and theoretical prediction, hinting that the luminous CUBS quasars may represent outliers in the BH-galaxy-halo relations or have less reliable virial BH masses.
Figure 7 (left panel) shows the relation between the halo mass and the BH mass, along with the theoretical prediction of M BH ∝ M 1.55 h from Booth & Schaye (2010) and best-fit relations from Ferrarese (2002) and Gaspari et al. (2019).For the Gaspari et al. (2019) relation, we follow Voit et al. (2023) and convert M 500 derived from the halo gas temperature to M h assuming a Navarro-Frenk-White mass profile with a concentration parameter c 200 ≈ 4. Despite the limited dynamical range in BH mass and large scatter, the BH mass and halo mass for the majority of the sample are roughly consistent with the local observed relations and the predicted M BH ∝ M

Figure 8 .
Figure 8. Galaxy overdensity as a function of BH mass (left), bolometric luminosity (middle), and radio loudness (right).The orange triangle shows the 3σ upper limit of the radio-loudness.The vertical and horizontal dashed lines indicate the 2σ overdensity and the radio-loud limit.We do not see a correlation between overdensity and BH mass, quasar luminosity, or radio-loudness in our sample.

Table 1 .
Summary of the quasar propertiesNote-The table columns are quasar name, alternative quasar names in the literature, quasar right ascension, declination, redshift, BH mass, luminosity at 3000 , and radio loudness (R).

Table 2 .
Summary of the galaxy groups (within the MUSE field of view) . JIL is supported by the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program.EB acknowledges support by NASA under award number 80GSFC21M0002.HWC and MCC acknowledge partial support from NSF AST-1715692 grants.FSZ acknowledges support of a Carnegie Fellowship from the Observatories of the Carnegie Institution for Science.Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program 104.A-0147.This paper includes data gathered with the 6.5-meter Magellan Telescopes located at Las Campanas Observatory, Chile.

Table 3 .
Summary of galaxy properties in each quasar fieldNote-The table columns are quasar name/galaxy ID, galaxy right ascension, declination, redshift, r-band magnitude from MUSE, B-band magnitude from SED fitting, B − I color, stellar mass, relative velocity to the quasar redshift, angular distance and projected distance to the quasar.

Table 3 .
(Continued) a failed to decompose galaxy and quasar light with PSF subtraction