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

TL;DR

This study investigates the chemical enrichment and kinematics of proximate molecular quasar absorbers using multi-wavelength spectroscopy, revealing insights into metallicity, dust, and gas dynamics near quasars.

Contribution

It provides the first detailed analysis of chemical and kinematic properties of proximate molecular absorbers, highlighting differences from intervening systems and the role of the quasar environment.

Findings

Gas-phase metallicities range from 2% to 40% of Solar.

High-ionisation phase traced by NV is almost ubiquitous.

No strong deficit of neutral argon despite hard UV fields.

Abstract

Proximate molecular quasar absorbers (PH2) are an intriguing population of absorption systems recently uncovered through strong H2 absorption at small velocity separation from the background quasars. We performed a multi-wavelength spectroscopic follow-up of thirteen such systems with VLT/X-Shooter. Here, we present the observations and study the overall chemical enrichment measured from the HI, H2 and metal lines. We combine this with an investigation of the neutral gas kinematics with respect to the quasar host. We find gas-phase metallicities in the range 2% to 40% of the Solar value, i.e. in the upper-half range of HI-selected proximate damped Lyman-alpha systems, but similar to what is seen in intervening H2-bearing systems. This is likely driven by similar selection effects that play against the detection of most metal and molecular rich systems in absorption. Differences are however seen in the abundance of dust (from [Zn/Fe]) and its depletion pattern, when compared to intervening systems, possibly indicating different dust production or destruction close to the AGN. We also note the almost-ubiquitous presence of a high-ionisation phase traced by NV in proximate systems. In spite of the hard UV field from the quasars, we found no strong overall deficit of neutral argon, at least when compared to intervening DLAs. This likely results from argon being mostly neutral in the H2 phase, which actually accounts for a large fraction of the total amount of metals. We measure the quasar systemic redshifts through emission lines from both ionised gas and CO(3-2) emission, the latter being detected in all 6 cases for which we obtained 3-mm data from complementary NOEMA observations. For the first time, we observe a trend between line-of-sight velocity with respect to systemic redshift and metallicity of the absorbing gas. [truncated]

Tables

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

Quasar (SDSS)a𝑎aa𝑎a$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.	$z$	Observing dates	Seeingb𝑏bb𝑏b$b$Delivered seeing corrected for airmass, as measured at the telescope by the Shack–Hartmann wavefront sensor (fits keyword HIERARCH ESO TEL IA FWHM). (arcsec)	airmass	Ref.
J001514.82+184212.3	2.63	29/07, 31/08, 31/08, 28/09 2019	1.67, 0.75, 0.65, 0.53	1.48, 1.39, 1.41, 1.42	(1)
J001930.55-013708.4	2.53	06/11, 07/11, 08/11 2021	0.91, 0.72, 0.71	1.30, 1.26, 1.16	(2)
J005917.64+112407.7	3.03	09/11, 09/11, 09/11 2021	0.77, 0.80, 1.14	1.24, 1.30, 1.47	(2)
J012555.11-012925.0	2.66	21/07, 01/09, 01/09, 11/10, 25/11 2021	1.14, 1.00, 0.98, 1.22, 0.92	1.11, 1.10, 1.17, 1.10, 1.17	(3)
J013644.02+044039.1	2.78	27/11, 28/11 2014	1.15, 1.15	1.16, 1.26	(4)
J085859.67+174925.1	2.65	24/12 2014, 15/02 2015	0.73, 1.38	1.36, 1.37	(4)
J123602.11+001024.5	3.03	07/02, 07/02, 07/02 2022	0.73, 0.83,0.87	1.20,1.12,1.11	(2)
J124829.51+063925.6	2.53	08/02, 07/03, 07/03 2022	1.21, 0.61,0.65	1.19, 1.18, 1.19	(2)
J125917.31+030922.5	3.23	06/02, 06/02 2022	0.79, 0.93	1.33, 1.43	(2)
J133111.41+020609.0	2.92	08/03, 29/05, 29/05 2022	0.80, 0.84, 0.98	1.25,1.17,1.30	(2)
J135808.94+141053.2	2.89	25/02, 01/03, 07/03 2022	1.01, 0.79 , 0.68	1.29, 1.29, 1.30	(2)
J222807.36-022117.1	2.77	31/08, 15/09, 12/10, 26/10 2021	0.87, 0.71, 0.84, 0.85	1.10, 1.12, 1.09, 1.09	(2)
J232506.62+153929.3	2.62	29/08, 31/08, 30/09 2019	0.84, 0.77, 0.67	1.33, 1.33, 1.34	(2)

Table 2. Table 2: Measurements of the quasar emission lines redshift.

Quasar	Redshift
	[O ii]	[O iii]	Mg ii	C iii]		He ii		Template	$z_{syst}^i$	CO(3-2)
				meas.	corr.	meas.	corr.
