© The Authors 2024

1. Introduction

In the current paradigm of galaxy evolution, the circumgalactic medium (CGM) plays a major role in regulating star formation and the activity of active galactic nuclei (AGNs). Loosely defined to be the gas virially bound to a galaxy but beyond the visible stellar disk (Tumlinson et al. 2017), the CGM is believed to be heavily influenced by galactic outflows, heating and chemically enriching the gas reservoir, as well as inflows of cool gas (≈10⁴ K) from the intergalactic medium (IGM). While this cold accretion from the IGM filamentary gas structures, also known as the cosmic web, onto the CGM is ubiquitously seen in state-of-the-art cosmological simulations (e.g., Dubois et al. 2014; Schaye et al. 2015; Springel et al. 2018), it is difficult to directly observe due to the expected low densities.

Constraints on CGM properties (e.g., gas density, temperature, metallicity, and kinematics) have been achieved through line-of-sight absorption studies (e.g., Hennawi et al. 2006; Prochaska et al. 2013; Rudie et al. 2019; Péroux & Howk 2020; Donahue & Voit 2022; Faucher-Giguère & Oh 2023) or by directly observing CGM emission with narrowband filters or long-slit spectra (e.g., Hu & Cowie 1987; Heckman et al. 1991; Weidinger et al. 2005; Hennawi & Prochaska 2013). However, the new generation of sensitive integral field unit (IFU) spectrographs such as the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT) and the Keck Cosmic Web Imager (KCWI; Morrissey et al. 2012) have made it possible to routinely observe the emission of the cool CGM gas around individual galaxies by utilizing, for example, bright nebular lines such as Lyα. Most commonly, high-z quasars (2 < z < 6) are targeted for such studies (e.g., Borisova et al. 2016; Arrigoni Battaia et al. 2019a; Farina et al. 2019; Cai et al. 2019; O’Sullivan et al. 2020; Mackenzie et al. 2021; Fossati et al. 2021), often unveiling large portions of the host galaxies’ CGM and with the bulk of the Lyα emission extending for ∼80 kpc. Indeed, quasars are estimated to reside in relatively massive halos across cosmic times (∼10^12.5 M_⊙; Shen et al. 2007; White et al. 2012; Eftekharzadeh et al. 2015; Fossati et al. 2021; Pizzati et al. 2024; Costa 2024), resulting in a virial radius of R_vir ∼ 100 kpc at z ∼ 3, close to their cosmic noon.

While these AGNs produce copious amounts of ionizing photons, greatly outshining their host galaxies, the exact balance between the powering mechanisms of their Lyα nebulae remains disputed. The suggested possibilities include the combination of recombination radiation after the ionization of the material by the quasar ionizing photons (Cantalupo et al. 2005; Kollmeier et al. 2010; Costa et al. 2022), resonant scattering of Lyα photons from the quasar (Cantalupo et al. 2014; Costa et al. 2022), shocks (Taniguchi & Shioya 2000), gravitational cooling radiation (Haiman et al. 2000; Dijkstra et al. 2006). Additionally, in situ scattering of Lyα photons is likely important in shaping the nebulae morphology (Costa et al. 2022; Arrigoni Battaia et al. 2023a). Notwithstanding these uncertainties, both current observations and simulations emphasize the main role of the quasar radiation in powering the extended emission (Costa et al. 2022; González Lobos et al. 2023; Langen et al. 2023; Obreja et al. 2024).

The most extreme Lyα nebulae known, commonly referred to Enormous Lyα nebulae (ELANe), extend up to ≈500 projected kiloparsec and therefore exceed the predicted virial radius for quasar halos. They have been detected around systems hosting multiple quasars: the Slug nebula, situated between a radio-loud and a radio-quiet quasar at redshift z = 2.3, which extends for 460 kpc (Cantalupo et al. 2014); Hennawi et al. (2015) reported the detection of four quasars at z ∼ 2 embedded in a 310 kpc long nebula; Arrigoni Battaia et al. (2018) presented a Lyα nebula extending over 297 kpc and associated with one bright and two faint quasars at z ∼ 3; and Cai et al. (2018) detected a 232 kpc wide nebula surrounding a z ∼ 2.45 quasar pair. This phenomenon might be explained by the overdensity of UV luminous galaxies able to power nebular emission and the surplus of gas fed into the system through tidal interactions during close encounters, which in turn could also fuel AGN activity.

Subsequently, targeted quasar pair fields often show remarkable extended emission. Arrigoni Battaia et al. (2019b) discovered a Lyα nebula stretching between two z ∼ 3 quasars with a projected separation of ∼100 kpc. Li et al. (2023) targeted four pairs consisting of one bright (g < 18 mag) and one faint quasar at z ∼ 2.4 with the Palomar Cosmic Web Imager (PCWI; Matuszewski et al. 2010) and found a maximum projected extent upward of 200 kpc in three of them, indicative of cool gas outside the CGM being illuminated. A quasar pair at z = 3.23 with a projected separation of 500 kpc observed for 40 hours as part of the MUSE Ultra Deep Field (MUDF) program hosts multiple Lyα nebulae extending up to 100 kpc and preferentially along the connecting line between the quasars (Lusso et al. 2019). This is interpreted as possible evidence of a large-scale gaseous filament connecting the two galaxies.

Given the aforementioned expected quasar halo mass of ∼10^12.5 M_⊙, quasar pairs could indeed be suitable tracers of filament orientations. In the MilleniumTNG simulation (Pakmor et al. 2023), the mean halo mass at the nodes of the cosmic web is 10^12.47 M_⊙ at redshift three, with almost all halos above 10¹³ M_⊙ being identified as a node (Galárraga-Espinosa et al. 2024). Per definition, nodes indicate the endpoint of cosmic web filaments, and therefore, two neighboring nodes are expected to be connected to each other by a gaseous filament.

In this paper, we extend the research of previous works by searching for extended nebular emission around a sample of eight quasar pairs observed with VLT/MUSE. Our work has been conducted in the framework of the survey Quasar Snapshot Observations with MUse: Search for Extended Ultraviolet eMission (QSO MUSEUM; Arrigoni Battaia et al. 2019a), using a very similar strategy. The paper is structured as follows. Section 2 describes the observations and the data reduction. In Sect. 3, we explain how the data are analyzed to unveil the extended Lyα emission. Sections 4 and 5 present the results of this work and discuss them in comparison to previous studies. Finally, Sect. 6 summarizes our findings. Throughout this paper, we assume a flat ΛCDM cosmology with H₀ = 67.7 km s⁻¹ Mpc⁻¹, Ω_m = 0.31, and Ω_Λ = 0.69. At the mean redshift of our sample, z ≈ 3.2, one arcsecond corresponds to roughly 7.7 kpc. Reported quasar magnitudes are in the AB system, and absolute magnitudes are K-corrected relative to redshift two (Ross et al. 2013).

2. Observations and data reduction

The sample was selected from the SDSS DR12 quasar catalog (Pâris et al. 2017) complemented by AGN searches at high-z (e.g., Bielby et al. 2013) as described in Arrigoni Battaia et al. (2019b). In brief, to ensure the physical association of the quasar pairs, targets were selected to be within Δz ≤ 0.03 of each other (or about ≤2000 km s⁻¹) based on their systemic redshift provided by the aforementioned works and to have an angular separation of ≤1 arcmin (∼450–500 kpc). Moreover, the redshift of the quasar pairs was restricted to the range 3 ≤ z < 3.9. The lower limit ensures that the Lyα line is within the wavelength range of MUSE, while the upper limit prevents it from being heavily contaminated by sky emission lines. Although the sample selected in this way comprised 17 pairs covering a wide range of projected distances between the two quasars, only close (∼50–100 kpc) and wide (∼450 kpc) pairs were observed due to weather conditions. In the two widest pairs, only one quasar was targeted in each system due to the large angular separation. One quasar pair from the selected and observed sample (ID 8.1 and 8.2) has already been published in Arrigoni Battaia et al. (2019b). We included it here to increase the sample size, and re-analyzed it with the same procedures for comparison. The full spectra of the quasars are shown in Fig. 1. If the quasar is not in the field of view (FoV) of the MUSE observations, its SDSS spectra is instead shown. In the VLA FIRST survey (Becker et al. 1994), a radio source is detected 0.8″ away from quasar 5.1 with f_1.4 GHz = 1.4 mJy beam⁻¹. No other quasar from this sample is detected in VLA FIRST. Figure 2 displays the redshift and magnitude of quasars in our sample in comparison to previous surveys of both single quasar nebulae and quasar pairs.

Fig. 1. Eight targeted quasar pairs. Each panel shows the spectra for the two AGNs extracted either from the final MUSE data cubes in an aperture of 3″or available from SDSS. Vertical dotted lines indicate the location of key emission lines at the SDSS redshift of the brighter quasar. The spectrum of quasar 1.1 is scaled down by a factor of 10. The mean noise per channel in the MUSE data at the respective wavelength of Lyα is $0.08·{10}^{-17}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{\AA }}^{-1}$.

