Impact of Systematic Redshift Errors on the Cross-correlation of the Lyman-$\alpha$ Forest with Quasars at Small Scales Using DESI Early Data

The quasar catalog was constructed following the logic in section 6.2 of [29], specifically figure 9. We briefly describe the method here.

The quasar catalog is produced using Redrock [24] and the two afterburner algorithms: a broad Magnesium II (MgII) line finder and a machine-learning based classifier QuasarNET [30, 31, 32]. Objects are first run through Redrock to determine both the class of the object (galaxy, quasar, star) and the best-fit redshift. Objects that are classified galaxy are then run through the MgII algorithm. If the MgII algorithm finds significant broad MgII emission, then their classification is changed to QSO and the redshift from Redrock remains unchanged. Any remaining objects are then run through QuasarNET. If QuasarNET classifies an object as QSO, Redrock is re-run for that object with the QuasarNET redshift as a prior. If neither of the three algorithms classify an object as QSO, then it is not included in the catalog. While the redshifts from QuasarNET are included in the catalog (as value-added content), the official redshift for each object comes from Redrock.

To finalize our catalog, we select only the objects that were targeted as quasars, using the classifications provided by the DESI target-selection pipeline [33]. Some of these objects were observed multiple times during EDR and M2 and therefore have multiple entries, but we want to keep only one observation per object in the quasar catalog. To select the best observation for the duplicated quasars, we first remove those that have a redshift warning flags set (i.e, those with a bad fit). If there are still repeat observations of the same quasar at this point, we select the observation with the highest TSNR2_LYA value. TSNR2_LYA, which is defined in section 4.14 of [23], is a measure related to the signal to noise ratio in the blue part of spectra for a given observation. We also remove any BAL quasars that were identified by the BALFinder [34] as having their absorption index (AI) > 0.

After these cuts, a total of 290,506 quasars are present in the quasar catalog. However, only the 106,861 quasars with redshift z > 1.88 will be included in the cross-correlation discussed in section 3.2. We refer to these quasars as tracer quasars or the tracer quasar catalog. Figure 1 shows the redshift distribution of the entire quasar sample with the dashed grey line representing this cutoff for the tracer quasar sample. The red hatched histogram shows the redshift distribution for the Lyman-$\alpha$ sample discussed in the next section.

We use the same Lyman-$\alpha$ forest sample as described in [35]. We will briefly describe the sample here, but refer the reader to [35] for a more detailed description.

We define the Lyman-$\alpha$ forest as the region between the Lyman-$\beta$ and Lyman-$\alpha$ emission lines at 1025.7 Å and 1215.67 Å, respectively. To reduce any contamination from the wings of those emission lines, we reduce the rest-frame wavelength range of all the Lyman-$\alpha$ forests to 1040 Å < $\lambda$ < 1205 Å. An example quasar spectrum highlighting the Lyman-$\alpha$ forest is shown in figure 2. Other features highlighted are the CIV and semi-forbidden CIII emission lines. The region between these lines shaded in orange is used during the calibration step of the continuum fitting process discussed in section 3.1 and in [35].

The Lyman-$\alpha$ forest sample is created from the same quasar catalog described in section 2.1, but the BAL quasars are not excluded. We set a limit on the observed wavelength range of the Lyman-$\alpha$ forest between 3600 Å and 5772 Å. The lower limit is set to the lowest wavelength measured by the blue DESI spectrographs [23] and corresponds to a minimum Lyman-$\alpha$ pixel redshift of $z=1.96$, meaning the lowest redshift a quasar can have is $z=2.02$. The upper limit is set by the small number of high redshift quasars in EDR+M2. This limit corresponds to a maximum quasar redshift of $z$ = 4.29 (with a corresponding maximum Lyman-$\alpha$ pixel redshift of $z=3.75$).

With these redshift cuts, there are 109,900 valid Lyman-$\alpha$ forests from the quasar catalog. Prior to fitting the continuum (described in section 3.1), we first mask regions of the spectra that are contaminated. This contamination can come from sky lines, galactic absorption, BALs, and Damped Lyman Alpha (DLA) regions. We mask any DLA regions identified with the DLAFinder [36] and BAL regions identified with the BALFinder [34] of the quasar spectra. For the sky lines and galactic absorption, we mask four areas of the spectra: sky lines from [5570.5, 5586.7] Å; and the Ca H, and K lines from the Milky Way from [3967.3, 3971.0] Å and [3933.0, 3935.8] Å, respectively. We also perform a calibration using the CIII between 1600 < $\lambda_{\text{RF}}$ [Å] < 1850, which is shaded in orange in figure 2. This calibration step follows the same process described in section 3.1 and in [35]. We require each forest to have a minimum length of 40 Å. Forests that are too short (don’t have enough pixels), have low signal (< 0), or have problems fitting the continuum (negative continuum) will be rejected during the calibration and continuum fitting process (section 3.1). After these rejections, the final Lyman-$\alpha$ sample has 88,511 forests. Figure 1 shows the redshift distribution of Lyman-$\alpha$ pixels in the sample.

