Evidence of First Stars-enriched Gas in High-redshift Absorbers

The first stars were born from chemically pristine gas. They were likely massive, and thus they rapidly exploded as supernovae, enriching the surrounding gas with the first heavy elements. In the Local Group, the chemical signatures of the first stellar population were identified among low-mass, long-lived, very metal-poor ([Fe/H]<-2) stars, characterized by high abundances of carbon over iron ([C/Fe]>+0.7): the so-called carbon-enhanced metal-poor stars. Conversely, a similar carbon excess caused by first-star pollution was not found in dense neutral gas traced by absorption systems at different cosmic time. Here we present the detection of 14 very metal-poor, optically thick absorbers at redshift z~3-4. Among these, 3 are carbon-enhanced and reveal an overabundance with respect to Fe of all the analyzed chemical elements (O, Mg, Al, and Si). Their relative abundances show a distribution with respect to [Fe/H] that is in very good agreement with those observed in nearby very metal-poor stars. All the tests we performed support the idea that these C-rich absorbers preserve the chemical yields of the first stars. Our new findings suggest that the first-star signatures can survive in optically thick but relatively diffuse absorbers, which are not sufficiently dense to sustain star formation and hence are not dominated by the chemical products of normal stars.


INTRODUCTION
Cosmological simulations show that the first (PopIII) stars are likely more massive than present-day "normal" stars, with a characteristic mass of ∼ 10 M and a maximum-mass possibly extending up to ∼ 1000 M (e.g. Hosokawa et al. 2011;Hirano et al. 2014). Among such a variety of stellar masses there are many channels to produce supernovae (SNe), and thus to contaminate the surrounding environment with the heavy elements newly produced by Pop III stars. Very massive first stars, 140M ≤ M P opIII ≤ 260M , explode as Pair Instability SNe yielding chemical abundance ratios that exhibit a strong odd-even effect (Heger & Woosley 2002;Takahashi et al. 2018) and a unique lack of cobalt and zinc over iron (Salvadori et al. 2019). First stars with intermediate masses, 10M ≤ M P opIII ≤ 100M , also evolve as SNe but they can have a variety of explosion energies, thus yielding very different chemical element ratios depending upon the mass of the progenitor star and the SN explosion energy (e.g. Heger & Woosley 2010;Limongi & Chieffi 2018 for theory and Placco et al. 2021;Skúladóttir et al. 2021 for observations). Low-mass stars, born from the first stars-polluted gas (Bromm & Loeb 2003;Schneider et al. 2003), can survive until the present-day and retain in their atmospheres a record of the chemical elements produced by these first stars.
The search for the chemical signatures of Pop III stars has focussed on ancient metal-poor stars in our cosmic neighbourhood. In particular, stars in the Milky Way (MW) halo and Local Group dwarf galaxies are prime targets as we can uniquely study individual stars (e.g. Beers & Christlieb 2005;Frebel & Norris 2015;Simon 2019). Historically a new class of Carbon-Enhanced Metal-Poor (CEMP) stars was first recognized from the observations of Beers et al. (1992). A few years later two simultaneous studies, by Norris et al. (1997) and Bonifacio et al. (1998), identify for the first time an object pertaining to the class of CEMP-no stars, i.e. stars that are very metal-poor ([Fe/H]< −2), strongly enhanced in carbon with respect to iron ([C/Fe]> +0.7), and not enriched in neutron capture elements ([Ba/Fe]< 0.0). At the moment many CEMP-no stars have been discovered (e.g. Christlieb et al. 2002;Rossi et al. 2005;Bonifacio et al. 2015;Norris et al. 2013;Zepeda et al. 2022;Aguado et al. 2022) and very recently it has been confirmed by Aguado et al. (2023) that these objects are most likely the descendants of massive first stars that exploded as low-energy supernovae (e.g. Iwamoto et al. 2005;Marassi et al. 2014).
Indeed, when the explosion energy of a supernova is not high enough to expel Fe-peak elements from the innermost layers, a large fraction of them fall back onto the remnant (e.g. Heger & Woosley 2010). During these faint SN explosions, therefore, only the outermost layers rich of carbon and other light elements are ejected, yielding large values of [C/Fe]. The idea that CEMP-no stars formed in an environment polluted by low-energy primordial SNe is further supported by the increasing frequency of CEMP-no stars towards lower [Fe/H] (e.g. Marsteller et al. 2005;Beers & Christlieb 2005;Lucatello et al. 2006;Lee et al. 2013;Placco et al. 2014;de Bennassuti et al. 2017;Yoon et al. 2018;Liu et al. 2021).
Among the very metal-poor stars enriched with carbon another population exists: stars that exhibit an excess in heavy elements formed by slow (or rapid) neutron capture processes, dubbed CEMP-s (or -r) stars. The available data suggest that most CEMP-s stars dwell in binary systems (Lucatello et al. 2005;Starkenburg et al. 2014;Hansen et al. 2016a;Arentsen et al. 2019) and thus that their carbon-excess is not inherited from the natal cloud but acquired via mass transfer. The surplus of carbon likely comes from an evolved star that has passed through the asymptotic giant branch (AGB) phase (Abate et al. 2015), during which s-elements are also produced (Karakas & Lattanzio 2014). On the other hand, the carbon-excess in CEMP-no stars is expected to be representative of the environment of formation (e.g. Hansen et al. 2016b;Zepeda et al. 2022), even in the rare case in which CEMP-no stars are found to dwell in binary systems (Aguado et al. 2022(Aguado et al. , 2023. Based on these results from ancient nearby stars, we expect that at high redshifts it could be possible to find very metal-poor gaseous environments primarily enriched by the first stars (Pallottini et al. 2014), thus showing a carbon-excess.
Quasar absorption lines provide an important gateway to infer observational constraints on galaxy formation and evolution and to look for the signatures of the first stars in gas at high-redshifts. Detecting gas exhibiting similar abundance patterns as CEMP-no stars would open a new window to investigate the properties of the first stars and galaxies in the early Universe. Yet, despite long searches at z > 2 − 3, this distinctive chemical signature has not been discovered in dense absorbers (Cooke et al. 2011b;Dutta et al. 2014), such as damped Lyman-alpha systems (DLAs). These DLAs, which have neutral hydrogen column density log(N HI /cm −2 ) > 20.3, trace most of the neutral gas in the Universe together with the galaxies interstellar medium (ISM).
A claim of detection of a CEMP-DLA (QSO J0035-0918) at z abs = 2.340 with [Fe/H]= −3.04 and [C/Fe]= +1.53 was published by Cooke et al. (2011a). However, two subsequent works have disproved the result reporting a much lower carbon-to-iron ratio for the same DLA absorption system, i.e. respectively equal to [C/Fe]= +0.51 ± 0.10 and [C/Fe]= +0.45 ± 0.19 (Carswell et al. 2012;Dutta et al. 2014). Another DLA with [C/Fe]= +0.59, [Fe/H]= −2.84 at z = 3.07 has been reported by Cooke et al. (2012), which is however below the limit value of [C/Fe]= +0.7. Two recent work Welsh et al. (2019Welsh et al. ( , 2022 presenting a collection of all the very metal-poor DLAs in the literature confirm the absence of carbon-enhancement claims. Aiming at detecting the chemical evidence of gas enriched by the first stars, Lyman limit systems (LLSs) and sub-damped Lyman-α systems (sub-DLAs), with 17.2 ≤ log(N HI /cm −2 ) ≤ 20.3, represent promising gaseous environments to look for the fingerprints of PopIII stars. Indeed, they are less dense than DLAs and therefore metal poorer (Fumagalli et al. 2016) and likely not strongly contaminated by subsequent generations of normal (Pop II) stars (Salvadori & Ferrara 2012), which are expected to form early on in the ISM of Pop III enriched galaxies. On the other hand, these systems trace optically thick, relatively diffuse gas that is not sufficiently dense to self-shield the UV radiation; consequently, they are likely characterized by more complex ionization patterns.
In this work, we exploited the XQ-100 quasar legacy survey (López et al. 2016) to collect a sample of 54 absorption systems at redshift z ∼ 3 − 4 selected by the presence of the Mg II absorption doublet. Among these systems, we identified a sub-sample of 37 diffuse optically thick LLSs and sub-DLAs absorbers that we have studied in detail. We performed Voigt profile fitting of metal absorption features and hydrogen Lyman lines in the quasar spectra to measure column densities. To derive relative abundances of different elements, we applied photoionization-model corrections to the measured ionic abundances.
In Sect. §2 we present our data-set including the line profile fitting of the absorption systems and the determination of the chemical abundances. In Sect. §3 and §4 we present and discuss the results. The conclusions are drawn in Sect. §5.

DATA ANALYSIS
The Large Programme "Quasars and their absorption lines: a legacy survey of the high-redshift Universe with VLT/X-shooter" (López et al. 2016) has produced a homogeneous and high-quality sample of echelle spectra of 100 quasars (QSOs) with emission redshift z ∼ 3.5−4.5. The targets were observed with the X-shooter spectrograph (Vernet et al. 2011) mounted at the ESO Very Large Telescope (VLT, Cerro Paranal, Chile). X-shooter is characterised by three arms that allow to cover in one observation the full spectral range between the atmospheric cutoff at 300 nm and the near-infrared K-band at 2500 nm, at an intermediate resolving power. The full spectral coverage, along with a well-defined target selection and the high signal-to-noise achieved (median SN R = 30), clearly makes XQ-100 a unique data-set to study the rest-frame UV/optical spectra of high-z QSOs in a single, homogeneous, and statistically significant sample. The adopted slit widths were 1.0 in the UVB arm and 0.9 in the VIS and NIR arms, to match the requested seeing and to account for its wavelength dependence. These slit widths provide nominal resolving powers of 5400, 8900 and 5600 for the UVB, VIS and NIR arm respectively 1 . The XQ-100 survey was designed to cover many science cases: from the detailed study of the intergalactic medium, to the detection of galaxies in absorption Berg et al. 2016;Christensen et al. 2017), from the properties of QSO themselves (Perrotta et al. 2016(Perrotta et al. , 2018 to cosmology (Iršič et al. 2017a,b). The spectra, reduced from the collaboration, were delivered to the public 2 . Two types of reduced data are provided for each target: (i) a joint spectrum of the three arms together; (ii) individual UVB, VIS, and NIR arm spectra, which includes telluric correction and fitted QSO continuum. When a target is observed more than once at different epochs, there is one spectrum for each epoch and a combined spectrum, putting together all epochs. The telluric corrections are calculated on the spectra corresponding to the single epochs, while the continuum is calculated only on the combined spectrum.
Since for our purposes we needed telluric corrected spectra, we created a coadded spectrum when multiple observations were present and then we re-determined the intrinsic QSO continuum in the framework of the Astrocook 3 Python software package (Cupani et al. 2020). In Astrocook, the emission continuum is estimated by first masking the most prominent absorption features and then interpolating the non-masked regions with a univariate spline of chosen degree.
The final spectrum for each QSO was created by cutting the noisy edges of each arm and "stitching" the three arms together. Also these operations were carried out using Astrocook.
Before creating the final spectra, we estimated the "effective" resolving power determined by the atmospheric conditions during observations. Indeed, if the seeing during observations is smaller than the width of the slit, the "effective" resolving power of the obtained spectrum will be larger than the nominal one. We recomputed the value of the resolving power (R) for each spectrum (Table 1) based on the average value of the seeing during the observations (reported in the ESO archive as DIM M ) and assuming a linear relation between resolving power and slit width. For example in the VIS arm, if DIM M < 0.9 arcsec, the new resolving power is obtained as R ef f ∼ (0.9/ DIM M ) × R nom (see D'Odorico et al. 2022). Note however that these 1 The nominal resolutions are different from those of López et al. (2016) since at that time, due to a problem in X-shooter's data reduction, the resolutions were calculated incorrectly. https:// www.eso.org/sci/facilities/paranal/instruments/xshooter.html 2 https://www.eso.org/qi/catalog/show/73 3 https://github.com/DAS-OATs/astrocook determinations have uncertainties of the order of 10% (see D'Odorico et al. 2023, in prep.).

Line Fitting
The analysis of the spectra has been carried out with the Astrocook software package. In Astrocook, absorption lines are detected as prominent local minima in the flux density spectrum, then they are identified by crossmatching their wavelengths with a list of ionic transitions commonly observed in QSO spectra and finding coincidences among the obtained redshift values.
For our study, we need to identify absorption systems with H I column densities in the range 17.2 ≤ log(N HI /cm −2 ) ≤ 20.3, i.e. LLSs and sub-DLAs, which can trace diffuse gas, such as the one in the outskirts of galaxies and in cosmic filaments (Lofthouse et al. 2022). To this aim, we searched the XQ-100 spectra for singly ionised magnesium doublets, Mg II, which are good probes of the optically thick, low ionisation gas. Mg II is one of the best known examples of strong resonance-line doublets: it has rest-frame wavelengths (see Table 3) longer than H I Lyman-α (Lyα), and therefore it appears on the red side of the Lyα emission line in the QSO spectrum. For this reason, it is relatively easy to identify as it is not embedded in the thick Lyα forest. The search for Mg II doublets was carried out using an automatic recipe of the Astrocook software. Subsequently, the detection of the absorption lines was confirmed visually. In our analysis we have not make any attempt to reach or determine the completeness of the absorber sample since it was not relevant for our scientific purpose.
In each line of sight, we restricted our search to a redshift range that avoids the proximity region of the quasar (5000 km s −1 from the quasar emission redshift) and that allows to cover the H I Lyα and Lyman-β (Lyβ) transitions of a given absorption system, for a more reliable determination of the H I column density. Furthermore, we have excluded the interval 13500-14500Å affected by strong telluric lines. This wavelength range corresponds to the redshift range z 3.83 − 4.18 for Mg II λ2796Å.
Detected systems are then modelled with Voigt profiles in the context of Astrocook. The Voigt profile fitting provides the central line redshift, z, the column density, N , and the Doppler broadening parameter, b. After the identification of the Mg II doublets, we proceeded with the search of other low ionization lines (see Table 3) at the same redshift. Assuming that they originate from the same gas and that turbulent motion is dominant over thermal one, we fitted all with the same redshift components having the same Doppler parameters. On the other hand, C IV and Si IV absorption lines, if present, have generally a different velocity structure and they were fitted separately. To better constrain the chemical properties of our absorbers we also estimated column density upper limits, in particular for Fe II and O I if they were not detected. We used equations 2 and 3 of D' Odorico et al. (2016) adapted for 1σ limits assuming a Doppler broadening parameter equal to the one of the other low ionisation transitions in the system and the spectral velocity bin (UVB: 20 km s −1 ; VIS: 11 km s −1 ; NIR: 19 km s −1 ) from López et al. (2016). The neutral hydrogen column density was determined by considering in the fit all the lines in the Lyman series free from strong blending.
A well known issue related to absorption line fitting is the identification of saturated lines. Indeed, with the resolving power provided by X-shooter, in many cases the lines are not resolved and therefore it is difficult to understand when they are saturated. We adopted the following empirical procedure for metal lines: when the difference between the normalized flux and the depth of the line was >0.7 we assumed that the line was saturated and considered the determined column density as a lower limit. In this hypothesis we are neglecting the possible dependence on the Doppler parameter value.
We selected 51 Mg II absorption systems plus three previously identified LLS systems whose Mg II absorption fall in the telluric band (S. Cristiani priv. comm.), for a total of 54 absorption systems. We compared our sample with those of prior studies of the XQ-100 survey investigating DLAs and Sub-DLAs (Berg et al. 2016(Berg et al. , 2017(Berg et al. , 2019(Berg et al. , 2021, finding matches for 17 DLAs and 16 Sub-DLAs. We excluded the DLAs from our sample and used the HI column density determined in Berg et al. (2019) for the sub-DLAs. Our analysis was then focused on a sample of 16 Sub-DLAs and 21 LLSs.
Three of the QSOs in our sample (J0247-0556, J1111-0804 and J1723+2243) have a reduced UVES spectrum available from the SQUAD database . We verified with Astrocook that for our systems there is a good agreement between column densities measured in the X-shooter and UVES spectra.

Determination of chemical abundances
The gas in LLSs and Sub-DLAs is not fully neutral, therefore it is necessary to apply a ionisation correction  to translate the observed ionic column densities into element abundances and thus, into gas metallicity.
To infer the chemical composition and the physical state of the absorbing gas, we have used ionisation models based on radiative transfer calculations at equilibrium and for a single gas phase. These calculations are the input of a Bayesian formalism that exploits Markov Chain Monte Carlo (MCMC) techniques to derive the posterior probability distribution function for quantities of interest, such as the metallicity, Z, and the physical density, n H , of the absorbing gas. In particular, the MCMC method we have adopted (Fumagalli et al. 2016) uses a grid of Cloudy models. Cloudy (v. 17, Ferland et al. 2017) is an open-source photoionization code, that maps each set of input parameters into the corresponding column densities of all metal ions that we considered in this work. The MCMC sampler searches Z and n H for models with the column density pattern that best matches the neutral hydrogen column density, N HI , and the metal ion column densities fitted by the software Astrocook. The output of the MCMC modelling is a probability distribution of the sets of parameters Z and n H . In this work, we adopt the minimal model parameters (Fumagalli et al. 2016), which assumes a slab of gas at constant density illuminated on one side by both the UV background (Haardt & Madau 2012) and the cosmic microwave background. All metals are assumed to be in the gas phase with a solar abundance pattern (Asplund et al. 2009).
To obtain the ionization corrections ion by ion, we run again Cloudy system by system optimizing the values of Z and n H determined by the MCMC run based on the observed column densities. We vary the parameters in an interval corresponding to the tenth and ninetieth percentile of their posterior distribution function. We considered a narrow range because the optimization on the full parameter space had already been performed by the MCMC code. Finally, from the photoionization model we obtained the corrections for each ion of each element. Once the column densities obtained from the Voigt fit have been corrected for ionization, we derived the absolute abundances of the various elements. The corrected column density of an element X has been determined using the formula N (X) = N (X i )/IC(X i ), where N (X i ) is the column density of the ion fitted from the spectrum and IC(X i ) the ionization correction obtained from Cloudy. Then, we derived the relative abundances of each element as: where N H and N X are respectively the column density of hydrogen and of a specific element X and log(N X /N H ) is the solar abundance (Asplund et al. 2009). In Fig. 1  [Fe/H]< −2, where we include also those systems for which [Fe/H]−1σ < −2. We have verified that these [Fe/H]< −2 systems are also characterised by a low total metallicity with respect to the solar value (see Fig.  6 and Tab. 2) and therefore that they are not just ironpoor because of dust depletion. On the other hand, the shape of the MDF at [Fe/H]> −2 is less constrained due to the fact that the iron measurements could be affected by dust depletion: the effect would be to increase the value of [Fe/H].
Note that the chemical abundance pattern of very metal-poor stars, at [Fe/H]< −2, is well established to be different from the solar value (e.g. Cayrel et al. 2004;Bonifacio et al. 2009;Yong et al. 2013). For this reason, we recompute the ionisation corrections for our very metal-poor absorbers by assuming the chemical abundance pattern derived by Cayrel et al. (2004) for a sample of C-normal giant stars with −4 <[Fe/H]< −2 in the Milky Way halo (hereafter we will refer to this pattern as Cayrel). This average chemical abundance pattern is indeed characterised by a very small star-to-star scatter. Hence, it can be considered as a "reference" for very metal-poor environments. Table 4 reports the values assumed in the new Cloudy models which are those derived by Cayrel et al. (2004) corrected for 3D and/or Non-Local Thermodynamical Equilibrium (NLTE) effects and for stellar physical processes that can alter the measured abundances. These corrections are required to use the stellar abundance values for the gaseous component. In particular we adopted: (i) the [C/Fe] value derived for dwarfs (Bonifacio et al. 2009), since the photospheric carbon abundances in giant stars can be altered by convective mixing. In other words we assumed that the difference between giant and dwarfs is due to the first dredge-up; (ii) the original [O/Fe] value derived by Cayrel et al. (2004), since Bonifacio et al. (2009) demonstrated that in giant stars 3D effects for oxygen are not important; (iii) the magnesium value derived by Andrievsky et al. (2010) Cayrel et al. (2004), since aluminium was already corrected for NLTE effects and NLTE effects for silicon are expected to be very weak.

in order to account for NLTE effects; (iv) the [Al/Fe] and [Si/Fe] values reported by
We defined a modified chi square to establish which model (solar or Cayrel) reproduces better our observed column densities, χ 2 = ( i ((N i model − N i obs )/σ) 2 ) 1/2 , where N model and N obs are the modelled and observed ionic column density values and σ is the relative error in the observed value. More than 75% (11/14) of the very metal-poor absorbers are better modeled with Cayrel's abundance pattern than with solar (see Table 5).
The observed column densities are better reproduced for most of the ions. In addition, the relative abundances obtained with the new photoionization model are consistent with the abundances obtained with the solar model for most chemical elements.
The results of the Voigt profile fitting for the 14 very metal-poor absorption systems are available online as supplemental material. The tables reporting the fit parameters are available for all the 37 absorption systems. 4 An example figure of two metal-poor absorption systems is shown in the Appendix together with a portion of the fit parameters table for guidance. Figure 1 shows the new MDF obtained for our absorbers, adopting the Cayrel abundances for very metalpoor systems and the solar values for those with [Fe/H]> −2. We can first note that the shape of the MDF remains essentially unvaried with respect to the one which was obtained by assuming solar values for all absorbers. This is because most of the changes are small, with a [Fe/H] difference < 0.4 dex, i.e. smaller than the width of the MDF bin. By further inspecting Fig. 1 we see that the MDF of our gaseous absorbers is roughly bimodal: the first broad peak is around [Fe/H]∼ −0.8 while the second one, which is more pronounced, appears at [Fe/H]∼ −2. This result implies that these absorbers are likely a variegated population (e.g. Fumagalli et al. 2016).

Metallicity Distribution Function
We compare the normalized MDF of our high redshift absorption systems (Cayrel model) with those of ancient very metal-poor stars observed in the Galactic halo (Bonifacio et al. 2021) and in Local Group ultrafaint dwarf galaxies (UFDs) (Simon 2019), respectively in the upper and lower panel of Fig. 2. We are aware that we are comparing the chemical abundances of diffuse high redshift gas, which can be located in the outskirts of galaxies, with the ones of ancient local stars. Still, the main point behind our comparison is that stars are born from gas.
The MDF of our sample spans a wide range of [Fe/H], covering the values measured in stars of both the MW stellar halo and nearby UFDs. The distribution of the halo stars covers a similar range of [Fe/H] with respect to our absorbers, but has a different shape. Indeed the stars of the Galactic halo have a unique peak, which is more pronounced, and appears around [Fe/H]∼ −1.
Note that, although very rare and hence almost invisible in the normalised MDF, stars with [Fe/H]< −2.5 in the Galactic halo are those showing the chemical imprint of the first stars (see Fig. 4). The low-Fe peak of our absorbers, therefore, suggests similarities between these high-z absorbers and the gas that may have hosted the imprint of the first stars, like the birth environment of very metal-poor halo stars.
The MDF of UFDs covers a narrower range of iron abundances, −4 <[Fe/H]< −1, than the MDF of our absorption systems and shows a peak at [Fe/H]= −2.5, almost overlapping with the low-Fe peak of the MDF of our absorbers (Fig. 2, lower panel). These findings suggest that our absorption systems might represent the gas-rich counterpart of UFDs at high redshift (Salvadori & Ferrara 2012;Skúladóttir et al. 2018). The absence of ultra metal-poor absorption systems ([Fe/H]< −4) which are observed among Galactic halo stars, could be due to observational biases. On the one hand, we select our systems by the presence of the singly ionised magnesium doublet, on the other, the resolution and signal-to-noise ratio of the XQ-100 spectra do not allow us to put very stringent upper limits on the iron column density, when the ionic lines are not detected.

Carbon-enhanced systems: stellar relics vs high-z absorbers
In Fig. 3 we compare the carbon-to-iron ratio, [C/Fe] measured in our very metal-poor gaseous systems and in Local Group stars (halo and UFDs) as a function of [Fe/H]. The stellar data we used were taken from Salvadori et al. (2015) and updated with new measurements for stars in UFDs (Ji et al. 2016;Spite et al. 2018) and newly discovered extremely metal-poor stars in the Milky Way (Starkenburg et al. 2018;François et al. 2018;Bonifacio et al. 2018;Aguado et al. 2019;González Hernández et al. 2020). All [C/Fe] values are corrected to account for internal mixing processes (Placco et al. 2014). We see that our very metal-poor absorbers exhibit the same trend in [C/Fe]  We also notice that among the 14 very metal-poor absorbers, 3 are carbon-enhanced, with [C/Fe]> +0.7. These C-enhanced absorbers have an upper limit on Fe II, implying that their true [C/Fe] could be even larger. Note that there are other two absorption systems with [C/Fe]≈ +0.7. Still, since the definition of C-enhanced very metal-poor absorbers varies in the literature, e.g. [C/Fe]≥ +0.7 (e.g. Beers & Christlieb 2005) and [C/Fe]≥ +1.0 (Bonifacio et al. 2015) we decided to include them in the C-normal sub-sample. For the 3 carbon-enhanced systems we checked the reliability of the column density measurements against our hypothesis of turbulent broadening with respect to thermal broadening. To this end, we fit the detected lines (Si II, C II and Al II) assuming they arise in the same gas and they are thermally broadened and we find negligible variation of the column density ( ≤ 0.02) for all the transitions and ≤ 0.1 for the C II in J1658-0739 at z=3.54604.  (Salvadori et al. 2015;Placco et al. 2019).
When considering this complete stellar sample, the excellent agreement between the chemical abundances measured in stars and gaseous absorbers at [Fe/H]< −2 is even more evident, for both C-enhanced and C-normal systems. In particular, we see that very metal-poor Cenhanced absorbers never overlap with CEMP-s stars, which are shifted towards higher [Fe/H] at any given carbon-to-iron ratio. Conversely, we notice that systems with [Fe/H]> −2 cover a wider range of [C/Fe] abundances than what is observed both in very metal-poor absorbers and in stars at the same [Fe/H]. Among these [Fe/H]> −2 absorbers, we observe four C-enhanced systems, three of which nicely overlap with CEMP-s stars. We also see many (9) C-normal absorbers that reside in the same region of C-normal stars. At [Fe/H]> 0, where a few star measurements are available, we notice one C-enhanced absorber and two systems that are strongly C-deficient with respect to stars ([C/Fe]< −1). The large [C/Fe] scatter of the absorbers in this [Fe/H]≥ −1 regime, along with the presence of C-deficient systems, suggest that a non-negligible amount of carbon (and iron) is likely depleted onto dust grains. Ultimately our results, which are in line with previous studies (e.g. Vladilo 1998;Quiret et al. 2016;De Cia et al. 2018;Vladilo et al. 2018), confirm that the dust contribution makes the comparison between the chemical abundances of gas and stars challenging at [Fe/H]≥ −1, while it can be neglected in the very metal-poor regime. Armed with these new findings we can make further comparison between the chemical abundances of gas and stars in the very metal-poor regime.

Other chemical elements
To investigate more deeply the chemical enrichment history of our C-enhanced very metal-poor absorbers, in Fig. 4 we compare the relative abundances of oxygen, magnesium, and silicon with respect to iron, with those of ancient very metal-poor stars. The stellar data we used for these elements were taken from the SAGA database (Suda et al. 2008). The sample of stars is different from the one used for carbon-to-iron ( Fig. 3 and upper left panel of Fig. 4). For this reason we do not distinguish between CEMP-no and CEMP-s/r stars.
In Fig. 4 we see that there is a general agreement between the chemical abundances of our absorption systems and the abundances of the stars in the MW and dwarf galaxies. Furthermore, at [Fe/H]< −2, we note the same trend for all the chemical elements in both stars and gas: the relative abundances increase as the [Fe/H] decreases. In particular, we see that C-enhanced of all our gaseous absorption systems (filled circles) and of stars in ultra-faint dwarf galaxies (light blue squares), dwarf spheroidal galaxies (blue squares) and the Milky Way (grey/black squares). CEMP-no halo stars are shown as filled black symbols, C-normal halo stars as filled bordered light grey squares and open symbols are CEMP-s/r stars. Stellar measurements have typical 1σ errors of 0.2 dex. We distinguish among C-enhanced (red) and C-normal (yellow) very metal-poor absorbers. Orange circles are iron-rich absorbers. Other panels: Comparison between the chemical abundances of the very metal-poor gaseous absorption systems (circles) and of the stars (squares) in the Milky Way (grey), ultra-faint dwarf galaxies (light blue) and dwarf spheroidal galaxies (blue). Red and yellow circles are respectively C-enhanced and C-normal, very metal-poor absorbers. Orange circles are iron-rich absorbers. Stellar data are taken from the SAGA database (Suda et al. 2008). In all the panels the orange shaded area identifies the region where chemical abundances could be affected by dust depletion. Figure 5. Absolute carbon abundance (A(C)) as a function of [Fe/H] of all our gaseous absorption systems (filled circles) and of stars in ultra-faint dwarf galaxies (light blue squares), dwarf spheroidal galaxies (blue squares) and the Milky Way (grey/black squares). CEMP-no halo stars are shown as filled black symbols, C-normal halo stars as filled bordered light grey squares and open symbols are CEMP-s/r stars. We distinguish among C-enhanced (red) and C-normal (yellow) very metal-poor absorbers. Orange circles are iron-rich absorbers. Stellar measurements have typical 1σ errors of 0.2 dex. We used the same stellar sample of Fig. 4 (upper left panel). The horizontal line separates the low-C from the high-C band, while the dashed one indicates the value [C/Fe]=+0.7. The orange shaded area identifies the region where chemical abundances could be affected by dust depletion.
very metal-poor absorbers are always overabundant also in the other chemical elements. The same is observed in CEMP-no stars (Vanni et al. in prep).
Oxygen (Fig. 4, upper right panel) shows a similar trend to that of carbon even if we have fewer measurements. Indeed, we could measure the oxygen abundance only for 15 of the 30 systems for which we have the C abundance. In Fig. 4, however, we report only 12 measurements since for 3 absorbers we have both upper limits for iron and oxygen. For the remaining systems, the O I absorption lines fall in the Lyα forest and are heavily blended with other absorption lines. Magnesium (Fig. 4, lower left panel), which is produced by massive stars and released in the gas during supernova explosions, shows a similar trend with iron to that seen for [C/Fe], although the values are less extreme. In particular, we see that all our Cenhanced very metal-poor systems are rich in magnesium, [Mg/Fe]> +0.7, as observed in CEMP-no stars (e.g. Frebel et al. 2005;Keller et al. 2014, and Vanni et al. in prep. for a global view). The correlation between C and Mg tells us that these elements are probably produced by the same sources, i.e. most likely primordial low-energy supernovae. At higher iron abundances, [Fe/H]> −2, we see a large scatter. Silicon (Fig. 4, lower right panel) also shows increasing relative abundances with respect to iron as [Fe/H] decreases. Once again, the C-enhanced very metal-poor systems are also enhanced in silicon, showing super solar abundances. Also for silicon, we see that the system-tosystem scatter increases at [Fe/H]> −2, with abundance ratios spanning a wide range of values, many also in the sub-solar regime. This result suggests that silicon is likely one of the chemical elements that is most affected by dust depletion.
We recall that the chemical abundances measured in stars might suffer NLTE effects, whose precise estimate, at the moment, is only available at [Fe/H]< −2 for a small stellar sub-sample and for a few elements. Such effects can lower the [C/Fe] and [O/Fe] values measured in metal-poor stars by more than 0.3 dex (Amarsi et al. 2019) while they can increase the [Mg/Fe] value by +0.4 dex (Andrievsky et al. 2010), thus resulting in a better agreement between our C-enhanced very metal-poor absorbers and CEMP-no stars (see Fig. 4). For [Al/Fe] the correction is extremely high, +0.65 dex (Cayrel et al. 2004), and for this reason we decided not to make the comparison with the values measured in our absorbers.
Ultimately, Fig. 4 demonstrates that in the very metal-poor regime, [Fe/H]< −2, the chemical abundance ratios measured in our gaseous absorbers are in very good agreement with those of present-day stars and that our C-enhanced very metal-poor absorption systems reside in the same regions of CEMP-no stars and have the same chemical properties.

First star signatures
The results described in the previous Sections emphasize the idea that our C-enhanced very metal-poor absorbers are the gaseous analogues of CEMP-no stars, which have been likely imprinted by primordial lowenergy SNe. Indeed, in our absorption systems we see an overabundance of Mg and Si, which is also observed in CEMP-no stars. Magnesium and silicon are key elements since they are produced by primordial low-energy supernovae but not by AGB stars, which are the pollutant of CEMP-s stars. AGB stars yield high C, N and O, and also produce Ba via the slow-neutron capture process.
In CEMP-s stars the C and Ba excess is expected to be acquired via mass transfer from a binary AGB companion, while in CEMP-no stars to be representative of the environment of formation. For this reason barium is used to discriminate between CEMP-s stars ([C/Fe]> +0.7, [Ba/Fe]> +1.0), and CEMP-no stars ([C/Fe]> +0.7, [Ba/Fe]< 0.0, (e.g. Beers & Christlieb 2005)). Unfortunately, no barium measurements are (nor will be) available for our absorption systems. Thus, to further validate the link between our C-enhanced very metal-poor absorbers and CEMP-no stars we need to carry out additional tests.
In stellar archaeology, the absolute abundance of carbon, A(C), displayed as a function of [Fe/H], is used to discriminate among CEMP-no and CEMP-s stars when barium abundances are not available. Indeed, it has been shown (e.g. Spite et al. 2013;Bonifacio et al. 2015;Yoon et al. 2016) that these two different populations dwell in two well separated regions: the high carbon band, A(C) > 7.4, and the low carbon band, A(C) < 7.4. The stars belonging to the low carbon band are mostly CEMP-no stars 5 , while those of the high carbon band are CEMP-s stars. We can thus use this diagnostic to unveil the nature of our C-enhanced absorption systems.
In Figure 5 we display the absolute carbon abundance of Local Group stars and of our absorbers as a function of [Fe/H]. For the absorption systems we computed A(C) = log(N C /N H ) + 12, where N C and N H are the measured column densities corrected for ionisation. We see that our three carbon-enhanced very metal-poor absorbers are found in the low C-band, which is consistent with what is found for CEMP-no stars. These systems are clearly separated from the C-enhanced absorbers at [Fe/H]> −2, which populate the high carbon band where CEMP-s stars reside. In other words, the division between the low-and the high-C band corresponds to a division in [Fe/H] of our C-enhanced absorbers. Hence, despite the lack of barium measurements, we can conclude that our three C-enhanced absorbers at [Fe/H]< −2 are CEMP-no absorption systems, while the four at [Fe/H]> −2 are CEMP-s. Note that the C-enhanced absorber recently discovered by Zou et al. (2020) reside in the CEMP-s region ([Fe/H]= −1.6, A(C)= 9.0). Thus, like for our four CEMP-s absorbers, the C-excess of the absorption system studied by Zou et al. (2020) most likely arise from the contribution of AGB stars, and not from the chemical elements yielded by the first stars.
5 Stars belonging to the low carbon band can be subsequently divided into Group II and Group III, which might have different astrophysical origin (e.g. Yoon et al. 2016Yoon et al. , 2019Zepeda et al. 2022) possibly linked with various dust composition in Pop III supernovae ejecta (e.g. Chiaki et al. 2017). However, given the upper limit on [Fe/H] of our CEMP-no absorbers, we are not able to distinguish between the two groups.

DISCUSSION
All gathered evidence supports the idea that the newly discovered C-enhanced very metal-poor absorption systems are the gaseous z ∼ 3 − 4 analogues of present-day CEMP-no stars. Hence, we propose to define CEMPno absorbers, LLSs/sub-DLAs with [Fe/H]< −2 and [C/Fe]> +0.7. An increasing number of theoretical studies is supporting the idea that CEMP-no stars observed in different environments have been most likely enriched by the first stellar generations (e.g. Iwamoto et al. 2005;Salvadori et al. 2015;Liu et al. 2021). Thus, we are providing the first clues of Pop III stars-enriched gas in high-z absorbers. Still, two questions naturally arise: what is the nature of these high-z CEMP-no absorption systems? And why CEMP-no absorbers have so far escaped detection?

Origin of CEMP-no absorbers
Different semi-analytical models and cosmological simulations have investigated the origin of CEMP-no stars in the Galactic halo (e.g. de Bennassuti et al. 2017;Hartwig et al. 2018;Liu et al. 2021) and in ultra-faint dwarf galaxies (e.g. Salvadori et al. 2015;Jeon et al. 2021, Rossi et al. in prep). We can thus exploit these findings to interpret the properties of our CEMP-no absorption systems. The moderately high [C/Fe] and relatively high [Fe/H] values of CEMP-no absorbers can be explained as the result of two different enrichment mechanisms (Vanni et al. in prep). These mechanisms are: (i) a pollution solely driven by Pop III stars, which explode as low-energy supernovae with different masses (Hartwig et al. 2018;Welsh et al. 2021;Jeon et al. 2021); (ii) an enrichment still dominated by the products of low-energy Pop III supernovae, but which is transiting towards a C-normal pattern due to the contribution of subsequent generations of normal (Pop II) stars exploding as corecollapse supernovae (de Bennassuti et al. 2017;Salvadori et al. 2015;Jeon et al. 2021).
In both cases, the chemical enrichment should be dominated by the chemical products of primordial lowenergy supernovae, which produce the carbon-over-iron excess (Salvadori et al. 2023). However, CEMP-no relic stars likely form during the first Gyr of cosmic evolution (z > 6), while our CEMP-no absorbers are observed at lower redshifts, z ≈ 3 − 4, when the Universe was > 2 Gyr old. What are these CEMP-no absorbers and how can they preserve the chemical signature of Pop III stars?
Two main possibilities exist. The first, is that our CEMP-no absorption systems are associated to the diffuse circum-galactic medium (CGM) of recently born Pop III galaxies at z ≈ 3 − 4. Still, these objects should be extremely rare. Indeed, direct detection of Pop III galaxies is still lacking and at z ≈ 3 − 4 pristine galaxies are expected to be extremely rare (e.g. Pallottini et al. 2014;Jaacks et al. 2019). The second, is that our CEMP-no absorption systems are associated to lowmass "sterile" mini-halos (Salvadori & Ferrara 2012). Radiative feedback processes can indeed increase the gas temperature of metal-enriched mini-halos, which are then too diffuse to form stars but too metal-enhanced to be photo-evaporated. The chemical signature of the first stellar generations can then be preserved in these gasrich systems until they evolve in isolation. This physical mechanism, which has been proposed by Salvadori & Ferrara (2012), should be quite common in the early Universe and can turn low-mass UFDs at the beginning of their evolution into "failed" UFDs: the gas-rich and high-z dark counterpart of ultra-faint dwarf galaxies. Red and yellow squares are respectively C-enhanced and Cnormal, very metal-poor absorbers. Orange squares are ironrich absorbers. The shaded red region shows the metallicity range, according to our results, for gas enriched by first stars. On the right we show the metallicity posterior distribution function of all our absorbers.

Comparison with Literature Data
To examine our results in a global context and understand why CEMP-no absorption systems have not been identified yet, we can compare the total metallicity of all our 37 diffuse absorbers 6 with literature data for LLSs (Fumagalli et al. 2011(Fumagalli et al. , 2016Crighton et al. 2016;Robert et al. 2019). The comparison is shown in Fig. 6, where the total metallicity of the absorbers is displayed as a function of redshift (or cosmic time) and our systems are coloured to distinguish among [Fe/H]> −2 absorbers, CEMP-no systems, and very metal-poor C-normal absorbers (see caption and Fig. 4). Literature data include: (i) the four serendipitously discovered LLSs that have Z < 10 −3 Z and cover the same redshift range of our absorbers, 3 < z < 4.5; and (ii) the sample of LLSs at 0 < z < 4.5 collected by Fumagalli et al. (2016), for which we report the median metallicity computed by the authors using redshift bins containing at least 25 LLSs.
The right panel of Fig. 6 shows the total metallicity distribution function of our absorbers. Firstly, we see that our absorption systems cover a wide metallicity range, 10 −3 Z < Z < Z , and that > 90% of our very metal-poor absorbers, [Fe/H]< −2, have a total metallicity Z < 10 −2 Z (the only exception being one C-normal absorber). This result is a further confirmation that the effect of dust is negligible at [Fe/H]< −2, as suggested by previous works (e.g. Vladilo 1998;Vladilo et al. 2018). Secondly, we note that the peak of the distribution of our absorption systems, Z ≈ 10 −2.5 Z , is in perfect agreement with the average metallicity value of LLSs derived by Fumagalli et al. (2016) in the corresponding redshift bin, 3 < z < 4.5. Then, we note that among our absorbers there is a lack of extremely metal-poor absorption systems, Z < 10 −3 Z , although our systems cover exactly the same redshift range of the most pristine LLSs. This is certainly due to our initial selection which was based on the presence of metal absorption lines (see Sec. 2.1). Finally, we note that CEMP-no absorbers showing the chemical signature of the first stars are not the most pristine objects. This result, which is in line with what is found in the stellar halo, is a consequence of the key chemical signature left by primordial low-energy supernovae: an excess of C (O, Mg, and Si) over iron. To identify the footprint of primordial low energy supernovae in the gas at high-z we should then look for the most ironpoor absorbers, which are not the most metal-poor but rather have 10 −3 Z < Z < 10 −2.3 Z . Among the 12 absorption systems with available carbon and iron measurements at Z < 10 −2.5±0.1 Z , we have that three are CEMP-no, implying that the fraction of CEMP-no systems is F CEM P −no (Z < 10 −2.5 Z ) ≥ 25%. We note 6 Note that the number of absorption systems with iron measurements is smaller (30) than those with metallicity measurements.
that this is just an indicative fraction since we have not determined the completeness of our selection of absorption systems based on the automatic detection of Mg II absorption doublets. Ultimately, we suggest that metal-poor LLSs are the most promising candidates as CEMP-no absorption systems. Indeed, from the one hand, DLAs are likely too dense to stop forming stars and thus are naturally dominated by the chemical products of normal Pop II stars. From the other hand, the most pristine LLSs are too metal-poor to be enriched by the chemical products of primordial low-energy supernovae.

CONCLUSIONS
In this paper, we have studied the chemical abundances of 37 optically thick Lyman-α absorbers (LLSs and Sub-DLAs) at z ∼ 3 − 4.5 identified in the spectra of the XQ-100 quasar legacy survey (López et al. 2016). Column densities of the absorption features have been derived through fitting with Voigt profiles and ionization corrections have been applied in order to derive the chemical abundances. The main results of our study can be summarised as follows: • The MDF of our absorption systems covers a broad range in iron abundance and is bi-modal: it has a broad peak at [Fe/H]∼ −0.8 and a more pronounced one at [Fe/H]∼ −2.
• The low-Fe tail of the MDF of our absorption systems almost overlaps with the stellar MDF observed in UFDs.
• The chemical abundance ratios (C, O, Mg, and Si over Fe) measured in very metal-poor absorbers are in good agreement with those of very metalpoor stars in the Galactic halo and dwarf galaxies.
• Conversely, at [Fe/H]> −2, the absorbers' abundance ratios exhibit a larger scatter than presentday stars, suggesting that dust contribution is no longer negligible at these [Fe/H] values.
• The three C-enhanced very metal-poor absorbers also show an overabundance of Mg and Si, which are produced by first stars exploding as low-energy supernovae and not by AGB stars. These overabundances are also observed in CEMP-no stars.
• All C-enhanced very metal-poor absorption systems have A(C) < 7.4, i.e. they reside in the so-called low carbon band, as it is observed for CEMP-no stars.
• Conversely, all C-enhanced systems at [Fe/H]> −2 dwell in the high carbon band, A(C) > 7.4, consistent with CEMP-s stars.
• CEMP-no absorbers are the most iron-poor among our diffuse systems, but they are not the most metal-poor ones.
• Among the 12 absorption systems with available C and Fe measurements at Z < 10 −2.5±0.1 Z , we have that three are CEMP-no, implying that F CEM P −no (Z < 10 −2.5 Z ) ≥ 25%.
Ultimately, our CEMP-no absorption systems seem to be the gaseous high-redshift analogues of locally observed CEMP-no stars, which have been more likely imprinted by primordial low-energy supernovae. Our new discovery of CEMP-no absorbers suggests that optically thick, relatively diffuse absorption systems are the best environments to identify the missing chemical signature of the first stellar generations in the gaseous component. These absorbers are likely too diffuse to be star-forming, a key requirement to prevent further chemical pollution and thus to preserve the first-star signature. Theoretical investigations and further observational studies are required to fully understand their nature.
In the very near future, extremely large samples of high-redshift quasars will be available from several surveys (DESI Abareshi et al. 2022;WEAVE-QSO Pieri et al. 2016;4-MOST de Jong et al. 2019). They will allow us to select more of these CEMP-no absorbers to statistically characterize their properties and num-ber density. A significant step forward in the analysis of these systems will be represented by ANDES, the high-resolution (R ∼ 100, 000) spectrograph for the Extremely Large Telescope (ELT), foreseen for the beginning of the 30s (Marconi et al. 2022). The ANDES spectral coverage extending to the NIR and the collecting power of ELT will allow us to carry out detailed studies of these systems resolving the metal absorptions and determining significant constraints on key elements like Zn.
In the end, our work, which fully complements stellar archaeology, represents a fresh start for the searches of first-star pollution in high-z environments that can provide unique insight on both the early phases of reionization and on the physical processes shaping the evolution of the first galaxies (Pallottini et al. 2014;Salvadori & Ferrara 2012).