\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}


\section{Introduction}

The problem of estimating entropy of high dimensional binary configurations or signals is ubiquitous in many disciplines.
In physics, we very often have at our disposal a set of configurations of some physical system generated by a Monte Carlo simulation
at a given temperature $T_0$. This data is very much geared towards computing   expectation values of various operators or their correlation functions, however obtaining
the entropy or free energy of the system is far from trivial. Indeed, to the best of our knowledge, there is no known way to compute the entropy directly from
these configurations~\cite{COMMENT} even for a system of a quite moderate size (e.g. for a $20 \times 20$ lattice).
The goal of the present paper is to propose machine learning based methods
which would allow to make such a computation.

Due to the lack of a direct method to obtain entropy, various indirect approaches have been proposed.

One standard way requires to perform separate Monte Carlo simulations at a series of temperatures from close to $T=0$
up to $T_0$, evaluate the heat capacity $C(T)$ from the variance of
the energy $\sigma^2_E(T)$
\begin{equation}
C(T) = \f{\partial\! \cor{E}}{\partial T} = \f{\sigma^2_E(T)}{T^2}
\end{equation}
and obtain the entropy by a numerical integration over~$T$:
\begin{equation}
S(T_0) = \int_0^{T_0} C(T) \f{dT}{T}
\end{equation}
Clearly this process is time consuming and requires performing multiple auxiliary Monte Carlo simulations even if one is interested only in
the entropy of the system at a specific temperature $T_0$.

Another approach is to use Wang-Landau sampling~\cite{WangLandau} for the given hamiltonian to obtain the density of states $g(E)$, and subsequently compute the entropy.
This way does not use the original Monte Carlo configurations at all and requires to perform a quite separate (and conceptually different) Monte Carlo entropic sampling
to evaluate the entropy.
Needless to say, Wang-Landau sampling requires for the system to be described by an explicit hamiltonian,
which is of course clearly given in the physics context, but requires separate modeling in other contexts (e.g. for neuron spike recordings).
In the latter case, in order to use any of the above approaches, one would have to assume a specific form of an effective hamiltonian describing the system, which is in fact an additional prior input. Determining its free coefficients is an involved problem by itself, even before
trying to evaluate the entropy.

Indeed, within the field of neuroscience, one is often interested in the dynamics of spiking time series of populations of neurons and in particular in assessing their information content through an estimation of entropy. For small populations,
one can clearly evaluate the entropy directly (see section~\ref{s.entropy} below), however as the number of available neurons grows, the space of allowed
configurations increases exponentially.
A common approach, as indicated above, is to model the system through an effective
Ising model~\cite{IsingNeurons} (equivalently a Hopfield network~\cite{Hopfield})
\begin{equation}
\label{e.maxentising}
H = \sum_{i<j} J_{ij} s_i s_j + \sum_i h_i s_i
\end{equation}
which provides a maximal entropy description subject to the constraint of reproducing all 1-point and 2-point correlation functions.
Once one solves this ``inverse Ising'' problem of fitting the parameters $J_{ij}$ and $h_i$ \cite{InverseIsing}, one can then use a physics based approach
(like mean field or one of the more precise methods like temperature integration or the Wang-Landau sampling mentioned above) to evaluate the entropy. As emphasized previously, this is however not a
trivial problem even if the exact coefficients $J_{ij}$ and $h_i$ are known.

The aim of this work is to propose a viable method which computes the entropy directly from the original configurations (e.g. the
Monte Carlo configurations at a given temperature, or neuron spike trains, or any other relatively high dimensional binary signal).
In particular, we would like the method to be applicable when the dimensionality of the system is such that there are basically no
repeated configurations. This is in fact the generic situation in physics simulations, where e.g. the Ising model on even a small $20 \times 20$
lattice has of the order of $10^{120}$ possible configurations\footnote{For the $20000$ Monte Carlo configurations at each temperature for the $20 \times 20$ Ising model used in this paper, we checked that indeed all configurations are distinct for temperatures $T\geq 2.7$.}.

\section{Entropy and its estimation}
\label{s.entropy}

