{'fraction_non_alphanumeric': 0.03219329589926834, 'fraction_numerical': 0.022630593840394758, 'mean_word_length': 4.530867683041596, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral Evolution in High Redshift Quasars from the Final BOSS Sample

We report on the diversity in quasar spectra from the Baryon Oscillation Spectroscopic Survey. After filtering the spectra to mitigate selection effects and Malmquist bias associated with a nearly flux-limited sample, we create high signal-to-noise ratio composite spectra from 58,656 quasars (2.1 \le z \le 3.5), binned by luminosity, spectral index, and redshift. With these composite spectra, we confirm the traditional Baldwin effect (BE, i.e., the anticorrelation of C IV equivalent width (EW) and luminosity) that follows the relation W_\lambda \propto L^{\beta_w} with slope \beta_w = -0.35 \pm 0.004, -0.35 \pm 0.005, and -0.41 \pm 0.005 for z = 2.25, 2.46, and 2.84, respectively. In addition to the redshift evolution in the slope of the BE, we find redshift evolution in average quasar spectral features at fixed luminosity. The spectroscopic signature of the redshift evolution is correlated at 98% with the signature of varying luminosity, indicating that they arise from the same physical mechanism. At a fixed luminosity, the average C IV FWHM decreases with increasing redshift and is anti-correlated with C IV EW. The spectroscopic signature associated with C IV FWHM suggests that the trends in luminosity and redshift are likely caused by a superposition of effects that are related to black hole mass and Eddington ratio. The redshift evolution is the consequence of a changing balance between these two quantities as quasars evolve toward a population with lower typical accretion rates at a given black hole mass.


INTRODUCTION
Spectral properties of quasars depend on the physical conditions of the accreting black hole and on the geometry of the emitting regions. Quasars display striking similarities such as broad emission lines and continua that follow an approximate powerlaw relationship in the near UV and optical wavelengths (Merloni 2016, and refer-ences therein). However, the exact properties of the continuum and emission lines vary across the population of quasars in a manner that is poorly understood (Richards et al. 2011;Shen & Ho 2014;Baskin et al. 2015, for instance). Our understanding of quasar diversity has improved in recent years through applications of several analysis techniques to the sample of spectroscopicallyobserved quasars. Among these techniques, the identification of trends among observed parameters across the whole population of quasars, principal component analysis (PCA; Boroson & Green 1992, for instance), and the use of composite spectra are among the most commonly employed. These techniques have allowed characterization of the most significant correlations between the continuum and emission line properties of quasars. Constraining the trends across the sample of quasars allows us to better understand the underlying physics of these systems and informs the modeling of quasar spectra for better redshift determination in studies of large-scale clustering.
The application of PCA decomposition to large samples of quasars allows the empirical characterization of the most significant correlations between spectroscopic properties. The PCA technique builds a minimum set of basis vectors across the observed parameter space to capture the total variance in the sample. This method has revealed a striking relationship in which the relative optical to X-ray flux and the equivalent width (EW) of Fe ii are shown to be directly correlated with the EW of [O iii] and the full width at half maximum (FWHM) of H β . Known as Eigenvector 1, this trend was first identified by Boroson & Green (1992) and later confirmed in additional samples (Shang et al. 2003;Yip et al. 2004;Sulentic et al. 2009). The Eddington ratio is believed to be the physical driver of the Eigenvector 1 correlations (e.g., Boroson 2002;Marziani et al. 2003;Shen & Ho 2014).
The similarity of quasar spectra over wide ranges of luminosity and redshift has enabled the use of composite spectra to study their global properties. Composite spectra offer several advantages over the other methods by virtue of their high signal-to-noise ratios and the tendency to average out the intrinsic peculiarities of any single source. Weaker emission lines unidentifiable in single objects are detectable in composite spectra. This method has been employed over a wide range of quasar samples such as the Large Bright Quasar Survey (LBQS; Francis et al. 1991;Hewett et al. 1995), FIRST Bright Quasar Survey (FBQS; White et al. 2000;Brotherton et al. 2001), Sloan Digital Sky Survey (SDSS; York et al. 2000;Vanden Berk et al. 2001), Hubble Space Telescope (HST; Zheng et al. 1997;Telfer et al. 2002), HST-Far Ultraviolet Spectroscopic Explorer (FUSE; Scott et al. 2004), HST-Cosmic Origin Spectrograph (COS; Shull et al. 2012;Stevans et al. 2014;Tilton et al. 2016), Baryon Acoustic Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013;Harris et al. 2016), and VLT-XSHOOTER (Selsing et al. 2016).
Despite the general similarities, quasars differ in their luminosity, spectral index, and emission line features. Unrestricted composite spectra mask the intrinsic dispersion in these key parameters which would otherwise provide clues to underlying physics. Composite spectra constructed from a subsample of sources binned by common observed properties allow the study of correlated properties in a manner which does not necessarily impose the orthogonality constraints of the PCA approach. For example, Richards et al. (2011) constructed composite spectra from ≈30,000 quasars from the 7th Data Release of SDSS as a function of C iv equivalent width and blueshift to reveal the empirical relationships between emission features of line species and the continuum. They reported that the C iv emission parameters capture an important relationship between the disk and wind components of the broad emission lines that appears to be sensitive to the shape of the ionizing continuum. Ivashchenko et al. (2014) constructed composite spectra binned by spectral index and showed that their composite spectra resulted in higher average systemic redshift estimates than the SDSS templates.
One of the most prominent observed relations in quasar spectra is the anti-correlation between the rest-frame equivalent width of UV emission (particularly C iv) and continuum luminosity, first identified by Baldwin (1977). The dependence of equivalent widths on luminosity is commonly referenced as the Baldwin effect (BE). The BE is now well established for nearly all broad emission lines. Kinney et al. (1990) and Zamorani et al. (1992) expanded the original work to larger ranges of quasar luminosity and redshift, thus revealing similar effects in Mg ii (see Sulentic et al. 2000, for an outline). Green et al. (2001) discovered a significant BE in UV iron emission, although also found that UV iron emission may have stronger correlation with redshift, implying an evolution-ary effect rather than a direct dependence on luminosity. Croom et al. (2002) and Dietrich et al. (2002) revealed the effect in narrow emission lines using composite quasar methods. The BE in the X-ray Fe Kα feature has also been identified in a number of studies (e.g., Page et al. 2005;Wu et al. 2009;Ricci et al. 2013). Baldwin et al. (1989) and Dietrich et al. (2002) showed that the slope of the BE steepens with increasing ionization potential.
In quasars with multiple observations, the intrinsic BE associated with change in luminosity has also been identified for different broad emission lines. The same trend of a decreasing EW with increasing luminosity is seen, but with a steeper slope than the ensemble BE that is observed over a population of quasars (see Pogge & Peterson 1992;Gilbert & Peterson 2003;Cackett & Horne 2006). Among the several physical explanations proposed to account for the BE and its ionization dependence, the most popular is a prediction that the continuum softens at short wavelengths with increasing luminosity. Models predicting this relationship result in fewer ionizing photons in luminous quasars (Netzer et al. 1992;Korista et al. 1998).
Despite the myriad of studies and observations of the BE, there are no rigorous models that provide a complete explanation of the effect's origin (for a review, see Shields 2007, and references therein). The BE remains an empirical correlation, so it is not clear if the variations in line strengths are driven by luminosity or by a more fundamental parameter such as black hole mass, Eddington ratio, or orientation. Because these parameters are all correlated with luminosity, it is difficult to disentangle the origins of the effect (Bachev et al. 2004;Shields 2007;Bian et al. 2012). By creating composite spectra over different ranges in FWHM of C iv emission, Warner et al. (2004) provide evidence that black hole mass may be largely responsible for the EW variations associated with luminosity. In contrast, after binning on Eddington ratio, the composite spectra do not show BE-like behavior. However, Kollmeier et al. (2006) suggest that the interpretation of these results may be clouded by residual selection effects.
In this study, we explore spectral diversity using composite spectra over an unprecedented sample of 112,000 quasars. We mitigate biases in the distribution of colors and bolometric luminosities that are intrinsic to a fluxlimited survey in each of the three redshift intervals. We generate high signal-to-noise ratio composite spectra at each redshift in bins of constant luminosity and spectral index.
We use these composite spectra to study the changes in quasar spectral features as a function of redshift, luminosity, and spectral index. By characterizing diversity through a linear interpolation between composite spectra, this approach is similar to PCA. Unlike the PCA approach, we impose no orthogonality restriction which enables one to find common modes in diversity associated with each of these three parameters.
The paper is structured as follows: Section 2 describes the data, the binning and sampling methods, and the construction of composite spectra. Section 3 presents the variation in spectral signature associated with luminosity, spectral index, and redshift. In Section 4, we probe different observed quasar parameters that may be responsible for the variation in spectral signatures seen in luminosity and redshift. Section 5 summarizes our findings and discusses possible future work. Throughout, we assume a ΛCDM cosmology with H 0 = 70 km s −1 , Ω M = 0.3, and Ω Λ = 0.7.

