© The Authors 2023

1. Introduction

Investigations of high redshift quasars significantly contribute to our knowledge of the physical conditions in the early Universe, including the growth of supermassive black holes (SMBH) and their co-evolution with galaxies, and their effect on the re-ionisation history of the Universe. Quasars also serve as bright beacons to illuminate the baryonic material along their lines of sight since their high luminosity allow their detection to very large distances. Finding high-redshift quasars has been the focus of numerous scientific investigations. From colour selection of sources in large imaging surveys (e.g. Fan et al. 2001; Willott et al. 2010; Wang et al. 2019; Matsuoka et al. 2022; Bañados et al. 2023), hundreds of quasars beyond redshift 6 have been discovered and compiled in databases (Inayoshi et al. 2020; Fan et al. 2023), while beyond redshift 7 only eight quasars have been discovered to date (Mortlock et al. 2011; Bañados et al. 2018; Matsuoka et al. 2019a,b; Yang et al. 2019, 2021; Wang et al. 2018, 2021).

Spectra of quasars allow detailed studies of intervening baryonic material from absorption lines arising either in the environment near quasars as well as along random lines of sight through the Universe, and are invaluable probes to analyse the intergalactic medium (IGM) and the epoch of re-ionisation (EoR). Due to the rapid decrease in quasar number densities compared to galaxies at the highest redshifts, quasars contribute only < 7% of the re-ionising photon budget at z ∼ 6 − 6.6 (Jiang et al. 2022). The neutral IGM appears as a damping wing at the red part of the Lyα absorption in quasar spectra (Miralda-Escudé 1998), and the shape of this damping wing in quasar spectra at different redshift can be used to compute the neutral gas fraction X_H I as the re-ionisation progresses (Bañados et al. 2018; Davies et al. 2018, 2020). A complete absorption by the IGM in the Gunn-Peterson (GP) absorption trough (Gunn & Peterson 1965) detected in z > 6 quasar spectra (Fan et al. 2001; Mortlock et al. 2011) reveals an increase in the Lyα optical depth with increasing redshift (Becker et al. 2021).

In the otherwise dark regions of the GP troughs, transmission spikes in the IGM start appearing at z ≲ 6 where ionised regions allow the transmission of Lyα photons at those redshifts. Such spikes therefore need to lie near sources that produce ionising radiation. Kakiichi et al. (2018) and Meyer et al. (2020) investigated the correlation between Lyα forest transmission spikes and the location of Lyα emitters and galaxies within a distance of 60 Mpc from the line of sight towards z > 6 quasars, and they argue that re-ionisation is primarily driven by a population of galaxies that are fainter than the detected galaxies.

In addition to probing the EoR through studies of hydrogen absorption, quasar spectra can be used to measure the evolution of heavy element enrichment with cosmic time from intervening absorption line systems to z ∼ 6 (see Becker et al. 2015, for a review of metal absorption lines). Quasars probe random lines of sight through the Universe, and intersect galaxies either in their circumgalactic medium (CGM) or interstellar medium (ISM) out to the quasar redshift, all of which imprint a spectral signature in the background quasar spectrum.

While the metallicity on the average increases with time, metal absorption lines also reflect the ionisation state and are therefore sensitive to the ultraviolet background radiation. Early analyses of the Lyα forest clouds and their correlations with C IVλλ1548,1550 absorption lines revealed that the IGM at z = 3 is enriched to a metallicity level of ≈10⁻³ Z_⊙ (Cowie et al. 1995; Ellison et al. 2000). While metals seem to be prevalent, absorption line densities¹, defined as the number of absorption line systems detected per unit (comoving) redshift path length interval, can be used to map the evolution of the metal enrichment and ionisation level. For example, C IV absorber line densities remain roughly constant with redshift out to z ∼ 4 and then show a steep drop at z > 5 (Becker et al. 2011; Codoreanu et al. 2018; D’Odorico et al. 2022; Davies et al. 2023a). A similar trend is found for the other high-ionisation line Si IV (Cooper et al. 2019; D’Odorico et al. 2022), suggesting a rapid physical change in the ionisation state at these redshifts. The change in ionisation state rather than a change in metal enrichment is supported from the observation of a change in the C IV/C IIλ1334 ratio (Simcoe et al. 2020; Davies et al. 2023a). Similarly, the line density of weak Mg IIλ2796 absorbers remain constant over a large redshift range of 2 < z < 7 (Chen et al. 2017), whereas the line density of O Iλ1302 increases over the range 3.2 < z < 6.5 (Becker et al. 2011, 2019) and is also interpreted as a change in the ionisation state at the EoR, where the Universe becomes progressively more neutral.

Bright quasars also allow detailed analyses of individual absorption systems by measuring column densities of various atomic species, similar to studies of 2 < z < 3 strong hydrogen absorption line systems (e.g. Prochaska et al. 2003). Focusing on the highest redshift absorbers at z > 6, this epoch approaches the time when the first population of massive stars produce the initial chemical enrichment. However, for intervening absorption lines systems the absence of measurable H I column densities in the quasars heavily absorbed Lyα forest or GP troughs makes it impossible to derive absolute metallicities, and only a few studies of relative metallicities at z > 6 have been performed to date. One example is the z ∼ 7.5 quasar ULAS J1342+0928 whose absorption system at z = 6.84 with measured iron, carbon and silicon suggests no alpha-element enhancement nor any special signatures of enrichment from very massive stars (Simcoe et al. 2020). When the absorption system lies close to the quasar redshift, in the so-called associated or proximate absorption systems, detailed metallicity measurements can be performed since the intrinsic H I column density can be modelled from fits to the H I damping wing (Bañados et al. 2019; Andika et al. 2022). However, proximate absorption systems do not trace a random population, but rather probe the high-density environments near the quasars (Ellison et al. 2010).

The connection between galaxies in emission and gas in absorption inform us of the transverse metal enrichment in the CGM of galaxies, albeit only along a single line of sight to the background quasar. Hence, combining the observations with simulations will reveal the spatial distribution of metals in the CGM of galaxies and its evolution with redshift (e.g. Finlator et al. 2018; Peeples et al. 2019). At z ≥ 5 overdensities of line emitting galaxies either selected via Lyα emission (Díaz et al. 2021) or [O III] λ5007 (Kashino et al. 2023) have been found within 200 kpc of the metal absorbers and approximately at the same redshifts.

Large amounts of ground-based telescope time have in recent years been devoted to spectroscopic studies of z > 6 quasars (e.g. D’Odorico et al. 2023). However, with the successful launch and commissioning of the James Webb Space Telescope (Gardner et al. 2023; Rigby et al. 2023), the question arises of what new insights observations of high redshift quasars obtained with JWST’s near infrared spectrograph, NIRSpec (Jakobsen et al. 2022), may bring to bear on the subject matter. While the highest spectral resolution possible with NIRSpec of R ≃ 2700 is only just adequate for intervening absorption line studies, the benefits offered by the instrument’s very high sensitivity and wide (0.81 μm ≤ λ ≤ 5.3 μm) infrared spectral coverage free of the atmospheric telluric absorption and OH sky emission lines that hamper ground based near-IR observations are self-evident.

In this paper, we present high signal-to-noise ratio R ≃ 2700 JWST/NIRSpec spectra of four well-studied quasars at 6.5 < z < 7.5. We focus our analysis exclusively on the intervening absorption line systems detected. The paper is organised as follows. Section 2 describes the observations. The spectra are presented and analysed in Sect. 3. The measurements of metal absorption line systems and the redshift evolution of line densities of atomic transitions in oxygen, carbon, silicon and magnesium, calcium and sodium are discussed in Sect. 4. In Sect. 5 we compute metal abundance patterns and compare with the enrichment of Population III and Pop II massive star Supernova explosions, and in Sect. 6 we compare the metal absorption system redshifts with transmission spikes seen in the quasar GP troughs. In Sect. 7 we present a discussion and conclude in Sect. 8. Throughout the paper we assume a flat cosmological model with H₀ = 67.4 km s⁻¹ Mpc⁻¹ and Ω_m = 0.315 (Planck Collaboration VI 2020).

2. NIRSpec observations

In this paper we present high signal-to-noise ratio (S/N ≃ 50 − 200), R ≃ 2700 NIRSpec spectra spanning 0.86 μm ≤ λ ≤ 5.3 μm of the four high redshift quasars: ULAS J1342+0928 (Bañados et al. 2018, z = 7.54, J_AB = 20.3), VDES J0020–3653 (Reed et al. 2019, z = 6.83, J_AB = 20.4), DELS J0411–0907 (a.k.a VHS J0411–0907, Pons et al. 2019; Wang et al. 2019, z = 6.81, J_AB = 20.4) and UHS J0439+1634 (Fan et al. 2019, z = 6.51, J_AB = 17.4). These data were obtained as part of the NIRSpec GTO programmes 1219 (PI: N. Lützgendorf) and 1222 (PI: C. Willott).

The NIRSpec Band I (0.86 μm ≤ λ ≤ 1.88 μm) and Band II (1.70 μm ≤ λ ≤ 3.15 μm) spectra were taken with the NIRSpec Fixed Slits together with the G140H/F070LP and G235H/F170LP grating and filter combinations, whereas the long wavelength Band III (2.88 μm ≤ λ ≤ 5.3 μm) portions were obtained employing the NIRSpec Integral Field Unit (IFU) and the G395H/F290LP grating and filter combination (see Jakobsen et al. 2022 and Böker et al. 2022 for details of these NIRSpec observing modes). The dates on which the observations were carried out and the position angles of the slits employed are summarised in Table 1.

