Chemical enrichment and host galaxies of extremely strong intervening DLAs towards quasars (cid:63) Do they probe the same galactic environments as DLAs associated with γ -ray burst afterglows?

We present the results from VLT / X-shooter spectroscopic observations of 11 extremely strong intervening damped Lyman-α absorbers (ESDLAs) that were initially selected as high N (H i ) (i.e. ≥ 5 × 10 21 cm − 2 ) candidates from the Sloan Digital Sky Survey (SDSS). We conﬁrm the high H i column densities, which we measure to be in the range log N (H i ) = 21 . 6 − 22 . 4. Molecular hydrogen is detected with high column densities ( N (H 2 ) ≥ 10 18 cm − 2 ) in 5 out of 11 systems, 3 of which are reported here for the ﬁrst time, and we obtain conservative upper limits on N (H 2 ) for the remaining 6 systems. We also measure the column density of various metal species (Zn ii , Fe ii , Si ii , Cr ii , and C i ), quantify the absorption-line kinematics ( ∆ v 90 ), and estimate the extinction of the background quasar light ( A V ) by dust in the absorbing gas. We compare the chemical properties of this sample of ESDLAs, supplemented with literature measurements, to that of DLAs located at the redshift of long-duration γ -ray bursts (GRB-DLAs). We conﬁrm that the two populations are almost indistinguishable in terms of chemical enrichment and gas kinematics. In addition, we ﬁnd no marked di ﬀ erences in the incidence of H 2 . All this suggests that ESDLAs and GRB-DLAs probe similar galactic environments. We search for the galaxy counterparts of ESDLAs and ﬁnd associated emission lines in 3 out of 11 systems, 2 of which are reported here for the ﬁrst time (at z abs = 2.304 and 2.323 towards the quasars SDSSJ002503.03 + 114547.80 and SDSSJ114347.21 + 142021.60, respectively). The measured separations between the quasar sightlines and the emission associated with the ESDLA galaxy (for a total of ﬁve sightlines) are all very small ( ρ < 3kpc). Because our observations are complete up to ρ ∼ 7kpc, we argue that the emission counterparts of the remaining systems are more likely below the detection limit than outside the search area. While the small impact parameters are similar to what is observed for GRB-DLAs, the associated star formation rates are on average lower than for GRB host galaxies. This is explained by long-duration GRBs being associated with the death of massive stars and therefore pinpointing regions of active star formation in the GRB host galaxies. Our observations support the suggestion from the literature that ESDLAs could act as blind analogues of GRB-DLAs, probing neutral gas with high column density in the heart of high-redshift galaxies, without any prior on the instantaneous star formation rate.


Introduction
Damped Lyα absorbers (hereafter DLAs) that are observed in the spectra of luminous background sources, such as quasars (QSOs) or γ-ray burst (GRB) afterglows, are defined to be systems with column densities log N(H i) ≥ 20.3 (Wolfe et al. 2005).DLAs are found to probe most of the neutral gas content of the Universe (see Noterdaeme et al. 2009;Prochaska et al. 2009).
Simulations (e.g.Pontzen et al. 2008;Rahmati & Schaye 2014;Yajima et al. 2012;Altay et al. 2013) have shown that at high redshift, intervening DLAs mostly probe gas in galactic haloes or in the circum-galactic medium (CGM).However, these simulations also indicated that the systems with the highest Based on observations performed with Very Large Telescope of the European Southern Observatory under Prog.ID 095.A-0224(A) using the X-shooter spectrograph.
neutral gas column densities (N(H i) ∼ 10 22 cm −2 ) are expected to probe dense gas that is located closer to the centre of the associated galaxy, at impact parameters smaller than a few kiloparsec (kpc).Interestingly, systems at the high column density end of the N(H i) distribution provide strong constraints on models of stellar feedback and formation of H 2 in simulations (see e.g.Altay et al. 2013;Bird et al. 2014).Nevertheless, the incidence of these intervening systems along the quasar lines of sight is rare because the projected area of this dense gas is small.
One way to probe the gas in the central regions of galaxies is then to study damped Lyα absorbers that are located at the redshift of long-duration GRBs and are observed in their afterglow spectra.Hereafter, we use the term "GRB-DLA" to refer to systems at the GRB redshift as opposed to intervening GRB-DLAs with z DLA z GRB , while "QSO-DLA" refers to intervening DLAs along quasar sightlines, as opposed to proximate QSO-DLAs with z DLA ≈ z QSO .This is in agreement with the standard terminology in the literature.
Long-duration GRBs are thought to be associated with the death of massive stars (e.g.Woosley & Bloom 2006), and are indeed found to be located in the actively star-forming regions of their host galaxies (see Fruchter et al. 2006;Lyman et al. 2017).
As expected, the N(H i) distribution of GRB-DLAs is found to be skewed towards high N(H i) values (see Fynbo et al. 2009;Selsing et al. 2019).
However, the main difficulty with studying GRB-DLAs is that high-resolution spectroscopic observations of GRB afterglows are difficult to obtain due to the rapid decrease of their luminosity.The situation becomes further complicated in the presence of dust in the associated absorber.Dark GRBs have high dust attenuation (A V > 0.5 mag) and are therefore difficult to detect because the GRB afterglow becomes too faint for optical follow-up (see e.g.Ledoux et al. 2009;Bolmer et al. 2019).
The very large number of quasar spectra available in the Sloan Digital Sky Survey (SDSS; York et al. 2000) now enables building samples of intervening absorbers with very high H i column densities.Noterdaeme et al. (2014) have found 104 systems with log N(H i) ≥ 21.7 (they coined the name extremely strong DLAs, ESDLAs, for these systems) in the spectra of 140 000 quasars from the Baryon Oscillation Spectroscopic Survey (BOSS) of the SDSS-III Data Release 11 (Pâris et al. 2017).From Lyα emission detected in stacked spectra of ESDLAs, Noterdaeme et al. (2014) suggested that these systems indeed arise at small impact parameters (typically ρ < 2.5 kpc, where ρ is defined as the projected separation between the centroid of the nebular emission and the quasar sightline), a result that was also substantiated by two direct detections at ρ ∼ 1 kpc reported by Noterdaeme et al. (2012) and Ranjan et al. (2018).
Furthermore, ESDLAs are found to have a higher density and H 2 incidence rate than the "normal" DLA population (Noterdaeme et al. 2015a;Balashev & Noterdaeme 2018).In general, all these findings are in line with the expectations based on the Kennicutt-Schmidt law, which implies that probing neutral gas at the highest column densities can help reveal molecular gas and star-forming regions originating inside the optical disc of the associated galaxy.
Several works have accumulated evidence for a similarity between GRB-DLAs and strong QSO-DLAs.Guimarães et al. (2012) showed that the metallicity and depletion of high N(H i) QSO-DLAs (log N(H i) > 21.5) are similar to GRB-DLAs.Noterdaeme et al. (2012) showed that the properties of one detected ESDLA galaxy (SDSS J 1135−0010), at very small impact parameter from the quasar line of sight, is similar to some GRB host galaxies.Furthermore, Noterdaeme et al. (2014Noterdaeme et al. ( , 2015b) ) discussed the similarity in absorption characteristics of ESDLAs and GRB-DLAs.
More recently, Bolmer et al. (2019) showed that there is no apparent lack of H 2 in their GRB-DLA sample as compared to QSO-DLAs, supporting the previous findings of the similarity in chemical enrichment as well as providing indirect evidence of the similarity of the surrounding UV field between the two absorber sub-sets.The presence of H 2 indicates that the sightlines pass close to the centre of the associated galaxies, where the gas pressure favours the H i-H 2 transition (see Blitz & Rosolowsky 2006;Balashev et al. 2017).Furthermore, impact parameter measurements of GRB-DLAs by Lyman et al. (2017) and Arabsalmani et al. (2015) give direct evidence for the proximity of the absorbing gas to its associated galaxy.
It is, however, important to note that ESDLA observations are scarce, which has led to very small and heterogeneous samples of ESDLAs in the past.In order to further study the chemical enrichment, the molecular content, and the association with galaxies of ESDLAs, we have therefore started a follow-up campaign of ESDLAs using the Very Large Telescope (VLT) of the European Southern Observatory.As part of this campaign, 11 ESDLA systems have been observed with the X-shooter instrument.This more than doubles the previous sample size.
A detailed analysis of one of these systems has already been presented by Ranjan et al. (2018).We present here the analysis of the remaining 10 ESDLAs and compare the properties of ESDLAs in our sample (supplemented with literature measurements) to those of GRB-DLAs.The observations and data reduction are presented together with the compilation of literature data in Sect. 2. The analysis of the absorption lines and dust reddening is described in Sect. 4 and that of associated emission lines in Sect. 5. We then discuss our results in Sect.6 and summarise our findings in Sect.7. Column densities are always stated in units of cm −2 .A standard ΛCDM flat cosmology is used throughout, with H 0 = 67.8km s −1 Mpc −1 , Ω Λ = 0.692, and Ω m = 0.308 (Planck Collaboration XIII 2016).