DATA
We utilize quasars observed in the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013), a component of SDSS-III (Eisenstein et al. 2011). Using the Apache Point 2.5-meter Telescope (Gunn et al. 2006), SDSS-III and BOSS concluded in 2014. The final public data release is the 12th in the series from SDSS, comprising 2,401,952 galaxy and 477,161 quasar spectra (DR12; Alam et al. 2015).
BOSS was designed to characterize the Baryonic Acoustic Oscillation (BAO) feature first detected roughly a decade ago in the clustering of galaxies (Eisenstein et al. 2005;Cole et al. 2005). BOSS used luminous red galaxies and the Lyman-α forest observed in the absorption of quasar spectra to measure the BAO signal in the correlation of matter. Spectroscopic observations are performed using plates with a radius of 1.5 • on the sky (Smee et al. 2013). On each plate, there are 1000 holes drilled for 2 arcsecond diameter fibers, where approximately 600-700 fibers are used for galaxy targets and 160-200 fibers are used for quasar targets. A total of 20 fibers are reserved for standard stars with ugriz colors consistent with main sequence F stars. The most recent BOSS BAO measurements from galaxies are reported in Alam et al. (2016) and from the Lyman-α forest in Delubac et al. (2015).
The majority of quasar targets were selected from SDSS ugriz imaging data (Fukugita et al. 1996). Additional selections included imaging data from the Galaxy Evolution Explorer (GALEX; Martin et al. 2005) and UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007). Together, the selection methods (Richards et al. 2009;Yèche et al. 2010;Kirkpatrick et al. 2011;Bovy et al. 2011Bovy et al. , 2012 resulted in a total density of confirmed Lyman-α quasars exceeding 15 deg −2 (Ross et al. 2012). The quasars used in Delubac et al. (2015) span the redshift range 2.1 ≤ z ≤ 3.5. In the analysis presented in this work, we define a control sample based on spectroscopically observed properties. We are therefore less susceptible to target selection biases in the BOSS sample and we include all targets regardless of their selection algorithm.

Data Reduction
Once the observations are completed, there are two major steps that comprise data reduction.
Onedimensional spectra are extracted and calibrated from the raw, two-dimensional CCD images (Stoughton et al. 2002). The extracted spectra are flux-calibrated based on model spectra fit to the aforementioned F stars found on each plate. The second step of the data reduction pipeline performs a classification and a redshift estimate for each spectrum (Bolton et al. 2012).
Significant changes have been made to the spectral extraction and coaddition of individual exposures since the DR12 release. These changes are motivated by the challenges of classifying lower signal-to-noise ratio spectra found in the program succeeding BOSS. Like BOSS, the extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016) is the cosmological component of an iteration of the SDSS; this time SDSS-IV (Blanton et al. in prep.). The first major change corrects a known bias in the coaddition of individual exposures. This correction has significant impact on the classification of galaxy spectra and is described in Hutchinson et al. (2016). The second change pertains to flux calibration and directly impacts the interpretation of quasar spectra. Because it is relevant to this work, we provide a description below.
As described in Dawson et al. (2013), the atmospheric differential refraction (ADR) bears a different signature in quasar spectra than the standard stars due to a choice in the placement of fibers where quasars are drilled to maximize the flux at 4000Å instead of 5400Å for standard stars. Harris et al. (2016) provide an empirical, global correction to this effect by using standard stars observed in the same fashion as quasar targets. Likewise, Margala et al. (2015) derive a first principles estimate of ADR for each quasar target and apply corrections to the final coadded spectra. In the newest reductions, we have improved upon both of these analyses by applying the flux corrections derived from our ADR model to quasar spectra in individual exposures before any flux calibration is applied. This approach has the distinct advantage of correcting for ADR flux loss in 15-minute increments corresponding to the individual exposure time, rather than making a correction to the final coadded spectrum based on the averaged flux loss.
The new data reduction pipeline is complete and has been used to reprocess all quasar data taken during BOSS. Tagged as v 5 10, these reductions of BOSS data and all eBOSS data taken in 2014-2016 will be released as DR14 in Summer 2017. The results show significant improvement in the flux calibration of quasar spectra. We quantify spectrophotometric errors by examining the residuals in the spectra of stars that appear as contaminants in the selection algorithms for BOSS quasars because these were observed in the same manner as quasars. We compute synthetic magnitudes in filters corresponding to gri from these stellar spectra using the DR13 reductions and the new DR14 reductions. By comparing those synthetic magnitudes to the gri magnitudes from the SDSS imaging data, we can characterize the relative systematic errors in flux calibration. As shown in Fig For our analysis, we use the new calibrated spectra found in DR14 and the classifications from the quasar catalog released along with DR12 (DR12Q; Pâris et al. 2016). Several techniques were implemented in DR12Q to improve the redshift precision, classification, and identification of peculiar quasars relative to the BOSS pipeline classifications. We use the PCA technique as the estimate of the redshift.

