Down-the-barrel observations of a multi-phase quasar outﬂow at high redshift VLT/X-shooter spectroscopy of the proximate molecular absorber at z = 2.631 towards

We present ultraviolet to near infrared spectroscopic observations of the quasar SDSSJ001514 + 184212 and its proximate molecular absorber at z = 2 . 631. The [O iii ] emission line of the quasar is composed of a broad ( FWHM ∼ 1600kms − 1 ), spatially unresolved component, blueshifted by about 600kms − 1 from a narrow, spatially-resolved component ( FWHM ∼ 650kms − 1 ). The wide, blueshifted, unresolved component is consistent with the presence of outﬂowing gas in the nuclear region. The narrow component can be further decomposed into a blue and a red blob with a velocity width of several hundred kms − 1 each, seen ∼ 5pkpc on opposite spatial locations from the nuclear continuum emission, indicating outﬂows on galactic scales. The presence of ionised gas on kpc scales is also seen from a weak C iv emission component, detected in the trough of a saturated C iv absorption that removes the strong nuclear emission from the quasar. Towards the nuclear emission, we observe absorption lines from atomic species in various ionisation and excitation stages and conﬁrm the presence of strong H 2 lines originally detected in the SDSS spectrum. The overall absorption proﬁle is very wide, spread over ∼ 600kms − 1 , and it roughly matches the velocities of the narrow blue [O iii ] blob. From a detailed investigation of the chemical and physical conditions in the absorbing gas, we infer densities of about n H ∼ 10 4 − 10 5 cm − 3 in the cold ( T ∼ 100K) H 2 -bearing gas, which we ﬁnd to be located at ∼ 10kpc distances from the central UV source. We conjecture that we are witnessing di ﬀ erent manifestations of a same AGN-driven multi-phase outﬂow, where approaching gas is intercepted by the line of sight to the nucleus. We corroborate this picture by modelling the scattering of Ly- α photons from the central source through the outﬂowing gas, reproducing the peculiar Ly- α absorption-emission proﬁle, with a damped Ly- α absorption in which red-peaked, spatially o ﬀ set, and extended Ly- α emission is seen. Our observations open up a new way to investigate quasar outﬂows at high redshift and shed light on the complex issue of AGN feedback.


Introduction
Feedback from active galactic nuclei (AGN) is an essential element in modern models of galaxy formation and evolution (e.g., Silk & Rees 1998). It may quench star formation (e.g., Zubovas & King 2012;Pontzen et al. 2017;Terrazas et al. 2020), impact galaxy morphology , regulate the growth of the supermassive black holes (Volonteri et al. 2016), and so on. It may also have positive feedback on star formation through compression of the gas (e.g., Zubovas et al. 2013;Richings & Faucher-Giguère 2018). The most supportive evidence for AGN feedback is the presence of kiloparsecscale outflows, which result from the propagation of energy and momentum from accretion disc winds to the host galaxy interstellar medium (e.g., Fabian 2012; Costa et al. 2015). At low and intermediate redshifts, such outflows have been observed in ionised and molecular phases, with many studies focusing on Type 2 quasars, where the obscuration of the nuclear emission facilitates the observations (e.g., Sun et al. 2017). Outflows have also been observed in unobscured, luminous AGN (Type 1 quasars), but mostly at low redshift, such as the well-known Mrk 231 at z = 0.04 (Rupke & Veilleux 2011;Feruglio et al. 2015). Naturally, observational constraints on hot and cold outflows in high-redshift Type 1 quasars are much harder to obtain, not only because of the dimming of the light, but also because the coarser angular resolution impedes observations close to the bright nucleus.
In turn, the cosmological distances and high luminosity of the nuclear emission of high-z quasars make them excellent background sources against which gas in any phase can be studied in absorption. This includes not only gas at any place along A&A 646, A108 (2021)  Notes. (a) Best-fitting relative abundance compared to the default template in molecfit: N(H 2 O) = 2.725 × 10 4 ppmv and N(CO 2 ) = 368.5 ppmv. Typical uncertainties of ∼0.01. (b) Best-fit spectral resolution of the NIR spectrum from molecfit with typical uncertainties of 0.1 pixel.
the line of sight (from the intergalactic medium, giving rise to the Ly-α forest or closer to intervening galaxies, producing damped Ly-α systems (DLAs) in the quasar spectrum), but also gas close to the supermassive black hole, in the quasar host galaxy or in galaxies belonging to the quasar group environment. While broad absorption lines have long been associated with gas flows close to the central engine (e.g., Weymann et al. 1991;Arav et al. 2018), large-scale outflows, when intercepted by the line of sight to the nucleus, should also produce a detectable signature in absorption. Absorption spectroscopy has proven to provide very detailed information about the kinematics, chemical composition, and physical conditions in the absorbing gas, in particular when molecules are detected. Indeed, the formation, survival, and excitation of molecules is very sensitive to the prevailing physical conditions, such as abundance of dust, temperature, density, and ambient UV field (e.g., Reimers et al. 2003;Srianand et al. 2005;Cui et al. 2005;Guillard et al. 2009;Balashev et al. 2017;Wakelam et al. 2017, among many others).
Combining information from both emission and absorption line measurements could hence provide us with a fresh view of quasar outflows and feedback at high redshift. We recently embarked on a search for damped molecular hydrogen (H 2 ) absorption lines at the quasar redshift in low-resolution spectra from the Sloan Digital Sky Survey (SDSS) III (BOSS) Data Release 14. Through this targeted search, we discovered a population of strong proximate H 2 absorption systems at z > 2.5 with log N(H 2 ) > 19 . Detailed follow-up observations are currently on-going using the multi-wavelength spectrograph X-shooter on the Very Large Telescope (VLT).
In this paper, we report on the observations of the first object in our sample, the quasar SDSS J001514.82+184212.34 (hereafter J0015+1842) at z = 2.63. We interpret the observations as evidence of kilo-parsec scale, multi-phase, biconical outflows that are mainly orientated along the line of sight and also probed in absorption against the nucleus. We present the observations and data reduction in Sect. 2, the analysis of spatiallyresolved emission lines in Sect. 3, and that of absorption lines and dust extinction towards the nucleus in Sect. 4. We discuss our results in Sect. 5, along with the derivation of the chemical and physical conditions in the absorbing clouds and the modelling of leaking Ly-α emission. We summarise our findings in Sect. 6. Throughout this paper, we assume a flat ΛCDM cosmology with H 0 = 68 km s −1 Mpc −1 , Ω Λ = 0.69, and Ω m = 0.31 (Planck Collaboration XIII 2016).