Suppose that an $N$-dimensional binary signal $\{x_i\}_{i=1..N}$ comes from a probability distribution
\begin{equation}
p(x_1,x_2,\ldots, x_N)
\end{equation}
Then the (Shannon) entropy\footnote{In this paper we will always use base 2 logarithms, thus the entropy is measured in \emph{bits}.} is given by
\begin{equation}
\label{e.shannon}
S = - \!\!\!\!\!\! \sum_{\{x_i\}_{i=1..N}}\!\!\!\! p(x_1,x_2,\ldots, x_N) \log_2 p(x_1,x_2,\ldots, x_N)
\end{equation}
where the sum is over the $2^N$ allowed configurations.
In the physics context, the Shannon entropy is of course equivalent to the thermodynamic entropy, when
the probability distribution is given by a Boltmann distribution associated to some hamiltonian:
\begin{equation}
\label{e.boltzmann}
p_{Boltzmann}(x_1,x_2,\ldots, x_N) = \f{1}{Z} e^{-\f{1}{T} H(x_1,x_2,\ldots, x_N) }
\end{equation}
where $T$ is the temperature.

The difficulty in evaluating the entropy directly from the definition~(\ref{e.shannon}) has two sources. Firstly, even if the explicit functional form (\ref{e.boltzmann})
is known and the hamiltonian $H(x_1,x_2,\ldots, x_N)$ is easy to evaluate as is usually the case,
the computation of the normalizing constant $Z$ (i.e. the partition function in physics) and the summation
over configurations in (\ref{e.shannon}) is intractable by brute force. Hence the need for the theoretically much more involved approaches
mentioned in the introduction.

Secondly, especially outside the physics context, we do not have at our disposal the functional form (\ref{e.boltzmann}) with a simple \emph{known} hamiltonian, and the problem lies in reliably estimating
the probability distribution $p(x_1,x_2,\ldots, x_N)$ given a set of $n$ samples from the distribution.
For relatively small $N$ and a sufficient number of samples, one can directly count the number of times $k$
a given configuration $\{x_i\}_{i=1..N}$ occurred and set
\begin{equation}
\label{e.plugin}
p(x_1,x_2,\ldots, x_N) = \f{k}{n}
\end{equation}
Substituting this back into (\ref{e.shannon}) gives the plug-in (or Maximal Likelihood) estimator.
Taking into account the bias for small sample sizes/occurrence counts leads to a variety
of improved estimators like the Miller-Madow \cite{MillerMadow}, Chao-Shen \cite{ChaoShen}, Grassberger \cite{Grassberger}, James-Stein \cite{Hausser}
and others. Further refinements involve Bayesian estimators based on Dirichlet priors (see a summary and references in \cite{Hausser}),
and the most refined estimators obtained by fitting a mixture of Dirichlet distributions:
the Nemenman-Shafee-Bialek (NSB) estimator \cite{NSB} and the CDM estimator of \cite{CDM}.

Another quite different approach is based on an expansion of the Shannon entropy in terms of mutual information.
This comes from rewriting exactly the entropy in terms of entropies of individual subsystems
and then adding corrections due to first pairwise mutual information and subsequently their higher order
generalizations:
\begin{equation}
\label{e.mie}
S = \sum_{i=1}^N S(i) -\sum_{i<j} I_2(i,j) +\sum_{i<j<k} I_3(i,j,k)+\ldots
\end{equation}
where $S(i)$ is the entropy of the subsystem made up of neuron/spin $x_i$,
$I_2(i,j)$ is the mutual information between the subsystems $x_i$ and $x_j$:
\begin{equation}
I_2(i,j) = S(i) + S(j) - S(i,j)
\end{equation}
etc. Unfortunately, truncation of the expansion (\ref{e.mie}) at say $2^{nd}$ order
may yield huge errors -- even leading to a negative entropy estimate for the Ising model in the intermediate temperature regime.
In principle this should be improved by higher order terms, however evaluating
them explicitly becomes combinatorially intractable for larger $N$.

Another interesting proposed approach (with somewhat wider goals) \cite{MaxentropyNature}, aims at maximizing the possible entropy
satisfying a set of information theoretic constraints. Unfortunately this method
also seems to be limited to moderate dimensionalities $N$.

\section{The proposed method}
\label{s.method}

As indicated in the introduction, we would like our method of estimating entropy to work also in the case of dimensionalities $N$, where
the given samples involve basically only distinct configurations, thus in formula~(\ref{e.plugin}), $k=0$ or $k=1$
and we cannot use any of the subsequent refinements.

