MEAN-FLUX-REGULATED PRINCIPAL COMPONENT ANALYSIS CONTINUUM FITTING OF SLOAN DIGITAL SKY SURVEY Lyα FOREST SPECTRA

In this paper, we carry out continuum fitting for publicly available spectra from the final SDSS Data Release 7 (DR7) quasar catalog (Schneider et al. 2010), which is comprised of 105,783 spectroscopically confirmed quasars observed from the 2.5 m SDSS Telescope in Apache Point, NM. The spectra cover the observed wavelength range λ_obs = 3800–9200 Å with a spectral resolution of R ≡ λ/Δλ ≈ 2000.

From this overall catalog, we select a subsample suitable for Lyα forest studies. First, we require that some portion of the quasar Lyα forest region, λ_rest = 1041–1185 Å, be within the observed wavelength range. Since the extreme blue end (near λ_obs ≈ 3800 Å) of the SDSS spectra are known to suffer from spectrophotometric problems, we use λ_obs = 3840 Å as the lower wavelength limit. This sets a minimum quasar redshift of z_QSO = 2.3. For the quasars that satisfy this redshift criterion, we excise the portions of the spectra that lie below λ_obs = 3840 Å. In addition, broad absorption line (BAL) quasars have continua that are difficult to characterize (although see Allen et al. 2011 for a method to recover quasar continua from BALs), therefore we discard quasars flagged as BALs in the Shen et al. (2011) value-added quasar catalog.

There are 13,133 quasars in the DR7 quasar catalog that satisfy the above criteria. We make further quality cuts by discarding 962 spectra that have SPPIXMASK = 0–12 bitmask set (this signifies issues with the fiber; see Stoughton et al. 2002 for further details on the SDSS bitmask system) and two spectra where the S/N was too low to normalize the spectra, leading to negative normalizations. This leaves us with a sample of 12,069 spectra to which we will apply the MF-PCA technique. Within individual spectra, we mask pixels which have either zero inverse variance or the SPPIXMASK = 16–28 bitmask set. This avoids the use of problematic pixels, such as rejected extractions, bright sky-lines, or bad flats.

The median S/N in the sample is S/N = 3.0 per 69  km s⁻¹ SDSS pixel within the Lyα forest, and S/N = 6.2 pixel⁻¹ in the λ_rest = 1225–1600 Å wavelength region. The redshift and S/N⁵ distributions of our final quasar sample is shown in Figure 1. It is clear that the Lyα forest data from DR7 are noisy. Most of the sightlines have median S/N < 10 within the Lyα forest, which is too noisy to be fitted individually using existing techniques.

Zoom In Zoom Out Reset image size

In addition, we need to deal with Damped Lyα Absorbers (DLAs) within the spectra. These are absorbing systems with neutral hydrogen column densities of $N_{{\rm H\,\mathsc{i}}} \ge 2 \times 10^{20}\; \mathrm{cm^{-2}}$ that result in complete absorption over large portions (Δv ∼ 10³ km s⁻¹) of affected sight lines. Since the MF-PCA technique (Section 3) fits the amplitude of the quasar continuum based on the mean flux of the low column-density Lyα forest, the excess absorption of a DLA within a sightline would bias the continuum estimate.

To correct for this, we use a catalog of 1427 DLAs identified in the SDSS DR7 spectra by Noterdaeme et al. (2009). First, we mask the wavelength region corresponding to the equivalent width of each DLA (Draine 2011):

$$\begin{equation} W \approx \lambda _{\alpha } \left[ \frac{e^2}{m_e c^2} N_{{\rm H\,\mathsc{i}}} f_{\alpha } \lambda _{\alpha } \left(\frac{\gamma _{\alpha } \lambda _{\alpha }}{c} \right) \right]^{1/2}, \end{equation} \tag{ 1 }$$

where λ_α = 1216 Å is the rest-frame wavelength of the hydrogen Lyα transition, e is the electron charge, m_e is the electron mass, c is the speed of light, $N_{{\rm H\,\mathsc{i}}}$ is the H i column density of the DLA, f_α is the Lyα oscillator strength, and γ_α is the sum of the Einstein A coefficients for the transition.

However, the damping wings of each DLA extend beyond the equivalent width, providing a small but non-negligible excess absorption to the pixels close to the DLA. We correct for this by multiplying each pixel in the spectrum with exp (τ_wing(Δλ)), where

$$\begin{equation} \tau _\mathrm{wing}(\Delta \lambda) = \frac{e^2}{m_e c^2} \frac{\gamma _{\alpha } \lambda _{\alpha }}{4 \pi } f_{\alpha } N_{{\rm H\,\mathsc{i}}} \lambda _{\alpha } \left(\frac{\lambda }{\Delta \lambda }\right)^2 \end{equation} \tag{ 2 }$$

and Δλ ≡ λ − λ_α is the wavelength separation in the DLA rest frame.

In order to predict the shape of the Lyα forest continuum, we need a set of template spectra with clearly identified continua at wavelengths λ_rest < 1216 Å. For this purpose, we will use two different sets of quasar templates, derived from quasars observed in the HST, and SDSS itself.

Suzuki et al. (2005) derived PCA templates from 50 quasars that had been observed by the Far Object Spectrograph (FOS) on the HST in the ultraviolet. At the low redshifts (0.14 < z_QSO < 1.04) of these quasars, the line density of the Lyα forest is sufficiently small that the quasar continuum could be clearly identified. This enabled the creation of templates in the range λ_rest = 1025–1600 Å, covering Lyβ λ1025 to C iv λ1549.