Observations and data reduction
Observations of the 11 quasars were carried out in service mode under good seeing conditions (typically 0.7-0.8 ) between April 2015 and July 2016 under ESO program ID 095.A-0224(A) with the multiwavelength medium-resolution spectrograph X-shooter (Vernet et al. 2011) mounted at the Cassegrain focus of the Very Large Telescope (VLT-UT2) at Paranal, Chile.The observations were carried out using a two-step nodding mode with an offset of 4 arcsec between the two integrations.Because the atmospheric dispersion corrector failed during this period, we aligned the slit in all observations with the parallactic angle at the start of each exposure; the slit angle then remained fixed on the sky during the integration.The parallactic angle changed little between different observations of a given target (±15 • ), except for the quasar SDSS J 223250.98+124225.29.In the following, we use a short notation for the quasar names, for example SDSS J 223250.98+124225.29 is referred to as J2232+1242.The log of observations is given in Table B.1.We reduced the data using the X-shooter pipeline (Modigliani et al. 2010) and combined individual exposures by weighting each pixel by the inverse of its variance to obtain the combined 2D and 1D spectra we used.For absorption line and dust analysis, we used only the combined 1D spectra.For the emission-line analysis, we used not only the combined 2D spectra, but also the individual 2D spectra.

Literature sample
Given the small size of our sample (11 systems), we included measurements for other ESDLAs selected from the same parent SDSS sample (Noterdaeme et al. 2014).The system towards J 1513+0352 has been singled out by Ranjan et al. (2018), but it belongs to the same selection and observational program that we present in this work.Four additional systems (SDSS J0154+1395, SDSS J0816+1446, SDSS J1456+1609, and SDSS J2140−0321) have also been selected based on high N(H i) content from SDSS, but they were observed with UVES instead of X-shooter (Guimarães et al. 2012;Noterdaeme et al. 2015a).One of them (SDSS J2140−0321) is in common with the present X-shooter sample.All these constitute our homogeneous (indicating homogeneity in selection based on N(H i) measurement from SDSS) ESDLA sample (14 systems).We furthermore built A125, page 2 of 22 a "total" sample of ESDLAs that meet the N(H i) criterion, but have been selected based on different criteria.In this total sample, we include the "homogeneous sample" together with the following targets: SDSS J0843+0221, which was targeted for the clear presence of H 2 lines seen in the SDSS spectrum (Balashev et al. 2014), SDSS J1135−0010, which was targeted for its Lyα emission (Noterdaeme et al. 2012), as well as three additional ESDLAs that are not from the SDSS sample (HE0027−1836, Q0458−0203, and Q1157+0128), but have coverage of the H 2 lines (Noterdaeme et al. 2015b).
For the comparison between absorption properties, we used only the homogeneous sample of ESDLAs.However, for the comparison of impact parameters, we use the total sample of ESDLAs because this comparison is mainly qualitative in nature, given the small sample size.In total, five emission counterparts have been detected towards ESDLAs: J1135+0352 (Noterdaeme et al. 2012;Kulkarni et al. 2012), Q0458-0203 (Møller et al. 2004;Krogager et al. 2012), J0025+1145, J1143+1420 (both from this work), and J1513+0352 (Ranjan et al. 2018).
For the purpose of comparing the absorption properties of ESDLAs with GRB-DLAs, we used the GRB-DLA sample listed by Bolmer et al. (2019) because their sample size (22 systems) is similar to ours and their targets have also been observed using VLT/X-shooter.For the analysis of kinematics (∆v 90 ), we combined measurements from Arabsalmani et al. (2015Arabsalmani et al. ( , 2018) ) to build a larger sample of GRB-DLAs (13 systems, excluding low-resolution data).Impact parameter measurements for GRB-DLAs were compiled from the following works: Castro et al. (2003), Perley et al. (2013), Thöne et al. (2013), Arabsalmani et al. (2015), and Lyman et al. (2017).

