Metal Lines Associated with the Lyman-alpha Forest from eBOSS Data

We investigate the metal species associated with the Ly$\alpha$ forest in the eBOSS quasar spectra. Metal absorption lines are revealed in the stacked spectra from cross-correlating the selected Ly$\alpha$ absorbers in the forest and the flux fluctuation field. Up to 13 metal species are identified associated with relatively strong Ly$\alpha$ absorbers (those with flux fluctuation $-1.0<\delta_{\rm Ly\alpha}<-0.6$ and with neutral hydrogen column density of ~ $10^{15-16}$ $\rm cm^{-2}$) over absorber redshift range of $2<z_{\rm abs}<4$. The column densities of these species decrease toward higher redshift and for weaker Ly$\alpha$ absorbers. From modelling the column densities of various species, we find that the column density pattern suggests contributions from multiple gas components both in the circumgalactic medium (CGM) and in the intergalactic medium (IGM). While the low-ionization species (e.g., C II, Si II, and Mg II) can be explained by high-density, cool gas ($T \sim 10^4$ K) from the CGM, the high-ionization species may reside in low-density or high-temperature gas in the IGM. The measurements provide inputs to model metal contamination in the Ly$\alpha$ forest baryon acoustic oscillations measurement. Comparison with metal absorptions in high-resolution quasar spectra and with hydrodynamic galaxy formation simulations can further elucidate the physical conditions of these Ly$\alpha$ absorbers.