The Band I and Band II Fixed Slit observations all employed 65 group NRSIRS2RAPID (Rauscher et al. 2017) full frame detector sub-integrations. For ULAS J1342+0928, VDES J0020–3653, and DELS J0411–0907, spectra were obtained in both the S200A1 and S200A2 fixed slits in order to span the detector gap, with two integrations carried out in each of three nodded positions along the slit for G140H/F070LP (11 554 s total exposure time) and one integration at each nod in G235H/F170LP (5778 s total exposure time). The brighter UHS J0439+1634 was only observed in the S200A1 slit at three nodded integrations in both G140H/F070LP and G235H/F170LP (2889 s total exposure time). These data were reduced and combined using the NIRSpec GTO Team pipeline (Carniani et al. in prep.). The resulting flux-calibrated spectra refer to Barycentric vacuum wavelengths sampled at Δλ = 2.36 Å per wavelength bin in Band I and Δλ = 3.96 Å per wavelength bin in Band II, which by design is close to the average native pixel sampling for these dispersers.

One noteworthy aspect of the shortest wavelength Band I G140H observations is that they were taken with the F070LP order-separation filter, rather than the default F100LP filter. This was done to include wavelengths in the Gunn-Peterson troughs below emitted Lyα down to 0.813 μm with the S200A1 slit and down to 0.858 μm with the S200A2 slit. This is permissible in this case since there is little or no light from the quasar emerging in the trough and therefore no second order light (beyond that of the sky background, which is subtracted out in the reduction) to contaminate the primary first order spectrum at wavelengths λ > 1.26 μm. However this choice did require tricking the pipeline into reducing the full Band I spectrum in two steps. First by running it for the F070LP filter giving the correctly reduced spectrum over 0.85 μm ≤ λ ≤ 1.26 μm and then running it a second time pretending that the F100LP filter was used and then scaling the resulting red end of the spectrum covering the 1.0 μm ≤ λ ≤ 1.88 μm region by the known ratio of the transmission of the F070LP and F100LP filters. The latter correction amounts to less than 4% in flux.

The Band III G395H/F290LP IFU observations of ULAS J1342+0928, VDES J0020–3653, and DELS J0411–0907 employed 25 group NRSIRS2 sub-integrations at 6 dithered IFU positions (11 029 s total integration time) and UHS J0439+1634 18 group NRSIRS2 integrations at 8 dithered positions (10 620 s total integration time). The IFU observations of VDES J0020–3653, and DELS J0411–0907 are presented in Marshall et al. (2023), who along with Perna et al. (2023) also provide the details of how the IFU data were reduced. Here we merely employ the one dimensional G395H/F290LP spectra extracted from the four IFU data cubes over a 1$\stackrel{″}{.}$0 diameter aperture centred on the quasars to supplement the corresponding Band I and II Fixed Slit spectra of primary interest. The resulting Band III spectra are sampled at Δλ = 6.65 Å per wavelength bin.

The final combined one-dimensional spectra of the quasars along with their 1σ error spectra are shown in Figs. 1–4. The signal-to-noise (S/N) ratios achieved in these four spectra are very high. This is illustrated for the representative case of VDES J0020–3653 in Fig. 5. The lower panel plots the Band I and II signal-to-noise ratio per wavelength bin calculated as the ratio of the observed signal divided by the one sigma error spectrum output by the GTO pipeline after median smoothing over a 40 bin window to remove the absorption lines. It is apparent that S/N ≃ 50 − 100 per wavelength bin is achieved at nearly all wavelengths. The overall variation with wavelength seen in this plot partly reflects the blaze functions of the two gratings employed and the variation in the continuum flux due to the broad emission lines in the quasar spectrum, but also the fact that not all wavelengths are exposed equally due to the gap between the two NIRSpec detector arrays. When employing the G140H grating with the S200A1 fixed slit the wavelength span 1.3021 − 1.3391 μm falls on the detector gap. For the S200A2 fixed slits the wavelength range 1.3479 − 1.3849 μm is lost to the gap. Similarly, for the G235H grating respectively 2.1825 − 2.2445 μm and 2.2594 − 2.3215 μm project to the gap for the S200A1 and S200A2 slits. Moreover, with the S200A1 slit, the G235M spectrum falls off the red edge of the detector for λ > 3.0491 μm but for λ > 3.1262 μm for the S200A2 slit. These unavoidable wavelength regions of reduced S/N only covered in one of the two slit spectra are indicated in Fig. 5.

The high S/N ratio observations of the bright (J_AB ≃ 20) quasars presented here represent a radically different and rarer regime of applicability for NIRSpec compared to the perhaps more familiar lower spectral resolution observations of the much fainter (J_AB ≃ 26 − 29) highest redshift galaxies that the instrument was primarily designed for (cf. Curtis-Lake et al. 2023; Bunker et al. 2023). It is therefore worth pointing out a few further subtle quirks of NIRSpec revealed by these data that are usually masked at lower S/N ratio.

While NIRSpec grating observations of faint high redshift galaxies are invariably severely detector noise limited (cf. Appendix in Jakobsen et al. 2022), the photon noise in our quasar spectra typically exceeds the total detector read noise and dark current noise by a factor ≃4. Furthermore, the quasar photon signals typically exceeds that of the (zodiacal light dominated) sky background experienced by JWST by a factor ≃20. However, one important feature that these quasar spectra still share with their fainter counterparts is that they are rather coarsely sampled both in wavelength and spatially along the slit. As explained in Jakobsen et al. (2022), with the aim of optimising the faint-end sensitivity limit of NIRSpec, the 194 mas angular width of the NIRSpec S200A1 and S200A2 Fixed Slits projects on to just under two mean 104 mas pixels on the detector. As a result, the apparent line spread function in a given spectrum at a given wavelength will depend on how precisely the slit width is sampled by the between two and three pixels that its monochromatic image happens to project to. This variable sampling issue is further exacerbated by the slit tilt and optical distortion that NIRSpec spectra are subject to, which necessitate that significant resampling of the raw 2D spectra be performed in the NIRSpec reduction pipeline in order to rectify the spectra to a common wavelength and spatial scale that will allow accurate background subtraction and summation of the net quasar signal along the spatial direction. That the final spectra are assembled by co-adding several differently sampled dithered sub-exposures taken at different positions along the slit leads to further smearing of the apparent line spread function. Lastly, the fact that the targets observed here are point sources, whose monochromatic images at the shorter wavelengths are not expected to fill the full width of the slit, adds even further variation in the spectral sampling. The net impact of these effects can be readily seen in Fig. 6 which shows close-ups of four representative narrow absorption lines detected in our final spectra, overlaid by their best fit Gaussian profiles. It is evident that the coarse sampling of NIRSpec spectra causes the apparent widths of the absorption lines to vary randomly at the ≃1 wavelength bin level, even among lines having comparable strengths. This stochastic variation in the apparent shape of the pixelated line spread function due to sampling noise clearly precludes performing meaningful detailed profile fits to our unresolved absorption lines. However, it ought not bias measurement of the main parameters of interest – line equivalent widths and centroid wavelengths – provided the error propagation calculation for these measured parameters correctly mirrors the noise level in the resampled spectra.

The top panel of Fig. 5 shows the anticipated 3σ observed equivalent width sensitivity limit to detecting unresolved absorption lines calculated from the lower S/N curve per the expression

$$\begin{array}{c}\hfill W_{\lambda }>3\frac{\mathrm{\Delta }\lambda }{(S/N)}{n}_{\mathit{FWHM}}^{3/2}\end{array}$$(1)

assuming the FWHM of the NIRSpec line spread function to be a multiple of the wavelength bin size δλ = n_FWHMΔλ with n_FWHM (the number of spectral pixels sampling the line spread function), conservatively set to n_FWHM = 2.3 for both bands. It follows that we anticipate being able to measure absorption lines having observed equivalent widths W_λ > 0.2 Å in our spectrum of VDES J0020–3653. The equivalent S/N and line sensitivities achieved for DELS J0411–0907 and ULAS J1342+0928 are comparable to those of VDES J0020–3653 as shown in Fig. 5. However, our spectrum of the considerably brighter UHS J0439+1634 reaches S/N ≃ 90 − 200, albeit with actual gaps in the spectrum at the missing S200A1 wavelengths.

We have not attempted to quantify the S/N achieved in the IFU-based Band III G395H extensions to our spectra. One reason for this is that the STScI pipeline used to reduce the IFU data outputs an estimate of the error spectrum that significantly underestimates the actual noise level in the data. This does not affect our results, since we only detect a single Na I absorption feature in the Band III G395H data (Sect. 4.4.7).

A final consequence of the coarse sampling and optical distortion that NIRSpec spectra are subject to is worth mentioning. The NIRSpec GTO pipeline employs a strict photon conserving (pixel projection weighting) resampling scheme when rectifying the native 2D spectra to a common wavelength and spatial scale. Nonetheless, in the case of a point source, the coarse sampling of the native 2D spectrum gives rise to Moire-style numerical fringing in the resampled spectrum along any resampled spatial row (cf. Smith et al. 2007). The amplitude of this unavoidable numerical fringing (or ‘wiggle’) at any given wavelength is effectively averaged out provided the net signal is integrated over a sufficient number of spatial rows large enough to capture most of the light from the target. In the spectra presented here the GTO pipeline performed this summation over a fixed five spatial rows. This is adequate to integrate down the fringing at wavelengths λ < 2.8 μm, but gradually fails at longer wavelengths and can be seen in our Band II G235H spectra due to the FWHM of the PSF increasing linearly with wavelength in the diffraction-limited NIRSpec, causing a varying amount of light to be missed in the five row summation. This numerical fringing effect is particularly pronounced in NIRSpec IFU data (Marshall et al. 2023; Perna et al. 2023), but is less obvious in our 1D Band III spectra since these were extracted by summing over a large number of spaxels. In any event, since the effect causes a modulation of the entire spectrum it should not affect measured equivalent widths and line centroids as long as the local continuum is set to follow the modulation (Sect. 3.2).