Before we present the details, let us comment on the general philosophy of the proposed method.

The methods of entropy estimation based on occurrence counts like~(\ref{e.plugin}) and its variations treat each distinct configuration as a structureless ``atomic'' object. All details of its constituents are thrown away. Our method, in contrast, strives to analyze the internal structure of the configurations and thus can work in the regime where all occurring configurations are distinct. In order to extract the possibly quite complex internal relationships between the constituents, we will use machine learning methods, which, as we will show, appear here in a very natural way.

Our starting point is the exact rewriting of a multidimensional probability distribution in terms of
a product of conditional probabilities
\begin{equation}
\label{e.condprobfact}
p(x_1,x_2,\ldots, x_N) = p(x_1) p(x_2|x_1) p(x_3|x_1,x_2) \cdot \ldots
\end{equation}
Let us first evaluate Shannon's entropy using the above decomposition.
We clearly get a sum of $N$ terms:
\begin{equation}
\label{e.decomposition}
S = S_1 + \Delta S_2 + \Delta S_3 + \ldots \Delta S_N
\end{equation}
The first term is just the entropy of the first neuron/spin:
\begin{equation}
S_1 = -p_1 \log_2 p_1 - (1-p_1) \log_2 (1-p_1)
\end{equation}
where $p_1 \equiv p(x_1=1)$.

The second term is more interesting. We need to evaluate
\begin{equation}
\Delta S_2 = -\!\!\!\sum_{\{x_1,x_2\}} \!\! p(x_1, x_2) \log_2 p(x_2|x_1)
\end{equation}
Since $x_2$ is a binary variable, the conditional probability is completely specified
by the simple function
\begin{equation}
p_2(x_1) \equiv p(x_2=1 | x_1)
\end{equation}
Then an estimate of $\Delta S_2$ based on the dataset $\{x_1^{(k)}, x_2^{(k)}\}_{k=1 \ldots n}$ is
\begin{equation}
\label{e.deltastwo}
-\f{1}{n} \sum_{k=1}^n x_2^{(k)} \log_2 p_2(x_1^{(k)}) + (1-x_2^{(k)}) \log_2 (1- p_2(x_1^{(k)}))
\end{equation}
We observe, that this is exactly the cross-entropy loss in a standard supervised classification problem
where we treat the value of the neuron/spin $x_2$ as a class label $y$ (which can be either 0 or 1),
and we try to predict its probability based on the value of neuron/spin $x_1$.
Similarly for $\Delta S_3$ we get
\begin{eqnarray}
-\f{1}{n} \sum_{k=1}^n &&\hspace{-0.5cm} x_3^{(k)} \log_2 p_3(x_1^{(k)}, x_2^{(k)}) + \nonumber\\
 &&\hspace{-0.8cm} + \left(1-x_3^{(k)}\right) \log_2 \left(1- p_3(x_1^{(k)}, x_2^{(k)})\right)
\label{e.deltasthree}
\end{eqnarray}
where now we need to predict the probability of $x_3=1$ in terms of
the values of $x_1$ and $x_2$.

Thus the decomposition (\ref{e.decomposition}) corresponds exactly to a sum of cross-entropy losses
of a sequence of iterative supervised classification problems where we predict the probability of
$x_j=1$ given the values of the previous\footnote{According to the given fixed ordering $x_1, x_2, \ldots, x_N$.}
spins $x_1, x_2, \ldots, x_{j-1}$.

Note, however, that these auxiliary classification problems have a somewhat different flavour from the typical ones encountered in machine learning.
There, we expect usually the target class to be completely determined by the features (predictors), like whether a given image represents a cat or a dog. So the ideal classifier should attain zero cross-entropy loss.
In our case, we expect generically for the best possible classifier to achieve only some positive cross-entropy loss.
It is exactly this nonzero cross-entropy which summarizes the contribution of the $j$-th spin/neuron to the total entropy.

Indeed, suppose that spin $x_j$ is completely independent of $x_1, x_2, \ldots, x_{j-1}$.
Then the conditional probability
\begin{equation}
p_j(x_1, x_2, \ldots, x_{j-1}) \equiv p(x_j=1|x_1, x_2, \ldots, x_{j-1})
\end{equation}
will just be a constant equal to the marginal $p(x_j=1)$ and $\Delta S_j$ will
be equal to the entropy of the neuron/spin $x_j$ by itself.

