Neutral carbon in diffuse interstellar medium: abundance matching with H2 for DLAs at high redshifts

We present the study of CI/H$_2$ relative abundance in the diffuse cold neutral medium. Using the chemical and thermal balance model we calculated the dependence of CI/H$_2$ on the main parameters of the medium: hydrogen number density, metallicity, strength of the UV field, and cosmic ray ionization rate (CRIR). We show that observed relative CI and H$_2$ column densities in damped Lyman alpha systems (DLAs) at high redshifts can be reproduced within our model assuming the typically expected conditions in the diffuse cold neutral medium (CNM). Using the additional observed information the on metallicity, HI column density, and excitation of CI fine-structure levels, as well as temperature we estimated for a wide range metallicities in the CNM at high redshifts that CRIRs to be in the range from $\sim10^{-16}$ to $\rm few \times 10^{-15}\rm s^{-1}$, hydrogen number densities to be in range $\sim10 - 10^3$cm$^{-3}$, and UV field in range from $10^{-2}$ to $\rm few \times 10^2$ of Mathis field. We argue, that since the observed quantities used in this work are quite homogeneous and much less affected by the radiative transfer effects (in comparison with e.g. dissociation of HD and UV pumping of H$_2$ rotational levels) our estimates are quite robust against the assumption of the exact geometrical model of the cloud and local sources of the UV field.


INTRODUCTION
Carbon is one of the most important elements in the neutral phase of the local interstellar medium (ISM).Indeed, at metallicities close to solar, its singly ionized form, C ii, provides the principal source of the electrons (that determine the chemistry in the medium) and is one of the dominant cooling processes in diffuse regions, that determines the formation of the cold phase of neutral medium (CNM).In turn, in translucent and dense molecular regions carbon is mainly in C i and/or CO, which also dominate the cooling of the medium.Therefore, it is essential to study carbon in CNM of ISM both observationally and theoretically.
From observations, the most robust way to constrain the abundance and ionization state of carbon in CNM is an absorption line spectroscopy.However, while C i is typically a sub-dominant form of carbon, regarding measured column densities, its abundance can be much more confidently constrained than C ii.Additionally, C i in absorption always observed populating three fine-structure levels, which provide an important diagnostics of the physical conditions in the medium (see e.g.Jenkins & Tripp 2001;Silva & Viegas 2002;Srianand et al. 2005;Noterdaeme et al. 2007;Balashev et al. 2011;Jorgenson et al. 2010;Klimenko & Balashev 2020;Kosenko et al. 2023).C ii abundance (Jenkins & Tripp 2011;Roman-Duval et al. 2019) and excitation of its fine-structure level (Balashev et al. 2022) can also be used for constraints on the physical state of the gas.However, the usage of C ii is much more limited than C i due to the ★ E-mail: s.balashev@gmail.comavailability of less number of transitions and saturation of the absorption lines, although the incidence rate of C ii is much higher than C i.In turn, while rotational excitation of CO can provide independent diagnostics of the physical conditions (Klimenko et al. in prep), observations of CO are much more limited than C ii and C i, since the cross-section of the associated dense medium is quite small.
The C i chemistry in diffuse ISM was comprehensively considered theoretically by e.g.Liszt 2011;Wolfire et al. 2008.However, these studies considered the observations obtained in our galaxy, i.e. examined only the local ISM.In turn, to date, there is a quite large number of observations of C i obtained in high-redshift DLAs, i.e. remote galaxies.These observations are mostly associated with the low-metallicity gas ( ≲ 0.3), where one can expect the changes in the thermal (see e.g.Bialy & Sternberg 2019;Balashev et al. 2022) and ionization balance (see e.g.Balashev & Kosenko 2020;Balashev et al. 2021), which is important for the chemistry of C i.An important constituent of these changes are cosmic rays, which start to be a dominant heating source of cold ISM at  ≲ 0.3 (Bialy & Sternberg 2019) and govern its ionization state.Therefore C i is expected to be sensitive to the cosmic rays ionization rate (CRIR) and hence can provide a way to measure it.
Previously, some indirect methods were also used to estimate the cosmic ray ionization rate based on the observations of the abundance of a several molecules that are sensitive to the degree of ionization of the medium.Most of the molecules used, such as H + 3 , OH + , H 2 O + , HD (see e.g.Hollenbach et al. 2012;Indriolo & McCall 2012;Indriolo et al. 2015;Kosenko et al. 2021), represent diffuse and translucent molecular gas.These methods show that there is a large spread in the measured CRIR: from 10 −17 to ∼ few×10 −15 s −1 .Such a scatter can be caused both by the natural inhomogeneity of the CRIR in the interstellar medium, associated with the locality of sources and the effects of cosmic ray propagation, and by systematic effects in the methods of estimation used.Therefore independent estimations of CRIR are highly valuable.
In this work, we present an independent formalism to calculate C i abundance in the diffuse ISM (see Sect. 2) and its application on the sample of C i measurements obtained in the high-redshift H 2bearing DLAs.In particular, in comparison with previous studies (e.g.Wolfire et al. 2008;Liszt 2011) we focus on the low metallicity gas, how it affects the C i abundance, and the importance of the cosmic rays in C i production (see Sect. 4.1).We developed the method (presented in Sect.4.2) for estimating the main physical conditions in the medium, based on measured C i and molecular hydrogen, H 2 , column densities, the population of the C i fine structure levels, and kinetic temperature (derived through the ortho-para ratio of H 2 ).A new important feature of this work that discriminates it from our previous studies (e.g.Klimenko & Balashev 2020;Kosenko et al. 2023) is the consideration of C i abundance and constraints on CRIR that C i provides coupled with explicit calculations of the thermal state of the gas.We applied this method to a sample of known C i/H 2bearing DLAs at high redshifts and compared it with some previous estimates (see Sect. 4.3).

