Proximate molecular quasar absorbers: Chemical enrichment and kinematics of the neutral gas

Proximate molecular quasar absorbers (PH 2 ) are an intriguing population of absorption systems that was recently uncovered through strong H 2 absorption at a small velocity separation from the background quasars. We performed a multi-wavelength spectroscopic follow-up of 13 such systems with VLT / X-shooter. Here, we present the observations and study the overall chemical enrichment measured from the H i , H 2 , and metal lines. We combined this with an investigation of the neutral gas kinematics with respect to the quasar host. We ﬁnd gas-phase metallicities in the range 2% to 40% of the solar value, that is, in the upper-half range of H i -selected proximate damped Lyman-α systems, but similar to what is seen in intervening H 2 -bearing systems. This is likely driven by similar selection e ﬀ ects that play against the detection of most metal-and molecule-rich systems in absorption. Di ﬀ erences are seen in the abundance of dust (from [Zn / Fe]) and its depletion pattern when compared to intervening systems, however, possibly indicating di ﬀ erent dust production or destruction close to the active galactic nucleus. We also note the almost ubiquitous presence of a high-ionisation phase traced by N v in proximate systems. In spite of the hard UV ﬁeld from the quasars, we found no strong overall deﬁcit of neutral argon, at least when compared to intervening damped Lyman-α systems. The reason likely is that argon is mostly neutral in the H 2 phase, which accounts for a large fraction of the total amount of metals. We measured the quasar systemic redshifts through emission lines from both ionised gas and CO(3–2) emission, the latter being detected in all six cases for which we obtained 3 mm data from complementary NOEMA observations. For the ﬁrst time, we observe a trend between the line-of-sight velocity with respect to systemic redshift and metallicity of the absorbing gas. This suggests that high-metallicity neutral and molecular gas is more likely to be located in outﬂows, while low-metallicity gas could be more clustered in velocity space around the quasar host, possibly with an infalling component.


Introduction
The co-evolution of active galactic nuclei (AGNs) and their host galaxies is regulated by galaxy processes, including secular evolution, inflow rate of gas (from the intergalactic medium to the interstellar medium and ultimately to the supermassive black hole (SMBH)) as well as merging history (see Hopkins et al. 2008).The apparent properties of AGNs may further depend on the orientation and the evolutionary stage at which the observations are taken (e.g., Urry & Padovani 1995).This constitutes an entire field of modern astrophysics, motivating large efforts from the community, in both theoretical and observational grounds.Because of the large diversity of observed properties, it is common to select objects based on a few among their many characteristics (e.g., redshift, brightness, colour, radio-Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 103.B-0260(A) and 105.203L.001.loudness, variability, etc.) and study the relationships between other marginalised properties.Each of these studies then bring a different and very valuable lightening angle to the overall picture.
As cold gas constitutes the fuel for both star-formation and the growth of the supermassive black holes, many works have focused on detecting the dense molecular component of the brightest AGNs -quasars-through CO emission (e.g., Omont et al. 1996;Weiß et al. 2007;Wang et al. 2016), in order to explore the links between star-formation, total molecular mass in the host and growth and activity of the SMBH (e.g., Bischetti et al. 2017).
On the other hand, absorption studies enable detailed investigations of the gas properties along the line-of-sight to the central engine, over a wide range of densities.For example, ionised winds powered by the accretion disc can be observed as broad absorption line systems (BALs) in roughly 15% of quasars, and have velocities from a few thousands km s −1 up to ∼0.3 c (Hamann et al. 2018).Damped Lyman-α systems (DLAs) are ab-sorption systems with high enough H i column density (N(H i) > 2×10 20 cm −2 ) to be self-shielded from ionising radiation and are frequently seen along quasar lines of sight (e.g., Prochaska et al. 2005;Noterdaeme et al. 2009).While most DLAs are unrelated to the quasar environment and only intercepted by chance (called intervening DLAs), some DLAs are found at roughly the same redshift as the background quasar and hence dubbed proximate DLAs (PDLAs) where "proximate" refers to proximity in velocity space (typically < a few 10 3 km s −1 ).Because of the much smaller redshift path in which to identify PDLAs compared to that available for intervening DLAs, the former are much rarer than the latter.Nonetheless, Prochaska et al. (2008) showed the incidence of PDLAs (per unit redshift) is higher than expected from the intervening statistics meaning that they are likely associated to a clustered environment around the quasar.Ellison et al. (2010) performed the first systematic study of PDLAs at high-resolution and concluded that at least a fraction of them exhibit different characteristics relative to the intervening population.They however suggested that PDLAs are unlikely related to the quasar host, but sample preferably more massive galaxies in the highly clustered region around the quasar.Interestingly, Ly-α emission has been relatively frequently found in the core of several PDLAs (e.g., Møller & Warren 1993;Leibundgut & Robertson 1999;Ellison et al. 2002), when this is more rarely seen among the several 10 4 known intervening DLAs (but see Fynbo et al. 2010;Noterdaeme et al. 2012;Ranjan et al. 2018 for some example of direct detection and Rahmani et al. 2010;Noterdaeme et al. 2014;Joshi et al. 2017 for statistical detection through stacking).Finley et al. (2013) and Fathivavsari et al. (2018) further identified a population of PDLAs with Ly-α emission spanning a range of strengths, up to the point where the H i Ly-α absorption becomes unseen, but the presence of neutral gas is confirmed by strong low-ionisation metal lines.Such emissions are unlikely to arise from star-formation but rather from Ly-α from the AGN that is not fully covered by the absorber.Because Ly-α photons scatter easily and up to large distances around quasars (e.g., Arrigoni Battaia et al. 2019;Borisova et al. 2016), concluding about the location, scale and geometry of the absorbing gas is not straightforward.However, signatures of high excitation of the gas suggest an origin at rather small distance from the AGN, at least for the most extreme cases.Noterdaeme et al. (2019, hereafter Paper I) presented a new approach focusing on the presence of strong molecular hydrogen absorption directly detected at the quasar redshift in Sloan Digital Sky Survey (SDSS) spectra, without prior on the presence of H i or metal lines.Not only H 2 is sensitive to the physical and chemical conditions in the cold gas, but such an approach also differs from studies based on, e.g., CO molecular emission in the sense that no prior is made on the overall mass of dense molecular gas in the galaxy host, but rather on the interception of more diffuse cold gas by the quasar line of sight.The strong excess of proximate H 2 absorbers with respect to intervening statistics shows that most of these systems must be associated with the quasar environment, although it is not clear whether these arise from the quasar host, or if the excess is due to strong clustering of galaxies around the quasar.
To shed light on this issue, we embarked into a follow-up study with the multi-wavelength, medium resolution spectrograph X-shooter on the Very Large Telescope (VLT).A detailed study of a first system revealed a likely origin in outflowing material intercepted by the line of sight to the central engine (Noterdaeme et al. 2021b, Paper II). Remarkably strong and wide CO(3-2) emission as well as 3 mm continuum emission was de-tected in the same, optically bright (moderately reddened) system through NOEMA observations (Noterdaeme et al. 2021a, Paper III).All this suggests that the nucleus is seen through a channel of lower density but where some of the molecular content is carried out by outflows, when the host still contains large amounts of dust and molecular gas with wide kinematics, possibly indicative of a post-merger system.
In this work, we focus on the overall chemical enrichment and kinematics of the gas of our follow-up sample.In Sect. 2 we present the sample, observations, data reduction and postprocessing.In Sect. 3 we present the column density and dust extinction measurements.In Sect.4, we present the emission redshifts and absorption line kinematics.We discuss our findings and compare our results with other PDLAs as well as strong intervening H 2 systems in Sect. 5. We conclude in Sect.6.