We follow the method in [35] to measure the flux transmission field for the Lyman-$\alpha$ sample. We will briefly describe the method here but refer the reader to [35] for a more detailed description.

In this work, we adopt the same 0.8 Å linear pixel size used by the DESI spectroscopic pipeline [23] and other Lyman-$\alpha$ forest analyses using DESI EDR+M2 data [14, 35]. Previous analyses in eBOSS used logarithmic pixels [8]. We follow the method described in [35].

For each forest in each line of sight, $q$, we calculate the flux transmission field $\delta_{\text{q}}$ as a function of observed wavelength:

Here, $f_{\text{q}}(\lambda)$ is the observed flux, $C_{\text{q}}(\lambda)$ is the unabsorbed quasar continuum, and $\overline{F}(\lambda)$ is the mean transmitted flux at a specific wavelength. To measure each delta field, we estimate the mean expected flux, $\overline{F}(\lambda)C_{\text{q}}(\lambda)$, which is the product of the terms in the denominator in eq. 3.1. We assume it can be expressed according to:

where $\Lambda\equiv\log\lambda$, $\overline{C}(\lambda_{\text{RF}})$ is the estimate of the mean continuum for the considered quasar sample, and $a_{\text{q}}$, $b_{\text{q}}$ are per-quasar parameters.

The estimated variance of the flux $\sigma_{\text{q}}^{2}$ depends on the noise from the DESI spectroscopic pipeline and the intrinsic variance of the Lyman-$\alpha$ forest ($\sigma_{\text{LSS}}^{2}$) and we model it as

where $\tilde{\sigma}_{\text{pip,q}}=\sigma_{\text{pip,q}}/\overline{F}C_{\text{q}}(\lambda)$ and $\sigma_{\text{pip,q}}$ is the spectroscopic pipeline estimate of the flux uncertainties. It is further modified by the correction, $\eta(\lambda)$, which accounts for inaccuracies in the spectroscopic pipeline noise estimation. In DESI $\eta(\lambda)$ ranges between 0.97 and 1.05 between 3600 Å and 5772 Å, respectively, and is generally close to 1 due to first performing a calibration using the CIII region. The mean continuum $\overline{C}(\lambda_{\text{RF}})$, the function $\eta(\lambda)$, $\sigma_{\text{LSS}}(\lambda)$, and the quasar parameters, $a_{\text{q}}$ and $b_{\text{q}}$, are fit in an iterative process and stable values are obtained after about 5 iterations.

In this work we study the Lyman-$\alpha$ forest-quasar cross-correlation function. We require the use of the Lyman-$\alpha$ auto-correlation function to constrain some of the fit parameters (discussed in section 5). We use the Lyman-$\alpha$ forest auto-correlation results from [14] and, when necessary (in sections 4 and 5.3), run the auto-correlation following the method in [14]. The results of the cross-correlation are presented in section 5.1.

The separation between two tracers is determined both parallel and perpendicular to the line of sight, with $\vec{r}$ = ($r_{\parallel}$, $r_{\perp}$). For two tracers i and j with redshift $z_{\text{i}}$ and $z_{\text{j}}$, that are offset by an observed opening angle $\Delta\theta$, the comoving separation can be computed as

where $D_{\text{c}}(z)=c\int_{0}^{z}dz/H(z)$ is the comoving distance, $D_{\text{m}}$ is the angular comoving distance, and $H(z)$ is the Hubble function. To use the comoving distance we have to assume a fiducial cosmology and we adopt the flat $\Lambda$CDM cosmology of Planck Collaboration et al. (2018) [37]. With this cosmology we have $D_{\text{c}}=D_{\text{m}}$. In the above equations, the two tracers are either two Lyman-$\alpha$ pixels for the auto-correlation or a Lyman-$\alpha$ pixel and a quasar for the cross-correlation. The redshift of a Lyman-$\alpha$ pixel is calculated as $z_{\text{i}}=\lambda_{\text{obs}}/\lambda_{\text{Ly$\alpha$}}-1$, where $\lambda_{\text{Ly$\alpha$}}=1215.67$ Å and the quasar redshift is given from the tracer quasar catalog described in section 2.1.

The cross-correlation of the Lyman-$\alpha$ forest with quasars is defined as