Absorption-line analysis
We analysed the absorption lines using standard Voigt-profile fitting, taking the atomic and molecular wavelengths, oscillator strengths, and damping constants from the source website of VPFIT 1 (Carswell & Webb 2014), where most records are taken from Morton (2003).We present here the general method and measurements for our sample of ESDLAs, but most of the figures showing multicomponent Voigt profile fits to the data are presented in the appendix to keep the main text concise.

Neutral atomic hydrogen
The column density of atomic hydrogen was measured by fitting the damped Lyα line together with other lines from the H i Lyman series when covered by the spectra.Because the strong damping wings affect the apparent quasar continuum over a wide range of wavelengths, we first estimated an approximate continuum by eye taking reference from the quasar template by Selsing et al. (2016) and then used Chebyshev polynomials to model the remaining fluctuations, while fitting the H i lines simultaneously.
The best-fit Voigt profiles to DLAs are shown in Fig. 1.
In Fig. 2 we compare the column densities obtained from VLT follow-up spectroscopy with those originally obtained from the low-resolution, low signal-to-noise ratio (S/N) data by Noterdaeme et al. (2014).There appear to be no systematic differences between the two studies, and most values agree to within 0.1 dex.This shows that the SDSS-based column densities are quite robust at the high end of the H i column density distribution.Nevertheless, our X-shooter analysis shows that the 1 https://www.ast.cam.ac.uk/~rfc/vpfitFig. 1.Fits to the damped Lyman-α lines for the ten systems analysed in this work.The normalised X-shooter spectra are shown in black with the best-fit Voigt profile over-plotted in red.Shaded regions depict statistical uncertainties.The relative velocity scale is defined with respect to the absorption redshift quoted in each panel.column densities of a few of the systems (towards J0017+1307, J1143+1420, and J2322+0033) are slightly below the ESDLA threshold that was initially defined by Noterdaeme et al. (2014), N(H i) ≥ 5 × 10 21 cm −2 .We continue to call them ESDLAs in the following for brevity.

Molecular hydrogen
We searched for H 2 absorption lines in the Lyman and Werner bands redshifted to the wavelength range covered by the UVB arm of X-shooter for all systems in our sample.In addition to the systems towards J 1513+0352 and J 2140−0321, for which the detection of H 2 has previously been reported by Ranjan et al. (2018) and Noterdaeme et al. (2015a), respectively, we also clearly detect H 2 absorption lines towards J 1143+1420 and J 2232+1242.The spectrum of J 0025+1145 presents strong absorption at the expected position of H 2 , but only one band is available due to a Lyman-break caused by a system at higher redshift.Because J 0025+1145 also presents strong reddening and C i absorption lines, the detection of a high H 2 column density is expected and treated as a firm detection throughout the paper.Although there is an indication of damping wings and the data are best fitted with log N(H 2 ) ∼ 20, this column density remains highly uncertain.
For all these systems, we measured the H 2 column densities in the lowest rotational levels (typically J = 0, 1, 2) by simultaneously fitting the corresponding absorption lines together with the continuum modelled by Chebyshev polynomials.The H 2 redshifts were allowed to vary independently from those of H i because the observed H i redshift is just indicative of the middle point in the saturated Lyα profile.We find that H 2 always corresponds to a component seen in the profile of metals.Figure 3 shows one such example.While for some of the H 2 -bearing systems, we do see signs of absorption lines from high rotational levels as well, we caution that the column densities could be more uncertain than the uncertainties from the best fit suggest because the spectral resolution is insufficient.Notwithstanding, in the case of SDSS J 2140−0321, the values derived using X-shooter agree surprisingly well (within 0.2 dex) with those derived using available UVES data by Noterdaeme et al. (2015a) for all rotational levels up to J = 4 (see Table 1), although the statistical errors obtained here from the fitting of the X-shooter data appear to be underestimated.This is due to the need to assume a single b-parameter for all J levels when X-shooter data were fitted, while this was left free when the UVES data were fitted.
For the six remaining systems, the presence of H 2 is difficult to ascertain due to blending with the Lyman-α forest at the   spectral resolution of X-shooter and/or a small number of covered transitions.Instead, we measured conservative upper limits by identifying the highest column density model that was still consistent with the observed spectra.To do this, we created synthetic H 2 absorption models with a fixed excitation temperature, T ex = 100 K, as is typically seen in H 2 -bearing DLAs, a correspondingly low Doppler parameter b = 1 km s −1 , and the increasing total N(H 2 ) value convolved with the appropriate instrumental resolution to mimic the X-shooter spectrum.
We also tested redshifts within a 150 km s −1 window centred on the strongest metal component.For each model, we compared the distribution of positive residuals only (negative residuals are mostly due to absorption from the Lyα forest) with the expected distribution given the error spectrum and assuming intrinsically Gaussian-distributed uncertainties.The maximum allowed column density is obtained when the fraction of positive residuals above 1σ and 2σ both exceeds the 32% and 4.5% level, respectively.Our conservative upper limit is then given by the model with the highest column density that is still consistent with the data.The upper limits we obtained range from log N(H 2 ) ∼ 16 to 18.3.

Metal column densities
Absorption lines from metals in various ionisation stages are seen to be associated with the ESDLAs.We are mostly interested in the chemical enrichment in the neutral and the molecular gas phase, therefore we focus on the low-ionisation species, Fe ii, Si ii, Cr ii, and Zn ii, whose first ionisation energy is below 13.6 eV, but whose second ionisation energy is above this value.These species were fitted simultaneously assuming a common velocity structure, that is, the velocities and Doppler parameters for various elements were tied together.We also included neutral magnesium because Zn iiλ2026 is easily blended with Mg iλ2026.While being fitted simultaneously with these species, we allowed Mg i to have a different velocity structure because it does not correspond to the main ionisation stage of magnesium, which is mostly Mg ii.
The number and location of the velocity components were first identified visually from the numerous Fe ii and Si ii lines to serve as an initial guess.We then fitted the metal species, adjusting the number of components if necessary.During the fitting process, we noted that the main Si ii lines (including even the relatively weak Si iiλ1808) were intrinsically saturated because of the very high column densities involved.However, weaker lines such as Si ii λ2335 remain undetectable at the achieved S/N with X-shooter.The resulting Si ii column densities are therefore highly uncertain.