3. Analysis of the spectra

The full one-dimensional spectra of the quasars are presented in Figs. 1–4, with emission lines commonly seen in quasar spectra marked. We find a small shift in the absolute flux levels between the Fixed Slit and the IFU spectra, and have scaled the extracted one-dimensional Band III spectra by a factor of ≈0.7 to match the flux level at the overlapping wavelengths in the Fixed Slit Band II spectra.

3.1. Quasar systemic redshifts

In this paper, we focus on the analysis of intervening absorption line systems, and therefore re-visit the estimated quasar systemic redshifts in order to distinguish intervening from proximate metal absorbers arising less than 3000 km s⁻¹ from the quasar redshifts. This velocity cut-off is chosen since absorbers nearer than this limit show different physical properties such as metallicities and ionisation levels compared to genuine intervening absorption systems (Ellison et al. 2010). UHS J0439+1634 is a broad absorption line (BAL) quasar, and therefore could potentially affect its environment to even larger velocities of ∼10 000 km s⁻¹. In this spectrum we find two absorption systems that are offset by 9000 km s⁻¹, and both systems show narrow lines. We choose to include both absorption systems in the statistical analysis having a uniform cut at 3000 km s⁻¹ for the four quasar spectra.

Determining systemic redshifts for quasars is complex because of the dynamical effects that present themselves in quasar spectra. Broad emission lines can be blueshifted by several thousand km s⁻¹ relative to the systemic redshift, and can be composed of multiple components with peaks at distinct wavelengths, reflecting distinct kinematic components in the broad line regions. With the extended coverage of the rest-frame optical Balmer lines and narrow lines from [O II] λλ3727,3729 or [O III] λλ4959,5007 afforded by the NIRSpec spectra, more accurate systemic redshifts can be derived.

Fits to Hα lines from the Band III spectra of VDES J0020–3653 and DELS J0411–0907 and redshift measurements are presented in Marshall et al. (2023). The asymmetric [O III] line profile of VDES J0020–3653 is caused by excess emission from extended regions spatially offset by ∼1 arcsec from the quasar (Marshall et al. 2023). Our fit to the Hα emission line from VDES J0020–3653 is best reproduced by a combination of a broad Lorentzian and a broad Gaussian profile plus much weaker narrow components of Hα and [N II] λ6584 lines from the host galaxy, which gives a redshift of z = 6.8601 ± 0.0001, that is slightly higher than z = 6.855 based on a fit of two Gaussian components reported by Marshall et al. (2023).

Fits of Lorentz profiles to Hα and Hβ emission from DELS J0411–0907 suggest a systemic redshift of z = 6.8253 ± 0.0002 consistent with the initial redshift reported by Pons et al. (2019), and z = 6.8260 ± 0.0007 based on [C II] 158 μm emission (Yang et al. 2021), while multi Gaussian component fits to [O III] λ5007 and Hα emission lines suggest a smaller redshift of z = 6.818 as the systemic one (Marshall et al. 2023). The difference in the inferred systemic redshift is caused by a narrow component of Hα that is blueshifted by ∼380 km s⁻¹ relative to the broad Hα component (Marshall et al. 2023).

The [O III] λ5007 line in the spectrum of ULAS J1342+0928 can be fit by a double Gaussian line profile, where the narrow component suggests a lower redshift of z = 7.5317 ± 0.0007 compared to the centroid of the broad component at z = 7.5435 ± 0.0002. In comparison, a single Gauss component fit to the Hβ emission line provides an adequate fit with z = 7.5353 ± 0.0004. The redshift derived from Hβ is smaller by 211 ± 35 km s⁻¹ compared to z = 7.5413 ± 0.0007 measured from [C II] 158 μm emission from the quasar host (Venemans et al. 2017).

The Mg II line profile in the spectrum of UHS J0439+1634 is asymmetric with a significant blue wing. Fitting the broad line with a Lorentzian profile, we derive a redshift of z = 6.5102 ± 0.0005 consistent with a fit to the same line in a GNIRS spectrum from the Gemini telescope that reveals z = 6.511 ± 0.003 (Fan et al. 2019). The [C II] 158 μm emission from this object suggests a systemic redshift of z = 6.5188 ± 0.0004 (Yang et al. 2019). The Hα line is covered by the NIRSpec Band III data, and its line profile shows a broad component that is well fit by a single Lorentzian profile with narrow emission line components superimposed. The narrow emission lines are recognised as Hα and [N II] λλ6548,6584 from the host galaxy, suggesting a redshift of z = 6.5185 ± 0.0003, consistent with that from the [C II] 158 μm emission. The systemic quasar redshifts inferred from our spectra and used in the following analysis are summarised in Table 2.

3.2. Normalisation of spectra

For computations of absorption line equivalent widths, we need to determine the level of the continuum emission around the lines. First, we make a global fit of the continuum level across the entire spectrum using the Astrocook code (Cupani et al. 2022), using nodes with a velocity spacing of 1000 km s⁻¹, and reject nodes that deviate at the 5σ level within each velocity-window. Some of these outliers are caused by the presence of absorption lines. A spline function is then fit between the remaining node-points. A small window-size is necessary in order to also fit the broad emission lines of the quasars. The spectra in Figs. 1–4 show strong modulations in the continuum flux level in the fixed slit spectra around 1.8 μm and 3.0 μm as described in Sect. 2. We therefore choose to manually correct the continuum node points within Astrocook before the normalisation of the spectra, f_norm(λ), and their associated error spectra, σ(f_norm, λ) are computed. The normalised spectra are shown in Appendix A with the lines belonging to the identified absorption systems at different redshifts overlaid.

With normalised spectra the observed equivalent widths (W_obs) and their uncertainties are computed by

$$\begin{array}{cc}\hfill W_{\mathrm{obs}}& ={\displaystyle \sum _{{\lambda }_1}^{{\lambda }_2}}(1-f_{\mathrm{norm}}(\lambda ))\mathrm{\Delta }\lambda \hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \sigma (W_{\mathrm{obs}})& =\sqrt{{\displaystyle \sum _{{\lambda }_1}^{{\lambda }_2}}{\sigma }^2(f_{\mathrm{norm},\lambda })\mathrm{\Delta }{\lambda }^2}\hfill \end{array}$$(3)

where Δλ is the wavelength bin per pixel in the spectra.

In order to analyse the quasar spectra bluewards of their Lyα emission, we have to take a different approach to normalise the quasar spectra since the IGM has effectively absorbed most of the emission in the Gunn-Peterson trough at these redshifts. We use the python module PyQSOFit (Guo et al. 2018) to construct a model of the quasar continuum emission and broad Lyα emission line. We choose to model the continuum as a power-law function at wavelengths not affected by broad emission lines and superimpose best fits of the quasar Lyα emission line on top of the continuum, using a set of Gaussian functions with different widths, including both broad and narrow components. The observed spectra are divided by the constructed model spectra and the normalised spectra show the relative transmission in the IGM. We use these normalised spectra to explore the connection between the metal absorption line systems and transmission spikes in the IGM in Sect. 6.

4. Metal-line absorption systems

In the four quasar spectra we identify metal absorption line systems at several different intervening redshifts as well as in the vicinity of the quasar redshifts. Compared to the IGM where metal ions are typically highly ionised, intervening ISM or CGM lines are typically recognised via low-ionisation absorption lines from O I, Si II, C II, S II, Al II, Zn II, Cr II, Mn II, Fe II, Mg II, and Mg I, but also may have higher ionisation lines from Al III, C IV, and Si IV. We used the atomic line list in Krogager (2018), which is compiled from the Vienna Atomic Line Data Base (Piskunov et al. 1995; Ryabchikova et al. 2015).

The discovery and classification is done through visual inspection searching for systems of lines commonly detected in quasar damped Lyα systems (e.g. Prochaska et al. 2003) and Gamma-ray burst afterglow spectra that include strong damped Lyα systems (e.g. Christensen et al. 2011). Visual identifications, just as machine generated ones, are aided by the presence of line doublets, such as Mg IIλλ2796,2803 or C IVλλ1548,1550 whose wavelength separations and fixed relative strengths given by their oscillator strengths are well known. The relative line strengths are different by a factor of two in the optically thin regime, while in the optically thick regime, the line strengths are identical. In addition, absorption lines may have multiple kinematic components and blending with other lines makes the identification of line doublets more complex. In order to identify an absorption line system, we require that at least two or more lines are found for each system such that the identity of the lines and the redshift of the absorption system is uniquely determined.

In the spectrum of VDES J0020–3653 we identify 16 absorption systems, in DELS J0411–0907 we find 13 systems, in UHS J0439+1634 17 systems, and in ULAS J1342+0928 we find 15 individual absorption systems. In total, 61 identified systems are detected at 2.4 ≲ z ≲ 7.5. Appendix A presents the four normalised quasar spectra with identified absorption line systems overlaid. The systems and individual absorption lines are presented in both graphical and tabular form in Appendix B.

4.1. Distinguishing C IV doublets from O I/Si II absorbers

We need to pay particular attention to the identification of C IVλλ1548,1550 doublet and O Iλ1302-Si IIλ1304 absorbers since the two pairs of transitions have transition wavelength ratios of 0.99834 and 0.99831, respectively. These two separations cannot be distinguished at the R ≃ 2700 spectral resolution of the NIRSpec spectra, and weak systems may not be detected at the signal-to-noise level needed to unambiguously detect the C IVλλ1548,1550 doublet. To single out potential O Iλ1302-Si IIλ1304 absorption systems, absorption lines from Si IIλ1260 and Si IIλ1526 should be detected too, since they are stronger by a factor of 8.9 and 1.6, respectively, relative to the Si IIλ1304 line. The Si IIλ1260 line may, however, fall in the Lyα forest region and be completely absorbed in the GP troughs. In addition to other silicon transitions, we also search for other associated lines from C IIλ1334 or Mg IIλλ2796,2803. When none of these lines are detected, we classify the absorber as a C IV system.

