© F. Arrigoni Battaia et al. 2019

1. Introduction

The current paradigm of structure formation predicts the presence of gaseous filaments connecting galaxies (e.g., White et al. 1987; Bond et al. 1996), ultimately forming an intricate web known as the intergalactic medium (IGM; Meiksin 2009). Given the expected low densities for such gas (n_H  ≲  0.01 cm⁻²) and the budget of ionizing photons in the ultraviolet background (UVB; e.g., Haardt & Madau 2012), the direct observation of the IGM is anticipated to be very challenging. For instance, the surface brightness in Lyα emission is predicted to be SB_Lyα  ∼  10⁻¹⁹−10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² (Gould & Weinberg 1996; Bertone & Schaye 2012; Witstok et al. 2019). Indeed, a direct detection of the IGM thus far appears to be elusive even with current facilities (e.g., Gallego et al. 2018; Wisotzki et al. 2018), such as, the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) and the Keck Cosmic Web Imager (KCWI; Morrissey et al. 2012).

However, it was noticed early on that quasars could act as flashlights, possibly photoionizing the surrounding medium out to large distances. The ionized gas would then recombine, emitting numerous Hydrogen Lyman-α (Lyα) photons as main product (e.g., Rees 1988; Haiman & Rees 2001). This boosted glow, SB_Lyα >  10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², may possibly then be within reach of state-of-the-art instruments (Cantalupo et al. 2005; Kollmeier et al. 2010).

Following this idea, several works aimed for Lyα emission from halos (nowadays known as the circumgalactic-medium, CGM; Tumlinson et al. 2017) out to intergalactic scales around individual high-z quasars in order to constrain the physical properties of the diffuse gas phases (e.g., Hu & Cowie 1987; Heckman et al. 1991; Møller et al. 2000; Weidinger et al. 2004, 2005; Christensen et al. 2006; Hennawi et al. 2009, 2015; Cantalupo et al. 2014; Martin et al. 2014; Arrigoni Battaia et al. 2016; Farina et al. 2017). At z  ∼  3, observations can now easily (∼1 h on source) uncover the emission within 50 projected kpc, and reach an average maximum distance of ∼80 projected kpc from the targeted quasar (Borisova et al. 2016; Arrigoni Battaia et al. 2019). This Lyα emission usually shows SB_Lyα  ∼  4 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² with relatively quiescent line widths σ_Lyα <  400 km s⁻¹ (Arrigoni Battaia et al. 2019), which are intriguingly similar to the velocity dispersion expected for halos hosting quasars at these redshifts (σ  =  250 km s⁻¹; M_DM  ∼  10^12.5 M_⊙; e.g., White et al. 2012). The uncertainties in the determination of the quasar systemic redshift along with the haze possibly introduced by the Lyα radiative transfer still hamper, in most cases, a secure interpretation of the gas kinematics and/or configuration of the system as traced by the extended Lyα emission (Arrigoni Battaia et al. 2019). Notwithstanding these open issues, the Lyα nebulae are usually interpreted to be physically associated with the targeted quasar and trace either the gravitational motions due to structure assembly (Weidinger et al. 2004, 2005; Arrigoni Battaia et al. 2018) or the violent feedback of the central engine (Cai et al. 2017). Alternative interpretations explain that the extended Lyα emission is not strictly associated with the targeted quasars, but it is due to structures along our line-of-sight to the quasar, like portions of the CGM of massive halos in the Hubble flow aligned along our line-of-sight (Cantalupo et al. 2019) or proto-galactic disks (Martin et al. 2019) illuminated by the quasar.

Thanks to the aforementioned effort in the detection of the CGM around high-z quasars, it starts to become evident that even around individual quasars it is extremely hard to detect diffuse emission at intergalactic distances (> 100 kpc) unless additional companions (mostly active) are present in close proximity (Hennawi et al. 2015; Arrigoni Battaia et al. 2019, 2018), or much more sensitive observations are conducted. Dense environments seem to supply additional cool dense gas necessary for the detection of Lyα signal on very large scales (Hennawi et al. 2015; Cai et al. 2017; Arrigoni Battaia et al. 2018). Further, the unification model for active galactic nuclei (AGN; e.g., Antonucci 1993) and evidences of anisotropic ionizing emission from high-redshift quasars (e.g., Hennawi & Prochaska 2007) hint to the existence of shadowed regions around individual quasars. The presence of multiple quasars within the same structure thus increases the probability of large-scale gas to be illuminated by hard ionizing photons. For these reasons, scientific teams have started to change their approach in unveiling IGM emission, passing from the targeting of individual quasars to short (Cai et al. 2018) or extremely long integrations (> 40 h; Lusso et al. 2019) of multiple high-redshift quasars, or overdensities hosting quasars (Cai et al. 2017).

Here, we report on our effort within this framework. In particular, in 2015 we designed a survey of z  ∼  3 physically associated quasar pairs using the MUSE instrument on the Very Large Telescope (VLT) of the European Southern Observatory (ESO). We now have the first data of these observations, and here we present the results of the study of the first target. Our work is structured as follows. In Sect. 2 we explain how we selected the quasar pairs in our survey. Section 3 presents the observations and data reduction for the quasar pair here studied. We highlight our results for the extended Lyα emission and for the detected absorptions in Sect. 4. Section 5 discusses the possible scenarios for the powering of the extended Lyα emission, while Sect. 6 presents the results of the modeling of the absorbers. Finally, we summarize our findings in Sect. 7.

We adopt the cosmological parameters H₀  =  70 km s⁻¹ Mpc⁻¹, Ω_M  =  0.3, and Ω_Λ  =  0.7, and therefore 1″ corresponds to about 7.7 kpc at z  =  3.020 (z_QSO2; details in Sect. 2). All magnitudes are in the AB system (Oke 1974), and all distances are proper.

2. Selection of the quasar pairs

The quasar pairs to be observed in our program have been selected from the twelfth data release of the Sloan Digital Sky Survey (SDSS) quasar catalog (Pâris et al. 2017) using the following criteria. First, a quasar pair has to be at the lowest redshift for which the Lyα emission is detectable with MUSE, that is, 3.0 ≲ z <  3.9, where sky lines are not dominant; secondly, the two quasars should have a difference in redshift of Δz ≤ 0.03 (corresponding to ≤2000 km s⁻¹). This small difference in redshift should ensure that the two quasars are physically associated (e.g., Hennawi et al. 2006a,b); thirdly, the two quasars should have a projected separation ≤0.5 arcmin, so that both quasars sit within the MUSE field-of-view; lastly, a quasar pair should be well visible from VLT ESO, that is, Dec < 27°.

Importantly, in our selection we did not impose any constraint on the current luminosity of the quasars in the pair¹. Our effort is thus complementary to the approach of Cai et al. (2018), who selected pairs with at least one bright quasar (g <  19) visible from the Palomar and Keck sites. The aforementioned criteria resulted in the selection of a total of 17 quasar pairs visible during the ESO semester P100. We however obtained data only on seven of these targets due to weather conditions.

In this work we focus on the quasar pair SDSS J113502.03−022110.9 – SDSS J113502.50−022120.1 (henceforth QSO1–QSO2), separated by 11.6″ (or 89 kpc) and whose properties are summarized in Table 1. In particular, we double checked the redshift estimate of the SDSS catalog by using the known relation for the blueshift of the [C III] line emission (Shen et al. 2016), and obtained consistent redshifts within the uncertainties². For completeness, we list in the table both redshifts, but we use the SDSS redshifts in the reminder of this work. The current redshift estimates place the two quasars at Δv = 896 ± 316 km s⁻¹, which corresponds to a distance of 2.9 ± 0.9 physical Mpc if all the velocity difference is due to the Hubble flow. However, if we look at their spectra (e.g., Fig. 1), the observed Lyα emission peaks are only separated by Δv = 598 ± 98 km s⁻¹ (or 1.9 ± 0.3 physical Mpc)³. These two quasars are ≈3.8 mag fainter than the average M₁₄₅₀ = −27.12 of the QSO MUSEUM sample of Arrigoni Battaia et al. (2019), and sit in a portion of sky with low galactic extinction A_V = 0.08 mag (Schlegel et al. 1998)⁴.

Fig. 1. One-dimensional spectra (black) for the two quasars of the pair, QSO1 (top) and QSO2 (bottom), as extracted from the final MUSE data, using a circular aperture with radius 2′′. The red spectra indicate the error vectors. The vertical dashed blue (magenta) lines indicate the position of important line emissions at the systemic redshift of QSO1 (QSO2). For both objects, we show in the inset plots a zoomed version of the spectrum at the location of the Lyα and C IV lines to highlight the presence of interesting absorptions. The vertical dashed gray and green lines indicate the location of an H Lyα absorptions present within both QSO1 and QSO2 spectra, and C IV absorption along the QSO2 sight-line. The fit to these lines is shown in Sect. 4.2 and Figs. 6 and 7. Residuals due to frequent sky lines are evident at wavelengths >7000 Å.

3. Observations and data reduction

The quasar pair QSO1–QSO2 was observed during UT 19 of February 2018 with clear sky conditions for the program 0100.A-0045(A) with the MUSE instrument on the VLT 8.2 m telescope YEPUN (UT4). The observations consisted of three exposures of 880 s each, rotated with respect to each other by 90°, and with a dither of a few arcseconds between them. The data have been acquired with the nominal spectral range, and thus cover the wavelengths 4750.2−9350.2 Å.

The data were reduced using the MUSE pipeline recipes v2.2 (Weilbacher et al. 2014). In particular, each of the individual exposures have been bias-subtracted, flat-fielded, twilight and illumination corrected, sky-subtracted, and wavelength calibrated using the calibration data taken closest in time to the science frames. The flux calibration of each exposure has been obtained using a spectrophotometric standard star observed during the same night of the science observing block. The individual exposures were then combined into a single data cube. While we apply an initial sky subtraction using the MUSE pipeline, skyline residuals are further suppressed using the software ZAP⁵ (Soto et al. 2016)⁶. The seeing of the final combined data is measured from the star 2MASS J11350307−0220597 (see Appendix A), resulting in a Moffat function with β = 2.5 and FWHM = 1.66″. The coadded spectrum of QSO1 and QSO2 as extracted from the final MUSE datacube are shown in Fig. 1. Furthermore, we present the white-light image of the observations field of view, obtained by collapsing the final MUSE datacube in Fig. 2.

Fig. 2. White-light image of the observed 57″ × 57″ field of view. We indicate the position of QSO1, QSO2, the star 2MASS J11350307−0220597 used to compute the point spread function of our data (Appendix A), and an interloper galaxy “G” (tentatively at z = 0.457 ± 0.001; Appendix B). Additionally, we indicate the 2σ isophote for the extended Lyα emission discovered around the quasar pair (Fig. 3).

The MUSE pipeline produces a variance datacube which is known to underestimate the true noise because it neglects the correlated noise introduced during the resampling of the datacubes (e.g., Borisova et al. 2016). To correct for this effect, we rescaled the variance cube layer by layer so that the average of each layer in the variance cube matches the average variance computed from each science layer after masking objects.

The final MUSE datacube has a 2σ surface brightness limit of SB_Lyα = 7 × 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² (in 1 arcsec² aperture) in a single channel (1.25 Å) at ≈4872 Å (Lyα at the redshift of QSO2). Given the stability of MUSE, further smoothing can allow us to push this sensitivity to lower levels (Sect. 3.2).

3.1. Subtraction of the quasars point spread function

Quasars easily outshine the radiation produced by the surrounding gas distribution, and their emission is smeared out to larger scales due to the seeing. For these reasons, the study of large-scale gas around quasars requires the subtraction of the unresolved quasar emission, as characterized by the point-spread-function (PSF) of the observations. This problem has been empirically tackled in the literature by subtracting a wavelength-dependent PSF constructed from the data themselves in several ways (e.g., Møller 2000; Christensen et al. 2006; Husemann et al. 2014; Borisova et al. 2016). Given the presence of a bright star within the field-of-view of our observations, we were able to reconstruct the wavelength dependent PSF layer by layer at high signal to noise (S/N), as described in detail in the Appendix A. The reconstructed layer-by-layer PSF was then subtracted at each quasar position out to a 5″ radius after matching the quasar emission within a 1 arcsec² circle. Before proceeding with the extraction of the Lyα signal, we removed all the continuum-detected sources from the datacube using the median-filtering routine contsubfits in ZAP (Soto et al. 2016). We masked the location of very bright or extended continuum objects, like the star 2MASS J11350307−0220597 (Appendix A) and an interloper galaxy “G” tentatively at z = 0.457 ± 0.001 (more details in Appendix B) to avoid contamination from residuals.

3.2. Extraction of the Lyα emission

The Lyα signal was extracted from the final PSF and continuum subtracted datacube using custom routines written in the Python Programming Language⁷. First, we produced a sub-cube of the wavelength range where the extended Lyα line is expected, allowing for wide shifts of the line, that is, ±7500 km s⁻¹ from the two quasar systemic redshifts. This sub-cube covers 4750.2−5010 Å. In the next step we smoothed the sub-cube with a Gaussian kernel of FWHM = 1.66″ (i.e., similar to the seeing).

We further marked all the regions in all the layers of the smoothed sub-cube above a S/N = 2⁸, obtaining a segmentation map for each layer. Using these segmentation maps, we found the largest connected area with detection above S/N = 2 to be of 2077 spaxels (or 83 arcsec²), in the layer of 4872.7 Å (basically at the systemic redshift of QSO2). We walk through the cube starting from this layer, first in the direction of increasing wavelength, and then in the direction of decreasing λ, to obtain a three-dimensional (3D) mask describing the Lyα emission within the datacube. As we move from layer to layer in either only increasing or only decreasing λ, we recursively attach to the detection area (defined from the previous processed layers) the regions of the new layer which have at least one voxel in common with it. We consider the union of the two areas obtained by walking the cube in increasing and decreasing λ as the final detection area over the whole smoothed sub-cube. To avoid losing S/N >  2 regions at the central layers which are only slightly detached from the largest detection area, we found that one should start from an initial mask defined by the collapse of the three “central” segmentation maps, that is, the combined segmentation map of the central layer (containing the largest connected detection area) together with the maps of the two adjacent layers. The selection of the largest detection on this collapsed layer and the percolation at larger and smaller wavelengths following the simple aforementioned constraints allow us to obtain a 3D mask that can be used for the analysis of the extended emission. To avoid the inclusion of spurious signal, we restrict the mask to spaxels with at least three layers along the wavelength direction. The final 3D mask comprises 19253 voxels, extends for 21.25 Å (or 17 layers), and its flux-weighted center is at 4872.7 Å, thus close to the systemic redshift of QSO2 (4872 Å).

4. Results