Sample and observations
The sample studied here is entirely drawn from Paper I.Among the statistical sample of 50 quasars with high-confidence proximate molecular absorbers, 22 have δ < 20 • , that is, reachable from the Southern hemisphere.We obtained VLT/X-shooter data for thirteen of them: eleven based on their observability during the allocated periods of this specific programme (P103 and P105 1 ) plus the quasars SDSS J013644.02+044039.1 and SDSS J085859.67+174925.1 for which we use data collected in P94 as part of another programme (Balashev et al. 2019)  2 .We note that since the quasars in the present sample were selected from the statistical sample upon their observability only (and not any other property) they are likely a fairly good representation of the latter.
X-shooter (Vernet et al. 2011) simultaneously covers the whole optical/NIR wavelength range from 0.3 to 2.5 µm, splitting the incoming light beam into three spectrographs ('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.All observing blocks (OBs) consisted of 2980 s, 3010 s and 5×600 s exposures for the UVB, VIS and NIR arms, respectively obtained in 'stare' mode.We used 1.0 and 0.9 -wide slits in the UVB and VIS, with slow read-out mode and 1×2 pixel binning.The NIR observations were performed with 1.2 -slits for the two quasars observed in P103 (SDSS J001514.82+184212.3 and SDSS J232506.62+153929.3).The remaining observations were performed with a narrower NIR slit (0.9 arcsec) and a K-band blocking filter to avoid stray light in the instrument, effectively cutting out all wavelengths longer than 2.08 µm for all remaining quasars but SDSS J013644.02+044039.1 and SDSS J085859.67+174925.1 (in P94) for which the Kband filter was not yet available.We executed two to five OBs for each quasar at different position angles.Most of the observations were obtained under good seeing, somewhat sharper than the slit widths at the corresponding wavelengths, resulting in higher than nominal spectral resolution.The observing log is given in Table 1.For convenience, we will refer to these objects using a short name ver-sion in the following: SDSS Jhhmmss.ss+ddmmss.s will be abbreviated as Jhhmm+ddmm (e.g.J0015+1842 instead of SDSS J001514.82+184212.3).

Data reduction and telluric corrections
All X-shooter spectra were processed using the official esorex pipeline version 3.5.3(Modigliani et al. 2010) for 'stare' mode reduction maintaining the default parameters and including the correction for instrument flexure (recipe 'xsh_flexcomp').The raw images were pre-processed using the Python package astroscrappy (a reimplementation of the L.A.Cosmic algorithm by van Dokkum 2001) to interpolate over pixels contaminated by cosmic ray hits using the neighboring pixels.Each arm of each OB was reduced separately using the standard star observed closest in time in order to flux calibrate the given spectrum.The 1D spectra were then extracted using optimal extraction (Horne 1986).
Next, we used the ESO tool molecfit v.4.2.2a.1 (Smette et al. 2015) to correct each individual 1D VIS and NIR spectra for the absorption lines produced by the Earth's atmosphere.In the VIS, the main telluric absorption features are due to O 2 and H 2 O, while the NIR is strongly affected by H 2 O and CO 2 bands.Following the recommendation from the molecfit manual, the correction was done by choosing a few spectral regions with clear but not strongly saturated atmospheric lines for each quasartypically 0.68-0.69,0.71-0.73,0.76-0.77,0.81-0.82and 0.91-0.91µm in the VIS.In the NIR, we had to adjust more the selected spectral regions for each quasar, depending on the strength of absorption bands and location of quasar emission lines.During this process, we also masked deviant pixels (e.g., cosmic remnants) as well as astrophysical lines.
The best-fit model was used to calculate the overall atmospheric transmission and correct the individual VIS and NIR 1D exposures accordingly.The telluric-corrected spectra were then shifted to vacuum-barycentric wavelengths, rebinned to a common wavelength grid and combined using an inverse-variance weighting.The same vacuum-baryocentric correction and combination of 1D spectra was applied to UV data though without telluric correction.
We note that molecfit also provides the actual spectral resolution, that is a free parameter of the transmission model.As expected from a seeing sharper than the slit widths for most of the observations, the resolving power is typically 25% higher than the nominal values in the documentation (R nom = 5400 (UV) and 8900 (VIS)).The average spectral resolution obtained from individual VIS spectra is then adopted when fitting the lines.In the UV, the spectral resolution is first estimated by scaling the nominal UV value by the observed-to-nominal VIS resolution ratio.We cross-checked the values are in broad agreement with those expected for a slit width that would correspond to the observed seeing.It was then adjusted during the fitting process if necessary.While the exact values have little influence on the metal column densities, we provide the adopted values of the resolutions for the combined spectra along with the best-fit model parameters in Appendix B to ease reproducibility.

NOEMA observations
While the present work focuses on X-shooter data, we also obtained observations with the NOrthern Extended Millimeter Array (NOEMA) and the PolyFIX correlator in the 3 mm band between May 2020 and May 2021 for 6 out of the thirteen objects presented here.The central frequency was chosen in order to cover the expected position of the CO(3-2) emission line, which we detected in all 6 cases.A description of data handling and analysis of NOEMA observations (including other quasars from the same parent sample but not observed with X-shooter) will be presented in a future paper but we already use the redshifts from the detected CO(3-2) emission lines in the discussion (Sect.5).The Gaussian fits to the CO(3-2) emission lines are presented in Appendix F.

Column-density measurements
We analysed the absorption lines through standard, simultaneous, multi-component Voigt-profile fitting using vpfit v.12.3 (Carswell & Webb 2014) to measure the column densities, Doppler parameters and redshifts.The quasar continuum was evaluated using spline functions, with nodes constrained by the unabsorbed regions and adjusted over strong lines in an iterative way when necessary.In order to best take into account blending between metal, H 2 and H i lines, the final absorption models were obtained by fitting all lines simultaneously.However, for the sake of presentation, we describe here the analysis separately for H i, H 2 and metals.The results from the fits are presented in the Appendix, with the main derived properties summarised in Table .3.

Total HI and H 2 column densities
The total H i column densities were mostly constrained from the saturated core of the Ly-α line, which is little sensitive to the exact placement of the continuum on top of the Ly-α emission line of the quasars.Notwithstanding, thanks to the wide wavelength coverage of X-shooter and the relatively high redshifts of quasars (implied by the selection upon H 2 in SDSS), all lines from the H i Lyman series are covered by our spectra, allowing to constrain the H i column density further.While the formal statistical errors from fitting the lines are very small (typically 0.01-0.02dex), we note that the continuum placement as well as choice of fitting regions are necessarily somehow subjective, possibly introducing some systematic errors.Allowing for continuum variation using Chebychev polynomials indicate uncertainties less than 0.1 dex for similar DLA observations with Xshooter (e.g., Balashev et al. 2019;Ranjan et al. 2020).We hence added here a conservative 0.1 dex uncertainty into the error budget.
In all the systems, we confirm the presence of strong H 2 absorption lines from the Lyman and Werner bands.These are typically detected from rotational level J=0 to J∼6-7 in our Xshooter spectra.The high column densities imply damped lines at least in the first few rotational levels and hence accurate column densities, since lines in this regime do not depend on the Doppler parameter b.We do also take into account and fit the detected high rotational levels, assuming the same Doppler parameter and redshift for all H 2 lines.However, the derived column densities for the high-J levels are likely to be more sensitive to the assumptions on the Doppler parameters.In fact, b-values for H 2 lines have been found to present a dependence on the rotational level in several cases, see, e.g., Noterdaeme et al. (2007); Balashev et al. (2009).A more detailed study and exploration of systematics is hence needed in order to use these lines to constrain the excitation of the gas and the prevailing physical conditions.This is postponed to a future work, while we focus here on the total column densities, which are dominated by the low- Notes. (a) The actual object designation following the rules from the International Astronomy Union, includes the "SDSS" preceding the J2000 coordinates sequence provided in this column. (b) Delivered seeing corrected for airmass, as measured at the telescope by the Shack-Hartmann wavefront sensor (fits keyword HIERARCH ESO TEL IA FWHM). References.
(1) rotational levels.The fits to the H 2 absorption bands is presented in Appendix A.
In Fig. 1, we compare the total H 2 and H i column densities obtained with X-shooter, with those derived from SDSS.In the case of H i the agreement is very good, and always better than 0.3 dex.For H 2 , for which the involved columns are much lower, the agreement remains remarkable (apart for the lowest column density system with a difference of 0.85 dex), in particular given the fact that the SDSS values were obtained by visually scaling and matching a template (i.e., assuming single component with fixed Doppler parameter and excitation temperature) to the lowresolution, low S/N SDSS spectra (Paper I).This agreement is similar to what was found by Balashev et al. (2019) for intervening H 2 systems.

Low-ionisation metal species
In this section, we focus on the main low-ionisation metal lines that provide reliable column densities, that is, lines that have multiplets, are not strongly saturated and not in a region strongly affected by telluric lines.For example, we did not attempt to fit the only singly-ionised and very strong carbon line C iiλ1334.We fitted the low-ionisation metal lines assuming a common velocity structure with turbulent broadening (i.e., tied redshifts and Doppler parameters).We first identified the velocity components using the strong Fe ii and Si ii lines in the visual -that also has higher spectral resolution than UVB-and obtained a first fit, including weaker lines.We then added other metal species and adjusted the fitting regions and the number of components whenever necessary to obtain the final fit.
During this process, we also fitted simultaneously some other intervening systems that posses lines blended with that of the proximate systems, or masked the corresponding blended regions.Note also that Zn iiλ2026 is systematically blended with Mg iλ2026, which we did take into account assuming the same components, although Mg i is not a species discussed in this work.When ever possible, we also used Mg iλ1827 and/or Mg iλ2852 to further constrain the (small) contribution of Mg i to the 2026 profile.We also took into account partial blending of Zn ii with Cr ii at 2062 Å.Finally, the analysis of neutral carbon (incl.fine-structure lines), neutral chlo- rine and excited fine-structure lines (from C ii, Si ii and O i) is post-poned to a future work since we focus here on chemical enrichment where their contribution is negligible.That being said, we already note the clear detection of Si ii* towards J0015+1842, J0019−0137, J0059+1124, J0125−0129, J1259+0309, J1358+1410, J2228−0221 and J2325+1539.The fit to the different metal lines is presented in Appendix B.
We note that the metal column densities obtained at medium spectral resolution should always be considered with some caution since the modelling of the velocity profile mostly implied fewer but broader (b frequently as high as 20 km s −1 ) components than typically observed in high-resolution spectra, where b < 10 km s −1 .Therefore, in spite of our best efforts to model the profiles and use lines with a range of strengths, unresolved line saturation may still affect the column densities estimates for some of the species more than the formal uncertainty suggests.Anyway, metallicities are derived using weak Zn ii lines in all but one case, where the total column densities are less sensitive to exact velocity structure decomposition.

Dust
We quantify the amount of dust in the absorbers by measuring the extinction of the quasar, following a template matching procedure as done in many works (e.g., Srianand et al. 2008;Ma et al. 2015;Krogager et al. 2016;Zhang et al. 2015).In order to perform the dust reddening measurement using the full wavelength range provided by X-shooter, we combined the spectra from all three arms to obtain a joint spectra for each quasar.The combination was performed using an inverse-variance weighting in the overlapping regions.During the combination, all three arms were interpolated onto a common wavelength grid using a fixed sampling of 0.2 Å.The joint spectra were lastly corrected for Galactic extinction using the map by Schlafly & Finkbeiner (2011).
The individual joint spectra were then fitted using the quasar template by Selsing et al. (2016) assuming a fixed extinction law at the redshift of the quasar3 .Since we observe no evidence of 2175 Å extinction features, we used the extinction law inferred for the Small Magellanic Cloud (Gordon et al. 2003).In order not to bias the fit we exclude regions of the spectra around the strong, broad emission lines.We furthermore exclude regions that are affected by strong telluric absorption.
For each target, we derive the dust extinction as the average of the A V values for the individual observations.The statistical uncertainty is determined as the standard deviation of the individual measurements.This uncertainty includes any uncertainty on the flux calibration of the different spectra.Lastly, from previous work we assign a systematic uncertainty of 0.07 mag (Noterdaeme et al. 2017) due to any possible intrinsic mismatch of the template and the given target (i.e.accounting for scatter in the intrinsic spectral energy distribution of quasars).The results are summarised in the before-last column of Table 3 and the fits shown in Appendix C.

Absorption line kinematics
We quantify the kinematics of the neutral absorbing gas using the velocity width as defined by Prochaska & Wolfe (1997): , where λ 5 , λ 50 and λ 95 are the wavelengths corresponding to, respectively, the five per cent, fifty and ninety-five percentiles of the apparent optical depth distribution.
For each system, we followed the criteria from Ledoux et al. (2006) and selected a few low-ionisation absorption lines that are not too saturated -it would otherwise be impossible to measure the apparent optical depth-nor too weak -in which case part of the gas would remain untraced-.We then checked for consistency between results.The ∆v 90 for each PDLA is shown in Appendix D and provided in the last column of Table 3.We note that while the smearing of the lines at medium spectral resolution can in principle result in a small overestimation of ∆v 90 (Prochaska et al. 2008;Arabsalmani et al. 2015), this is mostly true for low velocity widths and becomes only an issue when comparing statistics of samples taken at different spectral resolutions.Moreover, the spectral resolution and ∆ v are large enough for the corrections to be here negligible, although, strictly speaking, the smallest values (for J1331+0206 and J1259+0309) could be considered as upper-limits.

Quasar redshifts
A precise measurement of the quasar's redshift is then required to study the motion of the absorbing gas with respect to the central engine along the line of sight, and more generally, to investigate dynamical processes in the circum-nuclear region and/or in the host galaxy.Velocity shifts between the different emission lines are seen in quasar spectra, reflecting the diverse kinematics of the emitting regions (e.g., Gaskell 1982;Tytler & Fan 1992;Vanden Berk et al. 2001;Shen et al. 2016).Overall, it has been found that high-ionisation lines tend to be strongly blue-shifted with respect to the systemic redshift (z sys ), while low-ionisation lines provide a better mean to measure z sys , the most reliable of which being the oxygen lines.
We hence first focused on the [O ii] and [O iii] lines, which are redshifted in the NIR.Unfortunately, telluric absorption bands frequently affect these lines, and in spite of our best efforts to correct for the atmospheric absorption, it is not always possible to recover the unabsorbed emission line properly.In the cases where the lines were reasonably seen, we modelled the region covering Hβ and [O iii] using a simple linear continuum on top of which we added Gaussian lines (with more than one component whenever necessary), but tied the parameters of the two lines of the [O iii] doublet together: same velocity, same width, and a flux ratio of 3. [O ii] was fitted independently.In order to remove spikes, skyline or cosmic-ray residuals as well as absorption lines, we iteratively fitted the data rejecting deviant pixels (at more than 3 σ) at each iteration until convergence (i.e. until no more pixels are rejected).We then visually compared the derived model with the rebinned data, obtained from median averaging in bins of typically 25 pixels.Unfortunately, for two [O ii] lines, sky line residuals are close to the peak and hence the [O ii] redshifts may not be as reliable as expected.This effect is mitigated for [O iii], which is generally stronger and has two lines.
In principle, Mg ii is the following choice at shorter wavelengths (Shen et al. 2016).We note that owing to our selection on the presence of a proximate absorber, the line is affected by strong Mg ii absorption, which we masked during the fitting process.We then also considered the strong C iii] and the weaker He ii lines at 1908 and 1640 Å rest-frame, respectively, which are free from strong absorption.We followed the prescription from Shen et al. (2016) and measure the lines redshifts from the peak of the emission, obtained from a model that reproduces the data.Multiple Gaussians provided convenient models thanks to their flexibility and the relatively small numbers of free parameters.
The redshift obtained this way are given in Table 2 using the same rest-frame wavelengths as in Shen et al. (2016).In the case of [O ii], [O iii] and Mg ii, the authors found small systematic shifts (∼ 8, −45 and −57 km s −1 , respectively) compared to the systemic value, which they determined from Ca ii absorp-tion.This is less than our measurement uncertainties so that we do not attempt any correction here.In the case of the C iii] complex (which includes Si iii] and Al iii), they found a luminosityindependent systematic shift of −229 km s −1 and in the case of He ii the systematic shift depends on luminosity.Here as well, we followed the prescription from Shen et al. (2016) and applied the correction accordingly (see their Eq. 1 and best-fit parameters in their Table 2).We note that these empirical shifts include any possible physical shifts of the emitting regions, plus the fact that several species actually contribute to the complex.We provide also the corrected redshifts in Table 2, in an attempt to remove the average trends seen in quasars.As the average shifts from Shen et al. (2016) are computed based on low redshift quasars observed in the SDSS, we finally did a last test by matching the quasar composite from Selsing et al. (2016) with our data in a region of the visual spectrum covering the range 1620-2450 Å rest-frame, i.e., covering the features from He ii to [Ne iv].This range is slightly reduced in a few cases to avoid bad quality regions, generally at long wavelengths.Absorption lines, bad pixels and sky-line residuals are masked during the fitting process.The only four free parameters are a scaling value, an additional continuum (power-law+constant) and the redshift.
The different measurements (see Appendix E) then allow us to cross-check the emission redshift and decide on the best estimate based on ionised lines (which we call z i syst ), as given in the before-last column of Table 2.In general, we prioritised the lines following the prescription from Shen et al. (2016) although in some cases, we preferred a second-choice line with high S/N over a line with more uncertain measurement.This is the case for J0136+0440 and J0858+1749, where we preferred the corrected C iii] over Mg ii, and J1331+0206 where we preferred Mg ii over [O ii].At the end, He ii and template measurement were not used, but provided for comparison and to give more confidence to the measurement.Finally, in the last column of this table, we indicate the redshift of the CO(3-2) emission line derived from fitting the 3-mm NOEMA data with single Gaussian profiles.This is expected to be an even better indicator of the quasar host systemic redshift, since it traces the bulk of the molecular gas (Wang et al. 2016).The corresponding fits are shown in Appendix F. Interestingly, the CO-based redshifts appear to be systematically higher by about 250 km s −1 (average based on the 6 systems where CO(3-2) measurements are available) with respect to the systemic redshifts based on rest frame UV/optical transitions.This is discussed in Sect.5.3.