Neutral carbon
We searched for the presence of neutral carbon lines in our Xshooter spectra of ESDLAs because this species is known to be a good tracer of molecular hydrogen (e.g.Srianand et al. 2005;Noterdaeme et al. 2018).C i lines are detected in four systems of our sample, towards J 0025+1145, J 1258+1212, J 1513+0352, and J 2140−0321.All of them have clear H 2 lines, except J 1258+1212, for which we could only obtain a conservative limit of log N(H 2 ) < 18.3.We note that while the spectra J1143+1420 and J2232+1242 feature clear H 2 lines, we do not detect C i lines.However, detectable C i lines at the achieved spectral resolution and S/N are not necessarily expected.C i is also dependent on the metallicity (Noterdaeme et al. 2018;Zou et al. 2018;Heintz et al. 2019), and strong H 2 systems with low C i column densities have been identified in the literature (e.g.Balashev et al. 2017).
When we detected them, we simultaneously fitted the C i lines in various fine-structure levels of the electronic ground state, denoted as C i (J = 0), C i* (J = 1) and C i** (J = 2).
The velocities and Doppler parameters for each component are assumed to be the same over the different fine-structure levels, consistent with what is usually seen in C i absorbers, even at high spectral resolution.However, the fitting procedure did not allow us to obtain accurate Doppler parameters because the C i lines are intrinsically much narrower than the instrumental line spread function.The resulting column densities are thus to be considered with care and may well represent lower limits to the actual values.

Dust extinction
The extinction of the quasar light by dust in the absorbing system was estimated by fitting a quasar composite spectrum to the data with various amounts of reddening, A V , applied in the restframe of the absorber.For this purpose, we used the quasar template spectrum by Selsing et al. (2016).The appropriate reddening law is unknown, therefore we fitted various reddening laws for each target in our sample.The extinction laws for the average Small Magellanic Cloud (SMC), the average Large Magellanic Cloud (LMC), and the average LMC super-shell (LMC2) were described using the parameters obtained for different environments by Gordon et al. (2003).For each target, we shifted the template spectrum to the given quasar redshift and re-sampled the template onto the observed wavelength grid.The template was smoothed with a 20 ixel top-hat filter in order to avoid biasing the fit by noise peaks in the empirical template spectrum.We then applied the reddening in the rest-frame of the absorber with a variable amount of rest-frame V-band extinction, A V .This is the only free parameter in the fit because the absolute flux-scaling is obtained by normalising the reddened template to match the observed spectrum.The best-fit values of A V were obtained by a standard χ 2 minimisation, and the resulting best-fit values are provided in Table 3.
The main uncertainty for the extinction measurement using this method comes from the fact that the intrinsic spectral shape of the quasar is not known and is in fact highly degenerate with the amount of dust extinction.We quantify this systematic uncertainty based on the observed variations of quasar spectral indices of 0.19 dex (Krawczyk et al. 2015).This translates into an uncertainty of 0.033 mag on E(B − V), which should then be scaled by R V to obtain the systematic uncertainty on A V for the different extinction curves, corresponding to roughly 0.1 mag in A V .These systematic uncertainties are given in Table 3.For two targets, we were not able to fit the reddening because the spectra are bluer than the observed template, which would yield a nonphysical negative value for A V .For these targets, we therefore just report upper limits on A V < 0.1 in Table 3 with no extinction law.

Emission-line analysis
We searched for emission counterparts of the ESDLAs by searching for Lyα emission in the UVB spectra and nebular emission lines in the near-IR (NIR) spectra: the [O ii] λλ 3727, 3729 doublet, the [O iii] λλ 4959, 5007 doublet, Hα, and Hβ.
In the case of Lyα, the quasar light is naturally removed by the intervening DLA, so that we can directly search for emission in the DLA core.For emission lines in the NIR spectra, we first need to subtract the quasar continuum emission in order to search for the weak lines from the ESDLA counterparts.The quasar trace in the 2D data was subtracted using a procedure similar to that described by Møller (2000), where we used a local estimate of the spectral point spread function (SPSF) instead of relying on a parametrised model.This empirical SPSF was constructed using a median combination of the SPSF from the A125, page 5 of 22 quasar dominated trace on either side of the expected location of the emission line.We discarded the SPSF from columns in the spectrum within 40 pixels around the line centroid in wavelength space.
For each quasar, we looked for emission using not only the individual 2D spectra, but also the combined spectrum.While the latter is composed of individual 2D spectra with different position angles (hereafter PAs) (hence the spatial axis loses some meaning), the PAs differed by only 24 deg.The combined 2D spectra therefore allow us to obtain a higher S/N.
We detected emission lines associated with three ESDLAs in our sample: J0025+1145, J1143+1420, and J1513+0352.The details regarding the emission counterpart of J1513+0352 detected in Lyα are presented by Ranjan et al. (2018), and the details regarding the two new detections are presented in the following subsections.In these cases, we measured the line fluxes in the quasar-subtracted NIR spectra.For each detection, we extracted a 1D spectrum in an aperture centred at the position of the emission line with a width of twice the full width at halfmaximum (FWHM) of the spatial extent of the emission.Based on the extracted 1D spectra, we fitted the emission line with a Gaussian profile in order to obtain the integrated line flux and the observed line width.The uncertainty on the line flux was estimated by randomly sampling the Gaussian profile according to the best-fit values and their uncertainties.We verified that the estimated uncertainties were consistent with the noise in the residuals from the best fit.
For systems without detected emission lines, we estimated upper limits to the line flux based on 100 random apertures in the quasar-subtracted frame.In all cases, we used an aperture size in velocity space of 300 km s −1 , corresponding to twice the average line-width expected for these emission lines (e.g.Krogager et al. 2013).In the spatial direction, we used an aperture size of twice the measured FWHM of the trace.The noise in the aperture was the added in quadrature, and the resulting 3σ limits are given in Table 2.