Observations and data reduction
J0015+1842 was observed in stare mode with the X-shooter spectrograph at the VLT at the Paranal Observatory in Chile  between July and September 2019. The log of observations is provided in Table 1.
The X-shooter spectrograph simultaneously covers the wavelength range from 0.3 to 2.5 µm in three separate spectrographs, which are the so-called arms: UVB ranging from 0.3 to 0.6 µm; VIS ranging from 0.6 to 1.0 µm; and NIR ranging from 1.0 to 2.5 µm. For all observations, slit widths of 1.0, 0.9, and 1.2 arcsec were used for the UVB, VIS, and NIR arms, respectively. The slit was aligned to the parallactic angle at the start to minimise slit losses of the exposure, and it was kept fixed on sky during the ∼1 h long exposure (see Fig. 1). The observation from August 30 failed due to an error in the atmospheric dispersion corrector (ADC), and exposures from this date were therefore not included in the analysis. For one other observation (July 29), we noted a strong chromatic suppression of flux in the UVB arm 1 , hence the UVB spectrum from this observation has been discarded from our analysis. As the loss of flux in the VIS and NIR arms is much less severe, we used these spectra in our analysis but weighted them appropriately (using the inverse variance) when combining the individual spectra.
The raw spectra were processed using the official esorex pipeline for X-shooter version 2.6.8. Before passing the raw spectra through the pipeline, we corrected cosmic ray hits using the code Astroscrappy (McCully et al. 2018), which is a Python implementation of the algorithm by van Dokkum (2001). The pipeline then performed subtraction of bias and dark levels on the CCD, spectral flat fielding, tracing of the curved echelle orders, and computation of the wavelength solution. The individual curved 2D orders were then rectified onto a straightened grid. The sky background was subsequently subtracted before the individual orders were combined. The full spectrum was flux-calibrated using a sensitivity function calculated for each observation using a spectroscopic standard star observed during the same night. The three spectra (UVB, VIS, NIR) were stitched together in the overlapping regions, correcting for any small offsets between the absolute flux calibration of the individual arms. Lastly, in order to correct for possible slit loss, we scaled the calibrated spectra to match the SDSS i band photometry. Because of the possible variation of the quasar brightness between the time of observation for the SDSS photometry and our data, the absolute flux calibration has a systematic uncertainty of order 20%, which corresponds to the typical long-term variability of quasars (Hook et al. 1994).
From the 2D spectra produced by the pipeline, we performed an optimal extraction (Horne 1986) to obtain 1D spectra. The 1D and 2D spectra were subsequently corrected to vacuum wavelengths and shifted to the heliocentric reference frame before being corrected for Galactic extinction using the maps by Schlafly & Finkbeiner (2011). Lastly, we combined the individual observations for each arm into final 1D spectra using a weighted mean.
Telluric correction. The NIR spectra were corrected for telluric absorption using the software tool molecfit Kausch et al. 2015). The code fits a synthetic absorption profile to the telluric absorption features in the data. For the NIR spectra, the dominant molecular species giving strong telluric absorption lines are H 2 O and CO 2 . The relative abundances of these two species were fitted together with a fifth-order polynomial model for the intrinsic spectrum in order to reproduce the shape of the wings of the quasar emission lines. We fitted the telluric model using seven narrow bands of the spectra dominated by telluric absorption features only. These seven regions were chosen to encompass regions of interest, that is, the regions around the quasar emission lines from oxygen. The wavelength intervals used were as follows: 1.13-1.14, 1.37-1.39, 1.44-1.45, 1.75-1.76, 1.82-1.83, 1.94-1.97, and 2.02-2.03 µm. For each region, we iteratively masked out pixels affected by skyline residuals or bad pixels. The optimal parameters for each spectrum are given in Table 1. The optimised telluric absorption model for each spectrum was then used to correct the corresponding 1D and 2D spectra.

Emission line analysis
We analysed the 1D and 2D together spectra to extract information about the spatial and spectral properties of the various emission lines of interest. In addition to the 1D spectrum obtained through the optimal extraction method, which mostly samples the unresolved nuclear emission (due to the weighting profile), we also constructed a total 1D NIR spectrum directly extracted from the combined 2D data over a wider spatial range of ±9 pixels (or equivalently 1.9 arcsec) around the quasar trace without weighting. This is specifically constructed for the analysis of spatially extended emission lines such as [O iii]. In summary, we therefore consider three kinds of 1D spectra in the following: (i) the on-trace spectrum, which is obtained from the trace by optimal extraction; (ii) the overall spectrum, which is extracted within wide spatial range; and (iii) the off-trace spectrum, which is the difference between the two others.

[O ii]
The [O ii]λλ3727, 3729 doublet is redshifted into a strong telluric band and not directly visible in the original data. However, applying the telluric correction described in Sect. 2 allowed us to recover the lines, which, given their large velocity width, are strongly blended and appear as a single line in Fig. 2. To estimate the intrinsic line width we fitted the [O ii] doublet using the sum of two Gaussians plus a linear continuum in that region. The relative positions of the Gaussians were kept fixed to the doublet wavelength ratio, and their widths were forced to be the same. In principle, the amplitude ratio for the [O ii] doublet can vary between 0.35 and 1.5, depending on the kinetic temperature of the gas and the electron density (e.g., Kisielius et al. 2009). However, because this region of the spectrum is quite noisy, we assumed a fixed ratio of 1.3. This assumption has almost no influence on the total line flux and width. The obtained parameters are given in Table 2 Sect. 3.4), we considered the overall 1D spectrum instead of the on-trace spectrum when fitting this line. We also used the off-trace spectrum to obtain an initial guess for the parameters of the narrow component free from contamination by the broad A108, page 3 of 19 A&A 646, A108 (2021) Notes. The apparent peak and FWHM of the leaking Ly-α emission were estimated visually from the non-Gaussian spectral profile and the luminosity from direct integration. For the other lines, these parameters were obtained from fitting the Gaussian profiles as described in the text. Statistical uncertainties on the line luminosities are about 10% for NIR lines and 5% for UVB lines. However, the uncertainty on the absolute flux calibration is about 20%. This has no effect on the line ratios. spectrum collapsed over its full spatial extent at wavelengths 1.790 < λ(µm) < 1.825 and the optimal extraction elsewhere, for presentation purposes. The green line shows our estimate of featureless continuum emission. Individual Gaussians constrained during the fit are shown in different colours, with the total fitted profile shown in red. nuclear component (top panel of Fig. 3). We then fitted the overall spectral profile relaxing all parameters and rejecting deviating pixels iteratively. The corresponding line parameters are given in Table 2.
The peak of the narrow [O iii] emission provides us with an estimate of the quasar host's systemic redshift, which also matches that obtained from [O ii]. Given the associated uncertainty of about 100 km s −1 , we rounded the value to z = 2.631, which we take as a convenient reference redshift to define the zero point of our velocity scale throughout this work.

C iv emission
The region of 1D spectrum around the C iv emission line is shown in Fig. 4. Strong absorption from the proximate DLA (PDLA) is observed in the red wing of the emission line. We fitted the emission line using the sum of two Gaussians plus a linear continuum in that region, excluding the region affected by absorption. The C iv emission peaks at ∼5611 Å, while the C iv absorption is centred at ∼5625 Å. Furthermore, the C iv emission is shifted by about 700 km s −1 bluewards of the systemic redshift but matches the location of the broad [O iii] component well. Interestingly, while the shape and strength of the C iv absorption profile indicate strong saturation, a residual emission is seen in the trough, with a possible peak at the systemic velocity. Even though the value for each pixel remains relatively close to the noise level, the residual emission is consistently seen above zero for almost all pixels, with integrated flux F resi ∼ (1.3 ± 0.1) × 10 −17 erg s −1 cm −2 , that is, 13σ significance. In principle, we cannot rule out that this is due to the absorption profile being composed of a blend of many narrow unsaturated components. We note, however, that reproducing the sharp edges and the width of the profile while keeping the non-zero absorption in the centre is very difficult, even more so with the condition that the C iv column density be bounded by the metallicity. A more likely explanation is that the C iv emission is not fully covered by the C iv absorbing gas. This is corroborated by the spatial extent of the residual emission, which is discussed in the following sub-section.