Chemical enrichment
In the following, we will discuss the chemical properties of our sample based on the measured column densities and further compare them with other samples of absorption systems: (1) Proximate H i-selected DLAs observed at high spectral resolution by Ellison et al. (2010); (2) a compilation of intervening DLAs (De Cia et al. 2016) with robust abundance measurements and (3) a compilation of H 2 systems from Balashev et al. (2019) and confirmed from follow-up observations.This compilation contains H 2 -bearing systems identified with a number of different selections: in regular DLA systems (Noterdaeme et al. 2008;Balashev et al. 2014Balashev et al. , 2017)), extreme H i absorbers (Noterdaeme et al. 2015;Ranjan et al. 2020) and C i-selected systems (Ledoux et al. 2015;Noterdaeme et al. 2018).A discussion on the selection function is presented in Balashev et al. (2019), so we here only consider the sample as a whole, but keeping in mind that its dis-tribution in parameter space is a sum of several selections.We removed from the intervening samples some systems that have z abs ≈ z QSO and could actually be considered as proximate.In order to avoid possible evolutionary effect in our comparison, we also restrict these samples to z > 2. Finally, we removed known H 2 -systems from the sample by De Cia et al. ( 2016) to avoid duplicates with Balashev et al. (2019).
We use the standard notation of gas-phase abundances of the different species X expressed relative to the Solar values as [X/H] = log(N(X)/N(H tot )) − log(X/H) , and where we take the Solar abundances from Asplund et al. (2021) following the recommendations of Lodders (2003) on whether to take photospheric, meteoritic or average values (see, e.g., table 1 of Konstantopoulou et al. 2022).We do the usual assumption that ionisation corrections are negligible for DLAs thanks the selfshielding from > 13.6 eV photons, i.e., that the considered species are in the main ionisation stage in the neutral gas (e.g., N(Zn ii) = N(Zn)).We do however include H 2 in the total hydrogen column density, i.e., N(H) = N(H i) + 2N(H 2 ).We remind that we consider the gas-phase abundances only when quoting abundances as [X/H].In other words, we do not attempt any correction for dust depletion.We will further investigate this effect by considering volatile and refractory elements in the following sections.A summary of the total abundances of the main species of interest is provided in Table 3.