J0025+1145
We detected both Hα and [O iii]λ5007 in all individual 2D spectra as well as in the combined 2D spectrum, see Figs. 4 and 5.The Hα line flux in the combined frame was measured to be (4.63 ± 1.43) × 10 −17 erg s −1 cm −2 at an impact parameter of 0.23 ± 0.05 arcsec.
The line flux of [O iii] λ 5007 was measured to be (4.16 ± 0.48) × 10 −17 erg s −1 cm −2 in the combined 2D spectrum.The other emission line of the doublet, [O iii] λ 4959, falls exactly on top of a sky emission line; we were therefore not able to detect this line.We measure the impact parameter along the slit to be ρ = 1.8 ± 0.2 kpc for the PA = −168 deg observations and 2.8 ± 0.8 kpc for PA = 168 deg.Because the two PAs do not differ much, the true impact parameter and relative PA were not constrained from triangulation.In turn, we used the average impact parameter as measured from the combined 2D data, ρ = 1.9 ± 0.1 kpc.This value also corresponds to the most likely impact parameter when a uniformly distributed random angle is assumed between the slit and quasar-galaxy direction (see Ranjan et al. 2018).

J1143+1420
For J1143+1420, we detected the [O iii] λ 5007 emission line in the 2D NIR spectra, as shown in Fig. 6.The integrated flux in the combined spectrum is (2.74 ± 0.35) × 10 −17 erg s −1 cm −2 with the emission centred at z = 2.3232.Because the slit PAs of the two observations differed by only 21 deg, it is not possible to triangulate the exact sky position relative to the quasar.However, we find that in both spectra the emission is observed at a very small impact parameter of ρ = 0.7 ± 0.3 kpc.

Discussion
In this section we discuss the absorption and emission properties of ESDLAs and compare them to those of GRB-DLAs.The samples used for the comparison are described in detail in Sect.3.

Gas properties
In the following, the gas-phase abundance of a species (X) is expressed relative to the solar value as where   a component-wise estimation of neutral hydrogen from saturated Lyα profiles (nor from other Lyman-series absorption).We adopted the same solar abundances as in De Cia et al. (2016), that is, from Asplund et al. (2009), following the recommendations of Lodders et al. (2009) about the choice of photospheric, meteoritic, or average values.We used zinc as a reference element for metallicity for all the systems in our sample because this species is known to be volatile and little depleted into dust grains.In turn, iron is known to deplete heavily, and we therefore used the iron-to-zinc ratio, [Fe/Zn] = log N(Zn ii)/N(Fe ii)−log(Zn/Fe) , as a measure of depletion, as is commonly done in the literature.
A125, page 7 of 22 3), even though the latter was obtained from analysing a composite spectrum of SDSS data with low resolution and low S/N.In the left panel of Fig. 7 we compare the distribution of metallicities for ESDLA to that for GRB-DLAs.The metallicity distribution of ESDLAs is found to be very similar to that of GRB-DLAs.The mean GRB-DLA metallicity is [Zn/H] = −1.23 ± 0.02 and the two-sided KS test between ESDLAs and GRB-DLAs metallicities give P = 0.73.This strengthens the similarity between ESDLAs and GRB-DLAs that has been discussed in the literature, in particular, as recently reported by Bolmer et al. (2019).
The distribution of depletion factors in our ESDLA sample is also compared to that of GRB-DLAs in Fig. 7 (middle panel).The two distributions are found to be similar, with average values of [Fe/Zn] = −0.50 ± 0.05 (ESDLAs) and [Fe/Zn] = −0.56 ± 0.04 (GRB-DLAs).The two-sided KS test gives P = 0.95.Therefore we cannot reject the null hypothesis that the two distributions are the same.