Spatial extent of line-emitting regions
The spatial profiles of the used a region centred at the peak position within one FWHM. However, to assess whether the photons from the broad component come from compact (unresolved) or spatially extended regions, we restricted the region to 1.805−1.815 µm. This avoids contamination by the narrow component. In the case of C iv, we collapsed the data over the 5580−5650 Å region, which roughly corresponds to the FWHM centred at the line peak. We also collapsed the data in the 5619−5628 Å region, corresponding to the region of the C iv absorption trough where some residual emission is seen (Fig. 4).
We found that the spatial profiles of [O ii], H β the broad [O iii], and the overall C iv emission are unresolved, perfectly matching the spatial profile of the continuum and hence mostly originating from the nuclear region. In turn, the narrow [O iii] line is significantly extended, over several kpc, on both sides of the trace of the continuum. The C iv absorption trough acts like a natural coronagraph, removing the bright quasar glow at these wavelengths. This provides easier access to extended but fainter C iv emission. While it is hard to appreciate the spatial profile of the residual (weak) C iv emission from the 2D data (top panel of Fig. 4), when collapsing along the wavelength axis, and albeit still noisy, this profile clearly presents a spatial kpc-scale extension below the trace.
In order to further characterise the extension of the narrow We remark that the [O iii] emission extends above and below the trace for roughly negative and positive velocities (according to our reference redshift), respectively.

The Ly-α emission
The spectrum of J0015+1842 in the Ly-α region is the sum of the quasar emission absorbed by the proximate DLA and the leaking Ly-α emission superimposed on the DLA trough, as already mentioned in our discovery work . One difficulty with the low-resolution, low signal-to-noise ratio (S/N) BOSS fibre spectrum was to ascertain the combined model of absorption and emission since the quasar redshift was not easily determined and the wings of the DLA were not clearly detected. The higher quality X-shooter spectrum confirms our previous decomposition (see Fig. 8). The absorption profile is well constrained by the damping wings of the DLA as well as other lines from the Lyman series, which were fitted simultaneously with H 2 lines (Fig. 10). The intrinsic quasar emission profile was reconstructed simultaneously using a smooth spline profile, aided by the quasar template matching (see Sect. 4.5); however, the exact shape of the peak is highly uncertain since it is positioned at the same place as the DLA trough. Nonetheless, we note that the exact shape and strength of the intrinsic, broad Ly-α emission line has little influence on the analysis presented in this paper.
In Fig. 8, we also show the original BOSS spectrum. The residual Ly-α emission appears nearly two times stronger in the BOSS spectrum than in the X-shooter spectrum. The reason is that the optimal extraction method used for the X-shooter spectrum gives lower weight to the extended and spatially offset Ly-α emission. Instead, the 2D spectra clearly show that the residual Ly-α emission is spatially offset and extended ( Fig. 9), and the overall (non-optimal) X-shooter extraction does match the SDSS data (top-right panel of that figure) much better. This will be discussed in the following.
Since the system was identified through its H 2 absorption lines imprinted on the quasar continuum, the absorbing gas intercepts the line of sight to the central engine, removing the unresolved emission and leaving us with a perfect PSF subtraction of the nuclear emission over the saturated region (∼4405−4420 Å A108, page 5 of 19   observed) of the DLA trough. In other words, the DLA acts as a natural coronagraph (see e.g., Finley et al. 2013;Jiang et al. 2016). However, the leaking Ly-α emission remains partly blended with light from the quasar nucleus in the wings of the DLA (mostly in the red wing, λ obs ∼ 4420−4430 Å). We therefore subtracted the nuclear emission over a wider region in order to isolate the extended emission. For this purpose, we constructed a 2D model of the nuclear emission using a Moffat profile to describe the spatial point spread function (SPSF). The shape and position of the SPSF is constrained by the continuum emission far from the Ly-α region. The amplitude of the spatial  profile is then multiplied by the combined 1D spectral model for the nuclear emission and Ly-α absorption (i.e. the red line in Fig. 8). The resulting core-subtracted 2D spectra are shown in the right-hand panels of Fig. 9. The right-most panels indicate the corresponding spatial profiles for the nuclear (red) and residual (black) emission. In the top-right panel of Fig. 9, we show the collapsed 1D spectra of the extended Ly-α emission over the full spatial range presented in the 2D panels (i.e. ∼6 arcsec). The individual spectra for the different position angles are all consistent with one another, indicating that differential slit losses are negligible. Given the observed north-south extension of the emission, the narrow slit width (1 arcsec in UVB), and the position angles of the slits (see Fig. 1), this suggests that most of the extended Ly-α emission falls within the region covered by all three slits, with an offset of ∼0.7 arcsec from the nuclear emission, as seen from the side panels of Fig. 9. Most of the extended Ly-α emission then corresponds to the same location as the red component of the narrow [O iii] emission. However, the Ly-α flux remains slightly higher in the BOSS spectrum (2 arcsec diameter fibre), in particular bluewards of the Ly-α peak. We note that shifts in the zero level have been observed in BOSS spectra, but on much lower levels. Furthermore, other saturated lines in the BOSS spectrum do not show any residual flux. This means that the excess emission at ∼4400−4420 Å in the BOSS spectrum compared to the X-shooter one is likely real, suggesting that some extended emission is not covered by the X-shooter slits. The emission could then extend beyond 5 kpc from the nucleus in the direction perpendicular to the slits.

Absorption line analysis
We analysed the absorption lines using standard multi-component Voigt profile fitting and vpfit v10.3 (Carswell & Webb 2014) to obtain redshifts, Doppler parameters, and column densities of different species. We used both the UVB and VIS spectra with resolving powers of about 5400 and 8900, respectively. While we fitted all absorption lines simultaneously, the description of the analysis is split into sub-sections for convenience.

Neutral atomic hydrogen
The total H i column density was obtained not only from fitting the damped Ly-α absorption together with the unabsorbed quasar continuum (see Sect. 3.5), but also by including other lines from the Lyman series, the central redshift of the overall H i absorption also being constrained by the metal lines. We obtained a total column density log N(H i) = 20.71 ± 0.02, consistent with the previously estimated value from the lowresolution, low-S/N BOSS data. The region where extended emission is seen was ignored during the fitting process. The result of the fit is shown in Fig. 8.