Overall metallicity from volatile species: zinc and sulphur
Zinc and sulphur are generally considered as the best reference elements for the intrinsic metallicity (i.e., log Z/Z ≈ [Zn/H] ≈ [S/H]) since they are volatile and hence little depleted onto dust grains (see, e.g., Welty et al. 1999;De Cia et al. 2016, among many other works).We confirm that for all PDLAs where both zinc and sulphur abundances can be measured, there is good agreement between both values, with an average measured sulphur-to-zinc abundance ratio [S/Zn] ≈ 0 and dispersion 0.1 dex (see Fig. 2).We note that, in principle, this does not exclude depletion of both species to similar amounts.However, we could expect any depletion to start differing at higher [Zn/H], which is not apparent.We do not see any trend of [S/Zn] with [Zn/Fe] either.In any case, for comparison with other measurements in the literature, we will continue assuming log Z/Z ≡ [Zn/H] in the following.For the only system without Zn measurement (J0136+0440, where the corresponding lines are too heavily blended with atmospheric lines), we hence use Sulphur instead.Overall, we found metallicities in proximate H 2 systems to span a range from log Z/Z ≈ −1.6 to -0.4.This corresponds to the upper half range of H i-selected proximate DLAs, as investigated by Ellison et al. (2010), but similar to the metallicity of intervening H 2 -bearing DLA systems (e.g.Petitjean et al. 2006;Balashev et al. 2019).

Dust depletion from relative abundances
Silicon and iron are refractory elements known to be respectively mildly and strongly depleted onto dust grains and their gas-phase abundances are typically lower than those of zinc and sulphur.In Fig. 3, we show the abundances ratios [Zn/Fe] and [Zn/Si] as a function of metallicity.There is clearly an increase of these ratios with increasing metallicity -as has been observed in a number of works on intervening DLAs (e.g., Ledoux et al. 2003;Noterdaeme et al. 2008;De Cia et al. 2016;Balashev et al. 2019) - (km s −1 ) J0015+1842 20.71±0.1219.73 ± 0.02 -0.41±0.13-0.54±0.14-0.78±0.12-1.79±0.120.40±0.10450 J0019-0137 21.02±0.1120.29 ± 0.01 -0.90±0.12-1.00±0.11-0.96±0.12-2. Notes. () Statistical error only, from dispersion between measurements.The systematic error on A V is about 0.07 mag (see text). (b) Measured from the low-ion metal lines (see text). (c) Balashev et al. (2019) quoted A V = 0.16 and 0.29 for J0136+0440 and J0858+1749, respectively.However, these values were obtained with a varying intrinsic power law, as usually done in attempt to separate intrinsic slope from that due to an intervening absorber.Here, as for other quasars in our sample, we consistently re-measured A V using a fixed power-law to quantify amount of dust in the overall systems.
which can be attributed to the increasing abundance of dust with increasing metallicity.
Iron appears to be significantly depleted even down to the lowest metallicities while [Si/Zn] becomes positive only for metallicities higher than 1/10 th Solar.Interestingly, [Zn/Fe] in proximate H 2 systems are clearly higher for a given metallicity than the values seen in regular intervening DLAs suggesting a higher depletion of iron in proximate systems.This remains true when comparing only to intervening H 2 absorbers, even if the latter tend to have high depletion among the overall (mostly non-H 2 bearing) intervening population at a given metallicity.It is also interesting to see a lower dispersion for proximate systems around the main trend (dotted blue line), likely suggesting a more homogeneous population.
On the other hand, the difference between proximate H 2 and intervening systems (regardless of their H 2 content) is less clear in the case of [Zn/Si].Since silicon is an α-element predominantly produced in massive stars, it is possible that a stronger intrinsic depletion (with respect to intervening systems) is com-pensated by recent star-formation, or that the depletion sequence for Si and Fe also differs.
In short, the high depletion (from [Zn/Fe]) at a given metallicity suggests a mechanism for build-up and/or destruction of dust grains that is similar across the population of proximate H 2 systems but different from that seen in intervening DLAs.Any correction to account for some depletion of zinc and recover the intrinsic metallicity could then also be different for proximate systems.We note that assuming the same depletion sequences as in intervening systems, the corrected metallicities, following De Cia et al. ( 2016) would be a factor two higher, on average.
The abundance of dust can be independently quantified using the extinction per hydrogen atom, A V /N(H).While the 0.07 mag uncertainty on A V measurements remains quite large compared to their values, there is a clear tendency for higher A V /N(H) at higher [Zn/Fe], [Zn/Si] and metallicities, strengthening a dustdepletion origin for the observed underabundance of iron and silicon.The trend between depletion and value of A V /N(H) turns out to be a strong correlation with Pearson coefficient r = 0.9 and probability of chance coincidence 0.07%, see    (2016).These values appear to be in relatively good agreement with the reddening-based measurements, albeit the former being on average higher by 0.25 mag than the latter.This may be another indication for a different type of dust close to the quasar.However, a comparison of extinction derived from gas-phase abundances with that directly obtained from the spectral slope should be performed first in intervening systems to validate the prescription by De Cia et al. ( 2016) at such precision level.

