\section{Introduction}
\label{intro}

Molecular hydrogen (H$_2$) is by far the most abundant molecule in the Universe and predominantly resides in the cold neutral gas of galaxies. This cold phase may be subject to Jeans instability and is hence thought to be tightly linked to the onset of star formation. 
While the association between star formation and molecular gas 
is well established on kpc-scales for normal galaxies in the local Universe \citep[e.g.][]{Bigiel2008, Schruba11}, observations in emission in the distant Universe have much coarser spatial resolution and are mostly limited to the bright-end (and much probably high-metallicity end) of the galaxy population. In addition, the mass of molecular gas is estimated indirectly via tracers, like CO or dust emission, which are metallicity-dependent \citep[e.g.][]{Genzel12}.
While the link between star formation and molecular gas has been observed in a variety of environments, there are some observational indications that star formation could proceed directly from the atomic gas when the metallicity is low \citep[e.g.][]{Michalowski2015}. Indeed, in the diffuse ISM H$_2$ forms on the surface of dust grains, whose abundance scales with metallicity. Therefore, it is expected that at low metallicities, the H$_2$ formation timescale can be larger than the gravitational collapse time \citep[see e.g.][]{Krumholz2012}.

The star formation process can be studied not only in emission (where the relation between star formation and atomic gas is usually derived over the whole extent of the galaxy), but also through absorption-line studies. 
This is particularly useful for high-redshift galaxies, since the detection in absorption is independent of the luminosity of the associated object. 
Blind absorption-line studies indeed select systems owing to their cross-section only. 
The neutral gas in and around high-redshift galaxies in absorption studies is mostly detected via the \ion{H}{i}-Ly$\alpha$ absorption line, producing the so-called damped Lyman-$\alpha$ systems \citep[DLA,][]{Wolfe2005}, when the \ion{H}{i} column densities reaches $\log N(\ion{H}{i}) \gtrsim 20$, ensuring the gas is predominantly neutral. Statistical studies have shown that DLAs contain the bulk of the neutral atomic gas in the Universe \citep[e.g.][]{Noterdaeme2012c} and present a significant metallicity evolution with redshift \citep[e.g.][]{DeCia2018}, i.e., the smoking gun of star formation.
However, despite recent progress, it remains technically challenging to detect galaxy counterparts of DLAs in emission \citep[see][and references therein]{Krogager2017}, and therefore to understand the origin of DLAs and constrain the associated star-formation rates.
However,
absorption studies of cold and molecular gas already provide very useful information to understand the star-formation processes in high-redshift galaxies. Indeed molecular hydrogen is found in a fraction of high-$z$ DLAs \citep[$\sim 4\%$ for $\log N(\rm H_2) \gtrsim 18$,][]{Balashev2018} meaning that cold gas is intercepted by the line-of-sight. Such H$_2$-bearing DLAs potentially provide information on the transition from atomic to molecular gas and subsequently on the onset of star formation.

Locally, observations towards nearby stars have revealed that the H$_2$ column density significantly increases when the \ion{H}{i} column density reaches some threshold value.
This threshold directly traces the \ion{H}{i}-to-H$_2$ transition and is well-described by analytic and numerical modelling \citep[see][and references therein]{Sternberg2014}. The \ion{H}{i}-to-H$_2$ transition is also found to occur at higher \ion{H}{i}\ column density in the Magellanic clouds, owing to lower metallicities and stronger UV fields, in agreement with expectations from the models.
In turn, a global understanding of the \ion{H}{i}-to-H$_2$ transition at higher redshifts is still missing
\citep[see][]{Noterdaeme2015b}, mostly because of small statistical samples, which actually probe a wide range of physical conditions.
To date, there are $\sim 40$ confirmed intervening H$_2$-bearing DLAs detected at redshifts $z>2$.

Before the advent of the Sloan Digital Sky Survey \citep[SDSS;][]{York2000}, H$_2$ at high redshifts has been blindly searched, mostly in high-resolution DLA spectra from archival data. The detection rate is found to be quite low: \citealt{Ledoux2003, Noterdaeme2008} found H$_2$ in about 10\% of the DLAs surveyed using high-resolution spectra, while \citealt{Jorgenson2014} found only one H$_2$ in 86 blindly-selected DLAs, using mid-resolution data. \citet{Petitjean2006} proposed that 
above a metallicity $Z>0.2Z_{\rm \odot}$, the probability to detect H$_2$ increases 
to about 50\%. However, this is of little practical use since the metallicity is known accurately only once a high-resolution spectrum is obtained.
With a huge increase in the number of quasar spectra available, the SDSS spectroscopic database has presented a unique opportunity to search for H$_2$ systems. Different techniques have been proposed to pre-select these from low-resolution, low-S/N spectra. 

A first efficient pre-selection was based on the presence of neutral carbon lines \citep[e.g.,][]{Srianand2008, Noterdaeme2010}. Indeed, \ion{C}{i}\ is a fragile species that survives more easily in the same conditions as H$_2$ requires \citep[see e.g.][]{Liszt1981,Ge2001}. Since the strongest \ion{C}{i} absorption lines are usually redshifted outside of the Lyman-$\alpha$ forest, it was possible to systematically search for \ion{C}{i}\ in SDSS spectra without any prior on \ion{H}{i}\ content \citep{Ledoux2015}. 
While the sample of \ion{C}{i}-selected molecular systems is continuously increasing thanks to the SDSS \citep[e.g.][]{Ma2015, Noterdaeme2017}, the majority of the systems have $z\lesssim2$, so that H$_2$ (and sometimes \ion{H}{i}) lines are inaccessible from the ground, due to the atmospheric cutoff at $\sim3000~${\AA}. For those systems where H$_2$ is covered, \citet{Noterdaeme2018} find a 100\% detection rate of H$_2$ molecules.
However, it has been found that a selection based on \ion{C}{i}\ only identifies metal-rich systems (with $Z\sim Z_{\odot}$, see \citealt{Zou2018} and Ledoux et al., in prep.).

Another method to pre-select possible H$_2$-bearing systems is to consider extremely-strong DLAs (ESDLAs, with $\log N(\ion{H}{i}) > 21.7$, \citealt{Noterdaeme2014}), since a large critical column density for atomic-to-molecular conversion is expected for clouds at the typical metallicity of DLAs.
Such a conversion could furthermore explain the steeping slope of the observed $N(\ion{H}{i})$-distribution function at high $N(\ion{H}{i})$ \citep{Zwaan2006,Altay2013}. From a galactic-scale point-of-view, we also expect 
the corresponding lines-of-sight to probe gas at small galacto-centric distances, as confirmed by the detection of Lyman-$\alpha$ emission in composite ESDLA spectra \citep{Noterdaeme2014} and in individual systems \citep{Noterdaeme2012a, Ranjan2018}.
This picture is corroborated by the possible increase of thermal pressure with increasing \ion{H}{i}\ column density \citep{Balashev2017}, indicating that such lines-of-sight probe more central regions of galaxies.
An enhanced H$_2$ detection rate in ESDLAs has indeed been confirmed through follow-up of ESDLAs \citep[][]{Noterdaeme2015b,Ranjan2019} as well as in stacked SDSS spectra \citep{Balashev2018}. 

Finally, \citet{Balashev2014} proposed to directly identify H$_2$ Lyman-Werner bands in SDSS spectra. However, this is a complex task because of the low spectral resolution, the typically low signal-to-noise ratio, and the presence of Ly$\alpha$ forest absorption lines.
Due to these limitations, only a few percent of SDSS spectra are suitable for the detection of H$_2$, and only systems with high H$_2$ column densities ($N(\mbox{H$_2$}) \gtrsim 10^{19}$~cm$^{-2}$) can be detected.
In our original search of SDSS DR9, we have found about 25 good H$_2$ candidates (with low false detection probability (FDP), see \citealt{Balashev2014}); A number that has since then been increased to about 100 using SDSS DR14.

