© The Authors 2024

1. Introduction

Thanks to their tremendous luminosity, quasars have been used since their discovery in the 1960s as powerful cosmic probes. Their spectra carry a wealth of absorption signatures from the various intervening gaseous environments encountered along the line of sight to the observer, ranging from the diffuse ionized intergalactic medium (e.g., Rauch 1998) to large column densities of neutral gas that imprint characteristic damped Lyman-alpha systems (DLAs; with N(H I) ≥ 2 × 10²⁰ cm⁻²; see Wolfe et al. 2005).

Associated systems are believed to originate within the quasar environment. Due to their physical proximity, a low velocity difference between the quasar and the associated system – dominated by the kinematics of the gas rather than the Hubble flow – is expected. This motivated the first systematic studies of associated systems through exploring samples of proximate DLAs (PDLAs), defined as systems with velocity differences |Δv|< 5000 km s⁻¹ (see Ellison et al. 2010). The detectability of such systems in H I depends on the competition between clustering in the quasar environment and ionization by its intense radiation field of the later (Prochaska et al. 2008).

More recently, Noterdaeme et al. (2019) uncovered a population of proximate molecular absorption systems. Because the presence of H₂ requires a shielding H I layer (see, e.g., Sternberg et al. 2014), proximate H₂ absorbers are expected to be PDLAs. The incidence rate of proximate H₂ absorbers was found to be an order of magnitude higher relative to the intervening H₂ systems, indicating that proximate H₂ systems are indeed tracing the environments surrounding the quasars. The velocity of the foreground absorber relative to the background quasar is likely due to the kinematics of the gas, where positive and negative velocity differences suggest accretion and outflow processes, respectively, in line with an observed trend between this velocity and the gas metallicity (Noterdaeme et al. 2023). Detailed analyses of optical and millimeter data would thus permit further investigations of the outflowing (Noterdaeme et al. 2021) or merging (Balashev et al., in prep.) processes in these systems.

The above observations have not yet been connected to the quasar emission properties. Integral field spectroscopy (IFS) has proved to be invaluable in providing 3D maps of the emission in quasar fields, even reaching a striking 100% detection rate of extended Lyα emission around bright z ≈ 3 quasars, with scales reaching up to hundreds of kiloparsecs (Borisova et al. 2016).

However, the outshining continuum and Lyα emission from the central engine, even if expected to arise from a tiny region, poses a severe challenge to investigations of structures at galactic scales (∼10 kpc), closer to the quasar itself, where outflows, inflows, and galaxy interactions are expected to play a crucial role. The ability to recover the resolved emission depends on the efficiency of the quasar point spread function (PSF) subtraction methods used. Empirical PSF methods generally perform well at large distances from the central quasar but leave large residuals in the inner regions (see examples in Arrigoni Battaia et al. 2019; Farina et al. 2019). Other methodologies involve modeling the quasar spectra and then subtracting the model for each spaxel of the data cube (North et al. 2017), but the accuracy of the extracted emission strongly depends on the modeling, and biases by involving the use of quasar template spectra.

In this Letter we present the first IFS observations of a proximate H₂ system, in the quasar SDSS J125917.31+030922.5 (hereafter J1259+0309). A natural coronagraph is provided by the corresponding high-column density PDLA, enabling Lyα nebula extraction without the necessity of modeling and subtracting the quasar light, and making it the most precise mapping of nebular emission at low projected distances to date. This allows us to trace the gas structure in a complex, likely interacting environment.

2. Observations

Integral field spectroscopy observations of J1259+0309 were carried out by the Multi Unit Spectroscopic Explorer (MUSE) instrument (Bacon et al. 2010), mounted on UT4 at the Very Large Telescope (VLT). Observations were executed under good seeing conditions (average 0.7″). Six 10-minutes-long exposures were taken with different position angles (0°, 22.5°, 45°, 90°, 112.5° and 135°) for better spatial sampling.

Raw exposures were reduced using ESO MUSE Data Reduction Software (pipeline version 2.8.9; Weilbacher et al. 2020). We ran two iterations of the scipost recipe, the first to obtain data images for producing a sky mask to be used in the second iteration. Finally, the Zurich Atmosphere Purge (Soto et al. 2016) software was used to remove any residual pattern from the sky subtraction procedure. To account for correlated noise, the final variance cube was scaled to match the measured variance of the background. The PSF’s full width at half maximum (FWHM) is measured to be 0.7″ in the processed datacube.