Sample Selection
There are 175,294 quasars in the redshift range 2.1 ≤ z ≤ 3.5 included in the DR12Q catalog. We restrict the population to that redshift range as it has the most com- Stars g < 20 λ eff = 4000Å Figure 1. Residual bias between SDSS imaging and BOSS spectrophotometry computed for 1330 stars brighter than g = 20 that were targeted as quasars. The histogram for each panel shows the difference between the DR13 reductions (red) and the new improved DR14 reductions (blue) with ADR corrections. Residuals are presented for gri and r − i.
plete selection and will allow better understanding of the quasar population used in BOSS Lyman-α forest studies. We remove broad absorption line quasars (BALs) that are identified in DR12Q with the flag BAL FLAG VI due to their influence on emission lines. Damped Lyman-α (DLA) quasars identified by the same technique as those in the catalog from Noterdaeme et al. (2012) are removed to suppress contamination of quasar continua from metal absorption lines. After application of these constraints, 143,341 quasar spectra remain. Additionally, all quasars observed at an airmass above 1.2 are removed, as the spectrophotometric correction (see Harris et al. 2016) may be less certain at high airmass. The remaining sample is large enough that removing the ∼30,000 high airmass quasars will not significantly affect the final results. In total, 112,000 quasars were included in this study. The sample size at each phase of filtering is reported in Table 1.
In any flux-limited sample, the minimum luminosity increases with increasing redshift, so the range of allowed black hole masses and Eddington ratios able to account for a given luminosity at a higher redshift is then narrowed relative to a lower redshift (e.g., Fig. 2 of Sulentic et al. 2014, AAp, 570, A96). To inspect spectral variations, it is important to choose parameters that can be derived from both the low and high redshift ranges of the sample in identical fashion. While it is difficult to directly measure black hole mass and Eddington ratio, luminosity and spectral index impose selection effects and can be easily measured. In this work we explore spectral diversity against three parameters: bolometric luminosity (L bol ), spectral index (α λ ), and redshift (z). These are all direct spectroscopic observables; L bol is the parameter which is often used to explain variations in emission line properties, while α λ plays a major role in determining quasar colors. Because color measurements are key to the selection of quasar targets, variations in the spectral index with variations in emission line strengths can cause different degrees of completeness in different redshift ranges (e.g., Richards et al. 2002).
The quantity L bol is calculated from the monochromatic luminosity at 1350Å, where L bol = c×L 1350 , L 1350 is the continuum luminosity at 1350Å, and c = 4, consistent with the values in Xu et al. (2008) and Richards et al. (2009). The spectral index is defined by fitting the powerlaw f (λ) = bλ α λ to the estimated quasar continuum. We choose a range of wavelengths to model the quasar continuum in which visual inspection of the high signal-to-noise ratio composite spectrum from Harris et al. (2016) reveals no clear features. We use two regions: 1440Å to 1480Å and 2160Å to 2230Å. We adopt FWHM measurements from DR12Q, and record L bol and α λ and their associated measurement errors for every one of the 112,000 quasars. We divide the sample of 112,000 quasars according to distinct observable properties to generate a series of composite spectra. From these composite spectra, we explore the systematic changes in spectral features associated with changes in luminosity, spectral index, and redshift.

Composite Spectra with Controlled Selection
Because the quasar sample used in our analysis spans a large redshift range, and is selected primarily by ugriz colors, the distribution of L bol and α λ is subject to artificially-induced redshift dependence. Malmquist bias introduces a selection effect in which lower luminosity quasars are preferentially selected at low redshift and are typically rejected at high redshift. Likewise, the color se-lection performed in the observer frame samples quasars with different spectral indices at different redshifts due to the presence of emission lines. When generating composite spectra, we must control diversity within each bin such that there is no redshift evolution in L bol or α λ .
We first bin over three redshift intervals from 2.1 ≤ z ≤ 2.35, 2.35 ≤ z ≤ 2.6, and 2.6 ≤ z ≤ 3.5, chosen to ensure a similar number of objects in each bin. We next require the distribution of L bol and α λ in each of these three redshift bins to have the same high-order statistics by first constructing a two-dimensional histogram spanning these two parameters. We clip the distributions so that only objects within two standard deviations from the mean are included. At each bin (L bol , α λ ) in the histogram, we identify the redshift interval with the minimum number of objects and randomly down-sample from the other two redshift intervals until all three have the same number of objects.
The reshaping of the quasar population results in the removal of ∼ 45% of the spectra (see Table 1). We divide each of the three redshift bins into thirds of log(L bol ), denoted L low , L mid , and L hi for the lowest, middle, and highest luminosity bins, respectively. We divide each of these nine bins into thirds of α λ , denoted α low , α mid , and α hi for the bluest, intermediate, and reddest spectral indices. These cuts are indicated in Figure 2. We do not require the sample to be symmetric in both dimensions, so the ranges for α λ in the luminosity bins vary. The lowest luminosity bin's α λ ranges are shifted by −0.04 from the middle luminosity bin's α λ ranges, and the highest luminosity bin's α λ ranges are shifted by 0.04. This asymmetry between bins is exacerbated when considering the tails of the parameter's distributions and is particularly prevalent in spectral index. Thus when we explore the variations in composite spectra along one dimension, we utilize the middle bins (mid) of the fixed parameter because they are more consistent throughout. By design, the variability across redshift bins is much smaller and neglected. The median values of redshift, L bol , and α λ that comprise each of these bins is summarized in Table 2. In addition, we provide the number of spectra, median FWHM of the C iv emission line, and median signal-to-noise ratio per pixel in the region used to determine bolometric luminosity. The composite spectra representing either the middle bin in luminosity or the middle bin in spectral index are each derived from more than 1,000 individual spectra. The signal-to-noise ratio only varies by a factor of 1.5 between the lowest redshift bin and the highest redshift bin for any given luminosity bin.
We generate composite spectra from the quasars in each of the 27 bins. Spectra were first de-redshifted and resampled to a common wavelength grid. These spectra were normalized by the median flux in the same wavelength range used to find L bol (1340 ≤ λ ≤ 1360Å). Pixels with corrupted flux (such as those from poorly subtracted sky features) are masked as described in Harris et al. (2016). At every wavelength pixel, we simply take the median from all quasars that contribute unmasked flux after removing outliers more than three standard deviations from the median. Statistical errors for the measured flux are estimated through bootstrap methods. The resulting high signal-to-noise ratio composite spectrum in each bin represents the specific subset of the quasar population from which we study spectral variation and evolution.
The composite spectra, ordered from low luminosity at the top to high luminosity at the bottom, are displayed in Figure 3. The traditional BE is evident in the presence of strong emission features at the top of the figure and the weaker emission features at the bottom. Figure 4 represents an expanded view of the composite spectra in three wavelength regions featuring interesting emission features. There is clear redshift evolution in the line strengths even when fixing luminosity and spectral index. Table 2 The median value of all single measurements in each composite spectrum bin for redshift, log(L bol ), α λ , signal-to-noise ratio (S/N), and C iv FWHM, and the number of objects in each bin. C iv FWHM are obtained from the DR12Q catalog.