In order to study the H$_2$ absorbers selected directly from SDSS in greater detail, we have initiated a follow-up program using the Very Large Telescope (VLT). The higher spectral resolution and signal-to-noise ratio allow us to confirm the presence of H$_2$ in the candidate systems and to measure the gas properties.
One of the H$_2$-bearing DLAs at $z=2.786$ towards J\,0843+0221 has already been published owing to its unique and extremely high H$_2$ column density of $\log N(\rm H_2)\sim 21.2$ \citep{Balashev2017}.
In the present paper, we describe the full analysis of the remaining 7 systems from this follow-up sample.

The paper is structured as follows. First, we describe the observations and the methodology that we have used to analyze the H$_2$-bearing DLAs in Sect.~\ref{sect:observations} and \ref{sect:analysis}, respectively. In Sect.~\ref{sect:discussion}, we discuss the properties of the sample, which we compare to other H$_2$ samples arising from different selection methods.
Lastly, in Sect.~\ref{sect:conclusion}, we offer a summary of the findings presented in this paper.

\section{Observations and data reduction}
\label{sect:observations}
Since the main goal of our pilot observations was to confirm the reality of the H$_2$ absorption systems detected by \citet{Balashev2014}, we first selected the targets with lowest false H$_2$-detection probability from our list of DLA systems (extended to DR10) and that were accessible from the Paranal observatory during the semester of observations.
The sample was observed with the X-shooter spectrograph \citep{Vernet2011} mounted at the Unit Telescope 2 of the VLT (VLT-UT2, Kueyen). The observations were performed in service mode under program ID 094.A-0362 (PI: Balashev) between the end of 2014 and beginning of 2015.  The log of observations is given in Table~\ref{table:jour_obs}. For all observations, we used a slit width of 1 arcsec in the UVB arm, and 0.9 arcsec in the VIS and NIR arms, leading to a resolving power of $R\sim5400$, $\sim$8900, and $\sim$5600\footnote{The actual value of the spectral resolution depends on the atmospheric seeing value, which in some cases was smaller than the slit width, see Table~\ref{table:jour_obs}.} in the UVB, VIS, and NIR arms, respectively. For each target, we obtained two 1-hour exposures except for SDSS\,J234730.76$-$005131.68, which was observed only once before the end of the observations.

\begin{table}
	\centering
	\addtolength{\tabcolsep}{-2pt}
	\caption{Journal of X-shooter observations. 

	}
	\label{table:jour_obs}
	\begin{tabular}{lcccc}
		\hline
		Quasar & Date & Seeing & Airmass & Exp. times$\dagger$ \\
		       &      & ($\arcsec$) &    & (ks)                \\
		\hline
		J\,013644.02+044039.00 & 27.11.2014 & 1.15 & 1.16 & 3.6, 3.6, 3$\times$1.2\\
        J\,013644.02+044039.00 & 28.11.2014 & 1.15 & 1.26 & 3.6, 3.6, 3$\times$1.2\\
        J\,085859.67+174925.32 & 24.12.2014 & 0.73 & 1.36 & 3.6, 3.6, 3$\times$1.2\\
        J\,085859.67+174925.32 & 15.02.2015 & 1.38 & 1.37 & 3.6, 3.6, 3$\times$1.2\\
        J\,090609.46+054818.72 & 23.12.2014 & 0.79 & 1.26 & 3.6, 3.6, 3$\times$1.2\\
        J\,090609.46+054818.72 & 25.12.2014 & 0.67 & 1.22 & 3.6, 3.6, 3$\times$1.2\\
        J\,094649.47+121628.56 & 18.02.2015 & 0.83 & 1.27 & 3.5, 3.5, 3$\times$0.9\\
        J\,094649.47+121628.56 & 18.02.2015 & 0.99 & 1.29 & 3.5, 3.5, 3$\times$0.9\\
        J\,114638.95+074311.28 & 18.02.2015 & 0.79 & 1.19 & 3.6, 3.6, 3$\times$1.2\\
        J\,114638.95+074311.28 & 24.02.2015 & 0.77 & 1.25 & 3.6, 3.6, 3$\times$1.2\\
        J\,123602.11+001024.60 & 14.02.2015 & 0.87 & 1.13 & 3.6, 3.6, 3$\times$1.2\\
        J\,123602.11+001024.60 & 18.02.2015 & 1.15 & 1.12 & 3.6, 3.6, 3$\times$1.2\\
        J\,234730.76$-$005131.68 & 20.12.2014 & 0.52 & 1.39 & 3.6, 3.6, 3$\times$1.2\\
		\hline
	\end{tabular}
		\addtolength{\tabcolsep}{+2pt}
	\begin{tablenotes}
	\item $\dagger$ Exposure times given for the UVB, VIS, and NIR arms, in this order.
	\end{tablenotes}
\end{table}

We first reduced the spectra using the ESO Reflex pipeline \citep{Freudling2013}.
However, the unavailability of atmospheric dispersion correction at the time of observations led to a significant wavelength-dependent trace offset in the UVB, which was not properly taken into account by the pipeline.
Therefore, we used our own code to extract the 1-dimensional (1D) spectra along the trace starting from the 2D pipeline product (which incorporates dark/bias correction, flat-field correction, 2D geometry and arc-line solutions).
In addition, we improved the removal of cosmic-ray hits by flagging the corrupted pixels on the 2D exposures interactively. The 2D sky model was constructed using polynomial fitting along the spatial axis for each dispersion column and taking into account the rejection of outliers and correction of 2D geometrical issues at the edge of bright sky lines.
We used the optimal-extraction method with a Moffat profile to extract the 1D spectrum \citep{Horne1986}. Finally, we corrected the flux error in each pixel of each spectra by a wavelength-dependent factor so that it matches the observed dispersion relative to the quasar continuum. To do this, we identified unabsorbed regions in the continuum and estimated a smooth average of the dispersion of pixels using 5-10~\AA-wide sliding windows. This extraction resulted in an overall $\sim$25\% improvement in the S/N ratio and a better handling of lower quality regions, which were otherwise hard to manage through the pipeline.

\section{Analysis}\label{sect:analysis}

Before presenting the analysis of the different systems, we note that two out of seven DLAs studied in this paper have redshifts close to that of the quasar ($\Delta v \lesssim 3000$\,km/s) and could hence be considered as proximate. This is quite surprising at the first glance since the parent sample used by \citet{Balashev2014}, i.e., the DLA catalogue from 
\citet{Noterdaeme2012c}, was based on a search optimised for intervening DLAs and hence likely biased {\sl against} 
proximate DLAs. 
Recently, \citet{Noterdaeme2019} found that the incidence rate of proximate H$_2$-bearing DLAs (within $\sim 2000$~km\,s$^{-1}$ from the quasar) per redshift path unit can be a factor of 5 higher than that of intervening systems. However, the observed number of systems depends on the total redshift path probed by each sample and therefore intervening \mbox{H$_2$}\ systems should still dominate in numbers over proximate ones, due to the much longer total pathlength. The relatively high number of proximate DLAs here is more likely to originate from selection method used for this follow-up sample. Indeed, the objects were selected upon their probability of false \mbox{H$_2$}-detection. This probability is strongly correlated with the number of the band used and therefore, the closer DLA to the quasar, the more H$_2$ bands are available in the spectrum before blue cutoff, and the smaller can be the false detection probability. In this work, as we concentrate on the chemical and physical properties of the systems and found no striking difference between those within a few thousands km\,s$^{-1}$\ from the quasar and those located further in velocity space, therefore we consider them indistinctly.