Fig. 2. Quasar magnitudes and redshift covered by different studies of Lyα nebulae. Single quasars are shown as dots (Borisova et al. 2016; Arrigoni Battaia et al. 2019a; Cai et al. 2019; Mackenzie et al. 2021; Fossati et al. 2021; González Lobos et al., in prep.), while physically associated quasar pairs are connected with a line (our work as purple stars; Cai et al. 2018; Lusso et al. 2019; Li et al. 2023).

The sample presented in this paper was observed from November 2017 to March 2018 as part of the program 0100.A-0045(A) (PI: F. Arrigoni Battaia) in service mode with the MUSE instrument in Wide Field Mode on the VLT telescope YEPUN. Quasar properties and information about the observations are summarized in Table A.1. The observations were conducted at an average seeing of 1.16″ (or ∼8.9 kpc at the mean redshift of the sample) with clear weather conditions or thin cloud coverage. For each target, three exposures with 880 seconds of on source time each were taken, rotated by 90 degrees with respect to each other and a few arcseconds dithering pattern. The only exception is quasar pair 7, for which only two exposures are available.

Each cube covers a wavelength range of 4750–9350 Å with a spectral resolution at the mean wavelength of Lyα of R = 1815 (corresponding to FWHM = 2.81 Å or 165 km s⁻¹) and at a spectral sampling of 1.25 Å and a field of view of approximately 1′ × 1′ with a spatial sampling of 0.2″. The data was reduced using the MUSE pipeline v2.8.7 (Weilbacher et al. 2014, 2016, 2020) as described in Farina et al. (2019) and González Lobos et al. (2023), applying bias-subtraction, flat-fields, twilight and illumination correction, sky-subtraction, wavelength and flux calibration. After running the pipeline, there are still residual sky emission lines affecting the data. Therefore, we applied an additional skyline subtraction using the Zurich Atmospheric Purge tool (ZAP, Soto et al. 2016), minimizing both the residuals of sky emission lines and the influence of this procedure on the flux of emission at the expected wavelength and position of the Lyα glow. Subsequently, we masked artifacts due to the edges of the MUSE integral field units via visual inspection and determined the offset between individual exposures by centroid fitting of a bright point source in the FoV. Then, we median combined the science exposures and variance cubes for each object. To determine a more accurate estimate of the noise by taking into account the correlated noise introduced by the pipeline interpolations, we rescaled the associated variance cubes layer by layer to the RMS spectrum of a sigma-clipped background region of the science cubes as usually done in the literature (e.g., Borisova et al. 2016; Bacon et al. 2017). This final variance cube was subsequently used to calculate the surface brightness limit (SBL) and errors. The average SBL in a 30 Å pseudo narrowband (PNB) image centered on Lyα and after masking continuum sources via sigma-clipping is 1.54 ⋅ 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² in a 1 arcsec² aperture, while the average SBL in the same aperture and in the respective central slice is 3.15 ⋅ 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻². This is comparable, but lower than the layerwise SBL reported in Arrigoni Battaia et al. (2019a) (4.4 ⋅ 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²) and Borisova et al. (2016) (5 ⋅ 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²). Both surveys have comparable observational setups, but include more targets at z ∼ 3.3, where sky emission lines can lower the sensitivity at the Lyα line.

3. Data analysis

To reveal the extended emission, the quasar point spread function (PSF) has to be subtracted from the science data cube as it outshines the extended emission (Møller et al. 2000; Christensen et al. 2006; Husemann et al. 2014). Here we used the same empirical method described extensively in González Lobos et al. (2023). Specifically, the size of the PSF and therefore the required subtraction window are dependent on the seeing, which we estimated from the quasar continuum emission. We created narrowbands of 25 Å width and fitted a 2D Moffat function onto the quasar profile while conservatively excluding spectral regions where extended emission might occur, that is, slices around the quasar emission lines typically spanning around 500 Å. Extrapolated to the Lyα wavelength at the systemic redshift of the respective quasar, the Moffat FWHM provides a good estimate of the seeing of the observations. In order to obtain more robust results, we performed this estimation on the brighter quasar in the field, although the results are consistent for all pairs if performed individually. Additionally, we determined the seeing at the wavelength of Lyα using a star in the FoV, if available. Especially for faint quasars, the host galaxy contribution can be high enough that the quasar no longer appears as a point source, leading to an overestimation of the seeing. However, the seeing determined from a star is on average only 0.3% smaller than the one determined from the respective quasar, with the biggest deviation being 6%.

In previous works the PSF subtraction was performed in circular windows with a diameter ranging from 6–10 times the seeing depending on the luminosity of the quasar (Borisova et al. 2016; Mackenzie et al. 2021), with brighter quasars requiring a larger subtraction region. Given the diversity of luminosities in our sample, we accordingly varied the window size between objects in a similar range as previous studies (circular aperture with diameter 4–9 times the seeing), aiming to have virtually no contribution from the quasar PSF to the flux value at the edge of the circular region. If the chosen window size is too big, oversubtraction can occur, and sometimes lead to spurious emission detection in the PSF and the continuum subtracted data cube if there are absorption features in the Lyα line of the quasar.

We constructed an empirical PSF layerwise by creating narrowbands in the subtraction window and scaling it to the flux value in the central 1″ × 1″. In this normalization area, residual values are very error-prone and we therefore masked the pixels for subsequent analysis (as usually done in similar analyses, e.g., Borisova et al. 2016). To avoid including possible extended nebular emission in the PSF construction, wavelength slices containing quasar emission lines were skipped when constructing narrowbands and instead, the next narrowband redward of the line was used to create the PSF. This is especially important for the Lyα line as the narrowband at shorter wavelengths than the line can be heavily affected by absorption. The width of narrowbands has to be varied depending on the magnitude of the quasar in order to achieve sufficient signal-to-noise ratios to constrain the PSF shape, with a typical width of 350 channels or 437.5 Å for this sample. We subtracted the empirical PSF constructed in this way from the respective quasar in the data cube to reveal the large scale emission around it. Subtraction of two PSFs was not necessary in every field: In the quasar pairs 4 and 5, only the brighter quasar is within the MUSE FoV due to the large angular separation, and quasar 1.2 is very dim with an r-band magnitude of 24.6 at the time of discovery with VLT/VIMOS (Bielby et al. 2013) and is detected as a faint continuum source in the MUSE observation (i-band 22.75). Therefore, we subtracted only one PSF in these fields.

Also, two systems (quasar 5.1 and 3.2) required “non-standard” PSF subtractions. In quasar 5.1, continuum sources very close to the PSF disturb the subtraction and we therefore masked them before constructing the empirical PSF. Examining the nature of these sources will be the focus of a companion paper Herwig et al. 2024. The PSF of quasar 3.2 has a very asymmetric spectral shape, leading to severe oversubtraction in one half of the subtraction window when using our standard approach, with a mean value within a 30 Å narrowband around the Lyα peak equivalent to −2.3σ, and undersubtraction in the other half. This is likely due to a very bright narrow nebular line compared to the fainter and broader Lyα from the quasar. Therefore, we masked the quasar spectrum in the slices where nebular emission is suspected and replaced the values with an interpolation. During the construction, we then scaled the PSF by the masked spectrum. The mask encompasses 12 slices, or 15 Å, and was chosen to be the narrowest mask capable of removing the oversubtraction.

Lastly, we utilized the ZAP package to estimate the continuum using the function contsubfits, which applies a median filter. The obtained continuum cube was then subtracted from the PSF-subtracted data cube, yielding a data cube that should only contain nebular emission lines and noise. In the subsequent analysis, continuum sources, excluding the quasars, were masked.

4. Results

In this section we report all the observational results of our analysis. We discuss the most important implications in Sect. 5.

4.1. Surface brightness levels of extended Lyα emission

As described in the previous section, the PSF and continuum subtracted data cubes should only contain nebular emission lines and noise. In this work we only focus on the extended Lyα emission surrounding the quasar pairs. To facilitate comparisons with previous studies, for the subsequent analysis we chose PNBs with a width of 30 Å in observed frame around the Lyα nebula line. To select the central wavelength of the PNBs, we proceeded as follows. After constructing a PNB centered on the peak of the nebula emission line determined by visual inspection, we fitted a Gaussian function to the spectrum extracted within the largest 2σ isophote close to the respective quasar, determined from the SBL in this PNB. We then constructed a new PNB centered on the peak of the fitted line and iterated the process. We selected those that maximize the flux of the fitted line as optimized PNBs. Table 1 summarizes the properties of such optimized PNBs for each pair, and Appendix B provides an example of a PNB before PSF subtraction as well as the corresponding example of SB profile extracted from the PNB before PSF subtraction, after PSF subtraction, and of the empirical PSF.