VARIATION IN SPECTRAL SIGNATURE
The 27 composite spectra are used to compare changes in spectral features associated with distinctly different quasar populations. Inspecting the traditional parameter used in BE analyses, we investigate equivalent widths (W λ ) for lines of interest in each of the composite spectra. We calculate W λ for C iv because it was most commonly used in previous analyses. In addition, we calculate W λ for the Si iv & O iv blend located at approximately 1400Å; this feature is relatively isolated from neighboring lines, and, at 33.5 and 113.9 eV, the ionization potential is significantly different than the 47.9 eV potential of C iv (see Table 3). We could have chosen other lines, such as C iii; however, the other strong lines tend to have complex blends where it is difficult to fit the continuum.
When computing the equivalent width, we first fit a powerlaw to the estimated continuum in the vicinity of the relevant lines. For C iv, we use the interval [1445,1465]Å and [1700,1705]Å to determine the continuum; for Si iv & O iv, we use the interval [1360,1370] A and [1425,1440]Å. We determine the line flux via a discrete integral computed after continuum subtraction. For C iv we integrate the relative flux density over the region 1500 to 1600Å, and for Si iv & O iv we use 1370 to 1425Å. Figure 5 presents log(W λ ) as a function of log(L bol ), α λ , and z. Because of the remaining diversity in the tails of the spectral index distribution described at the end of Section 2.3, only the results for L bol and redshift using the intermediate bin in spectral index (α mid ) are shown. When displaying the results of W λ versus α λ , we use the intermediate bin in luminosity (L mid ). Thus the equivalent widths of only nine of the 27 composite spectra are presented in each panel.
Following Baldwin (1977), we fit the functional form W λ ∝ L βw to the observed log(W λ )/log(L bol ) relationship. The results of the fit for each redshift bin are shown in Figure 5 and found in Table 4. The slope of the relationship is roughly a factor of three lower for the Si iv & O iv blend than it is for C iv. This result is expected if Si iv is the dominant line in the complex, as the higher ionization potential species display a greater anticorrelation with luminosity (see review by Shields 2007). The C iv slopes here are steeper than the results from the application of the BCES(Y |X) method used in Bian  (2012), where they find β w = −0.238 ± 0.040. The slope of equivalent width with luminosity grows steeper by ∼ 20% as the redshift increases from z = 2.25 to z = 2.84. The same trend is seen in C iv and Si iv & O iv. Bian et al. (2012) report hints of an opposite trend with redshift; however, they acknowledge that the dynamic range in luminosity is quite different across the redshift range of their sample, thus possibly affecting the linear regression. Because we have carefully controlled the sample to have the same range of luminosity in each redshift bin, our analysis should be less susceptible to such a bias. We next explore the W λ correlation with α λ . As seen in the second panel of Figure 5, the slope of the relation log(W λ ) ∝ α λ m α is much shallower for spectral index than it is for luminosity, even though the range of both parameters spans nearly the entire population of quasars in our sample. From Table 4, it is also evident that the equivalent width of Si iv & O iv has a stronger dependence on spectral index than does the equivalent width of C iv.
Surprisingly, in the third panel of Figure 5, there is redshift evolution in equivalent width at fixed luminosity and spectral index. The equivalent widths appear correlated with redshift following the relation log(W λ ) ∝ zm z , whereby equivalent width is increasing with increasing redshift. The slope of the fit for each luminosity bin is listed in Table 4. As with the equivalent widthluminosity relationship, the trend is much steeper for C iv than it is for Si iv & O iv. Also evident in the table, the trend steepens with decreasing luminosity. Preliminary findings of a similar nature were reported in Vanden , but no conclusive findings have been subsequently reported.
We expand upon the correlations with equivalent width by examining spectral diversity and evolution as a function of wavelength in the high signal-to-noise ratio composite spectra. We explore diversity through differential spectra used to provide a pixel-by-pixel linear interpolation between composite spectra. The differential spectra are computed ∆f ∆x = f (λ)−f (λ) mid x−x mid for each observed parameter x = log(L bol ), α λ , and z. As above, we determine ∆f ∆log(L bol ) and ∆f ∆z only for the composite spectra constructed with the intermediate spectral index (α mid ≈ −1.52). For ∆f ∆α , we use the intermediate luminosity bin (L mid ≈ 46.26). We discuss the results for each parameter in the following subsections. The filled red circles correspond to the equivalent widths measured for C iv and the dotted red line corresponds to the best linear fit. The blue color scheme corresponds to Si iv & O iv. The Si iv & O iv equivalent widths are much weaker and are scaled by a factor of six for illustrative purposes. Each series is offset by 0.5 in log(W λ ) for clarity. Error bars are included for all measurements, but are smaller than the symbols that represent the data.