Dust reddening
With both similar column densities and chemical enrichment, we can expect the dust reddening induced by ESDLAs to be similar to that induced by GRB-DLAs.The two distribution are compared in the right panel of Fig. 7.While the two distributions overlap (the KS test gives P = 0.45), ESDLAs tend to have higher dust reddening than GRB-DLAs, with average values of A V = 0.2 and A V = 0.1.Differences between the Notes. (a) Instrument used for the observations: X-shooter (XS) or UVES. (b) Metallicities and depletion factors from the literature have been corrected to use our adopted solar abundances. (c) Extinction measurements are available for systems observed with X-shooter only, except for J0843+0221, for which an estimate has been obtained directly from the SDSS spectrum (Balashev et al. 2017).The best-fit extinction law is indicated in parentheses.Uncertainties are dominated by systematics due to intrinsic quasar-shape variations and hence are all of the order of 0.1 mag. (d) UVES-based measurement from Noterdaeme et al. (2015a). (e) No H 2 line is covered for this system due to the Lyman break of a system at higher redshift; see Noterdaeme et al. (2012) and Kulkarni et al. (2012). (f ) These quasars are bluer than the composite spectrum (see Sect. measured extinctions arising from GRB-DLAs and ESDLAs, if any, are likely not of physical origin.Instead, these can be explained by differences in the measurement methods and observational biases.We note about the extinction measurements themselves that the intrinsic spectral shape of quasars is more complex than those of GRB afterglows and can vary strongly from one quasar to another.In addition, there is a degeneracy between dust at the quasar redshift and dust in the absorber, while for GRB-DLAs, the light source and the absorber virtually have the same redshift.The higher reddening measured for ESDLAs compared to GRB-DLAs is also mostly driven by the absence of GRB-DLAs with A V > 0.3 in the sample from Bolmer et al. (2019).As discussed by these authors (and previously by Ledoux et al. 2009), the extinction by dust in so-called dark bursts (with A V > 0.5) is high enough for the optical afterglow to become too faint for spectroscopic follow-up, implying a bias against dimmed systems in the spectral follow-up campaigns of GRB afterglows.
In principle, a similar bias could also apply to quasars when they are identified based on their colours in a flux-limited sample, as demonstrated by Krogager et al. (2019) for the SDSS-II.However, the following SDSS stages are likely less inclined to dust bias because of improved selection methods and deeper observations.In addition, unlike GRB afterglows, quasars are not transient and therefore do not suffer from a strong targeting bias that affects the decision for follow-up.In conclusion, within the current limitations and taking the uncertainty on the dust measurements of ESDLAs and the absence of dark GRB-DLAs into account, there is no evidence for an intrinsic differ-ence between the dust extinction from ESDLAs and that due to GRB-DLAs.

Molecular hydrogen
Theoretical H i-H 2 transition models show that molecular hydrogen is mostly dependent on three parameters: the total gas column density, the abundance of dust, and the UV flux to density ratio (Krumholz 2012;Sternberg et al. 2014;Bialy & Sternberg 2016).We have shown that ESDLAs and GRB-DLAs have similar column densities, chemical enrichment, and dust properties.Comparing the molecular content of both populations could then provide indications of any possible difference in the physical conditions.For example, from an early small sample of GRB-DLAs without H 2 , Tumlinson et al. (2007) have suggested that the proximity to the explosion site in the case of GRB-DLAs could influence the abundance of H 2 .Ledoux et al. (2009) estimated a distance of 0.5 kpc between the absorber and the GRB explosion by photoexcitation modelling of Fe ii.They also showed that at these distances, the GRB afterglow photons cannot influence the abundance of H 2 and that X-shooter observations might find molecular hydrogen in GRB-DLAs.Bolmer et al. (2019) showed that there is no apparent lack of H 2 in GRB-DLAs, with 6 firm detections in their sample of 22.Moreover, because dusty systems are more likely to harbour H 2 , the absence of dark burst systems may lead to an underestimate of the actual H 2 incidence rate in GRB-DLAs.Quantifying the number for the supposed dark bursts, Bolmer et al. (2019) concluded that the H 2 detection rate in GRB-DLAs could be as high as 41%.
In our total sample of ESDLAs (Table 3), H 2 is detected in 10 out of 18 systems (no limit at all could be placed towards J 1135−0010), which is consistent with the detection rate in GRB-DLAs given the small number statistics and caveats about detection bias in GRB-DLAs.However, the detection limits are different from one system to another.For example, H 2 is firmly detected with log N(H 2 ) ∼ 17.4 towards HE 0027−1836, which has extensively been observed with UVES (see Noterdaeme et al. 2007;Rahmani et al. 2013), while we could only place a weak limit (N(H 2 ) < a few × 10 18 cm −2 ) for several systems observed with X-shooter.A fair comparison therefore cannot be made considering the detection rates at face value.When we only consider systems with log N(H 2 ) > 18.3 (our weakest detection limit), there are 6/14 (∼43%) strong H 2 systems in our homogeneous sample of ESDLAs (7/18 (∼39%) for the total sample).These values agree well with the ∼35% incidence rate of high H 2 column densities among ESDLAs as derived from a stack of SDSS spectra (Balashev & Noterdaeme 2018).Regarding the GRB-DLA sample, 4 (possibly 5, when we consider that GRB 151021A as an H 2 detection) are seen in the sample of 22 GRB-DLAs by Bolmer et al., that is, a lower rate of ∼20%.However, we note that while most GRB-DLAs tend to have very high N(H i) content, there are still 8 GRB-DLAs with log N(H i) < 21.5.When we compare only those GRB-DLAs with log N(H i) > 21.5, we find that GRB-DLAs have a strong-H 2 detection rate of 4-5 out of 14, i.e., about 30%.This is consistent with the ESDLA detection rate given the low number statistics.
To illustrate this further, we plot the location of the GRB-DLAs and ESDLAs in the column density versus metallicity space and identify the systems with strong H 2 lines, see Fig. 8.The two populations seem to be indistinguishable: most of the strong H 2 -bearing systems are located in the upper right region of the figure (high N(H i) and high metallicity, beyond the limit that was originally proposed by Boissé et al. (1998) as due to dust bias).The similarity between the location of H 2 -and non-H 2 -bearing ESDLAs with respect to GRB-DLAs strengthens the conclusions by Bolmer et al. (2019) and Ledoux et al. (2009) that the effect of the GRB explosion itself is not strong enough to significantly dissociate H 2 in the observed GRB-DLAs.In other words, GRB-DLAs most likely arise from gas located within the same galaxy as that of the GRB itself, but they are far enough away from the UV radiation of the explosion site or the associated star-forming region (with distances d > 0.5 kpc) to form H 2 .
Finally, we note that C i lines are detected in about 36% of our sample, while Heintz et al. (2019) observed about 25% in GRB-DLAs.However, this difference is not significant given the small number statistics of both samples, and even more importantly, the detection of C i depends strongly on the resolution and the exact S/N ratio that is achieved.

Kinematics
We compare the kinematics of the ESDLA and GRB-DLA samples last.We chose a commonly used way to quantify the kinematics of the DLAs by means of the velocity width, ∆v 90 , which is defined as the velocity interval comprising 90% of the line optical depth (see Prochaska & Wolfe 1997).We ignored the low-resolution (110 < FWHM < 480 km s −1 ) data for GRB-DLAs because a low spectral resolution leads to a broadening of the absorption lines, which in turn leads to an over-estimation of ∆v 90 .Arabsalmani et al. (2015) showed that QSO-DLAs and GRB-DLAs are indistinguishable in the ∆v 90 -metallicity plane.The authors also showed that in contrast to the general spatial distribution of QSO-DLAs around their host galaxy, the GRB-DLA sightlines pass from a deeper part of the dark matter halo potential well, indicating close proximity to the centre of their associated galaxy.We proceed to perform a similar comparison of GRB-DLAs with our subset of ESDLAs.
For our ESDLA sample, we estimated ∆v 90 using the metal absorption lines with peaked absorption strength ∼0.3 of the normalized flux.The exact lines we used for the measurements are given in Table A.2 in the appendix.At the spectral resolution of X-shooter, this corresponds to lines that are in principle neither saturated nor weak.We used the same correction for the spectral resolution as applied by Arabsalmani et al. (2015).The comparison of these samples is shown in Fig. 9, along with the values for the overall QSO-DLA population by Ledoux et al. (2006).The ESDLA sample is located in the same region of the ∆v 90metallicity plane as the GRB-and QSO-DLA samples.The distribution of the values of ∆v 90 for ESDLAs and GRB-DLAs are shown in the top panel of Fig. 9.The two-sided KS test with the null hypothesis that these two samples are drawn from the same distribution gives P = 0.86.

Associated galaxies
The similarity between the gas properties of GRB-DLAs and ESDLAs (i.e.metallicities, depletion, H 2 content, H i column densities, and kinematics) suggests that the two classes of absorbers probe a similar galaxy population and similar environments within their respective host galaxies.Because the longduration GRBs (which we consider here) are linked to the death of a massive star (e.g.Woosley & Bloom 2006), they are expected to arise at relatively small impact parameters in their host galaxies where young, massive stars are formed.This is indeed confirmed by observations (e.g.Fruchter et al. 2006;Arabsalmani et al. 2015;Lyman et al. 2017).
The QSO-DLAs, on the other hand, are observed to span a wide range of impact parameters (Fynbo et al. 2008;Krogager et al. 2017).However, observations and simulations both indi- For comparison, the grey points show the values for the overall intervening QSO-DLA population from Ledoux et al. (2006).The histograms in the top panel provide the ∆v 90 distributions of the two samples using the same colour-coding as in Fig. 8.
cate an anti-correlation between the typical impact parameter and the hydrogen column density (e.g.Rahmati & Schaye 2014).
It is therefore expected that ESDLAs and GRB-DLAs probe a similar range of impact parameters.Nonetheless, because the very highest column density QSO-DLAs (i.e.ESDLAs) are rare, the impact parameters of ESDLAs and the overlap with GRB-DLAs has been poorly studied until now.
To test this expectation further, we compare the impact parameters for our sample of ESDLAs to those observed for GRB-DLAs, see Fig. 10.For this comparison, we disregarded the ESDLAs that have only been observed with UVES because the UVES data do not provide useful constraints on the associated emission: the UVES slit is not fixed on the sky during integration.In total, emission counterparts are detected for five ESD-LAs (see Sect. 3).For the remaining systems, we only obtained upper limits to the line luminosities within the effective slit aperture.While we cannot formally reject the possibility that the emission from associated galaxies falls outside the region covered by the slit, we note that the impact parameters are very small (< 2 kpc) for all detections.This is well within the range of impact parameters for which we have uniform coverage, determined by half the slit width (i.e.≈6.6 kpc for the median redshift of this sample).However, because the slit is much longer, we do cover impact parameters out to ∼40 kpc along the random slit direction.We therefore conclude that the non-detections most likely arise because the associated emission is below the detection limit (see also Krogager et al. 2017).This is further bolstered by the fact that all five detections have ρ < 3 kpc and none has impact parameters in the range 3 < ρ < 6.6 kpc, where our effective sky coverage is still complete.Our findings are finally also consistent with the suggestion by Noterdaeme et al. (2014) that ESDLAs typically probe impact parameters ρ 2.5 kpc, as evidenced by the detection of Lyα emission in a stack of SDSS fiber spectra and the similarity between the emission properties and Lyα emitting galaxies.
In Fig. 10 we also show the results from simulations by Rahmati & Schaye (2014) in which the very high N(H i) regime is found at very small impact parameters.This is in excellent agreement with our observations.This is not the case at lower column densities, however, where most of the observed data points are located well above the expected median values.This is likely the result of biasing towards high-metallicity systems, which were shown to be statistically located at larger impact parameters than low-metallicity systems (Krogager et al. 2017) because galaxy masses and sizes are correlated with metallicity.
It is also possible that many of the absorber-galaxy associations suffer from an identification bias in which the absorber is associated with the brightest galaxy in a group.Because we detect systems at very small impact parameters, this bias is much less likely to occur in our ESDLA sample.
Considering indeed that the non-detections of galactic emission for several systems is more likely because the line fluxes are lower than our detection limits and not because the galaxy is located outside the area that is covered by the slit, we can estimate upper limits on the star formation rates using the calibration by Kennicutt (1998), SFR/M yr −1 = 7.9 × 10 −42 L(Hα)/erg s −1 . (2) The median upper limit to the star formation rate from Hα ranges from 6 to 30 M yr −1 .Using Lyα instead and assuming standard case-B recombination theory (Brocklehurst 1971;Osterbrock 1989), assuming a low Lyα escape fraction of 5% (Hayes et al. 2011), we obtain more stringent limits of typically 2−5 M yr −1 , which are lower than typically seen for high-redshift GRB-DLA host galaxies.This is expected, however, because GRBs are by selection associated with active star formation (see discussion by Vergani et al. 2015).

Summary
We have presented the analysis of a sample of 11 intervening quasar absorption-line systems selected for their very high H i column densities and observed with the X-shooter spectrograph at the VLT.One system has been presented by Ranjan et al. (2018), and we here measured the metal abundances, molecular and dust content, and gas kinematics, and searched for the associated galaxy counterparts for the remaining 10 systems.We clearly detect molecular hydrogen in 5 out of 11 systems, 3 of which are reported here for the first time.For the remaining systems, the presence of H 2 is hard to ascertain, and we instead provide conservative upper limits from the maximum H 2 column density that are still consistent with the X-shooter data.We supplemented our ESDLA sample with measurements from the literature and compared their absorption properties with those of GRB-DLAs.
We find that the two populations are indistinguishable in terms of chemical enrichment (metallicity and dust depletion) and there is no marked difference in the conditions for the presence of H 2 .This supports the suggestion by Noterdaeme et al. (2014) that ESDLAs likely arise from similar environments as GRB-DLAs and the suggestion by Bolmer et al. (2019) that GRB-DLAs are likely located far enough from the explosion site for their H 2 content to be little affected by the GRB explosion itself.
The similarity between ESDLAs and GRB-DLAs is further bolstered by the similar absorption-line kinematics and by the  Bouché et al. 2013;Fynbo et al. 2011Fynbo et al. , 2013;;Hartoog et al. 2015;Kashikawa et al. 2014;Krogager et al. 2012Krogager et al. , 2013Krogager et al. , 2017;;Ma et al. 2018;Møller & Warren 1993;Møller et al. 2004;Neeleman et al. 2018;Noterdaeme et al. 2012;Srianand et al. 2016;Ranjan et al. 2018;Rudie et al. 2017;Warren et al. 2001;Weatherley et al. 2005;Zafar et al. 2017), and the two new detections presented here are shown as filled circles.Values corresponding to GRB-DLAs are plotted as red squares.The grey line with the shaded region shows the predicted median and the 68% region from simulations by Rahmati & Schaye (2014).Right panel: zoom-in around the high column density region.The extent of the y-axis corresponds to the impact parameter out to which we uniformly cover the area around the quasar, which is determined by the slit width of our observations.The dashed green rectangle shows the region of typical impact parameters derived by Noterdaeme et al. (2014) for log N(H i) > 21.7 systems (see text).observed impact parameters.We detect emission lines associated with 3 systems out of our sample of 11 systems, two of them being reported here for the first time This brings the number of emission counterparts obtained at the very high N(H i)end of the quasar-DLA column-density distribution to a total of five.All emission counterparts of ESDLAs are located at very small impact parameters (ρ < 3 kpc), when much larger impact parameters are generally reported for the overall population of DLAs.Because our observations are complete out to ρ ∼ 7 kpc, we argue that the counterparts for the remaining ESDLAs are more likely below the detection limit and are not located outside the area covered by the slit.
Long-duration γ-ray bursts have been shown to be tracers of star formation throughout the history of the Universe (Vergani et al. 2015).While GRB-DLAs then probe neutral gas inside high-redshift, highly star-forming galaxies, ESDLAs may act as completely blind analogues, probing neutral gas in the heart of high-redshift galaxies without any prior on the instantaneous star-formation rate.In the near future, blind radio absorption line surveys such as the MeerKAT Absorption Line Survey (MALS, Gupta et al. 2016) will further explore the population of neutral gas absorbers without the effects of dust biasing that currently affect DLA samples.