If, on the other hand, $x_j$ would be completely determined by $x_1, x_2, \ldots, x_{j-1}$,
we would expect the ideal classifier to give vanishing cross-entropy and $\Delta S_j=0$.

Let us emphasize that the above procedure of computing the entropy depends explicitly on the predefined ordering of the particular
neurons/spins $x_1,x_2,\ldots, x_N$, as the auxiliary classification problems for two orderings are completely different.
On the other hand, the final outcome, the entropy of the system,
is clearly independent of any ordering. This property may be a nontrivial cross-check of the
quality of the employed machine learning classifier, hence it would be beneficial to evaluate
the entropy for a couple of permutations of the neurons/spins.

Let us note that in the limit of infinite data, the entropy estimate~(\ref{e.decomposition}) using some concrete classification algorithm will approximate the true entropy from above. This can be easily seen e.g. for the difference $\Delta S_2^{true} - \Delta S_2^{classifier}$ equal in the infinite data limit
\begin{equation}
\!\!\!\sum_{\{x_1,x_2\}} \!\! p(x_1, x_2) \log_2 \f{p_{classifier}(x_2|x_1)}{p(x_2|x_1)}
\end{equation}
Using Jensen's inequality this is smaller or equal than
\begin{equation}
\log_2 \!\!\!\sum_{\{x_1,x_2\}} \!\! p(x_1, x_2) \f{p_{classifier}(x_2|x_1)}{p(x_2|x_1)}
\end{equation}
which in turn vanishes due to
\begin{equation}
\log_2 \!\!\!\sum_{\{x_1,x_2\}} \!\! p(x_1) p_{classifier}(x_2|x_1) = \log_2 1 = 0
\end{equation}
A similar argument holds for other terms.
Hence if various classifiers yield different answers, we should chose the one which predicts the
lowest entropy. Of course, for a finite number of samples these theoretical conclusions are not guaranteed to hold.

Finally, as emphasized in more detail at the end of the following section, we should take care how we evaluate the predictions of the fitted classifiers and never evaluate them on the data used for training.

\section{Machine learning considerations}
\label{s.ml}

In section \ref{s.method}, we reinterpreted the decomposition (\ref{e.decomposition}) for
the formula for the entropy (\ref{e.shannon}) in terms of cross-entropy losses of
a sequence of supervised classification problems. This is very general and comes from
the fact that a discriminative supervised classification problem is just
by definition a certain model of the conditional probability distribution
\begin{equation}
\label{e.condprob}
p(y|x_1,x_2,\ldots,x_j)
\end{equation}
where $y$ is a binary (target/class) variable and $x_1,x_2,\ldots,x_j$ are some predictors/features.

The reason, why this interpretation is fruitful is that we would like to estimate (\ref{e.condprob})
in a data driven way based on the given set of binary signals (representing e.g. Monte-Carlo spin configurations,
neuron spiking time series etc.) and the field of Machine Learning provides in fact a very wide variety
of algorithms which are specifically designed to model \emph{very complex} conditional probability distributions\footnote{A classical example would be the probability that an image shows a cat conditioned on the set of pixels of the image.},
and thus can be directly used now to estimate entropy.
The most notable examples are logistic regression, $k$-nearest neighbours, deep neural networks,  random forests, gradient boosted trees,  and others.

Note that the relevance of the machine learning algorithm for computing entropy is indicated not by classification accuracy, but rather by a good estimate of probability. Thus e.g. Naive Bayes is not expected to be a good choice, Support Vector Machines (SVM) by default do not provide probabilities at all  and deep neural networks may require some caution (cf.~\cite{Khan}).

The freedom in the choice of classifier means that we have at our disposal a whole range of methods for computing entropy. This variety may be very useful when we have some prior knowledge about the given binary signals. In the case when the number of samples is not so large, Bayesian models/inference may be very effective. One can also use any of the techniques
developed within machine learning for feature selection as an ingredient for the
procedure of entropy estimation.

In this paper we will not attempt to investigate any of such refinements, but rather, as a proof of concept of the method, use a standard \emph{off-the-shelf} classifier and use it with basically default settings