Relationship Between Luminosity and Line
Features The differential spectra representing ∆f ∆log(L bol ) are shown in Figure 6. Here, the traditional BE is presented as an inverse correlation between f (λ) and luminosity, i.e., ∆f ∆log(L bol ) is negative in the vicinity of all strong emission lines. The biggest effect is seen in Ly-α and C iv. Smaller effects exist around the C iii λ1908 blend and Mg ii λ2799. Although Ly-α displays the strongest effect, most likely due to its intrinsic strength, the smaller effects in C iii and Mg ii are consistent with the trend that lower ionization potential species display a smaller anti-correlation. N v is a possible exception to this general trend (as outlined in Shields 2007), but the linear variation shown in ∆f ∆log(L bol ) appears fairly strong relative to the intrinsic line strength.
When examining ∆f ∆log(L bol ) over the three redshift bins, one clearly sees more pronounced emission line features at higher redshift. In the differential spectra representing the linear change from log(L bol ) = 45.95 to 46.28, the peak in C iv increases by ∼ 60% from z = 2.25 to z = 2.84. Similar changes in structure are present in the other prominent lines.
There is an asymmetry in the wings of C iv in Figure 6; the red side of the line shows more integrated flux than the blue side. This asymmetry leads to a systematic suppression of flux on the red side of the emission feature with increasing luminosity. The resulting blueshift of the central wavelength of C iv as a function of luminosity is clear in the right panel of Figure 4. Richards et al. (2011) report a related trend that higher C iv blueshifts are correlated with higher luminosity sources.
Finally, there is an asymmetry in the differential spectra determined between log(L bol ) = 45.95 to 46.28 and log(L bol ) = 46.28 to 46.66. There is stronger variation in lines when interpolating over the lower luminosity range than over the higher luminosity range in each redshift bin.
The signature of ∆f ∆log(L bol ) bears striking resemblance to the first eigenspectra derived from principal component analysis (PCA) in the results of Shang et al. (2004), Suzuki (2006), and Pâris et al. (2011). Shang et al. (2004) report that their first eigenspectrum is correlated with luminosity in a similar fashion to the BE. In all of these works, the first eigenspectrum possesses the same sign of variation across all emission lines, as displayed in Figure 6.  ∆log(L bol ) differential spectra derived from composite spectra at fixed spectral index and redshift. Emission lines of interest are labeled and marked with dotted lines. The α mid bin is used for the spectral index with a mean of α λ = −1.52. To reduce confusion, the differential spectra of each redshift interval are vertically offset by 0.5 units.

Signature of Spectral Index
At the beginning of Section 3, we found a trend of equivalent width with spectral index, where the red quasars show stronger emission relative to continuum than the blue quasars. The analysis of the differential spectra presented in Figure 7 reveals the same trend that all the emission lines vary in flux in a correlated fashion, increasing in flux with redder spectra (∆α > 0). Deeper investigation shows that the lines in the differential spectra are shifted blueward of the central wavelength. For example, we measure an average ∼ 1700 km s −1 outflow within C iv per unit change in α λ . However, the center of the lines displayed in Figure 4 do not suggest any noticeable change in the peak wavelength associated with changing spectral index. The trend seen in Figure 7 therefore indicates that the strength of a second outflowing component associated with these emission lines may vary with spectral index. Also evident is a stronger blueshifting with higher ionization potential species. N v is blueshifted on average by ∼ 2400 km s −1 per unit change in spectral index, and Ly-α only by ∼ 1400 km s −1 , although estimating the latter off-set may be more difficult due to the Ly α forest. O iv is heavily blended with Si iv, making it even harder to decouple the two.
A caveat in the analysis of ∆f ∆α is the appearance of some discrepancies between redshift bins. For example, ∆f ∆α computed over the reddest composite spectra at z = 2.25 shows an absorption profile in C iv that does not appear in any other ∆f ∆α interval. Likewise, the signature of ∆f ∆α varies considerably around N v. These levels of discrepancy do not appear in the modes of variation on luminosity. It is possible that the spectral index analysis is susceptible to systematic errors in the redshift estimation (see Denney et al. 2016), calibration, or modeling the continuum. We briefly explored the signature of ∆f ∆α with other techniques of removing or distorting the continuum and found similar discrepancies. For now, we highlight the interesting features in ∆f ∆α , but refrain from deeper analysis until a later date. Figure 7. Nine ∆f ∆α differential spectra derived from composite spectra at fixed luminosity and redshift. Emission lines of interest are labeled and marked with dotted lines. The L mid bin is used for the luminosity with a mean of log(L bol ) = 46.26. To reduce confusion, the differential spectra of each redshift interval are vertically offset by 2 units.

