ALMA Reveals a Large Overdensity and Strong Clustering of Galaxies in Quasar Environments at z ∼ 4

The quasar sample used in this study is the same sample presented by García-Vergara et al. (2019), who observed 17 quasar fields at optical wavelengths to search for LAEs in z ∼ 4 quasar environments. We refer the reader to the aforementioned work for details of the quasar selection strategy. Briefly, the 17 quasars were selected from the Sloan Digital Sky Survey (SDSS; York et al. 2000) and the Baryon Oscillation Spectroscopic Survey (BOSS; Eisenstein et al. 2011; Dawson et al. 2013) quasar catalog (Pâris et al. 2014) to lie within a redshift window of z ∼ 3.862–3.879 (corresponding to δv = 1066 km s⁻¹ at z = 3.87). This redshift window was chosen based on the coverage of the optical narrowband filter used to detect LAEs. The quasar redshifts were determined based on one or more rest-frame UV-emission lines (S iv λ 1396, C iv λ 1549, and C iii] λ 1908) and using the calibration of emission-line shifts from Shen et al. (2007) to estimate the systemic redshift. All the quasars have 1σ redshift uncertainties ≲800 km s⁻¹. We tabulate the quasar positions, redshifts, and i-band magnitudes in Table 1.

The 17 quasar fields were previously observed (Program ID: 094.A-0900) using the FOcal Reducer and low-dispersion Spectrograph 2 (FORS2; Appenzeller & Rupprecht 1992) instrument on the VLT in 30 hours of observing time. Optical deep observations were acquired using the narrow band HeI (λ_c = 5930 Å, FWHM = 63 Å) and the broad bands g_HIGH (λ_c =4670Å, hereafter g) and R_SPECIAL (λ_c = 6550 Å, hereafter R) to detect LAEs around the quasars. The total area observed per field was 6.8 × 6.8 arcmin², with an image pixel scale of 0.25''/pix. The medians of the 5σ limiting magnitudes reached in the fields were 24.45, 25.14, and 25.81 for the HeI, R, and g images, respectively.

A catalog of 25 LAEs within R ≲ 7 h⁻¹ cMpc in the quasar fields is available (García-Vergara et al. 2019). LAEs were selected based on a color selection criteria such that they have a (rest-frame) EW_Lyα > 28Å and a solid detection of S/N ≥ 5.0 in the HeI band. All the LAEs lie approximately within ±1598 km s⁻¹ (corresponding to ± 18.7 cMpc at z = 3.78) from the quasar systemic redshift, as given by the coverage of the used NB filter. Only two out of 25 LAEs lie within a radius of 59'' from the central quasar, which is the area covered by the ALMA observations presented in this study.

Our observations were carried out during ALMA Cycle 7, as part of the observing program #2019.1.00411.S (PI: C. Garcia Vergara). The data were taken in Band 3 (84–116 GHz) between 2019 October 4 and 2020 March 19. Most targets were observed in a single scheduling block with 45–50 minutes on-source, except for SDSS J0040+1706, SDSSJ1138+1303, SDSSJ1205+0143, and SDSSJ1224+0746, which were observed in two scheduling blocks.

The 17 pointings were centered at the optical coordinates of each quasar (specified in Table 1). The spectral setup consists of four spectral windows (SPWs). One of the SPWs (hereafter SPW0), with a 1.875 GHz bandwidth, was centered on the expected observed frequency of the CO(4–3) emission line from the central quasar, given the quasar systemic redshifts reported in Table 1. The 1.875 GHz bandwidth correspond to δv ∼ 6000 km s⁻¹ (or δ z = 0.098) at the observed frequencies in SPW0. The SPW0 was placed into the lower sideband and configured in Frequency Division Mode (FDM), which is usually used for spectral-line observations, and provides a spectral channel width of 0.488 MHz. Channels were then averaged together resulting in a final resolution of 7.8125 MHz (equivalent to ∼25 km s⁻¹ at the observed frequency).

The remaining three SPWs were located at higher frequencies, with a 2.0 GHz bandwidth. One of them (SPW1) was adjacent to SPW0, and the other two (SPW2 and SPW3) were in the upper sideband (at ∼12 GHz from the lower sideband). SPWs 1, 2, and 3 were configured in Time Division Mode (TDM), which is usually used for continuum observations because it optimizes the continuum sensitivity. This provides a spectral channel width of 15.625 MHz.

Most observations were taken in the C43-4 array configuration with a maximum baseline length of 783 m, but several observations were taken in a more extended configuration, with a maximum baseline length of 1231 or 2517 m. Table 2 summarizes the details of individual observing blocks. The shortest baseline was always 15 m, and the maximum recoverable scale ranged between 10'' and 14''. The number of 12 m antennas used ranged between 32 and 49.

Table 2. ALMA Observations and Imaging Summary

Target ID	Dates observed	Nant	Max baseline	Beam size (PA)	σcont	σSPW0	log ${L}_{\mathrm{CO}(4-3)}^{{\prime} }$
			[m]	[arcsec × arcsec], [deg]	μJy beam−1	mJy beam−1	[K km s−1pc2]
J0040+1706	2020 Mar 17	32	1231 m	1.21 × 1.00 (36)	12.8	0.32	9.78
	2020 Mar 18	39	1231 m
J0042-1020	2019 Oct 19	49	783 m	1.42 × 1.19 (−85)	8.9	0.23	9.64
J0047+0423	2019 Oct 14	48	783 m	1.38 × 1.26 (−37)	9.0	0.24	9.65
J0119-0342	2019 Oct 19	43	783 m	1.41 × 1.22 (−78)	8.4	0.23	9.64
J0149-0552	2019 Oct 19	49	783 m	1.39 × 1.19 (−65)	8.8	0.24	9.66
J0202-0650	2019 Oct 19	49	783 m	1.49 × 1.17 (−62)	8.5	0.23	9.64
	2019 Oct 14 a	43	783 m