where $\delta_{\text{Q}}(\vec{x})$ is the fluctuation of the number density in quasars, and $\vec{r}=(r_{\parallel},r_{\perp})$ (with $r^{2}=r_{\parallel}^{2}+r_{\perp}^{2}$) for bin A. If we assume that quasars are shot-noise dominated, then the cross-correlation is just the average Lyman-$\alpha$ transmission around quasars at $\vec{x}_{\text{Q}}$ [38].

For the Lyman-$\alpha$-quasar cross-correlation we use the same estimator as previous studies [8, 14, 39]:

where $i$ indexes a Lyman-$\alpha$ pixel and $j$ a quasar, and $\delta_{\text{i}}$ is the flux transmission for the Lyman-$\alpha$ pixel defined in eq. 3.1. The weights $w$ are defined in [14] and proportional to $(1+z)^{\gamma-1}$, where $\gamma=1.44$ for quasars and 2.9 for a Lyman-$\alpha$ pixel. The sums are over all quasar-pixel pairs with separations that are within the bounds of bin A, with the exception of Lyman-$\alpha$ pixels from their own background quasar which are excluded. The bounds of bin A are defined in comoving coordinates and are determined based on the selected bin width.

The correlation is binned along and across the line of sight, but the cross-correlation itself is not symmetric (i.e, you can’t swap indices $i$ and $j$ for $\delta_{i}$ in eq. 3.6). We define a positive separation along the line of sight when the quasar is in front of the Lyman-$\alpha$ pixel ($z_{qso}<z_{Ly\alpha}$) and likewise define a negative separation along the line of sight when the quasar is behind the Lyman-$\alpha$ pixel ($z_{qso}>z_{Ly\alpha}$). This asymmetry in the cross-correlation is what allows us to study systematic quasar redshift errors.

In this work we are studying small-scale effects on the cross-correlation, so we do not need to extend to large values of $r_{\parallel}$ and $r_{\perp}$, unlike the studies of BAO and the full-scale correlations in [14, 8]. We instead use a maximum separation of $80\,h^{-1}\,\mathrm{Mpc}$ in both $r_{\parallel}$ and $r_{\perp}$. For the cross-correlation, this corresponds to a minimum tracer quasar redshift of $z=1.88$ given the minimum Lyman-$\alpha$ pixel redshift of $z=1.96$ given in section 2.2.

Since our maximum separation in this study is $80\,h^{-1}\,\mathrm{Mpc}$, for the cross-correlation we restrict $r_{\parallel}$ $\in[-80,80]$ $h^{-1}$ Mpc and $r_{\perp}$ $\in[0,80]$ $h^{-1}$ Mpc. Previous measurements of 3D correlations to study BAO with eBOSS or DESI [14, 8] used a bin width of $4\,h^{-1}\,\mathrm{Mpc}$. Here, the focus of this paper is to measure a small-scale effect, we use a bin width of $1\,h^{-1}\,\mathrm{Mpc}$, along both transverse and longitudinal separations. This gives us 80 bins across the line of sight, 160 bins along the line of sight, and N = 80 $\times$ 160 = 12800 bins in total. We validate our choice to use a $1\,h^{-1}\,\mathrm{Mpc}$ bin width with mock data in section 4.

We follow the same method as in [14] to estimate the covariance matrix for the cross-correlation. The covariance matrix is estimated by splitting the sky into $s$ regions (according to HEALPix [40] pixels) and then computing their weighted covariance. For our analysis we use 440 HEALPix pixels with nside = 16. Assuming that correlations between regions are negligible, the covariance is written as

where A and B are two different bins of the correlation function $\xi$, $W=\sum_{\text{s}}W^{\text{s}}$, $W^{\text{s}}$ is the summed weight for one of the regions $s$, and $\xi_{\text{A}}=\sum_{\text{s}}W_{\text{A}}^{\text{s}}\xi_{\text{A}}^{\text{s}}/W_{\text{A}}$.

The cross-covariance matrix has N = $12800\times 12800=163$ million entries and the auto-covariance (from [14]) has N = $2500\times 2500=6.25$ million entries. The covariance matrices are dominated by the diagonal elements of the matrix. The estimates of off-diagonal elements are noisy, so they are modeled as a function of the difference of separation along the line of sight ($\Delta r_{\parallel}=|r^{\text{A}}_{\parallel}-r^{\text{B}}_{\parallel}|$) and across the line of sight ($\Delta r_{\perp}=|r^{\text{A}}_{\perp}-r^{\text{B}}_{\perp}|$) and then smoothed [8]. This smoothing is necessary to obtain a positive-definite, or invertible, covariance matrix. Correlation coefficients that have the same $\Delta r_{\parallel},\Delta r_{\perp}$ are averaged.