To account for narrower or wider lines and different nebula redshifts in cubes with two quasars in the FoV, we performed the same iteration individually for each quasar while also changing the width of the PNB to two times the FWHM of the Gaussian fit. This does not significantly change the resulting PNB or sensitivity, as the narrowbands optimized in this way are mostly already 30 Å wide, encompassing both nebulae. The only exception is quasar pair 2 (Fig. A.1), where the determined peak wavelengths λ_Neb are offset from each other by 7.5 Å. The obtained PNBs for all quasar pairs are shown as surface brightness maps in the first panel of Figs. 3 and A.1–A.7.

Fig. 3. Lyα maps and profiles for quasar pair 1 (ID 1.1 and ID 1.2). All maps have been smoothed by a 2D Gaussian kernel with a standard deviation of 2 pixels or 0.4″ and were obtained from the PSF-subtracted optimized PNBs. Shown is a 26″ × 26″ (or about 200 × 200 kpc) cut-out centered between quasar 1.1 and quasar 1.2. Top left: Lyα surface brightness map. The position of the brighter quasar is indicated with a light blue, diagonally dashed circle with a diameter of 1 arcsec, the dimmer quasar is indicated by the darker blue, horizontally dashed circle with the same diameter. Contour lines indicate a detection of 2σ (dark gray), 4σ (light gray) and 10σ (white). The surface brightness maps of all pairs are shown with the same colorbar range to ease comparison. Top middle: Velocity offset with respect to the PNB center of the Lyα line within the 2σ isophote, calculated within the PNB wavelength range. Quasar positions are marked with circles color-coded by their systemic redshift on the same colorbar scale as the map. Top right: Velocity dispersion of the Lyα line within the 2σ isophote. Bottom left: Spectrum of the brighter (QSO1) and dimmer (QSO2) quasar integrated within an aperture of 3″ and spectrum of the nebula integrated within the 2σ isophote indicated in the surface brightness map. The spectra are centered on the Lyα peak of the nebula chosen as the PNB center (black dashed line). The slices used for the PNB are marked in gray. All spectra are smoothed with a Gaussian kernel with a standard deviation of 1 pixel or 1.25 Å. The spectrum of quasar 1.1 is scaled down by a factor of 10. Bottom right: Cosmologically corrected surface brightness profiles calculated in logarithmically spaced annuli around the brighter quasar. Quasar positions are indicated by the dots, color coded by Mi(z = 2). The profile represented by red dots is calculated in the full PNB after masking continuum sources and subtracting the background, the errorbars in y-direction show the 1σ error. For comparison, the median stacked profile of all pairs is plotted as black dashed line. The green curve is calculated within the 2σ isophote indicated in the surface brightness map, with shaded regions showing the 1σ error calculated in the same annuli in the variance cube.

Additionally, we computed the mean surface brightness in logarithmically spaced annuli around the brighter quasar in each pair, extending out to 450 comoving kiloparsecs or 110 physical kiloparsecs at the typical redshift of our sample. During this profile calculation, residual background emission was subtracted by calculating the mean surface brightness within a sigma-clipped background region of each image. This calculation was performed both for the full PNB and for the PNB after masking regions outside the 2σ isophote shown in Figs. 3 and A.1–A.7. The second set of profiles are typically brighter, as they do not include regions with no signal in the calculation.

To get the typical surface brightness profile of physical quasar pairs, we stacked the individual profiles, putting the brighter quasar in the pair to 0 ckpc and taking the median value for each radius (Fig. 4, Table 2). We made this choice because we always find brighter emission closer to the brighter quasar of the pair even though in quasar pair 3 the two AGNs have similar magnitudes and one does not show associated extended emission. For comparison, the median surface brightness profiles of other z ∼ 2–3 quasar samples are shown. These can typically be described with a power law SB_Lyα ∝ r^α with α ≈ −2 (e.g., Borisova et al. 2016; Arrigoni Battaia et al. 2019a). However, in the case of quasar pairs, the surface brightness starts out with relatively low values close to the brighter quasar and falls off more gradually. A fitted power law to the profile weighted by the symmetrized 25th and 75th percentiles yields α = −1.57 ± 0.17. We performed fitting using the nonlinear least squares function curve_fit in scipy (Virtanen et al. 2020). As the slope might be influenced by the presence of dense gas in the host halo of the companion quasar, we repeated this analysis by only stacking the three pairs at wide separation. The obtained profile is not contaminated by a second quasar halo and is best fitted with a power law index of α = −1.17 ± 0.28 for the median stacked profile and α = −1.33 ± 0.16 for the mean profile.

Fig. 4. Cosmology corrected, circularly averaged Lyα surface brightness profiles of different quasar samples. The obtained average Lyα SB profile for z ∼ 3 quasar pairs is compared with (i) samples of individual z ∼ 3 (bright quasars, Arrigoni Battaia et al. 2019a; Fossati et al. 2021; faint quasars, Mackenzie et al. 2021; González Lobos et al., in prep.) and (ii) with z ∼ 2 quasar pairs (Li et al. 2023).

4.2. Lyα nebulae spectra and kinematics

We obtained a global nebula spectrum for each detection by adding the emission from spaxels within the largest 2σ isophote in the PSF-subtracted surface brightness maps close to the respective quasars. These spectra, together with the corresponding quasar spectra extracted within a 3″ aperture centered on each quasar, are displayed in the fourth panel of Figs. 3 and A.1–A.7. As the area over which the nebula spectra are integrated can be quite large and the targeted quasars are relatively faint, the peak of the nebular emission line is in some cases brighter than the peak of one of the quasar’s emission at the nebula peak wavelength (i.e., quasar pair 1, 3, 4, and 8).

From these Lyα nebula spectra we can obtain a first probe of the cool gas kinematics, even though such integrated spectra encompass large portions of the quasars’ surroundings, mixing kinematics on different scales. We calculated the second order of the nebular emission line, σ_Line, in a window of 2 × FWHM obtained from a Gaussian fit to each spectrum. Errors were obtained by resampling the spectra with the associated noise spectra 10 000 times and calculating the second moment. The 32nd and 68th percentiles of the obtained distribution were assumed as lower and upper error. The mean value of σ_Line of all nebulae, 434 km s⁻¹, is indicative of a relatively quiescent cool CGM gas, as opposed to violent kinematics with typical linewidths and velocity shifts above 1000 km s⁻¹ (Villar-Martín et al. 2003).

To further analyze the Lyα kinematics in each nebula, we obtained the line velocity shift and velocity dispersion by calculating the first and second moment from the subcubes (i.e., wavelength range) used to build the PNBs. In particular, for each pair, we set the central slice of the PNB wavelength range as reference velocity of the nebula and computed the moments within the 2σ isophote. The obtained maps are shown in the second and third panels of Figs. 3 and A.1–A.7. We find that the detected nebulae have velocity shifts of a few hundred km s⁻¹, and average velocity dispersions of 390 km s⁻¹.

4.3. Lyα nebulae morphology

To quantify the morphology of the detected nebulae, we determined multiple metrics using the largest 2σ isophote associated with the respective quasar: the integrated Lyα nebula luminosity L_Neb, the maximum extent of the nebula d_max as well as the maximum extent between quasar position and nebula 2σ isophote, d_QSO, max (visualized in Fig. A.3), the offset between quasar and nebula centroid d_{QSO − Neb}, the enclosed area A_2σ, the ratio between major and minor axis α, and the angular offset between the nebula semi-major axis and the vector connecting the two quasars in the pair, ϕ. To determine α and ϕ, as described in Arrigoni Battaia et al. (2019a), we first calculated the second order moments of the flux distribution as

$$\begin{array}{cc}\hfill M_{\mathrm{xx}}& =\langle \frac{{(x-x_{\mathrm{Cen}})}^2}{r^2}\rangle ,M_{\mathrm{yy}}=\langle \frac{{(y-y_{\mathrm{Cen}})}^2}{r^2}\rangle ,\hfill \\ \hfill M_{\mathrm{xy}}& =\langle \frac{(x-x_{\mathrm{Cen}})(y-y_{\mathrm{Cen}})}{r^2}\rangle ,\hfill \end{array}$$(1)

with the flux-weighted centroid coordinates (x_Cen, y_Cen) and the distance of the point (x, y) to the centroid, r. From this, the Stokes parameters follow as

$$\begin{array}{c}\hfill Q=M_{\mathrm{xx}}-M_{\mathrm{yy}},U=2M_{\mathrm{xy}}\end{array}$$

and were used to calculate the axis ratio, or asymmetry of the nebula,

$$\begin{array}{c}\hfill \alpha =\frac{1-\sqrt{Q^2+U^2}}{1+\sqrt{Q^2+U^2}}\end{array}$$

and the angle of the semi-major axis to the next x- or y-axis,

$$\begin{array}{c}\hfill \gamma =\arctan \left(\frac{U}{Q}\right).\end{array}$$