J0240+0357	2019 Oct 4	48	2517 m	1.24 × 0.91 (−53)	14.1	0.35	9.82
J0850+0629	2020 Mar 2	44	783 m	1.38 × 1.22 (49)	10.1	0.25	9.67
J1026+0329	2019 Oct 19	46	783 m	1.37 × 1.24 (−58)	9.8	0.23	9.63
J1044+0950	2020 Mar 2	43	783 m	1.41 × 1.30 (23)	9.5	0.23	9.64
J1138+1303	2020 Mar 16	43	783 m	1.30 × 1.17 (42)	9.3	0.24	9.66
	2020 Mar 19	45	1231 m
J1205+0143	2019 Oct 19	46	783 m	1.33 × 1.23 (−76)	8.5	0.23	9.64
	2020 Mar 3	43	784 m
J1211+1224	2020 Mar 15	43	783 m	1.47 × 1.20 (−32)	11.8	0.29	9.74
J1224+0746	2020 Mar 14	40	783 m	1.38 × 1.24 (−40)	8.7	0.23	9.64
	2020 Mar 15	40	783 m
J1258-0130	2020 Mar 3	44	783 m	1.46 × 1.18 (−46)	10.7	0.30	9.75
J2250-0846	2019 Oct 18	47	783 m	1.39 × 1.18 (−86)	10.4	0.25	9.67
J2350+0025	2019 Oct 14	48	783 m	1.35 × 1.24 (−50)	9.0	0.25	9.67

Notes. For each target, we list the dates of observations, the number of antennas used, and the baseline range. The synthesized beam sizes are for naturally weighted images for the final data products. We include the rms noise measured in the dirty continuum images and the rms noise per channel (channel width of 7.8125 MHz, equivalent to ∼25 km s⁻¹ at the observed frequency) measured in the dirty spectral cube created from SPW0 (which is targeting the expected CO(4–3) emission). We also report the limiting luminosity of the CO line at 5.6σ (the S/N threshold used in this work), assuming a line width of 331 km s⁻¹ (see details in Section 4.1). Since the rms is measured in the non-primary beam corrected image, the reported limiting luminosity corresponds to the center of the pointing for each field, but this increases with increasing radius.

^a Semi-pass observation.

Download table as:  ASCIITypeset image

The primary beam FWHM is 66'' at 94.5 GHz, and the typical size of the synthesized beam is 14 × 12 (see Table 2).

All data processing and imaging was performed using Casa (Common Astronomy Software Applications package, McMullin et al. 2007), version 5.6.1. We process the data using the standard ALMA pipeline, however, some data sets required further manual flagging of problematic antennas. Specifically, we flagged the antennas DA60 and DA58 in J1258-0130, DV04 and DV11 in J0240+0357, DV19 in J1138+1303 (2020 March 19 observation), and DV17 for J1044+0950. The sources were too faint for a successful self-calibration.

For each quasar field, we used Casa's task tclean to create both dirty-image and clean spectral cubes at native frequency resolution. We use the Högbom cleaning algorithm and recalculate the noise per visibility (setting fastnoise=False). All images were created using natural weighting to maximize the surface brightness sensitivity. For the cleaned images, we used a stopping threshold of 1.5σ without manual masking. We set the pixel size to 02, which ensures that the synthesized beam is sufficiently subsampled. Note that the line search was performed on the dirty cubes, and the clean cubes are only used to extract the detected sources (see Section 3.1).

J0240+0357 was observed with a rather sparse UV-plane coverage, resulting in a noticeably non-Gaussian dirty beam with large sidelobes. To improve the image quality, we reimaged this source using an outer Gaussian taper of 15.

We use the dirty spectral cubes created from SPW0 (the cube at the expected observed frequency of the quasar emission line) without a primary beam correction to compute the rms noise per 7.813 MHz channel throughout the scanned frequency range in a rectangular region that covers most of the region within the primary beam FWHM. The rms noise varies only slightly per channel, (typically ∼3%), so we compute the median rms noise measured in all the channels of the spectral cube for each field and report these values in Table 2.

Finally, we create a continuum image by combining all four SPWs. For each quasar field, we compute the rms noise in the dirty continuum images (without primary beam correction) over the same region as for the SPW0 cube and tabulate the values in Table 2. The achieved rms is in agreement with expectations from the ALMA sensitivity calculator. Note that, in our analysis, we do not subtract the continuum from the data.

We perform a blind search for emission lines in the quasar fields (i.e., on the 17 spectral cubes extracted from the SPW0). Since we do not expect to detect extremely bright sources in our cubes, we prefer to use the "dirty" cubes to perform the line search, allowing us to preserve the intrinsic properties of the noise in the cubes. We use findclumps, an IRAF-based routine originally developed for the blind search of CO lines performed in the ALMA spectroscopic survey in the Hubble Ultra Deep Field (ASPECS; Walter et al. 2016; Decarli et al. 2019), a 3D survey of gas and dust in distant galaxies. We refer the reader to these works for further details about the algorithm implemented in the findclumps routine. Briefly, findclumps performs a top-hat convolution using N_chn consecutive frequency channels of the cube to create a combined image, with N_chn ranging from 3 up to 19 (corresponding to 74–470 km s⁻¹ at z = 3.87) in steps of two channels. For each combined image, the rms is computed ⁹ and a search for sources is performed using SExtractor (Bertin & Arnouts 1996), resulting in one source catalog. The signal-to-noise (S/N) ratio of each source in this catalog is computed using the peak flux of the source and the rms computed for the combined image. We only keep sources with S/N ≥ 3.5 in the catalog. The search is performed not only within the primary beam FWHM but over the entire area covered by our observations, which is given by a circular area with a radius 59'' (at which the telescope sensitivity is 10% of the maximum).

Since the source search is performed repetitively using different numbers of channels, we obtain several individual source catalogs, and we combine these into one final catalog per cube; however, the resulting final catalog contains duplications. For example, if the FWHM of an emission line is equivalent to 19 channels, the algorithm will detect such a line with the maximum S/N in the combined image created by using 19 consecutive frequency channels. However, that source may also be detected in the combined image created by using fewer consecutive channels, although with lower S/N because in this case only part of the total flux is encompassed in the combined image.

Since the S/N of a line is maximized when detected combining a number of channels equivalent to the exact width of the line, we look for all the sources that fall within 15 (which is similar to the size of the synthesized beam in our images) and 650 km s⁻¹ (0.21 GHz at z = 3.87), and we only keep in our final catalogs the source with the highest S/N. This procedure effectively removes the duplications, but we could also be missing some real close pairs of sources (for example mergers). We repeat this line search process on all our cubes, resulting in 17 source catalogs.

