An Optical Spectrum of the Diffuse Galactic Light from BOSS and IRIS

We present a spectrum of the diffuse Galactic light (DGL) between 3700 and 10,000 A, obtained by correlating optical sky intensity with far-infrared dust emission. We use nearly 250,000 blank-sky spectra from BOSS/SDSS-III together with IRIS-reprocessed maps from the IRAS satellite. The larger sample size compared to SDSS-II results in a factor-of-two increase in signal to noise. We combine these data sets with a model for the optical/far-infrared correlation that accounts for self-absorption by dust. The spectral features of the DGL agree remarkably well with features present in stellar spectra. There is evidence for a difference in the DGL continuum between the regions covered by BOSS in the northern and southern Galactic hemisphere. We interpret the difference at red wavelengths as the result of a difference in stellar populations, with mainly old stars in both regions but a higher fraction of young stars in the south. There is also a broad excess in the southern DGL spectrum over the prediction of a simple radiative transfer model, without a clear counterpart in the north. We interpret this excess, centered at ~6500 A, as evidence for luminescence in the form of extended red emission (ERE). The observed strength of the 4000 A break indicates that at most ~7% of the dust-correlated light at 4000 A can be due to blue luminescence. Our DGL spectrum provides constraints on dust scattering and luminescence independent of measurements of extinction.