The sulphur-to-silicon ratio in proximate absorption systems
The [S/Si] ratio is also an interesting quantity to look at since both species are α elements and their intrinsic abundances are expected to follow each other.According to  H i. This can be due to a bias against systems with high column densities of dust that would be either excluded from the parent quasar sample (Ledoux et al. 2009;Augustin et al. 2018) or result in low S/N ratio of the corresponding spectra.It can also be due to a higher abundance of metals/dust required for the presence of H 2 at low H i column densities (e.g., Bialy et al. 2017), thereby increasing the dust shielding, H 2 formation rate and gas cooling through atomic lines.In the low N(H i) regime, we would qualify the observed [S/Si] in (non-H 2 ) PDLAs as highly dispersed and possibly interpret this as a wide range of ionisation.
To investigate the effect of differential dust depletion further, we plot [S/Si] as a function of [Zn/Fe] in Fig. 6.Here, we see an apparent correlation between the two abundance ratio in the case of proximate H 2 systems (Pearson correlation coefficient 0.76 with 0.4% probability to be due to chance coincidence).This favours depletion of silicon as a more likely origin for the high [S/Si] values.2016) together with the associated intrinsic dispersion4 .However, proximate systems differ from intervening DLAs in the sense that they have [S/Si]≈ 0 for a wider range of [Zn/Fe] values.This is observed for PDLAs regardless of whether they contain H 2 or not, when the intercept for intervening DLAs is almost zero.
In summary, the PH 2 systems with high [S/Si] values appear to behave like intervening DLAs.The five systems with highest values also have low H i column densities, which is naturally explained if high-H i high-depletion systems are missing due to dust bias.At low [Zn/Fe] (and high N(H i)), the [S/Si] values tend to differ between intervening and proximate systems, the latter having typically lower observed [S/Si] -from S ii and Si iicloser to solar.We conclude that the main driver of the [S/Si] ratio is depletion onto dust -with possibly a different depletion sequence than that seen for intervening DLAs, again suggesting different dust type proximate to the quasar.However, ionisation effects are still possibly at play, meaning that part of these lowionisation metals could actually arise from ionised gas.Larger samples of proximate DLAs (with and without H 2 ) are necessary to investigate this further.

Argon as probe of hard UV ionising field
Argon can provide an additional diagnostic on the ionising field.Owing to its high first-ionisation energy (15.76 eV), argon is expected to be mostly neutral in DLAs, which are selfshielded against H i-ionising photons.While argon is predominantly ionised by cosmic rays in local ISM, high-energy pho-  (Lodders 2003), the argon deficit is likely dominated by ionisation from the quasar metagalactic radiation, with local modulation by H i column inside the absorbing galaxies.We constrain the Ar i column densities in our sample of proximate H 2 absorbers using the two available absorption lines (at 1048 and 1066 Å).Because these fall into the Lyman-α forest and can easily be blended with H 2 /HD lines, we used a slightly different procedure: we fitted the Ar i lines by fixing the redshifts and Doppler parameters to that previously obtained for other metal metal lines and fitting for the column densities.We included the H 2 absorption profile constrained by independent fit (see Section 3.1) and avoided blends with Ly-α forest as best as possible.In a few cases, however, we could get no meaningful constraint.The fit to argon lines are shown in Fig. B.10.
The derived [Ar/S] (obtained from Ar i and S ii) values are shown as a function of the total H i column densities and overall metallicities in Fig. 7. First, we do not see any dependence of [Ar/S] on N(H i), in agreement with the expected lack of shielding by neutral hydrogen.We do not see any dependence with overall metallicity either, contrary to the possible negative (but not easily explained) trend seen by Zafar et al. (2014b).The lack of dependence on overall metallicity is in line with the expected non-depletion of argon onto dust.Second, [Ar/S] is found to be very similar (-0.35 dex on average) to that seen in intervening DLAs, which could be surprising given the proximity of the quasar.However, our measurements can still be affected by Ly-α 0.00 0.25 0.50 0.75 1.00 1.25 Symbols and colours for data points are as in Fig. 5.
interlopers so that the Ar i columns should in principle be considered as upper-limits -in a conservative approach5 -but more importantly, the metal column density is generally dominated by the H 2 -bearing component, where we observe lower deficit of argon (-0.2 dex) with respect to components without detectable H 2 (-0.7 dex).In Fig. B.10, we also show the expected Ar i profiles obtained assuming Solar [Ar/S] ratio.This over-predicts the absorption in most cases, but we also see that such over-prediction seems less severe at the velocity of H 2 components.
The relatively low deficit of argon in H 2 components most likely is a result from an efficient conversion of Ar + into neutral argon through chemical reactions with H 2 (Ar + + H 2 → ArH + + H followed by ArH + + H 2 → Ar + H + 2 , see Duley 1980;Schilke et al. 2014;Neufeld & Wolfire 2016).The situation could be similar to that with neutral chlorine (Balashev et al. 2015) and all argon being in neutral form in the presence of H 2 .
In turn, in the regions where hydrogen is predominantly atomic, the ionisation fraction of argon should be very dependent on the ionising field.Indeed, at higher I UV the photo-ionisation rate increases, while the H 2 abundance decreases, impeding the molecular-ion recombination route.Hence the ionisation fraction of argon should rise at least proportionally to I UV .Therefore, the systematically lower abundance in non-H 2 -bearing components confirms that the medium is likely exposed to significantly enhanced (and hard) radiation field and/or cosmic ray ionisation.While quantitative statements would require a full modelling of the observed column densities and a better-assessed association between H 2 and metal components (i.e. at high spectral resolution) -out of the scope of this paper-we can estimate that the  et al. 2014a).In spite of our best efforts to take into account for blends, all points should in principle be considered as upper-limits.
−0.2 dex deficit of argon means that about half of the metal column seen at the velocity of H 2 components could come from their (H 2 -free) H i-envelope where argon would be mostly in ionised form.