From γ, the angle ϕ between semi-major axis of the nebula and the connecting line between the quasars can be determined. The determination of positions and as a result most aforementioned values are seeing-limited and this imprecision accounts for the majority of uncertainties. We therefore assumed an error of 1 × seeing on all distance measurements and obtained errors for ϕ and α by randomly varying the centroid position 500 times within a circle of radius 1 × seeing. The minimum and maximum values of the obtained distribution were taken as upper and lower errors. Physical area measurements are also influenced by uncertainties on the redshift. A conservative estimate for this error is 10 percent of the measured area (Arrigoni Battaia et al. 2023a). To avoid taking into account the same nebula twice, we only determined morphologies with respect to the brighter quasar in the pair if both of them are associated with the same nebula or only one quasar is in the FoV (ID 1.1, 4.1, 5.1, 6.1 and 8.1). If there are two clearly distinct nebulae or only one nebula, analysis was performed for each of those individually with respect to the quasar they are associated with (ID 2.1, 2.2, 3.2, 7.1). This analysis was not performed for quasar 7.2, as we require a minimum enclosed area of 2 arcsec² to be able to sensibly constrain the morphology. Table 3 lists all the aforementioned metrics computed to constrain the morphology of the detected Lyα nebulae. On average, the extended emission spans about 90 kpc with a nebula luminosity of 2.8 × 10⁴³ erg s⁻¹.

Fig. 5. Lyα nebulae axis ratio versus their luminosity. The values for the nebulae discovered around the studied quasar pairs (red stars, accompanied by the quasar ID) are compared with those for nebulae found around bright (cyan dots and cyan dash-dotted line; Arrigoni Battaia et al. 2019a) and faint (dotted gray line; Mackenzie et al. 2021) z ∼ 3 quasars. The cyan star indicates the value for the Fabulous ELAN in Arrigoni Battaia et al. (2019a), which is hosted by an AGN triplet (of which two are type-I quasars as in this work). The cyan shaded region indicates the 25th and 75th percentiles of the distribution of bright quasars.

Fig. 6. Angle between the Lyα nebula semi-major axis and the line connecting the two quasars, ϕ versus the maximum distance between the bright quasar in each pair and its associated Lyα nebula 2σ isophote, dQSO, max. The data points are color coded following Mi(z = 2) of the associated quasar, while their sizes indicate the projected distance between the quasars in each pair. The quasar IDs are displayed next to the respective data points.

We summarize our main findings on nebulae morphologies in Figs. 5 and 6. The first figure shows the axis ratio of each nebula, α, against its nebula luminosity L_Neb, with the dashed red line indicating the median value of α_Median = 0.60. For comparison, we also provide the values reported in Arrigoni Battaia et al. (2019a) for nebulae detected around more luminous single z ∼ 3 quasars (QSO MUSEUM I) together with their median value α_Median = 0.71 (turquoise dash-dotted) and their 25th and 75th percentiles, and the median value of α found for nebulae associated with dimmer z ∼ 3 quasars (dotted gray line, Mackenzie et al. 2021). For a better visualization of the data, the side plots of Fig. 5 provide the distributions of all the points as histograms.

We find that the second least luminous nebula considered for the morphological calculations, ID 2.1, has the most circular shape. The rest of the sample can be divided into two groups: The first one falls well within the distribution of QSO MUSEUM I with moderate asymmetries (α ∼ 0.64), while the second group shows exceptionally high asymmetries (α ∼ 0.36). Spanning almost two orders of magnitude in nebula luminosities and encompassing quasars 2.2, 7.1, and 8, and therefore a multitude of projected separations and magnitudes, there is no clear trend in what drives these asymmetries. Overall, the shape of quasar pair nebulae is more asymmetric than single quasar nebulae with a median value below the 25th percentile of bright quasar nebulae, even though in previous studies, there is a trend of increasing α with lower quasar luminosity. In other words, we find that quasars with luminosities as low as those in Mackenzie et al. (2021) but part of a pair are associated with more asymmetric nebulae. The integrated nebula luminosity L_Neb spans a similar range in both samples, but peaks at a higher value for single quasars, with the brightest nebula having a luminosity of 1.8 × 10⁴⁴ erg s⁻¹, and a higher average luminosity of 4.5 × 10⁴³ erg s⁻¹. Such more asymmetric morphologies in quasar pairs could be due to the presence of intergalactic structures connecting the two quasars in each pair. This tentative evidence of more asymmetric morphologies needs to be confirmed by a larger sample of quasar pair observations.

To investigate this further, Fig. 6 displays the relation between d_QSO, max (the maximum distance between the bright quasar in each pair and its associated nebula 2σ isophote) and ϕ (the angle between the nebula semi-major axis and the line connecting the two quasars). In the case of isolated/unrelated quasars we expect to find a random distribution of ϕ angles and no trend with distance. Instead, the two values are correlated with a Spearman correlation coefficient r = −0.7833 and a chance of coincidental correlation of p = 0.0125. We further tested the null hypothesis by randomly orienting the nebulae 50 000 times and calculating the Spearman coefficient for the obtained distributions. In 0.7% of realizations, the Spearman coefficient is at least as significant as for the data distribution, that is, below −0.7833, and in 0.2% of realizations, the 4 most extended nebulae all show alignment angles below 20 degrees.

The correlation is not driven by the magnitude of the quasars (color of the dots) or the projected distance between the pair (size of the dots), although at lower projected distances, the offset angle ϕ tends to be smaller. Possibly due to the small sample size, the significance of this trend is low, with a p-value of 0.5 determined by performing a 2D Kolmogorov–Smirnov test (Peacock 1983; Fasano & Franceschini 1987; Press et al. 2007) using the public code NDTEST¹.

In total, 4 out of 9 nebulae (∼45%) have ϕ ≲ 20 degrees, and all of these 4 nebulae extend toward the other quasar in the pair. Due to the small errorbar on large nebulae, this main finding holds true even when considering the upper limit for the alignment angle. We therefore expect that the Lyα nebulae with the largest d_QSO, max and smallest ϕ trace intergalactic bridges possibly connecting the two quasars.

4.4. Notes on each quasar pair

4.4.1. Quasars 1.1 and 1.2

Quasar 1.1 (Fig. 3) is the brightest object of the sample; it is bright as the quasars targeted in Borisova et al. (2016) or Arrigoni Battaia et al. (2019a). It is also associated with the brightest Lyα nebula with a luminosity L_Neb ∼ 9 × 10⁴³ erg s⁻¹. Accordingly, its Lyα surface brightness profile is an order of magnitude brighter than the median profile of all quasar pairs. Quasar 1.2, on the other hand, is the faintest targeted object and sits at a projected distance of 47.7 kpc from quasar 1.1. The Lyα nebula shows a slight elongation between the two quasars, also indicated by the low value of ϕ ≈ 2°, and the nebular line velocity shift evolves from ∼220 km s⁻¹ around quasar 1.1 to ∼−400 km s⁻¹ close to quasar 1.2. However, this is not reflected in the systemic redshifts of the AGNs (see top-middle panel in Fig. 3), likely due to the associated large uncertainties on these values as only broad lines are available for the redshift determination, leading to typical errors of up to 400 km s⁻¹ (Shen et al. 2016). The extended Lyα emission traces rather quiescent kinematics with an average velocity dispersion of 385 km s⁻¹.

4.4.2. Quasars 2.1 and 2.2

Quasar pair 2, displayed in Fig. A.1, is associated with the smallest Lyα nebulae in both area and extent. The extended Lyα emission around quasar 2.1 is also the only one in the sample that is almost circularly shaped according to the asymmetry parameter α, and therefore, the value of ϕ (evaluated here to be about 40°) is not well constrained. The two nebulae do not seem to be physically connected, supported by their redshift, which is close to the systemic redshift of the respective quasar. This might indicate that, although the projected separation is relatively small (∼100 kpc), the quasars could be at a significant distance from each other. In this framework the pair would be stretched along our line of sight out to a distance of 2.7 Mpc obtained by assuming that the velocity difference between the two quasars (∼870 km s⁻¹) is all due to the Hubble flow. Alternatively, the two quasars might not be able to illuminate gas in the transverse direction due to obscuration. In other words, their ionization cones do not illuminate the gas in between the two objects. While quasar 2.1 and 2.2. are relatively faint, there have been observations of individual quasars of similar luminosities that show similar or more extended associated nebulae (e.g., compare ID4, ID5, ID6 and ID7 in Mackenzie et al. 2021).

4.4.3. Quasar 3.1 and 3.2

Quasar 3.1 and 3.2 have similar absolute magnitudes M_i(z = 2) =  − 25.16 and are the second faintest objects (after ID 1.2) observed in this study, but are associated with the third-biggest nebula, which is displayed in Fig. A.2. We note that this is true even when taking into account the slightly differing PSF-subtraction method used for this source, which impacts the estimated A_2σ. Indeed, the standard subtraction of an empirical PSF leads to artifacts close to the quasar, reducing the area enclosed by the 2σ isophote, but does not influence different measures of nebula extent like d_max, which are sensitive to large scale emission. The discovered extended emission does not stretch between the quasars, akin to other pairs in the sample, but is concentrated around quasar 3.2 and shows relatively quiescent kinematics ($\overline{\sigma }=388$ km s⁻¹) and a shallower surface brightness profile with respect to the median profile of the sample. While the two quasars are similarly luminous, in the spectrum of quasar 3.2, a blueshifted absorption trough is visible in the CIV line (Fig. 1), which may indicate a nuclear outflow (Weymann et al. 1991) possibly helping in illuminating the surrounding gas (Costa et al. 2022).