4.2. Comparison with known absorption systems

As surveys detect more and more high-redshift quasars, an increasing number of studies of intervening absorption lines have been performed. The recent XQR-30 survey is based on a sample of 30 quasars at 5.8 < z < 6.6 observed with VLT/X-shooter at higher-resolution (D’Odorico et al. 2023). The only quasar in common with the NIRSpec sample is UHS J0439+1634, where Davies et al. (2023b) identify a total of 43 absorption line systems from z = 2.2 − 6.48. All of the 17 systems we detect, apart from the Mg II doublet at z = 6.208, are also detected in Davies et al. (2023b). In contrast, in our lower spectral resolution data we are not able to clearly identify and detect the other absorption line systems. In particular, Davies et al. (2023a) find several weak C IV doublets in the range from z = 5.0 − 5.3 and also weak rest-frame equivalent width (W_r < 0.05 Å) Mg II absorbers at z ≈ 2.3 − 5 that we cannot confirm in the NIRSpec data. The non-confirmation is either due to blending with other stronger absorption lines from different redshifts, or the fact that only one of the two lines in the doublet can be detected in the NIRSpec data, while the other component is clearly ruled out from the knowledge of the fixed line strengths between the two doublet lines. We are also not able to confirm several of the reported Si IV doublets. We note that around z ∼ 5 C IV falls in a wavelength region affected by (weaker) telluric absorption lines and also strong sky emission lines, that Davies et al. (2023b) do, however, correct for in their analysis.

Turning to absorption line systems at z > 5, Cooper et al. (2019) report two systems at z = 5.936 and 6.178 towards DELS J0411–0907, and at z = 5.889, 6.271 and 6.843 towards ULAS J1342+0928. These strong systems are all clearly detected in the NIRSpec spectra. In addition, we also detect a few weaker z > 5 Mg II or C IV doublets towards DELS J0411–0907, whereas towards ULAS J1342+0928, we also find several additional low-ionisation systems at z > 7 identified through multiple lines. ULAS J1342+0928 was observed with both X-shooter at the VLT and with FIRE at the Magellan Telescope at a higher spectral resolution as presented in a detailed case study of the z = 6.84 absorption system in Simcoe et al. (2020). We retrieved the processed 2-dimensional X-shooter spectra Advanced Data Products from the ESO archive in order to verify whether the high-redshift absorption line systems detected in the NIRSpec data are also present in the higher resolution spectra. For the higher redshift system at z = 7.368, the X-shooter data reveals a clear detection of Si IIλ1260 whereas Mg II falls in a noisy region of the spectrum and is not detected. None of the absorption lines from the z = 7.443 or z = 7.476 systems can be detected in the X-shooter data, consistent with the less sensitive detection-limit in the ground-based data.

4.3. Line densities

Theoretically we expect a detection limit of absorption lines observed equivalent widths at W_obs = 0.2 − 0.5 Å as described in Sect. 2 and Fig. 5. With the wide redshift range of the absorption systems, this implies that the detection limit on the rest-frame equivalent width W_r = W_obs/(1+z_abs) depends on the absorber redshift. The typical detection limit of W_r is 0.025–0.03 Å. When comparing line densities with those of from the literature, we adopt the same W_r cut-off.

To compute line densities, we combine all tables in Appendix B, extract the individual absorption transitions of interest, for example, C IVλ1548, O Iλ1302, or Mg IIλ2796, and compute their W_r. Line densities from metal-absorption lines should represent random lines of sight through the Universe, and we therefore exclude absorbers at redshifts closer than 3000 km s⁻¹ from the quasar systemic redshift, since these absorption lines represent proximate systems that may be influenced by the quasars.

We compute the redshift path length, which is defined as the full spectral range in each spectrum where it would be possible to detect the chosen absorption line. For each quasar we determine the minimum and maximum redshift where a chosen absorption line can be discovered, excluding regions in the spectra where there are gaps or where the signal-to-noise ratio is low. The sum of this gives the sensitivity function, g(z), which presents for each redshift the number of lines of sight probing the absorption line of interest with an observed W_obs > 0.3 Å. The result is presented in Fig. 7 which illustrates that the NIRSpec observations in this study are sensitive to intervening O Iλ1302 and C IIλ1334 absorbers at 6 ≲ z ≲ 7.5, Si IIλ1526 and C IV absorbers at 5 ≲ z ≲ 7.5, Mg II at 2.2 ≲ z ≲ 7.5, and Si IVλ1393 at 5.5 ≲ z ≲ 7.5.

By integrating the sensitivity function one can compute the line density dn/dz, defined as the number of absorbers n within a chosen redshift interval dz. However, when analysing the line density over a very large redshift range it is conventional to instead analyse the comoving density (Bahcall & Peebles 1969) and use the comoving absorption path length X defined as:

$$\begin{array}{c}\hfill \mathrm{d}X=\frac{{(1+z)}^2}{\sqrt{{\mathrm{\Omega }}_M{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}\mathrm{d}z.\end{array}$$(4)

The line density is then computed as

$$\begin{array}{c}\hfill \mathrm{d}n/\mathrm{d}X=\frac{{\sum }_{z_1}^{z_2}{n}_i}{{\int }_{z_1}^{z_2}g(z)\mathrm{d}X},\end{array}$$(5)

where n_i is the number of absorbers of the chosen species in quasar number i, summed over the redshift interval between z₁ and z₂ in all the quasar spectra. The denominator is the integral of the absorption path length in that redshift interval.

To derive uncertainties of the line densities we assume that the distribution of absorbers is Poissonian in nature. For the small number statistics relevant in this work we use the upper- and lower 68% confidence intervals tabulated in Gehrels (1986).

4.4. Absorption line density redshift evolution

Among the many different atomic species detected in the four quasar spectra, in this work we focus on the previously well-studied ones among the low-ionisation and high-ionisation lines.

For each atomic transition we compute the line density as a function of redshift and compare the JWST data with line densities at lower redshifts from the literature. Additionally, we compare the line densities with the predicted evolution of both low- and high-ionisations column density distribution from simulations in the redshift range 5 < z < 8 (Huscher et al., in prep.). The simulations are based on updates to the TECHNICOLOR DAWN cosmological hydrodynamic simulation (Finlator et al. 2018, 2020). By probing random lines of sight through the simulation volumes from redshifts z = 5 to z = 10, one can identify various absorption lines (e.g. Doughty & Finlator 2019, 2023) and compute their equivalent widths and column densities. Compared to previous simulations, the model of Huscher et al. includes predictions of equivalent widths and column densities of other atomic transitions for the modelled absorbers, allowing us to compare with the observations of a range of elements in different redshift bins from z = 5 to z = 10. For comparison with the NIRSpec observations, we restrict our attention to the relevant redshift interval 5 < z < 7.5.

4.4.1. O I λ1302 absorbers

We detect seven O Iλ1302 absorbers listed in Table 3 within the total absorber path length interval ΔX = 7.68. The highest redshift absorption system at z = 7.476 lies close to the ULAS J1342+0928 quasar redshift at z = 7.535, and is considered proximate according to the cut-off of at 3000 km s⁻¹, and is not included in the analysis of the statistics of line densities. The two other O I systems at z > 7 represent the highest redshift detections of neutral oxygen absorption along the line of sight to high-redshift quasars reported to date. In all cases, when O I is detected, we also detect more than one other line at the same redshift, so we can be confident to clearly identify the O I absorption systems (see Appendix B).

Figure 8 shows the line density of O I absorbers with W_r > 0.05 Å as a function of redshift for the new NIRSpec measurement compared with values at lower redshifts from Becker et al. (2019). The line density for the highest redshift point has a significant uncertainty, implying that even though the line density shows an increase from z = 6 to z = 7 the values are also consistent with being constant.

For O I as well as for other low-ionisation line species, the model predicts an increase in line density with increasing redshift up to z ∼ 7 and at higher redshifts (z > 7) the line density decreases. The reason for this behaviour is a combination of a decreasing metallicity of the CGM with increasing redshift, while the turnover from z = 7 → 6 is caused by the increasing level of ionisation of the CGM with time (Doughty & Finlator 2019). With increasing cosmic time, O I absorbers have a smaller covering fraction as the clouds’ outskirts become increasingly ionised while the total metallicity of oxygen and its distribution in the CGM of galaxies does increase with time (Doughty & Finlator 2019).

4.4.2. C II λ1334 absorbers

In the JWST spectra we identify eight intervening C IIλ1334 absorption systems listed in Table 4 over a full path length interval of ΔX = 11.10 at 5.9 < z < 7.5. Due to the available search path, some of the C II lines do not have detected O I absorption, since these lie bluewards of the quasar Lyα wavelength. The upper right hand panel in Fig. 8 illustrates the C II line densities with W_r > 0.03 Å.

In comparison, Becker et al. (2019) find that most of their O I absorbers also have C II lines, but this is primarily caused by selection since to identify O I absorption lines, the presence of other lines (e.g. from C II) at the same redshift is required. The line density of C II at z < 5 therefore shows a redshift evolution similar to that of O I. Using the absorption line catalogue of carbon absorbers in Davies et al. (2023b), we compute the C II line densities in three redshift bins at 5.0 < z < 6.5 as illustrated in Fig. 8. Including the line density at z > 6 from the NIRSpec data we find a general trend of an increasing line density with redshift.

4.4.3. Si II λ1526 absorbers

Several Si II transitions can be detected in the spectra with UV transitions at rest-frame wavelengths at 1260, 1304, 1526, and 1808 Å. Due to the wide range of oscillator strengths and possible contamination with other absorption lines, not all absorption lines are detected for any given system. When deriving the line density of Si II it is convenient to use a line which is both a strong transition and also has a large redshift path length. The 1808 Å line has the longest redshift path length, but it is the weakest transition of the four, whereas the strongest 1260 Å line can only be detected in a very short path length interval because its wavelength is close to Lyα. Therefore we choose to investigate the line density based on the 1526 Å line, which is also a strong transition, and has a larger redshift path length interval.

In the JWST spectra we identify 12 intervening Si II absorption systems listed in Table 5 over a full path length interval of ΔX = 29.0 at 5 < z < 7.45. We compare the inferred line density with that derived from the database in Davies et al. (2023b) in the lower left hand panel of Fig. 8. It is apparent that the evolution of the Si II line density could either be constant or show a mild increase with redshift to z ≈ 7.

In the models of Huscher et al. (in prep.), rest-frame equivalent widths are computed for the Si IIλ1260 transition. For non-saturated lines, the relation between two lines in the linear part of the curve of growth scale as $W_{r_1}/W_{r_2}={\lambda }_1^2{f}_1/({\lambda }_2^2{f}_2)$, where λ₁ and λ₂ are the rest-frame wavelength of the two transitions and f₁ and f₂ are the oscillator strengths. We can therefore convert the estimated Si IIλ1260 line densities to that for Si IIλ1526. The model predicts a constant line density of dn/dX ≈ 0.04 from z = 5 to z = 8 for absorbers with W_r < 0.025 Å, which is smaller by a factor of ∼12 compared to the NIRSpec observations. Alternatively, the model predicts Si II equivalent widths that are a factor of 2.5 smaller than the observed ones.

4.4.4. Mg II λ2796 absorbers

The classical Mg IIλλ2796,2803 doublet is one of the most commonly studied atomic transitions in quasar spectra (e.g. Bergeron & Boissé 1991; Churchill et al. 1999). It is easily distinguishable in absorption because of its known line separation and the fact that the 2796 line has twice the oscillator strength as the 2803 line. The line doublet enters the optical range at z = 0.3 and can be traced up to z ∼ 7 in ground-based near-IR spectra (Matejek & Simcoe 2012; Chen et al. 2017). In comparison to ground-based observations of Mg II absorbers at z > 2.5 that are heavily affected by the presence of sky emission- and telluric absorption lines, the NIRSpec spectra are much cleaner, allowing easy visual identification of Mg II absorbers even at the highest redshifts. It is therefore not surprising that we find some previously undetected absorption systems in these well-studied high-redshift quasars.

The 3σ detection limit of the W_obs for Mg II systems in the NIRSpec spectra is 0.3 Å. Above this limit, we find 49 Mg II absorbers from z = 2.4 to z = 7.45 including the highest redshift detection of an intervening absorption system to date (see Appendix B). As for the other lines, absorption line systems that lie closer than 3000 km s⁻¹ to the quasar redshift are omitted in the computation of line densities. The highest redshift proximate absorption system at z = 7.476 is also detected in Mg II, and is excluded in the line density computation.

The total redshift path length for Mg II absorption systems covered by the four spectra is ΔX = 73.3. To compare with other studies, we focus on the absorbers with W_r > 0.3 Å, where we detect 30 systems. The lower right hand panel in Fig. 8 illustrates that the line density for Mg II absorbers with W_r > 0.3 Å remains constant or at most exhibits a shallow decrease to the highest redshifts in agreement with previous findings at z ≲ 6 (Matejek & Simcoe 2012; Chen et al. 2017; Bosman et al. 2017). The observed distribution agrees with the model prediction of a shallow decrease in line density with redshift at z > 6.

Focusing on all the detected Mg II systems including both strong- and weak lines, we derive the rest-frame equivalent width distribution described by a powerlaw function

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}W}=\frac{n_{\ast }}{W\ast }exp(-W/W_{\ast })\end{array}$$(6)