\subsection{Absorption lines}
\label{sect:alf}
We analysed all DLA systems in our sample in the same manner. We used a custom routine\footnote{We used this routine in recent analyses of absorption-line systems \citep[e.g.,][]{Balashev2016, Balashev2017} and also, partly, in \citet{Noterdaeme2018, Krogager2018b}, where we found that it usually is in good agreement with other codes used for Voigt-profile fitting, like vpfit \citep{Carswell2014} and VoigtFit \citep{Krogager2018a}.} for Voigt-profile fitting to derive the redshifts, column densities and Doppler parameters of the velocity components for \ion{H}{I}, H$_2$, \ion{C}{I}, and  various low-ionization metal absorption lines (e.g. \ion{Si}{ii}, \ion{Fe}{ii}, \ion{S}{ii}, \ion{Zn}{ii}). 
For each absorption line (except \ion{H}{I} Ly$\alpha$, described below) we derived the local continuum using a B-spline interpolation matched on the adjacent unabsorbed spectral regions. 
The spectral pixels used to constrain the fit were selected visually and iteratively. For each quasar we fit available exposures jointly without combining them.
The number of velocity components was also defined visually and increased in case of remaining structure in the residuals.
We used Monte-Carlo Markov Chain calculations with an affine invariant sampler \citep{Goodman2010} to obtain the posterior probability density function (PDF) of fitting parameters (assuming flat uninformative priors). While we directly obtained the PDF, in the following we use the standard way of reporting the best fitting parameters and uncertainties. The best fitting value corresponds to the maximum posterior probability and the estimated uncertainties correspond to the 68.3\% confidence interval around this value (corresponding to formal 1$\sigma$ uncertainty for Gaussian PDF). 

We measured the \ion{H}{i}\ column density from the Ly$\alpha$ line. However, since this line is strongly damped and 
frequently located close to the Ly$\alpha$ emission line, we simultaneously fitted the absorption profile and the unabsorbed continuum modelled usually by 4-8 order Chebyshev polynomials, depending on the wavelength range.
The values all agree within 0.2~dex from those measured at lower S/N and resolution power from the SDSS data.



We confirm the presence of H$_2$ in all systems through the detection of Lyman- and Werner-band absorption lines. {We only detect lines from several low rotational levels of the $\nu=0$ vibrational ground state. Therefore, in the following, we use an abbreviated notation for those levels, as follows: $\rm H_2Ji$, where index $i$ is the rotational quantum number.} While an increase of the Doppler parameter with increasing rotational level has been observed for several H$_2$ systems 
\citep{Lacour2005, Noterdaeme2007, Balashev2009, Ivanchik2010}, the spectral resolution of X-shooter is not high enough to study this effect. Therefore, we tied the Doppler parameters of different rotational levels when fitting H$_2$ lines. 
This implies that the reported column densities for high ($J>2$) rotational levels should be considered with caution, in particular because the corresponding lines are generally in the intermediate regime. 
Notwithstanding, this has a little effect on the measurement of the total $N($H$_2)$ which is dominated by the first rotational levels whose lines are damped and hence not sensitive to the exact Doppler parameter. Additionally, for each exposure we included two additional independent parameters when fitting H$_2$ lines, to linearly correct for the wavelength distortions caused by the unavailability of the ADC \citep[see e.g.][]{Noterdaeme2017}. We then used the obtained values of the zero point and the slope of distortion to correct the spectrum before fitting the \ion{C}{i} and low-ionization metal lines. Finally, since the seeing was narrower than the slit width for some exposures, we also varied spectral resolution for each exposure as independent fitting parameter. 
We found lower H$_2$ column densities than initially estimated from the SDSS in all cases, see Fig.~\ref{compNH2}. The reasons for this is that H$_2$ profiles are unresolved in SDSS, the blends with Ly$\alpha$ forest lines are hard to take into account and the signal to noise ratio are low. For these reasons, the SDSS-based \mbox{H$_2$}\ column densities 
were considered as indicative only by \citet{Balashev2014}. 
Notwithstanding, the column densities agree within 0.3~dex for 4 out of 7 DLAs (5/8 including J0843+0221). From the X-shooter data, we could see that these are indeed the cases where the H$_2$ profiles dominated by one main component, while for the other three DLAs with larger differences, there are at least two \mbox{H$_2$}\ components with comparable column densities. 

\begin{figure}
    \centering
    \includegraphics[width=\hsize]{Ncomp.eps}
    \caption{Comparison of \mbox{H$_2$}\ column densities automatically determined from the low-resolution low S/N SDSS data with those obtained from medium/high-resolution  observations. The sample studied in this paper is shown as filled symbols (all using X-shooter but J\,0843+0221 observed with UVES), while open symbols represent other known \mbox{H$_2$}\ systems (all observed with UVES but towards J81240.69+320808.52, observed with HIRES on Keck, see \citealt{Jorgenson2009}).}
    \label{compNH2}
\end{figure}{}




To fit the low-ionization metal lines, we tied the redshifts and Doppler parameters of various species since the broadening is dominated by turbulence which is independent of the atomic mass. We also varied spectral resolution for each exposure as independent fitting parameter. 
For all DLAs we detected associated \ion{C}{I} absorption lines at the position of H$_2$ components. 
Therefore, we fitted the \ion{C}{I} lines (and the associated fine-structure levels) independently from other metals.

The detailed fitting results, the corresponding figures and specific comments for each DLA in our sample are presented in
Supplementary materials.
The total column densities of \ion{H}{i},  H$_2$, \ion{C}{i}, the average molecular fraction, metallicity, and Fe depletion are summarised in Table~\ref{table:fit_summary}.
We also add the H$_2$-bearing DLAs towards J\,0843$+$0221 to our sample, since it was selected from SDSS using the same technique. However it was obtained with separate high resolution UVES observations, which confirmed the exceptionally high H$_2$ column density \citep{Balashev2017}.

\subsection{Dust-extinction measurements}

We obtained the estimates of $E(B-V)$ by fitting a quasar spectral template to the data assuming a fixed extinction curve. Since we observe no indications of the prominent 2175~\AA\ dust feature, we used the extinction curve observed towards the Small Magellanic Cloud (SMC) as parametrized by \citet{Gordon2003}, however, we modified the far UV curvature part using $c_4 = 0$ since this provided a better description of the data. The dust reddening is assumed to be caused solely by the H$_2$-bearing DLA. We used the quasar template derived by \citet{Selsing2016} and fitted only the data in bona fide continuum regions, i.e., regions which are not strongly affected by absorption or broad emission lines. In the near-infrared, we performed a 5-$\sigma$ clipping in order to discard outlying pixels introduced by the removal of skylines in the data reduction. Moreover, we masked out pixels that are affected by telluric absorption.

Before fitting, the template by Selsing et al. was smoothed with a Gaussian kernel ($\sigma$ = 7~pixels) to prevent the noise in the template to falsely fit noise peaks in the real data. In order to take into account the uncertainty on the template, we subsequently convolved the errors on the observed data with the uncertainty estimate for the template.

We then fitted the template to the data using 3 free parameters: the amount of dust in the absorber, $A_V$, and the change in the intrinsic power-law index, $\Delta\beta$, as well as a flux scaling, $f_0$, since we do not know the intrinsic brightness prior to reddening. The best-fit $A_V$ are given in Table~\ref{table:fit_summary}. We found that power law indices, $\Delta\beta$, to be between $-0.8$ and $0.2$ for this sample. The typical statistical uncertainties on $A_V$ of such measurements are $\sim 0.01$\,mag. However, there is a strong degeneracy between $A_V$ and the intrinsic power-law shape of the quasar spectrum. 
We estimated the systematic uncertainty related to the unknown spectral shape by assuming the inferred Gaussian dispersion for $\beta$ of $\sigma_{\beta} = 0.19$ \citep{Krawczyk2015}. We drew a random $\Delta\beta$ from a Gaussian centred at $\Delta\beta = 0$ with dispersion $\sigma_{\beta}$ and change the quasar template accordingly. We then refitted the spectra with the newly altered, random template to obtain a distribution of $A_V$.  Based on this distribution, we inferred an estimate of the systematic uncertainty of $\sigma_{\rm sys} \sim 0.1$\,mag. We therefore consider this systematic uncertainty as a conservative estimate of the error on individual measurements.

\subsection{Physical conditions} 
\label{sect:phys_cond}
In this section, we derive the temperature and number density of the H$_2$-bearing gas in our sample as well as the intensity of the UV irradiation onto the clouds. In cases where more than one component is detected in H$_2$, we use only the strongest component since \ion{C}{i} is not always detected for the weaker component. We note however that for one absorber (toward J1146+0743) we are able to properly constrain \ion{C}{i} and H$_2$ in both components, and for both components, we infer consistent results for temperature, density and UV intensity. Also for J\,0843+0221 we used the total column densities of H$_2$ and \ion{C}{i}, since absorption lines from $J=0,1$ H$_2$ rotational levels damped and two main components are not resolved.  