4.4.4. Quasar 4.1 and 4.2

Figure A.3 presents the results for quasar 4.1, while quasar 4.2 is not in the FoV of MUSE due to the large projected separation. Its direction is indicated by a green arrow. As can be seen in Fig. 1, the spectrum of quasar 4.1 shows deep absorption troughs blueshifted of all broad lines and has been therefore classified as broad absorption line quasar in past literature (e.g., Gibson et al. 2009). The extended emission around this object is among the brightest and most extended, spanning almost 100 kpc between the quasar and the 2σ isophote. There is a clear asymmetry in the nebula, also evident in the big offset of flux-weighted nebula centroid and quasar position of ∼20 kpc and the relatively low value of α. A secondary peak of the Lyα surface brightness is visible to the south of the quasar. This substructure could be gas associated with an inflowing companion galaxy, but no continuum source is robustly detected at this position. The i-band forced magnitude in an aperture of 3″ for the substructure is 23.7, while the average background magnitude in the same aperture is 24.3. The increase in surface brightness is also apparent in the surface brightness profile measured within the 2σ isophote, peaking at 300 ckpc. The profile calculated from the full 2D map deviates from the typical quasar pair profile as the central surface brightness is higher, but the emission falls off more steeply than the median profile. Once again, the kinematics are relatively quiescent throughout the full extent of the extended Lyα structure ($\overline{\sigma }=395$ km s⁻¹).

4.4.5. Quasar 5.1 and 5.2

Displayed in Fig. A.4 are the results for quasar 5.1, while quasar 5.2 has not been observed due to the high angular separation. Quasar 5.1 is closely accompanied by two continuum sources first detected in the MUSE data cube and examined in more detail in a companion paper (Herwig et al. 2024). The detected extended emission is found at a significant offset from the quasar, d_{QSO − Neb} = 33.6 kpc, and is not directly connected to it. The velocity shift in this nebula spans a large range compared to the rest of the sample, with an evolution from 500 to −500 km s⁻¹ in the east-west direction. The cool gas traced by the Lyα emission might be expelled from the interstellar medium of the quasar and companion host galaxies due to tidal interactions.

4.4.6. Quasar 6.1 and 6.2

Results for the close quasar pair 6 can be found in Fig. A.5. It is the highest redshift pair in the sample and, possibly as a result of this, shows strong IGM absorption blueward of the Lyα peak. The extended Lyα emission may also be affected by the enhanced IGM absorption compared to the rest of the sample. The quasar pair is associated with a low surface brightness filamentary nebula stretching between the two quasars. The nebula velocity evolves from ∼−300 km s⁻¹ around the fainter quasar with lower redshift to 400 km s⁻¹ close to the brighter quasar with higher redshift, indicating that the nebula is indeed associated with both quasars, whose velocity separation is of 870 km s⁻¹. Therefore, the centroid of the nebula is roughly between the pair, 35 pkpc away from the brighter quasar. The surface brightness profile within the 2σ isophote stays remarkably constant and likewise decreases slowly when calculated in the full surface brightness map. The kinematics are among the most quiescent in this sample, with $\overline{\sigma }=342$ km s⁻¹.

4.4.7. Quasar 7.1 and 7.2

Quasars 7.1 and 7.2 is the only pair at wide (i.e., approximately 500 kpc) separation for which both quasars have been observed (Fig. A.6). The brighter quasar in the pair is associated with a relatively small and faint (L_Neb ∼ 10⁴³ erg s⁻¹) nebula with a significant axis asymmetry and a velocity dispersion around 350 km s⁻¹. Quasar 7.2, on the other hand, is not associated with extended emission exceeding our cut-off area of A_2σ > 2 arcsec² to ensure that the biggest 2σ isophote detected close to the quasar is not caused by spurious emission or noise. For this feature, the velocity shift is almost 0 km s⁻¹ and the velocity dispersion reaches atypically high values above 700 km s⁻¹, supporting the exclusion of this isophote from the morphological analysis.

4.4.8. Quasar 8.1 and 8.2

This object, shown in Figure A.7, has already been published in Arrigoni Battaia et al. (2019b) and is reanalyzed in this work to ensure consistency. While the overall shape of the extended emission stays constant in both works, the method employed here is more conservative (pseudo-NB vs channel-by-channel detection) and consequently does not detect low surface brightness emission features. While the nebulae associated with the two quasars are still connected through a gaseous filament, the emission is much less coherent and smaller in area. However, as the emission appears physically connected, we treat it as one nebula in the analysis.

Quasar pair 8 is embedded in the biggest nebula of the sample with an extent of almost 160 kpc. The velocity shift of the gas evolves from −250 km s⁻¹ around quasar 8.2–350 km s⁻¹ close to quasar 8.1, and it is therefore comparable with the velocity separation of the two quasars (896 km s⁻¹; see discussion in Arrigoni Battaia et al. 2019b). We measure an average velocity dispersion of ∼400 km s⁻¹, with higher dispersions around quasar 8.2 (see Sect. 5.1 for further details). The surface brightness profile of the emission mirrors the median profile of all quasar pairs remarkably well, starting out faint in the center and falling off gently.

5. Discussion

5.1. Comparison to previous quasar pair studies

To our knowledge, the number of previous works targeting extended emission on CGM/IGM scales around physically associated quasar pairs is limited to seven studies (Cai et al. 2018; Arrigoni Battaia et al. 2019b; Lusso et al. 2019; Li et al. 2023), three of which have been serendipitous detections since the faint companions were not known during observation planning (Cantalupo et al. 2014; Hennawi et al. 2015; Arrigoni Battaia et al. 2018). These studies unanimously reported Lyα nebulae with exceptional extents (up to 460 kpc). In our work, however, extended emission often does not cover similarly large areas (> 10⁴ kpc²) or reach such high luminosities L_Neb (> 10⁴⁴ erg s⁻¹).

Li et al. (2023) targeted quasar pairs at z ∼ 2.5 with one bright (g < 18 mag) quasar and included a known ELAN first published in Cai et al. (2018), resulting in two out of five quasar pairs meeting their criteria for ELANe. Compared to that work, nebulae studied here have much smaller extents, even when the integrated luminosity is comparable. This fact is indicative of the differing methods to determine what significant emission means. Here we decided to employ a more conservative approach in order to match previous narrowband study designs, leading to smaller 2σ isophotes. A comparison of these results is therefore not straightforward, but we can firmly state that our data do not unveil nebulosities as bright as those in prototypical ELANe, which have SB_Lyα ∼ 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻² on hundreds of kiloparsecs (e.g., Hennawi et al. 2015; Arrigoni Battaia et al. 2018). Therefore, the sample presented here does not confirm the overabundance of ELANe around the overall population of quasar pairs as we do not identify any ELAN within our sample of 8 pairs. Importantly, the targeted sample extends the study of large-scale emission around quasar pairs to additional systems in which both quasars are faint (Fig. 2). This fact is also one of the reasons why we do not detect bright and extended nebulae around all pairs (see more details in Sect. 5.2). In Fig. 4 we compare our stacked z ∼ 3 quasar pair Lyα profile to the one obtained using the same method as employed in this work by Li et al. (2023) for z ∼ 2.5 quasar pairs after correcting them for cosmological dimming. This exercise shows that overall the SB levels are consistent, contrary to the results found for Lyα emission around single bright quasars, with z ∼ 3 quasars having a factor of ∼3 brighter extended Lyα emission (Cai et al. 2019). In other words, naively, we would have expected quasar pairs at z ∼ 3 to be associated with brighter emission than z ∼ 2 quasar pairs for similar AGN magnitude ranges. Therefore we inferred that the emission in our sample is modest because of the average lower luminosity of the targeted quasars (Sect. 5.2).

Contrarily, the Lyα nebula published by Lusso et al. (2019) in the MUSE Ultra Deep Field (MUDF) around a z ∼ 3 quasar pair (m_r^QSO1 = 17.9 and m_r^QSO2 = 20.5) shows a similar morphology to some objects studied here, in particular to quasar 4.1. Both pairs have a projected separation of about 500 physical kiloparsecs, with the Lyα emission around the brighter quasar extending toward the fainter second quasar and showing an additional smaller nebula along the same direction (Fig. A.3). While the extremely long integration time of 40 hours in Lusso et al. (2019) allows for sensitivity to exceptionally faint emission, the nebula sizes still do not exceed the expected virial radius. The elongation between that quasar pair also appears to fit the relation shown in Fig. 6 for our systems. Quasar pairs like the one targeted in Lusso et al. (2019) and quasar pair 4 are therefore excellent systems for ultra-deep integrations to unveil intergalactic bridges and directly study the IGM.