A very versatile nonlinear classifier which will be used in the numerical experiments
in the present paper is the so-called \emph{gradient boosted tree} classifier.
Specifically, we chose {\tt xgboost} -- its very efficient\footnote{For computational speed, in this paper we always use the settings {\tt tree\_method='hist'} and {\tt n\_jobs=-1}.} variant and implementation \cite{xgboost}.
In some cases, we will also use logistic regression for comparison.

Before we close this section, we need to address a remaining important issue
concerning the evaluation of the predictions $p_j(x_1,x_2,\ldots,x_{j-1})$ used in formulas like
(\ref{e.deltastwo}) or (\ref{e.deltasthree}). In Machine Learning one has to avoid ``overfitting'',
where the classifier will fit random noise in the data. This is especially dangerous for the more complex
nonlinear classifiers, which are incidentally of the most interest in the present context. Moreover, predicting directly on the training data used to fit the classifier would strongly skew the predicted probabilities and thus would tend to underestimate the entropy if we have a finite amount of data
(as is always the case).

A standard way to mitigate this problem is to partition the dataset into $k$ (at least $k=2$ but usually $k=5$) parts (folds),
train the classifier on the sum of $k-1$ parts, but compute predicted probabilities only on the remaining  unseen $k$-th fold.
Then repeat the process $k-1$ times, holding out another fold, until one gets predictions for all datapoints.
Then the cross-entropy loss like (\ref{e.deltastwo}) or (\ref{e.deltasthree}) entering (\ref{e.decomposition}) should be computed only using these held-out predictions\footnote{In order not to be biased by a specific choice of the partition of the dataset into $k$ folds, we pick a different shuffled partition for each auxiliary classification problem.}.

\section{The importance of nonlinearity}

In this section, we will consider synthetic datasets where some binary features are
functionally dependent on others.

Let $X_1$ be given by 50 independent random binary variables, each with probability $0.5$, and $X_2$ be defined similarly. We will now form four variants of a third set
of binary variables $X_3$ which are given respectively by $\text{NOT}\ X_1$, $X_1\ \text{OR}\ X_2$, $X_1\ \text{AND}\ X_2$ and $X_1\ \text{XOR}\ X_2$.
The final datasets will be obtained by sampling the
concatenation of $X_1$, $X_2$ and one of the variants of $X_3$.
Hence in total
we have 150 binary variables, whose entropy is clearly equal to $100.0$ bits irrespective of the
choice of $X_3$.
We will now perform a random reordering of the variables, take $10000$ samples and apply
our proposed method for computing the entropy.

It is illuminating to look at the answers obtained by using two standard classifiers:
logistic regression\footnote{We use the implementation from {\tt scikit-learn} with default parameters.} and gradient boosted trees ({\tt xgboost}). The former is essentially a baseline linear classifier, while the latter is, as mentioned earlier, a fully nonlinear one. The results are shown in Table~\ref{tab.dependence}.

\begin{table}[h]
    \centering
    \begin{tabular}{|c|c|c|c|c|}
    \hline
         &  NOT & OR & AND & XOR \\
    \hline
    lr    & 100.68 & 101.23 & 101.26 & 151.05 \\
xgb(100)  & 100.74 & 100.74 & 100.71 & 122.27 \\
xgb(200)  & $\cdot$& $\cdot$& $\cdot$& 112.60 \\
xgb(400)  & $\cdot$& $\cdot$& $\cdot$& 105.79 \\
\hline
    \end{tabular}
    \caption{The results for estimating the entropy for the four datasets, with logistic regression,
    {\tt xgboost} with the default number of trees (100), and with 200 and 400 trees respectively. The exact answer for all datasets is $100.0$ bits.}
    \label{tab.dependence}
\end{table}{}

Unsurprisingly, the interdependence structure of the XOR dataset cannot be captured by a linear model, and logistic regression indeed very strongly overestimates the entropy. The nonlinear {\tt xgboost} classifier is definitely better, but the structure of the dataset requires increasing the complexity of the model (by increasing the number of trees from the default 100 to 400 and possibly more\footnote{Note, however, that increasing the complexity of the model is not always beneficial. See e.g. the case of the 2D Ising model in the following section.}) in order to correctly estimate the entropy.