MODEL
In this Section we describe a formalism that was used to calculate relative C i/H 2 abundances, thermal state, and the excitation of C i fine-structure levels.In the following, during the modelling, we paid attention to the following quantities: H i and H 2 column densities, metallicity, the ortho-para ratio of H 2 (providing the measurements of the kinetic temperature), and the column densities on three C i fine-structure levels.The latter provide the measurement of the total column density of C i and hence the relative abundance of C i to H 2 .These quantities from the observations will be used to derive constraints on the main physical parameters of the model.

Model global parameters
We considered a homogeneous (specified by the total hydrogen number density  tot H ) medium with metallicity  (denoted by default in decimal units) that is exposed by the UV field,  (normalized to Mathis field, Mathis et al. 1983) and cosmic rays with ionization rate  -the primary ionization rate per hydrogen atom (measured in s −1 ).These are the main global parameters of interest.
For each considered species X in the medium we denote the number and column density as  X and  X , respectively 1 .These quantities evidently are related through the integration along the line of sight,  X = ∫  X .The column densities, are also observed quantities, and hence are the subject of the model calculation.For convenience, in the following text, we will use H 2 column density,  H 2 , as a natural unit of measure of the cloud thickness, instead of physical length.The column density of the species X can be calculated as 1 In the following, we assume that the number and column densities are measured in units cm −3 and cm −2 , respectively.

Metal and dust abundance
The thermal balance depends on the metal and dust abundances.The dust abundance was specified by the dust-to-gas ratio, DTG.The observations suggest that it is scaled with metallicity as a power law or broken power law.However, the measured dispersion of DTG is quite large, therefore for simplicity, we assume that where DTG ⊙ is the dust-to-gas ratio at solar abundance and  is typically constrained between 1 and 3.While this dependence is simplistic, as was mentioned by Balashev et al. 2022, it matches with the observations from Rémy-Ruyer et al. (2014).Finally, when we used the model calculations we took into account the dispersion of DTG-metallicity dependence.
The abundance of the elements in the gas phase (which is also necessary to properly calculate the cooling rate) can be determined using the metallicity and DTG ratio as where X/H ⊙ is a solar photospheric abundance of some species denoted by  (we use values taken from the Asplund et al. 2009), and  X is a depletion of species at solar metallicity.For the latter, we used the following values 0.53, 0.41, and 0.0 for carbon, oxygen, and argon, respectively.We did not apply depletion for other species, since they participate in neither thermal balance nor basic chemical reactions, considered in this work.
We also considered Polycyclic Aromatic Hydrocarbons (PAH) separately, while, we assumed that PAH abundance is linearly scaled with DTG ratio, i.e. the distribution of the grain size does not change with metallicity.In principle it doesn't have to be true, since the production of the large grains can be non-linearly dependent on the small grain abundance, as well as the large grains destructed/formed into/from smaller grains, which may change their relative fractions.While other authors Wolfire et al. (2003Wolfire et al. ( , 2008)); Bialy & Sternberg (2019) used artificially added  PAH coefficient to rescale the collisional coefficients of the PAH measured in the laboratory, we did not consider this factor, i.e. used  PAH = 1.