Lastly, one object studied in this work, pair 8, has already been published in Arrigoni Battaia et al. (2019b). While most of the results are the same, extended emission is detected in that work using a 3D mask, making it more sensitive to faint and narrow emission. Therefore, the area detected above 2σ is bigger, and the nebula is more coherent than in this work. Since the velocity information in that work was computed using only the wavelength ranges within the 3D mask and included the nebula outskirts, the average velocity dispersion is much lower (162 km s⁻¹). We caution that the dispersion calculation is highly dependent on the spectral width considered, and the inclusion of spaxels spanning 1.25–2.5 Å in the 3D mask necessarily leads to very low dispersion measures. This comparison highlights the importance of considering the different methods of extended emission extraction when studying Lyα nebula kinematics.

5.2. Comparison to single quasar samples

A large number of single quasar fields around redshift three have already been targeted with IFUs to study extended Lyα emission, although mostly focusing on more luminous quasars than presented in this work (Fig. 2). Multiple surveys comprising more than one hundred bright z ∼ 3–4 quasars detected extended Lyα emission on scales of, on average, 80–150 kpc in maximum projected size around 100% of the targets (Borisova et al. 2016; Arrigoni Battaia et al. 2019a; Fossati et al. 2021). This is in stark contrast to the quasar pair sample presented here, where two of the 14 observed quasars are not associated with extended emission and the detected nebulae are comparatively small. This difference can be appreciated by comparing the stacked SB profiles presented in Fig. 4. The SB profile of our z ∼ 3 quasar pair sample is 5 times fainter than the SB profile around bright single quasars. This difference is remarkable as some of the targeted quasar pairs could trace ongoing mergers between massive galaxies. In such events several processes (e.g., tides) could boost the gas densities on CGM scales and hence the associated Lyα emission (Sect. 5.4). Further, the differences cannot be explained by the observational setup as the same method is adopted for most of the targets in the aforementioned studies (i.e., an exposure time of at least 1 hour in MUSE WFM) but it could possibly arise due to target selection as our quasar pair sample has not been selected from SDSS based on having bright quasars in the pair.

Indeed, a single quasar survey selecting 12 similarly low luminosity quasars at z ∼ 3.1 found extended emission in each targeted field with a minimum extent of 60 pkpc, but did show a connection between quasar Lyα peak luminosity and nebula luminosity (Mackenzie et al. 2021). We note that Mackenzie et al. (2021) adopted an optimal extraction method for constructing a 3D region for the study of the Lyα emission. This technique, applied channel by channel, is more sensitive to the fainter portions of the nebulae with respect to our method. For this reason we tested our analysis by applying it to the two extreme nebula cases in that sample (see Appendix C). We find that our analysis is able to reproduce the principal shape and luminosity of the brightest nebula reported in that paper, while our measurement of the extent deviates by 30%. The faint nebula is only detected in the peak of the nebula emission and is thus three times smaller in maximum extent than reported in Mackenzie et al. (2021). Therefore, the difference in maximum extent of 15% between the faint single-quasar sample and our sample of quasar pairs is likely in part due to the difference in methodology.

The stacked SB profile from that work has roughly similar SB values as the stacked profile for the quasar pairs (Fig. 4), but declines faster at large radii. It is therefore likely that the gaseous environments around quasar pairs at z ∼ 3 are not the same as around single quasars at similar redshifts, indicating the presence of more cool gas at large projected distances.

This is further confirmed by comparing our stacked SB profile with the stacked profile of 59 similarly faint (−27 < M_i(z = 2) <  − 24) z ∼ 3.1 quasars from the QSO MUSEUM III survey (González Lobos et al., in prep.). The surface brightness is remarkably similar at intermediate distances to the central quasar, but is brighter at close separation and falls off quicker, in agreement with the profile in Mackenzie et al. (2021).

We quantify the differences between SB profiles around single quasars and quasar pairs by fitting them with power laws. The SB profile for single z ∼ 3 quasars can be described by a power law with an index reported to be between −1.8 (Borisova et al. 2016) and −1.96 (Arrigoni Battaia et al. 2019a) for bright quasars, and −2.41 ± 0.19 for the faint sample (González Lobos et al. in prep.). In contrast, we find a slope of −1.57 ± 0.17 for the quasar pair sample. Depending on the Lyα powering mechanism, the shallower slope of the profile could be evidence of (i) presence of denser and/or a bigger reservoir of ionized gas at larger radii in the CGM of pairs in comparison to single quasars (if Lyα is dominated by recombination radiation; Hennawi & Prochaska 2013) and/or (ii) more neutral hydrogen on large scales in the CGM of pairs (if Lyα is dominated by resonant scattered photons). Intriguingly, Lusso et al. (2018) found evidence for moderate excess of optically thick absorbers (with column densities of HI log N_HI ≥ 17.2) in closely projected z ∼ 2 quasar pairs compared to single quasars, arguing that they are likely to trace mostly structures located in denser (partially neutral) regions within the CGM or IGM where both quasars reside.

Further evidence of the presence of large-scale structures comes from the morphologies of the discovered nebulae around quasar pairs. In section 4.3 we have shown a comparison of their morphologies (described by their axis ratios) with respect to nebulae around single quasars (Fig. 5). Morphologies of nebulae around single quasars are mostly circular (Borisova et al. 2016; Arrigoni Battaia et al. 2019a), with even higher symmetry for faint quasars (Mackenzie et al. 2021). ELANe, known to host multiple AGNs, are more asymmetric (Arrigoni Battaia et al. 2023b). Similarly, quasar pair nebulae have a median axis ratio lower than nebulae around single quasars (Fig. 5). This is akin to redshift ∼2 single quasar nebulae with a median value of α = 0.54 (Cai et al. 2019) hypothesized to originate in galaxy mergers, which peak at z = 2. While this interpretation can be supported in the context of quasar pairs, the spread of α in our small sample is substantial and α is not directly correlated with the projected distance. Asymmetries can therefore not only be driven by galaxy mergers, but also by intergalactic structures.

5.3. Direction of nebulae with respect to the quasar pairs

Assuming each quasar inhabits a node of the cosmic web, the connecting line between two members of a pair can be used as a proxy for the direction of an intergalactic filament. Following this argument, the value of the angle ϕ, introduced in Sect. 4.3, indicates the level of alignment between the discovered Lyα nebulae and the cosmic web surrounding their host galaxies. A small angle ϕ, together with a large elongation d_QSO, max (see Sect. 4.3 for definition) should correspond to the best candidates for intergalactic bridges as tested in Fig. 6.

This expectation is based on the assumption that a nebula semi-major axis (used to compute ϕ) is a good proxy for the direction of the underlying cool gas distribution. Cosmological simulations post-processed with a radiative transfer tool have shown that this is true if the quasar host galaxy is not edge-on (Costa et al. 2022). The more edge-on the host galaxy the more lopsided the extended Lyα emission. The fact that the majority of quasar nebulae are symmetric may therefore indicate that their host galaxies are far from being edge-on and are sampling the surrounding hydrogen gas distribution well (Arrigoni Battaia et al. 2023b).

These simple predictions are further complicated by the fact that the morphology of the Lyα nebulae (and hence d_QSO, max and the angle ϕ) can be influenced by a multitude of factors. Indeed it has been shown that the following factors all play a role in shaping the extents and luminosities of quasar nebulae: dust in the host galaxy or the CGM (Smith et al. 2022; González Lobos et al. 2023) as dust can efficiently absorb Lyα photons; the molecular gas mass of the host galaxy, being a proxy for denser and more dust rich galaxies that may prevent the escape of ionizing and Lyα photons (Muñoz-Elgueta et al. 2022); the Lyα luminosity and/or magnitude of the powering sources, with the trend of brighter quasars having brighter and often more extended nebulae (Mackenzie et al. 2021); the availability of cool, 10⁴ K gas that has the potential to be Lyα bright; and geometry effects (e.g., orientation of the quasar’s ionizing cones with respect to the gas distribution; Obreja et al. 2024; orientation of cosmic web filaments with respect to the observing angle, which could lead to absorption of Lyα). In addition to those, quasar pair nebulae might not trace exclusively inflowing CGM gas (see Sect. 5.4), especially in close pairs.

However, an anti-correlation between nebula misalignment and extent can be seen in our quasar pair sample irrespective of quasar magnitude and projected pair distance (Fig. 6), indicative of a direct link between filaments of the cosmic web and the distribution of cool gas in the CGM. This trend suggests that the extended nebulae discovered in our study with small angle ϕ and large d_QSO, max (e.g., quasar pair 4) probe gas falling into the quasar dark-matter halo after accretion from the cosmic web. Contrarily, less extended nebulae are perfect laboratories to test all the aforementioned processes that could strongly affect the morphology of Lyα nebulae. This means that any alignment found between large-scale galaxy distributions around quasars and their Lyα nebulae (Arrigoni Battaia et al. 2023b) should break for extremely small nebulae as the smallest discovered in this study.