J0015$+$1842	2.6310	2.6309	–	–	–	–	–	–	2.6310	2.6285
J0019$-$0137	–	2.5217	2.5276	2.5185	2.5212	–	–	2.5253	2.5217	–
J0059$+$1124	–	3.0369	–	3.0338	3.0369	3.0353	3.0406	3.0395	3.0369	3.0378
J0125$-$0129	–	2.6534	–	2.6466	2.6494	–	–	2.6540	2.6534	2.6618
J0136$+$0440	–	–	2.7905	2.7820	2.7849	2.7800	2.7814	2.7846	2.7849	–
J0858+1749	–	–	2.6643	2.6450	2.6478	–	–	2.6589	–	–
J1236$+$0010	–	–	3.0284	3.0208	3.0239	–	–	3.0318	3.0284	–
J1248$+$0639	–	–	2.5273	2.5206	2.5233	–	–	2.5221	2.5273	–
J1259$+$0309	–	–	3.2365	3.2277	3.2309	–	–	3.2369	3.2365	3.2405
J1331$+$0206	2.9179	2.9147	–	2.9099	2.9129	2.9100	2.9155	2.9153	2.9147	2.9168
J1358$+$1410	2.8992	–	2.8926	2.8971	2.9001	–	–	2.8996	2.8926	–
J2228$-$0221	–	–		2.7676	2.7705	–	–	2.7763	2.7705	2.7780
J2325$+$1539	–	–	–	2.6044	2.6072	–	–	2.6051	2.6072	–

Table 3. Table 3: Main properties of the proximate absorption systems

Quasar	$\log N\left(H\text{i}\right)$	$\log N\left(H_2\right)$	[Zn/H]	[S/H]	[Si/H]	[Fe/H]	AVa𝑎aa𝑎a$a$Statistical error only, from dispersion between measurements. The systematic error on $A_V$ is about 0.07 mag (see text).	$\Delta v_{90}$b𝑏bb𝑏b$b$Measured from the low-ion metal lines (see text).
	(${\mathrm{cm}}^{-2}$)	(${\mathrm{cm}}^{-2}$)					(mag)	($\mathrm{km}{s}^{-1}$)
J0015+1842	20.71$\pm$0.12	19.73 $\pm$ 0.02	-0.41$\pm$0.13	-0.54$\pm$0.14	-0.78$\pm$0.12	-1.79$\pm$0.12	0.40$\pm$0.10	450
J0019-0137	21.02$\pm$0.11	20.29 $\pm$ 0.01	-0.90$\pm$0.12	-1.00$\pm$0.11	-0.96$\pm$0.12	-2.03$\pm$0.11	0.08$\pm$0.01	190
J0059+1124	21.77$\pm$0.11	19.07 $\pm$ 0.02	-1.40$\pm$0.15	-1.37$\pm$0.16	-1.03$\pm$0.34	-2.02$\pm$0.13	0.04$\pm$0.02	170
J0125-0129	21.90$\pm$0.15	20.05 $\pm$ 0.01	-0.96$\pm$0.15	-0.96$\pm$0.16	-0.92$\pm$0.16	-1.66$\pm$0.15	0.19$\pm$0.01	430
J0136+0440	20.73$\pm$0.11	18.65 $\pm$ 0.07	–	-0.59$\pm$0.11	-0.80$\pm$0.12	-1.24$\pm$0.11	0.04$\pm$0.01c𝑐cc𝑐c$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.	190
J0858+1749	20.40$\pm$0.11	19.72 $\pm$ 0.02	-0.52$\pm$0.12	-0.57$\pm$0.12	-0.96$\pm$0.11	-2.26$\pm$0.11	0.12$\pm$0.01c𝑐cc𝑐c$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.	190
J1236+0010	20.65$\pm$0.11	19.75 $\pm$ 0.01	-0.55$\pm$0.13	-0.48$\pm$0.11	-0.71$\pm$0.11	-1.70$\pm$0.11	0.07$\pm$0.02	400
J1248+0639	20.55$\pm$0.11	19.76 $\pm$ 0.01	-0.55$\pm$0.16	-0.48$\pm$0.19	-0.89$\pm$0.12	-1.76$\pm$0.11	-0.05$\pm$0.02	240
J1259+0309	21.42$\pm$0.11	19.10 $\pm$ 0.01	-1.57$\pm$0.12	-1.39$\pm$0.11	-1.33$\pm$0.11	-1.85$\pm$0.11	-0.05$\pm$0.01	90
J1331+0206	21.41$\pm$0.11	20.11 $\pm$ 0.01	-1.56$\pm$0.14	-1.67$\pm$0.12	-1.61$\pm$0.12	-1.94$\pm$0.12	0.04$\pm$0.02	80
J1358+1410	21.19$\pm$0.11	19.85 $\pm$ 0.02	-0.96$\pm$0.14	-0.96$\pm$0.11	-0.87$\pm$0.11	-1.76$\pm$0.11	0.06$\pm$0.03	170
J2228-0221	20.98$\pm$0.11	19.44 $\pm$ 0.02	-0.52$\pm$0.11	-0.54$\pm$0.11	-1.00$\pm$0.11	-1.94$\pm$0.11	0.24$\pm$0.03	110
J2325+1539	21.74$\pm$0.11	20.25 $\pm$ 0.01	-1.15$\pm$0.11	-1.03$\pm$0.11	-1.14$\pm$0.11	-1.70$\pm$0.11	0.16$\pm$0.03	240