4.1. Extended emission connecting the quasar pair

We used the smoothed-cube masking described in the previous section to detect extended Lyα emission associated with the quasar pair. The optimally extracted SB map of this Lyα emission is computed by integrating only the signal within the 3D mask as usually done in the literature (e.g., Borisova et al. 2016). Because of the irregular 3D morphology of the mask, each spaxel location of the optimally extracted SB map thus represents the signal integrated along a slightly different range in wavelength. We show this SB map in the left panel of Fig. 3. The extended emission is detected at faint levels (average SB of SB_Lyα  =  1.8 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²) on an area of 191 arcsec², covering the region between the two quasars⁹. At the positions of QSO1 and QSO2 the emission shows slightly higher levels with up to SB_Lyα ∼ 8 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² in proximity to the brighter QSO1. The total luminosity of the extended emission is L_Lyα = 3.2 × 10⁴² erg s⁻¹. The emitting structure shows a projected morphology reminiscent of intergalactic bridges or filaments connecting the two quasars. One bridge extends in the direction connecting the two quasars, while the second passes through the location of the interloper galaxy G (Appendix B). These structures have an average projected width of ∼35 kpc (or 4.5″) at the current spatial resolution and depth. To enable the visualization of the significance of the detection and of the noise properties in our dataset, we overlay the S/N = 3, and 4 contours on the optimally extracted SB map in Fig. 3, while in Appendix C we show a pseudo narrow-band image and a smoothed χ image of the central portion of the wavelength range covered by the 3D mask.

Fig. 3. Lyα emission around the quasar pair in a field of view of about 200 kpc × 200 kpc (or 26″ × 26″). Left: “optimally extracted” Lyα surface brightness map obtained after subtraction of the quasars point-spread-function (PSF) and continuum in the MUSE datacube (details in Sect. 3.1). To highlight the significance of the detected emission, we indicate the contours for S/N = 3 and 4. This image reveals Lyα bridges extending between the quasar pair. Right: flux-weighted velocity-shift map with respect to the systemic redshift of QSO2 obtained from the first order moment of the flux distribution. A velocity gradient between QSO1 and the portion of the nebula southern than QSO2 is evident (Fig. 4). In both panels we indicate the position of the quasars QSO1 and QSO2 prior to PSF subtraction (white circles), and the masked interloper galaxy G, tentatively at z = 0.457 ± 0.001 (more details in Appendix B). Also, to guide the eye, we overlay a grid spaced by 10″ (or 77 kpc). We also highlight the location of a bright knot (white cross) whose SB value is relevant for the discussion in Sect. 5.1, the direction along which we trace the variations in SBLyα in Sect. 5.2.4, and the seeing circle for these observations (bottom left corner).

We also compute the first moment of the flux distribution within the 3D mask, or in other words the flux-weighted velocity shift with respect to the systemic redshift of QSO2. We show the map for the shift in the right panel of Fig. 3. The obtained shifts are in the range −400 km s⁻¹ ≲ Δv ≲ +400 km s⁻¹ and show a gradient along the direction connecting QSO2 to QSO1. Specifically, starting from the southern regions close to QSO2, we see a shift of ∼ − 200 km s⁻¹ which increases to ∼200 km s⁻¹ at the location of QSO1. The northern bridge, instead, shows a shift of ∼ − 400 km s⁻¹ in the vicinity of QSO2 which similarly increases to ∼200 km s⁻¹ at the location of QSO1.

To investigate these velocity gradients and visualize the line profile, we extract spectra along two pseudoslits spanning the two bridges. In particular, for the bridge along the direction connecting the two quasars, we focus on obtaining spectra in five rectangular boxes. For this operation, we simply sum the fluxes layer by layer within each box, without using the aforementioned 3D mask. The rectangles have sides of 1× and 2× the seeing of our observations, that is, 1.66″ × 3.32″, and are placed as shown in the left panel of Fig. 4, starting by centering the first region at the position of QSO1. The extracted 1D spectrum for each region¹⁰ is shown in the central panel of Fig. 4. Each spectrum is normalized at its peak to enable a better comparison of the line emission at the different locations. This panel confirms the presence of a flux-weighted velocity gradient of about 400 km s⁻¹ along the direction QSO1−QSO2, though with slightly different values (−100, 300 km s⁻¹) reflecting the uncertainties in this measure (vertical arrows in right panel of Fig. 4). This gradient is smaller, but comparable with the velocity difference and location of the peaks of the quasars Lyα emission (Δv = 598 ± 98 km s⁻¹; vertical dashed-dotted lines in the right panel of Fig. 4). Similarly, we placed three boxes to cover the second bridge, starting with the first box “a” in vicinity of QSO1 (left panel of Fig. 4). The normalized spectra of these three boxes are shown in the right panel of Fig. 4, confirming the velocity gradient seen in the velocity map, from about −400 km s⁻¹ (box “c”) to 200 km s⁻¹ close to QSO1 (box “a”). Also along this pseudoslit, the velocity gradient roughly spans the velocity difference between the peaks of the quasars Lyα emission.

Fig. 4. Left: “optimally extracted” Lyα surface brightness map as in Fig. 3 with the overlaid pseudoslits used to extract the spectra shown in the central and right panels. We assign an ID (blue) to each box of the pseudoslits. Center: normalized spectra of the Lyα emission along the pseudoslit shown in the left panel with solid lines. Each spectrum is color-coded following the color of its box on the left (details in Sect. 4). The dashed (dotted-dashed) vertical lines show the systemic (peak of the Lyα) redshifts for QSO1 (blue) and QSO2 (magenta). The respective shaded regions indicate the errors on the redshifts, as estimated by SDSS. The velocity shifts Δv are computed with respect to the systemic redshift of QSO2. We indicate with vertical arrows the flux-weighted centroid for each spectrum. The spectrum for box 4 (with no clear emission) is shown in Appendix D (Fig. D.1). Right: same as for the central panel, but for the second pseudoslit shown in the left panel with dotted lines.

The similarity between the two velocity gradients along the two bridges is not surprising. Indeed, (i) the current observations are not extremely deep, leaving space for the presence of more diffuse gas (hence lower levels of Lyman-alpha emission) connecting the currently observed bridges, with the observed gas being only the densest portion of the structure. (ii) Current cosmological simulations of structure formation usually show multiple dense filamentary structures embedded in more diffuse intergalactic gas along the direction of massive halos, or multiple dense structures around massive interacting systems (e.g., Rosdahl & Blaizot 2012; Mandelker et al. 2019).

Further, along both bridges, the difference between the current quasars’ systemic redshifts appear to be wider, Δv = 896 ± 316 km s⁻¹, and seemingly less linked to the observed velocity difference within the extended Lyα emission. Nonetheless, the difference between the two quasars systemic redshifts could be due to the large uncertainties in those measurements (intrinsic uncertainties of 233 km s⁻¹; Table 1). On top of this, the observed gradient and difference with respect to the uncertain quasars’ systemic redshift could encode a mixture of radiative transfer effects, CGM kinematics and intergalactic displacement along the line of sight. We discuss the possible configurations of the system in Sect. 5.

As a next step, we compute the flux-weighted velocity dispersion map as the second moment of the flux distribution for the voxels encompassed by the 3D mask. Figure 5 shows this map, which is clearly noisy due to the narrow spectral range of the detected emission. The obtained velocity dispersions are indeed relatively quiescent, with an average σ_Lyα = 162 km s⁻¹ (or FWHM = 380 km s⁻¹)¹¹. This value is comparable, though lower than the average value observed around individual brighter z ∼ 3 quasars (σ_Lyα = 265 km s⁻¹; Arrigoni Battaia et al. 2019).

Fig. 5. Flux-weighted velocity-dispersion map obtained as the second order moment of the flux distribution within the 3D mask described in Sect. 3.2. The map shows the same field-of-view and uses the same symbols and nomenclature as in Fig. 3.

Aside from the Lyα emission, we did not detect any other extended line emissions associated with the quasar pair down to the depth of the current observations. In particular, we checked the C IVλ1549 and He IIλ1640 expected wavelengths as these two lines can give informations on metallicity, volume density n_H, and speed of shocks (if any; Arrigoni Battaia et al. 2015a). Importantly, the He II/Lyα ratio is sensitive to n_H in a pure recombination scenario, with the ratio decreasing from the expected value of 0.34 (at a temperature T = 2 × 10⁴ K) if He II is not completely doubly ionized (i.e., at high enough densities; Arrigoni Battaia et al. 2015b). Here the observations achieved a 2σ surface brightness limit (in 1 arcsec² aperture) in a single channel (1.25 Å) of SB${}_{\mathrm{C}\mathrm{I}\mathrm{V}}=5.0\times {10}^{-19}$ erg s⁻¹ cm⁻² arcsec⁻² and SB${}_{\mathrm{He}\mathrm{I}\mathrm{I}}=4.6\times {10}^{-19}$ erg s⁻¹ cm⁻² arcsec⁻², respectively for C IV (at 6208.8 Å) and He II (at 6573.5 Å). These slightly deeper sensitivities than at the location of the Lyα are due to the overall system efficiency of the facility which peaks at about 7000 Å. Considering the region where Lyα is detected (191 arcsec²; Fig. 3), we obtain 5σ upper limits for C IV and He II emissions in 5 channels maps, that is, SB_CIV < 2.3 × 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² and SB_{He II} < 2.1 × 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻². We show the maps at these wavelengths in Appendix E. The average C IV/Lyα and He II/Lyα ratios are thus constrained to be < 0.13 (5σ) and < 0.12 (5σ), respectively. Therefore, following Arrigoni Battaia et al. (2015b), if the Lyα emission is only due to recombination, the He II component cannot be fully doubly ionized given the observed low constraints, implying relatively high gas densities (n_H >  0.1 cm⁻³). On the other hand, the low limit in the C IV/Lyα ratio translates to metallicities likely lower than Z_⊙ unless the densities are extremely high, n_H ≫ 1 cm⁻³ (Arrigoni Battaia et al. 2015b). These limits are similar and consistent with what has been usually found for extended Lyα nebulosities around quasars (e.g., Arrigoni Battaia et al. 2015b, 2018; Borisova et al. 2016) and in the so-called Lyman-Alpha Blobs (LAB; Prescott et al. 2009, 2013; Arrigoni Battaia et al. 2015a) down to similar depths. We take into account the limits on these lines in Sect. 5, where we discuss the possible system configurations.

4.2. Gas absorption traced by the quasar pair

As already shown in Fig. 1, we find various absorbers along the two quasar sight-lines. Here we focus on reporting the properties of three absorbers (ABS1, ABS2, ABS3), while we discuss in detail their origin in Sect. 6. In particular, we study a strong absorber (ABS1) along the QSO1 sight-line and close to the systemic redshift of QSO2, and two others (ABS2 and ABS3) found along both sight-lines to the two quasars.

We analyze the absorption features, proceeding as follows. We first model the continuum of each QSO by fitting low-order polynomials to spectral chunks that are free from absorption lines. After the continuum normalization, we model the Lyα absorption lines with Voigt profiles using VPFIT¹² v10.0. Given the spectral resolution of MUSE, we keep the Doppler b parameter fixed at reasonable values while performing the fit of the Lyα absorption lines. These absorptions are near the Lyα emission of the QSOs, and therefore the inferred column densities might be sensitive to the continuum placement. To take this into account we generate a few continuum models for each QSO and repeat the Voigt profile fitting. The dispersion in the resulted N_HI is incorporated in the quoted errors. We model each of the doublet absorption lines (C IVλ1548, 1550 and N Vλ1238, 1242) using a double Gaussian profile. The sigma values of the two Gaussians of a doublet are tied to be the same, and the wavelength ratio is fixed at the value given by the atomic tables. We also allow the equivalent widths (EW) ratio of the two lines in a doublet to vary to take into account the possible saturation effect. We note that such models are not sensitive to continuum placements since the lines are reasonably narrow. We then obtain the lower limits of the column densities using the EWs of the lines and assuming the linear part of the curve-of-growth. For the non-detected transitions we use the S/N at the position of the lines to calculate the upper limit on the EWs. We further convert such limits to upper limits on column densities using the linear part of the curve-of-growth. The fits performed along the two sight-lines are reported in Figs. 6 and 7, while all the derived parameters are listed in Table 2.

Fig. 6. Profiles of the absorbers along the QSO1 sight-line at the H I Lyα, N V, Si II, C II, and C IV lines. The black histograms show the continuum normalized data, while the orange lines are the sum of all the Gaussian components of the best fit. The locations of the absorbers considered in this analysis are highlighted with vertical dashed lines, that is, ABS1 (magenta), ABS2 (lime), and ABS3 (gray). The zero velocity is set to the redshift of the C IV strongest component. Table 2 gives all the fit parameters. For ABS1, we exclude the best fit solution for NHI (orange) following physical arguments (Sect. 6). We thus show the two extreme alternative fits with NHI and b values that are favored by the current data, that is, b = 100 km s−1, log(NHI/cm−2) = 17 as a dashed-dotted blue line, and b = 200 km s−1, log(NHI/cm−2) = 15 as a dotted green curve.

Fig. 7. Profiles of the absorbers along the QSO2 sight-line at the H I Lyα line, N V, Si II, C II, and C IV lines. The black histograms show the continuum normalized data, while the orange lines are the sum of all the Gaussian components of the best fit. The location of the absorbers considered are highlighted with vertical dashed lines, that is, ABS1 (magenta; not present along this sight-line), ABS2 (lime), and ABS3 (gray). The zero velocity is set to the redshift of the C IV strongest component of ABS1 along the QSO1 sight-line, as in Fig. 6. Table 2 gives all the relevant fit parameters.