Redshift Evolution
As shown in Section 3 and Table 4, the equivalent width of both C iv and the Si iv & O iv complex change with redshift in a manner that is quite similar to the change with luminosity. An examination of the differential spectra representing ∆f ∆z , shown in Figure 8, displays this trend in more detail. The behavior of equivalent width with redshift appears as a positive correlation between redshift and the surrounding flux density of each emission line. As with ∆f ∆log(L bol ) , the amplitude of the correlation is the strongest around Ly-α and C iv.
Just as ∆f ∆log(L bol ) displayed a higher amplitude of the correlation at high redshift, the ∆f ∆z differential spectra show a stronger correlation at low luminosity. Comparing the highest luminosity bin to the lowest luminosity bin of the higher redshift interpolation, reveals a ∼ 1000% increase in the peak line strength of C iv. The peak flux density in the correlation of the Lyman-α line increases by ∼ 1100% over the same change of luminosity. As shown in Table 2, the signature of ∆f ∆z appears to be related to the change in C iv FWHM. For the bin corresponding to the lowest luminosity and intermediate spectral index, the C iv FWHM is 4162, 4018, and 3482 km s −1 at z = 2.25, 2.44, and 2.76, respectively. At the bin corresponding to the highest luminosity and intermediate spectral index, the FWHM is 4866, 4871, and 4547 km s −1 at the same redshifts. The average FWHM increases by 19.5% from z = 2.76 to z = 2.25 for the low luminosity bin, but only increases by 7% in the high luminosity bin.
In another similar trend to ∆f ∆log(L bol ) , the differential spectra of ∆f ∆z demonstrate asymmetry in the linear interpolations. The amplitude of the correlation is stronger when interpolating over the higher redshift composite spectra than the lower redshift range. At both low and high luminosity, the average FWHM barely changes between z = 2.25 and z = 2.44, while there is significant change in the average FWHM between z = 2.44 and z = 2.76. The same trends with FWHM can be seen in the ∆f ∆log(L bol ) differential spectra; larger changes in C iv FWHM lead to higher amplitudes in the differential spectra.
The striking resemblance of the redshift signature to the luminosity signature motivates us to quantify the Figure 8. Nine ∆f ∆z differential spectra derived from composite spectra at fixed spectral index and luminosity. Emission lines of interest are labeled and marked with dotted lines. The α mid bin is used for the spectral index with a mean of α λ = −1.52. To reduce confusion, the differential spectra of each redshift interval are offset by 3 units. similarity between all differential spectra. To compare the three modes of diversity, we take the mean of all six differential spectra shown in Figure 6 (although inverted to reveal similarities), all six in Figure 7, and all six in Figure 8. In doing so, we assume a common mode between all realizations of ∆f ∆x and improve the signal-tonoise ratio for each of the three differential spectra. The high quality differential spectra are shown in Figure 9.
There is a strong correlation between ∆f ∆z and ∆f ∆log(L bol ) . The expanded view on Lyman-α shows the same ratio of the change in Lyman-α peak flux relative to N v. The same trend is present when examining the change in the peaks of Si iv and O iv emission. We quantify the correlation between the luminosity and redshift signature by computing the correlation coefficient between the differential spectra. The correlation coefficient is computed over the range 1216Å < λ < 2000Å to exclude the effects of evolution of the mean optical depth in the Ly-α forest and to avoid introducing noise in the featureless part of the spectra redward of C iii. This analysis produces cor( ∆f ∆z | ∆f ∆log(L bol ) ) = −0.98, implying that the observed redshift evolution has a nearly identical signature to that of the traditional BE.
The ∆f ∆α differential spectrum shows the same overall trend that all emission lines vary in the same direction. However, the expanded view of Figure 9 highlights the blueshifting of emission lines discussed in Section 3.2. The blueshift in ∆f ∆α does not appear in ∆f ∆log(L bol ) , signifying a fundamental difference between these two modes of diversity. Removing the broadband color by distorting the composite spectra to have the same spectral index yields cor( ∆f ∆α | ∆f ∆log(L bol ) ) = −0.57. Figure 9. Mean of all the differential spectra, where ∆f lux ∆log(L bol ) is divided by −2.5 for illustration. The top panel displays differential spectra vertically offset by 1 arbitrary unit in flux density with emission lines of interest labeled and marked with dotted lines. The bottom set of panels offers extended views of pertinent features with the dotted lines reflecting positions of the central wavelengths of emission lines. From left to right they are Ly-α λ1216, N v λ1240, Si iv λ1397, O iv λ1402, and C iv λ1549. Each composite in the bottom panel is offset by 1.0, 0.1, and 0.6 units from left to right, respectively, to reveal comparisons. These figures are designed only to reveal the similarities in the modes, note that the y-axis here is invalid.

DISCUSSION
Having inspected the signatures of luminosity, spectral index, and redshift imprinted on quasar spectra, we further investigate quasar parameters that may be responsible for the highly correlated ∆f ∆log(L bol ) and ∆f ∆z signatures. The most obvious remaining independent parameters to explore are those associated with emission lines. Of the strong lines, C iv is both visible at all redshifts covered by this analysis and not contaminated by absorption in the Ly-α forest. As shown in Table 2 and discussed in Section 3.3, the FWHM of the C iv emission line evolves toward smaller values with increasing redshift in each bin of fixed luminosity. The FWHM in C iv emission decreases by 658 km s −1 (16%) from z = 2.25 to z = 2.76 in the lowest luminosity bin. Examination of Figure 8 reveals that the highest amplitude of variation with redshift occurs for this low luminosity bin. Conversely, there is almost no change in the average FWHM over the redshift interval 2.25 < z < 2.47 for the highest luminosity bin; the ∆f ∆z computed between these two composite spectra possesses the lowest amplitude of variation. Overall, at a fixed luminosity, the spectroscopic signature of red-shift evolution is tightly related to the redshift evolution of C iv FWHM.
Examining the change in C iv FWHM associated with luminosity reveals a similar trend to that in redshift. At z = 2.25, the FWHM spans 728 km s −1 (17.6% increase) over the interval 45.95 < log(L bol ) < 46.66. At z = 2.77, the FWHM increases 1107 km s −1 (31.8%) over the same interval in luminosity. As shown in Table 4, the BE slope is significantly steeper at high redshift than it is at low redshift. Figure 6 shows a clear correlation in the amplitude of ∆f ∆log(L bol ) with the span in FWHM. Just as with ∆f ∆z , the emission lines in ∆f ∆log(L bol ) show a larger decrease in strength when there is a larger increase in C iv FWHM.