Fig. 2 .
Fig. 2. Comparison of H i column densities measured in our sample of 11 systems observed with X-shooter to those originally obtained from the low-resolution, low S/N SDSS data.The horizontal bar in the upper left corner shows the typical uncertainty of SDSS measurements.The dashed line shows the one-to-one relation, with dotted lines showing ±0.1 dex around this relation.

Fig. 3 .
Fig.3.Velocity structure of metal lines tracing the overall neutral gas (top panels) compared with that seen in H 2 (bottom panels), where the transition naming is "L" for Lyman, followed by the band number (vibrational level of the upper state), the branch (corresponding to the selection function), and the rotational level of the lower state in parentheses.Downward arrows indicate the positions of the different components in the best-fit model.The absorption seen at +150 km s −1 (−150 km s −1 ) in the L1R(0) (L2R(1)) panel is not an additional component but absorption from L1R(1) (L2R(0)).The zero of the velocity scale is set relative to the redshift of the strongest metal component.

Fig. 4 .
Fig. 4. Detection of Hα emission associated with the z abs = 2.304 ESDLA towards J0025+1145.Top panels: three individual 2D spectra and the combination of all of them.The quasar trace has been removed over the central region (∆v ∼ 3000 km s −1 ).The Hα emission corresponds to the remaining blob in the centre of the images, and is quite evident in the combined 2D spectrum.The side panels show the spatial extent of the trace (black) and of the Hα emission (red) obtained by collapsing the 2D data along the wavelength axis.The location of their centroid (horizontal dotted lines) provides a measure of the impact parameter along the slit direction.The fifth panel shows the quasar-subtracted 1D data (in blue) summed over the spatial axis together with a Gaussian fit (in red) to the emission line.The normalised absorption profile of Fe iiλ2260 is shown in the bottom panel for comparison.The absorption seen at v ∼ −1600 km s −1 is due to Fe iiλ2249.