H 2 /HI abundance
The abundances of H 2 and H i determine the ionization state and chemistry of the cloud and are hence important for the thermal state and chemical abundances of the species.In principle, to calculate the relative H 2 and H i abundances, and the H i/H 2 transition one needs to calculate full radiative transfer in H 2 Lyman and Werner bands.However, this is a very time-consuming task and here to be able to efficiently calculate the models, we used the analytical description of H i/H 2 transition proposed recently by Sternberg et al. (2014); Bialy & Sternberg (2016)  2 .Following this formalism, the number density of H 2 as a function of H 2 column density,  H 2 can be written as where  H is the atomic hydrogen column density, which can be derived as a function of the H 2 column density (see Bialy & Sternberg 2016),  H 2 is a self-shielding function,  g ≈ 1.9 × 10 −21 [DTG] (in cm 2 ) is the dust Lyman-Werner photon absorption cross-section per hydrogen atom, and  is the ratio of the rates of free space H 2 photodissociation and H 2 formation on the dust grains.In turn, atomic hydrogen density is simply defined as  H =  tot H − 2 H 2 .

Chemical abundances
For the purpose of this work, we considered a sample of species, that are abundant in diffuse clouds and affect the thermal and ionization balance (hence have an impact on the chemistry and excitation of C i).
The major species of the hydrogen network are H, H 2 , H + , H + 2 , H + 3 .For the oxygen molecular network we considered O, O + , OH + , H 2 O + , H 3 O + , H 2 O, and OH.We also calculated He, He + , HeH + , Ar, Ar + , and ArH + that impact the oxygen molecular chemistry and electron (denoted by ) abundance.In the carbon network, we considered only C and C + and did not append the network of C-bearing molecular species.This is reasonable since the model will be applied only for diffuse clouds, before the onset of CO, where the abundances of any other C-bearing molecules are very low and affect neither carbon chemistry nor thermal balance.We also explicitly considered PAH species at three ionization states: PAH 0 , PAH + , and PAH − , since they affect ion recombination, such as H + and C + .
For H and H 2 we used analytical description of the abundances (see Sect. 2.3).We also used the fixed He, O, and Ar abundance (as described in Sect.2.2, except He for which we used value 0.085 of hydrogen abundance, independent from metallicity), since they are dominant ionization form in cold diffuse medium.For other species, we calculated the chemical abundances by solving the matrix equations of the chemical network.We used the same reaction rates as provided by Balashev et al. 2021 mostly compiled from UMIST database McElroy et al. 2013.For reaction rates of PAHs we used compilation by Wolfire et al. 2008.We additionally considered the dependence of photoreactions on the dilution of the UV field, with coefficients taken from the recent compilation by (Heays et al. 2017).To solve chemical equations we additionally assume that the total PAH abundance is

Excitation of the fine-structure levels
The excitation of fine-structure levels is important for the calculation of the cooling functions and the column densities of the fine-structure levels, that was used to compare with observed ones.We calculated the excitation of fine-structure levels of C ii, C i and O i using standard matrix calculations (e.g.Silva & Viegas 2002), assuming excitation by collisions, radiative excitation by CMB and UV pumping, and de-excitation by spontaneous transitions.We considered collisions with the most abundant species in diffuse ISM: H, H 2 (separating ortho and para), He (if available), and electrons.The collisional coefficients for C + were taken from (Barinovs et al. 2005;Flower & Launay 1977;Keenan et al. 1986); for C from (Abrahamsson et al. 2007;Schroder et al. 1991;Staemmler & Flower 1991;Johnson et al. 1987); for O from (Abrahamsson et al. 2007;Jaquet et al. 1992;Monteiro & Flower 1987;Bell et al. 1998).We considered the excitation by Cosmic Microwave Background with the temperature set as 2.725(1+) K at the particular redshift of each studied absorber.This excitation can be important only for C i (which has the smallest separation of the fine-structure levels) and redshift ( ∼ 2 − 4) while provides subdominant contribution to the population of C i levels, that is typically observed.