Table 4. Table 4: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J0019 − - 0137. The adopted spectral resolutions are 44.3 and 27.0  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Mn ii
2.5282	19.4$\pm$0.2	15.03$\pm$0.02	15.50$\pm$0.04	14.53$\pm$0.01	12.54$\pm$0.05	12.52$\pm$0.04
2.5290	20.9$\pm$0.8	13.97$\pm$0.10	14.01$\pm$0.02	13.36$\pm$0.01	$\dots$	$\dots$
2.5299	5.7$\pm$0.2	14.89$\pm$0.06	15.29$\pm$0.09	13.48$\pm$0.03	12.60$\pm$0.05	12.06$\pm$0.24
Total	-	15.29$\pm$0.03	15.71$\pm$0.04	14.59$\pm$0.01	12.87$\pm$0.04	12.65$\pm$0.07

Table 5. Table 5: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J0059 + + 1124. The adopted spectral resolutions are 45.0 and 25.0  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Cr ii	Ni ii
3.0343	29.0$\pm$1.1	15.22$\pm$0.05	15.72$\pm$0.09	15.02$\pm$0.06	12.38$\pm$0.24	12.77$\pm$0.55	13.13$\pm$0.70
3.0345	7.8$\pm$2.2	15.26$\pm$0.23	16.09$\pm$0.45	14.77$\pm$0.14	12.86$\pm$0.11	13.57$\pm$0.10	14.10$\pm$0.10
Total	-	15.54$\pm$0.12	16.25$\pm$0.32	15.21$\pm$0.06	12.98$\pm$0.10	13.64$\pm$0.11	14.14$\pm$0.11

Table 6. Table 6: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J1236 + + 0010. The adopted spectral resolutions are 50.0 and 26.5  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Ni ii
3.0258	5.8$\pm$0.5	$\dots$	13.18$\pm$0.11	12.12$\pm$0.15	$\dots$	$\dots$
3.0264	17.3$\pm$1.8	$\dots$	13.70$\pm$0.04	12.99$\pm$0.03	$\dots$	$\dots$
3.0270	5.9$\pm$1.0	$\dots$	13.57$\pm$0.10	12.74$\pm$0.07	$\dots$	$\dots$
3.0283	35.0$\pm$1.0	14.71$\pm$0.04	14.61$\pm$0.03	13.73$\pm$0.02	12.02$\pm$0.29	$\dots$
3.0290	11.5$\pm$1.0	$\dots$	14.12$\pm$0.08	13.43$\pm$0.05	11.73$\pm$0.24	$\dots$
3.0304	25.1$\pm$1.0	$\dots$	13.79$\pm$0.02	12.90$\pm$0.04	$\dots$	$\dots$
3.0315	23.3$\pm$3.4	$\dots$	13.21$\pm$0.02	12.40$\pm$0.11	$\dots$	$\dots$
3.0323	18.3$\pm$0.3	$\dots$	13.71$\pm$0.03	12.83$\pm$0.05	$\dots$	$\dots$
3.0331	12.9$\pm$0.4	15.30$\pm$0.03	15.43$\pm$0.04	14.24$\pm$0.02	12.69$\pm$0.03	13.16$\pm$0.10
3.0335	18.0$\pm$0.9	$\dots$	13.95$\pm$0.04	13.48$\pm$0.04	$\dots$	$\dots$
Total	-	15.40$\pm$0.03	15.55$\pm$0.03	14.51$\pm$0.01	12.81$\pm$0.06	13.16$\pm$0.10

Table 7. Table 7: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J1248 + + 0639. The adopted spectral resolutions are 40.0 and 21.0  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Mn ii
2.5277	19.6$\pm$0.8	14.01$\pm$0.08	13.62$\pm$0.02	13.03$\pm$0.02	$\dots$	$\dots$
2.5284	8.5$\pm$0.6	$\dots$	13.35$\pm$0.04	12.72$\pm$0.03	$\dots$	$\dots$
2.5292	18.7$\pm$0.3	14.79$\pm$0.03	15.08$\pm$0.04	14.11$\pm$0.01	12.30$\pm$0.07	12.27$\pm$0.06
2.5295	2.7$\pm$0.4	15.12$\pm$0.26	14.45$\pm$0.28	13.52$\pm$0.10	12.25$\pm$0.14	11.53$\pm$0.33
2.5298	3.4$\pm$0.9	$\dots$	14.25$\pm$0.31	12.87$\pm$0.09	11.86$\pm$0.33	11.53$\pm$0.37
2.5300	44.3$\pm$1.2	14.02$\pm$0.18	14.37$\pm$0.01	13.69$\pm$0.01	11.96$\pm$0.53	$\dots$
Total	-	15.33$\pm$0.16	15.29$\pm$0.05	14.37$\pm$0.02	12.73$\pm$0.11	12.41$\pm$0.08