We used the $J=0$ to $J=1$ rotational excitation temperature of H$_2$, $T_{01}$, as a proxy for the kinetic temperature $T_k$ in the cold diffuse medium \citep[see][]{Roy2006}. The inferred temperatures for the individual systems are given in Table~\ref{table:fit_summary}. We note that these estimates of the temperature only corresponds to the temperature of the H$_2$-bearing gas and not the overall DLA gas, for which there most likely is a contribution from a warm neutral phase with $T \gtrsim 10^3$\,K \citep[e.g.,][]{Srianand2012, Neeleman2015}.

In order to constrain the number density and the UV flux in the H$_2$-bearing medium, we used the relative population of \ion{C}{I} fine structure levels \citep[see][]{Silva2002}, H$_2$ rotational levels and the analytical theory of \ion{H}{i}-H$_2$ transition from \citet{Sternberg2014}. For \ion{C}{I}, we used the same calculations and code as in \citep{Balashev2017, Ranjan2018}. Here we briefly outline the assumptions used in the analysis. For \ion{C}{i} excitation, we assumed a homogeneous medium, where \ion{C}{i} fine-structure levels are populated by the Cosmic Microwave Background (CMB) radiation, UV pumping and collisions. The CMB temperature at the absorbers redshift is fixed to be $T_{\rm CMB}(z) = T_{\rm CMB}(z=0)(1+z)$ \citep[see, e.g.,][]{Srianand2008,Noterdaeme2011}.
Since the observed H$_2$ column densities are high, we assume that the gas phase containing the bulk of \ion{C}{i} is fully molecularized \citep[e.g.,][]{Balashev2015}, i.e. the main collisional partners are H$_2$ and He. 
Due to the relatively low measured \ion{C}{i} column densities we neglect self-shielding in \ion{C}{i} lines in which UV pumping occurs. 
{Since the collisional coefficients for \ion{C}{i} depend on temperature, we fixed the temperature to the measured $T_{01}$ value, reducing the free parameters to be the hydrogen number density ($n_{\rm H}$) and the intensity of the UV field ($I_{\rm UV}$, that we express relative to the Draine field). The exact shape of the spectrum cannot be constrained, since it depends on the unknown galaxy's stellar population and ambient conditions. Therefore, we assumed Draine shape for the UV background spectrum. We also added the UV metagalactic background in the calculation. However, at redshifts $z\sim 2-3$, it is $\sim$10$\div$30 times lower than the Draine field in the considered wavelength range \citep[see e.g.][]{Khaire2019}.}

{The population of the H$_2$ rotational levels also provides constraints on the ($I_{\rm UV}$, $n_{\rm H}$)-plane, and were previously used to estimate the physical conditions in a few high-$z$ H$_2$-bearing DLAs \citep[e.g.,][]{Klimenko2016, Shaw2016, Noterdaeme2017, Rawlins2018}.}
While UV pumping plays a direct role in the excitation of high rotational level, we are not confident that $J\ge3$ rotational levels are well constrained (see note in Sect.~\ref{sect:alf}). We hence here focus on the excitation of the $J=0,1,2$ levels only. These levels are predominantly thermalised and their excitation typically close to the kinetic temperature \citep[see][]{LePetit2006}, which is set
by the thermal balance, itself being a function of density and UV field. 
Notwithstanding, UV pumping could also contribute to the excitation of these low-J levels,
so that the assumption of a homogeneous cloud is not correct anymore (the resonant H$_2$ lines --at which pumping occur -- quickly saturate inside the cloud). We therefore used the 
PDR Meudon code \citep{LePetit2006} which performs full radiative transfer calculations in H$_2$ lines. 
We used grids of pseudo-spherical constant-density models with the only parameters being the metallicity (assuming dust scales to metallicity), density and intensity of UV field. {The temperature profile in the cloud was calculated from the solution of thermal balance for each model.}
We then ran grids of models with ten logarithmically-spaced points for both the number density and UV field intensity in the respective ranges of $0 \leq \log \rm n_{H} / cm^{-3} \leq 4$ and $-1.5 \leq \log I_{\rm UV} \leq 2.5$ (where $I_{\rm UV}$ is measured in Draine field units). For each system, we chose the appropriate grid corresponding to the observed metallicity. 
We used models assuming a slab of gas with thickness corresponding to $A_V=0.5$ irradiated from both sides. To compare the modeled and observed column densities of different H$_2$ rotational levels, we cut the model at a depth in the cloud where half of the total H$_2$ column density equals half the observed one. This approach mimics a cloud illuminated on both sides with a total column density of H$_2$ equal to the measured one.

We then compared the modelled excitation of H$_2$ rotational levels to the observed ones, including 
a factor of 2 (0.3 dex) uncertainty for the observed H$_2$ column densities to take account for unknown geometry (the H$_2$ cloud is likely not exactly spherical). In order to obtain a smooth representation of the results, we interpolated the calculated H$_2$ excitation on a denser grid in the $n_{\rm H}$--$I_{\rm UV}$ parameter space.

\begin{figure}
\centering
\includegraphics[trim={15 0 0 0},clip,width=\hsize]{J0858+1749_0_phys.pdf}
		\caption{
		Constraints on the hydrogen number density $n_{\rm H}$ and UV field strength ($I_{\rm UV}$, in Draine field units) in the H$_2$-bearing medium toward SDSS\,J085859.67+174925.32. 
		The green and violet color gradients represent the obtained probability density function (with the contour lines encompassing 68.3\% of it) using excitation of \ion{C}{i}\ fine-structure levels and H$_2$ rotational levels, respectively. The red contour shows the joint constraint with the star located at the most probable value. The red lines on the top and right axis indicate the marginalized distributions of respectively, 
		$n_{\rm H}$ and $I_{\rm UV}$ from the joint fit. The black line with downward arrows show the constraint obtained from the theoretical 
		\ion{H}{i}-H$_2$ transition in the medium.
		}
		\label{fig:J0858_phys}
\end{figure}

In Fig.~\ref{fig:J0858_phys}, we show an example of the constrained regions from \ion{C}{i}\ and H$_2$ for the DLA towards J0858+1749\ using green and purple colors, respectively. The figures for the remaining targets in our sample are shown in the Appendix. Although independently \ion{C}{i}\ and H$_2$ provide degenerate constraints in the $\rm n_{\rm H}$-$I_{\rm UV}$ parameter space, their constraints are nearly orthogonal allowing us to break the degeneracy. Therefore, a joint fit (shown by the red contour in Fig.~\ref{fig:J0858_phys}) provides significantly tighter constraints as evidenced by the marginalized joint probability densities shown on the right and top axes of Fig.~\ref{fig:J0858_phys}.
The black line in Fig.~\ref{fig:J0858_phys} indicates the constraint on the parameters using the analytical theoretical expression from \citet{Bialy2017}, based on 
the \ion{H}{i}-H$_2$ theory from \citet{Sternberg2014}, which relates the \ion{H}{i}\ column density in the envelope of H$_2$ cloud to the ratio of the UV intensity to density ($I_{\rm UV}/n$) for a given dust abundance. The latter is assumed here to scale with metallicity. 
Since neutral gas not directly related to the H$_2$-bearing cloud also contributes to the observed $N(\ion{H}{i})$, the constraint following \citet{Bialy2017} results in an upper limit on $I_{\rm UV}$ for a given $n_{\rm H}$. For all DLAs in our sample, the \ion{H}{i}-H$_2$ transition theory provides constraints in agreement with the joint measurements from \ion{C}{i} and H$_2$.

The estimates of $n_{\rm H}$ and $I_{\rm UV}$ for our sample are given in Table~\ref{table:fit_summary}.





\begin{table*}
	\centering
	\caption{Summary of column densities, chemical enrichment, and physical conditions in the sample. 
	}
	\addtolength{\tabcolsep}{-1pt}
	\label{table:fit_summary}
	\begin{tabular}{lccccccccccr}
		\hline
		Quasar & \ensuremath{z_{\rm abs}} & $\log N(\rm H_2)$ & $\log N(\rm \ion{H}{i})$ & $\log f$ & $\log N(\ion{C}{i})$ & [X/H]$^{\rm a}$ & [Fe/X] & $A_{\rm V}$\,$^{\rm b}$ & $T_{01}$ & $\log n_{\rm H}$ & ${\log I_{\rm UV}}^{\rm c}$\\
		       &       &   [cm$^{-2}$]     &  [cm$^{-2}$]    &          &  [cm$^{-2}$]  &           &        &             &  [K]     & [cm$^{-3}$] & \\
		\hline
		{\Large \strut}J\,0136+0440 & 2.779 & $18.65^{+0.06}_{-0.07}$ & $20.73\pm0.01$ &  $-1.79^{+0.06}_{-0.07}$ & $14.67^{+0.17}_{-0.15}$ & $-0.58^{+0.03}_{-0.03}$ & $-0.70^{+0.03}_{-0.04}$ & $0.16$ & $77^{+7}_{-6}$ & $2.2^{+0.1}_{-0.2}$ & $-0.2^{+0.2}_{-0.1}$ \\
        {\Large \strut}J\,0858+1749  & 2.625 & $19.72^{+0.01}_{-0.02}$ & $20.40\pm0.01$ &  $-0.53^{+0.01}_{-0.02}$ & $14.34^{+0.04}_{-0.04}$ & $-0.63^{+0.02}_{-0.02}$ & $-1.70^{+0.03}_{-0.03}$ & $0.29$ & $102^{+2}_{-2}$ & $1.8^{+0.1}_{-0.1}$ & $0.1^{+0.2}_{-0.2}$  \\
        {\Large \strut}J\,0906+0548  & 2.567 & $18.88^{+0.02}_{-0.02}$ & $20.13\pm0.01$ &  $-1.00^{+0.02}_{-0.02}$ & $13.86^{+0.06}_{-0.05}$ & $-0.18^{+0.05}_{-0.08}$ & $-1.28^{+0.08}_{-0.06}$ & $0.21$ & $116^{+26}_{-4}$ & $2.6^{+0.1}_{-0.1}$ & $0.7^{+0.1}_{-0.2}$ \\
        {\Large \strut}J\,0946+1216  & 2.607 & $19.97^{+0.01}_{-0.02}$ & $21.15\pm0.02$ &  $-0.93^{+0.02}_{-0.03}$ & $14.43^{+0.05}_{-0.05}$ & $-0.48^{+0.01}_{-0.01}$ & $-0.77^{+0.01}_{-0.01}$ & $0.25$ & $118^{+9}_{-8}$ & $2.3^{+0.1}_{-0.1}$ & $0.5^{+0.2}_{-0.2}$ \\
        {\Large \strut}J\,1146+0743  & 2.840 & $18.82^{+0.03}_{-0.02}$ & $21.54\pm0.01$ &  $-2.42^{+0.03}_{-0.02}$ & $13.91^{+0.09}_{-0.07}$ & $-0.57^{+0.02}_{-0.02}$ & $-0.67^{+0.02}_{-0.02}$ & $0.17$ & $101^{+12}_{-9}$ & ${2.6^{+0.2}_{-0.2}}^c$  &  ${0.6^{+0.2}_{-0.2}}^d$ \\
        {\Large \strut}J\,1236+0010  & 3.033 & $19.76^{+0.01}_{-0.01}$ & $20.78\pm0.01$ &  $-0.80^{+0.01}_{-0.01}$ & $14.11^{+0.08}_{-0.07}$ & $-0.58^{+0.04}_{-0.03}$ & $-1.09^{+0.04}_{-0.04}$ & $0.07$ & $118^{+7}_{-6}$ & $1.5^{+0.6}_{-0.6}$ & $0.4^{+0.4}_{-0.4}$ \\
        {\Large \strut}J\,2347$-$0051 & 2.588 & $19.44^{+0.01}_{-0.01}$ & $20.47\pm0.01$ &  $-0.80^{+0.01}_{-0.01}$ & $14.05^{+0.37}_{-0.25}$ & $-0.82^{+0.05}_{-0.03}$ & $-0.71^{+0.03}_{-0.05}$ & $0.04$ & $76^{+2}_{-2}$ & $2.0^{+0.4}_{-0.5}$ & $-0.4^{+0.3}_{-0.4}$ \\
		\hline
		{\Large \strut}J\,0843+0221 & 2.786 & $21.21^{+0.02}_{-0.02}$ & $21.82\pm0.11$ &  $-0.48^{+0.08}_{-0.08}$ & $13.98^{+0.04}_{-0.04}$ & $-1.52^{+0.08}_{-0.10}$ & $-1.01^{+0.07}_{-0.05}$ & $<0.1$ & $123^{+9}_{-8}$ & $2.82^{+0.05}_{-0.06}$ & $1.18^{+0.24}_{-0.17}$ \\
		\hline
	\end{tabular}
		\addtolength{\tabcolsep}{+1pt}
    \begin{tablenotes}
    \item (a) For all systems the metallicity given here is based on the abundance of sulphur, except for that towards J\,1146+0743 for which we used zinc because \ion{S}{ii} lines are blended.
    \item (b) The typical statistical uncertainty on $A_{\rm V}$ is 0.01~mag, but the systematic uncertainty due to intrinsic spectral slope variations is 0.1~mag (see text).
    \item (c) In units of the Draine field.
    \item (d) For the system towards J\,1146+0743, the estimated number density in the two components is similar. Here we give the value corresponding to the strongest H$_2$ component.
    \end{tablenotes}
\end{table*}


\section{Discussion}\label{sect:discussion}

\subsection{Comparison between different H$_2$-absorber samples}\label{sect:properties}

In this section, we compare the chemical properties of various samples of H$_2$-bearing DLAs (and sub-DLAs) in the distant Universe: 
this includes (1) high-redshift systems from this paper based on direct H$_2$ identification in SDSS spectra of DLAs \citep{Balashev2014}; (2) systems pre-selected using \ion{C}{I} without any prior on \ion{H}{i}\ content \citep[][]{Ledoux2015, Noterdaeme2018}; (3) systems pre-selected based on very high \ion{H}{i} column densities \citep[see][]{Noterdaeme2015b, Ranjan2019}; and (4) regular DLAs systems without any further or particular selection (mostly from the VLT/UVES database, see \citealt{Noterdaeme2008}, \citealt{Balashev2017} and references therein, but also from \citealt{Muzahid2015, Muzahid2016} for low-$z$ systems). 
For convenience, in the following we will refer to 
these samples as, respectively, 
{\sl direct-H$_2$}, {\sl \ion{C}{i}}, {\sl strong-\ion{H}{i}} and {\sl blind}.



We note that the H$_2$-bearing DLAs towards J\,0843$+$0221\footnote{This system has the second highest H$_2$ column density to date ($\log N(\mbox{H$_2$})=21.2$, \citealt{Balashev2017}), only slightly below that towards J\,1513+0352, selected for its high $N(\ion{H}{i})$ ($\log N(\mbox{H$_2$})=21.3$,  \citealt{Ranjan2018}).} and J\,1146$+$0743 belong to both the high-$N(\ion{H}{i})$ and the direct-H$_2$ selections \citep{Noterdaeme2014,Balashev2014}. 

We also compare the properties of these samples with
those observed at $z=0$ in the local group \citep{Savage1977, Rachford2009, Gillmon2006, Welty2012}. By construction (due to various selection effects) the different samples probe different regions in terms of the redshift, metallicity, \ion{H}{i}, and H$_2$ abundance, which we discuss in the following.

\subsubsection{Hydrogen content of molecular absorbers}

In Fig.~\ref{fig:H2vsHI}, we plot the dependence of H$_2$ on \ion{H}{i} column densities measured in different samples. 
It can be seen that, for a given $H_2$ column density, \ion{C}{i}-selected systems tend to probe lower 
\ion{H}{i}\ column densities than the other samples. 
This can be explained by the absence of any prior on \ion{H}{i}\ for the detection of these systems \citep{Ledoux2015}, and the corresponding sightlines may statistically intercept less unrelated \ion{H}{i}\ gas compared to DLAs. Notwithstanding, when compared to the cosmological \ion{H}{i}\ frequency distribution,
the $N(\ion{H}{i})$-distribution of the \ion{C}{i}-selected sample remains skewed towards high values. In other words, the probability to find \ion{C}{i}\ remains higher at higher \ion{H}{i}\ column densities, but the overwhelming incidence of neutral absorption systems with low $N(\ion{H}{i})$ \citep[e.g.,][]{Noterdaeme2012c}
makes them over-represented in the \ion{C}{i}\ sample.

The {\sl blind} sample should in principle be more representative of the overall DLA population. However, while the corresponding H$_2$-bearing systems indeed cover a wide range of \ion{H}{i}\ column densities, there is an excess of high-\ion{H}{i}\ systems. This reflects not only the increasing H$_2$-detection probability with increasing $N(\ion{H}{i})$ but also a bias towards high-$N(\ion{H}{i})$ in the UVES DLA/sub-DLA database \citep{Noterdaeme2008}. Since these systems were mostly observed in high spectral resolution, high-S/N data, the corresponding detection limits can indeed be as low as $N(\mbox{H$_2$})\sim 10^{14}~$cm$^{-2}$, so that systems with low molecular content \citep[e.g.,][]{Noterdaeme2007} are also present in the figure.

The strong-\ion{H}{i}\ sample, by construction, probes the high end of \ion{H}{i}\ column densities. Although the incidence rate of H$_2$ in that sample is quite high, $\sim 50$ per cent, the actual H$_2$ column densities correspond to a very wide range\footnote{Note that most ESDLAs were followed up at intermediate resolution using the X-shooter spectrograph, thus formally reported non-detections of H$_2$ actually correspond to upper limits of $\log N(\rm H_2)\gtrsim 16-17$.}, including the two highest H$_2$ column densities found in ESDLAs, towards J\,0843+0221 and J\,1513+0352 \citep{Balashev2017, Ranjan2018}.

Finally, the direct H$_2$-selection in SDSS is, by construction, mostly limited to high \mbox{H$_2$}\ column densities, typically $\log N(\mbox{H$_2$})>19$, that produce detectable lines at low spectral resolution. In addition, the search was limited to DLAs and strong sub-DLAs ($\log N(\ion{H}{i})\ge 20$) from the SDSS catalogue \citep{Noterdaeme2012c}.
{We note that six out of eight H$_2$-selected systems have 
high average molecular fractions, as seen for \ion{C}{i}-selected systems. In the case of the latter, the presence of \ion{C}{i}\ depends on the abundance of carbon (hence the metallicity) and the shielding from UV photons, both favouring high molecular fractions, as discussed in next section. In the case of the former, by selection we restrict to high \mbox{H$_2$}-column density systems when the \ion{H}{i}\ frequency distribution decreases with increasing \ion{H}{i}\ column density, also favouring high-molecular fraction systems. 
}

\begin{figure}
\centering
\includegraphics[width=\hsize]{H2vsHI.pdf}
		\caption{The dependence of $\log N(\rm H_2)$ on $\log N(\ion{H}{i})$ measured in absorption at high redshift and locally. The green pentagons, blue diamonds and yellow squares represent H$_2$-bearing DLAs from the `blind', \ion{C}{I} and strong-\ion{H}{I} selections, respectively. The sample presented in this paper and concerned with direct-H$_2$ detection is shown by red circles. Note that two systems, towards J\,0843+0221 and J\,1146+0743, belong to both 
		the H$_2$ and the \ion{H}{i}-selections and the symbols are then half red circle, half blue diamond. 
		The measurements in the local galaxies (Milky Way, LMC and SMC) are shown by purple crosses with slightly different hue. The dashed and dotted lines depict constant average molecular fractions of 0.5 and 0.1, respectively. The histograms 
	    on top show the $N(\ion{H}{i})$-distributions for the different samples using the same colour coding.
		\label{fig:H2vsHI}}
\end{figure}


\subsubsection{Dependence on metallicity}

Aside from the obvious differences in H$_2$ and \ion{H}{i}\ column densities, the different pre-selection techniques also exhibit prominent differences in the probed metallicities and redshifts. In Fig.~\ref{fig:MevsHI} we show metallicity as a function of redshift and \ion{H}{i}\ column density. We here use metallicities based on volatile species, i.e. either sulphur or zinc \citep[e.g.,][]{Vladilo2011}  depending on the line coverage and possible blends. 

From the left panel of Fig.~\ref{fig:Mevsz}, it can be seen that different samples probe slightly different redshifts. The obvious reason for this is that each selection utilizes a different redshift cutoff depending on the instrument used, the atmospheric cutoff in the blue, and the rest-frame wavelength of the transitions used for detection. Indeed, \ion{C}{i}\ lines are detectable down to lower redshifts than \ion{H}{i}, while the observation of H$_2$ lines in SDSS spectra requires even higher-redshifts than for \ion{H}{i}. Note that the sensitivity functions, $g(z)$, of the samples increases with decreasing redshift, and therefore the density of points is enhanced  close to their respective redshift cutoff.


H$_2$-bearing DLAs mainly occupy the upper half of the overall observed DLA metallicity distribution, consistent with previous works \citep[e.g.,][]{Petitjean2006}. 
While the \ion{C}{i}-selected systems probe mostly the high-end of the H$_2$-bearing DLA metallicity distribution (and consequently of the overall DLA distribution, Ledoux et al. in prep.), the strong-\ion{H}{i}-selected sample corresponds to low metallicity H$_2$-bearing systems. The direct H$_2$ selection tends to probe metallicities in between those of \ion{C}{i}\ and ESDLAs and this remains evident even after taking into account the small redshift evolution of metallicity.

The metallicity actually plays a dominant role in the detectability of \ion{C}{i}\ lines for a given \ion{H}{i}\ column density.
Not only does a high metallicity imply a high carbon abundance, but it also promotes the survival of \ion{C}{i}\ through shielding by dust and H$_2$, both of which scale with metallicity. Therefore the \ion{C}{i}-selection predominantly picks out high metallicity systems, in particular at low $N(\ion{H}{i})$ \citep[][Ledoux et al. (in prep)]{Noterdaeme2018,Zou2018}.
Interestingly, for a given metallicity, the \mbox{H$_2$}\ column densities in \ion{C}{i}-selected systems increases with increasing $N(\ion{H}{i})$. Since theoretical models of \ion{H}{i}-H$_2$ conversion indicate that the \ion{H}{i}\ envelope of a single cloud has almost constant column density whatever the $N(\mbox{H$_2$})$ (i.e. more gas added to a cloud is converted to H$_2$), this could indicate the interception of multiple \ion{H}{i}-H$_2$ layers. This is consistent with a higher probability to find cold gas close to a given cold gas cloud, either along the same line of sight (appearing as multiple components, see Ledoux et al. in prep), or transversely \citep{Krogager2018}.

The chemical conditions for the presence of H$_2$ at detectable levels are less restrictive than those for \ion{C}{i}, allowing the direct H$_2$ and strong \ion{H}{i}-selections to probe lower metallicities. Furthermore, in the high \ion{H}{i}\ column density regime, the presence of \mbox{H$_2$}\ in a given system appears to be little dependent on the metallicity. Indeed, the DLA towards J\,0843+0221 has a very high \ion{H}{i}\ and H$_2$ content despite its low metallicity, and the metallicity distribution of H$_2$-bearing ESDLAs is not significantly different from that of non H$_2$-bearing ESDLAs (see Fig.~\ref{fig:MevsHI}). 
The local physical conditions (density and UV field) seem to play a more important role for these systems. Indeed, these systems are expected to probe small galacto-centric distances \citep{Noterdaeme2014,Ranjan2019}, where the pressure could be increased \citep[in theoretical models of warm and cold gas distribution in galaxies][]{Wolfire2003}. 
We also note that high-$N(\ion{H}{i})$ DLAs with high metallicities could be under-represented among them due to the associated dust reddening. This could exclude the background quasars from optically-selected samples \citep[see, e.g.,][]{Geier2019, Krogager2019} and lead to infeasible amounts of observing time required for spectroscopic follow-up.


The intermediate metallicities seen for the H$_2$-selected sample suggest that the intrinsic metallicity distribution of H$_2$-bearing systems is shifted compared to that of the overall DLA population.
However, the paucity of low-metallicity systems in the direct H$_2$ sample {\sl in comparison to H$_2$-bearing ESDLAs}, if real, can be due also to differences in physical conditions. The presence of H$_2$ depends not only on the metallicity and column density but also on the number density ($n$) and intensity of the UV field ($I_{\rm UV}$). To understand this further, we compare the observed $N(\ion{H}{i})$ and metallicities with the theoretical expectation from the work by \citet{Bialy2016}.
We caution that the theoretical predictions are only valid for a single cloud whereas the actual situation is more likely a concatenation of H$_2$ clouds along the line of sight embedded in a diffuse and warm \ion{H}{i}\ phase that also contributes to the measured \ion{H}{i}\ column density. In that sense, all measured values of $N(\ion{H}{i})$ are upper limits to the \ion{H}{i}\ column density associated with the envelope of the H$_2$ cloud.
The different theoretical lines should therefore be regarded as the minimal conditions for \ion{H}{i}/H$_2$ transition at a given metallicity and \ion{H}{i}\ column density. That is, the calculated \ion{H}{i}\ column densities are only applicable for systems where a transition from \ion{H}{i}\ to \mbox{H$_2$}\ has occurred. For the typical ranges of the physical conditions in the cold ISM (mainly, the metallicity, number density and UV flux) 
this corresponds to a \mbox{H$_2$}\ column density larger than $\log N (\rm H_2) \sim 18-20$.
Therefore we highlight the points correspond to different ranges of \mbox{H$_2$}\ column densities in Fig.~\ref{fig:H2vsHI} using color saturation. The blank points represent the lowest \mbox{H$_2$}\ column densities ($\log N (\rm H_2) < 18$) for which the model lines are most likely not applicable. Almost all saturated H$_2$-bearing DLAs, regardless of their selection, are located in between the $I_{\rm UV}/n_{\rm H}\sim 10^{-3}$ and $I_{\rm UV}/n_{\rm H}\sim 1$ lines. Systems with low metallicity and \ion{H}{i}\ column density would require very low $I_{\rm UV}/n_{\rm H}$ ratios (below $\lesssim10^{-3}$) in order to efficiently form H$_2$. High density systems with low total column density must have very small physical dimensions leading to a low absorption cross-section. This is consistent with the absence of H$_2$ among such absorption-selected samples.
For example, for a given UV field typical of our Galaxy (i.e. $I_{\rm UV} = 1$) this correspond to $n\gtrsim10^3$\,cm$^{-3}$. Therefore at column density at the DLA threshold, and metallicity typical of DLAs ([X/H]$\sim-1.3$), a cloud should be compressed less than $<0.1$\,pc in order to H$_2$ to form. For a temperature of $T\sim 100$~K, this also means a thermal pressure $P\sim10^5$\,K\,cm$^{-3}$, which is not very common, especially at low N(\ion{H}{i}). 

To summarize, the three preselection methods presented above complement each other and enable us to study the physical conditions and the H$_2$/\ion{H}{I} transition in a wide range of environments.

\begin{figure*}
\centering
\setlength{\tabcolsep}{2pt}
    \begin{tabular}{cc}
    \includegraphics[trim={0 0 0 0},clip,width=0.5\textwidth]{Mevsz.pdf}
    & \includegraphics[trim={0 0 0 0},clip,width=0.5\textwidth]{MevsHI.pdf} \\
    \end{tabular}
\setlength{\tabcolsep}{6pt}
\caption{The metallicity as a function of the redshift (left panel) and $\log N(\ion{H}{i})$ (right panel) measured in DLAs. The symbols are the same as in Fig.~\ref{fig:H2vsHI} with colour saturation representing the H$_2$ column density. The mean metallicity in DLAs as a function of redshift is shown on the left panel by dash \citep{Rafelski2012} and dotted \citep[$N(\ion{H}{i})$ weighted and depletion corrected, see][]{DeCia2018} lines. The different lines in the right panel correspond to the $N(\ion{H}{i})$ in the envelope of a single H$_2$ cloud \citep{Bialy2017} for the different ratios of the  UV field to the number density, $I_{\rm UV} / n$, where $I_{\rm UV}$ is expressed in units of the Draine field and the $n$ in units of cm$^{-3}$.} 
\label{fig:MevsHI}
\label{fig:Mevsz}
\end{figure*}


\subsubsection{Fe depletion and dust}

In previous discussions, we did not differentiate between the metallicity and dust content, assuming that the dust content directly scales with the metallicity. However, the presence of H$_2$ is directly linked to the dust content, since the formation of molecules occurs on the surface of dust grains that also participate to the shielding of UV photons.
Unfortunately, the large systematic uncertainties of the $A_V$ measurements, compared to the actual (low) values, prevent us from quantifying the dust content directly from the measured extinction. Instead, we can use the depletion of refractory elements (like iron) compared to volatile ones (zinc or sulphur).
All samples follow a similar trend of increasing Fe depletion with increasing metallicity, as illustrated in Fig.~\ref{fig:DeptoMe}, reflecting the scaling of dust with metallicity \citep[see][]{Ledoux2003}.
Since the H$_2$ formation rate is roughly proportional to the dust abundance, we might statistically expect a enhanced depletion for H$_2$-bearing DLAs at a given metallicity. 
However, we do not clearly observe such a behaviour: only two out of eight DLAs in our sample have a significant stronger depletion (i.e. lower [Fe/X]) than seen from the overall DLA population at their given metallicity. This can be 
due to a fact that the quantities used for comparison are averaged throughout the absorption complexes \footnote{See, e.g., the H$_2$-bearing DLA towards J\,0000+0048, where the detailed analysis of high-resolution data indicates that the H$_2$-bearing component is significantly more depleted than the average over all DLA components.}. 
More importantly, the dependence on average depletion shown in Fig.~\ref{fig:DeptoMe} doesn't take into account the column densities. Observationally, it is already known that the total column of the dust probably plays a major role in the presence of H$_2$ at detectable levels \citep[see][]{Noterdaeme2008}. The physical reason is that larger dust column imply larger shielding from UV photons, either directly from the dust grains, or indirectly provided by H$_2$ self-shielding (since the more dust, the higher H$_2$ formation rate and the larger H$_2$ column density). 
In Fig.~\ref{fig:H2vsFedust} we show the dependence of the H$_2$ column density on the column density of iron locked into dust, which provides a way to quantify the column density of dust \citep{Vladilo2006}. Overall, the different samples are not distinguishable in this Figure and confirm the neat increase in the H$_2$ column densities at $\log N(\rm Fe_{\rm dust})\gtrsim15$. This suggests that the physical conditions (ratio of UV to number density) should not be drastically different for different H$_2$ samples, since the critical $N(\rm Fe_{\rm dust})$ is mostly defined by this ratio \citep[see, e.g.,][]{Bialy2017}.
For given physical conditions, an increasing $N(\mbox{H$_2$})$ with increasing $N(Fe_{\rm dust})$ could then be expected at higher dust column, regardless of whether this is due to a concatenation of several transition layers seen along a line of sight, or single clouds with larger column densities. Such a correlation may possibly be seen in the \ion{C}{i}-selected sample which possibly share more homogeneous conditions.


\begin{figure}
\centering
\includegraphics[trim={0.0cm 0.0cm 0.0cm 0.0cm},clip,width=\hsize]{FetoZnvsMe.pdf}
		\caption{The depletion versus metallicity in DLA at high redshifts. Symbols and color encoding are the same as in Fig.~\ref{fig:Mevsz}.\label{fig:DeptoMe}}
\end{figure}

\begin{figure}
\centering
\includegraphics[trim={0.0cm 0.0cm 0.0cm 0.0cm},clip,width=\hsize]{H2vsFe_col.pdf}
		\caption{H$_2$ column density versus column density of Iron in the dust. Symbols and color encoding are the same as in Fig.~\ref{fig:Mevsz}. The grey crosses show the upper limits for DLAs taken from \citealt{Noterdaeme2008, Ranjan2019}.
		\label{fig:H2vsFedust}}
\end{figure}


\subsection{Physical conditions}

We find values of the temperature in the range 75-120~K, with a median of $T\sim 100$~K\footnote{We note that this value corresponds to the temperature of the H$_2$-bearing gas and not the overall DLA system, which is more likely characterized by a warm neutral phase with $T \gtrsim 10^3$\,K.}. 
These values are slightly higher than what is measured towards nearby stars \citep[$T\sim 80$~K, see, e.g.,][]{Savage1977}, but lower than previously seen in the overall population of high-$z$ DLAs \citep[$T\sim 150$~K;][]{Srianand2005}. However,
$T_{01}$ temperatures in H$_2$-bearing DLAs are known to decrease with increasing H$_2$ column density \citep{Muzahid2015, Balashev2017}, a trend also seen in the local ISM.
The systems in our sample have typically higher \mbox{H$_2$}\ column densities than those from \citet{Srianand2005}, but lower than towards nearby stars, that can explain the difference. However, from our sample alone, we do not see any trend in T$_{01}$ vs. $\log N (\rm H_2)$, though our measurements are still within the measured dispersion of T$_{01}$. 
Metallicity could also play a role to explain the differences in the temperatures since it plays an important role in 
the thermal balance \citep[e.g.,][]{Wolfire2003}. 
{Indeed, the cooling is dominated by line emission (in particular [C\,{\sc ii}] at $158\mu\,{\rm m}$) while the heating is dominated by photoelectric effect on dust grains. Both these processes are directly linked to the gas metallicity.}  
Observationally, a low H$_2$ temperature of $T\sim 50$~K was measured in the super-Solar metallicity DLA towards J\,0000+0048 \citep{Noterdaeme2017} while higher temperatures where measured at low metallicities, $\rm [X/H]\lesssim -1$ towards  J\,0843+0221 and J\,1513+0352 \citep{Balashev2017, Ranjan2018}. 
{However, since the photo-electric heating involves UV photons,
the intensity of the UV field directly plays an important role.} Indeed, one can see in Fig.~\ref{fig:T01}, where the estimated UV flux is coded by the color saturation, that our sample (for which most systems span a narrow range of the metallicity) supports such a dependence. 
We can further speculate that, knowing the gas metallicity, the relative $J=1/J=0$ population of H$_2$ could give a useful diagnostic of the UV to density ratio, without the need of detailed modelling.

\begin{figure}
\centering
\includegraphics[trim={0.0cm 0.0cm 0.0cm 0.0cm},clip,width=\hsize]{H2TKin.pdf}
		\caption{The excitation temperature T$_{01}$ versus column density of molecular hydrogen. Symbols are as per Fig.~\ref{fig:H2vsHI}, with colour-filling for the sample presented here (red-circles) representing the inferred intensity of the UV field (see enclosed colour bar). 
	
		\label{fig:T01}}
\end{figure}

\begin{figure}
\centering
\includegraphics[trim={0.0cm 0.0cm 0.0cm 0.0cm},clip,width=\hsize]{H2_Pn.pdf}
		\caption{The dependence of the  thermal pressure ${\rm P}/k=nT$ on $\log N(\rm H_{\rm tot})$ measured in DLAs (circles) and local galaxies (square and triangles). Green and red circles correspond to the known H$_2$/\ion{C}{i}\ bearing DLAs at high redshift and our sample, respectively.
		\label{fig:PHI}}
\end{figure}
We found a typical hydrogen number density of $n \gtrsim 100$ cm$^{-3}$ and typical UV field on average twice higher than the Draine field. The obtained number densities and kinetic temperatures agree well with the canonical values for the cold ISM in the two-phase model from \citet{Wolfire2003}. This shows that H$_2$ is a good tracer of the CNM in DLAs. This is particularly useful since 21-cm absorption studies remain limited at high redshift \citep{Srianand2012}, due to the paucity of strong compact radio sources observable with current facilities \citep[but see][for future prospects]{GuptaMALS}.

Using our estimates of the number density and temperature, we compare the thermal pressures (${\rm P}/ k = n T$) to those seen in other H$_2$-bearing DLAs as well as in the local group in Fig.~\ref{fig:PHI}. 
Our sample probes a wide range of hydrogen column densities, from values similar to that of the Magellanic clouds \citep{Welty2016} down to about an order of magnitude lower than in the Milky Way \citep{Jenkins2011}. The pressure in our sample is about $\sim10^4$\,K\,cm$^{-3}$ with 0.5 dex dispersion around this value, similar to what is seen locally. Our data further strengthen the suggestion from \citet{Balashev2017} that the thermal pressure increases with the hydrogen column density. This is likely the consequence of both the pressure and the column density to increase with decreasing galacto-centric radius.


\section{Conclusion}
\label{sect:conclusion}

We presented follow-up observations with the X-shooter spectrograph on the Very Large Telescope of seven strong H$_2$-bearing DLAs at high redshifts. These DLAs were selected from the SDSS DLA catalogue \citep[see][]{Noterdaeme2012c} directly upon the presence of H$_2$ Lyman bands seen in the SDSS spectra of their background quasars, see \citet{Balashev2014}.
Our X-shooter observations confirm the detection of H$_2$ in all systems and hence the high efficiency of the direct H$_2$ detection method even at low spectral resolution. We measured the neutral, molecular, and metal column densities of the systems and found typical H$_2$ column densities in the range $\log N(\rm H_2) = 18.5-20$, when the \ion{H}{i}\ column density span a range $\log N(\ion{H}{i}) \sim 20.1-21.5$. The metallicities, obtained from undepleted species, are in the range 
$\log Z = -0.8$ to $-0.2$. 
Finally, in all systems, we detected absorption lines from \ion{C}{i}\ in various fine-structure energy levels associated to the H$_2$ components.
We used the rotational excitation of H$_2$, the fine-structure excitation of \ion{C}{i}, and theoretical relations for the \ion{H}{i}-H$_2$ transition to estimate the physical conditions in the H$_2$-bearing gas. 
We found a typical kinetic temperature of $\sim100$\,K and hydrogen densities $n_{\rm H}\sim100$\,cm$^{-3}$, within clouds embedded in a UV field about twice higher than the Draine field.
This confirms that H$_2$ is a good tracer of the cold neutral medium\footnote{However, \citet{Srianand2012} noted there is no one-to-one correspondence between the population of cold gas traced by H$_2$ and that traced by \ion{H}{i} 21-cm absorption.}.

Comparing our sample with other samples of H$_2$-bearing DLAs and sub-DLAs selected differently, we discuss how the selection methods probe different regimes of column density and chemical enrichment and how
the presence of H$_2$ depends on the physical conditions. 
We show that all these methods complement each other and by the SDSS limitation tend to select high H$_2$ column-density systems that can be efficiently used to study the \ion{H}{i}-H$_2$ transition.

\section*{Acknowledgements}
We thank the referee for useful comments and suggestions.
This paper was supported by RSF grant 18-12-00301.
PN and JKK acknowledge support from the \textit{French Agence Nationale de la Recherche}, under grant ANR-17-CE31-0011-01 (Projet "HIH2", PI Noterdaeme) and from the CNRS-RFBR Russian-French grant PRC No\,1845 (PI Noterdaeme). We are also grateful to the Institut d'Astrophysique de Paris (SB and VK) and to the Ioffe Institute in Saint Petersburg (PN and JKK) for hospitality when part of this work was achieved.




\bibliographystyle{mnras}