Thermal balance
The temperature of the medium was obtained using the thermal balance between the cooling and heating processes in ISM.We considered the cooling by Ly (see e.g.Spitzer 1978), recombination (we used coefficient from Bialy & Sternberg 2019), and C ii, C i, and O i fine-structure emission.The latter was calculated using excitation of the fine-structure levels (see Sect. 2.5) and assuming the optically thin regime.The main heating sources considered in this work are the photoelectric heating on the dust grains and cosmic rays, which depend on the physical parameters and electron abundance.We follow the standard formalism to obtain the heating rates, most recently described in Bialy & Sternberg 2019.However, we additionally take into account the dependence of the photoelectric heating on the extinction of the UV field, using scaling of the rate by a dilution of the UV photons, calculated using the typical extinction curve (Fitzpatrick & Massa 2007) with We did not consider heating by X-ray background, since it is likely important only in the case of the presence of a strong X-ray emitting source.We also did not take into account the turbulence dissipation heating.While as was mentioned in (Balashev et al. 2022) it can be important at low metallicities and even can be comparable with cosmic ray heating, the intermittent nature of the turbulence heating, may indicate that the heating can be very local (see discussion in Klessen & Glover 2016), which complicates the calculations.We also did not consider heating and cooling by H 2 , since it is only important either at low metallicity gas or at dense compact regions exposed by large UV field (see e.g.Bialy & Sternberg 2019), while we considered mostly diffuse ISM.

Solution
To obtain the column densities of the species for each set of the physical parameters, defined in Sect.2.1 and 2.2, we used the following routine.We divided the cloud of some length (specified by  H 2 ) by a number of layers, with 0.1 dex step starting from log  H 2 = 14.
For each layer, we first determine  H and  H 2 using the analytical formalism of H i/H 2 transition (see Sect. 2.3).With these values and also fixed values of He, O, and Ar abundance (see Sect. 2.4) we calculate the cloud structure layer by layer (starting from the first) using iterative routine.At each iteration we solved a chemistry network to get abundances,  X and  X , of other species in consideration (see Sect. 2.4) together with the thermal balance to get temperature (see Sect. 2.6) and excitation of fine-structure levels (see Sect. 2.5).The iterations were stopped when next step does not change the derived values by more than 0.5 percent for each quantity, which set by a trade-off between accuracy and computational time.In most layers, several iterations were enough starting from the values, obtained at the previous layer.The typical calculation time of one model was about a few seconds per one core on the standard laptop, which we found quite enough for efficient sampling of the models presented in Sect.4.2.

DATA SAMPLE
To compare the model presented in the previous section with observations, we compiled observed C i and H 2 column densities from all known H 2 -bearing DLAs detected at high redshift ( ≳ 2) towards quasar sightlines.We note that there are a few additional detected systems at high-redshifts presented in (Noterdaeme et al. 2018) that have only total C i and H 2 column densities, i.e. there is no information on C i fine-structure excitation,  01 temperature, and metallicity.We also did not consider so-called proximate H 2 -bearing DLAs (e.g.Noterdaeme et al. 2019Noterdaeme et al. , 2023)), which are typically defined as systems that have redshift close to the quasar redshift.i.e.Δ ≲ 3000 km/s.These systems can be significantly impacted by the AGN emission (e.g.Balashev et al. 2020) or even be produced in the AGN outflowing gas (e.g.Noterdaeme et al. 2021), therefore they may be not representative for the diffuse medium, and require a separate study.The compiled sample is presented in Table 1.In the case where several H 2 /C i-bearing components were detected along the line of sight, we provide explicitly the individual H 2 /C i components within DLAs, since the described modelling allows us to consider them independently.