Table 8. Table 8: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J1259 + + 0309. The adopted spectral resolutions are 50.0 and 28.5  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Cr ii
3.2453	32.2$\pm$0.5	13.59$\pm$0.29	13.95$\pm$0.02	13.53$\pm$0.04	$\dots$	$\dots$
3.2461	14.5$\pm$0.2	15.15$\pm$0.02	15.59$\pm$0.02	15.01$\pm$0.02	12.46$\pm$0.05	13.46$\pm$0.03
Total	-	15.17$\pm$0.02	15.60$\pm$0.02	15.03$\pm$0.02	12.46$\pm$0.05	13.46$\pm$0.03

Table 9. Table 9: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J1331 + + 0206. The adopted spectral resolutions are 45.0 and 26.0  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Cr ii	Ni ii
2.9221	21.5$\pm$0.4	14.92$\pm$0.04	15.35$\pm$0.05	14.97$\pm$0.04	12.50$\pm$0.09	13.32$\pm$0.07	13.78$\pm$0.05
Total	-	14.92$\pm$0.04	15.35$\pm$0.05	14.97$\pm$0.04	12.50$\pm$0.09	13.32$\pm$0.07	13.78$\pm$0.05

Table 10. Table 10: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J1358 + + 1410. The adopted spectral resolutions are 43.0 and 26.0  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Cr ii	Ni ii
2.8921	15.8$\pm$0.6	14.79$\pm$0.04	15.10$\pm$0.08	14.38$\pm$0.03	12.12$\pm$0.19	12.26$\pm$0.67	13.14$\pm$0.16
2.8926	30.6$\pm$1.4	15.09$\pm$0.03	15.62$\pm$0.04	14.66$\pm$0.02	12.57$\pm$0.11	12.98$\pm$0.18	13.51$\pm$0.08
2.8931	48.5$\pm$2.8	14.80$\pm$0.05	15.28$\pm$0.07	13.83$\pm$0.06	12.42$\pm$0.16	12.66$\pm$0.44	$\dots$
2.8940	15.0$\pm$2.3	$\dots$	13.75$\pm$0.31	13.62$\pm$0.04	$\dots$	$\dots$	$\dots$
2.8943	28.4$\pm$1.1	$\dots$	13.02$\pm$0.28	13.46$\pm$0.03	$\dots$	$\dots$	$\dots$
2.8959	23.7$\pm$1.0	13.54$\pm$0.57	13.66$\pm$0.02	13.08$\pm$0.03	$\dots$	$\dots$	$\dots$
Total	-	15.40$\pm$0.02	15.87$\pm$0.03	14.93$\pm$0.02	12.88$\pm$0.08	13.20$\pm$0.18	13.66$\pm$0.07

Table 11. Table 11: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J2228 − - 0221. The adopted spectral resolutions are 45.5 and 28.0  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Cr ii	Mn ii	Mg ii	Ti ii
2.7691	18.5$\pm$0.2	15.59$\pm$0.02	15.49$\pm$0.03	14.50$\pm$0.01	13.09$\pm$0.02	12.64$\pm$0.13	12.64$\pm$0.05	15.36$\pm$0.10	12.44$\pm$0.17
2.7701	10.3$\pm$0.5	13.92$\pm$0.15	14.15$\pm$0.06	13.28$\pm$0.03	$\dots$	$\dots$	$\dots$	14.10$\pm$1.62	$\dots$
Total	-	15.60$\pm$0.02	15.51$\pm$0.02	14.52$\pm$0.01	13.09$\pm$0.02	12.64$\pm$0.13	12.64$\pm$0.05	15.38$\pm$0.12	12.44$\pm$0.17

Table 12. Table 12: Results of Voigt-profile fitting to singly ionised metal absorption lines proximate to J2325 + + 1539. The adopted spectral resolutions are 50.0 and 26.5  km ​ s − 1 km superscript s 1 {\rm km\,s^{-1}}  for the UVB and VIS arm, respectively.

$z$	$b$	$\log N$ [cm-2]
	(km s${}_{}{}^{-1})$	S ii	Si ii	Fe ii	Zn ii	Cr ii	Ni ii	Mg ii	Ti ii	Mn ii