where dW is the rest-frame equivalent width increment. We do not include a redshift dependence, since the line density is constant with redshift. The best fit parameters are W_* = 0.758 and n_* = 1.82. Figure 9 presents the distribution compared to a large quasar sample examined in Chen et al. (2017), and compared to ground-based observations, we see no evidence for a change in the W_r distribution in the NIRSpec spectra.

Strong Mg IIλ2796 absorption lines with rest-frame W_r > 1 Å have received particular attention since these absorbers are found to have the same redshift evolution as that derived from the integrated [O II] luminosity of galaxies and are proposed to trace the cosmic star-formation rate density at z ≈ 1 (Prochter et al. 2006; Ménard et al. 2011). It is therefore relevant to investigate the strongest Mg II lines at the highest redshifts with the new NIRSpec data. We identify 11 strong Mg IIλ2796 absorbers with W_r> 1 Å listed in Table 6. The line density dn/dX is illustrated in Fig. 10. Since we only have four JWST spectra, small-number Poisson statistics dominate the uncertainties (Gehrels 1986), but overall we see that the line density is consistent with other high-redshift quasar and gamma-ray burst afterglow samples that find a decrease in strong W_r > 1 Å systems at z > 3 (Chen et al. 2017; Christensen et al. 2017; Zou et al. 2021). None of the absorbers detected at z > 7 fall in the strong-line category.

4.4.5. C IV absorbers

Turning to the high-ionisation transitions, we identify 17 C IVλλ1548,1550 absorption systems at 5.20 < z < 6.86 in our spectra as listed in Table 7. The total path length interval for C IV is ΔX = 30.8. Two of the absorbers lie closer than 3000 km s⁻¹ from the quasar systemic redshifts, and are excluded in the computation of absorption line density shown in Fig. 11. The NIRSpec data reveal that the C IV line density continues to decrease with redshift beyond z > 6. Compared to the observed line densities, the model predicts a much lower line density. This discrepancy was also noted in D’Odorico et al. (2022), who discussed the mismatching C IV as a problem in the models in Finlator et al. (2020). If models are forced to match the line density of C IV, then the line densities for Si IV would be over-predicted.

4.4.6. Si IV absorbers

In the spectrum of ULAS J1342+0928 we detect a candidate Si IVλ1393 line at z = 6.476, which is also detected in C IV. The spectrum does not clearly reveal the fainter Si IVλ1402 line. In the spectrum of J0439+1634 we detect a single Si IV doublet, also reported by Davies et al. (2023b), but we note that both transitions in this doublet lie at the wings of two broad absorption line features in the quasar spectrum. These absorbers are listed in Table 8.

The total path length for Si IV systems in the JWST spectra is ΔX = 18.60. Figure 11 compares the Si IV line density from the highest redshift NIRSpec data to the lower redshift measurements of D’Odorico et al. (2022). As for the C IV systems, Si IV line densities are known to exhibit a rapid decrease in line density at z > 5 (D’Odorico et al. 2022), which we confirm with our very few detections of Si IV absorption lines.

4.4.7. Ca II and Na I lines

Most of the strong absorption lines arising in the ISM, CGM and IGM are from rest-frame UV lines, but the extended NIR wavelength coverage of the NIRSpec Band III G395H data allows us to search for rest-frame optical absorption lines associated with any of the absorption systems. Since the ionisation potential of Ca IIλλ3934,3969 and Na Iλλ5891,5897 is less than that of hydrogen, these absorption lines arise in neutral regions and are frequently associated with dusty absorbers at lower redshifts detected in quasar spectra (z ≲ 1; Heckman et al. 2000; Wild et al. 2006). Weak and strong Ca II absorbers, with a division at W_r = 0.7 Å, probe either galaxy disc gas or halo-components, respectively (Fang et al. 2023). We detect Ca II lines in three absorption systems at z = 3.17, 3.631, and 3.37, however, the signal is not sufficient to recover both lines in the doublet in any of the systems. Two of the absorbers fall into the weak line category with W_r, 3934 < 0.7 Å, while the z = 3.631 system towards ULAS J1342+0928 is a strong one.

The other rest-frame optical line of interest, the doublet Na Iλλ5891,5897 is covered by the Band III data for the highest redshift absorbers. However, as discussed in Sect. 2, the statistical errors on the Band III portions of the spectra are poorly defined. We therefore choose to only focus on absorption line systems that were identified in the Bands I and II Fixed Slit data, and test whether an associated Na I doublet is present in any of the z > 3.9 systems. Only a single high-redshift absorber at z = 5.6814 towards ULAS J1342+0928 has a potential Na I absorption line doublet detection.

4.5. Combined line density redshift evolution

To compare the evolution of the atomic species, all the line densities in previous subsections are combined in Fig. 12, which shows the evolution with redshifts for different elements derived in the literature with an extension to the highest redshifts (z > 6.0) detected in the NIRSpec data. All the computed line densities from the JWST data are listed in Table 9.

Both C IV and Si IV experience a drastic drop in line densities continuing at the highest redshifts (z > 6.5) detected in the NIRSpec data, while the line densities of the low-ionisation species O I and C II continue to increase at these redshifts. The evolution of Si II from low redshift (z < 5) is not constrained, and with the large uncertainties from the NIRSpec data, the Si II line density could be constant at the redshifts covered in this work, but could also trace the same increase in line density as detected for O I and C II. The transition in line densities of both carbon and silicon from predominantly low-ionisation at high redshifts to high-ionisation species at low-redshifts occurs at z ≈ 6.0. In contrast, the line density of Mg II remains constant throughout the range from z = 3 to z = 7, or may even experience a shallow decline towards higher redshifts as predicted by the model.

5. Metal abundances

5.1. Metal abundance ratios