RESULTS
In this section, we present the results of the modelling and its application to the observational data.
An example of the calculated temperature, abundances, and excitation of C i fine-structure levels, together with values of the chemical reaction and cooling/heating rates are shown in Fig. 1 for the model with  tot H = 100 cm −3 ,  = 1,  = 10 −15 s −1 , and  = 0.3.One can see that profiles of the species behave well as expected -where the main feature is H i/H 2 transition (see panel (b) in Fig. 1), which affected the thermal balance (panel (e)), temperature (panel (a)), ionization (panels (b)) and chemical structure (panel (g), (h)) and at the less amount PAH abundances (panel (d)).It is important to notice, that the temperature, C i abundance, and relative population of C i fine-structure levels are more or less constant and do not show strong variations throughout the cloud.Therefore these quantities are likely to have little bias regarding the geometrical model of the cloud, since they are little coupled with the radiative transfer.This will be important for the discussion of the results in the following subsections.It is interesting to note that there is a slight increase of the kinetic temperature in the H 2 dominated region in comparison with atomic hydrogen dominated region, due to lower particle density in H 2 dominated region and lower collisional excitation rate of excited C ii fine-structure level by H 2 , than by atomic hydrogen.

Variation of C i/H 2 abundance
We investigate the dependence of C i/H 2 relative abundances on each physical parameter in the model.Since our calculations were based on one-sided radiation field models (inherited from formalism by Sternberg et al. 2014), to compare model results with observational data, here and in the following text we calculated the column densities of each species at the half of specified H 2 column density and then multiplied by two, i.e.  obs X = 2 X ( H 2 /2).This approximately emulates the plain parallel cloud of the thickness  H 2 exposed by UV field on both sides 3 .
In Figure 2 we plot the comparison of the observed and calculated dependence of C i and H 2 column densities.We plot all the systems from the parent sample, presented in Sect.3, while we note that for several systems, where the column densities on high fine-structure were not measured, the total C i column densities can be slightly higher, but not more than ≈ 30 per cent.It is evident to expect that C i abundance may depend on the metallicity (at least carbon abundance almost linearly scaled with it), and since the observed metallicities at high redshifts systems span about two orders from log  ∼ −1.5 to 0.5, we divided the sample in the four metallicity bins ≈ 0.4 sizes (in dex) and calculated theoretical  obs C i ( H 2 ) profiles for the mean values in each bin.As we mentioned above, C i abundance strongly depends on three global parameters of the model: the number density,  tot H , CRIR, , and UV field strength, .To show how variations of these parameters affect C i/H 2 relative abundance we constructed the base model with  tot H = 100 cm −3 ,  = 10 −15 s −1 and  = 1 and varied independently each of the parameters within two dexes, that correspond to the typical measured variations, see Klimenko & Balashev (2020); Kosenko et al. (2021).The resulting profiles of C i/H 2 abundances are shown in Figure 2. One can see that except a few points at intermediate metallicities ⟨log ⟩ ∼ −0.8 and −0.4 observed data can be reproduced within chosen ranges of the physical parameters.
One can see that at different metallicities C i/H 2 abundance is sensitive to different parameters.At high metallicities,  ≳ 0.5, it is mostly sensitive to the variation of the UV field and number density.This is simply because the UV field directly scales the C i abundance by photoionization process.However, at low metallicities,  ≲ 0.2 (actually corresponding to the most DLAs at high redshifts), C i/H 2 abundance also becomes quite sensitive to CRIR, since at these metallicities, C ii recombination rate depends on the electron abundance, which in turn depends on the hydrogen ionization fraction, which is sensitive to CRIR.Once metallicity approaches solar value, the electron densities start to be determined solely by the carbon abundance and hence C i/H 2 becomes a little sensitive on CRIR (the bottom left panel in Figure 2).Interestingly, at low metallicities, C i/H 2 abundance becomes almost insensitive to the UV flux (the top central panel in Figure 2), at the typically observed values in the range between 0.1 and 10 of the Mathis field.One can additionally note, that at intermediate H 2 column densities, i.e. log  H 2 ≲ 20, C i/H 2 abundance have opposite dependence on the number density at low and high metallicities, i.e. decreases and increases with the number density increase, respectively.