Molecular hydrogen content
Bluewards of λ obs = 4100 Å, the spectrum is crowded with numerous H 2 lines from different rotational levels (Fig. 10), confirming the detection from Noterdaeme et al. (2019) based on a low-resolution, low-S/N BOSS spectrum. We modelled the H 2 absorption profile using three components (which we label A, B, and C from blue to red). The strongest two components (B and C) are located close to each other in velocity space, so their respective absorption lines are strongly blended. Notwithstanding, these two components are also seen in the neutral chlorine profile (with total log N(Cl i) ≈ 13.4 ± 0.1, see Fig. 11), which is known to be chemically linked to H 2 (see e.g., Balashev et al. 2015 and references therein).
We first used standard multi-component fitting using vpfit, including rotational levels up to J = 5. Because of the decreasing S/N and increasing line blending towards the blue, only the data at λ > 3600 Å was used to constrain the fit of the molecular hydrogen lines. We tied together the redshifts and Doppler parameters between rotational levels for a given component.  well mixed, which is not necessarily the case. Indeed, an increase of b with J has been observed with VLT/UVES in a few cases of relatively low H 2 column densities (e.g., Noterdaeme et al. 2007;Balashev et al. 2009;Albornoz Vásquez et al. 2014). However, in this work, the fit to the lower rotational levels is independent of the exact b value since the lines are damped.
Given the complexity of the H 2 absorption spectrum and to account for the possibility of different Doppler parameters, which can have an effect on the high-J levels, we also fitted the H 2 lines with a Markov chain Monte Carlo (MCMC) procedure. During this fit, we again tied the Doppler parameters of J = 0 and J = 1, as these levels are mostly co-spatial, but those from other rotational levels were allowed to vary independently. However, we added a penalty function to the likelihood function to favour the solutions where Doppler parameters increase with J. We also included the H 2 absorption from J = 6 and J = 7 levels: these are albeit marginally detected. To get posteriors on fitting parameters, we used an affine sampler (Goodman & Weare 2010) with a flat prior on log N, b, and z. The result from the MCMC fit is shown in Fig. 10, but we note that the standard vpfit procedure provides a visually almost indistinguishable profile. Indeed, higher spectral resolution remains necessary to ascertain the velocity structure and widths of the lines. Hence, we provide the best-fit parameters from both fitting methods in Table 3.
In all three components, the column densities of high rotational levels (J 2) are significantly enhanced compared to what is seen in typical H 2 -bearing DLAs at high redshifts. The observed T 01 temperatures (∼100 K for components B and C and closer to ∼200 K for component A) remain consistent with what is seen for intervening systems, albeit with a small tendency to be in the upper range for the observed column densities (e.g., Balashev et al. 2019). This could result from an enhanced heating of the gas through photoelectric effect of UV photons on dust grains. We discuss the observed H 2 excitation diagrams quantitatively in Sect. 5.2.

Chemical enrichment
Metal absorption lines are seen spread over 600 km s −1 mostly bluewards of the systemic redshift. Since the H i column density is large enough for self-shielding from ionising photons, we focused our analysis on the low-ionisation species in order to estimate the gas phase metallicity, since these are in their dominant ionisation stage in the neutral gas.
The initial guess for the velocity decomposition was obtained from the absorption lines of singly ionised silicon and iron that have several transitions with a wide range of oscillator strengths.
Singly ionised sulphur and zinc were subsequently added to the model, and finally, all parameters were fitted simultaneously under the usual assumption that Zn ii, S ii, Si ii, and Fe ii are located in the same gas (i.e. we tied their Doppler parameters and redshifts). Since Zn iiλ2026 is blended with Mg iλ2026, we also included this line in the fit with an independent Doppler parameter and redshift. We note that Si iiλ1526 is detected in the overlap between the UVB and VIS spectra. We therefore included both these spectral regions to constrain the fit.
The result of the fits are shown in Fig. 12 and the corresponding parameters are provided in Table 4. The total column densities are log N(cm −2 ) = 15.52 ± 0.03 (Si ii), 14.46 ± 0.02 (Fe ii), 15.39 ± 0.08 (S ii), and 13.34 ± 0.07 (Zn ii). Indeed, we measured depletion by a factor of two for silicon, while 96% of iron is locked into dust grains.

Excited carbon and silicon
Neutral carbon absorption lines are detected in the spectrum of J0015+1842. This is not surprising given the presence of H 2 and the high metallicity of the system (Noterdaeme et al. 2018). Three components are detected, matching those seen in H 2 , and we detect all three fine structure levels of the ground state triplet (2s 2 2p 2 3 P J=0,1,2 ). All C i (J = 0), C i* (J = 1), and C i** (J = 2) lines were first fitted together using the bands at 1277, 1280, 1328, 1560, and 1656 Å, avoiding regions affected by blends (in particular for the 1328 Å band) and tying the redshifts and Doppler parameters for each component. The result of the fit is shown in Fig. 13, and the corresponding parameters are summarised in Table 5. The Doppler parameters are found to be as small as 1 km s −1 . While such small values are similar to what is seen in H 2 , and not uncommon for C iabsorbers (e.g., Noterdaeme et al. 2017), here they are observed far below the spectral resolution, implying large uncertainties on the column densities. We therefore repeated the fit using the MCMC method, which suggested that the uncertainties on the b-values provided by vpfit were likely underestimated. Finally, we performed a final test assuming a fixed, somehow higher b = 3.5 km s −1 value. The column densities obtained from the different methods remain similar, although differing in some cases by up to ∼0.5 dex. Nevertheless, all results indicate higher excited-to-ground, fine-structure level ratios than usually seen in intervening systems. We caution, however, that high spectral resolution observations remain vital to accurately measuring the corresponding column densities.
Singly ionised carbon and silicon are also detected in excited fine-structure states (C ii* and Si ii*, respectively). While C ii* is in principle useful in estimating the cooling rate of the gas through [C ii]158 µm emission (Pottasch et al. 1979), the corresponding absorption lines here are saturated, impeding the measurement of the C ii* column density. The Si ii* absorption is in turn much weaker but clearly detected in its strongest line at 1264 Å. We identify three main A108, page 9 of 19 A&A 646, A108 (2021) Notes. (a) Doppler parameter tied for all rotational levels (vpfit) and for J = 0 and J = 1 (MCMC). (b) The absorption lines corresponding to these levels are little sensitive on the Doppler parameter.
components that correspond to those seen in the metal lines. We also used the weaker 1309 Å line to ascertain the detection and constrain the fit, shown in Fig. 14. The corresponding parameters are given in Table 4.