The spectral resolution and sampling (Sect. 2) of our NIRSpec spectra is insufficient to accurately model the absorption line profiles via Voigt profile fitting techniques with the aim of deriving column densities for the different ions. However in simple cases, we can extract the same information from the measurements of the rest-frame equivalent widths. In the optically thin regime when the optical depth is τ_λ ≪ 1, the column density N scales linearly with the rest-frame equivalent width (Spitzer 1978):

$$\begin{array}{c}\hfill N=1.13\times {10}^{20}\frac{W_r(\mathrm{\AA })}{f{\lambda }^2(\mathrm{\AA })}[{\mathrm{cm}}^{-2}],\end{array}$$(7)

where f is the oscillator strength, and λ the rest-frame transition wavelength.

We derive the relative abundances assuming that the singly ionised atoms and neutral oxygen are the dominant state of ionisation in the absorption systems. This is a reasonable assumption, since none of the systems with O I absorption show any associated high-ionisation transitions. Moreover at the highest redshifts, the metallicities are expected to be low, so the optical depths of the transitions are low too. Finally, the column density ratios are compared to the solar photosphere abundances in Asplund et al. (2021). The silicon-to-oxygen and carbon-to-oxygen ratios are computed as:

$$\begin{array}{cc}\hfill [\mathrm{Si}/\mathrm{O}]& =log\left(\frac{N_{\mathrm{SiII}}}{N_{\mathrm{OI}}}\right)-log{\left(\frac{N_{\mathrm{Si}}}{N_{\mathrm{O}}}\right)}_{\odot }\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill [\mathrm{C}/\mathrm{O}]& =log\left(\frac{N_{\mathrm{CII}}}{N_{\mathrm{OI}}}\right)-log{\left(\frac{N_{\mathrm{C}}}{N_{\mathrm{O}}}\right)}_{\odot },\hfill \end{array}$$(9)

where the solar values are $log{\left(\frac{N_{\mathrm{Si}}}{N_{\mathrm{O}}}\right)}_{\odot }=-1.18\pm 0.05$ and $log{\left(\frac{N_{\mathrm{C}}}{N_{\mathrm{O}}}\right)}_{\odot }=-0.23\pm 0.06$ (Asplund et al. 2021).

5.2. Pop III enrichment signatures in [Si/O]

To date, no secure discovery of metal-free Pop-III stars has been made and the nature of the first generation of stars remains unknown. Pop-III stars leave behind a signature in the form of distinct metal-abundance patterns created by the stellar explosions of the massive stars. These patterns make their imprint on the chemical abundances of extremely metal-poor halo stars in the Milky Way (e.g. Frebel & Norris 2015), but also absorption line systems at high redshifts can be used to trace the nature of the first stars, since the time span for metal-enrichment of gas in the Universe is small.

In this section we focus on the six systems at z_abs > 6.5 with identified O Iλ1302 absorption lines in Table 3. Since no high-ionisation lines are detected in these absorption systems, they are predominantly neutral and we do not require ionisation corrections to derive metal abundance ratios. We examine whether the absorption lines are likely to be intrinsically saturated. Even though the spectra in Figs. B.1–B.61 suggest that none of the lines are saturated, at the low spectral resolution some lines could have hidden saturation. Indeed some of the lower-redshift strong absorption lines such as the strong Mg II lines are clearly saturated. For each of the absorption line systems we construct a conventional curve-of-growth with a best fit Doppler parameter between 20 and 40 km s⁻¹ from unblended Si II lines, where two or three separate transitions are detected. For the uncontaminated O I, C II, and Fe II lines, we verify that the transitions in W_r/λ vs. Nfλ fall on a straight line in the linear, non-saturated part of the curve of growth. This confirms that most lines are not saturated and that we can reliably derive the metal abundance ratios listed in Table 10.

Figure 13 shows the abundance ratios in a diagram of [Si/O] vs. [C/O] for the six z > 6 O I absorbers in our sample. Since C II in two of the six absorbers is partly contaminated by other absorption lines at other redshifts, and one is not detected, we provide only upper limits for three [C/O] ratios.

To compare the observed abundance pattern with the chemical enrichment from various types of stellar explosions, Fig. 13 illustrates the predicted ratios where the combined yields are weighted by the number of stars in a Salpeter initial mass function (IMF), reflecting the combined enrichment pattern generated by a stellar population. Pop III pair-instability SNe (PISN) with initial masses of 140 − 260 M_⊙ are illustrated by red coloured contours, and occupy a distinct region with a high [Si/O] and low [C/O] ratio (Heger & Woosley 2002). Pop III hypernovae (shown in green contours) and Pop III supernovae (blue contours) with a range of energies and initial masses of 10 − 100 M_⊙ also display distinct abundance patterns (Heger & Woosley 2010). Hypernovae are defined as having explosion energies of > 10⁵² erg, roughly ten times more than that of core-collapse SNe.

Several theoretical works have computed yields from metal-enriched Pop II. Using the yields table in Nomoto et al. (2013), stars with initial stellar masses between 10 − 40 M_⊙ and 5%–100% solar metallicity (illustrated in grey colour in Fig. 13) produce an abundance pattern that overlaps with that of hypernovae. The yields in Kobayashi et al. (2006) gives more or less identical results. A more recent yields table that includes Pop II stars with higher initial stellar masses between 13 − 120 M_⊙ (Limongi & Chieffi 2018) (purple contours) produce an IMF weighted yield with a higher concentration around [C/O] ≈ 0.1 and [Si/O] ≈ 0.3.

Variations of abundance ratios from Pop III and Pop II stellar yields have been used to infer the enrichment patterns predicted in Gamma-ray burst hosts (Ma et al. 2017), observed in low-metallicity absorption-line systems (Cooke et al. 2017; Welsh et al. 2019, 2023), and used to predict the yields in absorbers at the highest redshifts (Jeon et al. 2019). Stars in the local Universe also retain chemical signatures originating from the first populations of stars reflecting distinct abundance ratio patterns (Vanni et al. 2023a,b; Salvadori et al. 2023). See also Ma et al. (2017), Vanni et al. (in prep.) and D’Odorico et al. (in prep.) for abundance ratio diagrams other than [Si/O] and [C/O].

Figure 13 illustrates that also other intervening absorption-line systems at 4.5 < z < 6 (orange circles; Becker et al. 2019; Cooper et al. 2019) fall in the same region as the JWST absorption systems. However, several of these other absorption line systems at z < 5 have detected Si IV and C IV absorption lines that display different kinematical profiles and therefore arise in a different medium than the low-ionisation lines (Becker et al. 2019). In such cases a direct conversion from W_r of Si II and C II to a total Si and C column density is not valid. In fact, some of the absorption systems at z < 4.2 in Becker et al. (2019) have hydrogen column densities N(H I) < 10²⁰ cm⁻² for which ionisation corrections are necessary, and when these are included, all the absorbers at 3.0 < z < 4.2 have [Si/O] < 0.4 (Saccardi et al. 2023). For seven absorption systems at 5 < z < 6.2 their [Si/O] ratios without ionisation corrections remain close to the solar value (Becker et al. 2012), suggesting that Pop III enrichment does not dominate these systems.

In higher H I column density Damped-Lyα systems (DLAs), where the hydrogen column density is N(H I) > 2 × 10²⁰ cm⁻², the clouds are self shielded, such that the dominant contributions to the metal column densities are from the low-ionisation species. Above metallicities of 1% solar, C,N,O absorption lines are saturated in DLA systems and only lower limits on the column densities can be derived (e.g. Berg et al. 2016). The grey circles in Fig. 13 illustrate that metal-poor DLAs at z_abs ∼ 3 with metallicities from 10⁻³ to 10⁻² Z_⊙ in oxygen (Cooke et al. 2011, 2017; Welsh et al. 2019) generally have abundance patterns similar to solar values and consistent with Pop II SNe enrichment patterns, albeit with an overlap from the high-mass Pop III hypernovae enrichment pattern. However, when analysing other elements (Fe and Al), one of the DLA systems has an enrichment pattern which is better matched by Pop III SNe (Cooke et al. 2017; Welsh et al. 2023).

For the three absorption systems detected by NIRSpec with limits on the carbon abundance the enrichment could possibly be created by Pop III PISN. Two of the other absorbers show abundance ratios with high values of [Si/O] that are consistent with either enrichment by primarily lower mass (< 15 M_⊙) Pop III hypernovae or even Pop III core-collapse SNe for one of them. Models predict that the abundance ratio 0.5 < [Si/O] < 1.0 is produced by energetic Pop III SNe (with explosion energies > 10⁵¹ erg) from stars with initial stellar masses in the range from 10 − 15 M_⊙ (Heger & Woosley 2010). However, also single 15 M_⊙ star Pop II explosions produce a high [Si/O] ≈ 0.5 ratio (Woosley & Weaver 1995; Thielemann et al. 1996; Nomoto et al. 2006, 2013), so the enrichment by a single Pop II star remains a possibility. We do have to consider that gas clouds seen in absorption may not only be polluted by a single star alone, but by several SNe explosions over a range of stellar masses. When integrating the yields over the IMF, stars with masses in the range from 13 − 40 M_⊙ with initial population metallicities less than solar all produce an integrated ratio of [Si/O] ≲ 0.1 (Woosley & Weaver 1995; Kobayashi et al. 2006; Nomoto et al. 2013).

6. Transmission spikes in the Gunn-Peterson trough

One of the most debated questions arising from studies of the high-redshift Universe in the past decades, is what sources re-ionised the IGM. To further investigate this, we examine whether the sources associated with metal absorbers have ionised their surroundings.