Spectral Diversity with FWHM
These trends in the amplitude of spectroscopic variations motivate an investigation of the change in spectral profile against FWHM at fixed luminosity. We construct composite spectra by splitting the sample (see Table 1) used in Section 3 by another factor of two. For each z, log(L bol ), and α mid bin, we split the sample between quasars with FWHM < 5000 km s −1 and with FWHM > 5000 km s −1 . To study the dependence on the line chosen for the split, we use both C iv and Mg ii (where available) FWHM. Co-addition of all spectra within each bin results in 18 composite spectra for C iv and six for Mg ii. Examples of the composite spectra in the lowest redshift bin are shown in Figure 10.
There is a stronger effect on equivalent width when splitting on C iv than when splitting on Mg ii FWHM. The reasons for this behavior are not obvious and will be the subject of a future work. Focusing on C iv FWHM, there is a correlation of decreased equivalent width with increased FWHM even when holding the luminosity fixed; this is seen strongly in both C iv and Ly-α.
For each of the nine bins of fixed luminosity and redshift, we take the differential spectrum representing ∆f ∆F W HM using the C iv FWHM measurements; the results are shown in Figure 11. The emission line features are negatively correlated with C iv FWHM, indicating that within bins of fixed luminosity, an increase in C iv FWHM results in a reduction of equivalent width. The differential spectrum for FWHM computed at the highest luminosity bin and the mean differential spectrum for luminosity are compared to the ξ 1 eigenvector from Pâris et al. (2011) in the bottom panel of Figure 11. The Pâris et al. (2011) ξ 1 eigenvector more closely matches our ∆f ∆F W HM line ratio between Ly-α and C iv than it matches that line ratio in ∆f ∆log(L bol ) . Computing the correlation coefficients between the ξ 1 eigenvector, the FWHM differential spectrum, and the luminosity differential spectrum, the correlations are above 92% for each of the three pairs. While the correlation coefficients are high for all comparisons, the subtle difference in line ratios (e.g., Ly-α/C iv) implies that ∆f ∆log(L bol ) is a different signature than ∆f ∆F W HM , whereas in the ξ 1 eigenvector and ∆f ∆F W HM are more closely related. Of equal importance in Figure 11, there is much smaller redshift evolution in the ∆f ∆F W HM differential spectra than in the other differential spectra. For example, the peaks in C iv in ∆f ∆log(L bol ) change by roughly 60% between redshift z = 2.25 and z = 2.84. In the highest luminosity bin, the peak in C iv in ∆f ∆F W HM only changes by 15%. The most prominent redshift evolution in the FWHM differential spectra appears at low luminosity, where the peak in C iv changes by 29% between redshifts z = 2.25 and z = 2.83.
While the traditional BE is often attributed to a dependence on luminosity, the spectral variations appear in almost the same fashion when varying C iv FWHM or redshift while holding luminosity constant. The similarity of these trends indicates that the BE is not due to luminosity, but instead more directly related to C iv FWHM. The appearance of the BE with luminosity is due to the simple fact that that luminosity and FWHM tend to be correlated. The fact that ∆f ∆F W HM and the ξ 1 eigenvector from Pâris et al. (2011) are so similar indicates that C iv FWHM is a larger driver of quasar diversity than any of the parameters we have explored. We have not navigated the full parameter space of quasar properties, so it is possible that the ∆f ∆log(L bol ) signature is a superposition of several effects that include ∆f ∆F W HM .

Physics in the Differential Spectra
The high correlation between the differential spectra associated with luminosity and those with redshift implies that these variations likely arise from the same physical mechanisms. In this context, the BE is actually driven by some other fundamental quasar parameters that are correlated with luminosity and evolve with redshift. Black hole mass, Eddington ratio and orientation of the accretion disk are likely physical drivers for quasar diversity, but are impossible to estimate directly from the high redshift BOSS spectra. As explained in the previous subsection, the primary driver for spectral diversity in the BOSS quasars is C iv FWHM rather than luminosity. To the extent that C iv FWHM can be used as a proxy for the mass of the black hole, these correlations can be used to infer the physical mechanisms driving the BE and redshift evolution in quasar spectra.
The C iv FWHM is known to lead to estimates of black hole mass with significant scatter in relation to masses derived from Hβ or Hα (e.g., Trakhtenbrot & Netzer 2012;Mejia-Restrepo et al. 2016). In addition, biases in the FWHM -mass relation arise due to the presence of outflowing gas which introduces asymmetry and a blueshift with respect to the systemic redshift (e.g. Sulentic et al. 2007;Netzer et al. 2007;Brotherton et al. 2015;Coatman et al. 2016). As shown in Section 3.2, we have accounted for spectral index which encapsulates some of the outflow signature of C iv. In addition, we controlled each subsample to have a large number of quasars and identical average luminosity across all redshift bins. While it is likely that some biases from outflows associated with C iv FWHM still remain, the following arguments depend only on the broad assumption that increases (decreases) in C iv FWHM at a fixed luminosity represent increased (decreased) black hole mass and decreased (increased) Eddington efficiency.
Keeping the caveats of the exact FWHM -mass relation in mind, we use the C iv line to evaluate the trends of spectral diversity in terms of changes in black hole mass, as our redshift range does not provide access to Hβ or Mg ii. The redshift evolution in average C iv FWHM at a fixed luminosity therefore represents an evolution toward larger black hole masses at lower redshifts. Because luminosity has been factored out of the analysis of redshift evolution, this trend is opposite to what might be expected for a flux-limited sample. In the ∆f ∆z differential spectra, this evolution appears as a reduction in line strengths relative to continuum at later times.
Using the C iv FWHM measurements as a tracer of the black hole mass, the ∆f ∆F W HM differential spectra can be interpreted as an anti-correlation between black hole mass and C iv equivalent width. The observed relationship is consistent with the standard thin accretion disk model (Shakura & Sunyaev 1973). The trend we observe can be attributed directly to the black hole mass (where the peak temperature of the disk scales as T max ∝ M − 1 4 ), or it can be attributed to the Eddington ratio (where the peak temperature of the disk scales as T max ∝ f 1 4 edd ). Because we have fixed the luminosity to be constant across Figure 10. The composite spectra binned by FWHM from quasars from the lowest redshift of our sample. Each panel shows the composite spectrum for the quasars with FWHM < 5000km s −1 (black) and for quasars with FWHM > 5000km s −1 (maroon). The top two panels represent selection according to C iv FWHM and the bottom two panels according to the Mg ii FWHM. Mean FWHM, luminosity and redshift are reported in the legend of each panel. the ∆f ∆F W HM and the ∆f ∆z differential spectra, it is impossible to determine whether the equivalent width is suppressed by a decrease in temperature arising from increased black hole mass or a decreased Eddington ratio.
Measurements of the intrinsic BE inform the degeneracy in our analysis between black hole mass and Eddington ratio. The intrinsic BE can be attributed to changes in Eddington ratio or other factors that are independent of black hole mass, however, it is not clear if the relationship stems from the same physical mechanism as the global BE. As with the global BE, an anti-correlation between C iv equivalent width and continuum luminosity is seen as the intrinsic BE in multi-epoch spectra of AGN. This signature could be attributed to a shielding-gas scenario proposed in Wu et al. (2011) (also see Luo et al. 2015). In this model, an increase in the accretion rate (Eddington ratio) is accompanied by an inflation of the accretion disk which reduces the number of ionizing photons reaching the C iv emitting clouds. In the scenario that the physics of the intrinsic and global BE are related, the intrinsic BE predicts a suppression of line strengths with increased Eddington ratio across large quasar samples. The fact that we observe an opposite trend suggests that the spectroscopic signature of black hole mass dominates over the spectroscopic signature of Eddington ratio. Perhaps in the global BE there is a competition between black hole mass and accretion rate, causing a shallower slope in the luminosity-equivalent width rela-tionship than what is seen in the intrinsic BE.
Following the logic above, the redshift evolution of quasar spectra at fixed luminosity carries the spectroscopic signature of the global reduction in accretion rate in supermassive black holes. Assuming a specific relation where the black hole mass at fixed luminosity scales as the square of FWHM, the average black hole mass in the low luminosity bin increases by about 40% from z = 2.77 to z = 2.25, implying a 30% global decrease in Eddington ratio at the range of masses probed in this analysis. The Horizon-AGN simulation (e.g., Dubois et al. 2014) suggests that the average 10 8 M black hole decreases by a factor of 4-5 in Eddington ratio between z = 3 and z = 2 (Volonteri et al. 2016). While the simulations cover black holes at all accretion rates rather than the population studied here, they do indicate that the quasar population is evolving quickly. This rapid change in Eddington efficiency of the quasar population in the redshift range covered by the BOSS Ly-α forest quasar sample is apparent in Volonteri et al. (2016) and provides convincing evidence that there are fundamental parameters other than luminosity responsible for quasar diversity. Figure 11. Top: Nine ∆f ∆F W HM differential spectra derived from composite spectra at fixed luminosity, spectral index, and redshift. The α mid bin is used for the spectral index with a mean of α λ = −1.52. The differential spectra computed for each luminosity bin are vertically offset by 0.001 units for clarity. Bottom: The average ∆f ∆F W HM corresponding to the highest luminosity bin, the average ∆f ∆log(L bol ) , and ξ 1 eigenvector from Pâris et al. (2011). Emission lines of interest are labeled and marked with dotted lines. A scaling has been applied such that the spectra are normalized to have the same peak in Ly α.