2.6144	9.5$\pm$0.3	15.19$\pm$0.04	15.24$\pm$0.06	13.63$\pm$0.02	12.80$\pm$0.03	$\dots$	$\dots$	15.14$\pm$0.24	$\dots$	$\dots$
2.6163	19.4$\pm$0.5	15.01$\pm$0.04	15.46$\pm$0.03	15.05$\pm$0.02	12.29$\pm$0.07	13.35$\pm$0.03	13.80$\pm$0.03	15.29$\pm$0.17	11.97$\pm$0.78	12.73$\pm$0.03
2.6169	23.3$\pm$0.3	15.68$\pm$0.02	15.96$\pm$0.02	15.34$\pm$0.02	12.93$\pm$0.04	13.68$\pm$0.02	14.17$\pm$0.02	16.14$\pm$0.03	12.90$\pm$0.10	13.09$\pm$0.02
Total	-	15.87$\pm$0.02	16.14$\pm$0.02	15.53$\pm$0.01	13.23$\pm$0.02	13.84$\pm$0.01	14.32$\pm$0.01	16.23$\pm$0.04	12.95$\pm$0.12	13.25$\pm$0.01

Table 13. Table 13: Results from Gaussian fits of CO(3-2) emission lines.

Quasar	$z_{CO}$	FWHM	$F_{\mathrm{CO}\left(3-2\right)}$	$M\left(H_2\right)$
		($\mathrm{km}{s}^{-1}$)	(Jy $\mathrm{km}{s}^{-1}$)	$\left({10}^{10}{M}_{\odot }\right)$
J0015$+$1842	2.6285	1010 $\pm$ 120	1.1 $\pm$ 0.2	3.4 - 17
J0059$+$1124	3.0378	500 $\pm$ 110	0.24 $\pm$ 0.07	0.9 - 4.6
J0125$-$0129	2.6618	460 $\pm$ 80	0.5 $\pm$ 0.1	1.5 - 7.4
J1259$+$0309	3.2405	450 $\pm$ 180	0.2 $\pm$ 0.1	0.8 - 3.9
J1331$+$0206	2.9168	350 $\pm$ 100	0.22 $\pm$ 0.08	0.8 - 4.0
J2228$-$0221	2.7780	330 $\pm$ 60	0.27 $\pm$ 0.07	0.9 - 4.5

Full text

1 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-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 ${\rm km\,s^{-1}}$ up to $\sim$0.3 $c$ (Hamann et al., 2018). Damped Lyman-$\alpha$ systems (DLAs) are absorption systems with high enough H i column density ($N({H\textsc{i}})>2\times 10^{20}$ ${\rm 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 103 ${\rm 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-$\alpha$ 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 104 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-$\alpha$ emission spanning a range of strengths, up to the point where the H i Ly-$\alpha$ 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-$\alpha$ from the AGN that is not fully covered by the absorber. Because Ly-$\alpha$ 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 H2 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 H2 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 detected 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 post-processing. 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 H2 systems in Sect. 5. We conclude in Sect. 6.

2 Sample, observations and data reduction

2.1 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 $\delta<20^{\circ}$, 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 P105111Because of delays related to the observatory shutdown during the Covid-19 pandemic, the observations were completed in P109.) 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)222We note that SDSS J123602.11$+$001024.5 had already been observed by Balashev et al. (2019). This quasar was then re-observed under better conditions and from various position angles thanks to the availability of the X-Shooter atmospheric dispersion corrector. We use here the new data.. 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 $\mu$m, splitting the incoming light beam into three spectrographs (’arms’): UVB ranging from 0.3 to 0.6 $\mu$m; VIS ranging from 0.6 to 1.0 $\mu$m; and NIR ranging from 1.0 to 2.5 $\mu$m. All observing blocks (OBs) consisted of 2980 s, 3010 s and 5$\times$600 s exposures for the UVB, VIS and NIR arms, respectively obtained in ’stare’ mode. We used 1.0 and 0.9$\arcsec$-wide slits in the UVB and VIS, with slow read-out mode and 1$\times$2 pixel binning. The NIR observations were performed with 1.2$\arcsec$-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$\leavevmode\nobreak\ \mu$m for all remaining quasars but SDSS J013644.02$+$044039.1 and SDSS J085859.67+174925.1 (in P94) for which the K-band 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 version 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).