As a tracer of the underlying matter density field, the Lyα forest is a powerful cosmological probe that has been successfully used to constrain cosmological parameters. In particular, the Lyα forest data from the spectroscopic observation of quasars in the the Sloan Digital Sky Survey (SDSS) has been used to measure the Baryon Acoustic Oscillations (BAO), adding unique high-redshift (z ∼ 3) data points to constrain cosmology (e.g., Slosar et al. 2013;Bautista et al. 2017;du Mas des Bourboux et al. 2020).
The Lyα forest also encodes valuable information about the circumgalactic and intergalactic media (CGM and IGM), including their density and temperature distribution. While the IGM and CGM are the source of material for star formation and galaxy formation, their properties are also affected by gas outflow from galaxies as a result of feedback processes. The outflow can bring heavy elements into the CGM and IGM. The metal content is connected to the stellar population, abundance pattern, and star formation history. The metal species at various ionization levels can constrain the ionizing ultraviolet (UV) background. The Lyα forest observation provides a means to reveal metals. As metal absorption in the Lyα forest is an important systematic component to be modelled in Lyα forest BAO measurement (Bautista et al. 2015), studying metals in the Lyα forest also aids the cosmological application of the Lyα forest.
High-resolution spectroscopic observations of quasars with long exposure times are usually the choice to study metal lines in the Lyα forest or associated with the Lyα forest (e.g., Ellison et al. 2000;Schaye et al. 2003;Simcoe et al. 2004;Lehner et al. 2019Lehner et al. , 2021. For example, Lehner et al. (2021) study 2.2 z 3.6 H I-selected absorbers with neutral hydrogen column density in the range of 14.6 ≤ log(N HI /cm −2 ) ≤ 20. The metallicity is found to have a broad distribution (e.g., an interquartile range from −3.6 to −1.8 dex with a median of −2.4 dex for log N HI = 14.5-16.2 absorbers), and the median metallicity increases and the interquartile range decreases with increasing neutral hydrogen column density.
High-resolution spectroscopic observations of quasars usually probe only the Lyα forest along the sight line of a single quasar or those of a small number of quasars. In contrast, the large spectroscopic survey of quasars from the various phases of the SDSS survey is able to probe hundreds of thousands of quasar sight lines with moderate resolution and short exposure times. Although the signal-to-noise ratio and the spectral resolution from the SDSS observation typically do not allow the identification of metal absorption lines associated with the Lyα forest in any individual quasar spectrum, the large number of quasar spectra can be used to perform a statistical study of metal lines. Pieri et al. (2010) and Pieri et al. (2014) develop a technique to study metal absorption lines using composite spectrum of Lyα forest absorbers identified in SDSS spectra of quasars. Lyα absorbers in the Lyα forest of each quasar are selected based on the transmitted flux fraction. For each selected absorber, the quasar's spectrum is shifted to the restframe of the absorber. The shifted quasar spectra for all the selected absorbers are then stacked to produce a composite spectrum. Features uncorrelated with the absorbers are highly suppressed in the composite spectrum, and the coherently stacked metal absorption lines associated with the selected absorbers can be detected with high significance.
Applying the technique to the data from the SDSS (Pieri et al. 2010) and from the Baryon Oscillation Spectroscopic Survey (BOSS) of SDSS III (Pieri et al. 2014) has revealed a series of metal absorption lines in the composite spectrum to probe the physical conditions and metal enrichment of the CGM and IGM at z ∼ 2-3.
In this paper, following the same basic idea in Pieri et al. (2010), we study the metal lines associated with the Lyα forest 1 by applying a cross-correlation method to obtain the stacked spectra using data from the extended BOSS (eBOSS; Dawson et al. 2016) survey of the SDSS IV. We study the dependence of the metal absorption lines on the Lyα fluctuation field and their redshift evolution over 2 < z < 4.
In Section 2, we describe the data and the method to compute the stacked spectra of Lyα absorbers. We then present the main results in Section 3, including the column densities of 13 metal species as a function of Lyα absorber strength and redshift. A simple model is discussed to understand the results. Finally, we summarize and discuss our findings in Section 4. We leave some details on column density measurements in Appendix A.

DATA SAMPLES AND REDUCTION METHODS
The Lyα forest data in this work are based on the quasar observation from the Sloan Digital Sky Survey (SDSS; York et al. 2000), gathered during SDSS-III by the Baryon Oscillation Spectroscopic Survey (BOSS; Eisenstein et al. 2011;Dawson et al. 2013), and during SDSS-IV by the extended BOSS (eBOSS; Dawson et al. 2016;Blanton et al. 2017), with a small fraction observed during SDSS-I and II (Schneider et al. 2010).
All quasar spectra used for the analysis here were obtained using the eBOSS spectrographs mounted on the 2.5 m Sloan Foundation telescope (Gunn et al. 2006;Smee et al. 2013;Myers et al. 2015) at the Apache Point Observatory, publicly available in the sixteenth data release (DR16; Ahumada et al. 2020). The spectral resolution varies with wavelength (3600-10400Å), from R = 1800 to 2200, and we adopt R = 2000 in our analysis, and this approximation has no effect on our results.
Following du Mas des Bourboux et al. (2020), we discard pixels flagged as problematic in the flux calibration or sky subtraction by the pipeline and and correct for small residual flux calibration errors from the eBOSS pipeline. We then mask out spectral intervals, in observed wavelength, where the variance increases sharply owing to unmodeled emission lines from the sky, as well the wavelength range of Ca II H&K absorption of the Milky Way. Damped Lyman-Alpha systems (DLAs) are identified (Chabanier et al. 2021;Lyke et al. 2020;Parks et al. 2018), and all pixels where the absorption due to the DLA is higher than 20% are not used 3 . Finally, a spectral region is used only if it contains at least 50 pixels to avoid overfitting the mean transmitted flux.

Lyα Absorbers and the Forests
In this work, Lyα absorption systems (Lyα absorbers) are selected based on the fluctuation field of the transmitted flux in the Lyα forest. With a selected Lyα absorber, the fluctuation fields in the Lyman series range and in the range longer than the Lyα wavelength are stacked for our study. We call these two ranges Lyman series forest and metal absorption forest, respectively. To measure the flux fluctuation field of a given forest, we make use of the publicly available Python "Package for 2 https://data.sdss.org/datamodel/files/BOSS SPECTRO REDUX/RUN2D/spAll.html 3 We perform a test of removing pixels with DLA-caused absorption higher than 5% and find that this more conservative cut has little effect on our final stacked results. Also, by default we mask Lyα and Lyβ lines associated with DLAs. We do a further test of masking the metal lines from DLAs. We find that, except weakening the shadow lines in the stacked spectrum, it has no effect on the measurements of metals in the Lyα absorbers we select, as the DLA metal lines do not contaminate them.
IGM Cosmological-Correlations Analyses" picca 4 (du Mas des Bourboux et al. 2021). The measurement of the flux fluctuation field for the Lyα forest, the Lyman series forest, and the corresponding metal forest follows an approach similar to the one established for the Lyα forest BAO analysis, as presented in du Mas des Bourboux et al. (2020). In short, for each quasar q, the fluctuation field δ q (λ) in the transmitted flux for each forest is computed as where f q (λ) is the observed flux density and C q (λ) is the continuum (the flux density that would be observed in the absence of absorption).The quantityF (λ) is the mean transmission. The product C q (λ)F (λ) is the mean expected flux for this quasar. It is modelled as a universal function in quasar's rest-frame (with a uniform forest spectral template), corrected by a linear function of log λ. The model parameters are determined by maximizing the likelihood between the model flux and the observed flux in the forest region for each quasar (see du Mas des Bourboux et al. 2020 for detail). As noted in Bautista et al. (2017) and du Mas des Bourboux et al. 2020, the above fitting method introduces a small systematic bias in δ Lyα . The difference in the correlation functions between the true δ Lyα field and the one derived from the above fitting method based on mock data (Fig.11 of Bautista et al. 2017) allows us to estimate this bias, and we find that the derived δ Lyα is about 0.02 lower than the true value, too small to have any significant effect on our results. Besides the systematic bias, the mean transmitted flux fitting for each individual quasar also introduces a statistical fluctuation in δ Lyα , which is included in the uncertainty of δ Lyα (e.g., du Mas des Bourboux et al. 2020). Such a fluctuation would result in a spread in the true δ Lyα values for a given selected δ Lyα . The high signal-to-noise ratio (above ten; see below) we impose to select Lyα absorbers limits such an effect to be small. In retrospect, this is supported by the clear trend we find in the dependence of metal column density on δ Lyα (to be presented later). Finally, it is possible that the mean transmitted flux fitting is not accurate at other wavelengths (e.g., around a metal line). As mentioned in Section 2.2, the stacked spectrum is normalized by a local linear fit around each metal line, which largely removes the effect caused by the inaccuracy in the above fitting.
With the measured Lyα fluctuation field δ Lyα in the Lyα forest, we construct a catalog of Lyα absorbers.  Figure 1. Illustration of different forests and mean flux fits. The spectra are in quasar's rest-frame. The colored curves are the fits for the mean expected flux [CqF in eq. (1)] in the Lyα forest (red in the left panel), in the Lyman series forest (blue in the left panel), and in the metal forest (red in the right panel). See the text. Note that in the right panel the wavelength range within ±15Å of the quasar's C IV emission is excluded for the CqF fit, as the emission line center are too sharp to produce a good fit. The IDs of the two quasars in the DR16Q catalog are labelled in each panel.
We select Lyα absorbers in the redshift range 2.0 < z abs < 4.0, where z abs is the redshift computed from the Lyα absorption. The lower redshift limit is a consequence of the Lyα forest exiting the observed wavelength coverage for quasars with z 2. The upper redshift limit is adopted as there are only a small number of z > 4.0 quasars in DR16Q. For the purpose of our analysis, the Lyα forest is defined to be in the wavelength range of λ RF =1041-1185Å in quasar's rest-frame (similar to that used by McDonald et al. 2006), where the proximity effects of the quasar's Lyα and Lyβ emission line profiles are small .
In this work, we concentrate on overdense regions and define Lyα absorber pixels as local minima among three adjacent pixels in the Lyα forest. We select absorbers in four bins of flux fluctuation, δ Lyα = −0.95, −0.85, −0.75, and −0.65, respectively, with bin width of 0.1. These are the ranges for absorbers with noticeable metal detection, which is the focus of this work. The analysis of hydrogen distribution properties for absorbers in other ranges of δ Lyα , including those corresponding to underdense regions, will be presented in a separate paper. We require that for each absorber pixel the signalto-noise ratio in δ Lyα is above ten 5 . Such a cut preserves as many absorbers as possible for the study of metal absorption and at the same time it keeps meaningful divisions of the four absorber samples. To investigate the redshift evolution of the metals in Lyα absorbers, we further divide the absorbers into 10 redshift bins, with z abs = 2.1 ± 0.1, 2.3 ± 0.1, ..., and 3.9 ± 0.1. In Table 1, the numbers of Lyα absorbers in different δ Lyα and redshift bins are provided.
Like the Lyα forest, we can define the Lyβ forest region to be in the wavelength range of λ RF =978-1014Å in quasar's rest-frame, and similarly for higher order Lyman series. In practice, we denote the wavelength range of λ RF =787-1185Å in quasar's rest-frame as a single forest, namely the Lyman series forest, even though there are metal lines in this range. For metal lines redward of quasar's Lyα emission, each can have its own forest region. For example, in du Mas des Bourboux et al. (2019), they adopt quasar rest-fame wavelength range of 1260-1375Å as the Si IV (1394/1403Å doublet) forest and 1420-1520Å as the C IV (1548/1551Å doublet) forest. In this work, we define a single metal forest as the range of λ RF =1236-2776Å in quasar's restframe. The range corresponds to ∼ 5000 km s −1 redward of the quasar's Lyα emission and ∼ 2600 km s −1 blueward of the quasar's Mg II emission, and we also exclude the pixels within ±15Å (in quasar's rest-frame) around the quasar's C IV emission, reducing the effect of the quasar's emission lines on the C qF fit.
Figure 1 illustrates the different forests and the C qF fits in equation (1) for each forest.
With quasars in the redshift range of 2.08 < z q < 4.82, applying the absorber redshift cuts of 2 < z abs < 4 and strength cuts of −1.0 < δ Lyα < −0.6, we end up with 47,118 Lyα forests (λ RF =1041-1185Å), 46,804 Lyman series forests (λ RF =787-1185Å), and 47,076 metal forests (λ RF =1236-2776Å). Lyα absorbers selected from the Lyα forests will be correlated with flux fluctuations in the Lyman series forests and metal forests to study the metal content associated with those Lyα absorbers.

Stacked Spectra Based on Lyα Absorbers
The selected Lyα absorbers here are connected to neutral hydrogen in the CGM or IGM, and they can be affected by metal enrichment caused by star formation and feedback processes. While the corresponding metal    Fig. 9) correspond to 10 absorber redshift bins, from z abs = 2.1 to 3.9 with step 0.2 and roughly with lower curve for higher redshift (at rest-frame 1500Å). absorption features may be too weak to be identified in the eBOSS quasar spectra for individual absorbers, they can be revealed by stacking the spectra associated with a large number of Lyα absorbers (e.g., Pieri et al. 2010Pieri et al. , 2014. With the catalog of Lyα absorber pixels and the fluctuation fields computed for the Lyman series and metal forest at different wavelengths, we perform a stacking analysis to reveal metal absorption associated with Lyα absorbers of different strengths (δ Lyα ). We follow the same technique as in Font-Ribera et al. (2014) and du Mas des Bourboux et al. (2017) to do the stacking. That is, we compute the cross-correlation function between Lyα absorber pixels and the fluctuation field δ of various forests, which is carried out with picca. In detail, the cross-correlation is computed as where i denotes the forest pixel and k the Lyα absorber pixel, δ i is the flux fluctuation in one of the forests, and w i and w k are the weights [as defined in equations (4) and (7) of du Mas des Bourboux et al. 2020] 6 . The wavelength bin A represented by λ is determined through λ = λ Lyα (λ i /λ k ), where λ Lyα = 1215.67Å is the restframe wavelength of the Lyα line and λ i and λ k are the observed wavelengths of the forest pixel and Lyα absorber pixel, respectively. Clearly, λ is the wavelength in the absorber's rest-frame. Later we also use δ(λ) to represent ξ(λ). The computation is performed for absorbers in each δ Lyα bin and in each redshift bin. Because of systematics in continuum fitting, the upper envelope of each stacked spectrum, 1 + ξ, can slightly deviate from unity, and we apply a low-pass pseudo continuum to rectify this (e.g., Pieri et al. 2014). Specifically, for illustrating the stacked spectrum, we perform a spline fit to the 1 + δ(λ) in regions excluding absorption lines and normalize the spectrum by the fit. For analyzing each absorption line, we perform a local linear fit to the continuum and normalize the absorption line profile by the fit.
The covariance matrix of the cross-correlation is calculated by subsampling the data, similar to the approach in du Mas des Bourboux et al. (2017). We divide the sky into HEALPix pixels (Górski et al. 2005) and compute the cross-correlation using data in each HEALPix pixel (subsample). Using a division of the sky with nside=16, we obtain the number of subsamples to be between 27 and 790, depending on δ Lyα and z abs bins. The covariance is then given by the following, neglecting the small correlations between subsamples: Here ξ α = ξ(λ α ) is the mean measured correlation for wavelength λ α , s denotes a subsample with summed weight W s α and measured correlation ξ s α , and W α = s W s α .
The signal-to-noise ratio of the stacked spectra depends on the number of Lyα absorber-forest pixel pairs. In Figure 2, we show the dependence of the number of such pairs on wavelength, δ Lyα , and z abs . Given our selection, the number has almost no dependence on δ Lyα and only a weak dependence on wavelength. The redshift dependence has two parts. First, the coverage shifts towards shorter wavelength at higher z abs , a reflection of redshift effect with the fixed observed wavelength range. At the lowest absorber redshift, z abs ∼ 2.1, we are not able to see Lyman series beyond Lyα while at the long wavelength end it covers up to ∼2900Å. In contrast, at the highest absorber redshift, z abs ∼ 3.9, at the low wavelength end the Lyman limit can be covered, but at the long wavelength end we cannot go above ∼2150Å. Second, at a given wavelength, the redshift distribution of quasars leads to a strong redshift dependence of absorber-forest pairs, varying from a few times 10 3 (z abs ∼ 2.1) to a few times 10 1 (z abs ∼ 3.9). Therefore, we expect to have less constrained and fewer metal lines at higher absorber redshift. Since δ Lyα is already selected to have signal-to-noise ratio higher than 10, the number of pairs also indicates that the stacked spectra for each δ Lyα and at each z abs have a much higher signalto-noise ratio, making it possible to identify weak metal absorption lines. We note that there is a break around 1280Å, since for the absorber-forest pixel pairs we use the Lyα and Lyman series forest (the metal forest) below (above) this wavelength for the forest pixels, which is also adopted in our stacked spectra.
As an illustration, in Figure 3, we show the stacked spectra for Lyα absorbers in the whole redshift range of 2 < z abs < 4 with −1.0 < δ Lyα < −0.9. The Lyman series, from Lyα to Ly11, can be easily identified. We leave the discussion of Lyman series (including those in under-dense regions) to a separate paper and focus our discussion here to metal lines that show up in relatively strong Lyα absorption systems.
The high signal-to-noise ratio of such stacked spectra clearly helps reveal an array of metal lines in the wavelength range of 900-2880Å, including low-ionization lines (e.g., Mg II and Fe II) and high-ionization lines (e.g., O VI and C IV). In particular, the lines with strong oscillator strengths are very prominent (e.g., Si IV 1394/1403Å, C IV 1548/1551Å, and Mg II 2796/2804Å doublets). There also appears to have a series of shadow lines. For example, the doublets at λ ∼ 1560Å are the shadow of the  are not labeled in the figure, and a similar plot with major shadow lines identified and labelled can be found in the Appendix B ( Figure 15). The lines contaminating Lyα are weak, and they play no significant role in our δ Lyα selection. For each identified Lyman series line and metal line, we perform Voigt profile fitting to obtain the column density of each species. In doing the fitting, for a set of model parameters (column density N and Doppler/broadening parameter b) the Voigt model is convolved with a Gaussian kernel of full-width-at-halfmaximum ∼ 150 km s −1 to account for the eBOSS spectral resolution. For a given line, we typically choose more than 5 pixels on each side to determine the continuum, and we visually inspect the pixels to remove those belonging to shadow lines. For example, for Si IV 1394Å, there is a shadow line of Si II 1527Å at λ ∼ 1391Å caused by the contamination of C II 1335Å in Lyα; for Si IV 1403Å, there is a shadow line of Si IV 1394Å at λ ∼ 1404Å caused by the contamination of Si III 1207Å in Lyα. In fitting the lines, we also exclude the pixels near the shadow lines associated with the contamination of transitions around the Lyα line (such as Si II 1190Å, Si II 1304Å, and C II 1335Å). For example, the absorption feature around 1404.5Å in the third panel of Figure 7 is the shadow of the Si IV 1393Å associated with the contamination of Si III 1207Å in Lyα. The pixels around this feature are not included when we perform Voigt fitting for the Si IV 1402Å line. Following previous work (e.g., Pieri et al. 2014), for doublet, we choose to fit each line separately, and the difference in the constraints can provide an assessment of possible systematics. For strongly blended lines that are hard to be measured independently, such as the blend of C II and O VI around 1037Å and that of O I and Si II around 1303Å, we perform a joint fit.
In general, we find that consistent with previous work (e.g., Pieri et al. 2014) the broadening parameter b appears to be high, typically at the level of a couple of hundred of km s −1 . This is related to the resolution of the eBOSS spectra. The relatively low resolution means that each identified Lyα absorber can consist of multiple absorbers, and the b parameter includes the contributions of thermal motion in each individual absorber and the relative motion among these absorbers. This work focuses on the column densities of various species. In Tables 2 in the Appendix A, the parameter constraints from Voigt fitting for each metal line as a function of δ Lyα and z abs are provided.

RESULTS
With the stacked spectra, we first derive the H I column density of each set of Lyα absorbers by analyzing the Lyman series lines. Then we study how the metal lines depend on the strength of Lyα absorbers and their redshift evolution. Finally, we present an illustrative model of metal column densities to provide insights on the physical properties of Lyα absorbers selected in this work. Figure 4 shows the stacked H I Lyman series absorption profiles as a function of δ Lyα . Compared with previous work (e.g., Pieri et al. 2014), the substantial increase in the signal-to-noise ratio allows us to reach Lyman series of much higher order. Four absorber redshifts are chosen for the illustration, which have clear detection up to Lyλ (i.e., Ly11). For each set of a given δ Lyα , the absorption depth decreases with decreasing line wavelength, resulting from the lower oscillator strength f for higher-order Lyman lines.

Stacked Lyman Series and H I Column Density
We perform Voigt profile fitting for each Lyman series line to obtain the column density N HI and the b parameter (and hence the equivalent width). An example is shown in Figure 5 for z abs ∼ 2.9. As pointed out in Pieri et al. (2014), lower-order Lyman series lines (e.g., Lyα and Lyβ) are challenging to be used to derive precise measurements of the equivalent width. Higher-order, unsaturated lines are more suited for such a task, and the equivalent width (normalized by f λ) is expected to approach a constant value. As we are able to reach Lyman series lines of much higher order, we indeed find that the equivalent width reaches a plateau for each set of lines of a given δ Lyα (as shown in Fig. 5). This can be seen more clearly in the middle panel, where the column densities from higher-order Lyman series lines are consistent with being at a constant value within the error bars.
In what follows, the H I column density is presented together with those of metal lines and adopted in modelling the absorber systems.
At each redshift, for each metal line, there is a clear trend that the line strength increases with increasing |δ Lyα |, implying that more metals are present in denser regions. At fixed δ Lyα , each metal line becomes weaker toward higher redshift.
At a given δ Lyα and z abs , Voigt profile fitting is carried out for each metal line, as long as the signal-to-noise ratio allows. Specifically, if the minimum pixel value of an absorption line lies below 1σ of the noise level of the "continuum", we perform the profile fitting. The bestfit is shown as the red solid curve. The derived column density N and Doppler parameter b can be found in Tables 2 in Appendix A. Overall, if the results are averaged over the redshift range, we find broad agreements with those in Pieri et al. (2010) and Pieri et al. (2014). Figure 9 summarizes the dependence of the column density on δ Lyα and z abs for various metal species as well as neutral hydrogen. In each panel, up to four sets of data points are shown, corresponding to the four values of δ Lyα (−0.95, −0.85, −0.75, and −0.65). Each set of data points are color-coded by the absorber redshift z abs . The column density ranges of the panels are all chosen to across 3 dex for an easy comparison between the strengths of the column density dependence on δ Lyα and z abs of different species.
At each redshift, the column density of each metal species drops as δ Lyα increases from −0.95 to −0.65. The column density of neutral hydrogen (top-left panel) drops by about one dex, from log(N HI /cm −2 ) ∼16 to ∼15. The change of column densities of high-ionization metal species (such as O VI, N V, C III, and C IV) seems to track that of neutral hydrogen, all with a ∼1 dex decrease over the δ Lyα range. In contrast, the column densities of low-ionization metal species show a stronger dependence on δ Lyα . For example, that of Mg II drops by almost 2 dex over the δ Lyα range. Such a difference implies that high-and low-ionization species arise from different environments, as will be shown with a simple model in the next subsection. At a given δ Lyα , the column density of neutral hydrogen show almost no dependence on redshift. In comparison, those for metal species at fixed δ Lyα show a clear dependence on redshift, typically more than 0.5 dex from z abs ∼ 2.1 to z abs ∼ 3.5, which may result from the evolution in temperature and ionizing background.
To help understand these results, we turn to a simple model of Lyα absorbers.

Illustrative Model of Metals in the Lyα Forest of High Neutral Hydrogen Column Density
The column densities of various species encode information about the physical properties of the absorption system. Similar to Pieri et al. (2014), we apply a simple, illustrative photoionization model to interpret the column density measurements.

Parameter Dependence
To guide our choice of models to compare to the measurements, in Figure 10, we show the dependence of the column density of various species on the stop neutral hydrogen column density (N HI ), total hydrogen number density (n H ), and temperature (T ), with gas in a planeparallel geometry photoionized by an external ionizing UV background. The nearly linear dependence of the column density of a species on metallicity (with other parameters fixed) is not shown here. Figure 14 in Appendix A further explores the parameter dependence.
In Figure 10, the species are ordered such that the ionization potential decreases from left to right. Calculations are performed with version 17.02 of photoionization code Cloudy, last described by Ferland et al. (2017). The gas was assumed to be a uniform slab of constant density in thermal and ionization equilibrium and exposed to the external UV background. The models explored are centered at a model with total hydrogen density log(n H /cm −3 ) = −3 and temperature log(T /K) = 4.5. The metallicity is fixed at [X/H] = −1. Solar abundance is assumed and the Haardt & Madau (2012) UV background at z = 2.75 is adopted in the calculation. With the central model, we vary the three parameters, N HI , n H , and T , one at a time. Each panel shows one set of change, and thicker lines are for higher values of the corresponding varied parameter.
As with the metallicity dependence, there is a clear trend that a species' column density increases with increasing stop neutral hydrogen column density (top panel). The dependence is stronger for low-ionization species. For example, for a 3 dex change in N HI , the column density of Mg II increases by about 4.5 dex, while the change is only about 1.5 dex in that of O VI. This is expected in the photoionization model. At higher hydrogen column density, the ionizing background is more attenuated within the absorption system. The subsequent lower photoionization rate causes a given species to stay at a low-ionization state than a high-ionization state.
As for the dependence on total hydrogen density n H , a higher n H results in a higher recombination rate, making it hard for a species to maintain a state of high-level ionization. As a consequence the column densities of high-9 7 7 9 7 7 9 7 7 9 7 7 9 7 7 (Å) 1 1 9 0 1 1 9 0 1 1 9 0 1 1 9 0 1 1 9 0 1 1 9 0 1 1 9 0 1 1 9 0 1 1 9 3 1 1 9 3 1 1 9 3 1 1 9 3 1 1 9 3 1 1 9 3 1 1 9 3 1 1 9 3 1 2 0 6 1 2 0 6 1 2 0 6 1 2 0 6 1 2 0 6 1 2 0 6 1 2 0 6 1 2 0 6 (Å) 1 2 6 0 1 2 6 0 1 2 6 0 1 2 6 0 1 2 6 0 1 2 6 0 1 2 6 0 1 2 6 0 1 3 0 3 1 3 0 3 1 3 0 3 1 3 0 3 1 3 0 3 1 3 0 3 1 3 0 3 1 3 0 3  Figure 6. Dependence of line profiles in the stacked spectra on δLyα and z abs for 13 metal species. In each panel, the relevant species are labelled at the lower-left corner, and the horizontal tick mark interval is 1Å. The δLyα and z abs dependences are shown to be along the vertical and horizontal direction, respectively. Offsets are added for different values of δLyα. The vertical lines mark the wavelengths of the relevant lines, and the red curves show the Voigt profile fits. In each row, the cases for z abs < 3.6 and 3.6 < z abs < 4.0 absorbers are plotted in the left and right panel, respectively, and the division is mainly driven by the potentially different dynamical range in the y-axis of the right panel to accommodate the larger uncertainties in the profiles at higher redshifts.  1 3 9 3 1 3 9 3 1 3 9 3 1 3 9 3 1 3 9 3 1 3 9 3 1 3 9 3 1 3 9 3     ionization species drop by a large factor (middle panel). For example, as n H increases from log(n H /cm −3 ) = −3 by 1 dex, the column densities of O VI, N V, and C IV decrease by about 3 dex, 2.5 dex, and 2 dex, respectively. The recombination from higher-level ionization, on the other hand, increases the column densities of lowionization species by a mild amount, e.g., ∼1 dex for Mg II. In contrast, lowering n H by 1 dex leads to a substantial drop in the column densities of low-ionization species (e.g., by ∼2.5 dex for Mg II) and a clear increase in column densities of high-ionization species (e.g., by ∼2.5 dex for O VI). As we will see, such an increase for high-ionization lines at low density provides a possible mechanism to explain the measurements. Increasing the temperature by 0.5 dex from log(T /K) = 4.5 (hence reducing the recombination rate) has the effect of enhancing the high-ionization species relative to low-ionization ones (bottom panel). Compared to the lower density case (middle panel), the relative abundances of high-ionization lines do not change much. For example, the column density of O VI still appears to be lower than that of C IV. The column density change in the low-ionization species is not as drastic as the lower density case. Compared with the case of increasing temperature by 0.5 dex, decreasing the temperature by 0.5 dex only has a small effect, shifting the column densities by at most 0.5 dex.

A Two-Component Model
To see how the measured column densities of various species compare with Cloudy models, we use the absorption system selected by δ Lyα = −0.95 and 2.6 < z abs < 2.8 as an example. The data points in the left panel of Figure 11 show the inferred column densities (see Table 2 for 2.6 < z abs < 2.8 absorbers). For those lines with column density measurements barely missed at this redshift but with measurements at the neighboring redshifts, such as C II 1036Å and O VI 1038Å, we apply interpolations to obtain the z ∼ 2.7 values. For the same species with multiple lines, we find that the (logarithmic) column density measurements from different lines can have a systematic difference up to 0.10-0.25 dex (e.g., C IV or Fe II), and we take the maximum column density as the measurement 7 .
Overall, the column densities of high-and lowionization species are within a relatively narrow range of about 2 dex. In contrast, the single component models presented in Figure 10 all show a large range in column density. In particular, these models are unable to make the column densities of high-ionization species to be at a comparable level as those of low-ionization species, as seen in the data. Motivated by Pieri et al. (2014), where the high-ionization and low-ionization species are compared to different models, we explore the scenario that the inferred column density pattern results from a superposition of the contributions from a high-metallicity and a low-metallicity component.
We build two sets of Cloudy models with gas in a plane-parallel geometry, one with metallicity [X/H] = −2 and the other with [X/H] = −0.3. For each set, we compute models on a grid, log(n H /cm −3 ) = −4 to 0 with step 0.05 and log(T /K) = 4 to 6 with step 7 We find that the discrepancies between column densities from the different lines (e.g., doublets) of an ion are related to the difference in their cross-sections (oscillator strengths). The line with the larger cross-section is more easily saturated. With the lowresolution spectra, the Voigt fit has the tendency of pushing the model toward the linear regime by fitting it with a large b value. This would lead to a smaller fitted column density. Therefore, the column density from the line with the minimum oscillator strength is a better estimate. We perform tests of using the inverse variance weighted mean column densities in doing the fitting with the Cloudy model, and the results remain similar.  log NHI=16, log nH=-3, log T=4.0 log NHI=16, log nH=-3, log T=4.5 log NHI=16, log nH=-3, log T=5.0 log NHI=16, log nH=-3, log T=5.5 Figure 10.
Dependence of column densities of various species on neutral hydrogen column density (top), total hydrogen number density (middle), and temperature (bottom), based on Cloudy models. The species are ordered from highionization to low-ionization potential from left to right. The metallicity is fixed to [X/H] = −1 for all the cases. 0.05, all with stop neutral hydrogen column density log(N HI /cm −2 ) = 16 assuming the solar abundance and the Haardt & Madau (2012) UV background at z = 2.75. Note that the mean hydrogen density at this redshift is about 10 −5 cm −3 and the hydrogen density in the set of models covers 10 to 10 5 times the cosmic mean. For each combination of the low-metallicity and high-metallicity models, the relatively fraction is solved to best fit the measurements. Our purpose here is to look for combinations that can provide a reasonable explanation to the overall trend of the column densities of various species. If we use the measurement uncertainties, the results would be completely driven by the few species with small error bars (e.g., 0.02 dex for C III and Si IV in Table 2 for 2.6 < z abs < 2.8). At such an uncertainty level, the model assumptions may also matter (e.g., the adopted solar abundance pattern). To focus on the overall trend and given the possible systematics in deriving the column densities (e.g., seen in the differences from multiple  Figure 11. Illustration of the combination of two absorption systems to compare to the observed column densities of various species with δLyα ∼ −0.95 and 2.6 < z abs < 2.8. Top: the best fit (black) from a combination of a low-metallicity system (red) and a high-metallicity system (blue). Bottom: the 1σ and 2σ constraints in nH and T for the combination of a low-metallicity and a high-metallicity system. lines), we adopt a uniform error bar of 0.2 dex in this exercise. The solid line in the top panel of Figure 11 is the bestfit model, consisting of a small contribution (2.3%) of low-metallicity systems with log(n H /cm −3 ) = −2.9 and log(T /K) = 5.4 (red dotted line) and a large contribution (97.7%) of high-metallicity systems with log(n H /cm −3 ) = −1.5 and log(T /K) = 4 (blue dotted line). To have a better idea of the possible combinations, we show the constraints on n H and T in the bottom panel. The high-metallicity component is well constrained to be around log(n H /cm −3 ) ∼ −1.5 and T ∼ 10 4 K. The low-metallicity component, on the other hand, is only loosely constrained. It can come from systems with either low density (∼ 10 −4 cm −3 ) or high temperature ( 10 5 K), both driven by the high column densities of high-ionization species (e.g., O VI, N V, and Open circles show the observed column densities of various species with δLyα ∼ −0.95 and 2.6 < z abs < 2.8. A high-columndensity system (black dotted line) is added, and the combination with the other two systems (red and blue dotted line) is shown as the solid line, with the combination fractions labelled in the legend (not a model fit). This is only for illustration purpose, as there is degeneracy among model parameters (e.g., column density and metallicity) and the system is unlikely to be limited to three components. See the text for more detail. C IV; see the parameter dependence in Fig. 10). Such physical conditions appear to be in line with the prediction for highly ionized metals from hydrodynamic simulations (e.g., Rahmati et al. 2016). The high-density component likely comes from cool gas in the CGM, while the other one is contributed by gas in the IGM. The two-component model is consistent with the finding based on high-resolution observation that metal absorbers show bimodal physical properties with low-and high-metallicity branches (Kim et al. 2016).
We note that we have performed a test with Cloudy by switching to the HM05 UV background model. We find the overall trend, e.g., as seen in Figure 11, remains almost the same, with a slight shift in the number density direction. For example, the bestfit number density for the low-metallicity component shifts by about 0.2 dex, from log(n/cm −3 ) ∼ −1.5 to -1.3.

Multi-Component Model
It should be emphasized that the exercise shown in Figure 11 indicates that the selected absorption systems cannot be described by a simple single component model. However, it does not necessarily mean that the right model is the two-component model, which is merely used here to illustrate the possible complexity of such systems. In fact, the fit shown in the top panel of Figure 11 misses a few data points, both in the highand low-ionization species.
For example, the pattern in the high-ionization species is not fully reproduced, with the predicted C IV and Si IV lower than measured. While the column densities of most low-ionization species can be reasonably described by the model, the model underpredicts the column densities of Fe II and Mg II by about 0.5 dex and that of O I by about 1.5 dex.
An adjustment of metallicity can bring the model to better match the column densities of low-ionization species.
With N HI = 10 16 cm −2 , we find that a model with super-solar metallicity ([X/H] = 0.3) and log(n H /cm −3 ) = −1 can lead to a better match to the column densities of Fe II and Mg II, but the predicted O I column density still falls short of the measurement by ∼0.75 dex. We note that in Pieri et al. (2014) a model with [X/H] = 0.3 is able to provide a match to their O I column density measurement (see their Fig.12). The difference may have two causes. First, we have a different selection of hydrogen absorption systems and the measured O I column density, log(N OI /cm −2 ) = 13.77, from our selection is higher than theirs (by about 0.37 dex). Second, the solar abundance pattern used in our work appears to be different from theirs. As it is unlikely to have metallicity much higher than [X/H] = 0.3, within our adopted model we find it is hard to explain the high O I column density measurement solely by metallicity effect.
The mismatch between the measured and model O I column density could suggest that O is more abundant in real systems, i.e., the abundance is different from the solar abundance adopted in the calculation. This may also help resolve the difference seen in the highionization species, which relates to the CNO enrichment history. The high O/N ratio can be caused by a series of major star formation episodes in galaxies with short interval between them (e.g., Kobulnicky & Skillman 1998;Pettini et al. 2002;Pettini 2004), given that O is produced in core-collapse supernovae (> 10 Myr after star burst) and N comes from intermediate-mass stars (4-8M ; > 100 Myr after star burst). Galactic winds may bring such an abundance pattern into the CGM/IGM, where the absorption systems likely reside.
There are certainly other possibilities. Given the resolution of the eBOSS spectra, the strong absorption systems we select, in this case, δ Lyα ∼ −0.95, may have a small fraction of contamination from high column density systems. In high-column-density systems (e.g., Lyman limit systems), the column density of O I, the one with the lowest ionization potential for species considered here, can become enhanced, owing to the low photo-ionization rate inside the systems from the attenuated ionizing background (see the bottom two panels of Fig. 14 in Appendix A). Such an enhancement in O I column density is seen in DLAs (see e.g., Cooke et al. 2011a,b). Figure 12 illustrates a threecomponent model (not from model fitting). Loosely motivated by Figure  , which is able to produce the high column density of O I and improves the match to those of Fe II and Mg II. In combination, the three components lead to a trend following the measurements reasonably well. Increasing the metallicity, n H , and N HI all have the effect of boosting the column density of O I, the species with the lowest ionization potential, and the latter two have the additional effect of causing a steeper jump between Fe II/Mg II and O I column densities. The contamination from the third component is assumed to be 5% in the illustration. While the metallicity [X/H] = −1 is motivated by the measurements in Lyman-limit systems (e.g., Fumagalli et al. 2016;Lehner et al. 2016), given the degeneracy among [X/H], n H , and N HI , the values used here are only for the purpose of illustrating the possible contribution from high-column density systems. In addition, there is hardly any reason to limit to three components for the selected absorbers.
If we broadly take the contamination from Lyman limit systems to be at a level of a few percent, is it reasonable? To answer this we must understand the systems that meet the selection function. Similar to the case of flux transmission F < 0.25 in Pieri et al. (2014), isolated lines without damping wings do not reach −1 < δ Lyα < −0.9 in moderate resolution spectra (such as those observed in eBOSS). The selection function requires lines that are both saturated in high resolution and clustered with one another on scales within the full width half maximum of the SDSS resolution element (see Fig.2 of Pieri et al. 2014). For this reason we cannot simply integrate the column density distribution function in order to assess the proportion of Lyman limit systems present in the selected absorbers. Put in logic terms, the presence of N HI 10 14 cm −2 is a necessary but not sufficient condition for selection, since clustering is also required; however, N HI > 10 17.5 cm −2 is neither neces-sary nor sufficient, since such systems are not necessarily sufficiently clustered with other N HI > 10 14 cm −2 lines on ∼ 150 km s −1 scales. Hence the column density distribution function is largely orthogonal to our selection function.
In order to understand the populations of individual lines selected here we must compare with results from line fits of high resolution data and verify what the selection function provides. This was done in Pieri et al. (2014) with the conclusion that N HI > 10 17.5 cm −2 systems can contribute no more than 3.7% of the systems with a true flux transmission of F < 0.25 with SDSS resolution. While our selection function is for stronger features than those in Pieri et al. (2014), the selection function is otherwise the same. We have no basis to conclude that the somewhat stronger blended absorption here differs from the line model test in Pieri et al. (2014). We therefore also assume their conclusion that systems with N HI > 10 17.5 cm −2 should be no more than 3.7% of our sample.
The discussion here focuses on δ Lyα ∼ −0.95 systems. The O I line becomes too weak to be measured robustly for δ Lyα −0.85 systems (bottom panel of Fig. 6), and this is nearly so for Fe II and Mg II as well. Such a behavior lends support to the scenario of contamination from high-column-density systems to the strongest absorption systems (δ Lyα ∼ −0.95) in this study with eBOSS spectra. In Appendix C, we perform further tests through selecting presumably stronger absorbers within the δ Lyα ∼ −0.95 systems. The change in the metal column densities and the feature at the corresponding Lyman limit also provide further evidence on the contribution from high-column-density systems.
To summarize, the column densities of various lowionization species in δ Lyα ∼ −0.95 absorption systems at 2.7 < z abs < 2.9 can be reasonably explained by a model cloud with n H ∼ 10 −1.5 cm −3 and T ∼ 10 4 K, photoionized by a UV background (Haardt & Madau 2012) to have a neutral hydrogen column density of N HI ∼ 10 16 cm −2 . The relatively high column densities of high-ionization species suggests an additional component, with either low n H or high T , or a deviation from solar abundance. The high column density of O I, the one with the lowest ionization potential for species considered in this work, may signal contamination from high column density systems in view of the spectral resolution of eBOSS.
Multiple components appear to also apply to other δ Lyα systems and to systems at different redshifts. In Figure 13, the column densities of various species at different δ Lyα values are shown in different panels, and in OVI NV CIV CIII SiIV SiIII AlIII CII AlII SiII FeII MgII OI 10% red + 90% blue Figure 13. Column densities of 13 metal species as a function of δLyα and z abs . The four panels are for different values of δLyα. In each panel, the column density of each species is color-coded according to absorber redshift, low z abs to high z abs from left to right within each species. The solid black line illustrates the model combination of a low-metallicity component (low nH and high T ; red-dashed line) and a high metallicity component (high nH and low T ; blue dashed line), and the difference in the models across the four panels is in the neutral hydrogen column density, as labelled.
each panel the redshift dependence is color-coded. The Cloudy model (not fit) with a high-metallicity and a low-metallicity component of different values of n H and T is used to illustrate the possible multiple components. Histogram, instead of lines connecting column density of each species, is used here for the model and its components to avoid the possible confusion with the redshift dependence. The models across the four values of δ Lyα are the same, except for the stop neutral hydro-gen column density. The O I line is not robustly measured, except for the strongest absorption system with δ Lyα ∼ −0.95, consistent with it being contaminated by systems of high column density. Given the uncertainties discussed above, we make no attempt to further model the redshift dependence, which can be related to metallicity and density evolution of the selected systems.

SUMMARY AND DISCUSSION
We study metal absorption lines associated with relatively strong Lyα absorbers in the Lyα forest using quasar spectra from the eBOSS survey. Various metal species are revealed in the stacked spectra, and their column densities are derived as a function of Lyα absorber strength and redshift. We find that multiple components, including low-and high-metallicity gas of different densities and temperatures, are needed to interpret the inferred column density pattern of the metal species. The results can be used to study the chemical enrichment of CGM and IGM. They also provide inputs to model the metal contamination in the Lyα forest BAO measurement.
We construct a catalog of Lyα absorbers in 10 redshift bins over 2 < z abs < 4, selected according to the Lyα flux fluctuation δ Lyα . We focus on the relatively strong absorbers, with δ Lyα ∼ −0.95, −0.85, −0.75, and −0.65. The Lyα absorbers in each redshift bin are crosscorrelated with the flux fluctuation in quasar spectra in the full wavelength range (shifted in accordance with the Lyα absborber redshift) to obtain the stacked spectrum. The neutral hydrogen column density of these absorbers, from analyzing higher-order Lyman series lines, ranges from ∼ 10 15 cm −2 (for δ Lyα ∼ −0.65) to ∼ 10 16 cm −2 (for δ Lyα ∼ −0.95), which has only weak redshift evolution.
In the stacked spectrum, up to 13 metal species can be identified, including high-ionization ones (e.g., O VI, N V, and C IV) and low-ionization ones (e.g., Al III, Al III, C II, Si II, Fe II, Mg II, and O I). The column densities of these species drop as δ Lyα increases and toward higher redshift.
We use a Cloudy model with gas cloud photoionized by an external UV background to help understand the results. In agreement with previous work (Pieri et al. 2010(Pieri et al. , 2014, we find that the column density distribution of the various metal species can not be explained by a single gas component. In the model we consider, the low-ionization species point to absorption systems of metallicity [X/H] ∼ −0.3 with high density (n H ∼ 10 −1.5 cm −3 ) and low temperature (T ∼ 10 4 K). The high-ionization species are likely low-metallicity ([X/H] ∼ −2) with either low density (n H ∼ 10 −4 cm −3 ) or high temperature (T > 10 4.5 K). The two compo-nents likely correspond to cool gas in the CGM and gas in the IGM, respectively. For absobers selected to have δ Lyα ∼ −0.95, there is also a possibility of contamination from Lyα absorbers of higher neutral hydrogen column density.
The model explored in this work is for illustration purpose, highlighting the need for more than one component. In reality, for systems selected by δ Lyα there may be more than two or three components with a distribution of gas properties (n H , T , abundance pattern, etc). Large-volume hydrodynamic simulations of galaxy formation (e.g., Oppenheimer et al. 2012) with photoionization computation for metal enriched gas (e.g., Oppenheimer & Schaye 2013) can help guide the modelling effort by identifying Lyα absorption systems in the synthetic spectra (matched with eBOSS resolution) and connecting them to the physical properties of the underlying gas (e.g., Turner et al. 2016). Meanwhile, highresolution observation of the Lyα forest (e.g., Cowie et al. 1995;D'Odorico et al. 2016) can be used to study the components in the Lyα absorbers identified under eBOSS resolution and their associated metal properties. The complementary use of the high-resolution spectra and stacked eBOSS spectra would provide an ideal approach to understanding the metal content in the Lyα forest.  Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.  In Tables 2, we list the column density measurements of all metal species in the ten absorber redshift bins over 2 < z abs < 4 and in the four δ Lyα bins at each red-shift. The measurements result from Voigt fits with the stacked spectra. The vacuum wavelengths and oscillator strengths used in the fits are taken from the Atomic Line List 8 (Krogager 2018).
In Figure 14, we extended the exploration of Cloudy models in the main text, overlaid with the measured column densities. The temperature is fixed at T = 10 4 K. From top to bottom panels, the density n H increases from 10 −4 cm −3 to 1 cm −3 . Note that in the fourth panel the neutral hydrogen column density N HI changes from 10 16 cm −2 to 10 18 cm −2 . From the top four panels, we see that as n H increases, the overall trend is an increase in the column density of low-ionization species and a decrease in those of high-ionization species, a result of the increase in the recombination rate. The dependence on metallicity is to change the amplitude. The bottom panel shows the dependence on neutral hydrogen column density N HI . At a higher N HI , species with the lowest ionization potential tends to have a stronger increase in column density.

B. SHADOW LINES
In the stacked spectrum shown in Figure 3, we label the major metal lines associated with the selected Lyα absorber. Many shadow lines can be seen in the stacked spectrum. These shadow lines are caused by metal lines in Lyα absorbers with wavelengths close to Lyα -the Lyα absorption used to select a type of absorbers of interest (e.g., with −1.0 < δ Lyα < −0.9) can be contaminated by metal absorptions from neighboring Lyα absorbers. Then in the stacked spectrum of the specific absorbers of interest, the metal lines in the neighboring absorbers show up with shifted wavelengths, emerging as shadow lines.
As the shadow lines are not directly related to what we intend to study in this work, they are not labelled in Figure 3 to avoid any confusion. In Figure 15, we provide a version of the stacked spectrum with major shadow lines identified and labelled to aid any possible investigation of shadow lines.
First, in the second row from the top we label the five metal lines in the Lyα absorber with wavelengths close to Lyα, including Si II 1190 & 1193Å (blue), Si III 1207Å (cyan), Si II 1260Å (magenta), and C II 1335Å (red). If one of these lines with restframe wavelength λ contamination contaminates the selected Lyα absorption, a metal line with restframe wavelength λ line associated with the absorber where the contaminating line originates will show up as a shadow line at the following OVI NV CIV CIII SiIV SiIII AlIII CII AlII SiII FeII MgII OI  Figure 14. Dependence of column densities of metal species from Cloudy model. The gas temperature is fixed at T = 10 4 K. Across the top four panels, the total hydrogen density nH increases, and in each of these panels, different lines show the dependence on metallicity. Note that in the fourth panel the neutral hydrogen column density NHI changes from 10 16 cm −2 to 10 18 cm −2 . The bottom panel shows the dependence on NHI, with other parameters fixed. For comparison, the measured column densities for δLyα ∼ −0.95 absorbers at different absorber redshift (color-coded) are overlaid.
wavelength, We perform further tests to investigate the contribution of high-column density systems in the −1 < δ Lyα < −0.9 absorbers.
By default, an absorber is selected to be in the sample by requiring the corresponding pixel with −1 < δ Lyα < −0.9 to be the local minimum. With such a selection, the δ Lyα values of its neighbouring pixels can be either below or above −0.9. If the absorber has a higher column density, its neighbouring pixels more likely have δ Lyα lie below −0.9. We therefore impose an additional selection criterion by requiring the δ Lyα values of the pixels to the left and that to the right of the above selected pixel to be less than −0.9. More exactly, we require δ Lyα,i < −0.9 + σ δ , where σ δ is the uncertainty in δ Lyα and i indicates the i-th pixel next to the selected absorber, negative (positive) i for pixels to the left (right).
We create three subsamples, subsample-1, subsample-2, and subsample-3, by requiring 1, 2, and 3 left/right neighbouring pixels to satisfy the above additional selection, respectively. These subsamples potentially have increasingly higher column densities. As illustrated in Figure 16, in the case of a single absorber under the SDSS resolution, the subsample-1 selection prefers to have systems of log(N HI /cm −2 ) > 18 or those of log N HI > 17 with large b values, and subsample-2 and subsample-3 selections would select absorbers of log N HI > 18.5. The absorbers in the three subsamples follow the same redshift distribution as the original sample, as shown in Figure 17. Figure 18 compares the stacked spectra of the original −1 < δ Lyα < −0.9 sample and the three subsamples. It is clear that the stacked Lyα line becomes increasingly broader from the original sample to the subsample-1/2/3, and the sequence also displays an increasingly stronger wing-like feature. The Lyman series lines appear to become broader and deeper. All of these are consistent with the expectation that the subsamples consist of absorbers with an increasing contribution from systems of high neutral hydrogen column density.
To further test the column density, in the gray shaded region in Figure 18, we normalize the stacked spectra to the mean continuum level within 960 ± 1Å in the absorber's rest frame so that the flux decrements in the Lyman continuum (< 912Å) can be revealed. This part of the stacked spectra comes from absorbers with z abs 2.95. Similar to that seen in the spectra of Lyman-limit systems (e.g., Fig.8 in Fumagalli et al. 2020), we see marked flux decrements blueward of 912Å. If the decrement were caused by a single absorber, the implied log N HI would be about 16.9, 17.2, 17.5, and 17.5 for the original sample and the subsample-1/2/3, respectively, estimated using the decrement around 912Å. We note that a similar exercise in Pieri et al. (2014) leads to an estimate of log N HI ∼ 16.7 (16.4) for absorbers with flux transmission −0.05 < F < 0.05 (0.05 < F < 0.15), a similar level to that inferred for our original sample. This is argued in Pieri et al. (2014) as an upper limit because of effects like the large-scale excess in forest absorption associated with the selected absorbers. Anyway, the increasing decrement over these samples leads support to the scenario that high-column-density systems contribute to the −1 < δ Lyα < −0.9 absorbers selected in this work. A more realistic model, which is beyond the scope of this work, should consist of absorbers with a distribution of column densities, and the decrement is related to the average absorption exp(−τ ) over such a distribution of absorbers. Without such a model, we could not assess the fraction of Lyman-limit systems and the relative importance of the high-column-density systems and the strongly-clustered low-column-density systems.
Next we compare the metal species. Following the original sample and the subsample-1/2/3 sequence, we see that metal absorptions become increasingly stronger. In addition, a few species not clearly detected in the original sample start to emerge in the subsamples, which are labelled in Figure 18 Figure 17. Redshift distributions for Lyα absorbers in the original sample (black) with −1.0 < δLyα < −0.9 and in the subsamples (blue/green/red) with the 2/4/6 neighboring pixels around the selected absorbers in the original sample satisfying δLyα < −0.9 + σ δ , with σ δ the uncertainty in δLyα.
The number distributions are shown in the left panel, and the normalized distributions (i.e., the probability distribution functions; PDFs) are in the right panel.
The measured column densities of the metal species for the subsamples are shown in Figure 19, along with those for the original sample. From the original sample to subsample-1/2/3, the column densities of lowionization species (like O I and Fe II) in general increase by more than one order of magnitude, a much stronger change than high-ionization species (e.g., O VI, C IV and Si IV; with N V being an exception caused by blending with the Lyα shadow lines from the contamination of Si II 1190/1193Å in Lyα as shown in Fig. 15). Such a relative change between column densities of lowand high-ionization species is a typical feature in highcolumn-density systems (see the bottom two panels of Fig. 14) and is seen in the observational study with high-resolution quasar spectra (e.g., Lehner et al. 2021). Figure 12, in the top panel of Figure 19, we show a three-component Cloudy model (not a fit) to illustrate the role of adding a high-column-density component in reproducing the trend in the low-ionization species. In the bottom panel, we boost the contribution of the high-column-density system, which provides a reasonable match to the column densities of the lowionization species in subsample-3, the one presumably with the highest contribution from high-column-density systems. While the fractions in the illustration are not to be taken seriously, the numbers are self-consistent: the ∼65% high-column-density systems in subsample-3 would make ∼2% of the original sample, based on the ratio of the numbers of absorbers in the two samples (Fig. 17), and this is treated as a lower limit of the contribution to the original sample.

Similar to
To summarize, the relative change in the column densities of low-and high-ionization metal species clearly favor the scenario that from the original sample to the subsample-1/2/3 the average neutral hydrogen column density of the absorbers increases (with possible modulations by metallicity), instead of a pure effect of absorber clustering. The change in the Lyman series absorption profiles and the decrement at the Lyman limit also support such a scenario. Future investigations based on high-resolution observations and hydrodynamic simulations will help verify and quantify the finding here.  Figure 18. Similar to Figure 3, but with the stacked spectra of the original sample (black) and three subsamples (blue/green/red). The original sample consists of absorbers with −1.0 < δLyα < −0.9 in the redshift range of 2 < z abs < 4. The subsamples (blue/green/red) are from selecting stronger absorbers by requiring the δLyα values of the neighboring 2/4/6 pixels around the absorbers in the original sample to lie below −0.9 + σ δ , where σ δ is the uncertainty in δLyα. In the shaded region, the flux is normalized to the mean continuum level within 960 ± 1Å to show the absorption near the Lyman break (912Å). Table 2. Measurements of metal lines in the 2.0 < z abs < 4.0 stacked spectrum. Column density N and Doppler parameter b from Voigt profile fitting are in units of cm −2 and km s −1 , respectively. The complete set of metal line measurements is available in the online Journal in a machine-readable format.
Species λ (Å) −1.0 ≤ δLyα < −0.9 −0.9 ≤ δLyα < −0.8 −0.8 ≤ δLyα < −0.7 −0.7 ≤ δLyα < −0.6  Combination ( 3% , 32% , 65% ) Figure 19. Column densities of metal species in the original sample (black symbols) and in the subsamples (blue/green/red symbols). The color of symbols follow the same notation as in Figure 18. In each panel, three Cloudy model curves are shown -a component with a low hydrogen number density and high temperature (red dotted line), a component with a medium hydrogen number density (blue dotted line), and a component with a high hydrogen number density and neutral hydrogen column density (black dotted line). Different combinations (not model fits) of the three components are shown as the solid black (top panel) and the solid red (bottom panel) lines to illustrate the difference between the absorbers in the original sample and those in the subsample of the strongest absorption. See the text for detail.