This example serves to illustrate that the flexibility of our proposed method in the choice of the specific machine learning classifier can yield interesting insights into the data. By comparing different classifiers we can e.g. assess the complexity of nonlinear interdependence of our configurations or signals. Alternatively, we can sometimes restrict ourselves to a simpler and faster model if the more complex model does not yield a significantly lower entropy estimate.

\section{Ising model entropy and free energy from Monte Carlo configurations}

As a nontrivial cross check of the proposed method we will evaluate
the entropy of the 2D Ising model on a $20 \times 20$ periodic lattice
directly from sets of 20000 configurations obtained using Monte Carlo sampling
for temperatures ranging from $T=1.0$ to $T=4.0$, thus spanning
both the low and high temperature phase of the infinite volume theory and the critical phase transition in between.

The 2D Ising model is especially interesting as a testing ground
for the evaluation of entropy for a number of reasons. Firstly,
there are analytical exact formulas for the entropy and free energy
of the model put on any $L \times L$ lattice (see the formulas in Appendix~A). Secondly, it is not a trivial system but exhibits
a phase structure with two distinct phases separated by
a second order phase transition. Thirdly, it has a huge number of possible configurations $2^{L^2}$, thus making it a challenging testing ground for entropy estimation. Fourthly, it is a quintessential example for Monte Carlo simulations, thus the
direct computation of the entropy and free energy just from the Monte Carlo configurations would be an interesting proof of concept of the proposed method. Finally, generalizations of the Ising model with arbitrary pairwise couplings and inhomogeneous magnetic field (\ref{e.maxentising})
are commonly employed as maximal entropy models of neuron spiking in neuroscience.

\begin{table}[t]
    \centering
    \begin{tabular}{|c|c|c|r|}
    \hline
      {\tt max\_depth}   &  {\tt n\_estimators} & $S(T=T_c)$ & time \\
      \hline
       3  & 25 &   0.44125 & 7.7\\
       3  & 50 &   0.42815 & 12.1\\
       3  & 75 &   0.42868 & 16.5\\
       3  & 100 &  0.42960 & 21.0\\
       3  & 200 &  0.43330 & 38.7\\
       \hline
       5  & 100 &  0.43885 & 37.2\\
       5  & 200 &  0.45196 & 71.2\\
       \hline
       \multicolumn{2}{|c|}{logistic regression} & 0.43518 & 7.9 \\
       \hline
      \multicolumn{2}{|c|}{\bf exact} & {\bf 0.42468} & \\
      \hline
    \end{tabular}
    \caption{The estimated entropy per spin computed from 20000 Monte Carlo samples at $T=T_c$ for various choices of hyperparameters of {\tt xgboost} as well as logistic regression for completeness. The running time is in minutes on a 6-core desktop. The exact answer is given for a $20\times 20$ lattice (see Appendix~A).}
    \label{tab.hyper}
\end{table}{}

Since the 2D Ising model has a clear spatial structure, and we would like to test our method in a more general context, we will randomly permute the $400$ spins to set the ordering for (\ref{e.condprobfact}).

Before we present the predictions for the whole range of temperatures it is instructive to look at the dependence of the outcome on various factors.

Firstly, two key hyperparameters of the gradient boosted trees algorithm are the number of trees ({\tt n\_estimators}) and their maximal depth
({\tt max\_depth}). In Table~\ref{tab.hyper}, we show the dependence of the entropy estimate on these choices at the
critical point $T=T_c \sim 2.2691853$.

We see that for 25 trees, the model is too simple, while for
increasing depth and number of trees it starts overfitting. Thus, for all the remaining experiments we choose the default {\tt max\_depth=3} and just decrease the number of trees from the default $100$ to~$50$. This is also clearly better than a logistic regression baseline.

It is also interesting to analyze the dependence of the entropy estimate on the number of Monte Carlo samples. The results for
$n=10000$, $n=20000$ and $n=40000$ are shown in Table~\ref{tab.n}.
The results get better as we gain access to more data.
This motivated us to compute a linear extrapolation in $1/n$, although there are
no known theoretical grounds for a specific functional dependence.
This is, in fact, an interesting problem for further study.