Constraints on the physical parameters
As we showed in previous section abundance C i/H 2 , C i fine structure excitation and ortho-para ratio of H 2 depend on physical conditions ( tot H , , ) one can constrain the latter using observed quantities together with metallicity measurement.Importantly, as we discussed in Sect.4, aforementioned quantities (that we will use for comparison with observations) do not drastically vary along the cloud, and hence these estimates are most probably less biased by unknown geometry of the absorbing region than most other species do, e.g.HD or O-bearing molecules.Since the dependency on the physical parameters likely has a relatively complex shape in the parameter space, and our model is quite a little time-consuming to constrain these parameters we followed the Bayesian approach to sample the posterior probabilities function of the  tot H , ,  using affine invariant Monte-Carlo Markov Chain sampler (Goodman & Weare 2010).The likelihood function was constructed by comparison of model values with observations of the following quantities: (i) the kinetic temperature, measured by column densities of  = 0, 1 levels of H 2 ; (ii) the column densities of fine structure levels of C i (this includes the relative excitation and total C i column density) (iii) the total H i column density associated with DLA.The latter was used as an upper limit only since we do not know how much of the observed H i attributed to the C i-bearing component, but we found it useful since it can shrink the posterior distribution functions.We considered H 2 column density as a model parameter since it determines the thickness of the cloud in the model.However, we take into account the uncertainty of the measurements of  (H 2 ) using a Gaussian prior, which is well within the Bayessian approach.Additionally, the metallicity was also considered as a model parameter with Gaussian prior corresponding to the measured interval estimate for each system.Since the C i formation is quite sensitive to the dust-to-gas ratio, we considered the DTG ratio as an independent parameter during the modelling, while adding a prior to follow dependence with metallicity with  = −1.5 (see Eq. 1) and additional 0.5 dex dispersion, to take into account the observed scatter in DTG-metallicity relation.For the other model parameters we used flat priors in log space, that conditionally emulate a wide distribution, in the absence of other physically motivated priors.In the following, to report the point and interval estimates on the model parameters we used maximum a posterior probability estimate and the highest posterior density 0.683 credible interval, respectively.
An example of the result of the modelling of the H 2 /C i absorption  (e) the rate of the heating (dashed lines) and cooling (solid lines); (f, g, and h) the rate of reactions of the formation (solid lines) and destruction (dashed lines) for C, PAH − , and , respectively.system (at z=2.626 towards Q 0812+3208) using described method is shown in Fig. 3 using the corner plot python package (Foreman-Mackey 2016).As we see in this particular system, we were able to reproduce the observational quantities and obtain reasonable constraints on the model parameters.In this particular system, we found that UV field, CRIR, and number density are relatively enhanced in comparison with typically found values for diffuse clouds.For the latter, it is consistent with the observed high excitation of C i fine-structure levels.