Pâris et al. (2011) recently carried out a similar study, applying the techniques in Suzuki et al. (2005) to a subsample of 78 SDSS DR7 quasars. These z_QSO ≈ 3 quasars were selected to have full coverage of the Lyα forest and relatively high S/N (≳  10 pixel⁻¹). The transmission peaks in the Lyα forest were hand fitted with a low-order spline function to provide a continuum estimate. PCA templates were then derived in the spectral range λ_rest = 1020–2000 Å, which included the C iii λ1906 line. While this process might give a biased continuum level due to the low resolution and S/N of the templates, it should provide a good description of the relative shape of the quasar continuum which is required for MF-PCA—the mean-flux regulation process is designed to correct for uncertainties in the overall continuum level arising from pure PCA fitting.

We do not expect the redshift differences between the template and the DR7 quasars to be a significant issue, even for the Suzuki et al. (2005) quasars (〈z_QSO〉 ≈ 0.6). This is because various studies (e.g., vanden Berk et al. 2004; Fan 2006) have suggested that there is little redshift evolution of quasar spectra. However, the shape of quasar spectra is known to have a significant luminosity dependence, such as the well-known Baldwin effect (Baldwin et al. 1978), which is the anti-correlation between the strength of the C iv λ1549 emission line and the quasar luminosity. It is therefore reasonable to suppose that there would be a significant difference in the spectral shapes represented by two templates and the overall SDSS sample. This is because the Suzuki et al. (2005) sample is comprised of relatively low-luminosity, nearby quasars as opposed to the Pâris et al. (2011) quasars, which were selected to have high S/N and are therefore a luminous subsample of the SDSS quasars.

In order to compare the relative luminosities, we calculate λL₁₂₈₀, the intrinsic monochromatic luminosity near λ_rest = 1280 Å, for the quasars in our SDSS DR7 sample as well as the two template samples. We assume a standard ΛCDM cosmology with h = 0.7, Ω_m = 0.28, and Ω_m + Ω_Λ = 1. The respective distributions of λL₁₂₈₀ are shown in Figure 2. The SDSS DR7 quasars have a typical luminosity of λL₁₂₈₀ ≈ 10^46.3 erg s⁻¹, while the Suzuki et al. (2005) and Pâris et al. (2011) quasars are about 0.5 dex fainter and brighter, respectively. However, the combined luminosity distributions of the Suzuki et al. (2005) and Pâris et al. (2011) template quasars significantly overlap the full range of SDSS DR7 quasars, which justifies the use of both templates in this paper.

Zoom In Zoom Out Reset image size

The basic concept of PCA is that a normalized quasar spectrum, f(λ), can be represented as

$$\begin{equation} f(\lambda) \approx \mu (\lambda) + \sum ^{m}_{j=1} c_{j} \xi _j(\lambda), \end{equation} \tag{ 3 }$$

where μ(λ) is the mean quasar spectrum, ξ_j(λ) is the jth principal component or "eigenspectrum," and c_j are the weights for an individual quasar. The formalism for deriving the eigenspectra and weights is described in Suzuki et al. (2005) and Pâris et al. (2011).

The standard PCA formalism for deriving the weights, c_j, does not take into account spectral noise, which renders it unsuitable for noisy SDSS spectra (see Figure 1(b)). Instead, we first carry out a least-squares fit to the red side of each spectrum using the full λ_rest ≈ 1000–1600 Å eigenspectra as a basis. Due to the correlation between the weak emission lines within the Lyα forest and in λ_rest ∼ 1300–1500 Å (Suzuki et al. 2005), we expect this to provide a reasonable prediction for the shape of the continuum.

As described in Section 2.2, we have two separate sets of PCA eigenspectra from Suzuki et al. (2005) and Pâris et al. (2011). In principle, one could combine the two sets of quasar templates to generate one set of eigenspectra that would encompass the diversity of both template samples. However, the template spectra from Pâris et al. (2011) are not available to us at the time of writing, therefore we will carry out our fitting procedure separately for the two sets of PCA eigenspectra. Suzuki et al. (2005) had found that out of their 10 principal component eigenspectra, only the first eight components appeared to describe physical features in the spectra, while the ninth and tenth components seemed to describe mostly noise. Therefore, we will use only eight components from each set of eigenspectra for our fits. In addition, for the sake of consistency we limit ourselves to the rest-wavelength range λ_rest = 1020–1600 Å of each eigenspectrum even though the Pâris et al. (2011) eigenspectra extend up to λ_rest = 2000 Å.

However, we have found that fitting the SDSS spectra with just the PCA weights c_j was insufficient to account for the large diversity of the sample. Therefore, we introduce two additional fit parameters: a power-law component, α_λ, and redshift-correction factor, c_z. The power-law component, α_λ, is necessary due to the large range of slopes found in the SDSS quasars. Even though the third through fifth principal components in the HST eigenspectra include the spectral slope, they also describe some emission-line features—the introduction of α_λ as a free parameter allows an additional degree of freedom and enables a better fit to the emission lines and slope simultaneously. Because of this degeneracy between the slopes within the eigenspectra and α_λ, we do not interpret the latter as the slope of the underlying quasar power-law continuum. The power-law parameter also helps account for low-order spectrophotometric errors as well as dust extinction in the spectrum.