In this study we also used spectroscopic data from VLT/X-shooter and adopted the systemic redshift of the quasar host galaxy as z_syst = 3.2405 ± 0.0011 based on CO(3−2) emission from Northern Extended Millimetre Array observations (Noterdaeme et al. 2023). We use z_syst as the reference for zero velocity in the following sections. This z_syst implies a spatial scale of 7.5 pkpc per arcsec for a flat Λ cold dark matter cosmology with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹.

3. Methodology and results

3.1. The PDLA as perfect natural coronagraph

Observing Lyα emission in quasar fields normally requires precise modeling of the quasar PSF as it heavily dominates the emission. For IFS data, this is generally done empirically by estimating the PSF from broadband images extracted from the same data cube (see Borisova et al. 2016 or Farina et al. 2019 for examples). However, the outshining quasar emission itself inevitably leaves residuals that prevent us from measuring faint nebular Lyα emission in the inner regions, impeding any characterization of the emission at the scale of the quasar host or nearby companion galaxies. For J1259+0309, the presence of a PDLA (z_abs = 3.2461 ± 0.0001) with a high H I column density (log N(H/I) ≈ 21.4) efficiently suppresses the quasar light over a wide wavelength range (optical depth higher than 5 over ≈40 Å), so, for the first time, there is no need to perform PSF nor continuum subtraction, and we can directly extract the extended Lyα emission, including the relatively narrow component at the core of the DLA profile (see Fig. 1). In short, the PDLA acts as a natural coronagraph that suppresses the point source emission well before any dispersive process due to the Earth’s atmosphere, the telescope, or the instrument takes place.

Fig. 1. J1259+0309 spectrum and MUSE narrow band images around the DLA core. Top row: J1259+0309 X-shooter 1D combined spectrum around systemic the Lyα wavelength. The solid red line shows the reconstructed quasar emission, to which the PDLA Voigt profile model (log N(H/I) ≈ 21.4) is then applied to obtain the solid purple line. Vertical dashed red and purple lines indicate the expected location of lines at the quasar systemic redshift from CO(3−2) emission (also the zero point in the velocity from the top axis) and the PDLA, respectively. The inner subpanels present a zoomed-in view in around wavelength regions containing the Lyα emission in the PDLA core and a low-ionization absorption line. Bottom row: MUSE narrowband images from the wavelength regions A and B indicated in the top panel by the vertical blue and gray stripes, respectively. The third panel shows a quasar-continuum-only region at 5400−5412 Å. All images were smoothed using a Gaussian kernel, and have the same scale, and coordinates are relative to the quasar central position in the sky. In panel A, horizontal dashed rectangles illustrate the position the two X-shooter slits employed to obtain the combined 1D spectrum. A scale of 10 proper kpc (≈1.3″) is shown as a horizontal solid black line.

3.2. Identification and subtraction of a low-redshift galaxy

The left panel of Fig. 2 shows a white-light image of the data cube, in which several sources are detected. Within this 16″ × 16″ region, we identified only one additional line-emitting source, which turns out to be a z = 0.187 galaxy whose Hγ emission, falling at 5153 Å, could in principle contaminate the high-z Lyα map at this position. We predicted the Hγ emission based on the strength of the Hα and Hβ lines, assuming recombination and the values reported in Table B.7 from Dopita & Sutherland (2003), and subtracted the emission from the cube on a spaxel-by-spaxel basis wherever we managed to detect Hα with S/N ≥ 3.

Fig. 2. Spectra of the southern extension of the nebula and a close low-redshift interloper. Left-most panel: White-light image of the MUSE data cube around the quasar (at (0, 0) position) prior to subtracting contamination sources. 5 and 10σ contours of the extended Lyα nebula are overlaid. Enclosed by a dashed elliptical contour is the zgal = 0.187 line-emitting galaxy. Other panels: Spectra of the southern portion of the nebula (top row) and the line-emitting galaxy (bottom row) shown in three different wavelength portions covering Hα, Hβ+[O III] and Hγ at zgal, the last of which coincides with Lyα at zquasar = 3.24.

3.3. Optimal extraction of the Lyα nebula

We employed similar methods as outlined in Borisova et al. (2016) and Farina et al. (2019): The only pre-detection process was background subtraction at spaxels where the quasar is not detected. A S/N cube was then created from the results of the last step, which was converted into a binary cube, where unit values indicate voxels with S/N ≥ 2.5 (zero otherwise). Finally, a 3D segmentation mask was obtained by running a friends-of-friends algorithm on this binary cube. The largest group of connected voxels is considered the final mask. An extended explanation of the methodology will be presented in Urbina et al. (in prep.). The optimally extracted Lyα nebula (zeroth moment) and its first and second moments are shown in Fig. 3. In addition to the main emission around the position of the quasar, we detect an elongated secondary emission centered ≈2 arcsec (≈15 kpc) to the east and a weaker southern emission connected through a filamentary structure (S/N ≈ 2). We also present four zoomed-in velocity-binned images of the nebula in the inner regions for a better visualization of the structures at those scales.