We investigate the narrow transmission spikes at z > 5.7 in the Lyα forest region and in the GP troughs covered by the NIRSpec data using the normalised quasar spectra bluewards of the quasars Lyα lines described in Sect. 3.2. To locate narrow transmission spikes we search for peaks with a 3-σ prominence above zero, using the normalised spectrum divided by the normalised noise spectrum, namely a spectrum of the S/N ratio. We require that the peaks extend over three or more pixels to avoid detections from a single noisy pixel. The identified spikes have relative transmissions of 2–5% or higher. Figure 14 illustrates the location of the Lyα wavelengths for the identified metal-absorption line system redshifts plotted on the normalised spectrum in the GP trough region of the four spectra. The blue and red coloured dashed lines mark high- and low-ionisation absorption systems, respectively. The green triangles mark the spikes in the IGM where there is a significant transmission of Lyα photons, and Table 11 lists the redshifts of the encountered Lyα transmission spikes.

We briefly compare the detected IGM transmission spikes with previously published results. Yang et al. (2020) examine transmission spikes at z > 5.5 based on ground-based spectra of a large sample of high-redshift quasars including two of the quasars in this paper (DELS J0411–0907 and VDES J0020–3653). They list four spikes in DELS J0411–0907, with three being multiple peaks at an almost identical redshift of z ∼ 5.7. We also detect these four spikes in the NIRSpec data but in addition we also detect spikes at even higher redshift. For VDES J0020–3653, Yang et al. (2020) do not report spike redshifts. Due to the higher S/N ratio of the JWST spectra and the absence of noise due to telluric features, we can detect weaker spikes in the GP troughs.

Only a single low-ionisation absorption system (at z = 5.928 towards UHS J0439+1634) out of the total of 13 intervening low-ionisation absorption systems (red dashed lines in Fig. 14) match the transmission spikes in the IGM to within a velocity offset of 2000 km s⁻¹, whereas four out of nine of the high-ionisation system lie within ≈2000 km s⁻¹ of transmission spikes, all of which at z ∼ 5.8. Not surprisingly, the two proximate absorption systems in the spectra of UHS J0439+1634 and VDES J0020–3653 are indeed high-ionisation systems. Conversely, from the total of 21 absorption systems at z > 5.7, 16 spikes do not match any metal-absorption system within ±2000 km s⁻¹. The fact that most of the low-ionisation systems are not associated with a transmission spike supports the hypothesis that the metal-enrichment arises from star-formation followed by heavy element productions by supernova explosions preferentially in higher density regions that may be the last regions to be re-ionised (Oh 2002). On the other hand, half of the high-ionisation systems at 5.7 < z < 6.8 appear in regions of the IGM that are ionised to a degree that allows the escape of Lyα photons.

7. Discussion

7.1. A forest of lines or discrete absorption line systems?

Oh (2002) has suggested that as the Universe becomes progressively more neutral, a forest of weak metal lines may appear at redshifts close to the EoR. Analogous to the H I Lyα forest arising from neutral clouds of hydrogen in the IGM that becomes more numerous towards higher redshifts, metals in the IGM would also progressively become more neutral. For example an O Iλ1302 forest could give rise to a depression of the quasar continuum at the 10–20% level (Oh 2002). Other atoms such as C II, or Si II have ionisation potentials close to that of H I and have oscillator strengths that also make the lines useful as tracers of the enrichment and the increasing neutral IGM towards the re-ionisation epoch. Likewise, the Mg IIλλ2796,2803 doublet is easily detected, and could give rise to a continuous Mg II forest provided the neutral IGM is globally enriched with metals (Hennawi et al. 2021). However, ground-based spectra of high redshift quasars have not yet revealed any signature from a Mg II forest (Tie et al. 2023). In the NIRSpec data we detect a line density of Mg II that remains constant out to z = 7.5. Mg II lines are identified as discrete absorption systems rather than a global decrease in the continuum in the form of a Mg II forest. Also in the case of O I, we are only able to detect discrete absorbers that likely arise in the CGM or ISM of intervening galaxies. Further analysis of JWST/NIRSpec data using more advanced methods such as auto-correlation techniques may be able to find a metal forest, but that is beyond the scope of the analysis of individual absorption systems in this work. Detecting a Mg II metal-line forest with a global decrease in the continuum flux at the few percent level is difficult due to the wiggles in the NIRSpec data described in Sect. 2, but also the broad emission-line composite from Mg II and Fe II in the quasar spectrum will make it difficult to accurately determine the level of the unabsorbed continuum emission.

7.2. High-redshift metallicities

In the absence of measurable H I column densities, we can only measure relative metallicities. In the redshift interval 2 < z < 6 absorption systems display a relatively constant relative abundance pattern in [Si/O] and [C/Fe] with redshift and do not reveal clear signatures of Pop III SNe enrichment (Becker et al. 2012; Kulkarni et al. 2014), with a couple of possible exceptions from DLA systems at z ≈ 3 (Welsh et al. 2023). A relevant question is therefore if we have overestimated the [Si/O] or [C/O] ratios for the z > 6 absorption systems. In comparison with ground-based spectroscopic studies, the NIRSpec data have low spectral resolution and could potentially have hidden saturation, or be composed of multiple components rather than a single one assumed for the O I systems in this work. We argue against the possibility of hidden saturation based on the analysis of the curve of growth for the absorption lines in any of the O I absorbers where we find that they are consistent with being on the linear, optically thin region of the curve of growth. In agreement with our finding, numerical simulations at z ∼ 7 also find that O I absorption lines along random sight-lines are optically thin (Doughty & Finlator 2019).

Alternatively, multiple absorption components in the observed spectra could have variations in abundance patterns between them, but in order to display a high [Si/O] ratio, some of the components will inevitably have to have an extreme value. Ultimately, we need higher spectral resolution measurements to reveal and decompose the absorption lines in distinct components but this may have to wait until the arrival of high-resolution spectrographs on future extremely large telescopes (e.g. Maiolino et al. 2013).

Another possibility leading to inaccurate abundance ratios arises if the atoms exist predominantly in ionisation states other than Si II or C II in which case ionisation corrections are needed to derive [Si/O] and [C/O] ratios accurately. However, we do not regard this as likely for following reasons. Firstly, the intervening absorption systems are primarily neutral because there are no high-ionisation C IV or Si IV lines associated with the O I systems. Secondly, the O I absorbers are not associated with a transmission spike in the GP troughs at these redshifts, and consequently the environment is predominantly neutral. O I is charge-exchanged locked with H I, so O I is the dominant state of oxygen. Possibly, Si III and C III could be present, but these atomic states do not have absorption lines in the spectral range covered by the NIRSpec data. However, if that were the case, then the total [Si/O] and [C/O] would be even higher than detected in the high-redshift absorbers. For similar z ∼ 6 absorption systems, Cooper et al. (2019) argue that the low-ionisation systems are equivalents to z ∼ 3 DLA systems in which metallicities derived without ionisation corrections are reliable.

7.3. Ionisation state of the CGM

Our finding that the number density of the high-ionisation ions C IV and Si IV continues to drop at z > 6.5 extends the trends observed at slightly lower redshifts (Cooper et al. 2019; D’Odorico et al. 2022; Davies et al. 2023a). The disappearing of these absorption lines cannot be attributed to a change in metallicity, since C II and Si II are frequently detected with column densities that are ten times higher than in the high-ionisation states (Becker et al. 2011; Bosman et al. 2017; Cooper et al. 2019). Moreover, we find that the Mg II line density appears to be constant throughout the redshift range probed out to z ∼ 7.5, or at least only display a modest decrease to the highest redshifts covered, z ∼ 7.5. This supports the contention that the change in the high-ionisation line densities is caused by a softening of the UV radiation field dominated by star-forming galaxies rather than by the harder radiation of active galactic nuclei at these redshifts.

The under-production of C IV at z > 6 is known to be a problem in cosmological simulations (Keating et al. 2016; Finlator et al. 2018; Doughty et al. 2018) although interestingly we find that the number density of C II absorbers agrees with observations. We detect only two Si IV absorbers in the JWST spectra, and confirm that the line density of this transition continues to decrease with redshift consistent with model predictions. The fact that the observed line density of Si IV at z > 5.5 matches the simulations (D’Odorico et al. 2022), whereas the observed Si II line density is higher than in simulations, is not likely due to the hardness of the UV background in the simulations, because the ratio of these lines does not depend on the spectral hardness (Doughty et al. 2018). In the epoch where re-ionisation has not yet been completed, the low-ionisation ions will be the dominant ones, and may be present out to large distances from the galaxies (Doughty & Finlator 2023).

7.4. Environments of the metal absorbers

To date, not many studies have attempted to match metal-absorption lines to galaxies at z > 6 primarily due to the lack of instrument sensitivity since the galaxies are very faint. From ground-based observations, targeting Lyα emitters has recently been the best option due to the increased contrast of the emission lines above the background. By studying correlations of high-ionisation C IV systems at z ∼ 5 with Lyα emitters Díaz et al. (2021) found emitters within 200 kpc in projection from the absorbers. The environments of the even higher redshift low-ionisation systems are less clear from observations. Using JWST/NIRcam slitless spectra of one z = 6.3 quasar field Kashino et al. (2023) identify galaxies with [O III] emission lines at 5.3 < z < 6.2 and find that most of the detected galaxy redshifts match either high- or low-ionisations absorption-line systems in the quasar spectrum. Conversely, out of the eight intervening absorption line systems, only four of the systems have detected galaxies within 200 kpc of the quasar sight line (Kashino et al. 2023), where most of the non-detections are from low-ionisation absorber systems.