The redshift-correction factor, c_z = λ^fit_rest/λ^pipe_rest, translates the spectrum along the wavelength axis with respect to the rest wavelength given by the pipeline redshift, λ^pipe_rest, to a best-fitting rest wavelength, λ^fit_rest. It is required as the SDSS pipeline redshifts are not completely accurate (see, e.g., Hewett & Wild 2010). However, due to the asymmetry and velocity shifting of quasar emission lines at different redshifts, we do not necessarily interpret c_z as a true redshift correction—it is merely an ad-hoc parameter to obtain the best-possible fit to the spectrum. The full list of free parameters for our continuum-fitting procedure is shown in Table 1 (a_MF and b_MF are free parameters for the mean-flux regulation step, described in Section 3.2).

We are now in a position to carry out the fitting procedure. First, the quasar spectrum is shifted to the quasar rest frame using the pipeline redshift and normalized at λ_rest = 1275–1285 Å. We then use the least-squares fitting routine MPFIT (Markwardt 2009) to find the best-fitting set of parameters, [c_z, α_λ, c_j], given the spectrum and its noise. The initial fits from this procedure is represented by the black dashed lines in the examples shown in Figure 3.

Zoom In Zoom Out Reset image size

However, while the intrinsic quasar spectrum is generally well defined redward of Lyα in the SDSS spectra, in many cases intervening metal absorption lines can be seen in the spectrum in the λ_rest ≈ 1216–1600 Å wavelength. To prevent these absorption features from biasing the PCA fitting, we carry out a simple iterative procedure to mask these absorption lines: using the continuum, C_init, obtained from the initial least-squares fit, we mask pixels in which f(λ) − C_init(λ) < −2.5 σ(λ), where f(λ) and σ(λ) are the observed spectrum and pipeline noise, respectively. We then make a new PCA fit redward of the Lyα emission line and repeat this process until the fit converges. In Figure 3, the final PCA fits, C_PCA, are shown as black solid lines while the masked pixels are denoted by crosses.

From Figure 3, we see that the least-squares PCA fitting procedure generally works well, even with noisy (S/N ∼ few) spectra. The fit to the Lyα λ1216 emission line is sometimes imperfect, but unsurprising considering that the fitted range (λ_rest = 1216–1600 Å) only takes partial account of the line. Comparing the initial (dashed line) and final (solid line) fits in Figure 3, we see that the red-side metal absorption lines usually have little effect on the fits, but in certain cases (e.g., Figure 3(b) and (c)) absorption line masking noticeably improves the fit.

We use the absolute flux error to quantify the goodness of the PCA fits redward of Lyα:

$$\begin{equation} |\delta F| = \frac{\int ^{\lambda _\mathrm{max}}_{\lambda _\mathrm{min}} \left| \frac{ C_\mathrm{PCA}(\lambda) - \tilde{f}(\lambda) }{\tilde{f}(\lambda) } \right| d\lambda }{\int ^{\lambda _\mathrm{max}}_{\lambda _\mathrm{min}} d \lambda }, \end{equation} \tag{ 4 }$$

where C_PCA(λ) is the fitted continuum and $\tilde{f}(\lambda)$ is the observed spectrum smoothed by a 15 pixel boxcar to avoid biasing |δF| in noisy spectra. λ_max = 1600 Å and λ_min = 1225 Å represent the range over which we calculate |δF|.

The distribution of |δF| in the fitted data is shown in Figure 5, which plots |δF| against the red-side signal-to-noise per pixel, S/N_red, for PCA fits to a subset of the SDSS spectra as well as the mock spectra described in Section 4.

Figure 5 provides a useful diagnostic for the quality of the PCA fits on the SDSS spectra. We expect the fits to the mock spectra (green crosses) to represent the case in which the PCA eigenspectra describe the spectra nearly perfectly (see Section 4), therefore they typically have smaller values of |δF| than the real spectra. Clearly, the SDSS spectra with large |δF| are most likely bad fits, but note that the presence of metal absorption lines and other artifacts in the real data can bias |δF| to larger values even for good fits (we did not mask any lines when calculating |δF|, as our metal-masking algorithm is rudimentary and sometimes masks legitimate pixels).