\begin{table}[h]
    \centering
    \begin{tabular}{|c|c|}
    \hline
        $n$ & $S(T=T_c)$ \\
        \hline
      10000   &  0.42989\\
      20000   &  0.42815\\
      40000   &  0.42748\\
      \hline
      extrapolated & 0.42661\\
      \hline
      {\bf exact} & {\bf 0.42468}\\
      \hline
    \end{tabular}
    \caption{The estimated entropy per spin at $T=T_c$ as a function of the number of Monte Carlo samples, together with a linear extrapolation in $1/n$.}
    \label{tab.n}
\end{table}

As emphasized in section \ref{s.ml}, it is important to estimate the probabilities by predicting on an unseen test set, especially for more complex classifiers, so we always use $k$-fold cross validation.
Increasing $k$, increases the amount of training data for the machine learning algorithm on each fold, thus leading to better results. The drawback, however, is that the running time increases as one has to run each classifier $k$ times to obtain predictions for all samples. In Table~\ref{tab.cv} we show this dependence on~$k$. In all other simulations we take $k=5$.

\begin{table}[h]
    \centering
    \begin{tabular}{|c|c|c|}
    \hline
        $k$ & $S(T=T_c)$ & time\\
        \hline
        2 & 0.42935  & 2.7\\
        3 & 0.42857  & 5.2\\
        4 & 0.42829  & 7.4\\
        5 & 0.42815  & 9.1\\
      \hline
      {\bf exact} & {\bf 0.42468} &\\
      \hline
    \end{tabular}
    \caption{The estimated entropy per spin computed from 20000 Monte Carlo samples at $T=T_c$ for various number of cross validation folds together with the running time in minutes.}
    \label{tab.cv}
\end{table}

As emphasized in section~\ref{s.method}, the procedure for estimating the entropy depends explicitly on the ordering of the variables (spins/neurons) as the auxiliary classification problems (i.e. which spin/neuron to predict based on which subset of spins/neurons) are quite different for different orderings. Yet the sum of the cross entropy losses for all the classification problems should not depend on the ordering. This is a nontrivial consistency check of the method and may serve to assess whether the chosen machine learning classifier is adequate for the given dataset and to estimate the lower bound on the error\footnote{Of course, the difference between the entropy estimate and the true entropy may be much higher than that, as the classifier may \emph{systematically} fail to identify all predictive regularities in the dataset, c.f. logistic regression for the XOR dataset in the previous section.}.

In Table~\ref{tab.order}, we show the results obtained for a natural row-wise ordering, five different random orderings as well as two orderings based on correlation. For the {\it max corr.} one, we iteratively pick a spin which is most correlated\footnote{As measured by the sum of absolute values of correlation coefficients $\sum_{j\in previous} |C_{ij}|$.} with the ones chosen earlier, while for the {\it min corr.} ordering we pick the least correlated spin. We see that the results for various orderings are consistent between themselves.

\begin{table}[t]
    \centering
    \begin{tabular}{|c|c|}
    \hline
    ordering & $S(T=T_c)$ \\
        \hline
    row-wise            &  0.42821\\
    max corr.           &  0.42818\\
    min corr.           &  0.42820\\
    \hline
    random (seed 0)     &  0.42815\\
    random (seed 137)   &  0.42816\\
    random (seed 555)   &  0.42810\\
    random (seed 1621)  &  0.42812\\
    random (seed 4567)  &  0.42817\\
      \hline
      {\bf exact} & {\bf 0.42468} \\
      \hline
    \end{tabular}
    \caption{Entropy estimates computed for various orderings of the spins.}
    \label{tab.order}
\end{table}

With the above exploratory analysis done, we compute the entropy for the range of temperatures $T=1.0$ to $T=4.0$ from 20000 Monte Carlo configurations (at each temperature) and compare with the exact answer for the Ising model on the $20 \times 20$ periodic lattice (see Appendix~A for explicit formulas).
We use 5-fold cross validation, {\tt max\_depth=3} and {\tt n\_estimators=50}.
We use a random ordering with seed 0.
The results are shown in Fig.~\ref{fig.entropy}.

\begin{figure}[t]
\hspace{-0.2cm}\includegraphics[width=8.25cm]{ising_entropy.png}
\caption{Entropy per spin computed using {\tt xgboost} classifier compared with the exact answer. The magenta data point at $T=1.0$ uses data augmentation (see text).}
\label{fig.entropy}
\end{figure}