Fig. 3. Lyα nebula around J1259+0309. First panel: Optimally extracted Lyα nebula (zeroth moment) from continuum-subtracted MUSE data cube. The spatial projection of the segmentation mask is shown as a thick black contour around the detected emission, and it is approximately a 2.5σ = 2.45 × 10−18 erg s−1 cm−2 arcsec−2 contour. Thin black contours enclose 5, 10 and 15σ regions. The black horizontal line indicates a scale of 10 pkpc. For reference, the PSF is depicted by a black circle with a diameter of 0.7″. Second panel: Velocity map of the emission with the zero point defined by zsys. Third panel: Second moment map to show the equivalent Gaussian FWHM, where the instrumental FWHMins ≈ 170 km s−1 has been subtracted in quadrature. Each image has a linear projected size of 14″ × 17″ and coordinates are relative to the quasar position. Last panel: Zoomed-in views of the inner region of the emission in four ≈290 km s−1 wide velocity bins. The central value is indicated in the top-left corner of each image. Isophotal contours of SB values of 15, 30, and 45 × 10−18 erg s−1 cm−2 arcsec−2 are shown.

3.4. Additional Lyα emitters in the field

To explore the wider environment around J1259+0309, we searched for potential Lyα emitters (LAEs) in the cube using the LSDCat weak line-emitting source detection software (Herenz 2023). LSDCat employs a matched filtering approach, which is applied in the spatial and then in the spectral direction. For J1259+0309, we used a 2D Moffat profile with a FWHM ≈ 0.7″ in the spatial direction and a Gaussian profile with a FWHM of 300 km s⁻¹ for the spectral direction (i.e., a compact source with a narrow line emission). The output of the software is a detection significance cube. In Fig. 4, we show the extracted spectrum and narrowband images of the two LAE candidates detected this way and summarize their properties in Table 1.

Fig. 4. LSDCat LAE candidates. (a.1, a.2): LAE candidate spectra obtained obtained using 1″ apertures centered at the coordinates reported in Table 1. Dashed red and purple lines show the expected wavelength of the emission at the quasar systemic and PDLA redshift. The solid black line shows a Gaussian model fitted to the data. Propagated 1σ errors are shown as vertical blue lines. (b.1, b.2): LSDCat output detection significance maps at the wavelength channel where the peak of Lyα emission is reached.

4. Discussion

In Fig. 5 we present the circularly averaged surface brightness (SB) profile of the Lyα emission around J1259+0309, together with that of SDSS J095253.83+011421.9 (hereafter J0952+0114 Marino et al. 2019) and the best exponential fit to the QSO MUSEUM sample (Arrigoni Battaia et al. 2019). Apart from a few noticeable cases, most Lyα emission in this sample is concentrated within 50 kpc of the quasar and has a relatively isotropic profile with a median minor-to-major axis ratio of 0.73. However, PSF contamination prevent us from constraining the Lyα SB profile within the inner 2 − 3″. This is unfortunate since the close environment of quasars (within a few 10 kpc) is where outflowing or infalling gas and galaxy interactions should have most noticeable effects and hence are key to understanding the baryonic cycle or galaxy interactions. Moreover, precisely determining the SB profile of the inner regions can provide some insights into the powering mechanism of the nebulae. Postprocessing cosmological simulations by Costa et al. (2022) suggest that is not possible to recover an approximately flat profile at low projected distances when only assuming recombination and collisional excitation. As the profile shown in Fig. 5 is indeed approximately flat up to ≈10 kpc, we can conclude that mechanisms such as resonant scattering may be relevant in this case.

Fig. 5. Circularly average Lyα SB profiles (corrected for redshift dimming) as a function of projected distance from the central quasar. Solid orange, blue, and green lines stand for J1259+0309 (this work, with the 1σ error shown as the shaded area), J0952+0114 (Marino et al. 2019), and the average of the QSO MUSEUM sample (Arrigoni Battaia et al. 2019, with 1σ scatter as the light green area), respectively. The dotted green line is an extrapolation of the QSO MUSUEM sample toward low radii. The dashed orange line shows the eastward average profile for J1259+0309 (see Sect. 4). The sharp drop at 30 kpc in the J1259+0309 profile is due to sensitivity limitations. The PSF shape, derived from the quasar continuum, is shown as a solid grey line (with 1σ error) and is scaled for visual purposes.