Dust extinction
We fitted the dust reddening using the formalism by Fitzpatrick & Massa (2007) assuming various fixed extinction laws (Gordon et al. 2003) to identify the overall shape of the dust reddening law applied to the composite X-shooter spectrum from Selsing et al. (2016). The details of the method are presented by Noterdaeme et al. (2017). We found a best match using a rather shallow extinction curve like the one observed towards the Large Magellanic super-shell (LMC2). However, we needed to vary the strength of the extinction bump at 2175 Å in order to properly match the observed spectrum. The spectrum exhibits significant variations of the iron emission lines in the region between the C iv and Mg ii lines. We included a model of these iron lines using the template given by Vestergaard & Wilkes (2001). While this template improves the overall fit, the exact variations of the myriad of iron lines cannot be captured by a fixed line template. Taking systematic uncertainties due to variations in the intrinsic quasar shape into account, we found a best-fit solution of A(V) = 0.4 ± 0.1 mag, Fe 2 = 0.5 ± 0.2, Fe 3 = −1.5 ± 0.5, c 3 = 0.7 ± 0.1, where A(V) is the amount of extinction and Fe 2 and Fe 3 multiplicative factors parametrising the strengths of the Fe ii and Fe iii lines with respect to the template. The inferred strength of the 2175 Å bump is A bump = πc 3 /(2γ) = 1.2 ± 0.2, where γ = 0.95 is the width of the bump (Gordon et al. 2003). The best-fit solution is shown in Fig. 15. Based on the inferred extinction law from the spectral fitting, we constructed a dust model to reproduce this extinction law by varying the dust-grain size distribution and the ratio of silicate to graphite grains (Si/C). The grain size distribution was assumed to be a power-law distribution with a slope of −3.5 and a maximum grain size of 0.25 µm similar to the standard ISM size distribution (e.g., Nozawa & Fukugita 2013). The best match to the extinction law was found for a minimum dust grain size of 0.01 µm (compared to 0.005 µm for standard ISM dust grains) and a silicate-to-graphite ratio of Si/C = 2 relative to the standard ISM abundance ratio.

Results and discussions
In this section, we first discuss the implications of the emission line properties and why this suggests that we are observing galactic-scale outflows, and we show that the absorbing gas is likely a different manifestation of the same galactic-scale outflows, intercepted by the line of sight to the nucleus. We then derive the physical conditions in the absorbing gas and its distance to the central engine. Finally, we substantiate the overall picture by modelling the Ly-α absorption-emission profile both spatially and spectrally.

The extended emission from ionised gas
According to the standard paradigm of AGN, the broad emission lines of quasars arise from dense regions with high velocities close to the central engine. Since these regions are far too small (∼1 pc, e.g., Kaspi et al. 2007) to be resolved at cosmological distances, information about their sizes is generally obtained through reverberation mapping (Blandford & McKee 1982). However, valuable constraints can also be obtained from partial coverage by a foreground absorber (e.g., Balashev et al. 2011).
The narrow emission lines are in turn expected to arise further away from the obscuring torus, from gas at sufficiently low density where the emission is dominated by forbidden-line transitions. While initially thought to be less than kpc-sized, evidence of the presence of extended narrow line regions (NLR) around diverse types of AGN has increased. This is particularly evident for low redshift Type 2 AGN, where the observations are facilitated by larger angular sizes and obscuration of the nucleus emission (e.g., Sun et al. 2017). However, evidence for kpc-scale emission is also seen in high-redshift quasars (e.g., Rupke et al. 2017). This also includes the ubiquity of giant ( 100 kpc) Lyα nebulae around luminous quasars at 3 < z < 4 (Borisova et al. 2016). The narrow emission lines therefore contain important information for the study of ionised gas on large scales.
If the energy from accretion disc winds propagates to the host galaxy ISM, then we expect to observe some connection between the properties of the different types of emission lines. For example, Coatman et al. (2019) observed a correlation between the blueshifts of the overall C iv and [O iii] emissions, as well as an anti-correlation between [O iii] equivalent width and C iv blueshift. This suggests that high-velocity accretion disc winds, arising in a compact region very close (<1 pc) to the supermassive black hole (Netzer 2015) can influence the host-galaxy on kpc-scales, assuming [O iii] is extended over such scales 2 . While our observation here concerns a single object, the observed kinematics and equivalent width are in agreement with those observed by these authors. More importantly, we also confirm that most of the C iv emission arises from the broad nuclear region, but we observe a weak C iv emission at least as extended as the narrow [O iii] emission, supporting the above picture.
Finally, Tombesi et al. (2015) found evidence for accretion-disc winds driving galactic-scale molecular outflows, as later confirmed by Veilleux et al. (2017). Taken together, these works suggest that the H 2 absorption system observed here could be linked to the same outflow process.
The As shown in Sect. 3.4, the narrow [O iii] emission presents two blobs on both sides of the nucleus (at a distance of ∼5 kpc projected on the plane of the sky) and with opposite line-of-sight velocities with respect to our reference redshift. The bulk velocities of the emission are of several hundred km s −1 . These are at the limit between what is expected from rotational and nongravitational velocities (e.g., Comerford et al. 2018). However, the emission is well detached from the nucleus (in particular the red lobe, see Fig. 7) which is better explained by an outflow. This is also consistent with the presence of blueshifted broad  Fig. 6 reveals that as the distance from the nuclear emission increases, the line width of the narrow emission decreases, and the centroid of the emission approaches the systemic velocity (see right panels of that figure). These kinematic signatures are consistent with a decelerating outflow and a possible compression at the front of the outflow. However, this interpretation remains tentative and requires confirmation through deeper  observations with higher spatial resolution, ideally using integral field spectroscopy.

Kinematics: Absorption and emission as different manifestations of the same wind
While the velocity profile of neutral hydrogen seen in absorption is not directly accessible because of saturation of H i lines, from Fig. 16 it is clear that the profiles of low-ionisation and high-ionisation metal species share the same overall structure, with several components spread over about 600 km s −1 , almost completely blueshifted with respect to the assumed systemic redshift. Such a broad, blueshifted structure is consistent with the expected signature of outflowing gas. The kinematics contrast what is usually seen in intervening DLAs, where the velocity spread of low-ionisation metals are generally much smaller (Ledoux et al. 2006) and systematically less extended than more ionised species (see Fig. 4 of Fox et al. 2007). The profiles observed here suggest that ionised and neutral phases could be well mixed up within the gross structure of the outflowing gas. The absorption line kinematics observed in the spectrum of J0015+1842 also corresponds well to the kinematics of the blueshifted, extended, ionised gas seen in emission (the blue component of the extended [O iii] emission), indicating that we are looking at the same phenomenon both in emission and absorption and that the kpc-scale outflowing gas likely contains a mixture of different phases. A similar interpretation has recently been put forward by Xu et al. (2020). Based on the kinematics of [O iii] emission and C iv absorption in seven quasars, the authors suggest that outflows seen in emission and absorption are likely different manifestations of the same wind.