Once we have the entropy, we can immediately compute the free energy from
\begin{equation}
F= \cor{E} -T S
\end{equation}
where the expectation value of the energy is trivial to compute from the Monte Carlo configurations. The results are shown in Fig.~\ref{fig.free}.

\begin{figure}[t]
\hspace{-0.2cm}\includegraphics[width=8.25cm]{ising_free_energy.png}
\caption{Free energy per spin computed using {\tt xgboost} classifier compared with the exact answer. The magenta data point  at $T=1.0$ uses data augmentation (see text).}
\label{fig.free}
\end{figure}

Note that there is a significant deviation for the lowest temperature $T=1.0$. There, the system is in the ordered phase where all spins are predominantly oriented in the same direction. This means that the classification problems are very strongly imbalanced. In fact, it is quite surprising that the {\tt xgboost} classifier works so well for slightly higher temperatures in the ordered phase.

One way to alleviate the problem is to use more data.
Instead of generating further Monte Carlo samples,
one can instead use the very standard machine learning procedure of data augmentation, namely constructing additional samples based on the original data (e.g. cropping and flipping images in image classification tasks). In the context of the Ising model, we can rotate the spin configurations by multiples of $90^{\circ}$ or
perform arbitrary (periodic) translations. Note that from the point of view of the auxiliary classification problems, this yields completely novel additional data. For the magenta points in the figures, we added multiples of $90^{\circ}$ rotations of the original 20000 samples, giving all together 80000 data points. We see that the entropy estimate is now much better.

Of course, data augmentation is only possible if the relevant system has some known symmetries. This is a very common situation for various physical systems but does not occur e.g. for recordings of spiking neurons. In the latter case data augmentation is impossible and one would have to acquire more real data.

\section{Discussion}

In this paper we translated the problem of computing the entropy of a set of binary configurations or signals into performing a sequence of supervised classification tasks, whose sum of cross-entropy losses provides an estimate of the entropy.
We showed that the method is powerful enough to reproduce quite well the entropy and free energy of the 2D Ising model directly from Monte Carlo configurations.

This framework is very general and allows to use the whole machine learning toolbox with a wide range of diverse machine learning classifier algorithms for computing entropy.
This flexibility may be used, on the one hand, to chose a classifier most suited to the general structure of the data of interest. And, on the other hand, a comparison of the entropy estimates due to different classifiers may lead to interesting insights into the data, like being an indication of the inherent complexity/nonlinearity of the signals.

We hope that these methods will be very useful in physics, where they may be used to estimate the entropy and free energy directly from a set of Monte Carlo configurations.
Within neuroscience, they may be an aid in estimating the entropy of larger populations of spiking neurons as well as other brain signals. Comparing the entropy of a maxentropy model like~(\ref{e.maxentising}) with the entropy obtained from the original signal may be an aid in estimating the applicability of the specific model to the observed data.

There are numerous directions for further study. It should be possible to extend this approach to other information theoretic quantities.
It would be interesting as well to develop extensions for differential entropy (of continuous probability distributions). Another direction would be to investigate the optimal machine learning algorithms and/or feature selection procedures in various specific contexts like the case of very high dimensionality.

\bigskip

\noindent{}{\bf Acknowledgments.}
I would like to thank Przemysław Witaszczyk for the ongoing collaboration which motivated this investigation and Piotr Białas for interesting discussions on entropy.
This work was done in preparation for the Foundation for Polish Science (FNP) project \emph{Bio-inspired Artificial Neural Networks} POIR.04.04.00-00-14DE/18-00.

\bigskip

\noindent{}{\bf Code and data availability.} All code for reproducing the numerical experiments performed in this paper is available at \href{https://github.com/rmldj/ml-entropy}{\tt github.com/rmldj/ml-entropy}. The 2D Ising Monte Carlo configurations were generated by \href{https://github.com/rmldj/ising}{\tt github.com/rmldj/ising}, a minimally modified fork of \href{https://github.com/zeehio/ising}{\tt github.com/zeehio/ising}. The generated Monte Carlo configurations are also available at \href{https://doi.org/10.5281/zenodo.3457123}{\tt doi:10.5281/zenodo.3457123}.