CONCLUSION
Using 112,000 quasars covering the redshift range 2.1 ≤ z ≤ 3.5, we sample 27 different bins spanning redshift, log(L bol ), and α λ . Each of these bins represents an observationally distinct subset of the quasar population. We shape the distributions of log(L bol ) and α λ to be nearly identical for all redshifts, resulting in a final sample of 58,656 quasars from which to construct high signalto-noise ratio composite spectra. We derive differential spectra ∆f ∆log(L bol ) , ∆f ∆α , and ∆f ∆z to describe the linear variations from one composite spectrum to another. There is a unique signature in ∆f ∆α but we focus on the signature of diversity associated with luminosity and redshift in this work. Motivated by redshift evolution in C iv FWHM, we generate differential spectra ∆f ∆F W HM for direct comparison to the variations seen with luminosity and redshift. Our primary results are: 1. A comparison of the 27 composite spectra confirms the well-documented anti-correlation between equivalent width and luminosity. Fitting the relation W λ ∝ L βw , we find a redshift evolution in the relationship such that β w = −0.35 ± 0.004, −0.35 ± 0.005, and − 0.41 ± 0.005 for C iv at z = 2.25, 2.46, and 2.84, respectively. Higher ionization species typically reflect a larger anticorrelation (e.g., C iv vs Mg ii).
2. There is a surprising redshift evolution in quasar spectra at constant luminosity. The slope of C iv equivalent width as a function of redshift (log(W λ ) ∝ zm z ) is m z = 0.16 ± 0.008, 0.11±0.004, and 0.06±0.004 for log(L bol ) = 45.95, 46.28, and 46.66, respectively.
3. The differential spectrum associated with luminosity correlates with the differential spectrum associated with redshift at 98%.
4. The average value of C iv FWHM decreases significantly with redshift for a given luminosity. The redshift evolution of quasar spectra at a fixed luminosity is directly related to the degree by which C iv FWHM changes between redshift bins. The amplitude of spectroscopic variation due to luminosity is also larger for larger spans of C iv FWHM between luminosity bins, thus leading to a redshift evolution in the slope of the luminosity -equivalent width relationship. 5. The differential spectra that describe the change in flux with respect to C iv FWHM reveal that C iv FWHM is negatively correlated with equivalent width. The C iv FWHM differential spectra show a much smaller redshift evolution than do the luminosity differential spectra.
6. The differential spectrum associated with C iv FWHM is not identical to that of the luminosity, but does correlate at the 93% level. The small differences in line features between the two differential spectra imply that there are other factors beyond FWHM that contribute to the diversity associated with luminosity. The C iv FWHM differential spectra bear more resemblance (via the correlation coefficient and comparison of the Ly-α/C iv line ratios) to the ξ 1 eigenvector derived from high redshift quasar spectra than to the luminosity differential spectra, implying that the FWHM is largely responsible for the variation in quasar spectra.
7. Under the basic assumption that C iv FWHM increases with black hole mass at a fixed luminosity, our analysis shows that the observed BE is not directly related to luminosity, but instead to a superposition of spectral diversity associated with black hole mass and Eddington efficiency. Models for AGN growth predict that average black hole mass increases and average Eddington efficiency decreases with time. The redshift evolution in the BOSS quasar sample at a given luminosity can therefore be explained by cosmic evolution in the relative contributions of C iv estimated black hole mass and mass accretion rate.
This work presents a perspective on quasar spectra that disentangles several effects responsible for quasar diversity. The basis chosen here of luminosity, spectral index, and redshift allows identification of clear redshift evolution and demonstrates that luminosity is not the direct cause for the significant variations in line features. However, the high correlation between the luminosity and redshift differential spectra indicates that the assumed basis is not optimal for capturing the full diversity in quasar spectra. We have made the 27 median composite spectra and four composite differential spectra available via download at a publicly accessible website 1 . Future work will focus on the effects of spectral index and other observable properties (such as FWHM) to isolate the more fundamental drivers of spectral diversity. The sample chosen here also covers only the redshift range sampled by the BOSS Ly-α forest quasar sample. The eBOSS program will include an even larger sample of quasars in the 0.9 < z < 2.2 redshift range (Myers et al. 2015), bringing the total sample to nearly 1,000,000 quasar spectra obtained in the four generations of SDSS. Future work will include lower redshift quasars and extended rest-frame wavelength coverage that includes the Balmer lines and the forbidden narrow lines over this much larger sample. A large number of these quasars also have multi-epoch spectra (e.g., Shen et al. 2015), allowing studies of the intrinsic variability of quasars, such as the intrinsic BE. Such studies may allow a separation