Physical conditions in the absorbing gas
If the absorbing gas is also part of the outflow, which can typically span around a few tens of kpc from the AGN, then we expect to see effects of the quasar proximity on the physical conditions. The presence and excitation of atomic and molecular species provides a unique opportunity to investigate these physical conditions. Singly-ionised silicon is an interesting species since its excited fine-structure level has been detected in only a few high-z intervening DLAs to date (e.g., Kulkarni et al. 2012;Noterdaeme et al. 2015) and most likely in regions of high  pressure (Neeleman et al. 2015). However, it appears to be more conspicuous in proximate DLAs, where the excitation of silicon fine structure could be enhanced by both high densities and strong UV flux (e.g., Ellison et al. 2010;Fathivavsari et al. 2018). Hence, its presence may indicate close proximity to the quasar central engine. However, in the absence of other information, it is hard to discriminate between the case of high density or high UV flux. Molecular hydrogen can provide important information, since both the presence of H 2 and the excitation of its rotational levels depend on the strength of the UV field, the gas density and its temperature (e.g., Balashev et al. 2019). Owing to formation on dust grains and destruction through line absorption in the Lyman and Werner bands, the presence of molecular hydrogen mainly depends on the ratio between the UV intensity and the gas density (e.g., Sternberg et al. 2014). This implies that the density required for an atomic-to-molecular transition to occur follows n ∝ d −2 , where d is the distance to the central UV source.
For the H i column density and metallicity (or equivalently dustto-gas ratio) observed for J0015+1842, the normalisation of this relation 3 is about log n (cm −3 ) ∼ 4 at d ∼ 10 kpc distances from the nucleus .
Combining the information from Si ii* and H 2 , as well as other species sensitive to the density and UV radiation (e.g., C i fine structure, whether CO molecules are present or not), can therefore provide constraints on the physical state of the gas and its distance to the central UV source. In Fig. 17, we compare the apparent optical depths of the Si iiλ1808 and Si ii*λ1264 lines 4 , which are both optically thin. The ground state and excited profiles of singly-ionised silicon have similar velocity profiles with three main clumps, two of which coincide with the positions of H 2 absorption components. Yet, the strongest bulk component of Si ii* does not align with the H 2 components. In the top panel of Fig. 17, we show the ratio of Si ii*/Si ii. This ratio is not enhanced at the position of H 2 absorption, where instead it tends to be consistently lower than average over several pixels. Since molecular hydrogen traces the high-density regions in the absorbing gas, this may suggest that the high excitation of Si ii* is not due to high density but rather due to strong UV irradiation. 3 We note that this simple estimate assumes the transition to occur and does not exclude the presence of H 2 at lower densities (or distances) if the transition is not complete. 4 We caution that Si ii*λ1264 is partly blended with Si ii*λ1265. However, the latter has an oscillator strength ten times smaller than the former, so the impact on the apparent optical depth remains small.  Fig. 12 for the absorption lines from excited finestructure of singly-ionised silicon.

Modelling warm and cold gas as successive slabs
Since the quasar itself is likely the dominant source of photons striking the medium at distances few tens of kpc, the UV field received by each slab of gas (absorption component) actually depends on the transmitted field through the previous slab of gas on the line of sight. While understanding successive filtering slabs of gas makes the study of proximate systems significantly more complex than intervening ones (where the components can generally be treated independently), advantages are that the observed column densities correspond to the actual layers through which UV photons propagate, and we have prior knowledge about the shape of this field 5 . However, the large number of components and their unknown respective physical locations along the line of sight does not allow for a full modelling of the system, which would have far too many parameters to vary.
We simplified the problem by considering two successive slabs of gas along the line of sight to get an idea of the physical conditions in the different phases. In this picture, one slab of warm gas, directly illuminated by the quasar, is likely responsible for most of the Si ii*, while the second slab, for which the incident field contains no more H-ionising photons, is in a cold state and contains the H 2 .
We used the photoionisation code Cloudy v17.1 (Ferland et al. 2017), running grids of constant density models with plane-parallel geometry illuminated on one side by the built-in cloudy AGN spectral energy distribution. The cosmic microwave and cosmic ray backgrounds were included. The gas phase abundances were set to match the observed ones, and we included dust using the grain mixture and size distribution described in Sect. 4.5. In the first runs, the abundance of dust was scaled to self-consistently (i.e. conserving mass) match the observed depletion factors of iron and silicon. This in turn implies depletion factors of about 50% and 25% for carbon and oxygen, respectively. We tested and found that the results concerning the first layer are actually little sensitive to the abundance of dust, except, naturally, for the calculated extinction, which tends to strongly exceed the observed one due to the large total column density in the model. However, as we see below, the first layer contains mostly ionised gas, as well as neutral gas with high temperatures (>10 4 K), so it is most likely that grains, in particular small ones, are destroyed. We then assumed no grains in that phase -but these are kept in the second (cold) slab.
The results concerning the excitation of silicon in the first slab are shown in Fig. 18, where the UV flux is converted into 4 × 10 3 6 × 10 3 2 × 10 4 3 × 10 4 Observed wavelength (Å)  Fig. 15. Combined X-shooter spectrum of J0015+1842 in black with best-fit SED corrected for dust extinction shown in red. For visual purposes, the data have been filtered (rejecting 4-σ outliers; light grey) and smoothed (black) in order to remove the strong skyline residuals in the NIR data. The position of the 2175 Å extinction feature is indicated by an arrow, and the best-fit model without the 2175 Å feature is indicated by the purple line. The unreddened quasar template is shown in blue. Photometric data from SDSS (u, g, r, i, and z bands) and WISE (band 1) are shown as orange squares. distance to the source using the observed quasar luminosity.
These show that the observed excitation of Si ii* can be reproduced by gas with densities of the order of ∼10 cm −3 at any distance less than ∼50 kpc from the source. Exciting Si ii at the observed level, if located further, would require higher densities. However, this would also push the second slab further out, where the cold slab of gas containing H 2 would become too cold and little excited compared to observations, as we discuss later. We note that the two slabs need not have the same velocities since the ionising cross-section is wide.
In Fig. 19, we show the variation of number densities of neutral hydrogen, electrons, and singly-and thrice-ionised silicon, as well as the excitation of silicon and the electron temperature as a function of the cumulative H i column density for an example model with n = 10 cm −3 and d = 10 kpc. In the ionised gas, the excitation of Si ii is high owing to efficient collisions with electrons. The local Si ii*/Si ii ratio then drops by about a factor of ten after crossing the dissociation front, when the gas becomes neutral and both the electron density and the temperature decrease. The Si ii*/Si ii ratio then remains constant afterwards, as do n e and T . However, since singly- ionised silicon is not the main ionisation stage upstream, the overall column density ratio, N(Si ii*)/N(Si ii) (which is what we observe), is actually N(H i)-weighted, explaining the slow dependence on the stopping H i column, as shown by the cumulative N(Si ii*)/N(Si ii) ratio 6 (dashed-dotted curve) in the middle panel.
Changing the UV field affects the depth of the dissociation front (i.e. the column density of upstream ionised gas), but the results remain almost unchanged with respect to the cumulative H i column density. This explains the low dependence of the N(Si ii*)/N(Si ii) ratio on the distance to the UV source. It is interesting, however, to note that the range spread by Si ii*/Si ii at different depths (from about 0.01 in the ionised gas to 0.001 in the neutral gas) corresponds roughly to that seen in Fig. 14, which is consistent with the fact that the absorption system also contains ionised gas. The main effect of increasing the volumic density is an upwards shift of the Si ii*/Si ii curve, as collisional excitation of Si ii increases both in the ionised gas (collisions with electrons) and in the neutral gas (collisions with H i), so Si ii*/Si ii is roughly proportional to n in the UV-independent regime. At much higher densities, the increasing dependence of Si ii* on the distance to the AGN is mostly due to the dependence of the gas temperature and ionisation fraction on the strength of the UV field. Finally, we note that H 2 is not yet produced in the first slab, which is consistent with the non-detection of H 2 in the strongest Si ii* component.
We modelled the second slab of gas again using constant density models at various distances from the quasar. However, the striking radiation field is now the transmitted field through the first slab. The main effect of the first slab on the incident field is a removal of hydrogen-ionising photons. This means that the assumed parameters for the first slab actually have little effect on the second slab. We stopped the calculations of this second slab when reaching log N(H 2 ) = 19.3, as a representative value for individual H 2 components. The results are shown in Fig. 20.
The excitation of Si ii* is now significantly lower compared to the first slab for a given distance and density as the temperature decreases considerably. The Si ii*-to-Si ii ratio reaches the same level as in the first layer for distances of only a few kpc, for densities up to log n ∼ 5, and it increases at higher densities. However, there is also indication that the excitation of silicon is less at the velocities of H 2 components, with log(Si ii*/Si ii) between roughly −3 and −2.5, than at other velocities (where, as we showed above, the high excitation is easily explained by the high pressure in the first layer). In that case, the distance is slightly larger, between 5 and 12 kpc.
The T 01 excitation temperature strongly depends on the UV flux, but more weakly on the density. The observed value of about 100 K corresponds to distances of the order of 5-15 kpc. At such distances, the high rotational levels are highly populated, even though they remain lower than what is observed for the total column in the model (equating that of component C). This is a common situation in the modelling of static slabs and depends much on the actual geometry and turbulence, so we do not consider the 0.5 dex as a strong mismatch between the data and model. The excitation of neutral-carbon fine-structure levels has also proved to be a useful probe of the physical conditions in H 2bearing DLAs (e.g., Silva & Viegas 2002), in which it is usually detected (e.g., Srianand et al. 2005;Noterdaeme et al. 2018).
Here, the C i*/C i ratio is notably high in all three components, reaching the limit of a statistical weight ratio, and it appears to be mostly populated by UV pumping. The C i*/C i favours distances smaller than the above 10 kpc, it but remains consistent within the large uncertainties of the observations over most of the parameter space probed. For component C, where the constraints are slightly better, the excitation of C i rejects distances larger than about 10 kpc at ∼2σ.
The strict upper limit on the H i column density is consistent with these distances as long as the density is as high as log n ∼ 4.5 or more. The dust extinction is then also within the allowed range. Finally, the absence of CO at a detectable level (log N(CO) < 13.4 (3σ) following the procedure presented by Noterdaeme et al. 2018) rejects the high density, large distance corner of the plot.
Our model is of course very simplified, as the actual system presents numerous components at various velocities and unknown respective spatial location along the line-of-sight. The chemical properties, gas-phase metal abundances, and dust grain properties can also vary from one component to another. We also explored a number of varying parameters, like turbulence, shape of incident radiation field, column density of the first layer, and the use of the simple removal of H-ionising photons instead of a full transfer calculation. All this typically affects the distances and densities further by a factor of two. We thus conclude that within a factor of a few, the distances of the clouds are of the order of ∼10 kpc, and H 2 is produced in clouds with densities as high as n ∼ 10 4.5 cm −3 , located downstream, after H-ionising photons have been removed by a warm atomic layer. The latter is also responsible for most of the Si ii fine-structure excitation, while C i is only present in the denser molecular gas.