High redshift sample
We applied the procedure described in the previous Section to the C i/H 2 -bearing DLAs at high redshifts,  ≳ 2, detected so far.The parent sample was described in Sect. 3 and shown in Table 1.We selected a subsample of systems, where C i column density on each fine-structure level was measured, to make the analysis homogeneous since excitation of C i fine-structure levels provides important constraints on the model parameters.The described modelling allows us to consider the individual H 2 /C i bearing components of DLAs, unless they have appropriate column density measurements of H 2 , C i fine-structure levels, and  01 temperature.However, since we can not determine the fraction of the H i associated with each component, we used an assumption that the metallicities of the components within DLA are the same and equal to the DLA average.
The results of the determinations of the physical parameters are summarized in Table 2.We found that estimates on CRIR are located within the range from 10 −16 to few × 10 −15 s −1 , while the UV estimates show a wider range of values.The number densities are found to be in the range of  tot H = 10 − 1000 cm −3 (note that this is the total density of hydrogen nuclei, so the actual number density can be ≈ two times lower than  tot H ), which is consistent with the values, expected for the cold neutral medium.
In Fig. 4 we show the relative dependence between obtained values of CRIR and UV fluxes.It seems that there is no evident correlation between these parameters unless one considers some fraction of the data.We note, that since these parameters compete on their contribution to the thermal balance (through cosmic ray heating or photoelectric heating) in some absorption systems we found that they are degenerated.However, usage of the C i abundance can relax this degeneration, since it provides independent constraints on CRIR.Additionally, we did not confirm the quadratic dependence of the CRIR on UV flux, which was previously reported by Kosenko et al. 2021 for the diffuse ISM, based on the modelling of HD/H 2 relative abundance coupled with C i fine-structure excitation.In comparison with the aforementioned values, we obtained relatively higher values CRIR on average, while the significant fraction of the points around log  ∼ −15 and log  ∼ 1 are consistent with each other.The higher values of the CRIR that we found from C i modelling likely arose from the matching of C i/H 2 abundance, since the C i abundance depends on the recombination rate, which is tightly related to the electron abundance (directly at low metallicities, or indirectly through the dependence on the PAH − abundance).However, we note that the constraints based on C i are more robust, regarding the radiative transfer issues than HD/H 2 -based measurements.Indeed, D i/HD transition which is important to the total HD abundance is typically sharp (see e.g.Balashev & Kosenko 2020), and hence can be sensitive to the exact geometrical model of the cloud, since the irregular shape, non-uniformity, and patchy structure may allow deeper penetration of the dissociating UV radiation in the medium and making HD abundance less than in case of uniform model.In that sense, one can consider that HD/H 2 likely provides a lower limit on the CRIR, since higher ionization fraction  (Balashev et al. 2019;Klimenko & Balashev 2020;Kosenko et al. 2021, and references therein) binned by the metallicity, where systems in each bin shown in rows, with mean metallicity within the bin provided in the right bottom corner of the panel.The calculated dependencies using the described formalism are presented by solid, dashed, and dotted lines, that correspond to the variation on one of the parameters from the base model with  tot H = 100 cm −3 ,  = 10 −15 s −1 and  = 1.The arrow in the middle panel at the second row marked the position of J 0812+3208, whose detailed analysis is presented in Figure .3 if we assume higher CRIR (and higher ionization fraction) D i/HD transition occurs at a lower depth in the cloud.In turn, as we note above, all, C i/H 2 abundance, C i fine-structure excitation and the temperatures used in our study are quite homogeneous within the H + abundance, and hence the HD formation rate scales with the hydrogen ionization fraction in the medium and CRIR.medium and do not reveal large variations that can significantly bias results.
As we discussed above, the ionization and the thermal state depend on the dust abundance, since the high metallicities, log  ≳ −0.5, the photoelectric heating is the dominant heating source and the electron and H + recombination are tightly coupled with the PAH abundances.
In turn, at low metallicities, log  ≲ −0.5, dust (as it scales with for the H 2 /C i absorption system at z=2.626 towards Q 0812+3208.The blue lines and stripes show the point and interval estimate from the observations.The kinetic temperature, and column density of C i fine-structure levels were used to calculate a likelihood function, while the total H 2 column density and metallicity were used as model parameters, with priors corresponding to the measured interval estimate.metallicity) is less important for both the ionization and thermal state and they are mostly governed by cosmic rays.Additionally, at low metallicities, the electron production rate becomes insensitive to metal abundance, since it is dominated by hydrogen ionization, rather than carbon photoionization as at high metallicity.Since our sample probes a large range of metallicities, we have an opportunity to explore this metallicity dependence.In Fig. 5 we plot the dependence of the derived CRIR on the metallicity.One can see that our results do not show any significant correlation with metallicity, which may indicate that the results are robust against artificially induced effects regarding the onset of the cosmic ray dominance at low metallicities.However, in principle, this cannot be separated from the possible dependence of the physical conditions on the metallicities, which is expected since different metallicities attributed to different galaxies and their evolutionary states.Indeed it is known, that the metallicity of the DLAs correlates with the mass of the galaxy and its halo.The previous constraints based on HD/H 2 abundance also do not reveal a strong correlation, although as we discussed above, have systematically lower values of the derived CRIR.
In Fig. 6 we plot the dependence of the derived UV flux on the hydrogen number density.One can see again, that there is no strong correlation between these quantities.However, one can notice that while UV field estimates are found to be consistent, we derive systematically higher values of the hydrogen number density in our study, than were previously obtained using only excitation of C i fine-structure levels and kinetic temperature (note that values found by Kosenko et al. 2021 generally agree with what was obtained by Klimenko & Balashev 2020).The main discrepancy can arise from the fact that we also used C i abundance measurements and C i formation rate scales with number density, and hence the higher number density is favored to reproduce the observed relatively high C i/H 2 ratios.Additionally, we note, the results obtained by Kosenko et al. 2021;Klimenko & Balashev 2020 based on the modelling by MEUDON PDR code, during which they used a fixed CRIR to minimize the number of the model variables.As we discussed above, the CRIR may significantly impact the thermal and chemical balance at low metallicities and therefore can elevate the characteristic number density of the CNM, since it affects the position of CNM branch of the neutral phase diagram (see e.g.Bialy & Sternberg 2019;Balashev et al. 2022).Additionally, as was noted by Slava Klimenko (private communication, Klimenko et al. in prep), the previous version of MEUDON PDR code that was used by Klimenko & Balashev 2020;Kosenko et al. 2021 had incorrect treatment of CMB field, which intensity was a factor of 2 higher than it should be.Higher values of the CMB field, will increase the contribution of CMB to the excitation of C i levels at high redshifts, where the CMB temperature is enhanced by (1+z) factor, hence lowering the number density estimates.The detailed comparison of C i/H 2 and HD/H 2 modelling certainly deserves an additional study.However, taking into account the aforementioned issues, we suppose that the results presented in that study based on the C i/H 2 abundance are more robust, than the previous ones (Klimenko & Balashev 2020;Kosenko et al. 2021) based on analysis with MEUDON PDR code.