H 2 in proximate versus intervening systems
We compare proximate and intervening systems according to their H 2 and total hydrogen column densities in Fig. 8.The absence of proximate systems at low H 2 columns is naturally due to our selection based on the presence of strong H 2 Lyman-Werner bands (Paper I), and in fact, all our systems have log N(H 2 ) > 18.5.In turn, proximate H 2 were selected independently of N(H i), but as can be seen in this figure, they are limited to log N(H) > 20.5, when intervening systems can have lower N(H) even at similar H 2 column densities (see distributions on top of this figure).However, this is driven by a few very metalrich intervening absorbers (as can be seen from the colour coding), two of them following a C i-based selection, skewed towards very high metallicities.Actually metallicities tend to increase with decreasing N(H) and increasing N(H 2 ).This could be naturally explained by the 1/Z dependence of the critical column density for H 2 to form (e.g., Krumholz et al. 2008), but also to a bias against high-metallicity systems at high-N(H i) due to stronger dust extinction.This is better seen in Fig. 9, where the location H 2 systems is shown in the N(H)-Z plane.The different lines show the ratio of UV intensity to number density as a function of the H i column density in the envelope of a single H 2 cloud illuminated on one side 6 .Three quarters (9/12) of the proximate H 2 systems are located between the I UV /n = 10 −2 and 10 −1 cm 3 lines, suggesting similar physical conditions.We caution however that the measured H i column correspond to that integrated over the profiles while the calculations provide the H i-column that participate in the shielding of H 2 clouds.Notwithstanding, any atomic gaswith some amount of dust-located upstream along the line of sight will still participate in shielding H 2 from the quasar field.
Determining the physical conditions from population levels is thus necessary to constrain the distances of the cloud from the AGNs.In turn, owing to the different selections, intervening H 2 systems tend to have a wider spread of values.We also note that the lines of constant I UV /n closely follow constant metal column density (for volatile species), and hence almost constant extinction (e.g., Vladilo et al. 2006).In other words, the conditions that facilitate the presence of H 2 are also those that can extinguish the background sources to the point that may be excluded from the optically-selected samples.However, since most of the intervening systems were also detected towards quasars from the SDSS, it is likely that more reddened quasars with more metalrich proximate H 2 absorbers are actually already present in the SDSS database -in particular since such systems could afford lower densities to survive-but the S/N achieved in the blue prevented their detection directly from H 2 lines.
We also show in Fig. 9 the location of proximate and intervening DLAs without strong H 2 as open symbols.In the case of intervening non-H 2 systems, we used here the sample from Noterdaeme et al. (2008) where the presence of H 2 has been systematically assessed.The fact that no H 2 is detected, even at high metallicity and H i column density means that the corresponding system do not fulfil the requirements on I UV /n.For example, for a UV field about 1/10 th of the Draine field, such systems would have low number densities n < 1 cm −3 (which is typical for the warm neutral phase), too low to build up molecular hydrogen.Alternatively, their total H i column density is spread over many different clouds.We finally note that the above implicitly assumes static equilibrium, which may not necessarily be the case, in particular if H 2 is found in outflowing gas.Indeed, for n = 100 cm −3 the characteristic timescale for H 2 formation will be ∼ 10 7 yrs, which is then comparable to the dynamical timescales.

Kinematics
The velocity extent of low-ionisation metal lines, ∆v 90 is considered as a good way to quantify the kinematics of the neutral gas (since this information is not available from the H i lines).In the case of intervening DLAs, Ledoux et al. (2006) have first observed a correlation between this velocity width and the gas metallicity, suggesting that this reflects an underlying relation between mass and metallicity, since the motion of low-ionisation lines is expected to be governed by gravity.This has been discussed in a number of subsequent works (Møller et al. 2013;Neeleman et al. 2013;Christensen et al. 2019), with attempts to connect this to the mass-metallicity relation seen in emission 6 Following Eq. 2 of Paper I, based on theoretical calculations by Bialy et al. (2017) and where we used Z/Z as a proxy for σg , the grain Lyman-Werner photon absorption cross section per hydrogen nucleon normalised to the fiducial Galactic value.Note that such calculations do not provide the column density of H 2 which just continues building up as more gas would be added to a given cloud.In practice, the total column is dominated by H i, so that the plot against N(H i)-only is almost identical.(e.g., Møller & Christensen 2020).While the assumption that ∆v 90 depends primarily on the galaxy mass seem to be consolidated at least in a statistical sense, there are evidences that outflows, tidal-streams and environment in general can affect the velocity spread in a more complex way (e.g., Zou et al. 2018;Nielsen et al. 2022).In the case of proximate systems, because of the presence of the quasar, one could expect the situation to be even more complex.For example, the case of J0015+1842 is strongly suggestive of outflowing material (Paper II).Surprisingly enough, the distribution of ∆v 90 in proximate H 2 -systems follows the trend and dispersion seen in intervening DLAs, although limited to high-values, as expected from the presence of H 2 (Noterdaeme et al. 2008).Overall, the velocity spread of the neutral gas in proximate-H 2 systems suggests an origin in massive galaxies, but without clear indication of perturbed kinematics, or at least, not more perturbed than similar high-metallicity systems in the intervening population.
In order to investigate now the kinematics of the gas with respect to the quasar, we represent in Fig. 11 the kinematics of the neutral gas as a function of its metallicity.There is an apparent trend for high-metallicity systems to be around the quasar systemic redshift as obtained from the ionised lines, while lowmetallicity systems tend to be redshifted.The Pearson correlation coefficient between v 50 = c(z QSO − z 50 )/z QSO (i.e., the velocity of the centroid of the absorption profile ) and the gas metallicity is r = −0.50 with p-value 0.08.In other words, there is  only about 8% probability of this anti-correlation to be due to chance alone.One system, towards J0858+1749, is significantly more blueshifted than the rest of the sample.This does not drive the correlation however as removing it from the correlation test only changes the Pearson values to r = −0.49,p = 0.10.In fact, the ∼ 2000 km s −1 blueshift of this system makes it intermediate between what could be considered as clearly intervening or proximate.For example, at this velocity, the incidence of strong H 2 absorbers already reaches that of intervening ones (see Fig. 3 Noterdaeme et al. 2019).In fact, the derived UV flux in this system (Balashev et al. 2019) suggests that this system is located further from the quasar, at a distance of at least several hundreds kpc.To account for the statistical uncertainties on metallicity and on the quasar systemic redshift, we calculated the Pearson's (r,p) values for 10 000 random samples where the metallicity of each system is taken using a normal distribution around the best value with the corresponding uncertainties.Similarly, the velocity is chosen using a normal distribution around the measured value with σ set to the intrinsic velocity dispersion of z QSO provided by Shen et al. (2016) (i.e., 46, 56, 205 and 233 Mg ii and C iii], respectively).From this exercise, we find an approximate normal distribution of r values, giving r = −0.47±0.09(r = −0.44±0.11ignoring J0858+1749).
The p-values have a median of 0.11 (0.14) but with some tail towards larger values.In short, this means that the trend remains even taking into account statistical uncertainties on metallicity and emission redshift measurements.
In principle, the projected line-of-sight velocity towards the quasar for low-metallicity systems could be interpreted as neutral gas infalling onto the quasar host galaxy.However, it is unlikely that gas would feed directly the quasar and we may instead  Ledoux et al. (2006).
expect any feeding gas to orbit the quasar host galaxy without significant line-of-sight velocity component.It is possible however, that gas is orbiting the galaxy, in non-circular orbits, after a galaxy companion interaction, where outer neutral gas is coming back as a fountain after having been dragged out in a tidal tail.Alternatively, the systemic redshift could be underestimated if the blueshifts of the emission lines are larger than observed by Shen et al. (2016).In their work, the authors used Ca ii absorption as reference to compute the (non)-shifts of [O ii] and [O iii] and then the shift from other lines (Mg ii, C iii], etc.) are computed with respect to the latter.The use of SDSS data (covering wavelengths shorter than 10 400 Å) implies that the computed shifts correspond to low to intermediate redshift quasars (and hence likely not particularly luminuous).At high redshifts (z ∼ 5 − 7), Schindler et al. (2020) and Eilers et al. (2021) found larger Mg ii blueshifts in luminuous quasars, using precise systemic redshifts from [C ii]µm158 and/or CO emission line measurements from ALMA and NOEMA as reference.These lines arise respectively from star-formation and molecular gas, and hence should better represent the systemic redshift of the quasar host, which is also very likely to dominate over other, less massive galaxies in the quasar group.The sample presented here has redshifts z ∼ 3, so we may naively expect line shifts in-between those obtained by Shen et al. (2016) and those obtained by Schindler et al. (2020) or Eilers et al. (2021).For half of our sample, we do have CO(3-2) redshifts from our NOEMA observations, which we represent as ellipses on Fig. 11.On average, the CO redshifts are about 250 km s −1 higher than those obtained from the (corrected) ionised emission lines.If z CO does represent better estimates of the host galaxy redshift, then the absorbing gas now appears to have no particular line-of-sight velocity at low-metallicities, while the outwards velocities at high metallicities is reinforced.This maintains the anti-correlation to r ≈ −0.6, p ≈ 0.2 (now computed on 6 quasars).In that case, high-metallicity gas (above ∼25% Solar) is more likely to represent gas enriched by intense star-formation activity and being expelled through associated winds, as clearly seen in the case of J0015+1842 (Paper II), which has the highest metallicity among our sample.This would be consistent with chemically-enriched outflowing gas coming from nuclear regions, contrary to outer gas, due to the radial Z-gradient.
Determining the distances between the absorbing gas and the active nucleus by modelling the excitation of molecular and atomic species should help understanding the origin of the gas with different metallicities further.Such a study using the present sample is devoted to a future paper.