In the next section we discuss more about the origin of the cool gas as probed in all the diversity of quasar pair nebulae in this study.

5.4. Origin of the cool gas probed by Lyα around quasar pairs

The gas traced by extended Lyα emission around z ∼ 3 quasars is usually assumed to be the cool (∼10⁴ K) phase of the CGM, supported by typical Lyα halo radii < 100 kpc for the bulk of the emission, just short of the expected virial radius for quasar halos (e.g., Arrigoni Battaia et al. 2019a). This is further confirmed by the levels of the Lyα SB profiles around quasars as a function of redshift. The absence of an evolution from z ∼ 6 down to z ∼ 3 (Farina et al. 2019; Fossati et al. 2021) and the subsequent decrease to z ∼ 2 (Arrigoni Battaia et al. 2019a; Cai et al. 2019) seems in agreement with the theoretical expectation (e.g., Dekel & Birnboim 2006) that high-z (z ≳ 3) massive halos are able to accrete gas in a cool phase keeping the Lyα emission roughly constant, while at lower redshift, the gas mass in the cool phase (and hence the Lyα emission) is reduced by shock heating. At moderate depths (∼10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²), a few rare exceptions show very extended emission that can exceed the projected virial radius of the host halo and therefore reach into the IGM (e.g., Cantalupo et al. 2014; Arrigoni Battaia et al. 2019a). However, most of these outliers have active companions.

Quasar pairs seem to follow the same trends as single quasars, but with, on average, enhanced Lyα emission at larger distances (Sect. 5.2). The Lyα emission in quasar pairs can therefore better trace IGM structures given the closer presence of an additional massive halo. Indeed, cosmological simulations report the presence of denser and thicker filaments in denser environments (Cautun et al. 2014), similar to the locations quasar pairs are expected to populate in the cosmic web (e.g., Boris et al. 2007; Onoue et al. 2018). Given the presence of the second quasar, they are at least in a group environment. For this reason, processes favored by dense environments could affect the cool gas distribution more than around single quasars. In particular, galaxy interactions may become especially important in pairs at close separation. Both tidal stripping and ram pressure stripping could significantly enhance the density of the CGM in a preferential direction aligned with an interacting pair (Salem et al. 2015; Samuel et al. 2023) and may therefore boost the observed surface brightness out to large radii. Evidence for gas stripping around high-z quasars is indeed mounting (e.g., Chen et al. 2021; Vayner et al. 2023). The effects of these environmental processes can, in part, explain the alignment of extended nebulae, specifically in close pairs (Fig. 6), and the comparatively flat surface brightness profile observed in the median stack of quasar pairs (Fig. 4) and especially apparent in the profiles of some individual quasar pairs, for example pair 6 (Fig. A.5). However, the properties of the extended emission studied here (Table 3) do not correlate with the angular separation between the quasar pairs, as would be expected if the galaxy interaction alters the CGM in an observable way. Therefore, the morphology of extended Lyα emission alone does not provide evidence for the aforementioned environmental effects. Summarizing, the diversity of extended Lyα emission around quasar pairs probes the combined effects of dense environments and close massive halos, possibly tracing the surrounding large scale structure.

6. Summary and conclusions

Motivated by previously discovered ELANe around quasar pairs and signs of filamentary gas connecting two quasars in multiple systems, we conducted the largest study of extended Lyα emission around quasar pairs to date. Our sample consists of 14 quasars in eight pairs observed with VLT/MUSE with a snapshot strategy (45 minutes/source). The targets span redshifts of z ∼ 3–4 and i-band magnitudes between 18 and ∼22.8, dimmer than most single-quasar samples. Our key findings are as follows:

• Twelve out of the 14 observed quasars are associated with extended Lyα emission spanning 40–160 kpc, and three of the five close pairs are embedded in the same Lyα nebula (Table 3, Figs. 3, A.1–A.7). In contrast to previously discovered quasar pairs, none of these nebulosities meet the criteria for ELANe in either L_Neb or extent, lowering the fraction of ELANe around known quasar pairs significantly (Sect. 5.1). On average, the nebulae in our sample are 15% smaller in maximum projected extent than nebulae around similar single-quasar samples, but they have comparably quiescent kinematics (Section 5.2).

• Larger nebulae are preferentially aligned with the expected direction of the large-scale structure between the associated quasar pair, while smaller nebulae tend to be misaligned (Fig. 6). This fact is indicative of a direct connection between filaments and the cool CGM (Sect. 5.3). Wide pairs highly aligned with their filaments, such as quasar pair 4, are excellent targets for ultra-deep follow-up observations to detect these faint IGM structures.

• The SB profile of quasar pairs is comparable in brightness to single-quasar profiles at similar redshift but is fainter close to the quasar and falls off more gradually with a power law slope of −1.57 ± 0.17 (Fig. 4). Due to galaxy interactions important in close quasar pairs, the higher surface brightness at larger radii can in part be explained by stripped gas increasing the densities in the CGM. The presence of a second large halo can also add to the Lyα emission through increased contribution from IGM structures (Sect. 5.4).

• Extended emission around this sample of quasar pairs is morphologically diverse, but on average, it is more asymmetric than single-quasar nebulae with large offsets between the quasar position and the flux-weighted nebula centroid (Table 3, Fig. 5). This can be explained by the presence of more substructures in the overdense regions likely populated by quasar pairs and by close pairs sometimes sharing the same nebula.

• We did not find a trend in brightness or morphology of the extended Lyα emission with the projected distance between the associated quasar pair, as would be expected if the host galaxies’ interaction influences the CGM in a way that is observable through the extended Lyα emission (Sect. 5.4). While the sample studied in this work only encompasses close (∼50–100 pkpc) and wide (∼450–500 pkpc) pairs, extending the sample of observed quasar pairs to intermediate projected distances has the potential to reveal the full picture.

This study highlights once again the promise of Lyα nebulae around quasars tracing the large-scale structure in dense environments. As quasar pair nebulae show a higher surface brightness level at larger radii than individual quasars, and they provide the most accessible laboratory to directly study the IGM and its connection to the associated galaxies. Indeed, our sample contains excellent targets for future deep and ultra-deep observations to study the large-scale gaseous structures of the IGM (e.g., Tornotti et al. 2024). Additional snapshot observations of quasar pairs are needed to increase the number of such suitable targets and therefore cover the diversity of IGM filament structures and active galaxy interactions at different separation stages.

Acknowledgments

We thank the referee for their useful comments that helped us improve the manuscript. Based on observations collected at the European Southern Observatory under ESO programme 0100.A-0045(A). AWSM acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through grant reference number RGPIN-2021-03046. E.P.F. is supported by the international Gemini Observatory, a program of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the United States of America. All figures have been made with matplotlib (Hunter 2007) and APLpy (Robitaille & Bressert 2012; Robitaille 2019). We have also used the Python libraries numpy (van der Walt et al. 2011), scipy (Virtanen et al. 2020), spectral-cube (Ginsburg et al. 2015), and astropy (Astropy Collaboration 2013, 2018, 2022).

References

Appendix A: Additional material

In this appendix, we summarize important information about the MUSE observations and the quasar properties in Table A.1 and show in Figs. A.1 to A.7 Lyα maps and profiles for the remaining seven quasar pairs not shown in Fig. 3.

Fig. A.1. Same as Fig. 3 but for quasar pair 2 (ID 2.1 and ID 2.2).

Fig. A.2. Same as Figure 3 but for quasar pair 3 (ID 3.1 and ID 3.2).

Fig. A.3. Same as Fig. 3 but for quasar pair 4 (ID 4.1). Quasar 4.2 is not within the MUSE FoV, but the direction toward it is indicated with a green arrow in the Lyα surface brightness map. The white arrow indicates the line corresponding to dQSO, max in this object.

Fig. A.4. Same as Fig. 3 but for quasar pair 5 (ID 5.1). Quasar 5.2 is not within the MUSE FoV, but the direction toward it is indicated with a green arrow in the Lyα surface brightness map.

Fig. A.5. Same as Fig. 3 but for quasar pair 6 (ID 6.1 and ID 6.2).

Fig. A.6. Same as Fig. 3 but for quasar pair 7. As this is a wide pair, maps are presented for each quasar individually. Top row: Lyα surface brightness map, line velocity and velocity dispersion for quasar 7.1. Middle row: Lyα surface brightness map, line velocity and velocity dispersion for quasar 7.2.

Fig. A.7. Same as Fig. 3 but for quasar pair 8 (ID 8.1 and ID 8.2).

Appendix B: Example of a field before quasar PSF subtraction

As explained in Sect. 3, we subtracted the unresolved emission from each quasar to search for the extended Lyα emission around the pairs. For completeness, we show in Fig. B.1 one example of one field (quasar pair 1) before subtraction of the quasar PSF. As can be seen from the colorbar, the quasar PSF is much brighter than the unveiled extended emission shown on the same image cut-out in Fig. 3.