2.2 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 O2 and H2O, while the NIR is strongly affected by H2O and CO2 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 quasar -typically 0.68-0.69, 0.71-0.73, 0.76-0.77, 0.81-0.82 and 0.91-0.91 $\mu$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_{\rm 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.

2.3 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.

3 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, H2 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, H2 and metals. The results from the fits are presented in the Appendix, with the main derived properties summarised in Table. 3.

3.1 Total HI and H2 column densities

The total H i column densities were mostly constrained from the saturated core of the Ly-$\alpha$ line, which is little sensitive to the exact placement of the continuum on top of the Ly-$\alpha$ 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 ${\rm 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.02 dex), 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 X-shooter (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 H2 absorption lines from the Lyman and Werner bands. These are typically detected from rotational level J=0 to J$\sim$6-7 in our X-shooter 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 H2 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 H2 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-rotational levels. The fits to the H2 absorption bands is presented in Appendix A.

In Fig. 1, we compare the total H2 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 H2, 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 low-resolution, low S/N SDSS spectra (see Paper I). This agreement is similar to what was found by Balashev et al. (2019) for intervening ${\rm H}_{2}$ systems.

3.2 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$\lambda 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\textsc{ii}}\lambda$2026 is systematically blended with ${Mg\textsc{i}}\lambda 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$\lambda$1827 and/or Mg i$\lambda$2852 to further constrain the (small) contribution of Mg i to the 2026 profile. We also took into account partial blending of ${Zn\textsc{ii}}$ with ${Cr\textsc{ii}}$ at 2062 Å. Finally, the analysis of neutral carbon (incl. fine-structure lines), neutral chlorine 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 ${\rm km\,s^{-1}}$) components than typically observed in high-resolution spectra, where $b<10\leavevmode\nobreak\ {\rm 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.

3.3 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 $\AA$. 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 quasar444This implicitly assumes that no other intervening absorber is causing the reddening. Such an assumption is reasonable since the proximate systems likely contain most of the dust present along the line of sight as indicated by the presence of ${\rm H}_{2}$ and the low probability to intercept by chance another rare dust-rich absorber. Indeed, we do not find signs for additional metal-rich intervening absorbers in any of the quasar spectra.. 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.

4 Absorption line kinematics and systemic redshifts

4.1 Absorption line kinematics

We quantify the kinematics of the neutral absorbing gas using the velocity width as defined by Prochaska & Wolfe (1997): $\Delta v_{90}=c[\lambda_{95}-\lambda_{5}]/\lambda_{50}$, where $\lambda_{5}$, $\lambda_{50}$ and $\lambda_{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 $\Delta 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 $\Delta 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 $\Delta_{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.

4.2 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_{\rm sys}$), while low-ionisation lines provide a better mean to measure $z_{\rm 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$\beta$ 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 $\sigma$) 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 ($\sim 8$, $-45$ and $-57$ ${\rm km\,s^{-1}}$, respectively) compared to the systemic value, which they determined from Ca ii absorption. 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 luminosity-independent systematic shift of $-$229 ${\rm 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_{\rm syst}^{i}$), 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 ${\rm 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.

5 Results and discussions

5.1 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 ${\rm H}_{2}$ systems from Balashev et al. (2019) and confirmed from follow-up observations. This compilation contains H2-bearing systems identified with a number of different selections: in regular DLA systems (Noterdaeme et al., 2008; Balashev et al., 2014, 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 distribution in parameter space is a sum of several selections. We removed from the intervening samples some systems that have $z_{\rm abs}\approx z_{\rm 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 ${\rm 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({\rm X})/N({\rm H_{tot}}))-\log({\rm X/H})_{\odot}$, 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 self-shielding 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\textsc{ii}}$) = N(Zn)). We do however include ${\rm H}_{2}$ in the total hydrogen column density, i.e., $N({\rm H})=N({H\textsc{i}})+2N({\rm 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.

5.1.1 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_{\rm\odot}\approx[$Zn/H$]\approx[$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$]\approx 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 $\rm\log Z/Z_{\rm\odot}\equiv[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 H2 systems to span a range from $\log Z/Z_{\rm\odot}\approx-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 ${\rm H}_{2}$-bearing DLA systems (e.g. Petitjean et al., 2006; Balashev et al., 2019).

5.1.2 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) – 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/10th Solar. Interestingly, [Zn/Fe] in proximate ${\rm 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 ${\rm H}_{2}$ absorbers, even if the latter tend to have high depletion among the overall (mostly non-${\rm 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 H2 and intervening systems (regardless of their H2 content) is less clear in the case of [Zn/Si]. Since silicon is an $\alpha$-element predominantly produced in massive stars, it is possible that a stronger intrinsic depletion (with respect to intervening systems) is compensated 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 ${\rm 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_{\rm V}/N({\rm H})$. While the 0.07 mag uncertainty on AV measurements remains quite large compared to their values, there is a clear tendency for higher $A_{\rm V}/N({\rm H})$ at higher [Zn/Fe], [Zn/Si] and metallicities, strengthening a dust-depletion origin for the observed underabundance of iron and silicon.

The trend between depletion and value of A${}_{\rm V}/N(\rm H)$ turns out to be a strong correlation with Pearson coefficient $r=0.9$ and probability of chance coincidence 0.07%, see Fig. 4. The low dispersion around the approximate linear relation ($\log A_{V}/N($H$)\approx 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.

5.1.3 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 $\alpha$ elements and their intrinsic abundances are expected to follow each other. According to Ellison et al. (2010), based on modelling by Rix et al. (2007), the observed negative [S/Si] values at low H i column densities in their sample of PDLAs 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\textsc{i}})=20.75$ in our strong H2-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 PH2 have close to Solar ratios at high column densities $\log N({H\textsc{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 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 H2 at low H i column densities (e.g., Bialy et al., 2017), thereby increasing the dust shielding, H2 formation rate and gas cooling through atomic lines. In the low $N({H\textsc{i}})$ regime, we would qualify the observed [S/Si] in (non-${\rm 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 ${\rm 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. Fitting a linear relationship to the proximate H2 data, that takes into account uncertainties on both axes, we obtain ${\rm[S/Si]\approx 0.39(\pm 0.08)\times[Zn/Fe]-0.21(\pm 0.09)}$. A correlation between the two abundances ratio has also been observed in intervening DLAs (Pearson correlation $r=0.7$), for which we show the best fit linear relation obtained by De Cia et al. (2016) together with the associated intrinsic dispersion666The differences in adopted Solar values between this work and De Cia et al. (2016) only affect the [Zn/Fe] ratio by 0.005 dex and hence can safely be ignored here.. However, proximate systems differ from intervening DLAs in the sense that they have [S/Si]$\approx 0$ for a wider range of [Zn/Fe] values. This is observed for PDLAs regardless of whether they contain H2 or not,

when the intercept for intervening DLAs is almost zero.

In summary, the PH2 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\textsc{i}}$ high-depletion systems are missing due to dust bias. At low [Zn/Fe] (and high $N({H\textsc{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 ii– closer 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 low-ionisation metals could actually arise from ionised gas. Larger samples of proximate DLAs (with and without ${\rm H}_{2}$) are necessary to investigate this further.

5.1.4 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 self-shielded against H i-ionising photons. While argon is predominantly ionised by cosmic rays in local ISM, high-energy photons significantly above 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 $\alpha$-elements in intervening DLAs, with a small dependency on H i column density. They concluded that, since argon is extremely volatile (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 H2 absorbers using the two available absorption lines (at 1048 and 1066 Å). Because these fall into the Lyman-$\alpha$ forest and can easily be blended with H2/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 H2 absorption profile constrained by independent fit (see Section 3.1) and avoided blends with Ly-$\alpha$ 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. 31.

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\textsc{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-$\alpha$ interlopers so that the Ar i columns should in principle be considered as upper-limits –in a conservative approach777This should be less an issue at high spectral resolution as in the work by Zafar et al. (2014b) based on UVES data.– but more importantly, the metal column density is generally dominated by the ${\rm H}_{2}$-bearing component, where we observe lower deficit of argon (-0.2 dex) with respect to components without detectable ${\rm H}_{2}$ (-0.7 dex). In Fig. 31, 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 H2 components.

The relatively low deficit of argon in ${\rm H}_{2}$ components most likely is a result from an efficient conversion of Ar+ into neutral argon through chemical reactions with ${\rm H}_{2}$ (${\rm Ar^{+}+{\rm H}_{2}\to ArH^{+}+H}$ followed by $\rm ArH^{+}+{\rm H}_{2}\to 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 ${\rm 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_{\rm UV}$ the photo-ionisation rate increases, while the ${\rm H}_{2}$ abundance decreases, impeding the molecular-ion recombination route. Hence the ionisation fraction of argon should rise at least proportionally to $I_{\rm UV}$. Therefore, the systematically lower abundance in non-H2-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 ${\rm H}_{2}$ and metal components (i.e. at high spectral resolution) –out of the scope of this paper– we can estimate that the $-$0.2 dex deficit of argon means that about half of the metal column seen at the velocity of ${\rm H}_{2}$ components could come from their (${\rm H}_{2}$-free) H i-envelope where argon would be mostly in ionised form.

5.2 ${\rm H}_{2}$ in proximate versus intervening systems

We compare proximate and intervening systems according to their H2 and total hydrogen column densities in Fig. 8. The absence of proximate systems at low ${\rm H}_{2}$ columns is naturally due to our selection based on the presence of strong H2 Lyman-Werner bands (Paper I), and in fact, all our systems have $\log N({\rm H}_{2})>18.5$. In turn, proximate H2 were selected independently of $N({H\textsc{i}})$, but as can be seen in this figure, they are limited to $\log N({\rm H})>20.5$, when intervening systems can have lower $N(\rm H)$ even at similar ${\rm H}_{2}$ column densities (see distributions on top of this figure). However, this is driven by a few very metal-rich 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({\rm H}_{2})$. This could be naturally explained by the $1/Z$ dependence of the critical column density for ${\rm H}_{2}$ to form (e.g., Krumholz et al., 2008), but also to a bias against high-metallicity systems at high-$N({H\textsc{i}})$ due to stronger dust extinction.

This is better seen in Fig. 9, where the location H2 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 H2 cloud illuminated on one side888Following Eq. 2 of Paper I, based on theoretical calculations by Bialy et al. (2017) and where we used $Z/Z_{\rm\odot}$ as a proxy for $\tilde{\sigma_{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 ${\rm 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\textsc{i}})$-only is almost identical.. Three quarters (9/12) of the proximate H2 systems are located between the $I_{\rm UV}/n=10^{-2}$ and $10^{-1}$ cm3 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 ${\rm H}_{2}$ clouds. Notwithstanding, any atomic gas –with some amount of dust– located upstream along the line of sight will still participate in shielding H2 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 H2 systems tend to have a wider spread of values. We also note that the lines of constant $I_{\rm 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 H2 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 metal-rich proximate H2 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 ${\rm H}_{2}$ lines.

We also show in Fig. 9 the location of proximate and intervening DLAs without strong ${\rm H}_{2}$ as open symbols. In the case of intervening non-${\rm H}_{2}$ systems, we used here the sample from Noterdaeme et al. (2008) where the presence of ${\rm H}_{2}$ has been systematically assessed. The fact that no ${\rm 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_{\rm UV}/n$. For example, for a UV field about 1/10th 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 ${\rm H}_{2}$ is found in outflowing gas. Indeed, for $n=100\rm\,cm^{-3}$ the characteristic timescale for H2 formation will be $\sim 10^{7}$ yrs, which is then comparable to the dynamical timescales.

5.3 Kinematics

The velocity extent of low-ionisation metal lines, $\Delta 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 (e.g., Møller & Christensen, 2020). While the assumption that $\Delta 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 $\Delta v_{90}$ in proximate ${\rm H}_{2}$-systems follows the trend and dispersion seen in intervening DLAs, although limited to high-values, as expected from the presence of ${\rm H}_{2}$ (Noterdaeme et al., 2008). Overall, the velocity spread of the neutral gas in proximate-${\rm 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 low-metallicity systems tend to be redshifted. The Pearson correlation coefficient between $v_{50}=c(z_{\rm QSO}-z_{50})/z_{\rm 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 $\sim 2000$ ${\rm 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 ${\rm 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 $\sigma$ set to the intrinsic velocity dispersion of $z_{\rm QSO}$ provided by Shen et al. (2016) (i.e., 46, 56, 205 and 233 ${\rm km\,s^{-1}}$for $z_{\rm QSO}$ based on [O ii], [O iii], Mg ii and C iii], respectively). From this exercise, we find an approximate normal distribution of $r$ values, giving $r=-0.47\pm 0.09$ ($r=-0.44\pm 0.11$ ignoring 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 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\sim 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]$\mu$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\sim 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 ${\rm km\,s^{-1}}$ higher than those obtained from the (corrected) ionised emission lines. If $z_{\rm 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\approx-0.6$, $p\approx 0.2$ (now computed on 6 quasars). In that case, high-metallicity gas (above $\sim$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.

5.4 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 ${\rm km\,s^{-1}}$, suggesting that proximity in velocity does not necessarily mean physical proximity. Indeed, 5000 ${\rm km\,s^{-1}}$ in the Hubble flow corresponds already to about 16 Mpc proper distance at $z\sim 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 ${\rm H}_{2}$-absorbers (see Fig. 12). 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 ${\rm H}_{2}$-selected proximate systems.

6 Conclusion

The chemical enrichment and velocity dispersion of the neutral gas in strong proximate H2 systems appears not to be significantly different from what is seen in intervening systems. While the metallicities in proximate H2 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 ${\rm 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 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-${\rm H}_{2}$-bearing components. In turn, argon remains predominantly in neutral form in the the presence of ${\rm 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 ${\rm H}_{2}$-absorption 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 H2 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 H2-bearing DLAs means that proximate H2 should have similar UV-to-density ratios. This means that proximate ${\rm 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_{\rm abs}\approx z_{\rm 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 ${\rm 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 ${\rm 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_{\rm abs}>z_{\rm 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 high-metallicity 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-${\rm 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 PDLAs based on saturated cores biases the samples against systems with leaking Lyman-$\alpha$ emission, that is actually seen in about half of our proximate H2 systems. If, as expected, the amount of gas and its distribution around the quasar affects the ability of Lyman-$\alpha$ 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.

Acknowledgements.

We thank the anonymous referee for useful comments and suggestions. The research leading to these results has received support from the French Agence Nationale de la Recherche under ANR grant 17-CE31-0011-01/project “HIH2”. RC gratefully acknowledges support from the French-Chilean Laboratory for Astronomy (IRL 3386). SL acknowledges support by FONDECYT grant 1191232.

Appendix A Fit to H i Ly-$\alpha$ and H2 Lyman-Werner absorption bands

Appendix B Low-ionisation metal absorption lines

Appendix C Measurement of dust extinction

Appendix D $\Delta v_{90}$ measurements

Appendix E Optical/NIR emission lines

Appendix F CO(3-2) emission lines