The assessment of the reliability of the sources in our catalogs is crucial to distinguish real detections from spurious detections only caused by noise peaks exceeding our S/N threshold. The reliability of a line can be determined by exploring the noise properties of each cube. Specifically, we apply the same line search algorithm described in Section 3.1 to the negative cubes, and we create final catalogs containing (unphysical) negative lines detected with S/N ≥ 3.5. Since negative lines are only produced by noise fluctuations, this catalog can be used to determine the probability that a positive line detection is due to noise as a function of their S/N. Note, however, that this probability also depends on the line width, such that at a fixed S/N, a source with a wider line has a higher probability to be real than a source with a narrower line (see a detailed discussion in González-López et al. 2019).

To estimate the noise distribution of each data cube, and under the assumption of a Gaussian distribution for the noise, we fit a Gaussian (centered at S/N = 0) to the S/N histogram of the negative lines for the different width lines. To improve the statistics, we have grouped detections with different widths together before performing the fitting; specifically, we grouped lines with channels width N_chn 3 and 5, 7 and 9, etc. To better take into account the low-number statistics in the tails of the noise distribution, we perform the fitting using a Poisson maximum likelihood estimator.

We compute the reliability (or fidelity F) of a detection as:

$$\begin{eqnarray}&&F({\rm{S}}/{\rm{N}},{N}_{\mathrm{chn}})=1-\displaystyle \frac{{N}_{\mathrm{neg}}({\rm{S}}/{\rm{N}},{N}_{\mathrm{chn}})}{{N}_{\mathrm{pos}}({\rm{S}}/{\rm{N}},{N}_{\mathrm{chn}})},\end{eqnarray} \tag{ 1 }$$

where N_pos(S/N, N_chn) and N_neg(S/N, N_chn) are the number of positive emission lines and the expected number of negative emission lines, respectively, in each S/N bin. The expected number of negative emission lines is computed using a Gaussian function with the best-fitted parameters found by the fitting process described above.

In some previous works (e.g., Walter et al. 2016; Loiacono et al. 2021), the fidelity has been computed using the cumulative number counts (i.e., N_neg( ≥ S/N) and N_pos( ≥ S/N)) instead of the differential number counts. However, if a cumulative approach is used, the computed fidelity at intermediate-S/N values (for instance at S/N ∼ 4.5) would be completely dominated by the higher-S/N detections (for instance, a detection with S/N ∼ 9.0), resulting in a boost of the fidelity for lower-S/N detections. To avoid this artificial boost of fidelity of intermediate-S/N sources, here we use N_neg and N_pos as the differential numbers of detections in bins of S/N values.

The noise properties varies among the 17 cubes, so we repeat this process for all our 17 cubes individually and determine the fidelity of each detected source in the cubes. For visualization purposes only, we computed the median fidelity over the cubes and show this in Figure 1. Since we compute fidelity using the differential number of detections, the fidelity cannot be constrained in S/N bins where no positive lines are detected. For computing the median, we only use the bins with constrained fidelity over the fields (data points in Figure 1) and we linearly interpolate between the S/N values.

Zoom In Zoom Out Reset image size

The choice of the fidelity threshold to keep a source in our catalogs is crucial and has to be taken carefully. This choice varies between different studies, depending on the aim of the analysis. Decarli et al. (2019) is aimed to constrain the CO luminosity function (used in our work as a reference for the computation of the background number counts as described in Section 4.1), and they choose a fidelity threshold of 0.2. However, they treat the fidelity as an upper limit and use a Monte Carlo simulation to generate several realizations of their final catalog, in which they randomly assign a fidelity between 0.0 and their upper limits, and they only keep sources if this random fidelity value exceeds their threshold of 0.2. This means that lower-fidelity sources are less likely to be kept in the catalog for each realization.

This approach is adequate for measuring the luminosity function because, even if sources with low fidelity (likely to be spurious) are included in some realizations, the overall statistics of the number counts will stay mostly unchanged (as explicitly checked in Decarli et al. 2019). Additionally, low-fidelity sources are mostly composed of faint sources, only affecting the faint end of the luminosity function, which has little impact on the overall fitting for the number counts.

In the case of clustering measurements, secure detections are required, since all the sources contribute with the same weight to the final measurement, no matter their intrinsic flux. Including spurious detections would have a strong impact on the final measurement, especially for sparse samples. For this reason, we prefer to choose a more conservative fidelity threshold of ≥0.8, and we only keep sources fulfilling this criterion in our final catalogs. We nevertheless explore the impact of the fidelity threshold choice on our clustering results in the Appendix.

From our median fidelity computation (see Figure 1), we find that fidelity 0.8 is typically reached at S/N ∼ 5.6, S/N ∼ 5.7, S/N ∼ 5.8, and S/N ∼ 6.0, for sources with line widths 15−19 (371–470 km s⁻¹), 11−13 (272–322 km s⁻¹), 7−9 (173−223 km s⁻¹), and 3−5 (74−124 km s⁻¹) channels, respectively.

Since the findclumps algorithm determines the source S/N using a global rms computed on the primary beam FWHM region of the ALMA pointing, the S/N of the detected candidates located close to the edges of the pointing may be overestimated due to noise fluctuation affecting these regions of the pointing. However, we still could detect reliable sources in these regions, so we include sources within the entire area covered by our observations, but we increase the fidelity threshold up to ≥0.9 for sources located in the edges of the pointing (where the telescope sensitivity is lower than 20% of the maximum—or equivalently, at radii larger than 50.14'' from the central quasar).

Considering all the fields together, our final catalogs contain nine sources (excluding the detection of the central quasar). Note that this is the total number of sources detected in all the SPW0 cubes, but not all of them are necessarily included for the clustering computation (see details in Section 4.1).

We extract the 1D spectra for all the selected candidates from the clean cubes. At the resolution of our observations, the sources are not significantly resolved. Therefore, we perform the spectral extraction on the brightest pixel of the source, and we perform a Gaussian fit with a flat continuum to constrain the line center (ν_line) and the line standard deviation (σ_line). We use the immoments task from the Casa package to create a 0th moment map for each source. For that, the data cube is integrated over a frequency width given by 2.8σ_line centered on ν_line, which recovers ∼84% of the line flux.