Modelling the excitation of H 2 components individually using Meudon PDR
In order to further test our results on the molecular phase, we used the Meudon PDR code (Le Petit et al. 2006) to model the population of all detected H 2 rotational levels in the different individual components separately. We ran grids of isobaric cloud models, with thermal pressure from 10 4 to 10 8 k B K cm −3 and illuminated by a UV field with intensity from 1 to 3 × 10 4 times the Mathis field (Mathis et al. 1983), corresponding to distances of ∼160 to 1 kpc from the central engine 7 . We modelled each H 2 component independently, that is, we did not take into account shielding in H 2 UV lines of one component by another component located closer to AGN. The column densities in the different rotational levels are all compared to the data using a simple least-square likelihood function. The constraints on the thermal pressure, P, and distance to AGN, d for all three H 2 components are shown in Fig. 21. The constraints are similar for vpfit and MCMC H 2 fitting methods. The most likely 7 Assuming the UV field is dominated by the AGN. We note that there is no option to set the external radiation field shape to AGN in the PDR Meudon code. However, in the UV range, the shape of the Mathis field is quite similar to that of typical AGN, so that uncertainties related to the shape of the field are negligible compared to those arising from the unknown geometry of the clouds.  Table 3), respectively. The marginalised probability distribution functions (along right and top axes) and the best-fit values (stars) are shown for MCMC results only for the sake of clarity. The black line with upwards arrows show the constraints from the H 2 -H i transition, assuming the kinetic temperature is T ∼ 100 K, which is consistent with the observations. values correspond to high pressures P/k B ∼ 10 6 −10 7 K cm −3 located at distances of 10−20 kpc for all three components. The corresponding modelled excitation diagram is shown in Fig. 22. The typical kinetic temperatures in the H 2 dominated cores of the clouds is then ∼200−300 K, which is slightly higher than the observed T 01 . The number density in the H 2 -dominated medium is in the range of 10 4 −10 5 cm −3 , which is in agreement with the results using Cloudy. We also ran isochoric models and found that while isobaric models reproduce the observed H 2 excitation diagram slightly better, the constrained values of d, P, and n remain very similar. We note that the derived distances are also consistent with the limits predicted by the static H 2 -H i transition theory from Sternberg et al. (2014), shown as a black line in Fig. 21: below these lines, the predicted H i columns in the envelope of H 2 clouds would be higher than the total N(H i) observed.
While the distances are constrained to be of the order of 10 kpc, these do not allow us to determine which component is located closer to the AGN. In any case, since cold clumps can drop out from the bulk motion (see e.g., Wu et al. 2020), the information of the respective locations of H 2 clouds would not likely tell us whether the bulk outflow is decelerating or not.
Interestingly, the high density derived in the H 2 -bearing medium implies very small longitudinal sizes, of the order of ∼100 AU only. This is actually comparable to the size of the accretion disc, assuming the relation from Morgan et al. (2010) and the black hole mass of log M BH /M = 8.6 derived from width of the C iv line (FWHM ∼ 3000 km s −1 ) and the restframe continuum luminosity at 1350 Å using the calibration from Vestergaard & Peterson (2006). However, we do not see any partial coverage of the background source by H 2 . This could be due to the nuclear emission arising from a smaller area than expected, or the transverse extent of H 2 -bearing gas being larger than the longitudinal one. The latter could be effectively reproduced by a large covering factor of H 2 -bearing clumps around the component velocity. Given the typical (longitudinal) velocities of the components of ∼400 km s −1 future observations on  Fig. 21). The data points are artificially shifted on the x-axis to ease representation.
timescales of a few years could reveal a variability of the H 2 absorption lines, if the transverse velocities are of the same order.