In practice, we carry out the PCA fitting procedure using the two different PCA templates described in Section 2.2, then for each SDSS spectrum we select the fit which gives the lower value of |δF|. Fits with |δF| values under the 95th percentile of those from the mock spectra (red line in Figure 5) are then automatically considered good fits, while the rest are visually inspected and flagged for goodness of fit on the red side of the spectrum—approximately 90% of the spectra were adequately fit. The spectra which are not well fitted by our procedure consist mostly of objects which have strong absorption systems at λ_rest > 1216 Å, such as metal absorption from DLAs and weak BAL quasars—while our procedure includes a simple procedure to mask metals, this does not effectively deal with very strong metal absorbers. An example of this is shown in Figure 4(a). There are also quasars with unusual spectral shapes which are not represented in the template spectra described in Section 2.2, such as quasars with weak emission lines (Figure 4(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The spectra now have been had PCA fits carried out on them redward of Lyα, but the predicted continuum, C_PCA, extends blueward of Lyα (λ_rest < 1216 Å). For the objects that are well fitted on the red side of the spectrum, we expect the predicted continua to provide a reasonable prediction for the shape of the Lyα continuum blueward of Lyα, but the overall amplitude is uncertain due to the EUV–NUV power-law break described in the introduction. We now turn to the next fitting step, mean-flux regulation, to constrain the continuum amplitude.

In the least-squares PCA fitting step described in the previous section, we have hitherto used no information blueward of the quasar Lyα emission line due to the absorption from the Lyα forest. However, since the absorption redshift at any point in the Lyα forest is known (z_abs = λ_obs/1216 Å − 1), the average absorption averaged over each sightline can be used to constrain the predicted continuum. In this section, we will describe the use of the mean-flux evolution of the Lyα forest, 〈F〉(z), to regulate the amplitude and slope of the predicted PCA continuum. We refer to this as the "mean-flux regulation" step.

Using the PCA continuum C_PCA(λ_rest) fitted to the observed spectrum, we first extract the Lyα forest transmission F^init(λ_rest) = f(λ)/C_PCA(λ_rest) in the range λ_rest = 1041–1185 Å. The extracted Lyα forest is then divided into bins, and the mean flux, $\bar{F}^\mathrm{init}_\mathrm{bin}(\lambda _\mathrm{bin})$, is evaluated for each bin, where λ_bin = [1070 Å, 1110 Å, 1150 Å] are the central rest wavelengths of each bin.

We now introduce a quadratic fitting function blueward of a pivot point, λ_rest = 1280 Å, to obtain the mean-flux-regulated continuum:

$$\begin{eqnarray} C_\mathrm{MF}(\lambda _\mathrm{rest}) &=& C_\mathrm{PCA}(\lambda _\mathrm{rest})\nonumber \\ & & \times \big(1 + a_\mathrm{MF} \hat{\lambda }_\mathrm{rest} + b_\mathrm{MF} \hat{\lambda }_\mathrm{rest}^2\big), \end{eqnarray} \tag{ 5 }$$

where a_MF and b_MF are free parameters for the fit, while $\hat{\lambda }_\mathrm{rest}\equiv \lambda _\mathrm{rest}/ 1280\ \mathrm{\AA}- 1$. Note that for the lower redshifts (z_QSO ≲ 2.4) in which only a limited portion of the Lyα forest is accessible, we use only the linear parameter, a_MF, in order to avoid overfitting.

We again use least-squares fitting to find the values of a_MF and b_MF that provide the best fit between the extracted mean-flux $\bar{F}^\mathrm{fit}_\mathrm{bin}(\lambda _\mathrm{bin})$ and the external mean-flux constraint 〈F〉(z). In this mean-flux regulation step, the parameters c_z, α_λ, and c_j fitted to λ_rest > 1216 Å are kept fixed. For the mean-flux constraint, we use the power-law-only fit from Faucher-Giguère et al. (2008) without metal correction:

$$\begin{equation} \langle F \rangle (z) = \exp [-0.001845 (1+z)^{3.924}]. \end{equation} \tag{ 6 }$$

Note that we used their measurement with metals included, because we are applying these constraints to the SDSS Lyα forest spectra that have not had metals removed from their sightlines (and would probably be impossible to remove, in the case of the noisier spectra).

In principle, the errors in the continuum fit should now be at the level of a few percent, arising from some combination of the large-scale variance in the Lyα forest and errors in the fitting. In the next section, we will use mock spectra to quantify the level of continuum errors in the MF-PCA technique.

The first step is to create synthetic quasar spectra from PCA eigenspectra by making Gaussian realizations of the PCA weights, c_j, (see Equation (3)) in the manner described in Suzuki (2006). Note that this is an approximation, as the distribution of the weights may not be fully Gaussian, but it does generate realistic-looking quasar spectra. In principle, the PCA eigenspectra used to generate the mock spectra and those used in the fitting procedure should be separate but drawn from the same distribution. While we do have two sets of PCA eigenspectra (Section 2.2), they represent quasars with different luminosities. Hence, we do not expect to be able to use eigenspectra from Pâris et al. (2011) to fit mock spectra generated from the Suzuki et al. (2005) eigenspectra and vice versa. However, since we only use eight principal components in our PCA fitting step (Section 3.1), we can generate our mock spectra using 10 principal components in order to increase the uncertainty in the fitting. Nevertheless, the tests described in this section will primarily apply to the limit in which the PCA eigenspectra are a good representation of the fitted quasars, which we have argued (Section 2.2) is a reasonable assumption.

Therefore, we generate mock quasar spectra in the spectral range λ_rest = 1020–1600 Å, using 10 principal components from the Suzuki et al. (2005) eigenspectra. Next, we need to introduce the Lyα forest absorption to the mock spectra. For this, we use the publicly available Roadrunner Lyα forest simulations⁶ of White et al. (2010). These are N-body simulations with a box size of (750 h⁻¹ Mpc)³ and a grid scale of 187.5 h⁻¹ kpc in which the Lyα forest flux was derived using the fluctuating Gunn–Peterson approximation. The simulations were released in the form of 22,500 Lyα forest sightlines per box, output at redshifts z_box ≈ 2.00, 2.25, 2.50, and 2.75.

For a given mock quasar at redshift z_QSO, we select the simulation box with the closest redshift, z_box. Using Equation (6), we then re-normalize the mean flux of the box to 〈F〉(z = (1 + z_QSO)1100/1216 − 1), i.e., using the absorber redshift corresponding to the λ_rest = 1100 Å in the quasar spectrum. We choose to normalize the mean flux across the entire box rather than in individual spectra in order to preserve the variance across different lines of sight, which is a source of error in the MF-PCA continuum fitting. A random line of sight is selected from the set of skewers, and the transmitted flux in each pixel, F_i, is rescaled to F'_i = F_i × 〈F〉(z_{abs, i})/〈F〉(z = (1 + z_QSO)1100/1216 − 1), where z_{abs, i} is the absorber redshift corresponding to the pixel. This introduces redshift evolution of the mean flux, 〈F〉(z), within the individual sightlines, which had hitherto had a fixed value of 〈F〉.

The simulated Lyα forest absorption is added to the mock quasar spectrum in λ_rest < 1216 Å and smoothed to the approximate SDSS resolution, R = 2000. Gaussian noise is then added to the mock spectrum using the noise array of a randomly chosen SDSS quasar spectrum with the same z_QSO  and S/N. The mock spectra are then run through the MF-PCA fitting process described above to obtain continuum fits.

In Figure 6, we show several examples of the mock spectra and the fitted MF-PCA continua. The first thing to note is that mock spectra look realistic. They look similar to the real spectra shown in Figure 3, apart from the lack of metal absorption redward of the Lyα emission line.

Zoom In Zoom Out Reset image size

As described in Section 3.1, we use the mock spectra as a benchmark for the PCA fit quality on the red side (λ_rest > 1216 Å) of the spectra—we automatically accept all fits with absolute flux error, |δF|, less than the 95th percentile of the |δF| distribution measured from the mocks (Figure 5). For the fits with larger |δF| that require visual inspection, we use the mocks as a visual guide for what constitutes a good fit.

It is clear from Figure 6 that the mean-flux-regulated continuum, C_MF, is a corrected version of the PCA fit, C_PCA. In several cases, the initial PCA continuum,C_PCA, appeared unphysical (e.g., dipping below the peaks of the forest). These were rectified by the mean-flux-regulated fit, C_MF.

We can place this on a more quantitative footing by comparing the fitted continua, C_fit, to the "true" continua, C_true, which is known by construction in the mock spectra. We define the continuum-fitting residual,

$$\begin{equation} \delta C(\lambda _\mathrm{rest}) \equiv \frac{C_\mathrm{fit}(\lambda _\mathrm{rest})}{C_\mathrm{true}(\lambda _\mathrm{rest})} - 1. \end{equation} \tag{ 7 }$$

We can then carry out MF-PCA fits on large numbers of mock spectra to obtain statistics on the continuum-fitting errors as a function of S/N and quasar redshift. In Figure 7, we show the residuals from fitting 1000 mock spectra in two bins of S/N and quasar redshift, z_QSO. The orange lines show the residuals as a function of wavelength, δC(λ_rest), binned into rest frame 1 Å bins from a subset of 100 mocks. The dashed lines show the 1σ dispersion of the residuals estimated from bootstrap resampling at each 1 Å wavelength bin, while the dotted lines show the 10th and 90th percentile at each wavelength. Figure 7(a) represents one of the worst case scenarios—the signal-to-noise (S/N = 2–4) is low at λ_rest > 1216 Å, making it difficult to obtain a good fit for the continuum shape. At the same time, the redshift is sufficiently low that the Lyα forest occupies the blue end λ_obs ≲ 4000 Å of the SDSS spectrographs where the S/N deteriorates rapidly with decreasing wavelength. This causes a large scatter in the flux of Lyα forest even when averaged over large segments, introducing more errors to the mean-flux regulation process.

Zoom In Zoom Out Reset image size

At moderate signal-to-noise (S/N = 6–10; Figure 7(b)), the situation is significantly better. The 1σ dispersion of the fit residuals are well under 10% and 6%–7% in the central portion of the fitted region. Many of the residuals are flat to within a few percent across the Lyα forest region, indicating that the PCA-fitting has successfully accounted for the shape of the Lyα continuum. About one-tenth of the fitted continua are badly fit, with badly predicted continuum shapes and/or residuals greater than 10%.

We also calculate the mean bias, $\overline{\delta C}(\lambda _\mathrm{rest})$, of the residuals in Figure 7. Averaged over ∼10³ spectra, the MF-PCA technique yields a low bias of <1%, although this does not include any systematic errors in the 〈F〉(z) measurement used to constrain the overall continuum levels. Note that this can be significant at 〈z〉 ≳ 3, where the uncertainties in 〈F〉(z) are at the several percent level, and increasing rapidly beyond 〈z〉 ∼ 4.

To quantify the overall fit quality on each mock spectrum, we use the rms of the continuum residuals evaluated over λ_rest = 1041–1185 Å:

$$\begin{equation} (\delta C)_\mathrm{rms} \equiv \left[ \frac{\int \left(\frac{C_\mathrm{fit}(\lambda _\mathrm{rest})}{C_\mathrm{true}(\lambda _\mathrm{rest})} - 1 \right)^2 d\lambda _\mathrm{rest}}{\int d\lambda _\mathrm{rest}} \right]^{1/2}. \end{equation} \tag{ 8 }$$

Figure 8 shows the median rms error from runs of 1000 mocks as a function of redshift, for four different S/N bins. At lower redshifts, the rms error is relatively high for the low-S/N spectra because the mean-flux regulation is affected by the increased noise levels near the blue end of the SDSS spectra at λ_obs ∼ 4000 Å. As the observed Lyα forest region clears the blue end of the spectra, the rms error decreases to a minimum at z_QSO ≈ 3.0. It then rises with redshift at z_QSO > 3 due to the increasing variance in the Lyα forest, which adds to the error in the mean-flux regulation. At fixed redshift, the median rms decreases with S/N as might be expected. Below z_QSO ≈ 3, it drops below 5% rms for moderate signal-to-noise (S/N ∼ 5) and asymptotes to ∼3%–4% rms for the S/N > 10 spectra.

Zoom In Zoom Out Reset image size

To estimate the contribution from various sources of error in the MF-PCA fitting, we also carried out the mean-flux regulation directly on the simulated Lyα forest skewers without introducing a quasar continuum, or adding pixel noise. In other words, we directly fit the function Equation (5) to the mean flux evaluated over three bins in each skewer. The rms error in the continuum from this estimate is shown as the dashed line in Figure 8. For z_QSO = 2.3 quasars, the rms error contribution from the Lyα forest variance is about 1.5%; this increases to 4% at z_QSO = 4. This suggests that even in the limit of high-S/N, errors in the continuum shape from PCA fitting contributes 1%–2% to the overall rms continuum error.

To place the above results in context, we carry out another set of continuum fits on the mock spectra using the continuum model

$$\begin{equation} C_{\rm mean}(\lambda _\mathrm{rest}) = f_{1280}\times \mu (\lambda _\mathrm{rest})\times \left(\frac{\lambda _\mathrm{rest}}{1280\ \mathrm{\AA}}\right)^{-\alpha } \end{equation} \tag{ 9 }$$

in which the mean spectrum, μ(λ_rest), is multiplied with a power law, λ_rest∝λ^−α_rest, and f₁₂₈₀ is the flux normalization at λ_rest = 1280 Å. Both f₁₂₈₀ and α are determined separately for each quasar, with the power law fitted to the regions near λ_rest = 1280 Å and λ_rest = 1450 Å. This model is highly similar to that implemented in Slosar et al. (2011).

In Figure 9, we show the continuum residuals from C_mean, as a function of rest wavelength in the Lyα forest region. Comparing this plot with Figure 7(b), which show the MF-PCA fitting residuals from the same [S/N, z_QSO] bin, it is clear that MF-PCA dramatically reduces the range of fitting errors by a factor of three. Indeed, the C_mean(λ_rest) residuals are significantly larger than even the worst-case scenario for MF-PCA continua, represented in Figure 7(a). The significant bias (black line) of the residuals from C_mean(λ_rest) is puzzling at first glance, as one would expect to recover the mean spectrum (and hence no bias) when averaging over large numbers of spectra. We suspect this bias most likely due to an asymmetry in the distribution of power-law spectral indices in quasars (see Desjacques et al. 2007), which we have not accounted for in our mock spectra. However, this does not affect the scatter and shape of the continua, which is the present quantity of interest. The median rms error from the continua shown in Figure 9 is (δC)_rms = 0.13, which is worse than any of the [S/N, z_QSO] bins evaluated in Figure 8.

Zoom In Zoom Out Reset image size

The rms continuum-fitting error is not the most important quantity for studies of the one-dimensional Lyα forest flux power spectrum:

$$\begin{equation} P_F(k) = 2 \pi \int ^\infty _0 \delta (k) \delta ^*(k) dk, \end{equation} \tag{ 10 }$$

where δ ≡ F/〈F〉 − 1, and k ≡ 2π/l is the Fourier wavenumber. Rather, it is the Fourier power from the continuum errors that is the troublesome systematic. For example, in their measurement of P_F(k) from SDSS data, McDonald et al. (2006) were limited to scales of k ⩾ 0.0014 km⁻¹ s, corresponding to comoving distances of r ≲ 50 h⁻¹ Mpc or λ_rest ≲ 18 Å in the quasar rest-frame wavelength. This was due to the increasing influence of continuum power at large scales. It is therefore pertinent to investigate the amount of residual Fourier power introduced by the various continuum-fitting methods.

In Figure 10, we show the mean power spectrum of the continuum residuals from 1000 mock spectra, δC(λ_rest), for the power-law+mean continuum method, C_mean, and MF-PCA continuum fitting, C_MF, shown for two signal-to-noise bins. All the mock spectra had quasar redshifts in the range z_QSO = 2.9–3.1. All the residual power spectra have the same overall shape, with a bump at k ≈ 0.0005 km⁻¹ s corresponding to the weak emission lines in the intrinsic quasar spectrum at scales of λ_rest ≈ 50 Å. This is unsurprising, since imperfections in the continuum fitting should give similarly shaped residuals. At fixed k, the amplitudes follow the same pattern we had already discussed above: the residual power from MF-PCA fitting is significantly lower than that from the power-law+mean continuum fits. Even with noisy S/N ≈ 3 spectra, C_MF continuum fitting reduces the residual power by ∼30% compared to C_mean fitted to higher-S/N spectra.

Zoom In Zoom Out Reset image size

In their measurement of P_F(k) from SDSS data, McDonald et al. (2006) had used a mean continuum shape for all their spectra, which is similar to C_mean except that they did not fit for the individual power-law indices. Using this method, McDonald et al. (2006) found that residual continuum power started becoming problematic at scales greater than k = 0.0014 km⁻¹ s, indicated by the solid red vertical line in Figure 10. We thus estimate the residual power in C_mean at this scale at which continuum power interferes with P_F(k) measurements. We can then look for the points in the C_MF residual power spectra with the same limiting power (red arrows in Figure 10). Note that this assumes that the Lyα forest power is constant whereas the Lyα forest power decreases with scale, but the change is gradual. For our rough estimates, we can approximate it as constant over small logarithmic intervals.

The corresponding limiting values of k (red dotted lines) are significantly smaller than for C_mean. This suggests that MF-PCA continuum fitting could allow the Lyα forest flux power spectrum to be measured at larger scales than previously possible. For S/N = 6–10 spectra, the lower k-limit is now k = 0.0007 km⁻¹ s. This corresponds to a doubling of the accessible comoving scales: r = 2π/k ≈ 90 h⁻¹ Mpc at z = 2.75 compared with r ≈ 45 h⁻¹ Mpc in the McDonald et al. (2006) study, where these distances are calculated assuming a standard flat ΛCDM cosmology with h = 0.7, Ω_m = 0.28 and with a w = −1 cosmological constant. Even for noisy (S/N = 2–4) spectra, the accessible scale has been increased significantly to r ≈ 65 h⁻¹ Mpc.

The new continuum-limited scales approach the ∼100 h⁻¹ Mpc BAO scale at moderate signal-to-noise (S/N ≳ 6). In Figure 10, the black vertical triple-dot-dashed line indicates the wavenumber of first BAO peak, calculated from the prescription in Eisenstein & Hu (1998). Even though future BAO measurements in the Lyα forest are expected to be carried out in three dimensions (McDonald & Eisenstein 2007), the increase in accessible modes along the lines of sight will improve the robustness and precision of the measurements.

We have carried out MF-PCA continuum fitting on 12,069 quasar spectra from the SDSS DR7 catalog. The continua in the spectra range 1030 Å < λ_rest < 1600 Å have been made publicly available and can be downloaded via anonymous FTP.⁷ The IDL fitting code took ∼0.5 s per spectrum (including file input/output) on a single processor core of a 3.0 GHz Intel Quad Core desktop with 2 GB of RAM, allowing the entire SDSS DR7 Lyα forest sample to be fitted in about two hours.

Since we expect the MF-PCA technique to provide a good continuum fit only when there is a good PCA fit redward of the quasar Lyα line, the fitted spectra have been visually inspected to verify the fit quality in the λ_rest = 1216–1600 Å wavelength region. Approximately 89% of the spectra had reasonable PCA fits, and these have been flagged as such in the publicly available continua, although we recommend that users of the continua should make their own cuts on the fit quality. Approximately 30% of the spectra were better fitted by the low-redshift Suzuki et al. (2005) templates, while the rest were better fitted by the Pâris et al. (2011) templates. This is qualitatively as expected, since there is greater overlap by the Pâris et al. (2011) templates in the luminosity distribution of the SDSS quasars.

From the mock spectra analysis of the fit quality in Section 4, we have estimated the continuum-fitting error at each pixel within the Lyα forest, as a function of quasar redshift and spectral S/N. However, the errors have significant covariances, therefore it would be too unwieldy to provide the full error estimates for each spectrum although we can provide them upon request.

It is worth emphasizing that our continuum estimation methods require a prior assumption on the mean flux of the Lyα forest in our fits (Equation (6)). This means that our continuum estimates cannot be used to measure 〈F〉(z); we expect our continua to be useful primarily for higher-order Lyα forest flux statistics that are less sensitive to uncertainties in 〈F〉(z), such as the flux power spectrum.

While we have studied the performance of MF-PCA continuum fitting on mock spectra in Section 4, it is difficult to empirically constrain the quality of the fits. One possibility is to compare a small subset of the data with high-resolution, high-S/N spectra of the same objects. However, even with high-resolution spectra, it is questionable whether there are sufficient transmission peaks in the forest to adequately constrain the continuum shape; Faucher-Giguère et al. (2008) have shown that accurate fitting of the quasar continuum is difficult at z ≳ 2.5—even though they had corrected for the overall bias using cosmological simulations for their subsequent analysis, it is difficult to recover the true continuum shape for any single spectrum. Furthermore, most high-resolution spectra are obtained from echelle spectrographs with uncertain spectrophotometry, so it would be tricky to directly compare quasar sightlines which have been observed in both SDSS and high-resolution echelle spectrographs.

However, it is possible to get a sense of the efficacy of our MF-PCA continua by stacking large numbers of spectra. This cancels out the Lyα forest power from individual sightlines and allows the underlying continuum shape to be seen, albeit lowered due to the mean Lyα absorption.

Recall that the PCA coefficients, c_j, parameterize the shape of the quasar spectra. Therefore, if our continuum-fitting technique works, quasars with similar c_j measured from λ_rest > 1216 Å should have similar-looking continua within the Lyα forest region. Thanks to the large number of spectra in our SDSS sample, it is possible to stack ∼10² spectra with similar values of c_j to recover the collective shape of their Lyα forest continua.

We can select subsamples of quasars based on their values of σ₁ ≡ c₁/λ₁ and σ₂ ≡ c₂/λ₂, where λ₁ = 7.563 and λ₂ = 3.604 are the standard deviations of c₁ and c₂, respectively, in the low-redshift HST eigenspectra (Suzuki 2006). These two principal components account for approximately 80% of the total variance in the low-redshift quasar templates. We limit ourselves to spectra with S/N > 3 pixel⁻¹, and which have been visually inspected to be decent fits redward of λ_rest = 1216 Å. We also select quasars with z_QSO > 2.6 in order to ensure reasonably complete coverage of the Lyα forest. Within a subsample, each spectrum is first normalized near λ_rest = 1280 Å and rebinned into a common wavelength grid with Δλ_rest = 1 Å bins before being stacked. The same procedure is carried out on the MF-PCA continuum fitted to each spectrum to obtain a mean MF-PCA continuum for the subsample.

In Figure 11, we show four subsamples from our SDSS sample with different [σ₁, σ₂] with respect to the low-redshift Suzuki et al. (2005) eigenspectra. Redward of 1216 Å, we see that the least-squared PCA procedure generally does a good job of fitting the emission lines, although there are inaccuracies in fitting N v λ1240 and Si ii λ1306. Blueward of 1216 Å, the stacked spectrum appears to have a similar shape to the fitted continua, although the overall flux level is depressed due to the mean Lyα absorption.

Zoom In Zoom Out Reset image size

We can make a more direct comparison between the stacked spectra and the fitted Lyα forest continua by correcting each observed Lyα forest pixel by its mean flux (using Equation (6)) before stacking. The mean-flux-corrected Lyα forest is shown as the disembodied red line in Figure 11. It is gratifying to see that the stacked MF-PCA continua generally agrees well with the stacked Lyα forest spectra. Our technique can clearly account for the diversity in quasar continua: Spectra with clear emission-line features (Figure 11(a)) and those with smooth continua (Figure 11(b)) are well differentiated. Because we have corrected each Lyα forest pixel by the same mean flux (Equation (6)) that we have used to carry out the MF-PCA fits, we expect the amplitude of the corrected Lyα forest stacks, in λ_rest = 1041–1185 Å, to match those of the stacked MF-PCA continua, but the tilt and shape of the continuum bears testament to the success of the technique. In addition, since the MF-PCA continua shown in Figure 11 were a subset which used the Suzuki et al. (2005) quasar templates, this suggests that it was appropriate to use low-redshift templates to fit some of the z_QSO ≳ 2 SDSS spectra.

We have introduced MF-PCA continuum fitting, a new technique for predicting the Lyα forest continuum in low-S/N spectra. In tests on mock spectra, we have found that MF-PCA can predict the continuum at the 8% rms level in SDSS spectra with S/N ∼ 2 at z = 2.5, and <5% rms in S/N ≳ 5 spectra. This is a significant improvement over the ∼15% rms continuum errors previously achievable in low-S/N spectra. We are making available MF-PCA continuum fits for 12,069 Lyα forest spectra from the SDSS DR7 quasar catalog. The MF-PCA technique also significantly reduces the Fourier power from continuum-fitting residuals by a factor of a few in comparison with dividing by a mean continuum. This will allow a concomitant increase in the accessible scales for Lyα forest flux power spectrum measurements.

This improved continuum-fitting accuracy will significantly increase the value of low-S/N Lyα forest data. For example, the ongoing Baryon Oscillations Spectroscopic Survey (BOSS; Eisenstein et al. 2011) will obtain Lyα forest spectra from ∼150, 000 quasars at z_QSO ≳ 2 (Ross et al. 2011), with the aim of measuring the BAO feature in the Lyα forest absorption across different quasar sightlines (e.g., Slosar et al. 2011). The typical signal-to-noise (S/N ∼ 2) of BOSS spectra will be even lower than that of SDSS (S/N ∼ 4), therefore we expect the MF-PCA technique to contribute significantly to the utility of the BOSS data.

There are several ways in which the current work could be improved. The Suzuki et al. (2005) PCA templates, which we have used to fit some of the spectra, which were derived from a low-redshift quasar sample may not be a perfect descriptor of the SDSS data (although the test in Section 5 shows that it does a reasonable job). In addition, while the z_QSO ∼ 3 Pâris et al. (2011) templates were indeed obtained from SDSS quasars, they used a high-luminosity subset that is not representative of the full SDSS luminosity distribution. Furthermore, the hand-fitting technique that they had used to obtain continua from these spectra cannot be used for the lower-luminosity (and hence lower-S/N) quasars.

However, with a large data set such as SDSS or BOSS, it is possible to regard the Lyα forest absorption within individual spectra as a noise term that cancels out with sufficiently large numbers of template spectra. This will allow new eigenspectra to be generated from the data itself, although each individual spectrum will need to be corrected by its mean flux before being included in the eigenspectrum solution. In the near future, we will work on this technique to generate new eigenspectra from the BOSS data.

The other issue with the MF-PCA fitting is that it requires an assumed mean flux, 〈F〉(z) for the Lyα forest. This is not ideal, as the evolution of the mean flux is an important observable of the Lyα forest. This could in principle be overcome by solving simultaneously for the mean flux of the Lyα forest and the continuum-fitting parameters for the individual spectra, using maximum-likelihood techniques. This could potentially allow large Lyα forest data sets to be continuum-fitted and studied in a fully self-consistent fashion.