We measured the total integrated line flux from the 1D spectra by adding the flux from all the channels within ±1.4σ_line from the emission line. We note that this may result in a slightly different S/N ratio than the reported from the line search algorithm, since the computation is made over a different number of channels in some cases (see the gray area versus the dotted lines in the 1D spectra shown in Figure 2). The fluxes are not corrected by the primary beam response of the telescope. Note that the integrated line flux of the individual objects is not relevant for the clustering analysis, since clustering is based only on the position of the sources. The only requirement is that the total integrated line flux of all the sources should be always greater than or equal to 5.6σ (the S/N threshold used in this work), with σ being the limiting flux of the survey. The limiting flux of the survey is the only relevant flux measurement used in the clustering analysis (see Section 4.1), and it is computed based on the rms of the cubes (see Table 2).

Zoom In Zoom Out Reset image size

In Figure 2, we show the moment 0 maps and the 1D extracted spectra for all the sources in our catalogs, and the best-fitted parameters for the lines are listed in Table 3. In Table 5, we summarize the number of emission lines in each individual field.

Table 3. Properties of the Emission Lines Detected in All the 17 Fields Properties of the Emission Lines Detected in All the 17 Fields

ID	R.A. (J2000)	Decl. (J2000)	νline	FWHMline	S/N	Fidelity	δ θ	$\delta {\rm{v}}$	F	Optical Counterpart
	[deg]	[deg]	[GHz]	[km s−1]			['']	[km s−1]	[mJy km s−1]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
J0042-1020.2	00:42:18.48	−10:20:36.48	95.06	129	6.01	0.90	33.0	−1767.19	78.08 ± 13.75	R > 25.6, g > 26.3
J0119-0342.1	01:20:00.37	−03:42:59.90	95.33	58	5.76	0.85	45.1	−2266.36	48.65 ± 9.14	R > 25.7, g > 26.4
J0119-0342.3	01:20:02.15	−03:42:16.96	93.90	128	5.71	0.85	39.0	2270.94	70.32 ± 14.13	R > 25.7, g > 26.4
J0149-0552.2 a	01:49:08.41	−05:52:41.34	94.79	47	5.91	0.80	31.3	−851.97	42.31 ± 9.73	R > 25.7, g > 26.5
J0202-0650.1 a	02:02:54.72	−06:50:51.88	94.82	173	13.11	1.00	16.1	−829.85	193.75 ± 15.79	R = 24.4, g = 25.7
										Lyman break detection
										with g − R = 1.3
										located at 12
J0850+0629.1 a	08:50:14.22	06:29:49.93	94.80	428	17.60	1.00	11.9	−710.84	556.37 ± 27.58	R > 25.7, g > 26.3
J1211+1224.1 a	12:11:45.06	12:25:06.25	94.65	75	6.26	0.93	54.9	559.31	90.54 ± 16.91	R > 25.3, g > 26.3
										LAE located at 16
J1258-0130.1	12:58:43.81	−01:31:09.35	94.85	299	5.97	0.91	53.1	−1220.26	120.69 ± 26.47	R > 25.7, g > 26.4
J2250-0846.2 a	22:50:50.92	−08:45:40.39	94.43	728	8.50	1.00	32.7	398.63	277.69 ± 34.97	R > 25.6, g > 26.4

Notes. This table contains all the emission lines detected with Fidelity $\ge 0.8$ (or $\ge 0.9$ if they are Located at Radius >5014 from the Pointing Center). Columns (1), (2), and (3) indicate the name of the candidate and the R.A., decl. sexagesimal coordinates. Columns (4) and (5) indicate the central frequency and the FWHM of the detected line as determined by the best Gaussian+continuum fitting. Column (6) shows the S/N of the detection as determined by the line search algorithm. Column (7) shows the fidelity computed using Equation (1). Columns (8) and (9) show the angular distance and the velocity distance, respectively, between the emission line and the central quasar, assuming that the observed line corresponds to the CO(4–3) transition. For the computation of the velocity distance, we use the redshift based on the ALMA detection of the quasar, when it is detected with fidelity ≥0.8; otherwise, we use the redshift as determined from the optical quasar spectrum. Column (10) indicates the (non-primary beam corrected) integrated line flux measured from the 1D spectra presented in Figure 2. Column (11) shows the optical color–magnitude g and R measured in a 2'' aperture in the position of the source (or at the position of a counterpart, if one exists within 2''). The lower limits for the magnitude correspond to 3σ limits. We provide comments about optical counterparts for the source.

^a Included for the clustering measurement.

Download table as:  ASCIITypeset image

One of the sources (J2250-0846.2) is identified by our algorithm as two close (δ θ = 04, δv ∼ 625 km s⁻¹) sources. This could be either a galaxy merger or a galaxy with a two-peaked emission line, with total width FWHM > 700 km s⁻¹. In this study, we consider this as a single object, since even in the merger scenario this would likely trace the same dark matter halo. We study this object further in a forthcoming paper (C. García-Vergara et al. in preparation).

Our candidates span a large range of line widths, ranging from 47 km s⁻¹ to 728 km s⁻¹, with three of the lines showing small FWHM values (<100 km s⁻¹). ¹⁰ To explore whether the existence of such narrow lines is expected and reasonable, we checked the line widths of the 16 CO-emitting galaxies blindly found in the ASPECS survey with fidelity >0.9 (González-López et al. 2019). We find that their FWHMs range from 40 km s⁻¹ to 617 km s⁻¹, including two very narrow lines with 40 km s⁻¹ (with S/N = 7.9) and 50 km s⁻¹ (with S/N = 9.5), respectively. Using HST counterparts for the CO emitters detected in ASPECS, Aravena et al. (2019) suggest that these galaxies are likely very face-on, which would explain the detection of such narrow line widths in CO-emitting galaxy surveys.

As we targeted quasars with known redshifts, which are expected to be surrounded by large overdensities of galaxies, and given that we are tracing relatively small areas in the sky, we assume that all the detected lines correspond to the CO(4–3) transition and that the serendipitous detection of emission lines at other redshifts is negligible. To support this assumption, here we quantify the probability that each of our detected lines is actually a galaxy at a different redshift, and we also use our optical images to look for a counterpart of the detected emission lines.