In a survey of eight quasars at z > 6, Meyer et al. (2020) detect a weak correlation between transmission spikes and z ≈ 6 Lyα emitters and Lyman-break galaxies (LBGs) within 1–10 proper Mpc in projection along the line of sight, but argue that fainter galaxies below the detection threshold are necessary to drive the re-ionisation. In a similar analysis of a single quasar at z = 6.4 Kakiichi et al. (2018) find that three out of six LBGs lie close in redshift space to the IGM transmission spikes, and also conclude that fainter galaxies clustered around the brighter LBGs are needed to drive re-ionisation, but also argue that luminous quasars do contribute to the ionising flux at the latest stages of the re-ionisation. In a recent analysis of JWST NIRcam slitless data, Kashino et al. (2023) measure a clear correlation between galaxy redshifts at 5.7 < z < 6.1 and transmission spikes in the Lyα forest, but find that metal-absorption systems are not connected to the transmission spikes. This is consistent with our finding of no clear association of low-ionisation absorbers with IGM transmission peaks.

Whether any of the absorption systems have detectable galaxy counterparts nearby can be investigated with the Band III JWST integral field spectra. The 3 arcsec field of view of the IFU limits the search to within 7–10 kpc in projection from the quasar, which is much smaller than the typical size of 100 kpc of the metal-enriched CGM around high-redshift galaxies (Tumlinson et al. 2017). The sizes and covering fractions of the metal-enriched CGM depend on the galaxy masses and redshifts, where simulations predict that the median separation of an O I absorber at z = 7 to its host galaxy is about 250 h⁻¹ comoving kpc or a physical size of about 45 kpc, and that the absorber may probe either the CGM or the ISM of the host (Doughty & Finlator 2019). Even though the IFU of NIRSpec has a small field of view, there is a non-negligible chance of detecting the associated galaxy through emission lines at the absorber redshift in the IFS data cube, especially if the strongest absorbers probe the ISM. This association of metal absorption lines with galaxies will be investigated in forthcoming papers.

8. Summary and conclusions

We have presented JWST/NIRSpec spectra of four z > 6.5 quasars covering the wide spectral range 0.8–5.3 μm. Targeting luminous quasars with magnitudes of J_AB = 17–20, the data shows an exquisite signal with a S/N ratio of 50–200 per spectral bin allowing for an in-depth analysis of each of the quasar lines of sight. We detect metal-absorption line systems from intervening absorbers and systems associated with the quasars. In each of the four spectra we identify between 13 and 17 absorption line systems with a total of 61 systems at redshifts ranging from z = 2.4 to z = 7.5. Since the line of sight towards the quasars probe random locations, the identifications include detections of intervening metal-enriched gas in the CGM or the ISM of intervening galaxies at the highest redshift to date. Most absorption line systems are identified in multiple atomic transitions due to the high quality of the NIRSpec data, while some weaker systems are only detected through Mg II or C IV doublets.

At the quasar redshifts, the Lyα transitions of each absorption system fall in the Gunn-Peterson trough, so deriving hydrogen column densities or metallicities based on the JWST data is impossible, and likewise is it impossible to derive the neutral hydrogen gas fractions. Instead, we use low-ionisation lines such as O I as a proxy to examine the changing ionisation state of the IGM (Oh 2002). Rather than detecting a forest of metal lines in the IGM at the end of the EoR, we detect instead a number of discrete O I (and Mg II) absorbers, as logically expected since winds from galaxies are responsible for distributing the heavy elements out into the CGM of the host galaxies.

From the number statistics of the low-ionisation absorbers O Iλ1302, C IIλ1334, Si IIλ1526, and Mg IIλ2796, we derive the redshift evolution of the line densities beyond z > 6. Combined with other studies of lower redshift quasar spectra (Becker et al. 2019; Chen et al. 2017; Davies et al. 2023b), we find that the line densities of O I, C II, and Si II derived from the NIRSpec data reveal a continued increase beyond z > 6, whereas Mg II shows a redshift evolution that is consistent with being constant or mildly decreasing with increasing redshift. In all atomic species, the number of detected absorbers remains small and the redshift path-length interval is also limited in our NIRSpec data set of just four quasars. Consequently, the Poisson uncertainties are large and the line density evolution could also be constant to within ∼1σ uncertainties for all atomic species. Nevertheless, compared to numerical simulations, the observations follow the predicted redshift dependence for the low-ionisation lines at z = 5 → 8, where the model predicts an increase with increasing redshift to z ∼ 7 due to the increased neutral content and a decreasing UV ionising background flux. At z > 7 the model predicts that line densities decrease with redshift as a consequence of the decreasing metal enrichment of the CGM.

For the high-ionisation transitions C IV and Si IV we derive a decrease in line density beyond z > 6 following the previously reported rapid decrease seen for these lines from z = 4 → 6 (D’Odorico et al. 2022; Davies et al. 2023a). Again, the drop is consistent with the model prediction. The transitions in line densities from a predominantly low-ionisation level to a high-ionisation level occurs at z ≈ 6.0 for both carbon and silicon. The observations therefore suggest that there is a continuous change in the ionisation state rather than a drastic change in the distribution of metals at z > 6.

We matched the redshifts of the absorption line systems at z > 5.7 with spikes of transmission of Lyα photons in the GP troughs, and find that half of the intervening absorbers that have high-ionisation absorption lines fall within ±2000 km s⁻¹ of the Lyα redshift of a matching transmission spike. For low-ionisation absorbers on the other hand, only a single system out of 13 intervening absorbers falls within ±2000 km s⁻¹ of a spike. This demonstrates that low-ionisation systems are found in environments that are not ionised sufficiently to allow Lyα photons to escape. Low-ionisation metal-absorption line systems therefore likely trace higher density environments in the Universe that are among the last ones to become re-ionised.

While we are not able to derive absolute metallicities, our relative metallicity ratios are accurate provided that the absorption lines are not saturated, and that ionisation corrections are insignificant. We argue that the systems traced through O I absorption lines are reliable tracers, since in none of the cases do we detect high-ionisation lines and they are not associated with IGM transmission spikes. Focusing on silicon, carbon and oxygen atoms, we derive high values of the metal line ratio [Si/O] in the range from 0.5–1.0 in five of the z > 6 absorption systems. Such high values of [Si/O] are consistent with the enrichment pattern expected from energetic Pop III hypernovae with initial stellar masses from 10 − 15 M_⊙, or even from Pop III core-collapse SNe (Heger & Woosley 2010). However, explosions of single Pop II stars with a stellar mass of 15 M_⊙ can also produce a high [Si/O] ratio (Woosley & Weaver 1995; Kobayashi et al. 2006; Nomoto et al. 2013), so this scenario cannot be ruled out if the absorption system happen to probe a gas cloud enriched by a single star alone. Absorbers with the highest [Si/O] ratios and upper limits on [C/O] could be enriched by Pop III pair instability supernovae. The high [Si/O] ratio and limit on the carbon abundance cannot be explained by the yields produced by Pop II SNe. We therefore suggest that the CGM of the neutral absorption systems at z > 6 are dominated by the enrichment from Pop III explosions.

In the near future, surveys from LSST and Euclid will detect an increasing number of (fainter) quasars at z > 7. Also with the increasing depth and completed wide-field surveys with JWST, new detections of quasars at z > 10 are being made (Goulding et al. 2023; Maiolino et al. 2023). Follow-up near-IR spectroscopic observations with NIRSpec and future ELT spectrographs will allow to extend this work to even higher redshifts into the first half of the re-ionisation epoch beyond z > 7.5. Higher resolution spectra compared to the R ≃ 2700 NIRSpec observations are needed to fit the absorption line profiles with conventional Voigt fitting techniques to verify metal column densities, and with an increasing number of high redshift quasar spectra the larger path length interval will provide line densities with smaller Poisson uncertainties. Future observations of the fields surrounding the four quasars presented in this paper may be able to identify the galaxies responsible for the intervening absorption line systems and single out the galaxies that are possibly responsible for the Pop III metal-enrichment at z > 6.

---

Acknowledgments

We thank Kristian Finlator for useful discussions, and Ezra Huscher for sharing results from their simulations. We also thank Stefania Salvadori for valuable insights into Pop III enrichment. L.C. is supported by DFF/Independent Research Fund Denmark, grant-ID 2032–00071. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant no. 140. AJB has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Advanced Grant 789056 ‘First Galaxies’. MP acknowledges support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI), and the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. The work is based on data obtained by the European Southern Observatory, Paranal, Chile, Program IDs: 098.B-0537, 0100.A-0446 and 0100.A-0898. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with GTO programs 1219 and 1222.

References

Appendix A: Normalised spectra and absorption line systems

For visualisation purposes the normalised Fixed Slit spectra of the quasars are shown here, with all the different intervening and proximate absorption line systems overlaid. Absorption lines are colour coded to mark lines belonging to the same redshift absorbers. Figures with zoom in on individually detected absorption lines for each system are presented in Appendix B. We note that the Lyα forest is ignored, and we do not include the Band III data from 3-5.3 μm, because only a very few absorption lines are detected at these wavelengths.

Appendix B: Individual absorption systems

The appendix presents each of the detected absorption line systems. In each of the quasar spectra we have identified between 13–17 intervening absorption line systems at different redshifts, ranging from lower redshifts (z = 2.4) to absorbers located at the quasar redshifts.

In the figures in this appendix all panels illustrate the normalised quasar spectrum around the detected absorption line species within ±600 km s⁻¹ for each of the absorption systems. Flux errors per spectral pixel are overlaid, and the y-axis of each absorption line is scaled to match the depth of the line. Equivalent widths in the observed frame are computed by integrating the lines from within ±200 km s⁻¹, while for the stronger systems that display broader lines, the lines are integrated within ±300 km s⁻¹. The reported wavelengths of the absorption lines in all tables are computed from weighted mean flux values rather than based on absorption line fits. Line profiles for narrow line systems for a single absorption system may appear visually different as explained in Fig. 6 in Section 2. Some lines are blended or partially contaminated with other identified lines from other absorption systems. We describe in the comments to the tables which lines are blended with a line from another system. Some weak lines or lines towards longer wavelengths are dominated by errors from the complex normalisation of the quasar spectrum.

All Tables

All Figures