MUSE observations of J0952+0114, where PSF subtraction is not required due to the presence of a PDLA¹, show a much higher flux in the central regions compared to what is expected from extrapolating the MUSEUM profile but did not reveal much structure either. This particular quasar originally attracted attention because it lacks the broad Lyα component expected from the broad metal line emission (Hall et al. 2004). Jiang et al. (2016) later showed the presence of a PDLA that covers only the continuum and the broad-line region but is filled with relatively narrow emission (∼1000 km s⁻¹), explaining the difficulty in identifying the PDLA as such. This PDLA exhibits strong absorption from silicon in its excited fine-structure state (Si II*), which may indicate a high compactness of the associated gas (Fathivavsari 2020). Additionally, the quasar spectrum shows a higher line emission-to-continuum ratio compared to quasar composite spectra, suggesting differential reddening, where dusty gas in the PDLA does not fully cover the metal emission regions either. Regardless of its exact origin, the extended Lyα emission in J0952+0114 may still be dominated by a central point-like source that is observationally smoothed by the PSF over a few arcseconds. This would explain the very strong Ly-alpha emission in the DLA as well as the smooth symmetric shape centered on the quasar position, with a lack of internal structure at r < 20 kpc (∼3″).

In the case of J1259+0309, the central point source emission is efficiently suppressed over a wide wavelength range, and the remaining emission has a SB closer to what is expected from extrapolating the peripheral emission profile. While we do not know the actual extent of the absorbing gas, our Lyα mapping is not blinded by the central emission, which allowed us to resolve the structures presented in Fig. 3. We indeed observe a clear asymmetric shape of Lyα emission even at small separations of 2″ (≈15 kpc) from the quasar, which extends eastward. For the inner regions of the nebula (without the southern extension), we measure a flux-weighted minor-to-major axis ratio of 0.48. Only ≈10% of Lyα nebulae have a minor-to-major axis ratio lower than J1259+0309 (Arrigoni Battaia et al. 2019).

The relative velocity of the PDLA with respect to z_syst is Δv = 430 ± 80 km/s – the PDLA is at a higher redshift than the quasar – which suggests that absorbing gas must be moving toward the quasar. The gas traced by the DLA could hence be tracing a feeding stream or leading to another galaxy in the group that may eventually merge with the quasar host. This system is also a low-metallicity proximate H₂ absorber (Z ≈ 0.03 Z_⊙, log N (H₂)≈19.10; Noterdaeme et al. 2023) with little to no dust reddening. Therefore, for the H I-H₂ transition to occur under the influence of the quasar radiation, we conclude that there must be a relatively high volume density (Noterdaeme et al. 2019)². We finally note that the quasar is also a broad absorption line (BAL) quasar, which suggests quasar-driven outflows may be present as well. Nevertheless, this is not necessarily related to the PDLA since fewer than 10% of quasars are expected to be BALs (Trump et al. 2006).

These PDLA properties suggest the quasar is located in a dense environment with ongoing galaxy interactions. The Lyα emission shows consistent signatures of this. A positive gradient in the FWHM toward the central quasar (see the third panel in Fig. 3) suggests disturbed kinematics occurring in regions with outflows or infalling material, although the Lyα emission can be further broadened by radiative-transfer effects only (e.g., Laursen et al. 2011). The inner nebula shows enhanced secondary emission to the East, possibly due to a galaxy companion³, which drives the asymmetric shape as well. Additionally, the nebula extends southward from the central quasar at 7″ (∼50 kpc). This region shows a different kinematic signature compared to the inner regions, with a velocity gradient that may indicate a rotating gas (Fig. 3). Albeit detected only with S/N ≈ 2, a filamentary structure may connect this emission to the material in the immediate vicinity of the quasar. We also note that the number of LAE candidates detected within the J1259+0309 field (2 LAE) is ≈4 times larger than expected for similar quasar fields at z ∼ 3 (assuming a 2σ detection limit of 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² and after completeness corrections; Arrigoni Battaia et al. 2019), suggesting a denser-than-average environment.

The rich environment unveiled by the Lyα emission around J1259+0309 suggests that galaxy interactions are actively occurring, on scales of ∼15 kpc up to ∼50 kpc. Such a dense environment is not typical for quasars at these redshifts. This suggest that, by selecting the system based on H₂ absorption, we have selected a richer-than-average environment. As large programs surveying the typical z ∼ 3 population of quasars are already available, the next step is to explore a large number of quasar fields with proximate absorbers (Urbina, in prep.) to assess whether selecting quasars with proximate absorbers dictates their environmental properties or whether the properties seen in J1259+0309 field can be attributed to the cosmic variance.