3.4.1. Probability of Low-z Contaminants

First, we compute the probability that each source in our catalog is a low-z galaxy. For this computation, we first explore what transitions at low-z could be detectable given the frequency setup of our survey, and we compute the comoving volume traced for them in each quasar field. We consider a circular area given by R ≤ 59'' (the same used to look for emission lines) and a redshift coverage given by (i) δv = ±3000 km s⁻¹, corresponding to the width of the whole cube, the same as where we detect our nine candidates, and (ii) δv = ±1000 km s⁻¹, corresponding to the width used for the clustering analysis, where we detect five candidates (see details in Section 4.1). We summarize this information in Table 4. The main possible contaminants in our sample are the CO transitions CO(1–0) at z ∼ 0.2, CO(2–1) at z ∼ 1.4, and CO(3–2) at z ∼ 2.7.

Table 4. Probability of the Presence of Low-z Contaminants in Our Survey

Transition	${z}_{\min }$	${z}_{\max }$	Volume3000	N3000	N1000
			[h−1cMpc ]
(1)	(2)	(3)	(4)	(5)	(6)
CO(1–0)	0.206	0.231	6.486	0.023	0.008
CO(2–1)	1.413	1.461	147.276	0.250	0.083
CO(3–2)	2.619	2.691	253.083	0.103	0.034

Note. For the different transitions, we report the redshift and the comoving volume traced in one quasar field at the redshift of the transition. The comoving volume is given by a circular area with radius R ≤ 59'' and a redshift coverage given by δv = ±3000 km s⁻¹. Columns (5) and (6) show the expected number of lines per quasar field in the comoving volume indicated in column (4), and in 1/3 of that volume (i.e., assuming δv = ±1000 km s⁻¹), respectively. Note that, in the absence of clustering, the number of expected lines only depends on the volume; therefore, the number of lines simply decrease by a factor of three in column (6) compared with column (5).

Download table as:  ASCIITypeset image

To compute the probability, we use the CO luminosity function of the different CO transitions (Decarli et al. 2019) at the corresponding redshift, and integrate it down to the 5.6σ limiting luminosity of our survey (this is the same S/N threshold that is used in this work for the detections of CO emitters; see Section 3.3) to compute the number density of sources for each transition. We multiply this quantity by the traced comoving volume to obtain the expected number of contaminants in one field, and we show the result in Table 4. The limiting luminosity at the center of the pointing is reported in Table 2. Note that the limiting luminosity varies over the ALMA pointing, due to the decrease of sensitivity with the radius, thus we include this in our computation, following the procedure described in Section 4.1.

The probability that each of the lines presented in our final catalog is CO(1–0), CO(2–1), and CO(3–2) is 2%, 25%, and 10%, respectively. If we only focus on the sample used for the clustering analysis composed of five sources, we find that the probability that each of these sources is CO(1–0), CO(2–1), and CO(3–2) is 0.8%, 8.3%, and 3.4%, respectively.

3.4.2. Optical Counterparts

We detect an emission line in 10/17 quasars, which we assume to be the CO(4–3) transition, as they all have an optical counterpart and spectroscopic redshift from SDSS spectra. All 10 lines have fidelity ≥0.8 and are detected with S/Ns ranging from 5.9 up to 24.5. We relax our S/N criterion down to S/N > 4, and we look for possible line candidates for the other seven quasars, restricting the search to detections with line widths ≥124 km s⁻¹. When more than one candidate was found at the expected sky projected position of their optical counterpart (within 3''), we selected the one with higher S/N. We do not impose a requirement in the velocity space position of the quasar counterpart, since the redshift of optical quasars could be associated with larger uncertainties than reported in previous studies. An uncertainty in the optical computed redshift of 2970 km s⁻¹ (one half of the bandwidth) could be still possible. This search results in the detection of all the other seven quasar counterparts, detected with S/Ns ranging from 4.1 up to 5.3, but with low fidelity (<0.4). Therefore, some of these detections could be noise fluctuations.

We compute the offset between the CO(4–3) redshift versus the redshifted based on rest-frame UV-emission lines for all the quasars, and we show our results in Figure 3 and Table 5. For the secure quasar detections, we find a median velocity offset of ∣$\delta {\rm{v}}$∣ = 738 ± 651 km s⁻¹. This is consistent with the typical reported uncertainties of the optical-based redshifts (see Table 1). The only exception is the quasar located in the field J1224+0746, which exhibits a large offset of ∣$\delta {\rm{v}}$∣ = 2386 km s⁻¹ and therefore is located at 1004 km s⁻¹ from the edge of the SPW0. For completeness, we also include the velocity offset computed for the quasar detected with low fidelity in our ALMA observations, but in this study, we only use the redshifts from the secure detections. We present a detailed analysis of the quasar properties in a forthcoming study (C. García-Vergara et al. in preparation).

Zoom In Zoom Out Reset image size