Note on the AGN proximity from the presence of N v
While we focus on the neutral phase (and low ionization metal species) in this paper, we briefly discuss here the presence of N v absorption at the quasar redshift.Photo-ionisation of nitrogen to that stage requires an energy of 77.5 eV, meaning a hard ionisation spectrum since stellar light falls rapidly after 54 eV where He ii ionisation occurs.This makes N v one of the best diagnostic for the physical proximity of the AGN and distinguish between intervening and associated system (e.g., Perrotta et al. 2016).Fox et al. (2009) found that N v is detected in only 13% of intervening DLAs but they also found a surprisingly similar fraction in PDLAs out to 5000 km s −1 , suggesting that proximity in velocity does not necessarily mean physical proximity.Indeed, 5000 km s −1 in the Hubble flow corresponds already to about 16 Mpc proper distance at z ∼ 3. Ellison et al. (2010) found a tentative excess of N v detection rate (2/7) at similar velocities in their PDLA sample, but note that this fraction would increase if they limited the statistics to smaller velocity separations from the quasar.Finally, using composite spectra of about a thousand C iv absorbers and interpreting their results with the help of ionisation modelling, Perrotta et al. (2018) confirmed that N v is only present in proximate systems and conclude it is a good tool to identify gas physically related to the quasar host.N v absorption is clearly seen at the quasar redshift for nine out of our thirteen proximate H 2 -absorbers (see Fig. A.1).Note that the non-detection of N v in the two systems with lowest metallicity in our sample may just be due to column density effect rather than lower ionisation.Overall, the very high incidence of N v in our sample supports an origin in the close environment of the quasars.Careful ionisation modelling, together with comparative kinematic (N v is typically present over wider velocity ranges than low-ionisation species) should help understanding the ionised component of H 2 -selected proximate systems.

Conclusion
The chemical enrichment and velocity dispersion of the neutral gas in strong proximate H 2 systems appears not to be significantly different from what is seen in intervening systems.While the metallicities in proximate H 2 systems are higher than the general population of intervening DLAs, and in the upper-half when compared to the sample of proximate DLAs from Ellison et al. (2010), they are consistent with those intervening systems that also have high H 2 content (see Balashev et al. 2019).The relative abundances of sulphur, silicon, iron and zinc follow the expected trends due to increasing depletion of refractory elements with increasing metallicity (e.g., De Cia et al. 2016), indicating an increased abundance of dust with increasing chemical enrichment.This is also confirmed independently by the increasing Ellipses are shown at the location of the CO(3-2) emission whenever observed, with major axis equal to the line's FWHM.z CO is shifted with respect to z i syst by about 250 km s −1 on average (dash-dotted line).
quasar extinction per hydrogen atom with depletion and metallicity.We however note higher depletion of iron in proximate systems compared to intervening ones at given metallicity, suggesting different dust production and/or destruction close to the AGN.The ubiquitous presence of N v absorption indicates indeed a strong and hard UV field, which also explains the deficit of neutral argon in non-H 2 -bearing components.In turn, argon remains predominantly in neutral form in the the presence of H 2 , aided by ion-molecule reactions.
The lack of proximate and intervening systems with both high-metallicity/depletion and hydrogen column density is likely due to a selection effect against highly reddened systems which would be more severe than just exclusion from optically-selected quasar samples.Actually, such systems may already be present in the SDSS quasar database, but not selected through H 2absorption because of their low S/N in the blue.Those systems would also be interesting to recover since they should be able to survive with relaxed densities for a given intensity of UV field (i.e., distance to the central engine).Upcoming surveys should also provide less biased samples of quasars and absorbers.For example, 4MOST (de Jong et al. 2019) will collect optical spectra of quasars selected on astrometry (The 4MOST-Gaia Purely Astrometric Quasar Sample -4G-PAQS; PI: Krogager), infrared or X-ray properties (Merloni et al. 2019).The higher spectral resolution and larger collecting area should push deeper into the dust-obscured population of absorbers.Similarly, the MeerKAT Absorption Line Survey (MALS; Gupta et al. 2016) also uncovers quasars in a dust-unbiased manner, owing to radio+IR selection of the targets (see Krogager et al. 2018;Gupta et al. 2022a).A complementary view of proximate cold gas absorbers can then be done through associated 21-cm absorption.Further, 21-cm absorption provides different (and multiple) lines of sight owing to radio-continuum emission regions distinct from that producing the optical continuum (Gupta et al. 2022b).
The significant excess of proximate H 2 systems with respect to the expected incidence from the intervening statistics showed that most of the former are very likely to be associated to the quasar environment (quasar host or galaxies in the quasar group).On the other hand, the similar column densities of neutral and molecular hydrogen at similar metallicities when compared to intervening H 2 -bearing DLAs means that proximate H 2 should have similar UV-to-density ratios.This means that proximate H 2 -systems should have significantly higher density than their intervening analogues if located within the dominating influence of the quasar UV field (i.e., within the host galaxy, or possibly several hundreds of kpc away).The higher densities if located close to the quasars would explain the higher fraction of systems with excited atomic fine-structure lines.Indeed, Si ii* is clearly detected in eight out of the thirteen systems presented here.Determining the physical conditions will be necessary in that respect too.
While our parent sample was selected by having z abs ≈ z QSO , the systems could potentially have been found within a wide velocity range around the quasar owing to the large uncertainties on the quasar emission redshifts in SDSS spectra, together with the low S/N and resolution allowing H 2 -detection over a relatively wide range.It is then almost surprising to find that all systems but one ended up having velocities significantly smaller than 1000 km s −1 from the quasar, after securing precise absorption and emission redshifts.In the absence of peculiar motions (i.e., pure Hubble flow), this would correspond to at most a few Mpc.If peculiar motions contribute significantly to the observed velocities, then this becomes an upper-limit to the actual distance.In fact, velocities of z abs > z QSO systems gives a hint on the dispersion due to peculiar velocities since this obviously cannot be due to the Hubble flow.The kinematics are then all consistent with no Hubble-flow component at all (except for J0858+1749), suggesting that most systems are indeed associated to the quasar host or its group environment.
We observe a trend for different line of sight velocities of the absorbing gas with respect to the quasar as a function of metallicity.Assuming systemic redshifts from ionised lines, this would suggest that low-metallicity gas tends to move towards the quasar.However, it is more likely that the low-metallicity gas is actually more clustered (in velocity space) around the quasar, as corroborated by the use of CO-based redshifts, and that the highmetallicity gas is actually more likely to arise from outflowing winds.While this remains speculative for now, if confirmed in a large sample of systems, including non-H 2 ones, the relation between motion of neutral gas with respect to the quasar host galaxy and/or the central engine and the chemical enrichment of this gas would open a new window to understand quasar evolution.
Selection effects seem however to be complex when it comes to proximate systems.For example, the classical selection of PD-LAs based on saturated cores biases the samples against systems with leaking Lyman-α emission, that is actually seen in about half of our proximate H 2 systems.If, as expected, the amount of gas and its distribution around the quasar affects the ability of Lyman-α photons to scatter at large distances, such a bias will also affect any conclusion on the nature of proximate systems.
We conclude that proximate systems have a strong potential to investigate many processes driving the evolution of quasars.Determining the physical conditions and distances of gas clouds from the central engine will be the next natural step on this matter.