Using the Lyman-$\alpha$ forest when estimating the quasar continuum biases the estimates of the delta field in eq. 3.1 because the measured delta for a given pixel is a combination of the absorption signals located at other pixels. This produces a distortion in the correlation functions that occurs throughout, but mainly appear for small $r_{\perp}$. The biases arise because

1. 1.

   fitting the quasar parameters $a_{\text{q}}$ and $b_{\text{q}}$ biases the mean $\delta_{\text{q}}$ toward 0, and

2. 2.

   fitting the mean transmitted flux, $\overline{F}$, biases the mean $\delta$ at each observed wavelength, $\overline{\delta(\lambda)}$, toward 0.

These effects can be modeled by transformations to $\delta_{\text{q}}(\lambda)$. With the assumption that these transformations to $\delta_{\text{q}}$ are linear, the correlation function is distorted by a distortion matrix, $D_{\text{AA'}}$, such that

For the cross-correlation, the distortion matrix is calculated as

where $\eta$ is the projection due to the distortion that occurs during the continuum fitting process and is defined in equation 3.3 of [14], and the weights $w$ are defined in eq. 3.6.

When fitting the data, (described in section 3.4), the physical model of the correlation function is multiplied by the distortion matrix.

This section presents the theoretical model for the cross-correlation. Though this work focuses on the cross-correlation, the auto-correlation is necessary to include during the fitting process as it helps break degeneracies between parameters. We follow the method from [14] for fitting the auto-correlation, and use the auto-correlations from [14] in the baseline analysis.

In this model, the expected measured (distorted) cross-correlation in the ($r_{\parallel}$, $r_{\perp}$) bin $A$ is related to the theoretical (true) correlation by the distortion matrix as shown in eq. 3.8. The theoretical cross-correlation, $\xi^{\text{qf,th}}$, can be broken down into the sum of its components:

Here, the first term gives the contribution from the cross-correlation between the Lyman-$\alpha$ forest and quasars, the second term gives the contribution from the cross-correlation between quasars and other absorbers, like metals and high column density systems (HCDs), in the Lyman-$\alpha$ region, and the last term gives the contributions from the transverse proximity (TP) effect, which is the effect of radiation from the quasar on the nearby surrounding gas. We describe these terms in depth throughout this section and the parameters of this model are described in table 5.

The Lyman-$\alpha$-QSO contribution to the correlation function is derived from the Fourier transform of the tracer biased power spectrum, $P(\textbf{k},z)$, which is written as

where k is the wavenumber with modulus $k$ and $\mu_{\text{k}}=k_{\parallel}/k$. The bias and linear redshift-space distortion parameters, $b_{\text{i,j}}(z)$ and $\beta_{\text{i,j}}$, are for tracer $i$ or $j$. $G$ is a correction that accounts for the averaging of the correlation function binning and $F_{\text{NL}}$ accounts for the non-linear effects on small scales (large k). $P_{\text{QL}}$ is the quasi-linear matter power spectrum defined in eq. 3.16. For the auto-correlation, $i=j$ and the only tracer is Lyman-$\alpha$ absorption. In this work, we are concerned with the cross-correlation, $i\neq j$ and the two tracers are Lyman-$\alpha$ absorption and quasars.

In eq. 3.11, the bias $b$ and redshift-space distortion $\beta$ parameters for each tracer $i,j$ appear with the standard Kaiser factor: $b_{\text{i,j}}(z)(1+\beta_{\text{i,j}}\mu_{\text{k}}^{2})$ [41]. The fit of the cross-correlation is only sensitive to the product of the quasar and Lyman-$\alpha$ biases. The quasar bias $b_{\text{q}}$ is redshift dependent and is given by

where $z_{\text{eff}}$ is the effective redshift of the sample and $\gamma_{\text{q}}$ = 1.44 [42]. Following the analysis of [14] we allow $b_{\text{q}}$ to remain free in the fit. On the other hand, the quasar redshift-space distortion (RSD) parameter, $\beta_{\text{q}}$, is not a free parameter. It is related to the quasar bias $b_{\text{q}}$ by:

where $f=0.9703$ is the linear growth rate of structure in our fiducial cosmology. We then derive the value of $\beta_{\text{q}}$ from the result of $b_{\text{q}}$.

We assume that the Lyman-$\alpha$ bias parameter, $b_{\text{Ly$\alpha$}}$, is redshift dependent but use the approximation that the Lyman-$\alpha$ RSD parameter, $\beta_{\text{Ly$\alpha$}}$, is redshift independent. Following [14], $b_{\text{Ly$\alpha$}}$ is given by

with $\gamma_{\text{Ly$\alpha$}}$ = 2.9 [43]. We allow $b_{\text{Ly$\alpha$}}$ and $\beta_{\text{Ly$\alpha$}}$ to remain free in the fit.

The Lyman-$\alpha$ absorption contains contributions from both the inter-galactic medium (IGM) and high-column-density (HCD) systems. HCD absorbers, such as Damped Lyman Alpha, trace the underlying density field. When these HCD systems are correctly identified and the appropriate parts of the absorption region are properly masked and modeled, they will have no effect on the measured correlation functions. However, if these systems are left unidentified, they will affect the measured correlation function as a broadening/smearing in the radial direction. Following the method in [44], this effect introduces a $k_{\parallel}$ dependence to the effective bias and can be modeled as

where $(b_{\text{Ly}\alpha},\beta_{\text{Ly}\alpha})$ and $(b_{\text{HCD}},\beta_{\text{HCD}})$ are the bias and redshift space distortion parameters associated with the IGM and HCD systems. Following [8], we use $F_{\text{HCD}}=\exp(-L_{\text{HCD}}k_{\parallel})$, which is approximated from the results in [45]. $L_{\text{HCD}}$ is the length scale for unmasked HCDs, and is degenerate with other parameters. For this work, we fix $L_{\text{HCD}}=$ $10\,h^{-1}\,\mathrm{Mpc}$ and allow $b_{\text{HCD}}$ to be free in the fit. We fix $\beta_{\text{HCD}}=0.5$.

While this work focuses on comoving separations smaller than the BAO scale, for consistency with other analyses [14, 46], we model the quasi-linear matter power spectrum, $P_{\text{QL}}(\textbf{k},z)$, as the sum of a smooth and peak components with an empirical anisotropic damping applied to the BAO peak component:

$P_{\text{sm}}$ in eq. 3.16 is derived from the linear power spectrum $P_{\text{L}}(k,z)$ via the side-band technique described in [47]. The redshift dependent linear power spectrum, $P_{\text{L}}$, is derived from CAMB [48], and we then define $P_{\text{peak}}=P_{\text{L}}-P_{\text{sm}}$. The BAO peak is affected by non-linear broadening [49] and is corrected for by introducing the parameters $(\Sigma_{\parallel},\Sigma_{\perp})$. $\Sigma_{\parallel}$ is related to $\Sigma_{\perp}$ and the growth rate such that

where $f=0.9703$ is the linear growth rate of structure in our fiducial cosmology. Given that we do not fit the BAO peak in this work, we fix ($A_{\text{peak}},~{}\Sigma_{\parallel},~{}\Sigma_{\perp}$) = (1, $6.37\,h^{-1}\,\mathrm{Mpc}$, $3.24\,h^{-1}\,\mathrm{Mpc}$) [14].

$F_{\text{NL}}$ in equation 3.11 corrects for additional effects at small scales. In the auto-correlation these are thermal broadening, peculiar velocities, and nonlinear growth structure, parameterized according to equation 3.6 in [50]. For the cross-correlation the most important small-scale correction is due to quasar velocities and the precision of quasar redshift measurements. Following [14, 51], we model this using the Lorentz-damping form:

where $\sigma_{\text{v}}$ is a free parameter that describes the precision of quasar redshift measurements. Alternative models for $F_{\text{NL}}$ use a Gaussian form.

The last term in eq. 3.11, $G(\textbf{k})$, accounts for the effects of the binning of the correlation function on the $(r_{\parallel},r_{\perp})$ separation grid. Assuming the distribution in each bin is homogeneous, and following the method used in [52],

where $R_{\parallel}$ and $R_{\perp}$ are the scales of the smoothing, i.e, are the radial and transverse widths of the bins, respectively. In this work we focus on bin widths of $1\,h^{-1}\,\mathrm{Mpc}$ for both $R_{\parallel}$ and $R_{\perp}$, but validate the choice between $1\,h^{-1}\,\mathrm{Mpc}$ and $4\,h^{-1}\,\mathrm{Mpc}$ in section 4.

The next component that contributes to the model cross-correlation in eq. 3.10 is the sum over the non-Lyman-$\alpha$ absorbers, or metals. The power spectrum for the cross-correlation with metals has the same form as that used for the Lyman-$\alpha$-quasar cross-correlation, though it is simplified by neglecting the effect of HCDs. Because there is little absorption by metals, it is further simplified by not separating the smooth and peak components. It is however, more complicated due to the fact that the ($r_{\parallel}$, $r_{\perp}$) bins in the correlation correspond to an observed $(\Delta\theta,\Delta\lambda)$ calculated assuming that the absorption is due to the Lyman-$\alpha$ transition.

This causes a shift in the model correlation function, and the contribution of each absorber is maximized in the ($r_{\parallel}$, $r_{\perp}$) bin that corresponds to vanishing physical separation. For Lyman-$\alpha$ absorption this corresponds to ($r_{\parallel}$, $r_{\perp}$) = (0,0). For the other absorbers, this maximum occurs at $r_{\perp}$ = 0 and $r_{\parallel}\approx\frac{c}{H(z)}(1+z)(\lambda_{\text{m}}-\lambda_{\text{Ly}\alpha})/\overline{\lambda}$, where $\overline{\lambda}$ is the mean value between the $\text{Ly}\alpha$ and metal absorption [14]. The values of the $r_{\parallel}$ separations for the metal absorbers considered in this work are given in table 1.

The shifted-model correlation function is calculated with respect to the unshifted-model correlation function, $\xi^{mn}$, for each absorber pair $(m,n)$ by introducing a metal matrix $M_{\text{AB}}$ as in [53]:

where $M_{\text{AB}}$, the metal matrix, is defined as

Here, $W_{\text{A}}=\sum_{\text{(m,n)$\in$A}}w_{\text{m}}w_{\text{n}}$ and $w_{\text{m}},w_{\text{n}}$ are the weights for each absorber and are redshift dependent such that $w\propto(1+z)^{\gamma}$. $(m,n)\in A$ refers to separations calculated using the reconstructed redshift of the absorber pair and $(m,n)\in B$ refers to separations computed using the redshifts $z_{\text{m}}$ and $z_{\text{n}}$ of the absorber pair $(m,n)$.

In the fits to the data each metal species has its own individual bias parameters $(b,\beta)$. Since the correlations are only visible in ($r_{\parallel}$, $r_{\perp}$) bins that correspond to small physical separations and the amplitudes for the SiII and SiIII metal species are determined by the excess correlations, we do not have enough signal to determine the $(b,\beta)$ parameters separately. As in previous works ([8, 51]), we fix $\beta=0.5$ for all metal species. Since we are using the auto-correlation from [14] we include correlations from the four metal species given in table 1 in the correlations.

The last term in eq. 3.10 is due to the transverse proximity effect. In the area surrounding a quasar, there is significantly less Lyman-$\alpha$ absorption because the radiation emitted from the quasar dominates over the UV background and increases the ionization fraction of the surrounding gas, thus making it more transparent to Lyman-$\alpha$ photons. The Lyman-$\alpha$ absorption of a background quasar will decrease in strength at the redshift of the foreground quasar when the Lyman-$\alpha$ forest is close to the quasar line of sight. This will therefore affect the correlation between quasars and the Lyman-$\alpha$ forest. As in [14], we model this effect in the correlation with

where $r$ is the comoving separation in units of  $h^{-1}$ Mpc, $\lambda_{\text{UV}}=$$300\,h^{-1}\,\mathrm{Mpc}$ [54], and $\xi_{0}^{TP}$ is an amplitude that will be fit.

Systematic errors on the measurement of the quasar redshift can cause a shift in the model estimation of the correlation function in $r_{\parallel}$. To account for this, we allow for a shift of the absorber-quasar separation in the $r_{\parallel}$ direction

The shift is mostly due to systematic quasar redshift errors. However, any asymmetries in $r_{\parallel}$ will also affect $\Delta r_{\parallel}$. The model of the cross-correlation function is not symmetric about $r_{\parallel}$ due to the contribution of metal absorptions and the variation of mean redshift with $r_{\parallel}$. The continuum-fitting distortion also introduces an asymmetry in $r_{\parallel}$. [55] has showed that while $\Delta r_{\parallel}$ is sensitive to relativistic effects, they are small. We do not model any relativistic effects in this work as [56] showed that they are partially degenerate with $\Delta r_{\parallel}$.

The measured Lyman-$\alpha$ forest and quasar cross-correlation and the result for the baseline fit are presented in figure 5. We show the average of the correlation and model in each $r_{\perp}$ slice listed. The correlations have a lower amplitude for larger $r_{\perp}$ slices and therefore appear noisier. The model for the baseline fit is calculated jointly using the auto- and cross-correlations and has 13 free parameters. The values for the fixed parameters used in this model are presented in table 5 in Appendix A.

We present the results of the baseline fit in the first column in table 2. The baseline fit has a joint $\chi^{2}$ value that corresponds to essentially zero probability. This is likely due to our model failing at small scale separations. We re-fit our model, increasing the minimum separation to $5\,h^{-1}\,\mathrm{Mpc}$ and $10\,h^{-1}\,\mathrm{Mpc}$ to check the results of the $\chi^{2}$ values. We found a very minimal improvement in $\chi^{2}$ values and a worsening of $\Delta r_{\parallel}$ values. We therefore decide to keep the minimum separation in our fit at $0\,h^{-1}\,\mathrm{Mpc}$ and leave for future studies any further tests of the model on $\chi^{2}$ values.

We find $\Delta r_{\parallel}$ $=-1.94\pm 0.15$ $h^{-1}$ Mpc and $\sigma_{\text{v}}=5.8\pm 0.5$ $h^{-1}$ Mpc. Comparing these results to table 1 of [14], we find that $\sigma_{\text{v}}$ is within 1.0$\sigma$ from the value reported in [14]. Our result for $\Delta r_{\parallel}$ is also consistent within 1.5$\sigma$ with [14], while our fit value for $\xi_{0}^{\text{TP}}$ is 2.0$\sigma$ with respect to the results reported in [14]. These differences are likely due to our smaller bin width, and the different ranges in $r_{\parallel}$ and $r_{\perp}$ used.

A negative value for $\Delta r_{\parallel}$ means the measured cross-correlation is shifted in the positive $r_{\parallel}$ direction, which is seen in all 9 $r_{\perp}$ slices in figure 5. Based on our definition of a positive separation, i.e, when the quasar is in front of the Lyman-$\alpha$ pixel, this indicates that the estimated quasar redshifts are less than the true quasar redshifts.

Our initial thought for this discrepancy is the disappearance of the MgII emission line from quasar spectra above redshift $z\approx 2.5$. The MgII line, a reliable indicator of systemic redshift [60, 39], is located beyond the observable range of the DESI spectrographs at these redshifts. We therefore rely on the high ionization broad lines, which often display significant velocity shifting [21], to obtain redshift estimates at $z>2.5$. This is also suggested by [29], who report a “kink" at $z\sim 2.5$ when comparing the redshift estimates from DESI to those from SDSS (their figure 16). However, this is unlikely the cause of the “kink" due to the observable range of the BOSS spectrographs being very similar to DESI [61], as eBOSS would lose access to the MgII line at similar redshifts to DESI.

The apparent underestimation of DESI redshifts persisted when updating the quasar templates in Redrock for processing the DESI year one quasar sample, as they did not account for the Lyman-$\alpha$ optical depth [27]. In section 6, we discuss the improvements in redshift performance for DESI year 1 data at $z>2$ when using an updated version of Redrock and quasar templates that properly model Lyman-$\alpha$ optical depth.

We test and study here the evolution and dependence of the parameters $\Delta r_{\parallel}$ and $\sigma_{\text{v}}$ on quasar redshift to investigate more in depth the cause of the measured biases. We split the tracer quasar sample described in section 2 into four redshift bins: z1: $z\leq 2.2$, z2: $2.2<z\leq 2.5$, z3: $2.5<z\leq 2.8$, and z4: $z>2.8$. As a reminder, quasars with redshift $z<1.88$ will not contribute to the correlation as the separation will be larger than our maximum separation of $80\,h^{-1}\,\mathrm{Mpc}$. We use the same sample of Lyman-$\alpha$ forests as the baseline, but we create new quasar catalogs with approximately 41,000 objects in the first catalog, 28,000 objects in the second, 17,000 objects in the third, and 19,000 objects in the fourth. We then re-calculate the cross-correlation and fit for each redshift bin following the methods in section 3. Since the auto-correlation only depends on the Lyman-$\alpha$ forest sample, we do not need to re-compute it and can use the auto-correlation from [14] for each redshift bin when computing the joint fit.

The fit results for each redshift bin are given in table 3. Each fit has the same configuration as the baseline fit presented in the first column of table 2. We plot these results vs the average redshift in each bin in figure 6, highlighting the baseline value in the shaded region. We find a strong evolution with redshift for $\Delta r_{\parallel}$ with the value becoming more negative as the average redshift of the bin increases for the first three bins. However, the fourth redshift bin does not clearly follow this trend. Though this bin does not have the least amount of quasars, the quasars in this bin are at a higher redshift and are generally fainter and have lower signal to noise. The fit for this redshift bin is also the worst of the four with a $\chi^{2}=1.053$, which corresponds to near-zero probability given the 11,631 degrees of freedom in the fit. Despite this, the overall trend indicates that the quasar redshifts are being underestimated by a larger amount as the quasar redshift increases, which is consistent with the quasar templates not accounting for the Lyman-$\alpha$ forest optical depth. If the loss of the MgII emission line in the quasar spectra was the dominant contribution to the redshift errors, we would expect a sudden jump in $\Delta r_{\parallel}$ and $\sigma_{\text{v}}$ as a function of $z$ around $z=2.5$. Such a feature is not observed, so the loss of MgII is unlikely to be the culprit.

The measured evolution of $\sigma_{\texttt{v}}$ as a function of $z$ is much less striking than that of $\Delta r_{\parallel}$. Still, there is an indication for a mild increase of $\sigma_{\texttt{v}}$ for higher redshift quasars, consistent with less precision in the measured redshift.

BAL systems can impact the shape of quasar emission lines and introduce shifts in the estimated quasar redshift, up to several hundred km s^-1 [19]. BAL quasars often have larger redshift errors compared to quasars with no BAL features because the BAL features impact the blue side of emission lines [18]. This is true in DESI, where Redrock tends to overestimate redshifts for BAL quasars [18]. The current pipeline used for most Lyman-$\alpha$ studies, Picca, has the option to mask out BAL regions during the continuum fitting process, but does not update the quasar redshift when these masks are applied so any systematic errors on the quasar redshift will still be present. [19] shows that when BAL regions of a quasar spectrum are masked and new redshifts are obtained, those new redshifts shift by 240 km s^-1 on average. We explore in this section the effect of including the BAL quasars in the tracer catalog as well as using the updated BAL redshifts from [19] once the BAL regions of the quasar spectra are masked.

We first test if including BAL quasars in the tracer quasar catalog will have an effect on $\Delta r_{\parallel}$ or $\sigma_{\text{v}}$. For this test we use the same Lyman-$\alpha$ forest sample as [35] which is the same as our baseline. The tracer quasar catalog is the same as described in section 2.1 except the BALs are not removed. It consists of 128,026 quasars with redshift $z>1.88$. We use the auto-correlation from [14] and follow the same method described in section 3 to measure and fit the cross-correlation. We find $\Delta r_{\parallel}$ = $-1.84\pm 0.13$  $h^{-1}$ Mpc and $\sigma_{\text{v}}=7.3\pm 0.6$  $h^{-1}$ Mpc. These results, along with those from the baseline configuration, are shown in figure 7 and presented in table 2 in the BAL column. The shift of $\Delta r_{\parallel}$ in the positive direction would seem to indicate that including BAL quasars in the tracer catalog reduces the quasar redshift errors. However, this is actually due to BAL quasars typically having overestimated redshifts, as the BAL features tend to impact the blue side of the emission lines [18]. Any overestimation of quasar redshifts would shift $\Delta r_{\parallel}$ in the positive direction. This is the opposite of what we observe when the BAL quasars are excluded, where the shift in $\Delta r_{\parallel}$ is negative due to the underestimation of redshifts. It is likely that the inclusion of the BAL quasars partially compensates for the negative shift.

To test if updating the BAL quasar redshifts after the BAL regions are masked improves the $\Delta r_{\parallel}$ parameter, we use the quasar catalog from [19]. This is the same catalog as described above, but with updated redshifts for the BAL quasars, of which 28,185 BAL quasars have an updated redshift. Because the redshifts are updated, the Lyman-$\alpha$ forest sample changes slightly, both in the number of forests and the redshift distribution of pixels. While the effect is small, we make the choice to not use the Lyman-$\alpha$ forest sample described in [35]. We follow the method in [35] and in section 3 for the continuum fitting process and our final Lyman-$\alpha$ sample consists of 88,432 forests. We use the tracer catalog with updated redshifts and follow the method described in section 3 to measure and fit the correlations. We can no longer use the auto-correlation from [14], so we measure those following the method in [14]. The results are shown in figure 7 and in table 2 in the ZBAL column. We find $\Delta r_{\parallel}$ = $-1.90\pm 0.11$  $h^{-1}$ Mpc and $\sigma_{\text{v}}=7.0\pm 0.6$  $h^{-1}$ Mpc. The value of $\Delta r_{\parallel}$ has shifted in the negative direction and is now equal ($0.35\sigma$) to the baseline value. This negative shift is consistent with what is found in [19], which found that the redshifts of BAL quasars after masking were shifted to the blue.

In both configurations when BAL quasars are added in at the catalog level the values of $b_{\text{Ly}\alpha}$ and $\beta_{\text{Ly}\alpha}$ change by a factor of two compared to the baseline. These parameters are highly correlated with each other and with $b_{\text{q}}$ ($\rho\sim 0.9$). By adding the BAL quasars at the catalog level we are changing the value of $b_{\text{q}}$, so it is not unexpected that $b_{\text{Ly}\alpha}$ and $\beta_{\text{Ly}\alpha}$ also change.

This shows that BAL quasars can affect the correlations when they are included in the tracer catalog and their redshifts are not updated after the BAL regions are masked. When the redshifts are updated the inclusion of BAL quasars in the tracer catalog has no effect on $\Delta r_{\parallel}$.