In Figure 4, we show the sky distribution of the detected companion lines around the central quasar, and in Figure 5 we show the distribution of the velocity offsets between the lines and the central quasar for all our fields together. For the computation of the velocity offsets, we use the redshift based on the ALMA detection of the quasar, when it is detected with fidelity ≥0.8; otherwise, we use the redshift as determined from the optical quasar spectrum (reported in Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We measure a volume-averaged projected cross-correlation function between quasars and CO(4–3) emitters defined by

$$\begin{eqnarray}&&\chi (R)=\displaystyle \frac{{\int }_{{V}_{\mathrm{eff}}}\xi (R,Z){dV}}{{V}_{\mathrm{eff}}},\end{eqnarray} \tag{ 2 }$$

where ξ(R, Z) is the real-space correlation function, assumed to have a power-law shape $\xi {(r)=(r/{r}_{0,\mathrm{QG}})}^{-\gamma }$, where r_0,QG and γ are the correlation length and the slope of the correlation function, respectively. The real-space separation between objects, r, has been written in Equation (2) as a function of their two components: the transverse comoving separation R, and the radial comoving separation Z. V_eff is the effective volume of the survey, which is a cylindrical shell centered on each quasar, with inner and outer radius given by ${R}_{\min }$ and ${R}_{\max }$, respectively, and with a height defined by the radial comoving distance ${Z}_{\max }-{Z}_{\min }$.

We compute χ(R) in logarithmically spaced radial bins centered on the quasar by using the estimator

$$\begin{eqnarray}&&\chi (R)=\displaystyle \frac{\langle {QG}(R)\rangle }{\langle {QR}(R)\rangle }-1,\end{eqnarray} \tag{ 3 }$$

where 〈QG(R)〉 is the number of quasar–galaxy pairs observed in our survey, within the cylindrical shell volume, and 〈QR(R)〉 is the number of quasar–galaxy pairs that we expect to observe within the same volume but in the absence of clustering (i.e., assuming the background number density of CO-emitting galaxies at z ∼ 4).

If we had ALMA observations covering large areas of the sky at random locations (i.e., regions not containing quasars), using the same setup as for our data, we could directly measure 〈QR(R)〉 from the data. However, this is not the case. We cannot use the three quasar–free SPWs (SPWs 1, 2, and 3) for this purpose, due to the following reasons. First, the central frequency (at least for SPW1) is relatively close (<6000 km s⁻¹) to the quasar location, and thus the background number counts could not necessarily be reached. ¹² This leaves us with only the other two SPWs, which are far enough away, reaching the background number density, although they would be tracing the number density of CO(4–3) at slightly lower redshift (z ∼ 3.2). Second, the volume covered by SPW2 and SPW3 is small, resulting in extremely low statistics for the number counts. ¹³ This would result in large uncertainty in the computed 〈QR(R)〉, and therefore in the χ(R) value.

Third, the SPW1, SPW2, and SPW3 were observed in TDM mode, whereas the SPW0 was observed in FDM mode (see Section 2.3). This results in different spectral noise properties, as well as different spectral resolutions (the spectral resolution in SPW1, SPW2, and SPW3 is two times lower than the resolution in SPW0), which complicates the comparison of the line number density in the different cubes (we would be sensitive to a different type of lines in SPW1, SPW2, and SPW3 compared to the lines detected in SPW0). Considering all the mentioned reasons, we thus decided to estimate 〈QR(R)〉 using the CO(4–3) luminosity function measured at z = 3.8 from the ASPECS survey (Decarli et al. 2019). We detail this in what follows.

The number of quasar–galaxy pairs that we expect to observe in the absence of clustering is given by $\langle {QR}(R)\rangle ={\int }_{{V}_{\mathrm{eff}}}{n}_{\mathrm{CO}}(L^{\prime} \geqslant {L}_{\mathrm{Lim}}^{{\prime} }){dV}$, and for the cylindrical shell volume defined above, we can write

$$\begin{eqnarray}&&\langle {QR}(R)\rangle ={\int }_{Z\min }^{{Z}_{\max }}{\int }_{R\min }^{{R}_{\max }}{n}_{\mathrm{CO}}(L^{\prime} \geqslant {L}_{\mathrm{Lim}}^{{\prime} })2\pi {RdRdZ},\end{eqnarray} \tag{ 4 }$$

where ${n}_{\mathrm{CO}}(L^{\prime} \geqslant {L}_{\mathrm{Lim}}^{{\prime} })$ is the number density of CO(4–3) emitters with line luminosity above the limiting luminosity of our survey. To compute the limiting luminosity, we follow Solomon et al. (1997),

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{Lim}}^{{\prime} }}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}{\mathrm{pc}}^{2}}=\displaystyle \frac{3.257\times {10}^{7}}{1+z}\displaystyle \frac{{F}_{\mathrm{Lim}}}{\mathrm{Jy}\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\left(\displaystyle \frac{{\nu }_{0}}{\mathrm{GHz}}\right)}^{-2}{\left(\displaystyle \frac{{D}_{{\rm{L}}}}{\mathrm{Mpc}}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

where ν₀ is the rest-frame frequency of the line (ν₀ = 461.04 for CO(4–3)), D_L is the luminosity distance at the redshift of the source, z, and ${F}_{\mathrm{Lim}}$ is the limiting integrated line flux.

The limiting integrated line flux for each quasar field is computed using the rms noise of the cubes from Table 2, as well as the S/N threshold defined by our fidelity threshold. Specifically, we assume a typical FWHM line width of 331 km s⁻¹, which is the median width measured in the high-fidelity CO line sample of the ASPECS survey ¹⁴ (González-López et al. 2019), and we use the corresponding S/N threshold for such a width (S/N = 5.6 for a detection with fidelity 0.8 as shown by the magenta curve in Figure 1). The limiting luminosity of the survey at the center of the pointing is reported in Table 2.

If the rms sensitivity of the survey were roughly flat, the term ${n}_{\mathrm{CO}}(L^{\prime} \geqslant {L}_{\mathrm{Lim}}^{{\prime} })$ in Equation (4) could be taken out from the integral and it could be simply computed by integrating the CO(4–3) luminosity function down to the limiting luminosity of the survey computed from Equation (5). However, since the ALMA sensitivity decreases when the radius within the pointing increases, the limiting luminosity actually depends on R, and therefore we have to model this dependency and numerically compute the integral over R in Equation (4).

We model the sensitivity variation over the pointing as a 2D Gaussian with FWHM given by the primary beam size (∼66''), and compute the primary-beam-corrected limiting integrated line flux as a function of R, which gives us the limiting luminosity as a function of R from Equation (5).

Finally, the limiting luminosity as a function of R is used to compute ${n}_{\mathrm{CO}}(L^{\prime} \geqslant {L}_{\mathrm{Lim}}^{{\prime} })$ at every dR integration step in Equation (4). Here, we use the CO(4–3) luminosity function at z = 3.8 measured by Decarli et al. (2019), which is given by the Schechter parameters $\mathrm{log}{{\rm{\Phi }}}_{* }[{\mathrm{Mpc}}^{-3}{\mathrm{dex}}^{-1}]=-3.43$, ${\mathrm{logL}}_{* }^{{\prime} }[{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}]=9.98,$ for a fixed α = −0.2.

In principle, in Equation (4) we should include a multiplicative term encapsulating the completeness of the sample. This completeness refers to the fraction of sources detectable by the search algorithm, and it is dependent on the line flux, line width, fidelity, and S/N of the sources. ¹⁵ Based on previous estimations of the performance of the line search algorithm findclumps using simulated mock sources (e.g., González-López et al. 2019; Decarli et al. 2020; Loiacono et al.2021), we expect that at the high fidelity and S/N threshold used in our study, the completeness is near 100%, so we have neglected this completeness correction. This could result in a slight underestimation of the clustering of sources around quasars, but our main conclusions will not significantly change.

The computation of 〈QG(R)〉 and 〈QR(R)〉 is performed individually for each of the 17 fields, and then they are stacked to obtain the final χ(R) value using Equation (3). We measure the clustering using ${Z}_{\max }=-{Z}_{\min }=8.19\,{h}^{-1}\,\mathrm{cMpc}$ (corresponding to δv = ±1000 km s⁻¹ at z = 3.87). We note that the quasar with the largest redshift offset is located at 1004 km s⁻¹ from the edge of SPW0 (see Section 2.1), and therefore the choice of δv = ±1000 km s⁻¹ allows us to trace a symmetric and complete volume around the quasar for all the fields. This choice results in a total of five sources included in the clustering analysis (see Table 3). Given the small size of our sample, we assume that Poisson error dominates our measurement, and we compute the one-sided Poisson confidence interval for small number statistics from Gehrels (1986). We show the measured cross-correlation in Figure 6 and tabulate it in Table 6.

Zoom In Zoom Out Reset image size

Table 6. Quasar–CO Cross-correlation Function

${R}_{\min }$	${R}_{\max }$	〈QG(R)〉	〈QR(R)〉	χ(R)
h−1 cMpc	h−1 cMpc
0.096	0.165	0.000	0.008	−1.00${}_{-0.00}^{+226.13}$
0.165	0.282	0.000	0.025	−1.00${}_{-0.00}^{+74.91}$
0.282	0.483	2.000	0.065	${29.89}_{-19.96}^{+40.75}$
0.483	0.827	2.000	0.119	${15.79}_{-10.85}^{+22.15}$
0.827	1.417	1.000	0.063	${14.88}_{-13.14}^{+36.53}$

Download table as:  ASCIITypeset image

Note that choosing a larger radial comoving range (δv) would result in more sources for the measurement, increasing the statistics, but the clustering signal would be diluted because we are integrating the signal in a larger volume and up to large distances from the quasar, where the background is close to being reached. We also note that keeping the volume small also decreases the probability of having a contaminant in the sample (see Section 4). We explore the effect of the chosen δv in the Appendix.

For the seven quasars not (securely) detected in our observations, we use the optical-based redshifts, which could have an offset from the ALMA-based redshift of ∼700–800 km s⁻¹ (see Figure 3). This would mean that, for these quasars, we may be including or excluding companions at larger and smaller redshift distances, respectively. The confirmation of the precise redshift for these quasars is the only way to correct for these uncertainties in our clustering measurement, but we explore how much the clustering signal changes if we use ALMA-based redshifts for all the quasars, even if they are detected at low fidelity. We find that the same five CO emitters are kept in this case, confirming the stability of the obtained correlation function.

To determine the real-space cross-correlation parameter r_0,QG that best fits our data, we use a Poisson maximum likelihood estimator. Given the noise of the data and the small physical scales proved in our study, we follow common practice and fit our data assuming a fixed slope γ = 1.8 (e.g., Ouchi et al. 2004; García-Vergara et al. 2017, 2019; Fossati et al. 2021). We find that the maximum likelihood and 1σ confidence interval for the correlation length is ${r}_{0,\mathrm{QG}}\,={8.37}_{-2.04}^{+2.42}\,{h}^{-1}\,\mathrm{cMpc}$. We use the best-fitted parameter in Equation (2) to compute the corresponding χ(R) value, which is shown as a black line in Figure 6.

Finally, we use our 〈QG(R)〉 and 〈QR(R)〉 binned values to compute the observed and expected cumulative number counts of CO(4–3) lines in our 17 quasar fields (〈QG( ≤ R)〉), and show this in Figure 6. In the whole volume survey (1751.7 h⁻³ cMpc³ for the 17 fields over δv = 1000 km s⁻¹), we find a total of five CO(4–3) lines, while we expect only 0.28 CO(4–3) lines from the background alone, resulting in a total CO(4–3) line overdensity of ${17.6}_{-7.6}^{+11.9}$ in quasar fields. We note that, due to the low-number statistics, the computed overdensity is significant at the 2.2σ level, and either deeper observations or observations of a larger sample of quasar fields would be required to confirm the overdensity with higher significance. We also compute the total overdensity per field as the ratio of 〈QG(R)〉 per field over 〈QR(R)〉 integrated over the radial bins and provide these values in Table 5. Although the individual overdensity is dominated by low-number statistics and affected by cosmic variance, we include this to study possible correlations between the overdensities and the quasar properties (C. García-Vergara et al. in preparation).

We compare our results with the clustering of LBG and LAE around z ∼ 4 quasars. Since the previous LBG and LAE clustering studies extend up to larger scales (R ≲ 9 h⁻¹ cMpc) than those traced in our study (R ≲ 1.5 h⁻¹ cMpc), for this comparison we assume that the small-scale quasar–CO cross-correlation function can be extrapolated toward larger scales following a single power-law shape.

We find that the QSO-CO cross-correlation length is slightly lower, but consistent within error bars with the QSO-LBG cross-correlation length, which is given by ${r}_{0,{\rm{Q}}-\mathrm{LBG}}={9.78}_{-1.86}^{+1.68}\,{h}^{-1}\,\mathrm{cMpc}$ for a fixed γ = 1.8 ¹⁶ (García-Vergara et al. 2017). This suggests that CO emitters and LBGs would inhabit dark matter halos of similar masses at z ∼ 4 (we further discuss this point in Section 4.3).

We can also compare our measurements with the clustering of LAEs around this same quasar sample, thus providing a direct comparison of optical and dusty galaxy populations in these fields. We find that the cross-correlation length for CO-emitting galaxies is ${3.0}_{-1.4}^{+1.5}$ times higher than the cross-correlation length for LAEs (${r}_{0,{\rm{Q}}-\mathrm{LAE}}={2.78}_{-1.05}^{+1.16}\,{h}^{-1}\,\mathrm{cMpc}$ with γ = 1.8; García-Vergara et al. 2019) around quasars, although we caution that the uncertainties are still large, and therefore the cross-correlation lengths are consistent within 2σ error bars. We note that the redshift window traced by both studies is similar (δv = ±1000 km s⁻¹ for CO emitters, and δv = ±1600 km s⁻¹ for LAEs), and thus we do not expect that this discrepancy is the result of a dilution in the signal due to differences in the traced volume. We have explicitly checked this in the Appendix, where we find that the quasar–CO cross-correlation length is only 1.1 times smaller if we assume a δv = ±1600 km s⁻¹ (see Figure 8).

Tracing a different quasar sample at z = 3–4.5, and using deeper optical observations, (Fossati et al. 2021) measure a slightly higher quasar–LAE cross-correlation length of ${r}_{0}={3.15}_{-0.40}^{+0.36}\,{h}^{-1}\mathrm{cMpc}$ at a fixed slope γ = 1.8. ¹⁷ This is still ${2.7}_{-0.8}^{+0.8}$ times lower than the cross-correlation length for CO-emitting galaxies. We note that the physical scale traced in this study (R ≲ 0.6 h⁻¹cMpc) is slightly smaller than that traced in our analysis.

The difference in the clustering of CO emitters and LAEs around quasars is unlikely to be caused by differences in the halo mass hosting both populations (see Section 4.3), therefore we suggest that this discrepancy is related to physical processes affecting the visibility of the LAEs around quasars. We further discuss this interpretation in section Section 5.

The autocorrelation of a galaxy population is a powerful tool, since it can be directly related to the dark matter halo mass in which that population reside (e.g., Cole & Kaiser 1989; Mo & White 1996; Sheth & Tormen 1999). The autocorrelation of CO-emitting galaxies at z ∼ 4 has never been measured before, mainly due to the lack of large and deep surveys of these galaxies at high-z. The largest samples of CO emitters at z ∼ 4 currently available are composed of only a few tens of sources (e.g., Decarli et al. 2016, 2019), which do not provide enough statistics for an autocorrelation measurement. Therefore, up to now it has been challenging to compute precise masses of CO-emitting galaxies and understand how they fit into different evolutionary scenarios for high-z galaxies.

However, under certain assumptions, the cross-correlation between CO emitters and quasars can be used to infer the clustering of CO-emitting galaxies in blank fields, providing initial constraints on the clustering of this population at z ∼ 4.

First, we assume that our small-scale cross-correlation can be extrapolated toward larger scales following a power-law shape given by $\xi {(r)=(r/{r}_{0,\mathrm{QG}})}^{-\gamma }$. Although the autocorrelation function of quasars and galaxies has been found to slightly deviate from a power law toward smaller (≲0.2 h⁻¹cMpc) scales (e.g., Ouchi et al. 2005; Hennawi et al. 2006, but see also Shen et al. 2010), likely due to the transition between the one-halo to two-halo terms, the autocorrelation is typically reasonably well-approximated as a power-law, and thus we assume that the one-halo term does not strongly boost the signal in our measurement.

Second, we assume a deterministic bias model, in which the QSO–galaxy cross-correlation function can be written as ${\xi }_{\mathrm{QG}}=\sqrt{{\xi }_{\mathrm{QQ}}{\xi }_{\mathrm{GG}}}$, where ξ_QQ and ξ_GG are the autocorrelation of quasar and galaxies, respectively. We also assume that ξ_QQ and ξ_GG have a power-law shape with the same slope γ = 1.8. Under these assumptions, the correlation lengths can be related by

$$\begin{eqnarray}&&{r}_{0,\mathrm{QG}}=\sqrt{{r}_{0,\mathrm{QQ}}{r}_{0,\mathrm{GG}}}.\end{eqnarray} \tag{ 6 }$$

Using the quasar autocorrelation length previously reported at z ∼ 4, given by r_0,QQ = 22.3 ± 2.5 h⁻¹cMpc (recomputed from Shen et al. 2007 with a fixed γ = 1.8), and our measured QSO–galaxy cross-correlation length, we obtain that the autocorrelation length of CO emitters at z ∼ 4 is given by r_0,GG = 3.14 ± 1.71 h⁻¹ cMpc.

This measurement is slightly lower than the LBG autocorrelation length at z ∼ 4 (${r}_{0,\mathrm{LBG}}={4.1}_{-0.2}^{+0.2}\,{h}^{-1}\,\mathrm{cMpc};$ Ouchi et al. 2004), and is slightly higher than the LAE autocorrelation length at z ∼ 4 (${r}_{0,\mathrm{LAE}}={2.74}_{-0.72}^{+0.58}\,{h}^{-1}\,\mathrm{cMpc}$, Ouchi et al. 2010). We caution that the uncertainties of our estimation are still large, which makes the CO autocorrelation length consistent within uncertainties with the autocorrelation length of both LAE and LBG.

Interestingly, the reported autocorrelation length of CO emitters is lower than the autocorrelation length of S_{870μ m} > 4.0 mJy submillimeter galaxies (SMGs) detected with ALMA at 1.5 < z < 3 (${r}_{0,\mathrm{SMGs}}={7.7}_{-2.6}^{+2.8}\,{h}^{-1}\,\mathrm{cMpc};$ Stach et al. 2021, but see García-Vergara et al. 2020). If the SMG clustering does not strongly evolve with redshift as suggested by Stach et al. (2021), then this would mean that CO emitters are lower-mass galaxies compared to the population of SMGs typically detected in continuum surveys.

While our measurements are still noisy, they provide a first rough constraint of the clustering of CO-emitting galaxies, allowing us for the first time to locate them within the context of evolutionary galaxy models. Larger and deeper surveys of emitting lines around quasars are still needed in order to constrain the clustering of CO emitters more precisely. Alternatively, large and deep surveys of CO-emitting galaxies in blank fields would offer an independent and more direct constraint of the CO clustering at these redshifts. However, such large surveys are very expensive and extremely challenging because of the required sensitivity and the low number density of CO emitters in blank fields.