Fig. B.1. Surface brightness in the PNB of quasar pair 1 before PSF subtraction and spatial smoothing. Shown is the same image cutout as in Fig. 3 but with a different color bar.

This becomes even more apparent when plotting the surface brightness profile of the data before subtraction, the constructed PSF and the data after subtraction in linearly spaced annuli as shown in Fig. B.2 for PNBs of quasar pair 1. At the peak of the PSF, the unveiled extended emission is roughly 1.8 orders of magnitude dimmer than the quasar emission and it only starts dominating the emission in the PNB at a distance from the quasar of roughly 200 ckpc. We emphasize that during the analysis, we subtracted the empirical PSF channel-by-channel, and Fig. B.2 therefore only functions as a simplified visualization of the subtraction.

Fig. B.2. Surface brightness profile of quasar pair 1 determined in linearly spaced annuli in the PNB before PSF subtraction (black), after PSF subtraction (blue), and the empirical PSF constructed (gray dashed line).

Appendix C: Reanalyzed data of faint single quasar nebulae

Fig. C.1. Lyα SB maps of the faintest (left) and brightest (right) nebula first reported in Mackenzie et al. (2021) (ID4 and ID6 respectively) and reanalyzed using the method presented in this work. The maps are plotted as described in Fig. 3 within a 28″ × 28″ window.

In this work, we employ a conservative approach of using PNBs to detect extended Lyα emission. Previously published samples of nebulae have often been optimally extracted within a 3D mask containing all pixels above a signal-to-noise ratio of 2 that are associated with the nebula. That method uses a less conservative detection criteria, and it could therefore result in more extended detections. To facilitate comparison, we reanalyze the two most extreme nebulae published in Mackenzie et al. (2021) using 30 Å PNBs: ID4 is reported to host the faintest and smallest nebula (Fig. C.1, left); ID6 is associated with the brightest and most extended Lyα emission (Fig. C.1, right). We find that our method is able to detect extended Lyα emission as shown in Mackenzie et al. (2021), with the only exception of thin structures at low S/N and at the edges of the nebulae in their optimally extracted maps. For this reason, we find that our measurements for the extents (20 kpc for the small nebula and 134 kpc for the large nebula) are smaller than those reported in that work (60 kpc and 190 kpc). These differences are clearly due to the different extraction criteria for the nebulae but do not affect our results.

All Tables

All Figures

Fig. 1. Eight targeted quasar pairs. Each panel shows the spectra for the two AGNs extracted either from the final MUSE data cubes in an aperture of 3″or available from SDSS. Vertical dotted lines indicate the location of key emission lines at the SDSS redshift of the brighter quasar. The spectrum of quasar 1.1 is scaled down by a factor of 10. The mean noise per channel in the MUSE data at the respective wavelength of Lyα is $0.08·{10}^{-17}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{\AA }}^{-1}$.

In the text

Fig. 2. Quasar magnitudes and redshift covered by different studies of Lyα nebulae. Single quasars are shown as dots (Borisova et al. 2016; Arrigoni Battaia et al. 2019a; Cai et al. 2019; Mackenzie et al. 2021; Fossati et al. 2021; González Lobos et al., in prep.), while physically associated quasar pairs are connected with a line (our work as purple stars; Cai et al. 2018; Lusso et al. 2019; Li et al. 2023).

In the text

Fig. 3. Lyα maps and profiles for quasar pair 1 (ID 1.1 and ID 1.2). All maps have been smoothed by a 2D Gaussian kernel with a standard deviation of 2 pixels or 0.4″ and were obtained from the PSF-subtracted optimized PNBs. Shown is a 26″ × 26″ (or about 200 × 200 kpc) cut-out centered between quasar 1.1 and quasar 1.2. Top left: Lyα surface brightness map. The position of the brighter quasar is indicated with a light blue, diagonally dashed circle with a diameter of 1 arcsec, the dimmer quasar is indicated by the darker blue, horizontally dashed circle with the same diameter. Contour lines indicate a detection of 2σ (dark gray), 4σ (light gray) and 10σ (white). The surface brightness maps of all pairs are shown with the same colorbar range to ease comparison. Top middle: Velocity offset with respect to the PNB center of the Lyα line within the 2σ isophote, calculated within the PNB wavelength range. Quasar positions are marked with circles color-coded by their systemic redshift on the same colorbar scale as the map. Top right: Velocity dispersion of the Lyα line within the 2σ isophote. Bottom left: Spectrum of the brighter (QSO1) and dimmer (QSO2) quasar integrated within an aperture of 3″ and spectrum of the nebula integrated within the 2σ isophote indicated in the surface brightness map. The spectra are centered on the Lyα peak of the nebula chosen as the PNB center (black dashed line). The slices used for the PNB are marked in gray. All spectra are smoothed with a Gaussian kernel with a standard deviation of 1 pixel or 1.25 Å. The spectrum of quasar 1.1 is scaled down by a factor of 10. Bottom right: Cosmologically corrected surface brightness profiles calculated in logarithmically spaced annuli around the brighter quasar. Quasar positions are indicated by the dots, color coded by Mi(z = 2). The profile represented by red dots is calculated in the full PNB after masking continuum sources and subtracting the background, the errorbars in y-direction show the 1σ error. For comparison, the median stacked profile of all pairs is plotted as black dashed line. The green curve is calculated within the 2σ isophote indicated in the surface brightness map, with shaded regions showing the 1σ error calculated in the same annuli in the variance cube.

In the text

Fig. 4. Cosmology corrected, circularly averaged Lyα surface brightness profiles of different quasar samples. The obtained average Lyα SB profile for z ∼ 3 quasar pairs is compared with (i) samples of individual z ∼ 3 (bright quasars, Arrigoni Battaia et al. 2019a; Fossati et al. 2021; faint quasars, Mackenzie et al. 2021; González Lobos et al., in prep.) and (ii) with z ∼ 2 quasar pairs (Li et al. 2023).

In the text

Fig. 5. Lyα nebulae axis ratio versus their luminosity. The values for the nebulae discovered around the studied quasar pairs (red stars, accompanied by the quasar ID) are compared with those for nebulae found around bright (cyan dots and cyan dash-dotted line; Arrigoni Battaia et al. 2019a) and faint (dotted gray line; Mackenzie et al. 2021) z ∼ 3 quasars. The cyan star indicates the value for the Fabulous ELAN in Arrigoni Battaia et al. (2019a), which is hosted by an AGN triplet (of which two are type-I quasars as in this work). The cyan shaded region indicates the 25th and 75th percentiles of the distribution of bright quasars.

In the text

Fig. 6. Angle between the Lyα nebula semi-major axis and the line connecting the two quasars, ϕ versus the maximum distance between the bright quasar in each pair and its associated Lyα nebula 2σ isophote, dQSO, max. The data points are color coded following Mi(z = 2) of the associated quasar, while their sizes indicate the projected distance between the quasars in each pair. The quasar IDs are displayed next to the respective data points.

In the text

	Fig. A.1. Same as Fig. 3 but for quasar pair 2 (ID 2.1 and ID 2.2).
In the text

	Fig. A.2. Same as Figure 3 but for quasar pair 3 (ID 3.1 and ID 3.2).
In the text

	Fig. A.3. Same as Fig. 3 but for quasar pair 4 (ID 4.1). Quasar 4.2 is not within the MUSE FoV, but the direction toward it is indicated with a green arrow in the Lyα surface brightness map. The white arrow indicates the line corresponding to dQSO, max in this object.
In the text

	Fig. A.4. Same as Fig. 3 but for quasar pair 5 (ID 5.1). Quasar 5.2 is not within the MUSE FoV, but the direction toward it is indicated with a green arrow in the Lyα surface brightness map.
In the text

	Fig. A.5. Same as Fig. 3 but for quasar pair 6 (ID 6.1 and ID 6.2).
In the text

	Fig. A.6. Same as Fig. 3 but for quasar pair 7. As this is a wide pair, maps are presented for each quasar individually. Top row: Lyα surface brightness map, line velocity and velocity dispersion for quasar 7.1. Middle row: Lyα surface brightness map, line velocity and velocity dispersion for quasar 7.2.
In the text

	Fig. A.7. Same as Fig. 3 but for quasar pair 8 (ID 8.1 and ID 8.2).
In the text

	Fig. B.1. Surface brightness in the PNB of quasar pair 1 before PSF subtraction and spatial smoothing. Shown is the same image cutout as in Fig. 3 but with a different color bar.
In the text

	Fig. B.2. Surface brightness profile of quasar pair 1 determined in linearly spaced annuli in the PNB before PSF subtraction (black), after PSF subtraction (blue), and the empirical PSF constructed (gray dashed line).
In the text

	Fig. C.1. Lyα SB maps of the faintest (left) and brightest (right) nebula first reported in Mackenzie et al. (2021) (ID4 and ID6 respectively) and reanalyzed using the method presented in this work. The maps are plotted as described in Fig. 3 within a 28″ × 28″ window.
In the text