Fig. 5 .
Fig. 5. Same as Fig. 4 for the [O iii]λ5007 emission line at the redshift of the z abs = 2.304 DLA towards J0025+1145.

Fig. 8 .
Fig. 8. Metallicity vs. H i column density for the sample of ESDLAs (homogeneous, blue) and GRB-DLAs (red).Filled symbols represent the systems with firm detection of strong H 2 lines (log N(H 2 ) ≥ 18.3).The dashed line shows the constant total metal column density corresponding to log N(Zn ii) = 13.15 from Boissé et al. (1998).

Fig. 9 .
Fig. 9. Comparison of the metallicities and ∆v 90 values of ESDLAs (blue circles) and GRB-DLAs (red squares) observed with X-shooter.For comparison, the grey points show the values for the overall intervening QSO-DLA population fromLedoux et al. (2006).The histograms in the top panel provide the ∆v 90 distributions of the two samples using the same colour-coding as in Fig.8.

Fig. 10 .
Fig.10.Impact parameter, ρ, vs. the neutral hydrogen column density, log N(H i), measured along the line of sight.Systems associated with the generic population of QSO-DLAs are plotted as blue circles (fromBouché et al. 2013;Fynbo et al. 2011Fynbo et al. , 2013;;Hartoog et al. 2015;Kashikawa et al. 2014;Krogager et al. 2012Krogager et al. , 2013Krogager et al. , 2017;;Ma et al. 2018;Møller & Warren 1993;Møller et al. 2004;Neeleman et al. 2018;Noterdaeme et al. 2012;Srianand et al. 2016;Ranjan et al. 2018;Rudie et al. 2017;Warren et al. 2001;Weatherley et al. 2005;Zafar et al. 2017), and the two new detections presented here are shown as filled circles.Values corresponding to GRB-DLAs are plotted as red squares.The grey line with the shaded region shows the predicted median and the 68% region from simulations byRahmati & Schaye (2014).Right panel: zoom-in around the high column density region.The extent of the y-axis corresponds to the impact parameter out to which we uniformly cover the area around the quasar, which is determined by the slit width of our observations.The dashed green rectangle shows the region of typical impact parameters derived by Fig. A.12. Same as Fig. A.11 for the z abs = 2.444 towards J 1258+1212.
Fig. A.14. Portion of the X-shooter spectrum of J 0025+1145 covering the L0-band absorption lines of H 2 from the DLA at z = 2.304(2).The rotational levels J are indicated above each blue tick mark.The normalised spectrum is shown in black, and the synthetic profile with total H 2 column density log N(H 2 ) = 20.1 is over-plotted in red.

Table 2 .
Ranjan et al. 2018) of ESDLAs from our X-shooter observing programme.All measurements are from this work, except J1513+0352 (fromRanjan et al. 2018).The upper limits correspond to the mean detection limit in the spatial area covered by the slit (see text). Notes.

Table 3 .
Main absorption properties of extremely strong DLAs with VLT spectroscopic follow-up.

Table A .
1. Total column densities of low-ionization metal species in extremely strong DLAs.