The leaking Lyα emission
Sensitive integral field spectrographs, in particular MUSE at the VLT (Bacon et al. 2010), have revealed a near ubiquity of Ly-α nebulae around high-redshift quasars (e.g., Borisova et al. 2016;Arrigoni Battaia et al. 2018). Extended Ly-α emission has been observed up to hundreds of kpc in some cases, and routinely over several tens of kpc. In turn, the strong glare from the quasar itself complicates the observations of the Ly-α emission at the small distances from the nucleus. In the case of J0015+1842, the presence of the DLA removes the nuclear emission, leaving us with free access to the Ly-α emission as soon as it is outside the BLR.
To verify if the conjectured configuration of clumpy, gaseous outflows from J0015+1842 is a likely representation of the system (see cartoon in Fig. 24), we constructed a 3D model in an adaptively refined mesh and carried out Ly-α radiative transfer (RT). The model is based upon our observational constraints, as discussed below; where we have no constraints, we assumed plausible values (in particular concerning the properties of the inter-cloud medium). In both cases, we then varied the physical parameters in certain ranges to get a feeling for the importance and degeneracies of these parameters.
The Ly-α RT was conducted using the Monte Carlo code MoCaLaTA 8 (Laursen et al. 2009a,b). Given a grid of cells containing thephysical parameters of the gas (densities, temperature, and velocity field), the code traces the paths of individual photons, producing a 3D (two spatial + one frequential) cube that can be compared to our observation. Given the large number of unknowns, rather than performing a stringent fit to the observed spectrum (as has sometimes been done, see e.g., Verhamme et al. 2008;Gronke 2017), we aimed to construct a model that would at the same time be consistent with our various observational constraints, and -most importantly -roughly reproduce both the observed Ly-α spatial and spectral shapes.
The gas is assumed to have been expelled in a biconical outflow down through one of which we observe the quasar nucleus. Given the similar transverse and longitudinal length, we set the opening angle of the cone to θ 45 • . The number Comparison between simulated 2D spectrum of photons escaping in the particular direction of the observer looking down the barrel of the jet (bottom) with the actual observed X-shooter 2D spectrum (top). The distance along slit is taken from the unresolved quasar continuum emission. We note that the observed spectrum is also affected by unrelated absorption, which we do not model here (e.g., at 4396 Å). of low-ionisation absorption lines seen along the line of sight (covering factor f c ∼ 10) can then be reproduced distributing N cl 750 000 effective clouds in the outflows, each of radii r cl 40 pc, T cl = 10 4 K and n HI,cl = 0.25 cm −3 . These are located between inner and outer radii of 16 and 20 kpc, respectively 9 and immersed into a mostly ionised inter-cloud medium with T ICM = 10 5 K and a residual neutral gas n HI,ICM = 0.005 cm −3 . Both phases have metallicity Z = 0.5 Z . With these numbers, the covering factor is f c = 10 +3 −2 H i column density is log N(H i) = 20.63 +0.09 −0.11 . The gas follows a velocity profile decreasing linearly from 250 km s −1 to 0 at the end of the jet. In addition to this, the clouds have a velocity dispersion of σ v,cl = 10 km s −1 . Finally, photons are emitted from the centre with an intrinsic spectrum given by the reconstructed unabsorbed emission (see Fig. 8).
We generated a 2D spectrum of the photons escaping (in the particular direction of the observer looking down the barrel of the jet) as would be seen with X-shooter, that is, with the same slit parameters, resolution, and seeing as observed. This is compared to the observed spectrum in Fig. 23. While far from a perfect match, our model reproduces (i) the damped absorption line towards the nucleus, (ii) the red-peaked leaking Ly-α emission, and (iii) the spatial offset of this emission down the trace. We note that Ly-α photons are also likely to be produced from recombination in the outflow as well as in the host galaxy, hence contributing to the overall observed profile.

Conclusion
Active galactic nucleus feedback is a very complex issue that motivates significant efforts from the community in terms of theory, simulations, and observations. While progress in understanding galactic-scale outflows as the main observable evidence of feedback is being made quickly, Cicone et al. (2018) recently noted that these remain largely unconstrained as they imply various scales and different gas phases. Furthermore, the observation of outflows, especially at high redshift, is notoriously difficult. 9 These values are taken to be representative of the distances derived from the physical conditions and the fact that clouds are more easily ionised close to the quasar. As we aim at a global agreement between the model and actual data and still have many unknowns, we chose to fix them instead of considering these radii as free parameters. can be more easily scattered out from the cones towards the observer from the receding outflow and appear red peaked.
Our VLT/X-shooter observations of the proximate molecular absorber towards J0015+1842 suggest that spatially resolved emission of ionised gas and kinematically wide absorption in various phases towards the nucleus are different manifestations of a single outflowing process. From detailed investigation of the physical conditions in the cold phase of the absorbing gas, we derived a distance to the central engine of the order of 10 kpc. This, together with the similar kinematics and transverse extent of the ionised gas seen in emission, indicates a picture of a multi-phase quasar outflow, down the barrel, that is, intercepted by the line of sight to the bright nucleus. This configuration is corroborated by modelling of the resonant Ly-α emissionabsorption spectral and spatial profiles.
Of course, it may be reasonable to recall that this is likely an extremely simplified version of the actual geometry. It is, for example, well possible that the ISM in the quasar host is disrupted by the ionised outflow, giving rise to a much more complex structure of the various phases in the outflow. While the present observations open up a new way to study outflow properties, they do not yet tell us what the effect of these outflows on the host galaxy are. Understanding how much molecular outflows affect the host galaxy's reservoir, and hence whether outflows are effectively quenching star formation, could be addressed from spatially resolved observations of CO emission lines (e.g., Carniani et al. 2017).
Since the present analysis includes only a single object, it is legitimate to ask whether the prevalence of (multi-phase) outflowing gas is larger in quasars with proximate H 2 absorbers than in the overall population. By spreading out H 2 gas, outflows may contribute to increasing the H 2 absorption cross-section, and hence proximate H 2 absorbers could represent an effective way of pre-selecting these. If, as shown here, outflowing gas facilitates the scattering of Ly-α photons on galactic scales, then the presence of leaking Ly-α emission in a fraction of them could be an additional effective criterion. Finally, the presence of excited atomic species could indicate proximity of the gas to the central engine. Outflowing gas being responsible for a fraction of proximate H 2 absorbers could provide an explanation to the at least five-fold enhancement of the incidence rate of H 2 -selected proximate DLAs compared to intervening statistics ) when H i-selected proximate DLA show only marginal enhancement (Prochaska et al. 2008) and likely missed systems with strong Ly-α leakage. Of course, a A108, page 18 of 19 fraction of proximate H 2 absorbers could also originate from galaxies in the quasar group.
We note that a correlation between a fraction of leaking Ly-α photons (not absorbed by the foreground cloud) and the strength of Si ii* in metal-selected proximate DLAs (without information of H 2 ) has been observed by Fathivavsari (2020). The author interprets this as a decrease of the covering fraction of the background Ly-α emission by clouds with increasing compactness (i.e. of higher density for a given column density). However, there may be no causality between the density inside individual (sub)-pc-scaled clouds and the fraction of leaking Ly-α photons (over kpc-scales) throughout the whole system. Instead, the prevalence of dense clouds and the Ly-α leaking fraction may both depend on the gas flow configuration, including the outflow viewing angle (see e.g., Gupta & Saikia 2006 for a similar discussion on the detectability of H i 21 cm absorption) and quasar properties. Similarly, den Brok et al. (2020) recently suggested that the morphologies of extended Ly-α nebulae can be used to understand the geometry of high-redshift AGN on circumnuclear scales. Nevertheless, the interpretation of single clouds partially covering the nucleus may still be valid for extreme cases where the Ly-α appears completely unabsorbed.