Acknowledgments

F.U. acknowledges support by Subdirection of Human Capital ANID (national MSc 2023/22231861). S.B. is supported by RSF grant 23-12-00166. S.L. acknowledges support by FONDECYT grant 1231187. Authors thank the referee for all comments and suggestions provided.

References

All Tables

All Figures

Fig. 1. J1259+0309 spectrum and MUSE narrow band images around the DLA core. Top row: J1259+0309 X-shooter 1D combined spectrum around systemic the Lyα wavelength. The solid red line shows the reconstructed quasar emission, to which the PDLA Voigt profile model (log N(H/I) ≈ 21.4) is then applied to obtain the solid purple line. Vertical dashed red and purple lines indicate the expected location of lines at the quasar systemic redshift from CO(3−2) emission (also the zero point in the velocity from the top axis) and the PDLA, respectively. The inner subpanels present a zoomed-in view in around wavelength regions containing the Lyα emission in the PDLA core and a low-ionization absorption line. Bottom row: MUSE narrowband images from the wavelength regions A and B indicated in the top panel by the vertical blue and gray stripes, respectively. The third panel shows a quasar-continuum-only region at 5400−5412 Å. All images were smoothed using a Gaussian kernel, and have the same scale, and coordinates are relative to the quasar central position in the sky. In panel A, horizontal dashed rectangles illustrate the position the two X-shooter slits employed to obtain the combined 1D spectrum. A scale of 10 proper kpc (≈1.3″) is shown as a horizontal solid black line.

In the text

Fig. 2. Spectra of the southern extension of the nebula and a close low-redshift interloper. Left-most panel: White-light image of the MUSE data cube around the quasar (at (0, 0) position) prior to subtracting contamination sources. 5 and 10σ contours of the extended Lyα nebula are overlaid. Enclosed by a dashed elliptical contour is the zgal = 0.187 line-emitting galaxy. Other panels: Spectra of the southern portion of the nebula (top row) and the line-emitting galaxy (bottom row) shown in three different wavelength portions covering Hα, Hβ+[O III] and Hγ at zgal, the last of which coincides with Lyα at zquasar = 3.24.

In the text

Fig. 3. Lyα nebula around J1259+0309. First panel: Optimally extracted Lyα nebula (zeroth moment) from continuum-subtracted MUSE data cube. The spatial projection of the segmentation mask is shown as a thick black contour around the detected emission, and it is approximately a 2.5σ = 2.45 × 10−18 erg s−1 cm−2 arcsec−2 contour. Thin black contours enclose 5, 10 and 15σ regions. The black horizontal line indicates a scale of 10 pkpc. For reference, the PSF is depicted by a black circle with a diameter of 0.7″. Second panel: Velocity map of the emission with the zero point defined by zsys. Third panel: Second moment map to show the equivalent Gaussian FWHM, where the instrumental FWHMins ≈ 170 km s−1 has been subtracted in quadrature. Each image has a linear projected size of 14″ × 17″ and coordinates are relative to the quasar position. Last panel: Zoomed-in views of the inner region of the emission in four ≈290 km s−1 wide velocity bins. The central value is indicated in the top-left corner of each image. Isophotal contours of SB values of 15, 30, and 45 × 10−18 erg s−1 cm−2 arcsec−2 are shown.

In the text

Fig. 4. LSDCat LAE candidates. (a.1, a.2): LAE candidate spectra obtained obtained using 1″ apertures centered at the coordinates reported in Table 1. Dashed red and purple lines show the expected wavelength of the emission at the quasar systemic and PDLA redshift. The solid black line shows a Gaussian model fitted to the data. Propagated 1σ errors are shown as vertical blue lines. (b.1, b.2): LSDCat output detection significance maps at the wavelength channel where the peak of Lyα emission is reached.

In the text

Fig. 5. Circularly average Lyα SB profiles (corrected for redshift dimming) as a function of projected distance from the central quasar. Solid orange, blue, and green lines stand for J1259+0309 (this work, with the 1σ error shown as the shaded area), J0952+0114 (Marino et al. 2019), and the average of the QSO MUSEUM sample (Arrigoni Battaia et al. 2019, with 1σ scatter as the light green area), respectively. The dotted green line is an extrapolation of the QSO MUSUEM sample toward low radii. The dashed orange line shows the eastward average profile for J1259+0309 (see Sect. 4). The sharp drop at 30 kpc in the J1259+0309 profile is due to sensitivity limitations. The PSF shape, derived from the quasar continuum, is shown as a solid grey line (with 1σ error) and is scaled for visual purposes.

In the text