Fig. 1 .
Fig. 1.Comparison of H 2 (purple diamonds) and H i (green hexagons) column densities obtained using X-shooter and SDSS spectra.The dashed line depicts the one-to-one relation, and the dotted lines show ±0.3 dex around it.

Fig. 2 .
Fig. 2. Observed abundance of sulphur against that of zinc in proximate H 2 systems.

Fig. 4 .
Fig. 4. The low dispersion around the approximate linear relation (log A V /N(H) ≈ 1.3[Zn/Fe] − 23.5) indicates that the quasar template matching provides reliable estimates, at least relatively to each other.Interestingly, it is also possible to indirectly derive the extinction from the gas phase abundances, assuming some scaling with respect to the Galactic values, which we also represent in Fig. 4, based on calculations following De Cia et al.(2016).These values appear to be in relatively good agreement with the reddening-based measurements, albeit the former being on average higher by 0.25 mag than the latter.This may be another indication for a different type of dust close to the quasar.However, a comparison of extinction derived from gas-phase abundances with that directly obtained from the spectral slope should be performed first in intervening systems to validate the prescription by De Cia et al. (2016) at such precision level.
Fig.3.Zinc-to-iron (bottom) and zinc-to-silicon (top) for our sample of proximate H 2 systems (circles with error bars).The colour-scale provide the value of A V /N(H), which is proxy of the dust-to-gas ratio.For two systems the derived A V values are actually negative (hence non-physical) and the corresponding points are unfilled.Green triangles represent measurements in intervening DLAs(De Cia et al. 2016) with filled ones corresponding to known H 2 intervening systems(Balashev et al. 2019).The solid blue lines show the best linear fit to the value for the proximate H 2 systems, taking into account uncertainties on both axis.The green lines and shadow region show the best linear fit and associated 1 σ uncertainty for the overall sample of intervening DLAs.

]Fig. 4 .
Fig. 4. Dust-to-gas ratio measured from the continuum spectrum (filled circles) and derived from the metal column densities (open diamonds, following De Cia et al. 2016), as a function of depletion factor.The blue (resp.purple) dotted line corresponds to a linear fit to ([Zn/Fe], log A V /N(H)) data points, where A v is based on continuum fitting (resp.metal column density).

Fig. 7 .
Fig.7.Observed [Ar/S] ratio (based on Ar i and S ii) as a function of overall metallicity (top) and H i column density (bottom) in our sample of proximate systems.Blue circles correspond to the overall value in each system, while purple diamonds and red crosses to the H 2 and non-H 2 components, respectively.The dashed lines correspond to the corresponding mean values, with the distribution shown on the side panel.The green line corresponds to the mean [Ar/S] value (corrected to same Solar reference) and dispersion for intervening DLAs(Zafar et al. 2014a).In spite of our best efforts to take into account for blends, all points should in principle be considered as upper-limits.

Fig. 8 .
Fig. 8. Column density of molecular hydrogen versus total neutral hydrogen column density N(H) = N(H i) + 2N(H 2 ) for our sample of proximate H 2 absorbers and z > 2 intervening strong H 2 -bearing DLAs fromBalashev et al. (2019).For a fair comparison, we removed systems that may potentially be proximate, keeping only systems with a conservative cut of c (z abs − z QSO )/(1 + z QSO ) < −10 000 km s −1 .Our selection of proximate H 2 systems based on strong H 2 lines directly seen in the low-resolution SDSS spectra naturally skews our sample to high H 2 columns.The blue horizontal line shows the approximate and empirical lower-limit of at log N(H 2 ) > 18.5 for the present sample.The top histograms show the N(H)-distributions for the proximate (blue) and intervening (green) systems with N(H 2 ) > 18.5.Finally, the different lines correspond to different averaged molecular fractions f H2 ≡ 2N(H 2 )/(2N(H 2 ) + N(H i)).

Fig. 11 .
Fig. 11.Summary of absorbing gas metallicity and kinematics.The different coloured segments represent the ∆v 90 regions, with dots indicating the location of H 2 components.The zero of the velocity scale is set to the best estimate of the systemic redshift z sys from ionised lines (Sect.4.2).Positive velocities mean here z > z i sys .Ellipses are shown at the location of the CO(3-2) emission whenever observed, with major axis equal to the line's FWHM.z CO is shifted with respect to z i syst by about 250 km s −1 on average (dash-dotted line).

FigFigFig
Fig. A.2. Portion of X-shooter spectrum of J0019−0137 featuring H 2 absorption bands.The solid red line shows the overall absorption model, while contribution from non-H 2 lines (H i, metals and lines from intervening systems) is shown as dashed profile.The blue segments connect lines from different rotational levels (shown here for convenience from J=0 to J=5) of a given band (indicated by the label above, 'L' for Lyman, 'W' for Werner and the following number being the upper vibrational level.

Table 1 .
Quasar sample and log of X-shooter observations

Table 2 .
Measurements of the quasar emission lines redshift.

Table 3 .
Main properties of the proximate absorption systems Rix et al. (2007)010), based on modelling byRix et al. (2007), the observed negative [S/Si] values at low H i column densities in their sample of PD-LAs is due to ionisation by a hard spectrum.We compare this ratio in our sample of proximate systems with those from Ellison et al. as a function of the H i column density in Fig.5.Contrary to what these authors observed in PDLAs, we do not see any drop in [S/Si] below log N(H i) = 20.75 in our strong H 2 -bearing PDLAs.While we do not probe H i columns as small as these authors do, it is interesting to note that both PDLAs and PH 2 have close to Solar ratios at high column densities log N(H i) > 21.In contrast with the PDLAs from Ellison et al., we see a tendency for increasing [S/Si] with decreasing H i-column, which is opposite to what is expected if such ionisation corrections are at play but rather reflects an increasing depletion with decreasing Ellison et al. (2010)licon abundance ratio versus H i column density.The horizontal dashed line marks the Solar ratio while the vertical line shows the N(H i) limit below which Ellison et al. (2010) observed under-Solar abundances in PDLAs (see text).Note that, while not in our sample of proximate H 2 , one PDLA fromEllison et al. (2010), namely J2340-00, also has strong H 2 and is marked with a blue circle.
Zafar et al. (2014b))ove 13.6 eV can still penetrate the cloud and ionise argon, that further has a low recombination to ionisation rate in diffuse atomic medium (see e.g.Sofia & Jenkins 1998).Zafar et al. (2014b)found a deficiency of argon by about 0.4 dex (with standard deviation about 0.2 dex) with respect to other αelements in intervening DLAs, with a small dependency on H i column density.They concluded that, since argon is extremely volatile Fig. 6.Sulphur-to-silicon vs zinc-to-iron abundance ratios in proximate H 2 and DLA systems.The green line and shaded area represent the best-fit for intervening DLAs by De Cia et al. (2016) and the 1 σ scatter around this relation.The dashed blue line is a linear fit to the proximate H 2 data.Dashed vertical and horizontal lines mark Solar values.
0 0 Metallicity vs total hydrogen column density.The different lines correspond to the hydrogen column density of a single H 2 cloud illuminated by a UV field for different UV to density ratios I UV /n as indicated above each line, with I UV in units of the Draine field and n in cm −3 .
Filled symbols correspond to H 2 -bearing systems.Proximate H 2 systems are from this work, proximate DLAs from Ellison et al. (2010) and intervening systems from Noterdaeme et al. (2008) and Balashev et al. (2019).
Fig. 10.Metallicity versus velocity extent of the low-ionisation metal lines for our sample proximate H 2 -systems and a sample of intervening DLAs (mostly from Ledoux et al. 2006 with some additions by Noterdaeme et al. 2008).The green line shows the best-fit relation for high-z intervening DLAs by