INTRODUCTION
No part of the sky is completely dark. The sky background includes airglow, scattered sunlight and moonlight, scattered light and thermal emission from dust in both the solar system and the Galaxy, and extragalactic emission from sources including unresolved galaxies and the cosmic microwave background. Each of these carries unique information. In this paper we measure the diffuse Galactic light (DGL): optical emission, both starlight and nebular emission from gas, scattered by diffuse interstellar dust.
The DGL was first detected at low Galactic latitudes as excess flux above the zodiacal light and airglow (Elvey & Roach 1937). Subsequent studies have observed the DGL from the ground (Elsässer & Haug 1960), from a sounding rocket (Wolstencroft & Rose 1966), and with data from satellites including Voyager (Murthy et al. 1991), PIONEER (Matsuoka et al. 2011), and HST (Kawara et al. 2017). Measurements of the DGL now span wavelengths from the infrared (Arendt et al. 1998;Sano et al. 2015Sano et al. , 2016 to the ultraviolet (Lillie & Witt 1976;Hurwitz et al. 1991;Seon et al. 2011).
The DGL can constrain properties of the interstellar radiation field (Murthy & Henry 1995;Witt et al. 1997), as well as properties of interstellar dust including grain size distribution and composition (Mathis 1973;Schimi-novich et al. 2001;Sujatha et al. 2005). Unfortunately, it always appears together with other sources of radiation. The DGL component can be disentangled, however, if another tracer of Galactic dust is available. Brandt & Draine (2012, hereafter BD12) measured the DGL spectrum, tracing optical background emission with blanksky calibration spectra taken by the second phase of the Sloan Digital Sky Survey (SDSS-II, York et al. 2000;Abazajian et al. 2009). To trace Galactic dust they used measurements of 100 µm intensity by the InfraRed Astronomy Satellite (IRAS, Neugebauer et al. 1984) as processed by Schlegel et al. (1998), hereafter SFD.
More data, especially optical spectra, have become available since the work of BD12. The Baryon Oscillation Spectroscopic Survey (BOSS, Dawson et al. 2013), part of SDSS-III, has observed nearly 250,000 blanksky spectra, providing new measurements of the diffuse optical background. This presents an opportunity to improve the precision of the measured DGL spectrum and to place additional constraints on dust and the interstellar radiation field. Other surveys including DESI (DESI Collaboration 2016) and GAMA (Driver et al. 2011) are similarly expanding the number of measured optical spectra.
In this paper we derive new measurements of the DGL using the combination of BOSS/SDSS-III optical spec-tra and the IRIS reprocessing of the IRAS sky maps (Miville-Deschênes & Lagache 2005). Our approach builds on that of BD12, generalizing their correlation model to account for self-absorption by dust. We structure the paper as follows. In Section 2 we describe the BOSS and IRIS data as well as the linear model of BD12 and our update to the model. Section 3 presents our DGL spectra, investigates the robustness of our model and spatial variation in the spectra, and compares to the results of BD12. In Section 4 we compare our DGL spectra to the predictions of a radiative transfer model, finding clear signatures of scattered starlight and looking for constraints on stellar populations and dust grain size distribution. In Section 5 we present evidence for extended red emission (ERE), namely an excess of intensity at ∼6500Å over the predictions of a radiative transfer model, and we estimate quantities including integrated ERE intensity, comparing to previous detections in the diffuse ISM. Then we derive an upper bound on luminescence at 4000Å. We conclude with Section 6.

DATA AND METHODOLOGY
In this section we summarize the data we use for optical spectra and far-infrared intensity, then detail our methodology for obtaining the spectrum of the DGL. The optical data are blank-sky spectra from BOSS/SDSS-III while the far-infrared data are the IRISreprocessed IRAS maps. Our methodology largely follows BD12 but with an important modification to account for self-absorption by dust.

BOSS Data
For optical light intensity, BD12 used blank-sky spectra from SDSS-II, the second phase of SDSS. In this paper we use spectra from BOSS, part of SDSS-III. The blank-sky spectra are taken for sky calibration. Of the 640 simultaneous spectra per plate in SDSS-II, typically 32 were taken from regions of the sky with no detectable sources. In contrast, the BOSS spectrographs accommodated 1000 simultaneous spectra, and BOSS observations used at least 80 blank-sky spectra per plate. Each BOSS plate has a radius of 1.5 deg and an area of 7.0 deg 2 .
The blank-sky spectra were reduced in the same way as the science spectra, with the same sky subtraction and spectrophotometric calibration. Since the calibration was intended for point sources, a correction must be applied when looking at extended sources. BD12 found that division by 1.38 is appropriate for SDSS-II, and we apply the same factor. Any difference in this value for BOSS is degenerate with the bias in our DGL spectrum (Section 2.6) and does not otherwise affect our results.
There are more data in BOSS than SDSS-II; DR12 contains 239,000 blank-sky spectra compared to the 92,000 of DR7. Although BOSS fibers are smaller than SDSS-II fibers (2 diameter instead of 3 ), the combined effect is to reduce the uncertainty in our results. As we show in Section 3, our uncertainties are dominated by sky coverage and sampling rather than formal statistical errors, and thus benefit fully from the larger number of spectra in BOSS.
BOSS also has improved sky coverage over SDSS-II, as seen in Figure 1. The upper two panels show sky fiber density for SDSS-II and BOSS in Galactic coordinates with (l, b) = (0, 0) at the center. BOSS generally has a higher sky fiber density than SDSS-II, with far more coverage in the southern hemisphere resulting in a more uniform patch. The lower two panels show fiber density weighted by the intraplate variance in 100 µm intensity, which determines the impact of fibers on our DGL spectrum. The variance is calculated for each plate using the locations of the sky fibers. From the lower panels we can see that BOSS has more fibers that contribute significantly to the spectrum.
BOSS covers a wider wavelength range than SDSS-II: 3600 to 10,400Å rather than 3800 to 9200Å. This means the BOSS uncertainties are much lower for wavelengths near the edges of the SDSS-II range. It also places the [O ii] λλ3727-30Å nebular emission doublet within the BOSS wavelength range. The resolving power of the BOSS spectrograph ranges from R = 1560 to R = 2270 for the blue channel and from R = 1850 to R = 2650 for the red channel. The wavelengths observed by SDSS-II varied slightly from night to night, so BD12 interpolated the spectra onto a common wavelength array. For BOSS, the same set of wavelengths is extracted for every spectrum, so no interpolation is necessary.
We mask all measurements from six fibers that were identified in DR12 as being affected by intermittently bad CCD columns 1 . We also perform an automated check for plates with an outsized effect on our results, looking at the impact of each plate on the average value of the DGL spectrum over 500Å intervals. Then, by visual inspection, we identify twelve plates with a nontrivial effect on the spectrum. We mask nine of these plates 2 after finding abnormalities in their optical spectra which appear to be caused by corrupted data and issues with sky subtraction. We do not mask the other three plates, but their impact is not large enough to affect our results.

100 µm Maps
BD12 used the SFD dust maps to trace far-infrared intensity; these are based on data from the IRAS and COBE (Boggess et al. 1992) satellites. We instead use the IRIS maps (Miville-Deschênes & Lagache 2005) obtained by reprocessing the IRAS and COBE images with improvements in zodiacal light subtraction, calibration, and destriping. IRIS keeps the ∼4. 3 resolution of IRAS. Our results are consistent with what we find using the Figure 1. Comparison of the sky coverage of SDSS-II (left) and BOSS (right) sky fibers in Galactic coordinates. BOSS has more uniform sky coverage, especially in the southern Galactic hemisphere, with well over twice as many sky fibers overall. The lower panels show fiber density weighted according to variance in I100 µm, which determines a plate's weight in our DGL correlation spectra. The variance is calculated for each plate using the fiber locations, and we multiply the I100 µm values by a wavelength-dependent optical depth correction factor (Equation (3)), setting λ = 7000Å for this plot. The handful of plates that we mask due to high I100 µm or abnormalities in their optical spectra are not included (Section 2.1).
SFD maps, and the uncertainties in our results decrease when the IRIS maps are used.
In Section 2.4 we introduce a new model for the relationship between optical and far-infrared intensity. This model requires the dust optical depth as a function of position and wavelength, which we obtain from the SFD maps. SFD performed additional processing to convert the 100 µm map to a map of extinction, including removing point sources and estimating the dust temperature. We convert from the SFD E(B − V ) values to wavelength-dependent optical depth using the R V = 3.1 dust model of Fitzpatrick (1999).

Linear Model from BD12
BD12 used a linear model to correlate sky fiber intensity with 100 µm intensity. Both quantities have a similar dependence on dust column density and dust temperature, as described in Section 2.2 of BD12. Briefly, scattered intensity is proportional to dust column density times the intensity of illuminating starlight, and this product is proportional to τ 100 µm T 5.5 dust if far-infrared emissivity ∝ ν 1.5 as estimated by Planck Collaboration XI (2020). Assuming thermal equilibrium and converting total emission to I 100 µm , we also have roughly I 100 µm ∝ τ 100 µm T 8 dust for dust with T ≈ 18K. If intraplate variations in DGL and I 100 µm are both primarily due to variations in dust column density, rather than variations in illuminating starlight, temperature variations will be small (Planck Collaboration Int. XLVIII 2016). Then the difference in temperature dependence, T 5.5 versus T 8 , can be neglected and the intraplate variations in DGL can be approximated as linearly dependent on I 100 µm .
The intensities from the SDSS pipeline, I λ,sky , have already been sky subtracted on a plate-by-plate basis. For consistency the mean value of (νI ν ) 100 µm , the farinfrared intensity, must also be subtracted from each plate. The model is then λI λ,sky,j,p = α λ [(νI ν ) 100 µm,j,p − (νI ν ) 100 µm p ] (1) for fiber j, plate p, and wavelength λ, where p represents an average across the fiber locations on plate p.

BD12 defines the quantity
x j,p = (νI ν ) 100 µm,j,p − (νI ν ) 100 µm p . (2) The correlation coefficient α λ and its associated uncertainty are calculated using the maximum likelihood estimator (Equations (5) and (6) of BD12). We use the term "correlation spectrum" to refer to the correlation coefficient as a function of wavelength.
This linear model is only accurate in the optically thin limit. As optical depth increases, so does the fraction of DGL photons that are absorbed by dust along our line of sight. BD12 therefore masked fibers and plates with high 100 µm intensity, which in turn implies a high dust column density. They applied a threshold of I 100 µm < 10 MJy sr −1 , corresponding to an optical depth at 5000Å of τ ≈ 0.7 for 18 K dust or τ ≈ 1.1 for 17 K dust (Hensley & Draine 2020).

Updated Model to Account for Self-Absorption
Here we generalize the linear model of BD12 to treat self-absorption along the line of sight. This results in a model that is applicable along sight lines with moderate optical depth. We introduce a correction factor, updating the definition of x from Equation (2) as follows: x λ,j,p = ((νI ν ) 100 µm β λ ) j,p − (νI ν ) 100 µm β λ p . (4) We derive the optical depth τ λ from 100 µm maps, as described in Section 2.2. Since τ λ is a function of wavelength, x becomes a function of both wavelength and sky location, while in BD12 it was a function only of sky location. Our model is then λI λ,sky,j,p = α λ x λ,j,p , and α λ can be found with the same equations as for the linear model except that x j,p is now replaced with x λ,j,p . In the optically thin limit (τ λ → 0) the correction factor becomes β λ = 1, and we recover the linear model of BD12. In the optically thick limit (τ λ → ∞) the correction factor goes to β λ = τ −1 λ . As stated in Section 2.3, 100 µm emission is roughly proportional to τ 100 µm T 8 dust , and if temperature variations are small enough, I 100 µm ∝ τ 100 µm is a good approximation. Then there are two factors of τ in Equation (4) which cancel so that x λ approaches a constant value for a given wavelength. From Equation (5), λI λ,sky also approaches a constant in this limit; only dust within the last τ ≈ 1 contributes to the visible scattered light. We expect the use of Equation (4) over (2) to have the largest impact at blue wavelengths, where dust optical depth is highest. Hβ λ4863 3.9 ± 0.3 4.5 ± 0.6 3.5 ± 0.4 To interpret α λ we consider the model definition of Equation (5). Both λI λ,sky and x λ are mean-subtracted quantities, but the proportionality constant α λ is the same if the means are not subtracted, and so we can write (in an averaged sense): Because a factor of β λ appears here, α λ is no longer proportional to λI λ . We instead choose to plot which has a straightforward interpretation as intensity and is therefore better for comparing to model spectra. We employ the same masking as BD12, with a threshold of 10 MJy sr −1 except where indicated otherwise. We mask based on SFD 100 µm values for the sake of comparison; the IRIS values are generally larger, so if we were to use them for masking we would mask more fibers. The details of masking should not impact our results because, as we will show in Section 2.4, the updated model is insensitive to the masking threshold.
For most analyses we increase the signal-to-noise ratio by binning the correlation spectra with a bin width of 50Å, and we mask each nebular emission line listed in Table 1 by removing the data in a region around the central wavelength corresponding to a velocity range of ±200 km s −1 .

Bootstrapping
BD12 obtained uncertainties on their DGL spectrum using the errors reported by the SDSS pipeline and the formal uncertainties of a linear fit. This will underestimate the true uncertainty if there are any correlated or non-Gaussian errors, and the effect is amplified when data are binned. To overcome this limitation we use bootstrap resampling ("bootstrapping") to compute the uncertainties on our DGL spectra.
Bootstrapping is a technique where many samples of size N are taken from a data set of size N , with replacement. This approximates sampling from the true underlying distribution responsible for the data. A statistic is computed for each sample and the variation in the statistic across samples can be used to compute uncertainties.
We use bootstrapping to generate 10,000 realizations of the sample of BOSS sky spectra. Each of these realizations draws from the sample of observed plates rather than fibers. We treat separate observations from the same physical plate (on different nights) as distinct. In this way we approximate repeating the survey many times, each time with a different placement of the BOSS plates. Then we compute a correlation spectrum for each of the 10,000 bootstrap samples and find the 16 th and 84 th percentiles of α λ at each wavelength.
With the exception of Figure 2, where we also include formal uncertainties, errors shown in figures are half of the 68% confidence intervals derived from bootstrapping. This directly corresponds to the standard deviation if the errors are Gaussian. In some figures we plot the 16 th and 84 th percentiles to form an envelope. For binned plots we first bin the correlation spectra before taking percentiles.

Bias Factor
Based on a comparison to the predicted results from stellar models and a radiative transfer model, it appears that α λ is biased low by close to a factor of two. In Section 4.3 we describe this comparison, and we adopt bias factors of C = 2.1 for the northern Galactic hemisphere and C = 1.9 for the southern hemisphere.
BD12 estimated a 20% uncertainty in the bias factor, based primarily on how the predicted DGL spectrum changes when dust optical depth is varied in the radiative transfer model and when a model of the local ISRF is used instead of stellar population synthesis models. We adopt the same 20% uncertainty, resulting in C = 2.1 ± 0.4 in the north and C = 1.9 ± 0.4 in the south. There is a 13.5% gain uncertainty in the IRIS 100 µm map (Miville-Deschênes & Lagache 2005), and while this affects the scaling of the measured DGL spectrum, it does not change the percent uncertainty in the bias factor. Any error in our flux calibration factor (Section 2.1) would have a similar effect.
An estimate of the bias factor should be applied whenever the true ratio of λI λ to (νI ν ) 100 µm is desired. We apply bias factors for the ERE calculations in Section 5.2, where we are concerned with finding actual intensity values and comparing to previous work. Be-cause of the high uncertainty in the bias factor, we do not apply it in plots except where explicitly indicated by a factor of C.
The source of the bias is not clear. Here we investigate two possibilities: the effect of neglecting 100 µm noise in our model and the difference in resolution between the optical and 100 µm data. We find that neither effect is large enough to explain the observed bias.
Section 4.2 of BD12 describes how uncertainties in the 100 µm map, which are not considered in our model, can lead to a bias. However, we estimate the size of this effect and find that the statistical uncertainties in the IRIS map are too small to explain a factor-of-two bias. Equation (14) of BD12 provides a way to estimate the bias if the measurement uncertainty and the variance in 100 µm intensity are both known. While the IRIS maps do not include uncertainty maps, Miville-Deschênes & Lagache (2005) estimates the median noise level to be 0.06 ± 0.02 MJy sr −1 . We also find that the median variance in 100 µm intensity across the BOSS plates is 0.06 MJy 2 sr −2 . The predicted bias factor is then 1.07 ± 0.04.
Another possible source of bias is the difference in size between the BOSS fibers and the IRAS pixels. The ∼4. 3 IRAS pixels are much larger than the 2 BOSS fibers, so there is some error in determining the 100 µm intensity at the location of a BOSS fiber. However, we apply Equation (14) of BD12 to find that the bias from this effect is also small. Projected dust emission closely follows a P (k) ∼ k −2.9 power law down to 15 (Miville-Deschênes et al. 2016), which means that the power on small scales falls off sharply. H i observations suggest that this k −2.9 power law continues down to 5 (Roy et al. 2010). We integrate P (k)k dk to estimate the contribution of different spatial scales to the variance. Assuming that the k −2.9 power law extends down to 2 , scales from 4. 3 to 3 • contribute about 25 times as much variance as scales between 2 and 4. 3; this alone would result in a bias factor 1.05.
Since random errors in the 100 µm map and the resolution mismatch both predict bias factors much smaller than two, something else must be responsible for the bias. Any deviations from our simple model relating optical intensity to 100 µm emission will contribute to the bias, and we suspect this is the primary cause. These deviations include physical sources of spatial variation such as contamination by zodiacal light and associated 100 µm emission, variations in starlight intensity which affect dust grain temperature and therefore 100 µm emission, and variations in the anisotropy or spectrum of incident starlight. In addition, there might be variations in the physical properties of dust grains arising from processes such as grain growth.
BD12 demonstrated that the bias factor is independent of wavelength for a linear model. Our nonlinear model introduces a possible wavelength dependence, but the nonlinearity correction is 20% at all wavelengths (Section 3.2).

INITIAL RESULTS
In this section we present DGL correlation spectra computed using the approach of Section 2.4. We describe the differences between our results and those of BD12 due to new data and an updated model. Then we show that our results are insensitive to the masking threshold for plates and fibers with high dust optical depth, and we explore the sensitivity of our results to the spatial footprint of the survey. Finally, we measure the strengths of nebular emission lines. Figure 2 shows spectra of the DGL computed using the approach described in Section 2.4, with uncertainties from bootstrapping (Section 2.5). The spectrum in the left panel is based on data from SDSS-II, and the right panel shows the new spectrum using data from BOSS.

BOSS Results
The SDSS-II and BOSS correlation spectra, overall, are very similar. However, the BOSS spectrum is slightly more peaked near 6500Å, and the optical intensity relative to 100 µm emission appears to be higher for the BOSS spectrum at most wavelengths. We show in Section 3.3 that these differences can be explained by sky coverage, and that some of the difference may be due to spatial variation in the uncertain bias factor (Section 2.6) rather than a true difference in emission.
The BOSS DGL spectrum represents a substantial gain in precision over the spectrum from SDSS-II, largely due to the increased number of optical spectra. Between 5000Å and 8000Å, with nebular emission lines masked and the spectra binned into 50Å bins, the median signal-to-noise ratio is 24 for BOSS compared to 12 for SDSS-II. This difference is even more apparent near the edges of the SDSS-II wavelength range because BOSS covers a wider range.
Bootstrap uncertainties in the binned spectra are larger than the formal uncertainties calculated from errors given by the SDSS-II and BOSS pipelines. The dominance of bootstrap uncertainties makes the larger sample size and more uniform spatial sampling of BOSS especially important.

Updated Model Results
In Section 2.4 we describe a new model for the DGL that accounts for dust self-absorption. Figure 3 shows that this updated model does indeed reduce the need to mask regions of high dust optical depth, as it makes our results insensitive to the choice of masking threshold. The left panel adopts the linear model of Section 2.3 with a range of masking thresholds, while the right panel applies the same thresholds for the nonlinear model of Section 2.4.
For low 100 µm masking thresholds (e.g. 10 and 15 MJy sr −1 ), the optical depth is small enough that the linear model is a reasonable approximation; the updated model gives similar results. As expected, the biggest correction is at the blue end of the spectrum where there is more self-absorption. The correction at these wavelengths is ∼20%, with the linear model underestimating the DGL intensity due to its neglect of self-absorption.
As the masking threshold increases, fibers with higher optical depth are included and effects from absorption become more significant. The results of the linear model change drastically, indicating that the model is breaking down. In addition, the newly-included plates have a large impact on the spectrum, increasing the bootstrap uncertainties. On the other hand, the results from the nonlinear model remain virtually identical even with a masking threshold of 25 MJy sr −1 . At this intensity the optical depth for 18 K dust ranges from τ ≈ 2 at 4000Å to τ ≈ 0.5 at 10,000Å.

Spatial Variation
The DGL spectrum is expected to vary across the sky due to differences in the interstellar radiation field and/or dust properties. In Figure 4 we plot DGL spectra for the footprints of BOSS in the northern (b > 0) and southern (b < 0) Galactic hemispheres. Subdividing the data further quickly results in spectra with bootstrap uncertainties too large for detecting spatial variations.
This division is a natural choice for comparisons with SDSS-II, which observed a similar region to BOSS in the north but had limited and patchy sky coverage in the south. In fact, the SDSS-II correlation spectrum closely matches the north BOSS spectrum, suggesting that the difference seen in Figure 2 between the SDSS-II and BOSS spectra can be attributed to sky coverage.
While the difference between the north and south spectra in Figure 4 initially appears highly significant, we find that much of it can be explained by spatial variation in the bias factor (Section 2.6). In Section 4.3 we adopt bias factors of C = 2.1 in the north and C = 1.9 in the south; when this scaling is applied, the north and south spectra agree well with each other and with plausible stellar models over the range 4200Å to 5000Å. For wavelengths above 5500Å, there are significant differences in the spectra that cannot be easily resolved with a scaling factor, and we explore these in Sections 4 and 5. Possible explanations include differences in stellar populations, dust properties, and ERE.
As in BD12, we observe significant variation in the DGL spectrum with Galactic latitude, but there does not appear to be a clear trend in the latitude dependence. Several effects could contribute including variation in dust properties and variation in the ISRF. Especially at high latitudes, where the signal is weaker, there may also be contamination by extragalactic light (Yahata et al. 2007) and systematics from an imperfect model. It is difficult to separate these effects without a more detailed model of the galaxy.  . Left panel: the linear correlation model of BD12 fails at high optical depth due to self-absorption by dust; the DGL spectrum depends on the 100 µm intensity threshold (given in MJy sr −1 ) used to mask plates and fibers. Right panel: our updated correlation model for the DGL (Section 2.4) accounts for self-absorption by dust and is valid in regions of modest optical depth. The resulting spectrum is insensitive to the threshold used to mask fibers. We plot α λ in the left panel, but in the right panel we plot α λ (Equation (7)), which is more appropriate for comparison. The dashed lines are half of 68% confidence intervals obtained from bootstrapping (Section 2.5). The data are binned with a bin width of 50Å and nebular emission lines are masked.  Figure 1. Much of the difference in the spectra seen here could be explained by a difference in the uncertain bias factor (Section 2.6), but there may also be signs of differences in stellar populations and ERE (Sections 4 and 5). The envelopes around the spectra are 68% confidence intervals found with bootstrapping (Section 2.5), and the errors plotted with dashed lines are half of this interval. The data are binned with a bin width of 50Å and nebular emission lines are masked.

Emission Lines
Peaks from nebular emission are clearly visible in Figure 5; they represent light emitted by interstellar gas and scattered by dust grains. We mask these lines in our analysis and in most plots, but the line strengths can be useful in their own right as tracers of ISM properties; they carry signatures of the characteristic densities, temperatures, and ionization states of the nearby ISM. Table 1 reports the equivalent widths of nebular emission lines for the full sky and for the northern and southern Galactic hemisphere. The listed uncertainties are half of the 68% confidence intervals found with bootstrapping (Section 2.5).
We calculate the equivalent widths by fitting Gaussians, assuming a flat continuum and using the maximum likelihood estimator to simultaneously fit the continuum level and height of each peak. We center the fits on the known rest wavelengths of the lines and assume widths corresponding to instrumental broadening from the BOSS spectrograph. This gives a width for the Hα line that is in good agreement with the value we find by including the width as an additional parameter in the fit and performing a 1D nonlinear optimization.
The one exception to this is the [O ii] λλ3727-30Å doublet, which is outside the SDSS-II wavelength range but inside the range of BOSS. In principle the relative line strength can serve as an indicator of electron density (Draine 2011), however the noise level is too high to provide meaningful constraints. Instead we assume a ratio of [O ii]λ3730/[O ii]λ3727 = 1.5, corresponding to the limit of low electron density, and fit a sum of two Gaussians. This limit is implied by the values we find for [S ii]λ6718/[S ii]λ6733, and it is consistent with our expectations for the Warm Ionized Medium (Madsen et al. 2006) and for low-density H ii regions such as the one around the nearest O-type star, ζ Oph (Maíz- Apellániz et al. 2004).
Hα and Hβ are present in absorption in starlight, so our fit underestimates their strength in emission. We correct for this using Equations (8) and (9) of BD12 which are approximate linear relations between the size of the 4000Å break and the equivalent widths of the Hα and Hβ absorption lines.
The equivalent widths are generally larger in the north, which appears to be due to a greater share of the interstellar radiation field coming from gas emission. The line ratios in the north and south are similar, suggesting similar physical ISM states.

STARLIGHT AND DUST
The interstellar radiation field at optical wavelengths is dominated by starlight, so we expect the DGL to be mainly composed of scattered starlight. Here we confirm this by comparing our measured DGL correlation spectra to model stellar spectra. We then attempt to determine the properties of the stars that are responsible and to provide some constraints on the properties of dust grains. Because we are concerned here with starlight, we mask wavelengths where we expect nebular emission lines from interstellar gas, as described in Section 2.4.

Models
The stellar population in the solar neighborhood covers a range of ages, but with a higher proportion of old stars than young stars (Casagrande et al. 2011). Our starlight templates are composite stellar populations from Bruzual & Charlot (2003, hereafter BC03): exponentially declining star formation over 12 Gyr with timescales of 5 Gyr and 9 Gyr ("t5e9" and "t9e9") and constant star formation over 6 Gyr ("cst"). We consider linear combinations of these model spectra with metallicity Z = 0.02.
To account for scattering and extinction by dust, we apply the same radiative transfer model described in Section 5.1 of BD12, summarized here. We assume an infinite plane-parallel galaxy and a sight line with latitude b = 40 • . The dust distribution is Gaussian: ρ dust = exp −z 2 /2σ 2 with σ = 250 pc. We assume a V -band dust optical depth of τ V = 0.15 csc b. There are two exponential stellar distributions with scale heights of 300 and 1350 pc, and the 300 pc component contains 90% of the stars. For scattering we use a Henyey-  Table 1 and are generally larger in the north due to a greater contribution from ISM light. The errors plotted with dashed lines are half of the 68% confidence intervals found with bootstrapping (Section 2.5).
Greenstein phase function. The total absorbed power is taken as an estimate of the total power radiated in the IR, which is then converted to the 100 µm bandpass using the Draine & Li (2007) dust model. The relation is (νI ν ) 100 µm ≈ 0.52I TIR . This allows us to find the ratio of scattered light to 100 µm emission, λI λ / (νI ν ) 100 µm , which can be compared to α λ (Equations (6) and (7)). The implementation of this radiative transfer model in BD12 included two errors which we correct. These turned out to be largely compensating errors; our final results are very similar to those of BD12. We consider dust models from Zubko et al. (2004) and Weingartner & Draine (2001), hereafter ZDA04 and WD01. Both are "bare" models consisting only of PAHs, graphite grains, and silicate grains. Values for albedo, wavelength-dependent cross section, and anisotropy parameter (g = cos θ) are taken from these models. The WD01 model contains more large grains, resulting in a redder spectrum as shown in Figure 3 of BD12.

Evidence of Starlight
We expect the DGL spectrum to be dominated by scattered starlight, and we see convincing evidence that this is the case, both in the spectral features of binned plots and in a more detailed comparison to model spectra. This serves as a check that our method is working properly and suggests that we should be able to determine properties of stars and dust grains by comparing to model spectra. In the right panel of Figure 2 (and in other binned plots) there are clear indicators of starlight, including the 4000Å break and absorption near 5170Å.
In Figure 6 we plot the correlation spectrum at the native BOSS resolution alongside a model from BC03 (exponentially declining star formation over 12 Gyr with a timescale of 5 Gyr, Z = .02). We divide each spectrum by a high-order polynomial to remove the continuum, using least squares to fit the polynomial coefficients. This takes into account the higher uncertainties toward the edges of the wavelength range. While not physically motivated, we find that a 10 th order polynomial is sufficient to remove the power from the continuum.
The observations agree with the model of scattered starlight better than with a flat continuum (χ 2 ν = 0.60 versus 0.69). The upper panel of Figure 6 (400Å wide) shows that the observed DGL spectrum matches the features of the stellar model quite well, down to the smallest scale resolved by the model. The lower panel (2000Å wide) shows that the DGL spectrum also closely follows larger-scale variations in the stellar model.
The resolution of the BC03 model spectrum corresponds to a velocity dispersion of roughly σ v = 70 km s −1 . We smooth with a Gaussian kernel corresponding to σ v = 90 km s −1 , which is the value that minimizes the χ 2 between the observed spectrum and the model. Adding 70 km s −1 and 90 km s −1 in quadrature results in an effective velocity dispersion of σ v ≈ 115 km s −1 , consistent with the resolution of the BOSS spectrograph and somewhat larger than the velocity dispersion of disk stars (Holmberg et al. 2007).

Fitting Procedure
What types of stars are responsible for the DGL spectrum? We consider linear combinations of three models: exponentially declining star formation over 12 Gyr with timescales of 5 Gyr and 9 Gyr ("t5e9" and "t9e9") and constant star formation over 6 Gyr ("cst"). Although we attempt to quantitatively constrain stellar populations, the results are inconclusive and so we perform a heuristic comparison between the observed spectra and stellar models. We first try simultaneously fitting several Lick indices which BC03 found to be sensitive to age and metallicity (Section 4.3 of BC03). Unfortunately, some are Balmer series indices which we are unable to use because we cannot accurately separate the contributions from nebular emission. The fits from all three spectra are similar enough that we are unable to determine a clear best-fit combination. We also try a fit to the whole spectrum after continuum normalization but are again unable to draw a meaningful conclusion.
Our best constraints, then, come from fitting by eye. The left panel of Figure 7 shows fits to the north and south BOSS spectra, determined in Section 4.4 assuming ZDA04 dust. We primarily consider wavelengths toward the edges of the range when evaluating fits because of a possible contribution from ERE (Section 5).
To estimate the uncertain bias factor (Section 2.6), we first assume that the predicted spectrum derived from the t5e9 model represents the true level of the DGL spectrum. Then we scale the observed spectra to match the average level of the t5e9 spectrum over the range 4200-5000Å. This suggests bias factors of C = 2.1 in the north and C = 1.9 in the south, in agreement with the value of C = 2.1 ± 0.4 adopted by BD12.

Stellar Populations in the North and South
Here we present our best-guess stellar populations based on the fitting method described in the previous section. For the region covered by BOSS in the northern Galactic hemisphere, we prefer the t5e9 model; it matches the observed DGL spectrum reasonably well when the ZDA04 dust model is used, as shown in the left panel of Figure 7. The t5e9 model fits better than either t9e9 (more gradually declining star formation) or cst (constant star formation). Both of these models, when combined with ZDA04 dust, predict the DGL to be bluer than observed.
The DGL spectrum in the southern hemisphere appears to be slightly bluer than the spectrum in the northern hemisphere. The left panel of Figure 7 shows that, with the scaling described in Section 4.3, the t5e9 model slightly overpredicts the southern DGL at the reddest wavelengths. There are several effects that could account for this, including random variation or differences in metallicity, dust properties, or stellar populations. ERE is expected to cause a broad peak, as we discuss in Section 5, but it should not cause a tilt.   Figure 7. Left panel: measured DGL correlation spectra, together with predictions derived from BC03 starlight models using a simple radiative transfer model (Section 4.1) and the ZDA04 dust model. Our choice of starlight models is described in Section 4.4. We scale all spectra to a common level over the wavelength range 4200-5000Å, which corresponds to applying bias factors (Section 2.6) of C = 2.1 in the north and C = 1.9 in the south. The southern DGL spectrum shows a clear excess over the models at ∼6500Å, and we interpret this as ERE, while the excess in the north is less clear. Right panel: we plot the excess directly to visualize the ERE peak. The envelopes around the spectra are 68% confidence intervals found with bootstrapping (Section 2.5). The data are binned with a bin width of 50Å and nebular emission lines are masked.
A search of Gaia EDR3 (Gaia Collaboration et al. 2021) suggests that the southern region of the BOSS survey contains a higher fraction of young stars within 500 pc. As a proxy for young stars, we search for stars that are bluer than B p − R p = 0.5 and have absolute magnitudes brighter than G = −1. The approximate northern and southern footprint of BOSS each contain 18 such stars within 500 pc of the Sun. However, the northern footprint is much larger, with well over twice as many stars overall. Therefore a difference in stellar populations is a plausible explanation.
We see in the left panel of Figure 7 that adding constant star formation makes the spectrum bluer and improves agreement with the measured DGL spectrum in the south. We add enough constant star formation to match the level of the south spectrum at 9000Å. Our preferred model in the south is then a linear combination of the t5e9 and cst models, with ∼40% of the present-day star formation due to the cst population.

Dust Models and Systematics
The synthetic spectra in Figure 7 are computed using the ZDA04 dust model. Here we consider the WD01 dust model, which includes more large grains and yields a redder scattered spectrum.
The ZDA04 model results in a slightly better fit, particularly in the south. The same stellar models with the WD01 dust model overpredict the DGL spectrum at red wavelengths. The level of the south spectrum at 9000Å can still be matched by adding more constant star formation, with this population contributing ∼60% of the present-day star formation. The resulting spectral shape, however, does not agree quite as well.
Unfortunately, the bootstrap uncertainties and the limited nature of our radiative transfer calculations prevent us from drawing strong conclusions. There are possible systematics in our model (as defined in Section 4.1), including effects from our choice of model parameters, which could change our dust model preference. This would have only a minor impact on our conclusions about stellar populations, with both the north and south spectra remaining best modeled with mostly old stars, but varying the spectral shape could tip the balance in favor of the WD01 dust model.
Changing the metallicity from Z = 0.02 to Z = 0.008 results in a 9% difference in the predicted correlation spectrum at 9000Å. Changing the value of τ V at b = 90 • from 0.15 to 0.05, a plausible value based on Figure 11 of BD12, results in a difference of 22%. These changes both lower the spectrum at red wavelengths (with the scaling described above, where the blue end is held fixed), and they are both comparable to the 15% difference caused by switching to the WD01 dust model. Other parameter choices have a smaller impact on the spectrum. These include changing the latitude of the sight line from b = 40 • to 20 • or 60 • and removing the thick-disk stellar distribution with a scale height of 1350 pc.

ERE DETECTION
Our DGL spectra, particularly in the southern Galactic hemisphere, show a broad excess over stellar models that is centered near 6500Å. A wide range of previ-ous works have found a similar excess, known as ERE, in reflection nebulae (e.g. Schmidt et al. 1980;Witt et al. 1984;Witt & Boroson 1990) and the diffuse ISM (e.g. Gordon et al. (1998, hereafter G98); Szomoru & Guhathakurta (1998, hereafter SG98)). The carrier responsible is unknown; see Witt & Lai (2020) for a recent review of constraints on ERE models. Observations of a variety of objects and environments enable comparisons that can help identify the ERE carrier.
The first measurements of ERE were made in reflection nebulae, and due to similarities in dust properties, it was predicted that ERE would also be found in the diffuse ISM of the Milky Way (Witt & Lai 2020). Since then there have been several detections. Most have relied on broad-band measurements of sky intensity in a small number of bands, including G98 and Witt et al. (2008), which looked at the diffuse ISM. SG98 made spectroscopic measurements of three high-latitude dust clouds. Our spectroscopy of the diffuse ISM complements these nicely.

Our Detection
Here we describe our detection of ERE. In the left panel of Figure 7 we plot our DGL spectra for the regions covered by BOSS in the northern and southern Galactic hemisphere, with estimates of the uncertain bias factor (Section 2.6) applied. The south spectrum appears slightly more peaked, suggesting the presence of ERE. In Section 4.4 we identified various effects that can introduce a tilt in the predicted DGL spectrum: differences in stellar populations, dust properties, and metallicity. However, none of these can reproduce a broad peak.
Next we show that a comparison to models of scattered starlight provides further evidence for ERE. Unlike in the case of reflection nebulae, we cannot directly measure the spectrum incident on the dust, so we instead make use of BC03 model spectra. We apply the radiative transfer model from Section 4.1 and fit by eye at ∼4500Å and 9000Å as described in Section 4.3. In the left panel of Figure 7 we plot our preferred model spectra (Section 4.4) together with the observed spectra. The right panel shows the excess in the observed spectra compared to the model predictions; this excess serves as our estimate of ERE.
For the DGL spectrum in the northern BOSS region, there does appear to be some excess at ∼6500Å, but the noise level is too high to claim a significant detection. For the spectrum in the southern region, on the other hand, there is a significant peak in the excess compared to our best-guess model (green curve; a combination of exponential and constant star formation) over the range 5000-8000Å. We also plot the peak for only exponential star formation (gold curve), but there is little difference in the intensity or shape of the inferred ERE spectrum. The higher significance of our ERE detection in the south is partly because the spectrum is more peaked, and partly because of smaller uncertain-ties. This spatial difference may explain why BD12 did not claim a detection of ERE: SDSS-II sampled mostly from the northern hemisphere.
In Figure 7 we include bootstrap envelopes for the observed spectra, but we neglect model systematics, which also contribute to the uncertainty. The width and size of the ERE peak depend to some extent on choices we make about scaling, the dust model, and the mix of stellar populations. The model fit is quite good at blue wavelengths (4200-5000Å), but the rest of the systematics we have considered (Sections 4.4 and 4.5) affect the level of the spectrum at red wavelengths. Any combination of dust model and stellar population that matches observations at the red and blue end underpredicts the spectrum in the middle, suggesting an ERE peak. It is possible that ERE could extend to longer wavelengths, in which case we would be underpredicting the ERE by choosing a model that matches the DGL spectrum at ∼4500Å and 9000Å.
The plots in Figure 7 use the ZDA04 dust model. If the WD01 model is used instead, the t5e9 spectrum predicts a smaller ERE peak but also overpredicts the DGL intensity at red wavelengths. When constant star formation is added to match the level of observations at the red end of the spectrum, the ERE estimate for WD01 dust is almost the same as for ZDA04 dust.
Previous studies have suggested that the ERE carrier is excited primarily by UV photons (e.g. Darbon et al. 1999;Witt et al. 2006;Lai et al. 2017). Our results are consistent with this, since we find stronger evidence for ERE in the southern region of the BOSS survey. As discussed in Section 4.4, Gaia EDR3 identifies a higher proportion of luminous (absolute G < −1 mag), relatively blue (B P − R P < 0.5 mag) stars within the southern BOSS footprint.

Numerical Estimates
In this section we quantify our ERE detection by subtracting the predictions of our radiative transfer model from the observed DGL spectra, using our preferred stellar populations and dust model.
Our correlation spectra are biased low (Section 2.6), and for the purpose of these calculations we multiply by factors of C = 2.1 in the north and C = 1.9 in the south (Section 4.3) to correct for this. These factors are important for our estimate of the integrated ERE intensity, but they do not affect the significance of our detection or the ratio of ERE intensity to the intensity of scattered starlight. BD12 estimated a 20% uncertainty in their adopted bias factor of C = 2.1, and we expect this to be one of the dominant sources of uncertainty except where the bias factor cancels. Bootstrap uncertainties are comparable: we estimate ∼15% in the south for the integrated ERE intensity.
Following the procedure of Witt & Boroson (1990), we divide the peak into quartiles of equal integrated intensity, then report the 2 nd quartile as the peak and the difference between the 1 st and 3 rd quartiles as the width. We find a peak of 6500Å and a width of 1100Å. To get the integrated ERE intensity in units of flux per solid angle, we rearrange Equation (6) and multiply by the bias factor C to find that the optical intensity is given by and for I 100 µm we take a mean over the unmasked BOSS sky fibers. Integrating the excess in I λ over the prediction of our preferred BC03 model gives I ERE ≈ 0.28×10 −5 erg s −1 cm −2 sr −1 in the south (with I 100 µm = 2.6 MJy sr −1 ). In the north we find I ERE ≈ 0.04 × 10 −5 erg s −1 cm −2 sr −1 (with I 100 µm = 1.4 MJy sr −1 ), although we reiterate that there is not a clear detection in the north. Next we compute the ratio of integrated ERE intensity to the intensity of scattered starlight over the range 5000-8000Å, roughly corresponding to the photometric R band. This range is chosen because previous works, both photometric and spectroscopic, have reported a similar quantity. We obtain I ERE /I sca ≈ 0.20 for the south and 0.06 for the north. We also compute a ratio of I ERE to total infrared (TIR) intensity, converting from (νI ν ) 100 µm to I TIR using the Draine & Li (2007) dust model; BD12 estimates a ∼10% uncertainty in this conversion. The result is I ERE /I TIR ≈ 0.018 in the south and 0.005 in the north, suggesting that ∼1% of the total power incident on the grains is radiated as luminescence.
Calculating the ratio of I ERE to the UV power incident on the dust allows an estimate of the energy conversion efficiency. We take γ = I UV /I TIR , roughly equal to the ratio of UV to total incident power. If ERE is excited by photons more energetic than 6 eV or 8 eV, we find using our best-guess BC03 model that γ ≈ 0.56 or 0.46 in the south, corresponding to 3% or 4% of the incident UV power converted to ERE.

Comparisons
Our DGL spectrum provides a measurement of ERE from the diffuse ISM, averaged over large regions of the sky. In this section we compare our ERE constraint with other measurements over smaller spatial regions.
A measurement of ERE requires subtracting the contribution of scattered light, and usually a variety of other backgrounds including airglow and zodiacal light (e.g. SG98). G98 avoided airglow and zodiacal light by using measurements from spacecraft outside the zodiacal dust cloud, but they still had to remove light from stars and galaxies. Our method has the advantage that we only need to remove the contribution of scattered starlight as described in Section 5.1. Airglow does not correlate with 100 µm intensity, although the SDSS spectra have been sky subtracted, which further helps to isolate the DGL. Zodiacal light has been explicitly removed from the SFD and IRIS dust maps. It is possible that residuals from the zodiacal light removal may be large at high latitudes (Lauer et al. 2021), but we find that masking data within 10 • of the ecliptic has a negligible effect on the DGL spectrum. The optical spectra we use have also been placed in regions of the sky without detected optical sources.
The moderate resolution and high signal-to-noise ratio of our DGL spectrum enable a relatively clean detection of the ERE. Our peak wavelength of ∼6500Å is higher than the value of 6000Å from SG98, but still consistent with the observed trend of ERE peak wavelength increasing with radiation field density (Smith & Witt 2002).
A comparison of our integrated ERE intensity with literature results is difficult. The bootstrap uncertainty is ∼15% in the south, and we multiply by a bias factor with an uncertainty of ∼20% (Section 2.6; we adopt C = 2.1 in the north and C = 1.9 in the south). In addition, our method determines a spatially averaged ratio of ERE intensity to 100 µm emission, leaving open the possibility of spatial variation in this quantity. With those caveats, the value of I ERE ≈ 1.2 × 10 −5 erg s −1 cm −2 sr −1 from SG98 is larger than our inferred value for the southern BOSS footprint by a factor of four, and much larger than our value for the northern footprint.
If we instead consider the ratios I ERE /I sca and I ERE /I TIR , as defined in Section 5.2, we can avoid some of the problems mentioned above. There should be less spatial variation in these quantities because ERE, optical intensity, and infrared intensity all depend on dust column density. In the case of I ERE /I sca , the uncertain bias factor cancels. This ratio is known to depend on scattering geometry (G98), which can complicate interpretation, but the values here are all from observations taken primarily at Galactic latitudes of b 25 • .
Our value of I ERE /I sca ≈ 0.20 is smaller than the value of ∼0.3 from SG98, although it falls in the range 0.05-2 reported by G98 who had much better spatial resolution. G98 estimates that I ERE /I TIR ≈ 0.03. For the filaments observed by SG98, the 100 µm intensities they report also imply I ERE /I TIR ≈ 0.03, although they note that zodiacal light has not been subtracted which could increase this fraction. Our determinations of I ERE /I TIR ≈ 0.018 in the south and 0.005 in the north are smaller than the values found by G98 and SG98 unless the correction for bias in determining our spectrum is as large as C = 3.

Blue Luminescence
Our DGL spectrum can constrain Blue Luminescence (BL), an emission process that occurs at blue wavelengths (Vijh et al. 2005). This constraint exploits the presence of stellar features in the DGL. The size of the 4000Å break, a prominent feature of old stellar populations, should be preserved by scattering but diluted if there is any luminescence.
In Figure 8 we plot the 4000Å break for the observed DGL spectra in the north and south and for stellar mod- . Measured DGL correlation spectra around the 4000Å break, plotted together with predictions derived from BD12 starlight models using a simple radiative transfer model (Section 4.1) and the ZDA04 dust model. Our choice of starlight models is described in Section 4.4. The wavelength range 3850-4150Å is used to quantify the break in Section 5.4; the break is stronger than expected based on these models. We scale all spectra to a common level over the wavelength range 4200-5000Å, which corresponds to applying bias factors (Section 2.6) of C = 2.1 in the north and C = 1.9 in the south. The observations are binned with a bin width of 10Å, while the model spectra are not binned. els; high noise levels and spectral features make visual comparisons of the strength of the break difficult. We therefore define δ 4000 as the ratio of I λ just below the break (3850 to 4000Å) to I λ just above the break (4000 to 4150Å), finding δ 4000 = 0.73 ± 0.04 in the south, 0.71 ± 0.05 in the north, and 0.72 ± 0.03 for the full sky.
The t5e9 BC03 spectrum and ZDA04 dust model predict δ 4000 = 0.80 for the scattered spectrum; Table 2 reports values of δ 4000 for other model combinations. Thus the observed 4000Å break is stronger than expected, even compared to old stellar population models. This is evidence against significant luminescence at 4000Å. If a stellar model correctly describes the 4000Å a Our best-fit models described in Section 4 b Contributes ∼40% of star formation c Contributes ∼60% of star formation break and the contribution from BL is constant in I λ , then the fraction f BL contributed by BL to the total dust-correlated radiation at 4000-4150Å is given by With the full-sky data we use bootstrapping to find a 99.7% confidence interval for δ 4000 of [0.62, 0.81], and from there we can place a 3σ upper limit of f BL = 0.07, assuming the t5e9 stellar model and ZDA04 dust. If the stellar population contains more young stars, this bound decreases.
6. CONCLUSION In this paper we have presented a new measurement of the spectrum of the DGL. This work represents a substantial improvement over the DGL spectrum derived by BD12 using a similar approach. BOSS provides more data with better and more uniform sky coverage, resulting in a spectrum with lower uncertainty, and our model is more robust, applicable in optically thick regions with moderate self-absorption.
The DGL spectrum appears to be largely consistent with a simple radiative transfer model for scattering by dust of starlight from the local stellar population. It reproduces detailed spectral features of starlight, allowing us to infer properties of stellar populations and dust grains. Comparing results from the regions covered by BOSS in the northern and southern Galactic hemisphere, there is a possible difference at the red end of the spectrum, consistent with a higher proportion of young stars in the southern region. It is difficult to constrain the dust model, as stronger scattering at red wavelengths may be compensated by a younger, bluer stellar population. With this caveat, we have a slight preference for the ZDA04 dust model over the WD01 model. The former has fewer large grains and weaker scattering at the red end of the spectrum.
Comparisons to model spectra also reveal a broad excess in the DGL spectrum centered near 6500Å, especially in the southern hemisphere. A natural explanation appears to be ERE, supporting previous measurements that have detected ERE in the diffuse ISM. ERE is believed to result from luminescence of dust grains excited by UV photons, and our stronger detection of ERE in the south is consistent with this; Gaia detects a higher fraction of luminous young stars within 500 pc in the southern BOSS footprint relative to the northern footprint. The observed 4000Å break is consistent with scattered starlight, with no evidence for BL by dust; we obtain an upper limit of ∼7% on the fraction of dustcorrelated light at 4000Å contributed by BL.
The large number of sky spectra and the extensive, relatively uniform sky coverage of BOSS provide a powerful new probe of the DGL. Future analysis could improve the DGL spectrum even further by adding data from surveys including eBOSS (Dawson et al. 2016), the continuation of BOSS, or DESI (DESI Collaboration 2016), which will collect an unprecedented number of optical sky spectra. The increase in data should improve spatial resolution, allowing DGL spectra to be computed with acceptable signal-to-noise ratios over smaller patches of sky. In addition, existing and future surveys may allow for dust maps with higher resolution and higher fidelity. These maps would provide better templates for correlating optical intensity, and they could further increase the signal-to-noise ratio and decrease the bias in our results.
Another potential avenue for improving this analysis is the use of H i maps as a tracer of Galactic dust. At the high Galactic latitudes central to the present analysis, H i closely traces the dust column density without complications induced by variations in the dust temperature (Lenz et al. 2017). The use of spectroscopic H i data from the HI4PI (HI4PI Collaboration et al. 2016) and GALFA-HI (Peek et al. 2018) surveys may improve upon our use of the IRAS 100 µm map as a dust template.
Future probes of the DGL with increasing precision and better spatial resolution could map the interstellar radiation field, dust properties, and dust luminescence across the local universe.