CONCLUSIONS
We present the modelling of the relative abundance of C i/H 2 in the diffuse cold neutral medium using the solution of the H i/H 2 transition, chemical network (with an explicit treatment of the PAH), and excitation of fine-structure levels coupled with the thermal balance.This allows us to accurately describe the dependence of C i/H 2 on the global physical parameters, metallicity, number density, UV field, and cosmic ray ionization rate.We show that typically estimated values of these parameters can reproduce the observation from C i/H 2 -bearing DLA systems at high redshifts.These systems mostly represent the low-metallicity gas, for which the cosmic rays start to dominate the heating and ionization state of the medium, and hence severely affects the chemistry of the gas.We used C i/H 2 abundances, the population of C i fine structure, and H 2 rotational levels measured in DLAs to get constraints on global physical parameters.We found for the sample of all known C i/H 2 -bearing DLAs at high redshifts that CRIR is in the range ∼ 10 −16 to few × 10 −15 s −1 , hydrogen number density is in the range ∼ 10 − 10 3 cm −3 , and UV field is in the range 10 −2 to few × 10 2 of Mathis field.We did not find any strong correlation between obtained parameters.We also compare measured values with previous constraints obtained using only excitation of C i and H 2 levels, and HD/H 2 relative abundance.We argue that the method proposed in this paper provides more robust estimates on the physical conditions, since it takes into account the importance of CRIR and is based on the observational quantities, which are quite homogeneous within the medium, i.e. less affected by the radiative transfer effect, as e.g.HD/H 2 do.

Figure 1 .
Figure 1.The results of the calculation of the cloud structure with  tot H = 100 cm −3 ,  = 1,  = 10 −15 s −1 , and  = 0.3.The panels show profiles of: (a) temperature; (b) abundances of the main elements; (c) number densities on the C fine-structure levels; (d) PAH adundances; (e) the rate of the heating (dashed lines) and cooling (solid lines); (f, g, and h) the rate of reactions of the formation (solid lines) and destruction (dashed lines) for C, PAH − , and , respectively.

Figure 3 .
Figure 3.Posterior distribution functions of the constrained model parameters  tot , ,  , and DTG (along with the distributions of the measured quantities) for the H 2 /C i absorption system at z=2.626 towards Q 0812+3208.The blue lines and stripes show the point and interval estimate from the observations.The kinetic temperature, and column density of C i fine-structure levels were used to calculate a likelihood function, while the total H 2 column density and metallicity were used as model parameters, with priors corresponding to the measured interval estimate.

Figure 4 .Figure 5 .
Figure 4.The estimations of cosmic ray ionization rate,  , and UV field,  in diffuse CNM at high redshift H 2 -bearing DLAs.The red circles and green squares represent the results obtained based on C i/H 2 (this work) and HD/H 2 (from Kosenko et al. 2021) abundances, respectively.The gray symbols depicted the upper limits.

Figure 6 .
Figure 6.The estimates of UV field, , and hydrogen number density,  tot H , in diffuse CNM at high redshift H 2 -bearing DLAs.The color code is the same as in Fig. 4.

Table 1 .
Observational data on C i/H 2 -bearing absorption systems at high redshifts ( ≳ 2) used in this work 4(and hence higher CRIR) is necessary if D i/HD transition is shifted deeper into the cloud.Conversely, The dependence of the column densities of C i on H 2 .The colored points represent observed abundances in high-redshift DLAs

Table 2 .
The constraints on the physical conditions in high redshift C i/H 2bearing DLAs obtained from the joint modelling of C i/H 2 /H i column densities, excitation of C i fine-structure, and lower H 2 rotational levels.