For ABS1 we estimate log(N_{H I} = 15−17, allowing the Doppler b parameter to vary uniformly between 100 and 200 km s⁻¹, with smaller b at higher N_{H I}. Allowing for even smaller b parameters down to 50 km s⁻¹ increases the goodness of the fit (χ² decreases from ∼9 to ∼4). However, these small b values require very large column densities (log(N_{H I}/cm⁻²) > 18), which are disfavored by the lack of absorption at the location of low-ion transitions, like Si IIλ1260 or C IIλ1335. As an additional test, we check for the presence of an associated Lyman limit system (LLS; log(N_{H I}/cm⁻²) > 17.2) by looking for its 912 Å break in the LRIS data of the Quasars Probing Quasars database (QPQ; Findlay et al. 2018). We find no clear evidence for the break, confirming that the values log(N_{H I} = 15−17 are favored. For completeness, in Table 2 and in Fig. 6 we report examples for b = 50 (solid orange line), 100 (dashed-dotted blue line), and 200 km s⁻¹ (dotted green line).

ABS1 has associated absorption in N Vλ1240 and C IVλ1549. The absorption at the C IV wavelength is best fit by two components (log($N_{\mathrm{C}\mathrm{I}\mathrm{V}}^{\mathrm{strong}}/{\mathrm{cm}}^{-2})>14.9$; log($N_{\mathrm{C}\mathrm{I}\mathrm{V}}^{\mathrm{weak}}/{\mathrm{cm}}^{-2})>14.5$), while the one at N V can be fit by a single Gaussian line (log(N_{N V}/cm⁻²) > 15.5) at the current spectral resolution. The fit of the H I absorption places ABS1 at z = 3.005 ± 0.001. This redshift is at Δv = −230 ± 140 km s⁻¹ from the systemic redshift of QSO2. The two C IV components show velocity shifts of +173 km s⁻¹ and −66 km s⁻¹ with respect to the H I absorption, respectively for the strong and weak components. The N V is redshifted by +45 km s⁻¹. These shifts justify the large b parameter allowed during the fit of the H I absorption.

The metal absorptions show relatively quiescent widths with the two C IV components being characterized by a velocity dispersion σ^strong = 76 ± 4 km s⁻¹ (or 0.39 ± 0.02 Å) and σ^weak = 58 ± 10 km s⁻¹ (or 0.30 ± 0.05 Å), and the N V by σ = 215 ± 10 km s⁻¹ (or 0.89 ± 0.04 Å), after correcting for the MUSE spectral resolution. The larger value for N V could be partially due to the superposition of a second unresolved component. However, if we assume that the two unresolved components share roughly the same σ, we would get σ ∼ 150 km s⁻¹, still larger than C IV. Keeping in mind the large uncertainties in the determination of the two quasars redshifts, our fit overall suggests that the strong absorber ABS1 could be associated with QSO2. We cannot completely exclude that this absorber is due to intervening associated gas to QSO1 (at Δv = −1120 ± 140 km s⁻¹), but the narrowness of its metal absorptions rules out the scenario in which this gas is outflowing at small distances from QSO1. Therefore, ABS1 is most likely produced by gas at least on CGM scales (around QSO1, QSO2 or the system comprising the two quasars; Sect. 6).

ABS2 is placed at Δv = −900 km s⁻¹ from QSO2, and, conversely, has log(N_{H I}/cm⁻²) ∼ 14 along both sight-lines with no other absorption lines detected. This absorber is thus very similar to Lyα forest clouds (Meiksin 2009).

The fit of ABS3, located at Δv  =   − 2800 km s⁻¹ from QSO2, shows high log(N_{H I}/cm⁻²) ∼ 19 along both sight-lines (b  =  50 km s⁻¹). However, as for ABS1, the absence of low ion transitions possibly implies smaller values of N_{H I}, that is, log(N_{H I}/cm⁻²) ∼ 15−17 (b = 200−100 km s⁻¹). Also in this case we look for an associated 912 Å break in the QPQ database (Findlay et al. 2018), finding no evidence for a LLS. The values log(N_{H I}/cm⁻²) ∼ 15−17 are thus favored also in this case. Further, ABS3 shows strong C IV absorption only toward QSO2 (log(N_{C IV}/cm⁻²) > 14.7). For both ABS2 and ABS3 we do not find associated galaxies at their corresponding redshifts (Sect. 6).

Higher spectral resolution observations are required to firmly constrain the properties of all these three absorbers. Nevertheless, in Sect. 6 we show that already our current data allow us to roughly infer their nature.

5. The powering of the extended Lyα emission

Three powering mechanisms could be responsible for the extended Lyα emission detected around quasars: photoionization by the quasar (e.g., Heckman et al. 1991; Haiman & Rees 2001; Weidinger et al. 2005), scattering of Lyα photons from the quasar (e.g., Dijkstra & Loeb 2008), or shocks due to the quasar activity (e.g., Cai et al. 2017). These mechanisms do not exclude each other, and could possibly act together. We explore in turn their contributions, if any, in the system studied here, using analytical considerations. The modeling of these mechanisms in a cosmological context (e.g., Gronke & Bird 2017) is beyond the scope of this work.

First, we focus on fast quasar winds. This phenomenon has been so far traced in emission only out to a few tens of kpc from the central engine (e.g., Harrison et al. 2014), and is usually manifested in emission lines with FWHM ≳ 1000 km s⁻¹ and velocity shifts of at least few hundreds of km s⁻¹ (e.g., Mullaney et al. 2013; Cai et al. 2017). The extended Lyα emission detected in our data differs substantially as it shows a relatively quiescent line profile with an average velocity dispersion σ_Lyα = 162 km s⁻¹ (or FWHM = 380 km s⁻¹). Considering that this value is not corrected for the instrument spectral resolution and that resonant scattering of Lyα photons could broaden the line, it is safe to say that fast winds do not play a major role in shaping the Lyα extended structure and its emission level that we observe. This is in agreement with what has been routinely found with short exposures for extended Lyα emission around z ∼ 3 quasars (Arrigoni Battaia et al. 2019).

Secondly, we consider a photoionization scenario from both quasars. Indeed, the case in which only one of the quasars shines on the gas seems to be ruled out by: (i) the higher Lyα fluxes in proximity of each quasar, and (ii) the absence of emission on large scales in the NW direction from QSO1 and SE direction from QSO2, where most likely only the contribution of one quasar (modulo its opening angle and presence of gas) is relevant. We thus explore the quasar pair photoionization scenario in the two limiting regimes for the recombination emission: optically thin (N_HI  ≪  10^17.2 cm⁻²) or optically thick (N_HI  ≫  10^17.2 cm⁻²) gas to the Lyman continuum photons. We do this in two steps. First, we show some expectations by following the model for cool gas around quasars introduced by Hennawi & Prochaska (2013), and then we model the system using the photoionization code Cloudy (version 17.01), last described in Ferland et al. (2017).

5.1. Analytical estimates for the extended Lyα emission

In the framework of Hennawi & Prochaska (2013), the cool (T ∼ 10⁴ K) gas is organized in clouds characterized by a single uniform hydrogen volume density n_H, a cloud covering factor f_C, and a hydrogen column density N_H. Knowing these quantities and the luminosity of a quasar, one can estimate the Lyα emission at a distance R from it.

Specifically, in the optically thick case, the Lyα SB scales with the luminosity of the central source and should decrease as R⁻² with increasing distance from a quasar (see Hennawi & Prochaska 2013 for the derivation of the formula):

$$\begin{array}{cc}\hfill {\mathrm{SB}}_{\mathrm{Ly}\alpha }^{\mathrm{thick}}& =5.7\times {10}^{-17}{\left(\frac{1+z}{4.014}\right)}^{-4}\left(\frac{f_{\mathrm{C}}}{1.0}\right){\left(\frac{R}{50\mathrm{kpc}}\right)}^{-2}\hfill \\ \hfill & \times \left(\frac{L_{{\nu }_{\mathrm{LL}}}}{7.6\times {10}^{29}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}}\right)\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{arcsec}}^{-2},\hfill \end{array}$$(1)

where L_νLL is the specific luminosity at the Lyman edge. To obtain this luminosity for the two quasars, we assume a spectral energy distribution (SED) which follows the form L_ν  =  L_νLL(ν/ν_LL)^α_UV blueward of the Lyman limit ν_LL, and adopt a slope of α_UV  =   − 1.7 (Lusso et al. 2015). As done in Arrigoni Battaia et al. (2015b), L_νLL is computed by integrating the Lusso et al. (2015) composite spectrum against the SDSS filter curve to give the correct i-band magnitude of the two quasars (as listed in Table 1). We find $L_{{\nu }_{\mathrm{LL}}}^{\mathrm{QSO}1}=7.6\times {10}^{29}$ erg s⁻¹ Hz⁻¹ and $L_{{\nu }_{\mathrm{LL}}}^{\mathrm{QSO}2}=5.7\times {10}^{29}$ erg s⁻¹ Hz⁻¹ for QSO1 and QSO2, respectively. The two quasars have a bolometric luminosity of $L_{\mathrm{bol}}^{\mathrm{QSO}1}=1.5\times {10}^{46}$ erg s⁻¹ and $L_{\mathrm{bol}}^{\mathrm{QSO}2}=1.1\times {10}^{46}$ erg s⁻¹, when using a standard quasar SED template as described in Sect. 5.2.1.

We promptly demonstrate that the optically thick scenario is unlikely to be in place in this system. We can indeed explore different configurations (e.g., different distances between the quasars), and add up the contribution to ${\mathrm{SB}}_{\mathrm{Ly}\alpha }^{\mathrm{thick}}$ given by Eq. (1) for each quasar. For this discussion, we focus on the region of the bridge indicated by a white cross in Fig. 3, which is characterized by SB_Lyα = 3.5 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻², and assume in Eq. (1) an average redshift of 3.014 and f_C = 1, unless specified. We first consider the case in which the redshift difference is mainly tracing peculiar velocities and thus the distance between the quasars and the region considered is roughly given by the projected distance, R_{x − QSO1} = 57 kpc and R_{x − QSO2} = 31 kpc, respectively. Following Eq. (1), the sum of the contributions due to the two quasars would then give SB${}_{\mathrm{Ly}\alpha }^{\mathrm{thick}}=1.5\times {10}^{-16}$ erg s⁻¹ cm⁻² arcsec⁻², which is about 40× higher than the observed value. Even if we consider a factor of two larger distances, we would obtain an SB${}_{\mathrm{Ly}\alpha }^{\mathrm{thick}}$ 11× higher than observed. This can be reconciled by invoking a very low covering factor (f_C ∼ 0.02−0.09), obscuration of the quasars in the direction of the emitting gas, or larger distances between the two quasars and the observed gas. Low covering factors for the emitting clouds are disfavored as the emission would have looked much clumpier than observed (e.g., Arrigoni Battaia et al. 2015a). Conversely, with the current dataset we cannot firmly verify if the two quasars are strongly obscured (by e.g., dust on small scales or their host galaxy) so that only few percent of their luminosity shines on the gas. However, we obtain a crude estimate of the intrinsic extinction E(B − V) affecting the two quasars by fitting their spectra with a reddened version of the expected power law of the composite SDSS spectrum (α_opt = −0.46, Vanden Berk et al. 2001; Sect. 5.2.1), after normalising it to the continuum at 8200Å. The power law is reddened using an SMC extinction curve (Pei 1992), in which E(B − V) is a free parameter and R_V  =  2.93 is fixed. For both quasars we found E(B − V)< 0.06, indicating that the spectra of these quasars do not show significant intrinsic extinction along our line of sight. Nonetheless, the lack of strong obscuration and of dust has to be directly explored with follow-up observations, for example, in the near-infrared (e.g., Banerji et al. 2015) and submillimeter (e.g., Venemans et al. 2017) regimes. We next explore larger distances.

The uncertain redshift difference between the two quasars could reflect their distance within the Hubble flow. In this configuration, the zone considered for our estimates would then sit at much larger distances than previously considered. If we assume the region to be at half way between the two quasar systemic redshifts, R_{x − QSO1}  =  R_{x − QSO2}  =  1.45 Mpc, we obtain SB${}_{\mathrm{Ly}\alpha }^{\mathrm{thick}}=1.2\times {10}^{-19}$ erg s⁻¹ cm⁻² arcsec⁻². This value is about 30× smaller than the observed SB. Considering shorter distances given by the peak of the Lyα emission, R_{x − QSO1}  =  R_{x − QSO2}  =  0.95 Mpc, would only increase the SB to SB${}_{\mathrm{Ly}\alpha }^{\mathrm{thick}}=2.8\times {10}^{-19}$ erg s⁻¹ cm⁻² arcsec⁻². To match the value at the considered position, the two quasars should lie at a distance R_{x − QSO1}  =  R_{x − QSO2} = 267 kpc, which would translate to a very small velocity or redshift difference, that is, Δv  =  223 km s⁻¹ or Δz  =  0.002 (comparable to the redshift error). Even in this configuration, a fully optically thick scenario is ruled out for small distances from each quasar (if they shine on the gas).

We then focus on the optically thin case, which has been shown to only depend on the gas physical properties (e.g., n_H, N_H) provided the radiation is intense enough to keep the gas sufficiently ionized to be optically thin to the Lyman continuum photons (Hennawi & Prochaska 2013):

$$\begin{array}{cc}\hfill {\mathrm{SB}}_{\mathrm{Ly}\alpha }^{\mathrm{thin}}& =1.8\times {10}^{-18}{\left(\frac{1+z}{4.014}\right)}^{-4}\left(\frac{f_{\mathrm{C}}}{1.0}\right)\left(\frac{n_{\mathrm{H}}}{0.24{\mathrm{cm}}^{-3}}\right)\hfill \\ \hfill & \times \left(\frac{N_{\mathrm{H}}}{{10}^{20.5}{\mathrm{cm}}^{-2}}\right)\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{arcsec}}^{-2}.\hfill \end{array}$$(2)

As shown in Eq. (2), if we assume the median N_H value from absorption studies of quasar halos (log N_H = 20.5; Lau et al. 2016) and a plausible n_H for CGM gas¹³, the optically thin scenario can match the observed average SB_Lyα. This first order calculation holds only if the two quasars are able to keep the gas ionized enough to be optically thin to the ionizing radiation. As we demonstrate in the next section, this is not the case for distances R ≳ 100 kpc, and so a fully optically thin scenario holds only if the system extent is similar to or slightly larger than the observed projected distance.

5.2. Photoionization models for the extended Lyα emission

In the following sections, we construct photoionization models assuming different configurations of the quasar pair to test which one is more likely given the constraints on the different extended line emissions reported in Sect. 4. Specifically, we will base our investigation on the estimates presented in the previous section, and thus focus on three configurations: (i) the quasars sits at a separation similar to the projected distance, (ii) the quasars are within the Hubble flow with a separation of 2.9 Mpc, and (iii) the quasars are placed at an intermediate distance between the two aforementioned cases. Before describing the Cloudy calculation, we first describe the parametrization of the two quasars spectral energy distributions (SEDs), and discuss how we consider the impact of resonant scattering.

5.2.1. The assumed SED for the two quasars

For the quasars’ SEDs we adopt the same assumptions as in Arrigoni Battaia et al. (2015b) because we do not have complete coverage of the quasars’ spectra. The only exception to the modeling of Arrigoni Battaia et al. (2015b), is the assumption of the simple power-laws measured by Vanden Berk et al. (2001) and Lusso et al. (2015) in the rest-frame optical and UV, respectively. In the mid-IR part of the SED we assume the composite spectra by Richards et al. (2006). In Fig. 8 we show the shape of the assumed SED for both QSO1 and QSO2, together with their MUSE spectra. We provide below a summary of the power-laws used to parametrized the SEDs:

$$\begin{array}{c}\hfill f_{\nu }\propto \{\begin{array}{cc}{\nu }^{{\alpha }_{\mathrm{opt}}},\hfill & \mathrm{i}f0.11\mathrm{Ryd}\le h\nu \le 1\mathrm{Ryd}\hfill \\ {\nu }^{{\alpha }_{\mathrm{EUV}}},\hfill & \mathrm{i}f1\mathrm{Ryd}\le h\nu \le 30\mathrm{Ryd}\hfill \\ {\nu }^{\alpha },\hfill & \mathrm{i}f30\mathrm{Ryd}\le h\nu <2\mathrm{keV}\hfill \\ {\nu }^{{\alpha }_{\mathrm{X}}},\hfill & \mathrm{i}f2\mathrm{keV}\le h\nu <100\mathrm{keV}\hfill \\ {\nu }^{{\alpha }_{\mathrm{HX}}},\hfill & \mathrm{i}fh\nu \ge 100\mathrm{keV},\hfill \end{array}\end{array}$$(3)

where α_opt = −0.46 (Vanden Berk et al. 2001), α_EUV = −1.7 (Lusso et al. 2015), α = −1.65 (i.e., obtained to match an α_OX = −1.5; Strateva et al. 2005), α_X = −1, and α_HX = −2. These assumptions are regarded as standard in photoionization modeling of active galactic nuclei (AGN; e.g., Baskin et al. 2014).

Fig. 8. Spectral energy distribution (SED) of QSO1 (blue) and QSO2 (orange), used as incident radiation in the Cloudy calculations. We compare the models with the available MUSE data (lighter color for each quasar). Left panel: the vertical lines indicate the energies used to define the different power-laws (Sect. 5.2.1). Right panel: zoomed in version of the box highlighted in the left panel. The emission from the wavelength ranges around the C IV and Lyα line locations are masked to prevent contributions from resonantly scattered photons, as done in Arrigoni Battaia et al. (2015b).

5.2.2. Approximating the impact of resonant scattering

Because of the large optical depth at line center (e.g., Gould & Weinberg 1996), Lyα photon propagation should be affected by substantial resonant scattering under most astrophysical conditions. Even at very close separation from a quasar, the gas can be found to be optically thick to the Lyα transition (i.e., N_{H I} ≳ 10¹⁴ cm⁻²; e.g., Gallimore et al. 1999). Hence, a Lyα photon typically experiences a large number of scatterings before escaping the system or clouds in which it starts to interact (e.g., Neufeld 1990; Dijkstra et al. 2006).

However, it is usually found and assumed that the scattered Lyα line photons from the quasar do not contribute significantly to the SB_Lyα surrounding quasars on large scales, that is, ≳100 kpc (e.g., Hennawi & Prochaska 2013; Cantalupo et al. 2014; Arrigoni Battaia et al. 2015b). Indeed, the quasar’s Lyα photons very efficiently diffuse in velocity space. Consequently, the vast majority of these photons escape the system at very small scales (≲10 kpc), without propagating to larger distances (e.g., Dijkstra et al. 2006; Verhamme et al. 2006).

In this work we thus use a twofold approach. First, we neglect the contribution due to resonant scattering in the Cloudy calculations, such that we can mimic the expected negligible contribution of scattering on large scales and have “clean” predictions. To achieve this, we mask the quasars’ input SEDs at the Lyα line location as done in Arrigoni Battaia et al. (2015b) (Fig. 8). This method still allows us to account for the scattered Lyα photons arising from the diffuse continuum produced by the gas itself, which, however, appear to be negligible in our calculations¹⁴. Second, we introduce an approximate estimate for the contribution from resonant scattering of quasar Lyα photons, which is found to be more relevant on small scales. To compute this estimate, we need to know: (i) the fraction of the quasar’s Lyα photons seen by a parcel of gas in the nebula, or the probability that the quasar’s Lyα photons scatters in the direction of a portion of the nebula, and (ii) the probability of scattering and escaping the nebula in the direction of the observer. For each photon, both these probabilities can be written as P = W(cos(θ))e^−τ_esc, and are thus governed by the phase function W(cos(θ)) (or angular redistribution function, which parametrizes the probability of a photon to be scattered in a certain direction) and by the optical depth for the Lyα photons τ_esc ∼ N_HIσ_α(ν, T), where σ_α(ν, T) is the cross section for the Lyα scattering (e.g., Stenflo 1980; Laursen et al. 2009; Dijkstra 2017). For simplicity, we assume: (i) W(cos(θ)) ∼ 0.5 as it corresponds to the most probable value of cos(θ), and (ii) a similar optical depth between quasar and nebula, and nebula and observer. Also, as the cross section depends on the gas motions, we assume the gas to be in infall toward the quasars with a velocity of 200 km s⁻¹ as shown in cosmological simulations (e.g., Goerdt & Ceverino 2015). We then compute the estimate for the SB_Lyα due to scattering as

$$\begin{array}{c}\hfill {\mathrm{SB}}_{\mathrm{Ly}\alpha }^{\mathrm{scatt};\mathrm{QSO}}=\frac{f_{\mathrm{conv}}}{4\pi {(1+z)}^4}{R}^{-2}{\displaystyle {\int }_{\lambda =4865\mathrm{\AA }}^{4888\mathrm{\AA }}}{L}_{\mathrm{Lya};\mathrm{QSO}}(\lambda )P{(\lambda ,N_{\mathrm{HI}},T)}^2\mathrm{d}\lambda \end{array}$$(4)

where f_conv is the conversion from steradians to arcsec², R is the distance from the quasar, and L_Lya; QSO(λ) is the observed quasar luminosity spectrum. We convolve this with the aforementioned probability to observe a quasar Lyα photon after scattering, P(λ, N_HI, T)², and use the observed wavelength range [4865, 4888] in which we see extended Lyα emission (e.g., Fig. 4). As reference, if we integrate the quasars’ spectra in this range without applying the probability we get L_Lya; QSO1 = 1.17 × 10⁴³ erg s⁻¹ and L_Lya; QSO2 = 6.74 × 10⁴² erg s⁻¹ for QSO1 and QSO2, respectively. We note that these luminosities are similar to the luminosity of the extended structure (L_Lyα = 3.2 × 10⁴² erg s⁻¹). The use of the observed spectrum L_Lya; QSO(λ) is conservative because a non-negligible fraction of the quasars’ photons could have been absorbed in the system and along our line-of-sight before reaching the observer. We use the N_HI and T of the Cloudy calculations in the formula of P.

This treatment is very crude and has to be regarded as an indicative reference, given that we use a fixed set of parameters for θ, the relative gas velocity, and P. Only a Monte Carlo simulation of Lyα radiative transfer applied to cosmological simulations of quasar pairs could properly handle this problem and give more detailed insights. However, Monte Carlo simulations of Lyα radiative transfer are beyond the scope of this work, and, in any case, none of our results should depend strongly on this effect given the large extent of the system studied.

We stress that similar considerations also apply to the resonant C IV line (e.g., Berg et al. 2019). However, in this work we neglect the contribution of resonant scattering to the C IV line, since this process should be less efficient for the C IV photons than for the Lyα, due to the much lower abundance of metals. Taking into account that resonant scattering is important mainly at small distances from the quasars, neglecting this effect does not affect the main results of this work.

5.2.3. Photoionization models for a single quasar

To have a reference for the subsequent modeling of the quasar pair, we first show the results of photoionization of gas illuminated by a single faint quasar, QSO1. On top of the assumption for the quasar SED and for the resonant scattering already presented, we select the model parameter grid for this visualization as follow. We assume (i) a standard plane-parallel geometry for the slab, (ii) a fixed volume density n_H = 1, 0.1, 0.01 cm⁻³ whose values should encompass possible values in the quasar CGM, (iii) a fixed metallicity Z = 0.1 Z_⊙ close to the value seen in absorption studies around z ∼ 2 quasars (∼0.3 Z_⊙; Lau et al. 2016), and (iv) we stop our calculations when a total Hydrogen column density N_H = 10^20.5 cm⁻² is reached. This value is the median N_H estimated for absorbers around z ∼ 2 quasars out to an impact parameter of 300 kpc (Lau et al. 2016). We then place the slab of gas at increasing distance from the quasar to show how this would affect the predicted Lyα emission. Specifically, we place the slab at 30 different distances spaced in logarithmic bins between 20 and 1500 kpc.

Figure 9 shows the results of this calculation for the N_HI (top panel) and the SB_Lyα (lower panel) as a function of distance from QSO1. The two regimes described in Sect. 5.1, optically thin and optically thick to the ionizing radiation, are readily evident (the dotted gray line in the top panel indicates N_HI = 10^17.2 cm⁻²). A slab can be optically thin further away from the quasar than a denser slab, following the relation

$$\begin{array}{c}\hfill R_{{\mathrm{n}}_{\mathrm{H}}^{\mathrm{smaller}}}=\sqrt{n_{\mathrm{H}}^{\mathrm{larger}}/n_{\mathrm{H}}^{\mathrm{smaller}}}{R}_{{\mathrm{n}}_{\mathrm{H}}^{\mathrm{larger}}}.\end{array}$$(5)

Fig. 9. Cloudy predictions for plane parallel slabs with log(NH/cm−2) = 20.5, photoionized by a single quasar with the characteristic SED of QSO1. The slabs – characterized by uniform nH – are placed at increasing distance from the quasar. Top: column density of H I as a function of distance from the quasar. The horizontal dotted line indicates the threshold between the optically thin and thick regimes. For nH = 0.1 cm−3 the quasar is able to keep the gas ionized up to ∼200 kpc. Bottom: predicted SBLyα as a function of distance from the quasar. The dotted lines indicate Cloudy models which take into account the presence of the z = 3 UVB by Haardt & Madau (2012). The gray shaded region shows the range of observed SBLyα.

This can be easily obtained by comparing the number of ionizing photons at the two different distances or, in other words, by finding at which distance the ionization parameter U¹⁵ is the same for models with different densities. We note that the N_HI saturates to the total gas content on short distances after the models transition from optically thin to optically thick.

The prediction for the optically thin regime does not follow exactly the aforementioned relation SB_Lyα ∝ N_Hn_H as Cloudy takes into account both temperature changes of the recombination coefficients and the contribution to the Lyα emission from cooling. Both these phenomenon increase with distance from the quasar as the temperature drops increasing the recombination efficiency (e.g., Storey & Hummer 1995) and collisional coefficients (T  ∼  10^4.2 K; e.g., Raymond et al. 1976; Wiersma et al. 2009). In the optically thick regime, SB_Lyα ∝ L_νLL, scaling with the distance following R⁻², as expected. For > 100 kpc, the presence of additional ionizing photons from the metagalactic ultraviolet background (UVB; e.g., Haardt & Madau 2012) introduces mild differences in the predicted SB_Lyα, and very slight changes in the ionized fraction. This is illustrated by the deviation from the predicted R⁻² relation toward higher SB_Lyα of the dotted curves, which show the Cloudy models run with the UVB from Haardt & Madau (2012) at z  =  3. We do not show the scattering contribution here since it seems irrelevant at these distances (e.g., the dashed lines in Fig. 10). From Fig. 9, it is already clear that relatively dense gas (n_H >  0.1 cm⁻³) is needed to produce the high levels of SB_Lyα detected around the observed quasar pair in the short exposures with MUSE. Heckman et al. (1991), Cantalupo et al. (2014), Hennawi et al. (2015), Arrigoni Battaia et al. (2015b) have already shown that such dense gas is needed to explain the emission around individual quasars.

Fig. 10. Cloudy predictions for plane parallel slabs with total log(NH/cm−2) = 20.5 and nH = 0.5 cm−3, illuminated by the quasar pair QSO1 and QSO2 placed at a separation equal to their observed projected distance (89 kpc). Top left: column density of H I, NH I, as a function of distance. For each model the data-points are color-coded by their respective ionization parameter. In this scenario all the models are optically thin to the ionizing radiation, that is, NH I < 1017.2 cm−2. Top right: comparison of the observed (black line with shaded 1σ error) and predicted SBLyα (blue line). The brown and dashed lines indicate the contribution due to scattering of Lyα photons from the quasars, as explained in Sect. 5.2.2. The red dotted line indicates the total SBLyα summing up the Cloudy prediction and the scattering contribution. The thin black dotted lines are the observed SBLyα along the directions NE for QSO1 and SE for QSO2, at angles 52 and 142° east from north, respectively (details in Sect. 5.2.4). Bottom left: comparison of the observed (black line is the 2σ upper limit) and the predicted (blue line) He II/Lyα ratio as function of the distance from the quasars. The green shaded area represents the parameter space allowed by the observations. Bottom right: comparison of the observed (black line is the 2σ upper limit) and the predicted (blue line) C IV/Lyα ratio as function of the distance from the quasars. The green shaded area represents the parameter space allowed by the observations. The vertical dotted lines in each panel indicate the position of the two quasars, while the striped yellow regions represent the zones used to normalize the quasar PSF, characterized by large uncertainties and, therefore, not considered in the analysis. The dotted red line represents the ratio corrected for the presence of Lyα scattering.

5.2.4. Photoionization models for a quasar pair at the observed projected distance

In this section we present the modeling of a photoionization scenario in which the two quasars sit at a separation similar to the observed projected distance (89 kpc), and both illuminate the gas responsible for the extended Lyα emission. In this framework, the two quasars are likely in a merger phase which would explain the observed velocity shift of the Lyα emission and the uncertain difference in velocities between the quasar systemics. Large peculiar velocities are thus in play.

Our model grid covers distances [10, 80] kpc from each quasar in steps of 10 kpc. This is achieved by normalizing each quasar spectrum at different values of f_νLL depending on its distance from the slab. To be conservative, we do not consider distances smaller than 10 kpc, because of uncertainties due to the quasars’ PSF subtraction, and because in such close proximity to the quasars we expect density variations and effects due to, for instance, the interstellar medium of the host galaxy. For simplicity, we assume (i) a plane parallel geometry, (ii) a fixed volume density n_H = 0.5 cm⁻³, (iii) Z = 0.1 Z_⊙, and (iv) we stop our calculations when N_H = 10^20.5 cm⁻² is reached.

In Fig. 10 we show the prediction of this set of models in terms of: N_HI (top-left), SB_Lyα (top-right), and line ratios He II/Lyα (bottom-left) and C IV/Lyα (bottom-right). For the observational data, we extract the average Lyα emission along the direction connecting the two quasars, using a slit with width 2× the seeing of our observations (solid black line). It is important to note that, of course, in proximity of the quasars there are variations in the SB_Lyα depending on the direction along which we place the slit. To appreciate the difference in Lyα profiles along different directions close to the two quasars, we show how the SB_Lyα behaves along the NE direction from QSO1 and the SE from QSO2 at angles of 52 and 142° east of north (black dotted lines in top-right plot). These two directions have been chosen because they are perpendicular and parallel, respectively, to the direction connecting the two quasars (142° east of north). For the ratios, we divide the 2σ SB limits per layer at the He II and C IV locations (Sect. 4.1) by the aforementioned SB_Lyα within the two quasars. The allowed parameter space is indicated by the green shaded region.

Figure 10 shows that our model grid can reproduce the roughly flat SB_Lyα and the line ratios observed, with the He II/Lyα ratio possibly starting to show some tension with our simple modeling. In this configuration, the gas emitting Lyα emission is highly ionized and thus optically thin to the ionizing radiation. This is due to the relatively high ionization parameter U at each location, log U >  −1.2. The presence of Lyα resonant scattering appears to be non-negligible on scales R <  20kpc from the quasar, and could help in explaining the low He II/Lyα ratios observed at these locations. Finally, we stress that the assumption of a constant N_H value along all of the emitting bridge, implies a total mass of cool (T ∼ 10⁴ K) gas of M_cool = f_CAN_Hm_p/X = 3.9 × 10¹⁰ M_⊙, where A is the area covered by the bridge, m_p is the proton mass, and X = 0.76 is the Hydrogen mass fraction (Hennawi & Prochaska 2013). If the two quasars are hosted by a halo of M = 10^12.5 M_⊙ (average halo hosting quasars at these redshifts; White et al. 2012), the detected cool gas mass would represent 9.4% of the total gas mass within the halo, after removing the mass expected to be in stars, M_* = (8.3 ± 2.8)×10¹⁰ M_⊙ (Moster et al. 2018). Once taken into account that our observations are not sensitive to very diffuse gas (n_H ∼ 10⁻² cm⁻³), this estimate seems surprisingly close to the fraction of cool gas seen in similar massive halos in current cosmological simulations (15%; e.g., Cantalupo et al. 2014; Arrigoni Battaia et al. 2018).

5.2.5. Photoionization models for a quasar pair within the Hubble flow

Now, we assume a photoionization scenario in which the two quasars sit at their systemic redshifts, and thus at a distance of 2.9 Mpc. As this stretched configuration would place the filamentary emission along our line of sight, we again assume that both quasars illuminate the gas responsible for the observed extended Lyα emission. Because of the finite speed of light, this assumption requires that QSO2 has been active for at least 18.9 Myr, while QSO1 for 9.4 Myr. These values seem reasonable given the current estimates for quasars’ lifetimes (e.g., Martini 2004; Eilers et al. 2017; Schmidt et al. 2018; Khrykin et al. 2019)¹⁶. In this framework, the two quasar halos are not yet strongly interacting and the velocity shift of the Lyα line (Fig. 3) would be a mixture of complex radiative transfer effects and velocities tracing the Hubble flow.

Our model grid covers distances [10, 2890] kpc from each quasar in steps of about 100 kpc. As in Sect. 5.2.4, this is achieved by normalizing each quasar spectrum at different values of f_νLL depending on its distance from the slab. For simplicity, we assume (i) a plane parallel geometry, (ii) a fixed volume density n_H  =  0.5 cm⁻³, (iii) Z = 0.1 Z_⊙, and (iv) we stop our calculations when N_H = 10^20.5 cm⁻² is reached. These N_H, n_H, and Z are likely too high for the average cloud in the IGM (e.g., Meiksin 2009), but here we are interested in conservatively high values which should produce the highest signal observable as the models remain optically thick (Sect. 5.1 and Hennawi & Prochaska 2013). As we explore large distances from the quasars, we also run models with the UVB of Haardt & Madau (2012).

In Fig. 11 we show the predictions of this set of models in the same observables as in Fig. 10. We plot the predictions of the clean Cloudy models (blue solid line), the Cloudy models with the UVB as additional source (blue dotted line), and the Cloudy models plus our approximated contribution of the Lyα scattering (red dotted line). As expected from the single source model presented in Sect. 5.2.3, the Lyα emission shows its maximum levels SB_Lyα ≈ 3 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² at the transition between the optically thin and thick regimes (R ∼ 100 kpc from each quasar), with the expected decline as R⁻² in the optically thick regime, reaching the minimum (SB_Lyα ≈ 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²) at the half distance between the quasars. The contribution of ionizing photons from the UVB almost precisely doubles the Lyα emission at this location. In this “Hubble-flow” scenario, the He II and C IV line emissions will be extremely faint and already at the limit of current facilities capabilities for close separations (R ∼ 100 kpc) from each quasar. In this regard, it is interesting to note that our approximate treatment of scattering creates a region with peak He II/Lyα at a distance of about 100 kpc from the quasar. This effect remains to be verified with detailed radiative transfer simulations. Also, the ratio C IV/Lyα peaks at the same location as He II/Lyα. Its trend, however, is not mainly driven by the assumption on the Lyα scattering, but by the higher excitation of Carbon on smaller scales.

Fig. 11. Cloudy prediction for plane parallel slabs with total log(NH/cm−2) = 20.5 and nH = 0.5 cm−3, illuminated by the quasar pair QSO1 and QSO2, placed at a separation of 2.9 Mpc, as derived from their systemic redshifts. Top left: column density of H I, NH I, as a function of distance. The data-points are color-coded by the ionization parameter. In this scenario most of the models are optically thick to the ionizing radiation, NH I ≫ 1017.2 cm−2, with only the CGM regions of the two quasars characterized by optically thin gas. Top right: predicted SBLyα for the Cloudy models without (solid blue line) and with (dotted blue line) the UVB. The red dotted line indicates the total SBLyα, summing up the Cloudy prediction and the scattering contribution estimated following Sect. 5.2.2. The gray shaded region shows the range of observed SBLyα. Bottom left: predicted He II/Lyα ratio as function of the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. Bottom right: predicted C IV/Lyα ratio as function the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. In each of the bottom panels, the dotted red lines represent the ratios corrected for the presence of Lyα scattering as modeled in Sect. 5.2.2, while the green horizontal line indicates the local 2σ upper limit on each ratio (Fig. 10). In all four panels, the vertical dotted lines indicate the position of the two quasars, while the striped yellow regions show the zones used to normalize the quasars PSFs.

In this configuration, it is difficult to directly compare our photoionization models with the observations as complex projection effects can drastically change the predicted curves. For this reason, we do not attempt to plot our data in Fig. 11, but we only show the observed range of SB_Lyα and the local 2σ upper limit on the ratios He II/Lyα and C IV/Lyα. Nevertheless, from our models, it is clear that the observed emission would be dominated by gas at small distances from the two quasars, that is, in their CGM. In this scenario, we would thus expect to see two nebulae sitting at the systemic redshift of the quasars, and thus to find at least some overlapping emission showing double peaks, with each peak sitting at the systemic of the two quasars or at the redshift of the Lyα peak of the two quasars. We inspect our data for such a signature, finding a signal at the systemic of QSO1 only in close proximity to its location (2″ or 10 projected kpc) and along the direction connecting the two quasars (within box 2), as shown by the blue line in Fig. 4. This signature is very concentrated spatially (<1″), and for this reason we suspect that it is due to a compact object. Also, there is tentative evidence for a double peak in close proximity of QSO2 (black line in Fig. 4). This double peak is also extremely localized and could be due to radiative transfer effects at this location. We thus conclude that there are no obvious signatures of a superposition of two nebulae at different redshifts.

Finally, we stress that, in this framework, the direction of the discovered bridges of Lyα, stretching between the two quasars, would be due to a very improbable chance alignment of dense structures in the two distinct CGMs. This alignment is quite unlikely also because of the absence of additional extended emission in other directions. We thus argue that this scenario is not able to reproduce the observations.

5.2.6. Photoionization models for a quasar pair at an intermediate distance

As already discussed in Sect. 5.1, an interesting configuration places the two quasars at an intermediate distance with respect to the two extremes considered so far. Specifically, we consider a distance of 600 kpc. Indeed, our analytical estimates suggest that optically thick models would be able to reproduce the observed emission if the two quasars sit at ∼300 kpc from the central region of the observed bridges. As this configuration also stretches considerably the bridges along our line of sight, we assume that both quasars shine on the gas. In this framework the two quasar halos are probably approaching and the velocity shift of the Lyα line (Fig. 3) should be interpreted as a mixture of complex radiative transfer effects, velocities tracing the approaching quasar halos, and extent of the structure along our line of sight.

Our model grid covers distances [10, 590] kpc from each quasar in steps of about 50 kpc. We make the same assumptions as in Sects. 5.2.4 and 5.2.5, when computing this grid of models. We assume (i) a plane parallel geometry, (ii) a fixed volume density n_H = 0.5 cm⁻³, (iii) Z = 0.1 Z_⊙, and (iv) stop the calculations when N_H = 10^20.5 cm⁻² is reached. These N_H, and Z are likely too high for the average cloud in the IGM (e.g., Meiksin 2009), but they represent well the properties of absorbing gas seen around high-z quasars (e.g., Lau et al. 2016). Following the modeling of a single quasar (Sect. 5.2.3), the n_H is chosen large enough to allow for a match of the observed SB_Lyα.

We show the predictions of this set of models for N_HI (top-left), SB_Lyα (top-right), He II/Lyα (bottom-left) and C IV/Lyα (bottom-right) in Fig. 12. The color scheme is the same as in Fig. 11. As expected from the analytical modeling and from the single source calculation in Sect. 5.2.3, the Lyα emission is predicted to be roughly at the same level of SB_Lyα ≈ 2.5 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² throughout all the extent of the bridges. This happens even though the model transitions between the optically thin and thick regimes at around R ∼ 100 kpc from each quasar. At small distances (R ≲ 50 kpc), our calculation shows that the contribution from scattering could be important. Furthermore, in this scenario, the contribution of ionizing photons from the UVB of Haardt & Madau (2012), is irrelevant (all three curves fall on top of each other for distances larger than 50 kpc). As already seen in the Hubble-flow scenario (Sect. 5.2.5), the He II and C IV line emissions are extremely faint, and basically barely observable with current instruments. Interestingly, our approximate treatment of scattering creates also in this scenario a region with peak He II/Lyα at a distance of 50−100 kpc from the quasars. Detailed radiative transfer simulations will be able to verify this effect. The C IV/Lyα line ratio peaks also at the same location as the He II/Lyα one. As we discussed in the previous scenarios, this trend is driven by the higher excitation of Carbon on smaller scales, and not by the assumption on the Lyα scattering.

Fig. 12. Cloudy prediction for plane parallel slabs with total log(NH/cm−2) = 20.5 and nH = 0.5 cm−3, illuminated by the quasar pair QSO1 and QSO2, placed at an intermediate separation of 600 kpc. Top left: column density of H I, NH I, as a function of distance, color-coded by ionization parameter. In this scenario, a central region of about 300 kpc is optically thick to the ionizing radiation, NH I ≫ 1017.2 cm−2, while the CGM regions of the two quasars are optically thin. Top right: predicted SBLyα for the Cloudy models without (solid blue line) and with (dotted blue line) the UVB. The red dotted line indicates the total SBLyα, summing up the Cloudy prediction and the scattering contribution estimated following Sect. 5.2.2. The predicted SBLyα is basically flat at a value of SBLyα = 2.5 × 10−18 erg s−1 cm−2 arcsec−2, well within the observed range (gray shaded region). Bottom left: predicted He II/Lyα ratio as function of the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. Bottom right: predicted C IV/Lyα ratio as function of the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. In each of the bottom panels, the dotted red lines represent the ratios corrected for the presence of Lyα scattering as modeled in Sect. 5.2.2, while the green horizontal line indicates the local 2σ upper limit on each ratio (Fig. 10). The vertical dotted lines in all four panels indicate the position of the two quasars, while the striped yellow regions give the zones used to normalize the quasars PSFs.

As for the Hubble-flow scenario (Sect. 5.2.5), it is difficult to compare our photoionization models with the observations, as complex projection effects (e.g., absorption from the structure itself for both emitted and impinging radiation) can drastically change the predicted curves. For this reason, we do not plot our data in Fig. 11, but show the available information as done in Fig. 11. Nevertheless, from our models, it is clear that the observed flux can be equally due to emission from dense CGM and IGM surrounding the two quasars, with the central region of the bridge possibly characterized by optically thick gas. If this is the case, we would expect to see differences in the Lyα line shape as we move along the bridge, with the presence of double peaks or strong asymmetries (e.g., Neufeld 1990; Laursen et al. 2009) in its central region. We cannot exclude the presence of these features below the current MUSE spectral resolution (FWHM ≈ 2.85 Å or 175 km s⁻¹ at 4870 Å). Deep observations at higher spectral resolution with available IFUs, e.g., MEGARA (Gil de Paz et al. 2016) or KCWI (Morrissey et al. 2012), or longslit spectroscopy could help to clarify the shape of the Lyα emission, and assess if optically thick gas is present in this structure.

Finally, we can calculate a rough estimate for the gas mass in this extended structure by assuming a cylindrical geometry for each bridge with extent 600 kpc and diameter 35 kpc. In this case, we can compute the total cool gas mass as M_cool = Vf_V(n_H/0.5 cm⁻³)m_p/X = 8.8 × 10¹²(f_V/1.0) M_⊙, where V is the volume covered by one of the bridges, m_p is the proton mass, X = 0.76 is the Hydrogen mass fraction, and f_V is the volume filling factor (e.g., Hennawi & Prochaska 2013). Multiplying by the number of bridges, we thus obtain a total cool gas mass of M_cool = 1.8 × 10¹³(f_V/1.0) M_⊙. As the parcels with high densities n_H = 0.5 cm⁻³ are only the tracer of the structure, that is, the volume filling factor of parcels with n_H = 0.5 cm⁻³ is expected to be much lower than unity for gas on such large scales¹⁷, this estimate has to be regarded as an upper limit for the total cool gas mass along the observed structure. Confirming the presence of high densities n_H = 0.5 cm⁻³ within the IGM would imply finding parcels of gas similar to interstellar medium densities spread out along filaments. This scenario sounds plausible if we are tracing emission close to faint undetected galaxies, but quite unrealistic at the moment for “pure” IGM (n_H ≲ 10⁻² cm⁻²; e.g., Meiksin 2009).

6. Modeling the absorbers

In Sect. 4.2 we presented the observed properties of absorbers ABS1, ABS2, and ABS3. In this section we discuss their nature in different system configurations.

6.1. ABS1: metal enriched CGM or IGM absorber

The redshift of ABS1 suggests a link with QSO2. However, peculiar motions could mimic such an association, and ABS1 could be related to QSO1, QSO2 or be in the IGM. Furthermore, the two similar components seen at the C IV line could be due to different structures along the line of sight (Fig. 6). We constrain the nature of ABS1 by constructing photoionization models with Cloudy under three different system configurations. We briefly outline here the models and the results, while we present them in detail in Appendix F.

The photoionization models need to match the observed column densities reported in Table 2, and the Lyα emission at the location of the absorber (i.e., QSO1). The joint constraints from absorption and emission are key in assessing the physical properties of the gas (e.g., n_H), and thus its configuration (Hennawi et al. 2015). Unfortunately, as our PSF subtraction algorithm is not reliable within the 1 arcsec² region around QSO1, we can only assume conservative limits for the Lyα emission, that is, below the 5σ value per channel, which is equivalent to a SB_Lyα below 1.75 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻². The same applies also to QSO2. As we show in Appendix F, the loose constraint on the Lyα emission does not allow a firm evaluation of the absorber’s location.

Specifically, the three system configurations here considered are as follows: (i) ABS1 is only illuminated by the QSO1’s radiation, (ii) ABS1 sees the radiation from both QSO1 and QSO2, and the two quasars lies at a separation similar to the observed projected distance, and (iii) ABS1 is illuminated by both quasars, with QSO1 and QSO2 separated by 600 kpc. For all the models we consider the presence of the UVB as additional ionizing source. Importantly, we focus on these three configurations as they are allowed by the modeling of the extended Lyα emission shown in Sect. 5.

As explained in detail in Appendix F, in all the three configurations we find that ABS1 is a cool (4.1 ≲ log(T/K) ≲ 4.4), metal enriched (Z >  0.3 Z_⊙) absorber, located on CGM or IGM scales around the quasar pair. The relatively high metallicity is constrain by the presence of the strong N_{N V} absorption. Further, our analysis suggests that the location of ABS1 should be characterized by −1.7 ≲ log U ≲ −0.6. The current data, however, do not allow us to put stringent constraints on its precise position due to its loosely constrained n_H and N_H. Finally, we note that the resulting characteristics of ABS1 are similar to the absorbers usually studied along background sightlines piercing the halo of a foreground quasar (e.g., Lau et al. 2016). Those absorbers, however, show lower U than ABS1 (log U <  −1.7; e.g., Fig. 6 in Lau et al. 2016). Indeed, in those cases the quasar pairs are not physically related and the absorbers should not receive much of the radiation from the background quasar. Observations at higher spectral resolution together with deeper IFU data have the potential to firmly constrain the physical properties of ABS1, and thus its position.

6.2. ABS2 and ABS3: CGM or IGM coherent structures along the quasar pair sight-line

As reported in Sect. 4.2, ABS2 (log(N_{H I}/cm⁻²) ∼ 14) and ABS3 (log(N_{H I}/cm⁻²) = 15−17) appear on both quasars sight-lines, suggesting they trace coherent structures on ∼100 kpc (the projected separation between the two quasars). At the current depth of the observations, these absorbers are not associated to any continuum source in the MUSE field-of-view, nor to Lyα-emitting galaxies at the absorption redshift. We evaluate a 5σ upper limit for the counterpart (if any) in a seeing aperture, finding L_Lyα <  3.0 × 10⁴¹ erg s⁻¹ ($\sim 0.1L_{\mathrm{Ly}\alpha }^{\ast }$ of Ciardullo et al. 2012). Intriguingly, all these properties are very similar to the absorber detected at Δv = −710 km s⁻¹ along the line of sight to the quasar pair observed by Cai et al. (2018) with KCWI.

The wider MUSE wavelength range allowed us to detect the presence of strong C IV absorption for ABS3. This C IV detection is only visible along the QSO2 sight-line (log(N_{C IV}/cm⁻²) > 14.7). The presence of this relatively strong high-ionization metal line absorption might indicate that this portion of ABS3 is located at a closer distance to QSO2 (or strong ionizing sources, e.g. a shock front) than the remainder of the structure. The absence of absorption at the N V wavelength might indicate a low metallicity for ABS3. The values log(N_{H I}/cm⁻²) = 15−17 require a relatively large Doppler b parameter (200−100 km s⁻¹), which could be due to turbulences in expanding shells around the quasar pair. Data at higher spectral resolution are needed to explore this occurrence and to firmly constrain the properties of ABS2 and ABS3, which are likely CGM or IGM structures coherently extending in front of the quasar pair.

7. Summary and conclusions

Recent observations of extended Lyα emission around individual quasars suggest that multiple quasar systems are surrounded by more extended and rich structures (Hennawi et al. 2015; Arrigoni Battaia et al. 2018, 2019). In an effort to characterize the Lyα emission from CGM and IGM scales, we have initiated a “fast” survey (45 min on source) of z  ∼  3 quasar pairs with MUSE on VLT, complementing the work by Cai et al. (2018). In this study we focus on the first targeted faint z  ∼  3 quasar pair, SDSS J113502.03−022110.9 – SDSS J113502.50−022120.1 (z  =  3.020−3.008; i  =  21.84, 22.15), separated by 11.6″ (or 89 projected kpc).

We discovered the presence of filamentary Lyα emission connecting the two quasars at an average surface brightness of SB_Lyα  =  1.8 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻². Using photoionization models constrained with the information on Lyα, He IIλ1640, and C IVλ1548 line emissions, we show that the emitting structures could be explained as intergalactic bridges with an extent between ∼89 up to 600 kpc. The faintness of the two quasars and the high levels of Lyα emission seem to rule out a 2.9 Mpc extent for the bridges along our line-of-sight, as it could be inferred from the difference between the systemic quasars redshifts. The intergalactic nature of the emission is also supported by the narrowness of the Lyα line (σ_Lyα  =  162 km s⁻¹). At the current spatial resolution and surface brightness limit, the projected average width of the bridges is ∼35 kpc.

Additionally, we studied three absorbers found along the two quasar sight-lines. We detect strong absorption in H I, N V, and C IV along the background quasar sight-line, which we interpret as due to at least two components of cool (T  ∼  10⁴ K), metal enriched (Z >  0.3 Z_⊙), and relatively ionized circumgalactic or intergalactic gas characterized by an ionization parameter of −1.7 ≲ log U ≲ −0.6. Two additional H I absorbers are detected along both quasars sight-lines, at   ∼   − 900 and −2800 km s⁻¹ from the system. The H I absorber at −2800 km s⁻¹ has associated C IV absorption along only the foreground quasar sight-line. These two absorbers are not associated to any continuum or Lyα emitters within the MUSE field of view, possibly tracing large-scale structures or expanding shells in front of the quasar pair.

The observations presented in this study confirm that intergalactic bridges can be observed even with short exposure times, if peculiar or overdense systems are targeted (e.g., multiple AGN systems). This is likely due to the presence of dense (n_H  ∼  0.5 cm⁻³) gas on large scales coupled with the ionizing radiation originating from multiple sources. Deep high spectral resolution observations of such systems could firmly constrain the physical properties of the emitting gas and impinging ionizing continuum, providing a new leverage to improve current cosmological simulations of structure formation.

Acknowledgments

We thank Guinevere Kauffmann for providing comments on an early version of this work. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 0100.A-0045(A). We acknowledge the role of the ESO staff in providing high-quality service-mode observations, which are making this project feasible in a shorter time-scale. A.M. is supported by the Dunlap Fellowship through an endowment established by the David Dunlap family and the University of Toronto. This work made use of matplotlib (Hunter 2007).

References

Appendix A: The point spread function of our MUSE data

To model the PSF of our data, which is needed to subtract the unresolved emission from the two quasars (Sect. 3.1), we rely on the only bright star (i = 16.2; r = 16.4) within our observations field of view, 2MASS J11350307−0220597 (Cutri et al. 2003). This star has been so far classified as single point source in all the available catalogs we explored, that is, the 2MASS All Sky Catalog of point sources (Cutri et al. 2003), the AllWISE Source Catalog (Wright et al. 2010), the 14th Data Release of the Sloan Digital Sky Survey (SDSS DR14; Abolfathi et al. 2018). This star is not saturated in our data.

A faint red source (r ≈ 22) at about 4.4″ is present on the right side of this star, which is clearly visible in the white-light image (Fig. A.1). Given the red spectrum, this faint source does not contribute significantly at the wavelength of interest for the Lyα emission. We however remove this low-level contaminant by replacing in each layer the values at its position with the values at the symmetrical position with respect to the star centroid. This is to avoid the introduction of any systematic in the PSF subtraction, and in the subsequent extraction of the Lyα signal that we seek. The result can be visually inspected in the right panel of Fig. A.1, while we show the normalized profile of this “corrected” star out to five arcseconds in Fig. A.2 (open green squares).

Fig. A.1. Five arcseconds circle cutout of the white-light image of the bright star used as PSF in this work. Left: white-light image without post-processing. A faint source on the right part of the star’s PSF is clearly evident. Right: white-light image after replacing the faint source values with the symmetric portion of the MUSE dataset. After this correction, the PSF is well behaved at any wavelength out to five arcseconds (Fig. A.2 and Appendix A for details).

Fig. A.2. Normalized profile of the bright star used as PSF in this work. The open green squares show the normalized profile derived from the white-light image shown in the right panel of Fig. A.1, while the small dots show the normalized profiles for the star within the 17 MUSE layers encompassed by the obtained 3D mask for the extended Lyα emission (Sect. 3.2). The red dashed line is the best-fit Moffat profile to the white-light image data (β = 2.5 and FWHM = 1.66″). Given the brightness of the used star, there is very good agreement between the profile obtained from the white-light image and the individual layers.

The star profile is well fit by a Moffat function with β = 2.5 and FWHM = 1.66″. The value for β is in agreement with the usually assumed value for the MUSE instrument (β = 2.8; e.g., Bacon et al. 2017). For completeness, in the same plot we also show the star profile for each of the 17 layers of the 3D mask of the extended Lyα emission (small gray dots) built in Sect. 3.2. It is clear that the star PSF is defined at high S/N out to five arcseconds even in the individual layers, showing a profile consistent with the white-light image. In our analysis we adopt a normalized version of the star data layer-by-layer (after removal of the faint source) as empirical PSF.

Appendix B: The superposed galaxy at a lower redshift

The MUSE observations unveil the presence of a faint galaxy located in projection between the quasar pair. In Fig. B.1 we show the i-band image extracted from the MUSE datacube using the transmission curve of the corresponding SDSS filter (Fukugita et al. 1996). The image encompasses the same field of view of Fig. 3. We indicate the position of the two quasars, of the faint galaxy and of the 2σ isophote of the Lyα emission. The galaxy has a i magnitude of i = 23.82 ± 0.03 when extracted in a circle with radius of 1″.

Fig. B.1. SDSS i-band image extracted from the final MUSE datacube for the same field-of-view shown in Fig. 3. We indicate the position of QSO1, QSO2 and of a faint galaxy unveiled in projection between the targeted quasar pair. This galaxy is an interloper at a different redshift, tentatively at z = 0.457 ± 0.001 (Appendix B and Fig. B.2). For comparison purposes, we also show the 2σ isophote for the extended Lyα emission (black). To guide the eye we overlay a grid spaced by 10″ (or 77 kpc), as done in Fig. 3.

We show the spectrum of this faint galaxy in Fig. B.2. It is evident a relatively strong line emission at λ = 5435 Å, F = (2.1 ± 0.1)×10⁻¹⁷ erg s⁻¹ cm⁻². Given the absence of any other signature useful to identify the galaxy redshift, we cannot firmly place this galaxy in a cosmological context. Its redshift, however, is surely not close to the quasar pair as there are no known strong line emissions at a rest-frame wavelength of ∼1800 Å. Further, the galaxy morphology seems resolved even with the large seeing of these observations, possibly hinting at a low-redshift nature for this object. For reference, we compute its redshift by assuming the line emission to be [O II]λ3729. We find z = 0.457 ± 0.001. If this galaxy is indeed a foreground object, its dust and gas could absorb the higher redshift Lyα photons of interest to us. Deeper spectroscopy could unveil the nature of this galaxy and quantify its effect on the extended Lyα emission.

Fig. B.2. One-dimensional spectrum (black) of the faint galaxy shown in Fig. B.1, as extracted using a circular aperture with radius 1″ from the MUSE data. The red spectrum indicates the error vector. The vertical dashed blue line indicates the position of the only detected line emission (λ = 5435 Å). This galaxy is not associated to the quasar pair as there are no known strong line emissions at a rest-frame wavelength of ∼1800 Å. If this line is [O II]λ3729, this galaxy would be at z = 0.457 ± 0.001.

Appendix C: Narrow-band and χ maps of the Lyα bridge

In Sect. 4.1 we show the optimally extracted map of the extended Lyα emission connecting the quasar pair. For completeness and comparison purposes, we present here also a pseudo narrow-band image. Specifically, we collapsed the five layers (or 6.25 Å) of the final MUSE datacube centered at the wavelength of 4872.7 Å. This wavelength corresponds to the central layer of the 3D mask obtained in Sect. 4.1. To avoid introducing too large of a sky noise, the wavelength range of the pseudo narrow-band is chosen to be small and comparable to the width of the Lyα line in the central part of the observed structure. We caution that the chosen width does not encompass the whole velocity range spanned by the aforementioned 3D mask. The top panel of Fig. C.1 shows the SB map obtained in this way after a smoothing with a Gaussian kernel with FWHM = 1.66″ (i.e., the seeing of the observations). The extended Lyα emission connecting the two bridges is readily visible.

Fig. C.1. Top: pseudo 6.25 Å (5 layers) narrow-band map centered at the central wavelength (4872.7 Å) of the 3D mask of the Lyα emission obtained in Sect. 3.2. The map, obtained after PSF and continuum subtraction, shows a 57″ × 57″ (or 438 kpc × 438 kpc) FoV and it is color coded following the Lyα surface brightness. Bottom: χsmth map for the same wavelength range as in the top panel, and obtained using a Gaussian kernel with FWHM = 1.66″ (i.e. similar to the seeing), as explained in Appendix C. To guide the eye, in both panels we overlay a grid spaced by 15″ (or 115 kpc) and we indicate the position of QSO1 and QSO2 prior to their PSF subtraction. In both panels, the interloper galaxy G is indicated with its contour.

Further, we visualize the noise properties of this map and the significance of the detection by constructing a smoothed χ image of the same dataset following the recipe in Hennawi & Prochaska (2013) and Arrigoni Battaia et al. (2015a), for a Gaussian kernel with FWHM  =  1.66″. The smoothed χ image is obtained by dividing the smoothed data shown in the left panel of Fig. C.1, I_smth, by the smoothed sigma image σ_smth computed by propagating the variance image of the unsmoothed data (details in Arrigoni Battaia et al. 2015a). The bottom panel of Fig. C.1 shows this smoothed χ image after masking a circular region of radius 4″ around the bright star 2MASS J11350307−0220597. This map reveals that the extended Lyα emission is detected at relatively high significance, and that the noise behaves quite well throughout all the field of view.

Appendix D: Spectrum of box 4 along the pseudoslit

Here we present the spectrum of box 4 along the pseudoslit used in Sect. 4.1 (Fig. 4). Figure D.1 shows this spectrum in physical units. We omitted these data from Fig. 4 as it would have made that normalized plot harder to read. The faint level of emission at this location is in agreement with the optimally extracted map presented in Sect. 4.1.

Fig. D.1. Spectrum of box 4 along the pseudoslit used in the analysis of Fig. 4. The dashed (dotted-dashed) vertical lines show the systemic (peak of the Lyα) redshifts for QSO1 (blue) and QSO2 (magenta). The respective shaded regions indicate the error on the redshift as estimated by SDSS.

Appendix E: χ maps at the C IV and He II wavelengths

In Sect. 4.1 we quoted upper limits for the C IV and He II extended line emissions. Here we show a cut of the final MUSE datacube at their expected observed wavelengths given the flux-weighted center of the Lyα emission, 6208.8 Å and 6573.5 Å respectively. In particular, we construct smoothed χ images following the method described in Fig. C.1. These maps have the potential of better visualizing the presence of extended emission.

Figure E.1 presents the two smoothed χ maps obtained using a Gaussian kernel with FWHM = 1.66″. The white circles indicate the position of the two quasars prior to their PSF subtraction. We mask a circular region of radius 4″ around the bright star 2MASS J11350307−0220597. As mentioned in Sect. 4.1, there is no evidence for extended emission at these wavelengths.

Fig. E.1. Top: χsmth map of 57″ × 57″ (or 438 kpc × 438 kpc) FoV for 6.25 Å (5 layers) centered at the wavelength expected for the C IV line emission (6208.8 Å) given the center of the 3D mask for the Lyα emission. The map is obtained after PSF and continuum subtraction using a Gaussian kernel with FWHM = 1.66″ (i.e., similar to the seeing) as explained in Appendix C. Bottom: same as the top panel, but centered at the wavelength expected for the He II line emission (6573.5 Å). Both maps, covering the corresponding velocity range of Fig. C.1, do not reveal the presence of any clear detection of extended C IV or He II line emission associated with the extended Lyα emission. In both panels, we overlay a grid spaced by 15″ (or 115 kpc) and we indicate with circles the position of QSO1 and QSO2 prior to their PSF subtraction. In both panels, the interloper galaxy G is indicated with its contour.

Appendix F: Modeling the absorber ABS1

In Sect. 6.1 we summarize the results of our photoionization models concerning ABS1 in three different system configurations. In this appendix we present in detail the assumptions and the predictions of these calculations.

For simplicity, we assume the following for all the Cloudy models here discussed: (i) a plane parallel geometry, (ii) three values of fixed volume density n_H  =  0.1, 0.01, 0.001 cm⁻³, (iii) three values of fixed metallicity Z  =  0.1, 0.5, 1 Z_⊙, and (iv) a column density stopping criteria (N_H  =  10^20.5 cm⁻²). We do not consider higher values for n_H as these would result in higher SB_Lyα than the assumed upper limit for the emission (e.g., Fig. 9). Therefore, all models are already in agreement with the limits on the emission at the absorber position (SB_Lyα <  1.75 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻²; Sect. 6.1). We further note that all the models presented in this section include the presence of the UVB at z ∼ 3 (Haardt & Madau 2012). The three system configurations probed are as follows.

First, as ABS1 is only seen along the QSO1 sight-line, we assume the absorber to be illuminated only by QSO1. In this framework, QSO2 is obscured in the direction of ABS1, that is, ABS1 is not within the ionizing “cones” of QSO2. We thus run Cloudy models assuming as input only the continuum of QSO1 and the UVB, and consider distances in the range [20, 1500] kpc.

Figure F.1 shows how the column densities for the different ions change as a function of distance from QSO1 in the grid of models at solar metallicity. From left to right, we show the results for n_H = 0.1, 0.01, 0.001 cm⁻³, respectively. We note that the decrease in n_H (and thus increase in U) causes the predicted curves to be shifted toward larger distances (as photoionization models are self-similar in U). This shift follows Eq. (5) so that, for example, the curves for n_H = 0.01 cm⁻³ are at ∼3.2 times larger distances than the ones for n_H = 0.1 cm⁻³. Allowing for higher n_H values (> 0.1 cm⁻³) would require the absorber to be at small distances (<40 kpc) from QSO1. This seems to be ruled out not only by the Lyα levels implied by higher n_H, but also by the relatively quiescent kinematics of the metal absorptions (Table 2).

Fig. F.1. Predictions of photoionization models for the absorber ABS1 with Z = Z⊙, in the case it is illuminated by QSO1 and the UVB (see Appendix F for details on the models assumptions). The predicted column densities for H I, C IV, N V, C II, and Si II are plotted as a function of distance from QSO1. The horizontal green lines with arrows indicate the observational limits for the same metal ions, while the blue shaded regions show the observational limits for H I. The green hatched boxes indicate the regions where the models matched the observations. The model curves are color-coded following the predicted temperature. The gray dashed lines represent the curves for N V for Z = 0.5 Z⊙. The models in agreement with the observations are characterized by 4.1 ≲ log(T/K) ≲ 4.4 and −1.7 ≲ log U ≲ −0.7. To match the N V absorption, the metallicity should be Z >  0.3 Z⊙.

Fig. F.2. Predictions of photoionization models for the absorber ABS1 with Z = Z⊙, in the case it is illuminated by QSO1, QSO2 (placed at a projected distance of 89 kpc) and the UVB (details on the models assumptions in Appendix F). All symbols and colors are explained in the caption of Fig. F.1. Similarly to Fig. F.1, the models in agreement with the observations are characterized by 4.1 ≲ log(T/K) ≲ 4.4, −1.7 ≲ log U ≲ −0.7, and Z >  0.3 Z⊙.

In each panel, the curves are color-coded by their temperature, and we indicate the observed limits on the metal ions column densities as horizontal green lines with arrows, and the limits on Hydrogen as blue shaded regions. The green hatched boxes indicate the regions where the models match the observations. It is evident that these hatched regions encompass models with the same temperature in all three panels of Fig. F.1, that is, 4.1 ≲ log(T/K) ≲ 4.4. These also translate to the same ionizing parameters −1.7 ≲ log U ≲ −0.7. Furthermore, the observed lower limit on N_{N V} requires the models to be relatively metal enriched (Z >  0.3 Z_⊙). As example, we show the predictions for N_{N V} at Z = 0.5 Z_⊙ as a dashed gray line in each panel. These photoionization models thus predict that ABS1 is a cool (4.1 ≲ log(T/K) ≲ 4.4) absorber, already enriched, and located at a distance $40\sqrt{0.1{\mathrm{cm}}^{-3}/n_{\mathrm{H}}}$ kpc ≲ $R\lesssim 130\sqrt{0.1{\mathrm{cm}}^{-3}/n_{\mathrm{H}}}$ kpc from QSO1, where the ionization parameter is constrained to be −1.7 ≲ log U ≲ −0.7. It is thus clear that, in this configuration, ABS1 could be located from the CGM of QSO1 out to the IGM (even at Mpc distances from QSO1).

As we found good agreement between our models and the observed Lyα emission for a configuration in which the two quasars sit at their projected distance (Sect. 5.2.4), we assume that the distance along the line-of-sight between QSO1 and QSO2 is negligible and that both illuminate ABS1. Therefore, in the next step we run Cloudy models assuming as input the continua of both QSO1 and QSO2, scaled accordingly to their distance from the absorber. In particular, we consider distances in the range [20, 1500] kpc from QSO1, and distances $d_{\mathrm{QSO}2}=\sqrt{d_{\mathrm{QSO}1}^2+{89}^2}$ kpc from QSO2.

Figure F.2 shows how the column densities for the different ions change as a function of distance from QSO1 in the grid of models at solar metallicity. The addition of QSO2 at a small projected distance from QSO1 only slightly changes the predictions of the previously considered configuration, with the absorber now positioned at slightly larger distances from QSO1. The location of the hatched boxes is shifted roughly following Eq. (5). Indeed, on scales comparable to the distance between the two quasars, Eq. (5) is no longer strictly valid. Specifically, for n_H = 0.1 cm⁻³ ABS1 sits at 40 kpc ≲ R_QSO1 ≲ 170 kpc and 98 kpc ≲ R_QSO2 ≲ 190 kpc. For n_H = 0.01 cm⁻³, we find 150 kpc ≲ R_QSO1 ≲ 570 kpc and 170 kpc ≲ R_QSO2 ≲ 580 kpc, while for n_H = 0.001 cm⁻³, we get 490 kpc ≲ R_QSO1 ≲ 1600 kpc and 500 kpc ≲ R_QSO2 ≲ 1610 kpc. For the lowest-density grid, larger distances are also allowed. All the selected ranges where the models agree with the observations correspond to temperatures 4.1 ≲ log(T/K) ≲ 4.4 and ionizing parameters −1.7 ≲ log U ≲ −0.7. We again find that the models require metallicities Z >  0.3 Z_⊙ in order to match the observed absorptions in the metal ions, especially N_{N V}. For this configuration, ABS1 could be located from the CGM of the system comprising QSO1 and QSO2 out to the IGM.

Last, a configuration in which the two quasars sit at an intermediate distance of ∼600 kpc, can also explain the observed levels of Lyα emission (Sect. 5.2.6). Therefore, we also model this case for the illumination of ABS1. In particular, we run Cloudy models assuming as input the continua of both QSO1 and QSO2, considering distances 20 ≤ d_QSO1 ≤ 1500 kpc, and accordingly $R_{\mathrm{QSO}2}=\sqrt{{(d_{\mathrm{QSO}1}-600)}^2+{89}^2}$ kpc.

Figure F.3 shows how the column densities for the different ions change depending on the distance from QSO1, in the grid of models at solar metallicity. In each panel we indicate the locations of QSO2, at 600 kpc from QSO1. However, the absorber is never closer than 89 kpc from QSO2, as per the formula above. From left to right, we show the results for n_H = 0.1, 0.01, 0.001 cm⁻³, respectively. It is clear that the presence of QSO2 at 600 kpc creates a more complex behavior of the curves with respect to the previous two configurations. A high ionized region can now be seen around QSO2 as well. For this reason, there is not a straightforward formula to derive the distance at which the absorber is located depending on its density. We can derive d_QSO1 by looking for the position where the same U is achieved in all three different density cases. For n_H = 0.1 cm⁻³, ABS1 would be located at 40 kpc ≲ R_QSO1 ≲ 140 kpc and 470 kpc ≲ R_QSO2 ≲ 570 kpc, or at 500 kpc ≲ R_QSO1 ≲ 680 kpc and 90 kpc ≲ R_QSO2 ≲ 140 kpc. These distances correspond to 4.1 ≲ log(T/K) ≲ 4.4 and −1.7 ≲ log U ≲ −0.7, or 4.1 ≲ log(T/K) ≲ 4.2 and −1.7 ≲ log U ≲ −1.5. For n_H = 0.01 cm⁻³, ABS1 would be located at 120 kpc ≲ R_QSO1 ≲ 530 kpc and 114 kpc ≲ R_QSO2 ≲ 488 kpc, or at 710 kpc ≲ R_QSO1 ≲ 1000 kpc and 145 kpc ≲ R_QSO2 ≲ 500 kpc. These locations result in 4.2 ≲ log(T/K) ≲ 4.4 and −1.2 ≲ log U ≲ −0.6, or 4.1 ≲ log(T/K) ≲ 4.3 and −1.8 log U ≲ −0.9. Finally, for n_H = 0.001 cm⁻³, ABS1 would be located at 900 kpc ≲ R_QSO1 ≲ 1500 kpc and 400 kpc ≲ R_QSO2 ≲ 900 kpc. Larger distances (not modeled here) are allowed in the lowest density case. All the aforementioned models give similar ranges for the temperature 4.1 ≲ log(T/K) ≲ 4.4, and ionization parameter −1.2 ≲ log U ≲ −0.6. For each panel, the selected distances thus reflect similar T and U as the two configurations previously discussed. It is thus clear that, in this configuration, ABS1 could be located from the CGM of QSO1 or QSO2 out to the IGM. As in the two previous cases, the observed lower limit on N_{N V} requires the models to be relatively enriched, with metallicities Z >  0.3 Z_⊙.

Fig. F.3. Predictions of photoionization models for the absorber ABS1 with Z = Z⊙, in the case where it is illuminated by QSO1, QSO2 (placed at 600 kpc from QSO1) and the UVB (details on the models assumptions in Appendix F). All symbols and colors are explained in the caption of Fig. F.1. Similarly to Figs. F.1 and F.3, the models in agreement with the observations are characterized by 4.1 ≲ log(T/K) ≲ 4.4 and −1.7 ≲ log U ≲ −0.7 (in Appendix F we report the detailed ranges for each plot). The metallicity is constrained to be Z >  0.3 Z⊙.

All Tables

All Figures

Fig. 1. One-dimensional spectra (black) for the two quasars of the pair, QSO1 (top) and QSO2 (bottom), as extracted from the final MUSE data, using a circular aperture with radius 2′′. The red spectra indicate the error vectors. The vertical dashed blue (magenta) lines indicate the position of important line emissions at the systemic redshift of QSO1 (QSO2). For both objects, we show in the inset plots a zoomed version of the spectrum at the location of the Lyα and C IV lines to highlight the presence of interesting absorptions. The vertical dashed gray and green lines indicate the location of an H Lyα absorptions present within both QSO1 and QSO2 spectra, and C IV absorption along the QSO2 sight-line. The fit to these lines is shown in Sect. 4.2 and Figs. 6 and 7. Residuals due to frequent sky lines are evident at wavelengths >7000 Å.

In the text

Fig. 2. White-light image of the observed 57″ × 57″ field of view. We indicate the position of QSO1, QSO2, the star 2MASS J11350307−0220597 used to compute the point spread function of our data (Appendix A), and an interloper galaxy “G” (tentatively at z = 0.457 ± 0.001; Appendix B). Additionally, we indicate the 2σ isophote for the extended Lyα emission discovered around the quasar pair (Fig. 3).

In the text

Fig. 3. Lyα emission around the quasar pair in a field of view of about 200 kpc × 200 kpc (or 26″ × 26″). Left: “optimally extracted” Lyα surface brightness map obtained after subtraction of the quasars point-spread-function (PSF) and continuum in the MUSE datacube (details in Sect. 3.1). To highlight the significance of the detected emission, we indicate the contours for S/N = 3 and 4. This image reveals Lyα bridges extending between the quasar pair. Right: flux-weighted velocity-shift map with respect to the systemic redshift of QSO2 obtained from the first order moment of the flux distribution. A velocity gradient between QSO1 and the portion of the nebula southern than QSO2 is evident (Fig. 4). In both panels we indicate the position of the quasars QSO1 and QSO2 prior to PSF subtraction (white circles), and the masked interloper galaxy G, tentatively at z = 0.457 ± 0.001 (more details in Appendix B). Also, to guide the eye, we overlay a grid spaced by 10″ (or 77 kpc). We also highlight the location of a bright knot (white cross) whose SB value is relevant for the discussion in Sect. 5.1, the direction along which we trace the variations in SBLyα in Sect. 5.2.4, and the seeing circle for these observations (bottom left corner).

In the text

Fig. 4. Left: “optimally extracted” Lyα surface brightness map as in Fig. 3 with the overlaid pseudoslits used to extract the spectra shown in the central and right panels. We assign an ID (blue) to each box of the pseudoslits. Center: normalized spectra of the Lyα emission along the pseudoslit shown in the left panel with solid lines. Each spectrum is color-coded following the color of its box on the left (details in Sect. 4). The dashed (dotted-dashed) vertical lines show the systemic (peak of the Lyα) redshifts for QSO1 (blue) and QSO2 (magenta). The respective shaded regions indicate the errors on the redshifts, as estimated by SDSS. The velocity shifts Δv are computed with respect to the systemic redshift of QSO2. We indicate with vertical arrows the flux-weighted centroid for each spectrum. The spectrum for box 4 (with no clear emission) is shown in Appendix D (Fig. D.1). Right: same as for the central panel, but for the second pseudoslit shown in the left panel with dotted lines.

In the text

	Fig. 5. Flux-weighted velocity-dispersion map obtained as the second order moment of the flux distribution within the 3D mask described in Sect. 3.2. The map shows the same field-of-view and uses the same symbols and nomenclature as in Fig. 3.
In the text

Fig. 6. Profiles of the absorbers along the QSO1 sight-line at the H I Lyα, N V, Si II, C II, and C IV lines. The black histograms show the continuum normalized data, while the orange lines are the sum of all the Gaussian components of the best fit. The locations of the absorbers considered in this analysis are highlighted with vertical dashed lines, that is, ABS1 (magenta), ABS2 (lime), and ABS3 (gray). The zero velocity is set to the redshift of the C IV strongest component. Table 2 gives all the fit parameters. For ABS1, we exclude the best fit solution for NHI (orange) following physical arguments (Sect. 6). We thus show the two extreme alternative fits with NHI and b values that are favored by the current data, that is, b = 100 km s−1, log(NHI/cm−2) = 17 as a dashed-dotted blue line, and b = 200 km s−1, log(NHI/cm−2) = 15 as a dotted green curve.

In the text

Fig. 7. Profiles of the absorbers along the QSO2 sight-line at the H I Lyα line, N V, Si II, C II, and C IV lines. The black histograms show the continuum normalized data, while the orange lines are the sum of all the Gaussian components of the best fit. The location of the absorbers considered are highlighted with vertical dashed lines, that is, ABS1 (magenta; not present along this sight-line), ABS2 (lime), and ABS3 (gray). The zero velocity is set to the redshift of the C IV strongest component of ABS1 along the QSO1 sight-line, as in Fig. 6. Table 2 gives all the relevant fit parameters.

In the text

Fig. 8. Spectral energy distribution (SED) of QSO1 (blue) and QSO2 (orange), used as incident radiation in the Cloudy calculations. We compare the models with the available MUSE data (lighter color for each quasar). Left panel: the vertical lines indicate the energies used to define the different power-laws (Sect. 5.2.1). Right panel: zoomed in version of the box highlighted in the left panel. The emission from the wavelength ranges around the C IV and Lyα line locations are masked to prevent contributions from resonantly scattered photons, as done in Arrigoni Battaia et al. (2015b).

In the text

Fig. 9. Cloudy predictions for plane parallel slabs with log(NH/cm−2) = 20.5, photoionized by a single quasar with the characteristic SED of QSO1. The slabs – characterized by uniform nH – are placed at increasing distance from the quasar. Top: column density of H I as a function of distance from the quasar. The horizontal dotted line indicates the threshold between the optically thin and thick regimes. For nH = 0.1 cm−3 the quasar is able to keep the gas ionized up to ∼200 kpc. Bottom: predicted SBLyα as a function of distance from the quasar. The dotted lines indicate Cloudy models which take into account the presence of the z = 3 UVB by Haardt & Madau (2012). The gray shaded region shows the range of observed SBLyα.

In the text

Fig. 10. Cloudy predictions for plane parallel slabs with total log(NH/cm−2) = 20.5 and nH = 0.5 cm−3, illuminated by the quasar pair QSO1 and QSO2 placed at a separation equal to their observed projected distance (89 kpc). Top left: column density of H I, NH I, as a function of distance. For each model the data-points are color-coded by their respective ionization parameter. In this scenario all the models are optically thin to the ionizing radiation, that is, NH I < 1017.2 cm−2. Top right: comparison of the observed (black line with shaded 1σ error) and predicted SBLyα (blue line). The brown and dashed lines indicate the contribution due to scattering of Lyα photons from the quasars, as explained in Sect. 5.2.2. The red dotted line indicates the total SBLyα summing up the Cloudy prediction and the scattering contribution. The thin black dotted lines are the observed SBLyα along the directions NE for QSO1 and SE for QSO2, at angles 52 and 142° east from north, respectively (details in Sect. 5.2.4). Bottom left: comparison of the observed (black line is the 2σ upper limit) and the predicted (blue line) He II/Lyα ratio as function of the distance from the quasars. The green shaded area represents the parameter space allowed by the observations. Bottom right: comparison of the observed (black line is the 2σ upper limit) and the predicted (blue line) C IV/Lyα ratio as function of the distance from the quasars. The green shaded area represents the parameter space allowed by the observations. The vertical dotted lines in each panel indicate the position of the two quasars, while the striped yellow regions represent the zones used to normalize the quasar PSF, characterized by large uncertainties and, therefore, not considered in the analysis. The dotted red line represents the ratio corrected for the presence of Lyα scattering.

In the text

Fig. 11. Cloudy prediction for plane parallel slabs with total log(NH/cm−2) = 20.5 and nH = 0.5 cm−3, illuminated by the quasar pair QSO1 and QSO2, placed at a separation of 2.9 Mpc, as derived from their systemic redshifts. Top left: column density of H I, NH I, as a function of distance. The data-points are color-coded by the ionization parameter. In this scenario most of the models are optically thick to the ionizing radiation, NH I ≫ 1017.2 cm−2, with only the CGM regions of the two quasars characterized by optically thin gas. Top right: predicted SBLyα for the Cloudy models without (solid blue line) and with (dotted blue line) the UVB. The red dotted line indicates the total SBLyα, summing up the Cloudy prediction and the scattering contribution estimated following Sect. 5.2.2. The gray shaded region shows the range of observed SBLyα. Bottom left: predicted He II/Lyα ratio as function of the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. Bottom right: predicted C IV/Lyα ratio as function the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. In each of the bottom panels, the dotted red lines represent the ratios corrected for the presence of Lyα scattering as modeled in Sect. 5.2.2, while the green horizontal line indicates the local 2σ upper limit on each ratio (Fig. 10). In all four panels, the vertical dotted lines indicate the position of the two quasars, while the striped yellow regions show the zones used to normalize the quasars PSFs.

In the text

Fig. 12. Cloudy prediction for plane parallel slabs with total log(NH/cm−2) = 20.5 and nH = 0.5 cm−3, illuminated by the quasar pair QSO1 and QSO2, placed at an intermediate separation of 600 kpc. Top left: column density of H I, NH I, as a function of distance, color-coded by ionization parameter. In this scenario, a central region of about 300 kpc is optically thick to the ionizing radiation, NH I ≫ 1017.2 cm−2, while the CGM regions of the two quasars are optically thin. Top right: predicted SBLyα for the Cloudy models without (solid blue line) and with (dotted blue line) the UVB. The red dotted line indicates the total SBLyα, summing up the Cloudy prediction and the scattering contribution estimated following Sect. 5.2.2. The predicted SBLyα is basically flat at a value of SBLyα = 2.5 × 10−18 erg s−1 cm−2 arcsec−2, well within the observed range (gray shaded region). Bottom left: predicted He II/Lyα ratio as function of the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. Bottom right: predicted C IV/Lyα ratio as function of the distance from the quasars without (solid blue line) and with (dotted blue line) the UVB. In each of the bottom panels, the dotted red lines represent the ratios corrected for the presence of Lyα scattering as modeled in Sect. 5.2.2, while the green horizontal line indicates the local 2σ upper limit on each ratio (Fig. 10). The vertical dotted lines in all four panels indicate the position of the two quasars, while the striped yellow regions give the zones used to normalize the quasars PSFs.

In the text

Fig. A.1. Five arcseconds circle cutout of the white-light image of the bright star used as PSF in this work. Left: white-light image without post-processing. A faint source on the right part of the star’s PSF is clearly evident. Right: white-light image after replacing the faint source values with the symmetric portion of the MUSE dataset. After this correction, the PSF is well behaved at any wavelength out to five arcseconds (Fig. A.2 and Appendix A for details).

In the text

Fig. A.2. Normalized profile of the bright star used as PSF in this work. The open green squares show the normalized profile derived from the white-light image shown in the right panel of Fig. A.1, while the small dots show the normalized profiles for the star within the 17 MUSE layers encompassed by the obtained 3D mask for the extended Lyα emission (Sect. 3.2). The red dashed line is the best-fit Moffat profile to the white-light image data (β = 2.5 and FWHM = 1.66″). Given the brightness of the used star, there is very good agreement between the profile obtained from the white-light image and the individual layers.

In the text

Fig. B.1. SDSS i-band image extracted from the final MUSE datacube for the same field-of-view shown in Fig. 3. We indicate the position of QSO1, QSO2 and of a faint galaxy unveiled in projection between the targeted quasar pair. This galaxy is an interloper at a different redshift, tentatively at z = 0.457 ± 0.001 (Appendix B and Fig. B.2). For comparison purposes, we also show the 2σ isophote for the extended Lyα emission (black). To guide the eye we overlay a grid spaced by 10″ (or 77 kpc), as done in Fig. 3.

In the text

Fig. B.2. One-dimensional spectrum (black) of the faint galaxy shown in Fig. B.1, as extracted using a circular aperture with radius 1″ from the MUSE data. The red spectrum indicates the error vector. The vertical dashed blue line indicates the position of the only detected line emission (λ = 5435 Å). This galaxy is not associated to the quasar pair as there are no known strong line emissions at a rest-frame wavelength of ∼1800 Å. If this line is [O II]λ3729, this galaxy would be at z = 0.457 ± 0.001.

In the text

Fig. C.1. Top: pseudo 6.25 Å (5 layers) narrow-band map centered at the central wavelength (4872.7 Å) of the 3D mask of the Lyα emission obtained in Sect. 3.2. The map, obtained after PSF and continuum subtraction, shows a 57″ × 57″ (or 438 kpc × 438 kpc) FoV and it is color coded following the Lyα surface brightness. Bottom: χsmth map for the same wavelength range as in the top panel, and obtained using a Gaussian kernel with FWHM = 1.66″ (i.e. similar to the seeing), as explained in Appendix C. To guide the eye, in both panels we overlay a grid spaced by 15″ (or 115 kpc) and we indicate the position of QSO1 and QSO2 prior to their PSF subtraction. In both panels, the interloper galaxy G is indicated with its contour.

In the text

	Fig. D.1. Spectrum of box 4 along the pseudoslit used in the analysis of Fig. 4. The dashed (dotted-dashed) vertical lines show the systemic (peak of the Lyα) redshifts for QSO1 (blue) and QSO2 (magenta). The respective shaded regions indicate the error on the redshift as estimated by SDSS.
In the text

Fig. E.1. Top: χsmth map of 57″ × 57″ (or 438 kpc × 438 kpc) FoV for 6.25 Å (5 layers) centered at the wavelength expected for the C IV line emission (6208.8 Å) given the center of the 3D mask for the Lyα emission. The map is obtained after PSF and continuum subtraction using a Gaussian kernel with FWHM = 1.66″ (i.e., similar to the seeing) as explained in Appendix C. Bottom: same as the top panel, but centered at the wavelength expected for the He II line emission (6573.5 Å). Both maps, covering the corresponding velocity range of Fig. C.1, do not reveal the presence of any clear detection of extended C IV or He II line emission associated with the extended Lyα emission. In both panels, we overlay a grid spaced by 15″ (or 115 kpc) and we indicate with circles the position of QSO1 and QSO2 prior to their PSF subtraction. In both panels, the interloper galaxy G is indicated with its contour.

In the text

Fig. F.1. Predictions of photoionization models for the absorber ABS1 with Z = Z⊙, in the case it is illuminated by QSO1 and the UVB (see Appendix F for details on the models assumptions). The predicted column densities for H I, C IV, N V, C II, and Si II are plotted as a function of distance from QSO1. The horizontal green lines with arrows indicate the observational limits for the same metal ions, while the blue shaded regions show the observational limits for H I. The green hatched boxes indicate the regions where the models matched the observations. The model curves are color-coded following the predicted temperature. The gray dashed lines represent the curves for N V for Z = 0.5 Z⊙. The models in agreement with the observations are characterized by 4.1 ≲ log(T/K) ≲ 4.4 and −1.7 ≲ log U ≲ −0.7. To match the N V absorption, the metallicity should be Z >  0.3 Z⊙.

In the text

Fig. F.2. Predictions of photoionization models for the absorber ABS1 with Z = Z⊙, in the case it is illuminated by QSO1, QSO2 (placed at a projected distance of 89 kpc) and the UVB (details on the models assumptions in Appendix F). All symbols and colors are explained in the caption of Fig. F.1. Similarly to Fig. F.1, the models in agreement with the observations are characterized by 4.1 ≲ log(T/K) ≲ 4.4, −1.7 ≲ log U ≲ −0.7, and Z >  0.3 Z⊙.

In the text

Fig. F.3. Predictions of photoionization models for the absorber ABS1 with Z = Z⊙, in the case where it is illuminated by QSO1, QSO2 (placed at 600 kpc from QSO1) and the UVB (details on the models assumptions in Appendix F). All symbols and colors are explained in the caption of Fig. F.1. Similarly to Figs. F.1 and F.3, the models in agreement with the observations are characterized by 4.1 ≲ log(T/K) ≲ 4.4 and −1.7 ≲ log U ≲ −0.7 (in Appendix F we report the detailed ranges for each plot). The metallicity is constrained to be Z >  0.3 Z⊙.

In the text