Data-driven Discovery of Diffuse Interstellar Bands with APOGEE Spectra

Data-driven models of stellar spectra are useful tools to study nonstellar information, such as the diffuse interstellar bands (DIBs) caused by intervening interstellar material. Using ∼55,000 spectra of ∼17,000 red clump stars from the APOGEE DR16 data set, we create second-order polynomial models of the continuum-normalized flux as a function of stellar parameters (T eff, logg , [Fe/H], [α/Fe], and age). The model and data show good agreement within uncertainties across the APOGEE wavelength range, although many regions reveal residuals that are not in the stellar rest-frame. We show that many of these residual features—having average extrema at the level of ∼3% in stellar flux on average—can be attributed to incompletely removed spectral lines from the Earth’s atmosphere and DIBs from the interstellar medium (ISM). After removing most of the remaining contamination from Earth’s sky, we identify 84 absorption features not seen in unreddened sightlights that have <50% probability of being noise artifacts—with 25 of these features having <5% probability of being noise artifacts—including all 10 previously known DIBs in the APOGEE wavelength range. Because many of these features occur in the wavelength windows that APOGEE uses to measure chemical abundances, note that characterization and removal of this nonstellar contamination establish an important step in reaching the precision required for chemical tagging experiments. Proper characterization of these features will benefit Galactic ISM science and the currently ongoing Milky Way Mapper program of Sloan Digital Sky Survey V, which relies on the APOGEE spectrograph.


INTRODUCTION
Stellar spectra capture the parameters of a star's evolutionary state and record the chemical composition of the material in which it formed.Small samples of high resolution stellar spectra have been used to describe the individual element abundance distributions of the Milky Way (MW) in the local neighbourhood (e.g.Edvardsson et al. 1993;Feltzing & Gustafsson 1998;Prochaska et al. 2000;Bensby et al. 2003).With the advent of large surveys -such as RAVE (Steinmetz et al. 2006), SEGUE (Yanny et al. 2009;Rockosi et al. 2022), APOGEE (Wilson et al. 2012;Majewski et al. 2016Majewski et al. , 2017)), Gaia-ESO (Gilmore et al. 2012), GALAH (De Silva et al. 2015;Martell et al. 2017), LAMOST (Zhao et al. 2012), and H3 (Conroy et al. 2019) -has come the ability to map abundances across the disk, bulge, and halo of our Galaxy (e.g.Bergemann et al. 2014;Rojas-Arriagada et al. 2014;Nidever et al. 2014;Hayden et al. 2015;Buder et al. 2019Buder et al. , 2022;;Wylie et al. 2021;Eilers et al. 2022).These large data ensembles have also enabled new, statistically-motivated questions to be tackled about topics such as the underlying dimensionality of individual abundance distributions and the information content of stellar spectra (e.g.Mitschang et al. 2014;Ting et al. 2015;Price-Jones & Bovy 2018;Ness et al. 2019Ness et al. , 2022;;Ting et al. 2019;Feeney et al. 2021;Weinberg et al. 2022;Griffith et al. 2022).The answers to these questions are key to understanding the origin of individual elements and the utility of those elements to reconstruct the assembly history of the MW.
Chemical tagging -the ability to distinguish co-natal stars based on chemical abundances derived from spectra -is one of the foundational ideas of stellar surveys.Understanding the conditions that create particular populations of stars informs our stellar physics models and puts constraints on models of galaxy formation and evolution.In theory, stars that are born together were formed from the same gas cloud and thus share a chemical signature in their atmospheres.In practice, the level of precision required for chemical tagging is not currently feasible (< 0.02 dex; Ness et al. 2019).
The difficulties around chemical tagging become even more severe if there are unknown or unmodeled features in a spectrum, especially if those features impact wavelength regions used for measuring chemical abundances.In the visible and infrared (IR) regimes, the largest and most obvious source of non-stellar signal comes from the Earth's atmosphere, including sky line emission from OH and telluric absorption from CH 4 , CO 2 , and H 2 O.Because of detailed measurements of the night sky's effects as well as knowing the rest-frame that spectral features are produced in, astronomers are able to account for and remove the bulk of Earth's atmosphere's signature.However, many spectra suffer from imperfect sky line and telluric removal, which leaves residual features capable of confusing spectral analysis pipelines.
Another (often ignored) source of contamination comes from intervening dust and gas along the line-ofsight (LOS) to a star.Due to the velocity offset between gas/dust clouds and stars, spectral features from the Interstellar Medium (ISM) can appear at different wavelength locations in a set of observations at different LOS in the Galaxy.This issue is complicated further when the identification or central wavelength of an ISM-based feature is unknown or poorly constrained.Without a complete and accurate model of a star's light, it is often difficult to know a priori whether a particular residual feature is caused by non-stellar sources or is simply unknown physics/missing chemical species in the model.
One common detection and characterization method for diffuse interstellar bands (DIBs) is to measure a feature's presence in multiple spectra of different stars and then to show correlations between ISM properties (e.g.extinction from dust) and that feature's strength.Typically, O-and B-type stars have been used as the preferred background "lighthouses" to study interven-Table 1.
Most precise measurements of rest-frame wavelengths for currentlyknown DIBs that fall inside of the wavelength regions covered by the APOGEE spectrograph a .DIBs that fall between the wavelength coverage of the three APOGEE detectors have been omitted.2011) 15272.42 ± 0.04 Zasowski et al. (2015) 15616.13± 0.07 Elyajouri et al. (2017) 15651.38 ± 0.07 Elyajouri et al. (2017) 15671.82± 0.03 Elyajouri et al. (2017) 15990 ing structure because their spectra exhibit weak or rotationally-broadened features that make DIB identification much easier.However, accurate spectral modelling of the more-numerous later stellar types can reveal previously-unknown DIBs that were hidden by stellar features while also probing DIB properties along many more lines of sight (e.g.Ebenbichler et al. 2022).
Efforts to detect, characterize, and map these DIBs have historically been focused on the optical regime (e.g.Herbig 1995), though a growing number of studies have been exploring the near-IR (e.g.Joblin et al. 1990;Geballe et al. 2011;Cox et al. 2014;Zasowski et al. 2015;Elyajouri et al. 2016Elyajouri et al. , 2017;;Tchernyshyov & Peek 2017;Tchernyshyov et al. 2018;Ebenbichler et al. 2022;Smoker et al. 2023).For instance, the ten currently-known DIBs that fall in the near-IR H-band (1.51 − 1.7 µm) wavelengths seen by the APOGEE spectrograph are summarized in Table 1, which is in stark contrast to the hundreds of known optical DIBs (e.g.559 DIBs in 5000 − 9000 Å as measured by Fan et al. 2019).It is particularly important to understand sources of IR features as this regime is able to peer through the dusty regions of our Galaxy's disk.
Astronomy's burgeoning "Big Data Era" has facilitated the development of novel data-driven approaches to understand stellar spectra that are less reliant on un-derlying physical models.A few successful techniques to characterize stellar light include using deep learning (Leung & Bovy 2019), polynomial models of stellar labels (e.g.The Cannon; Ness et al. 2015), and non-Gaussian Processes (e.g.Feeney et al. 2021).One significant benefit of data-motivated models is that they can describe stellar features -and correlations between featuresin spectra that are currently unknown to physics-based models.Additionally, the data models do not rely on many of the simplifying assumptions that are common in synthetic models (e.g.local thermal equilibrium, 1D radial stellar models and atmospheres).Finally, datadriven models are especially useful when the physics is not well constrained, such as in the low-density environments of the ISM that are currently impossible to recreate on Earth.As a star's light passes through intervening gas and dust on its way to our telescopes, it is imprinted with ISM signatures from many chemical species whose identities and properties are generally unknown.Detailed characterization of all the signatures in a spectrum are therefore important in disentangling the origin of various spectral features -which furthers the science goals of both abundance measurements and ISM studies.
Currently, the Milky Way Mapper (MWM) program of SDSS-V is using the APOGEE spectrograph to collect millions of stellar across all regions of the MW to understand its formation history and the physics of its stars (Kollmeier et al. 2017); any improvements in the APOGEE analysis pipeline will therefore have compounding effects on the MWM science goals.Finally, better constraints on near-IR DIBs can be studied with tomography techniques (e.g.Tchernyshyov & Peek 2017;Tchernyshyov et al. 2018) to develop a more complete picture of our Galaxy's ISM.
In this paper, we describe the APOGEE spectra and stellar parameters used in our analysis in Section 2 and then present a data-driven model of those spectra in Section 3. In Section 4, we study the spectral residuals to show that DIBs, tellurics, and sky lines are responsible for many of the relatively-large remaining features.We remove the Earth-based residuals to detect and characterize the remaining DIB features in Section 5. Finally, we summarize our results in Section 6.

DATA
This work makes extensive use of stellar spectra, abundances, and parameters of Red Clump (RC) stars in the MW as measured by the APOGEE spectrograph (Wilson et al. 2012;Majewski et al. 2016Majewski et al. , 2017) ) on the Sloan Telescope at the Apache Point Observatory as a component of the Sloan Digital Sky Survey (SDSS; York et al. 2000; Eisenstein et al. 2011;Blanton et al. 2017).The RC sample -defined by Bovy et al. (2014) using stellar parameters and simulated stellar evolution -boasts high spectral signal-to-noise ratios (SNR) as well as precise stellar parameters and abundances.These properties make the RC sample an ideal population for data-driven modelling and for studying non-stellar residuals.
The APOGEE spectra cover the H-band (∼ 15000 − 17000 Å) with high spectral resolving power (R ∼ 22500) and fine pixel spacing (∆λ ∼ 0.2 Å • pixel −1 ).The publicly-available spectral data were reduced using the APOGEE data reduction pipeline as described in Nidever et al. (2015) and are given in the rest-frame of each star.Our analysis uses the individual visit spectra instead of the combined visit spectra to account for changes in the LOS velocity -and, therefore, the location of non-stellar features -of each observation.Distributions of the number of visits per star and the median  2018), which masks use of the bitmask definitions of Holtzman et al. (2015) spectral SNR of the individual visit spectra are given in Figure 1.
The spectra were analysed by the the APOGEE Stellar Parameter and Chemical Abundances (ASPCAP García Pérez et al. 2016) pipeline.We use the ASPCAP T eff , log g, [Fe/H], and [α/Fe] measurements and uncertainties from Data Release 16 (DR16; Jönsson et al. 2020), while stellar ages come from the catalogue of Sit & Ness (2020).As exemplified in Figure 2, the stars in our sample occupy a relatively narrow range in stellar parameters and abundances.After noticing a minor secondary peak near ([Fe/H], [α/Fe]) = (−0.6,+0.2) dex in the [α/Fe] versus [Fe/H] panel, we removed stars with abundances above the red line to ensure that our modelling only focuses on a single chemical population.
For the individual visit spectra in our sample, we remove pixels that have SNR < 50 pixel −1 .We also set the maximum pixel SNR to 200 pixel −1 as recommended by ASPCAP, which suggests that the "uncertainty floor floor is at the level of 0.5%" 1 .Because of known superpersistence issues in the blue detector, we mask out spectral observations where the fiber number is ≤ 100; this selection removes approximately 6500 RC stars that do not have a single observation with a fiber number greater than 100.Finally, following the approach of Price-Jones & Bovy (2018), we remove data from pixels that have any of the the bitmask flags listed in Table 2.
These APOGEE individual visit spectra in the stellar rest-frame, along with the Sit & Ness (2020) ages and ASPCAP parameters and abundances, are used in combination to build data-driven models in Section 3, which are the basis of our analysis.

Preprocessing Spectra
First, we continuum normalize the individual visit spectra using iterative B-spline fitting.At each iteration, the B-spline is defined by 50 Å-spaced knots.For the first iteration, all flux measurements in a spectrum are used to define the initial spline.For subsequent iterations, the new spline is measured using only flux values that are within 3σ (in flux uncertainty) for fluxes below the spline or 5σ for fluxes above the spline.This masking of fluxes is done to ensure that strong absorption features do not overly impact the continuum measurement.
We iterate the spline fitting up to 100 times per spectrum, but stop iterating if the current iteration's spline is very similar to the previous one: where c i,j,k is the continuum spline value at pixel i from the previous iteration and c ′ i,j,k is the current iteration's continuum spline value at pixel i for observation/visit number k of star j.This condition is such that the summed absolute fractional change is smaller than 0.001%, which tends to occur when the subset of fluxes being masked hasn't changed from one iteration to another.This process usually only takes a handful of iterations (i.e.≤ 5), and virtually all of the spectra converge on a spline well before the 100 maximum iterations.
To capture any remaining continuum, we repeat this continuum-fitting process after the first model fit to define a continuum-adjustment B-spline.We divide each individual visit spectrum by the previously-defined continuum spline and the best-fit model to get an approximate noise spectrum that may still have some continuum trends in it.We then use the same iterative B-spline fitting process as above, but using 25 Å-spaced knots and a 3σ threshold above and below for masking.The finerspaced knots and the narrower threshold are because the residual spectrum ideally only consists of noise and any trends larger than ∼ 20 Å likely arise from an incomplete initial continuum removal.
The continuum-normalized flux in pixel i for spectral observation k of star j is then y i,j,k = f i,j,k /c i,j,k with corresponding uncertainty of σ i,j,k = σ f,i,j,k /c i,j,k , where f i,j,k and σ f,i,j,k are the raw flux and uncertainty values.The continuum-normalized fluxes are then used as the inputs for the data-driven modelling.

Modelling Flux using Stellar Labels
Following the approach of Ness et al. (2016), we define the continuum-normalized flux in each pixel to be a 2ndorder polynomial of stellar labels (see also Price-Jones & Bovy 2018).We note that a key difference between these bodies of work is the stellar labels used in the fitting: Price-Jones & Bovy (2018) uses T eff , log g, and [Fe/H], Ness et al. (2016) uses the same, but also includes [α/Fe] and mass, while our model uses T eff , log g, [Fe/H], [α/Fe], and age.Though a detailed comparison of these models is outside the scope of this work, we expect that our model will produce extremely similar results to the Ness et al. (2016) case because of the tight relationship between mass and age for RC stars (e.g.Martig et al. 2015).The vector of stellar labels in our analysis for star j is therefore given as (1) such that the vector contains all the parameters to the first and second powers, cross terms, and a constant.
Our model of the continuum-normalized fluxes is defined as where ⃗ θ i are the coefficients for the stellar label terms in Equation 1 for pixel i and describes the noise as a result of the uncertainty in a pixel's flux.This functional form implies that the data likelihood at pixel i is We then see that the likelihood at a given pixel is maximized when (6) which defines the best-fit coefficient vector at pixel i for a set of normalized fluxes, uncertainties, and stellar labels.
To propagate the uncertainties on the stellar parameters to the model coefficients, we repeatedly draw realizations of the stellar parameters for each star and remeasure the best-fit coefficients at each pixel.After 500 iterations, we take the median of the best-fit coefficients to be the coefficients of final model.The residual flux for pixel i of observation k of star j is defined to be for realization l of the the stellar parameters samples, ⃗ x j,l , which yields a best-fit coefficient measurement of θi,l .To propagate the uncertainty on the stellar labelsand therefore the uncertainty on the best-fit coefficients -to the residuals, we also repeatedly measure the residual flux values for the 500 realizations, giving samples of r i,j,k,l .The final residual flux, ri,j,k , is taken to be the median of these realizations, with an uncertainty that is given by where var (r i,j,k,1 , . . ., r i,j,k,500 ) is the variance of the 500 residual measurements at a given pixel for a given spectrum.In words, the residual uncertainty is a quadrature sum of the normalized flux uncertainty and the propagation of the best-fit coefficient uncertainty.
A comparison of the the best fit model and data for a single observation of one star in our sample is shown in Figure 3.The uncertainty-scaled residual distribution for this star is close to the expected unit Gaussian distribution (mean of 0, standard deviation of 1), which shows that the model does a good job of capturing the information in the spectrum.If we assume that  the fluxes and uncertainties reported by APOGEE accurately describe the data, then the best-fit Gaussian to the residual distribution (orange histogram) having a width of ∼ 1.15 ± 0.1σ implies that there is information in this star's spectrum that the model does not capture at the level of ∼ 15% of the flux uncertainties.
If we instead look at the uncertainty-scaled residual distribution across all observations and all pixels, we get the distribution shown in Figure 4.As before, the distribution is quite similar to the unit Gaussian, which sug- By highlighting a single pixel, we focus on the residual information from all of the spectra at a particular stellar rest-frame wavelength.Left: A well-modelled pixel, where the uncertainty-scaled flux residuals (blue in top panel) follow the expected unit normal (orange in top panel).In the middle panel, the residuals are binned by heliocentric velocity (number of spectra in each bin given by N * in the legend), and we see no differences between the velocity bins.The bottom two panels show the median and standard deviations of the velocity-binned histograms, again showing no noticeable trend in the residuals with velocity.Right: A pixel where the residuals do not follow a unit Gaussian, or any Gaussian for that matter.After binning by heliocentric velocity in the middle panel, and plotting the binned medians and standard deviations in the bottom two panels, we see a noticeable trend in the median and width of the residual distribution with velocity.
gests that the model performs well across all the spectral observations in our sample.Again, that the width is greater than 1σ reveals that there may be more information in the spectra than the model is able to describe.
Using the residual histogram and best-fit Gaussian distribution in Figure 4, we vertically align the distributions so that they have the same height at their peaks so that we can quantify their difference in the wings (i.e.excess in the data over the expectation).First, we find that the data excess in the wings corresponds to ∼ 3.3% of the total pixels.Next, we draw many realizations of residual flux measurements from each bin in the blue data histogram where the data counts exceed the expected distribution.The results of this process are summarized in Figure 5, where the blue curve shows the median distribution of the size of flux residual in excess of the expectation.The blue point above the histogram shows the median and 68% region of the distribution, and the grey histograms show individual realizations of residual flux samples.In summary, the residuals that are not explained by the model have a median size of ∼ 2.9% of the stellar flux, with a 68% region covering 2% to 5% of the stellar flux.

STRUCTURE IN RESIDUALS
To further explore the residuals, we begin by looking for trends in the flux residuals as a function of wavelength (in stellar rest-frame) and other explanatory variables.In many APOGEE pixels, we find that the model almost perfectly describes the data; for instance, the pixel summarized in the left panels of Figure 6 shows a residual distribution (scaled by the residual uncertain- with each panel consisting of 21 pixels; these cutouts include 15 elemental features, one region we've identified as being mostly continuum, one region of high telluric contamination, and one region near the strongest known DIB in APOGEE (15272 Å feature).The residual spectra have been Gaussian-smoothed in the horizontal direction with a kernel width of 5 km/s.Where applicable, ionization levels are shown.For the C window, the dominant feature is from CO, and in the O window, the dominant feature is from OH.Some elemental regions (e.g. the DIB region, 3rd from the top) show a residual minimum that moves through the element window as a function of heliocentric velocity; this can be explained as an unmodelled residual feature that is not in the same reference frame as the star.
ties) that agrees very well with the unit Gaussian.When the stars are split up into different heliocentric velocity (V HELIO ) bins -as is shown in the left middle panel -we see little difference between the residual distributions, and no obvious trend in those distributions' medians or standard deviations (left bottom panels).
On the other hand, there are some pixels where the residual distributions are distinctly non-Gaussian and have widths that are much larger than 1σ.An example of this is shown in the right panels of Figure 6 for a pixel near the central wavelength of the strongest known DIB in the APOGEE wavelength range (e.g.λ 0 = 15272.42Å, from Zasowski et al. 2015).Breaking the stars up into V HELIO bins in the middle right panel, we now see significant differences between the residual distributions' shapes as well as their medians and widths (right bottom panels.In general, the extreme velocity bins have more positive residual medians and smaller standard deviations while the moderate velocity bins have negative residual medians and larger widths. We next look at the trends in the residuals with heliocentric velocity across neighbouring pixels at different wavelength cutouts.First, we sort the stars by their heliocentric velocity and then smooth the residuals (using a combination of inverse variance weighting and Gaussian smoothing based on V HELIO at each pixel for each residual spectrum) to produce Figure 7.Each of the 18 panels show the smoothed residuals using 21 APOGEE pixels (i.e.width of ∼ 4 Å) centered on a particular feature, with the central wavelength of that feature listed on the right edge of the cutout.These regions correspond to 15 elemental features used by ASPCAP to measure abundances (these are the same central wavelengths as used by Feeney et al. 2021), one region we identify as having minimal visible features (i.e.continuum), one region where we notice a large amount of residual Earth-based contamination, and one region around the strongest DIB in the APOGEE wavelength range.
The y-axis of each panel in Figure 7 corresponds to wavelength relative to the central wavelength label on the right of each plot, with blue wavelengths at the top and red wavelengths at the bottom.The x-axis is the sorted list of heliocentric velocities of the stars such that a vertical column in this figure corresponds to a single residual spectrum of a particular observation with that velocity.
By eye, some of the wavelength cutouts (e.g.Fe I, Ni I, C, Cr I, continuum) show small amplitude in smoothed residuals and not much correlation with heliocentric velocity.Other cutouts (e.g.K I, DIB, P I, sky, Na I, Yb II) show relatively strong trends in residuals with V HELIO .
Focusing on the DIB panel in particular, we see a residual minimum (i.e. an absorption feature) move through the cutout as a function of heliocentric velocity.This is because the strong DIB feature, while present in many of the APOGEE spectra, appears at different wavelength locations in the stellar rest-frame spectrum as a result of the offset between the velocity of the star and the velocity of the DIB-producing source.The wavelength of a DIB feature in a star's rest-frame spectrum is given by: where λ rest,DIB is the wavelength of the DIB in its own rest-frame, v * is the star's heliocentric velocity, and ∆v is the LOS velocity offset between the DIB source and the star.Because DIB features will show up at different wavelength locations in the stellar rest-frame spectra, our model is not able to describe the feature and thus leaves DIB residuals behind.Similarly, features from the Earth's sky are also likely to show up as residuals because they too appear at different wavelengths in the stellar rest-frame spectra: For the Earth-based features, sorting by heliocentric velocity means that the features will occur at similar wavelength locations, which is why they appear strengthened in the velocity-sorted, smoothed residual  panels.Comparing the features in the sky and the DIB panels, we also notice a difference in the horizontal width of the features; DIB features are known to be quite broad (e.g.FWHMs of 3.5 Å to 6 Å for a handful of 15272 Å DIB detections in Table 1 of Zasowski et al. 2015), so it agrees with expectations that the DIB panel residual has a visually larger width than, for example, the width of the local maximum diagonal stripe in the left half of the sky panel.Based on this difference, the narrow width of features in the P I and Ce II panels suggest these residuals are likely caused by the Earth instead of DIBs, while the large width of the residual minima stripes in the K I and Na I panels may originate from DIBs.
The particular pattern seen in the DIB panel of Figure 7 is a result of the average RC heliocentric velocity having a correlation with the average ISM heliocentric velocity in the MW disk.Because the RC stars in our sample are generally quite old (average of ∼ 2.8 Gyr, from Figure 2), they have had a longer time to kinematically decouple from the gaseous disk they were likely born in.In the Galactocentric cylindrical radius range of our stellar sample -95% of the APOGEE RC sample is within 6.4 < R < 13.4 kpc with a median at 10.1 kpc -the gaseous disk has been seen to have relatively flat rotation at ∼ 220 km/s (e.g., from the HII measurements of Brand & Blitz 1993).Using RC stars from APOGEE and GALAH, Khanna et al. (2019) measure rotation velocities that are also relatively constant (210 < V ϕ < 230 km/s) in the 6 < R < 11 kpc region for |Z| < 0.75 kpc.In particular, their R > 9 kpc measurements -where approximately 75% of our APOGEE RC stars fall -shows remarkably stable V ϕ ∼ 215 km/s.This suggests that RC stars can be thought of as belonging to a rigidly rotating disk that is spinning ∼ 5 km/s slower than the gaseous disk.
When we sort the DIB-based features by heliocentric velocity, we are largely sorting by azimuthal angle, as can be seen in Figure 8.If we assume both the RC stars and the ISM are described by rigid-rotating disks, then there is an average stellar V HELIO and an average DIB V HELIO for each bin in azimuthal angle; it is the relationship between the average stellar V HELIO and average DIB V HELIO that produces the pattern in the DIB panel.We can explore this relationship explicitly by tracing the location of the DIB minimum as a function of stellar rest-frame wavelength.This is done in Figure 9, where the orange best fit line to the DIB minimum location and the known rest wavelength of the DIB feature allow us to estimate the DIB source velocity as a function of stellar V HELIO .
We next compare the relationship we measure for gas V HELIO as a function of stellar V HELIO to our expectations.To that end, we assume the RC stars all belong to a rigidly rotating disk with V ϕ ∼ 215 km/s.Similarly, we assume there is also a gaseous disk rotating at V ϕ ∼ 220 km/s.For each RC star in Figure 8, we use the star's Galactic position to obtain the expected velocity vector from the rotating stellar disk, which is then transformed into a heliocentric velocity after accounting for the sun's position and motion2 .We follow a similar procedure for generating the simulated V HELIO of the gas, but we now need to assume a distance from the sun to the DIB source along each LOS.To be agnositic of the exact DIB source distance, we choose a  fraction of the total LOS distance to each star, and use those new implied Galactic coordinates and the gas disk to measure V HELIO for the intervening gas.We repeat this process for distance fractions of 100%, 50%, 10%, and 1% of the total distance to each star.A comparison of the velocity difference between simulated gas and simulated stars is shown in Figure 10, where the data points are colored by the different fractional distances of the gas.The orange line is the result of using the best fit line in Figure 9 with the known rest-frame wavelength of the DIB to measure heliocentric velocity of the DIB as a function of stellar V HELIO ; to be clear, the orange line in Figure 10 is not a fit to the simulated data.We see good agreement between the orange line and the simulated cases when the DIB sources are, on average, between 1% and 50% of the distance to the stars.This implies that the velocity offset function we measure from the location of the DIB minimum in the stellar rest-frame is caused by a source that is in the foreground of the stars, as we would expect for an ISM-based absorption feature.
Additionally, we see that the gas velocity offset in the simulated data has a relatively tight correlation with stellar velocity.This is ultimately what causes the pattern we see in the DIB panel of the smoothed residual in Figure 7; for the RC sample, stellar velocity has an approximately linear relationship with the gas offset velocity, so the wavelength location of a DIB in the stellar rest-frame, Equation 9, effectively becomes a function of stellar velocity alone.This explains why sorting by V HELIO amplifies DIB signals and causes their location to move smoothly in wavelength.
By comparing the apparent strength of residual spectral features with various parameters, we found a slight correlation in residual strength with spectral SNR.In particular, residual spectra with lower median SNR show higher levels of contamination, especially in wavelengths regions where tellurics and sky lines are known to occur.We shift the residual spectra to the observer reference frame using each observation's LOS velocity to allow for Earth-based features to align in wavelength.We then inverse-variance combine residual spectra in different spectral SNR bins, which yields the coadded spectra shown in Figure 11; the combined uncertainty in each of the SNR bins are approximately the same, with median combined SNRs of ∼ 550 pixel −1 .The top panel In even the highest SNR bins, the residual features are commonly on the order of 1% − 2% of the stellar flux, and as high as 4% − 5% in the extreme cases, which is similar to the peak absorption depths seen in 15272 Å DIB detections (e.g. the bottom panel of Figure 3 in Zasowski et al. 2015).The structure we observe in our residual spectra reveals that a complicated combination of Earth-based features and DIBs are common across the APOGEE wavelength range.These trends of residuals with spectral SNR and heliocentric velocity highlighted in this work are what prompted the authors of Ness et al. (2022) to include SNR and V HELIO as explanatory variables in their regression fitting to measure the relationship between abundances and stellar evolutionary state.
We posit that the non-stellar residuals seen in Figure 7 may be responsible for the larger-than-expected scatter in the some of abundance distributions measured by AS-PCAP (e.g.Zasowski et al. 2019;Jönsson et al. 2020); some elemental abundances (e.g.Na) are determined by flux measurements at a small number of pixels, so the presence of DIBs at those wavelengths may have a relatively large impact.If we attribute a large portion of the residual excess in Figure 5 to non-stellar features, then it wouldn't be unexpected to see a large effect from features with local minima that are ∼ 3% of the stellar flux.

REMOVING SKY LINES AND TELLURICS TO IDENTIFY DIFFUSE INTERSTELLAR BANDS
To identify the DIBs present in the residual spectra, we first need to characterize and remove the more common contamination from Earth-based features.One advantage of isolating the Earth-based residuals versus the DIBs is that we are easily able to shift to the rest-frame of the sky line and telluric residuals.
Because almost all the spectra have residual Earthbased features, but not all spectra have DIB contamination, it is useful to split the residual spectra up into two groups: those with and without the strongest DIB feature.Once we have these two groups, we can study the Earth-based contamination using the low-DIB-strength group.Those results can then be used to remove the Earth-based contamination from the high-DIB-strength group.
To accomplish this, we need to rank the spectra by their DIB strength.We fit an inverted Gaussian to the strong DIB feature that is present around 15272 Å DIB in each residual spectrum, after shifting to the heliocentric frame using the V HELIO for each observation; an example of this is show in Figure 12 3 .This allows us to measure a LOS velocity and equivalent width for the strong DIB in each residual spectrum.We then use the EW measurements to define a low-and a high-DIBstrength group; residual spectra that failed to fit an absorption feature near 15272 Å DIB were automatically assigned to the low-DIB-strength group, and then a threshold of EW = 1 × 10 −2 Å was chosen to divide the spectra with DIB detections to ensure an approximately equal number of spectra (∼ 21000) in the lowand high-DIB-strength groups.For the residual spectra in the high-DIB-strength group, the 15272 Å DIB EW is compared to the K-band reddening4 in Figure 13 with data points colored by the EW bins used in our later analysis.
For the low-DIB-strength group, we fit a 2nd-order polynomial model to the observer-frame residual fluxes at each pixel using a process similar to what is described in Section 3.2.In this case, however, the stellar label vector in Equation 1 is replaced by and the data we are fitting are the residuals in the observer frame.This model enables us to quantify the predictive power of SNR alone to describe the spectral residuals that are not captured by our five label model from Section 3.2.
Examples of the resulting best-fit model for different SNR bins are shown in Figure 14.Overall, this model captures much of the behaviour we see in the coadded residuals of Figure 11.There are a few regions that stand out visually (e.g. in the blue shaded regions where CO 2 telluric features are particularly strong) that have a saw-toothed shape to the residuals; this is a shape that can be created by subtracting two Gaussians with the same width but a slight offset in mean.We suggest that these Earth-based residuals are caused by a wavelength mismatch between the raw spectra and the telluric models being subtracted.To be clear, these stacked residual spectra are produced by stacking a couple thousand spectra together, which means that any issues of random error in the individual wavelength solutions should be minimal.Instead, these stacked spectra are revealing systematic offsets in the wavelength solution from the telluric model.
As shown in Nidever et al. (2015), the APOGEE data reduction pipeline's telluric and sky removal performs very well in most cases.Holtzman et al. (2018) further shows improvement in APOGEE's Earth-feature removal; however, they also point out that the reduction pipeline's approach is to flag and ignore pixels near particularly strong telluric and sky lines instead of a more computationally expensive approach that might achieve smaller residuals.Of particular interest, Figure 3 of Holtzman et al. (2018) reveals the same saw-tooth residuals that we see in our Figure 14, with their figure focusing on the 15680 − 15800 Å region that is sensitive to tellurics from CO 2 .While the APOGEE reduction process works well for their science goals, our analysis has shown that detailed accounting of Earth's atmosphere's light can reveal previously-hidden information.
By looking at the wavelengths and shapes of the Earth-based features in Figure 14, we can see that some of the smoothed residual patterns in Figure 7 are likely explained by the night sky.The Ce II panel of Figure 7, for example, shows a repeating alternation of residual fluxes from low to high, which is visually dissimilar to the DIB panel which shows a single feature moving through the cutout window.The Ce II panel's cen-tral wavelength of 15789.1 Å places it in a region of saw-toothed features in Figure 14, so this mottled pattern that appears after sorting by heliocentric velocity is likely the result of these saw-tooth shapes constructively and destructively combining.
While this SNR-dependent Earth-residual model performs well in describing the observer-frame residuals, we choose to not use this model in our subsequent DIB analysis without extensive further testing.Instead, we remove sky line and telluric features by pairing up residual spectra in the low-and high-DIB-strength groups based on spectral SNR.That is, for each high-DIB-strength residual spectrum, we identify a low-DIB-strength residual spectrum with median spectral SNR that is within 5 pixel −1 of the high-DIB-strength spectrum's SNR.We ensure that each low-DIB-strength residual spectrum is used only once.
Once each high-DIB-strength residual spectrum has a corresponding low-DIB-strength residual spectrum, we subtract each pair and add their uncertainties in quadrature in the observer frame.This has the effect of removing the Earth-based residuals from the high-DIB-strength residual spectra while leaving the possible DIB features relatively untouched.The subtracted spectra are then shifted from the observer frame to the 15272 Å DIB rest-frame using the velocity measured in the inverted-Gaussian fitting for the high-DIB-strength spectrum in the pair.Because each pair of spectra are likely not taken at the same exact time or with the same observing conditions, their subtraction may not perfectly remove Earth-based features in each resulting spectrum.However, averaging these subtracted spectra together does remove the vast majority of obvious telluric features.Because our subsequent analysis is based on combinations of spectra, we believe this approach is a simple and model-free method of dealing with telluric contamination.Furthermore, we find that our final results are insensitive to the particular pairings of lowand high-DIB-strength spectra, so long as the two spectra are similar enough in SNR.
The next step is to break the spectra up into different bins based on the strength of the 15272 Å DIB.We define 10 DIB-strength bins such that approximately 2100 spectra are in each bin; these are shown by the different colors in Figure 13.We then combine the residual spectra in these bins using population fitting to measure a population mean and uncertainty/standard deviation at each pixel.The precise statistics used for this combination are explained in Appendix A, but the key takeaway is that the population fitting takes into account both the uncertainty in each residual measurement as well as the spread in those measurements to give back realistic means and uncertainties at each pixel.
The results of this population-fitting process are shown in Figure 15, where the colors denote the same DIB-strength bins as in Figure 13.The vertical red lines show the locations of the previously-known DIBs in Table 1; by eye, the increasing DIB-strength bins show increasing depth of features at most of these wavelengths.For the 15225 Å feature, the large literature uncertainty in central wavelength may suggest that the true central wavelength of that DIB should correspond to one of the features near 15250 Å instead.We also notice additional absorption features (e.g.slightly redward of 15700 Å) that have shapes similar to the previously-known DIBs.
We recognize that all the DIB features are likely not in the same rest-frame as the 15272 Å DIB, as they might be produced by different intervening clouds along the sight-line to each star and therefore have different relative LOS velocities.However, the velocity difference between DIB sources along a single LOS is not likely to be large enough to shift the wavelength by more than a few APOGEE pixels.For instance, the vast majority of the 15272 Å DIB velocities we measure are in the ±30 km/s range, so at ∼ 16000 Å, a 10 km/s velocity offset would manifest as a ∼ 0.5 Å offset in wavelength.When we combine spectra in each DIB-strength bin, we may slightly reduce the signal of DIB features from other sources because of a combination of wavelength offsets and their strength not correlating with the 15272 Å strength, but we still expect to detect their presence.Our decision to use the 15272 Å DIB restframe will, however, preserve the strength of features that are being produced by the 15272 Å DIB source, so this technique is particularly useful for identifying DIBs that correlate with the strongest DIB in the APOGEE wavelength range.For future work where extremely precise central wavelengths are needed or the goal is to find every possible DIB in the APOGEE wavelength range, one approach would be to cross-correlate the residual spectra using different wavelength cutouts to find the velocities that amplify the signals of all DIBs.
To identify possible locations of DIBs, we smooth the highest DIB-strength residual spectrum (i.e. the darkest purple line in Figure 15) with a Gaussian kernel width of 5 pixels; after some visual vetting, this yields 133 local maxima and minima as possible DIBs in emission and absorption.Of the possible DIBs, we highlight the locations of 35 features in Figure 15 that we will ultimately find are likely produced by the same source that produces the 15272 Å DIB; we choose to not show all 133 possible features for visual clarity.As expected, this step finds features nearby to all of the previously-known DIBs.Many of the local extrema we find are likely spurious detections, but we choose to err on the side of testing too many extrema versus applying a more complicated thresholding at the detection step.For each local extremum, we then step redward and blueward from that wavelength in the smoothed spectrum until the slope is near 0 to define a useful wavelength region around each feature.We then manually check these regions to confirm that they cover the entirety of an apparent feature, tweaking the boundaries where necessary.This defines the wavelength cutouts we will later use for measuring the EW strength of the possible DIBs and for quantifying DIB detection probabilities.The list of possible features agree well with the features we find when we smooth a population-combined spectrum using all ∼ 21000 residual spectral pairs; we choose to use the feature locations from the highest DIB-strength spectrum instead of the total population because the signal from true DIBs will be strongest in the high DIB-strength bin.
As mentioned in the last few paragraphs, the central wavelengths we measure in the 15272 Å DIB rest-frame for DIBs that are produced by a different source may be offset from their true value by a few pixels, but this is still an improvement on the wavelength location for a few of the previously-known DIBs in Table 1.Additionally, because the spectra in each DIB-strength bin are a combination of ∼ 2100 individual residual spectra, a systematic velocity offset between different DIB sources would be required to produce a significant shift in the central wavelength location.This implies that our central wavelength measurements are likely within a handful of pixels of the true rest-frame wavelength of each DIB feature.To be conservative, we estimate the uncertainty on the central wavelengths of the possible DIBs to be ∼ 1 Å.
We next measure the equivalent width of the local extrema in each of the DIB-strength bins of Figure 15.Using the wavelength cutouts defined above, we empirically integrate the flux.Specifically, we use the 5 pixels at the edges of each cutout (3 inside, 2 outside of the cutout) to define an average flux and wavelength on the blue and red edge; these blue and red average fluxes define a line for each DIB-strength bin that we use at the local continuum measurement.After subtracting off the local continuum line for each bin, we perform a trapezoid numerical integration to measure the area of the feature.To propagate uncertainties in the flux, we repeated sample flux measurements -including the fluxes used to measure the local continuum line -to get a distribution on the empirical EW for each bin.The EW measurements we report are the median and standard  deviation of those realizations.An example of the resulting EW measurements is shown in the left panel of Figure 16 for the possible DIB near 15706 Å, with a cutout of the feature shown in the right panel for the different DIB-strength bins.We note that the nearest Earth-based feature is a sky emission line at 15720.6 Å, which is sufficiently far away from the new 15706 Å DIB to have a negligible impact on our results.
To characterize how much the EW measurements change between the bins, we perform a population fit to the EW measurements of the possible DIB feature (e.g. the y-axis values of the left panel); this population fit again follows the approach detailed in Appendix A, and an example of the resulting medians on the population mean and width are shown with a black line and grey region in the left panel of Figure 16.We are particularly interested in the population width because a large width indicates that the EW measurements between the bins are quite different.
Of course, each possible DIB feature may only be a spurious detection; that is, a feature may only show up in the highest DIB-strength bin by random chance alone.To determine the significance of a detection, we repeat the binning, combining, and EW-measuring process, but this time, we bin randomly.To be clear, instead of using the strong-DIB-strength ranking to bin our spectra, we assign spectra randomly into 10 different bins, with the same ∼ 2100 spectra per bin.Using random binning, we expect that the features in each bin will be quite similar to each other, as will their EW measurements.
To quantify exactly how strong a detection is (i.e. its significance above randomness), we also measure the population width of the EW measurements from random binning.In particular, we are interested in how the population width distributions compare between the random binning and the binning by the strong-DIBstrength cases.To account for the variation in the random binning that might occur by chance alone, we repeat the random binning 5 times.The resulting population width distributions of the populations agree quite well between the realizations for each possible feature.This implies that our results aren't overly sensitive to the particular choice in random binning.To be careful, however, we average the population width distributions from the 5 random realizations together.Measurements from this point on that refer to "random" come from the "averaged random" results.
Recognizing that strong-DIB-strength might not be the only/best parameter to explain the change in a feature's strength between bins, we explore a few additional binning options: [Fe/H], median spectral SNR, and A K reddening.Respectively, these can be thought of as testing for features that are the result of residual chemical information that the model did not capture, remaining Earth-based contamination or other SNR effects, and DIBs that correlate with the amount of LOS dust but Comparison of the population width distributions of EW measurements around the possible DIB absorption feature shown in Figure 16 when sorted and binned by different parameters.In each case, we measure the EW of the possible DIB feature from the combined residual spectra in different parameter bins, and then measure the population width of the EW measurements (e.g. the shaded grey region of Figure 16).Left: Distributions on the EW population widths for different parameter bins; the black histogram is the resulting average of 5 realizations of random sorting and binning.Right: Distribution on the difference between the population widths in the left panel and the random population width.In this case, the largest population width is achieved after binning by DIB-strength, suggesting that the absorption feature in Figure 16 is best correlated with 15272 Å DIB-strength.
don't originate from the 15272 Å DIB source.As before, we sort by a given parameter and then bin the spectra so that there are approximately 2100 spectra per bin.In every case, we measure EWs in each bin at each possible feature, and then use those EWs to measure a population width distribution, which are then compared to the results from random binning.An example of these population width distributions for the different sorting parameters is shown in the left panel of Figure 17; in the right panel, samples from the random distribution have been subtracted from samples in the other distributions to get a distribution of σ SORT − σ RANDOM .
We can then directly integrate these population width difference distributions to measure the probability that the population width for a given sorting parameter is greater than the population width from random sorting.In cases where none of the probabilities are greater than 50%, we decide that randomness alone likely produced the observed feature.We also use the medians of these distributions to determine which sorting parameter yields the greatest difference from random, implying that a particular feature is best explained by that parameter.Because the relationship between each feature's EW measurements and a given sorting parameter is potentially complicated (i.e.not a simple functional form), we argue that our approach is a cautious, statistically-motivated way of making detections without assuming a parameterized relationship.We are simply answering the question: "Which of the sorting parameters produces the largest difference in EW measurements between the bins?"If a feature has its best sorting from DIB-strength, it may suggest that feature is produced by the same source as the 15272 Å DIB.If the best sorting is A K , it may be that feature is still truly a DIB, but that it has a different source than the 15272 Å DIB.
We summarize the probabilities from all the sorting parameters for all 133 possible features in Figure 18.At our best estimate of the wavelength for each feature, a vertical line of 4 points are plotted to show the probability above random for each of the different sorting parameters; the point with the highest median population width is enlarged compared to the others, provided that the maximum probability above random is greater than 50% (i.e.there is only one enlarged point per feature).These probabilities reveals that 0 of the features are best explained by random chance alone (i.e.spurious detections), while the other 55, 25, 29, and 24 features correlate best with the 15272 Å DIB strength, [Fe/H], A K , and spectral SNR respectively.Comparing to the previously-measured DIBs (vertical red lines), we recover > 93% detection probabilities in each case.Apart Figure 18.Summary of the probabilities that a particular sorting parameter's population width of EW measurements is greater than random chance can explain (e.g. the legend values in the right panel of Figure 17) for all possible DIB features.The possible DIB wavelengths are our best estimates of each feature's central wavelength; the exact locations of the features that are best described by AK or DIB-strength sorting are given provided in Table 4.There are 0 possible DIB features where the maximum probability of all sorting parameters is less than 50%, meaning that none of the features are best explained by chance alone.At each possible DIB wavelength, the four sorting parameter's probabilities are plotted in a vertical line; the sorting parameter with the maximum median population width is enlarged compared to the others.Of the 133 possible DIB features, 55 show the largest difference from random when binned by 15272 Å DIB strength, another 29 features are potentially DIBs produced by sources other than the 15272 Å DIB, 24 features may originate from SNR effects, and 25 features are best explained by [Fe/H] binning which may be a result of the model imperfectly describing the data.
from the 15225 Å DIB -which finds best sorting using A K -all remaining previously-measured DIBs in Table 1 have best sorting from DIB-strength.These results may suggest that most of the previously-known DIBs in APOGEE are produced by the same source/species along each LOS while the remaining DIB has a different origin.
While we may not have perfectly captured the central wavelength or cutout region around each of the possible features, we emphasize that our process is likely biased towards calling a true DIB spurious than the other way around.A more rigorously-defined cutout region and fitting of the local continuum would likely yield stronger detection signals than we report, suggesting there may be even more DIB features than we discover.Highlydetailed analyses in the future may reveal DIBs that are not included in our results, and possible DIBs we classify as having non-ISM origin may be found to have significant signal.
In general, the average probability above random for detected features is higher for the DIB-strength sorting compared to the other sorting parameters (e.g.mean Pr DIB = 0.89 versus mean Pr AK = 0.83).This is likely a consequence of our choice to bin in the 15272 Å DIB restframe; as discussed previously, DIBs from other sources with slightly different rest-frames will experience a reduced signal, so their difference from random is not as strong.Table 3 lists the number of features in each bin with detection probabilities (i.e.probability above random) above various thresholds.This is equivalent to asking how many of the enlarged points persist above changing heights of the shaded grey region in Figure 18.We notice that the number of detected features in the A K sorting drops off faster than the number of features in the DIB-strength bin, and this again is likely caused by binning in the 15272 Å DIB rest-frame.
Linking the resulting DIBs back to the features seen in Figure 7, many of the element windows occur near newly-detected DIBs.The Na I panel, for instance, is near a detected feature at 16382.1 Å that is best explained by A K sorting with 93% probability above random.This feature likely explains the diagonal streak in the smoothed residuals of the Na I panel.It may also be responsible for increasing the scatter in Na abundances that ASPCAP reports, emphasizing the need to account for DIBs in future stellar abundance pipelines.
To summarize the patterns in the V HELIO -sorted, smoothed residuals, the Fe I, Ni I, C, and Cr I panels show minimal residual features suggesting these windows are free of significant non-stellar light.Like the As an initial step towards identifying the chemical species that producing each DIB, we compare the central wavelength locations to known hydrogen transitions in the APOGEE wavelength range.The newlydetected DIB in Figure 16, for instance, is within 1 Å of the Brackett series n = 15 to n = 4 transition (λ 0 = 15705.0Å), suggesting that atomic hydrogen in the ISM is a probable source of this feature.That we are able to detect features at many known hydrogen wavelengths without a priori searching there bodes well for our general methods.
A complete summary of the 84/133 DIBs whose best sorting parameter is either DIB-strength or A K is listed in Table 4. Eight of the features are in emission and the remaining seventy-six are in absorption.All wavelengths are given in the 15272 Å DIB rest-frame, and the Wavelength Range column gives the window we used to measure the EWs; the continuum wavelengths are defined by the 5 pixels nearest to the values in the Wavelength Range column.By summing the probabilities above random, we expect that 73.3 of the 84 DIB detections are truly DIBs.For features that have central wavelengths within 2 Å of a known hydrogen recombination line, we give the name of the series and level.Many of these newly-discovered DIBs occur in the same wavelength regions that were obscured by incomplete sky line and telluric removal; this, combined with their sub-percent sizes, explains why our analysis has been able to reveal these features for the first time.

SUMMARY
We have created data-driven models of RC stellar spectra using ∼ 5.5 × 10 5 individual observations of ∼ 1.7 × 10 5 stars from the APOGEE dataset.The modeling uses five parameters -T eff , log g, [Fe/H], [α/Fe], and age -to predict the spectra, and the resulting models agree quite well with the data (residual Gaussian mean near 0 and width of ∼ 1.16σ) across the APOGEE wavelength range.This implies that these five labels are sufficient to explain the majority of information present in the spectra.Consequently, there is not a substantial amount of residual information that may be leveraged for pursuits such as chemical tagging.Though the residuals of the data-minus-model are relatively small (∼ 3% of stellar flux on average), it is very important to understand and isolate their astrophysical origin.We discover that there are many pixels where the residuals suggest that a significant number of non-stellar features are also present in the stellar spectra.We identify which of these features are likely Earth-based and which are likely Diffuse Interstellar Bands.We find 84 possible DIBs in APOGEE spectra that have less than 50% probability of appearing by chance alone, including all previouslydiscovered DIBs in this wavelength region.
Our key results include: 1.The residuals of our data compared to the model show correlations with heliocentric velocity, which we show is evidence that many of these residual features are not in the stellar rest-frame, such as sky lines, tellurics, and DIBs.These residual features appear at the level of 3% of the stellar flux on average (Section 3.2, Figure 5; Section 4, Figures 6, 7, and 11); 2. The size of Earth-based residuals appear anticorrelated with spectral SNR.The shape of these residuals suggest that they may be removed by correcting for a wavelength offset between the sky model and raw observations (Section 5, Figure 14); 3.After removing the Earth-based residuals, we combine residual spectra in the rest-frame of the strongest DIB in the APOGEE wavelength range (λ 0 = 15272 Å).We detect 84 DIB features in absorption (including all of previously-measured DIBs) and emission that show highest correlation in strength with either K-band reddening (A K ) or EW strength of the 15272 Å DIB feature (Section 5, Figures 15 and 18, Tables 3 and 4).
Future work based on our results will focus on measuring the impact that unaccounted-for DIB and Earthbased features have on ASPCAP-measured abundances.It would also be worthwhile to compare the DIB wavelengths to lines produced by chemical species (i.e. more than just atomic hydrogen) that are known components of the ISM, as well as to correlate the new DIB strengths with additional tracers of ISM density.Another useful undertaking that would advance this work is detailed joint modelling of the DIBs, sky lines, tellurics, and stellar spectra in APOGEE, such as employed by the MADGICS method of Saydjari et al. (2023) using Gaia spectra.A complete understanding of all sources of spectral features is necessary for chemical tagging experiments, and this will have the added benefit of improving our understanding of the ISM.
The authors thank the anonymous referee for comments that helped improve the clarity of this paper.KM, CMR, and PG were supported by NSF Grant AST-2206328.We thank Julianne Dalcanton and the Center for Computational Astrophysics Astrophysical Data Group for useful discussions.KM thanks his PhD advisors, CMR and PG, for enthusiastically supporting his decision to pursue a Designated Emphasis in Statistics, without which, this work would not have been possible.Nearest detected DIB to the previously-known DIBs in Table 1.
Gaussian for p (r i |σ r,i )).The full posterior distribution is then: p(r i ,σ 2 r,i , r ′ i,1 , . . ., r ′ i,n * |r i,1 , σ 2 i,1 , . . ., r i,n * , σ 2 i,n * ) ∝ p (r i , σ r,i ) • n * j N r ′ i,j |r i , σ2 r,i • N r i,j |r ′ i,j , σ 2 r,i,j (A1) We next see that the posterior full conditional on r ′ i,j is given by: p(r ′ i,j |r i , . . . ) ∝ N r ′ i,j |r i , σ2 r,i • N r ′ i,j |r i,j , σ 2 r,i,j = N r ′ i,j |µ r ′ ,i,j , σ 2 r ′ ,i,j where σ 2 r ′ ,i,j = σ−2 r,i + σ −2 r,i,j −1 and µ r ′ ,i,j = σ 2 r ′ ,i,j • σ−2 r,i • ri + σ −2 r,i,j • r i,j .We can then use these results to integrate over r ′ i,j in the full posterior distribution of Equation A1 to find the marginal posterior of (r i |σ r,i , data): With the functional forms of the distributions in Equation A3 and A4 in hand, we are able to draw samples of (r i , σr,i |data) fairly quickly.First, we evaluate p(σ r,i |data) for a reasonable range and number of σr,i values, then we use those probabilities to draw σr,i samples.Next, we use those σr,i samples to draw samples from the (r i |σ r,i , data) Gaussian distribution, which is an relatively easy and efficient step.Once we've repeated this process a sufficient number of times, we can take the median of the (r i , σr,i |data) samples as our best estimate of the underlying population mean and width.Sometimes, like in the main text of this work, we are also interested in the distribution on (r i , σr,i |data) itself, and this technique allows us to study this distribution.

a
TheCox et al. (2014) values have been converted from their reported Air wavelengths to Vacuum.

Figure 1 .
Figure 1.Top: Histograms of the number of spectroscopic visits per star.The bin width has been set to 1.The median and mean number of visits per star are shown with vertical red lines.There are a total of 55028 individual visit spectra for the 17104 RC stars in our sample.Bottom: Histogram of the median spectral SNR for each visit used in our analysis, with mean and median shown with vertical red lines.

Figure 2 .
Figure 2. Distribution of the stellar parameters for the 17104 RC stars in our sample.The titles above each histogram on the diagonal summarize the median and 68% region of each distribution, which are displayed graphically with the vertical dashed black lines.The red line in the [Fe/H] versus [α/Fe] plot near the center of the figure shows where we mask out approximately 2600 stars that lie above this relationship; we do this to keep the 2D abundance distribution singly-peaked.

Figure 3 .
Figure 3.Comparison of best fit model (red line) and normalized flux (blue line) for one RC star with median SNR of 105 pixel −1 .The residuals in the lower left panel and the right panel have been scaled by their corresponding uncertainties.The orange histogram in the right panel is the best fit normal distribution, which shows that this star has good agreement with the data (mean near 0, standard deviation near 1) and that the difference between the data and model are larger than expected by the flux uncertainties at the level of ∼ 15%.The gaps in the data around 15850 Å and 16450 Å correspond to the wavelength gaps between the three APOGEE CCDs.

Figure 4 .
Figure 4. Distribution of uncertainty-scaled residuals for ∼ 3 × 10 8 pixels from all ∼ 5 × 10 5 spectral observations of ∼ 17000 stars in our sample (blue histogram).The orange histogram is the best fit Gaussian, which shows that the residuals are centered near zero and have a standard deviation near 1, as would be expected if the model and flux uncertainties perfectly describe the observed fluxes.Because the best fit Gaussian width is 1.164±0.004,the model doesn't perfectly describe the data within their uncertainties, as is expected if there is non-modelled information remaining.

Figure 5 .
Figure 5. Distribution of the size of residual fluxes that are in excess of the best-fit Gaussian in Figure 4.The grey histograms show different realizations of measuring the excess, while the blue histogram shows the median of these realizations.The blue point at the top of the figure summarizes the median and 68% region of the distributions, revealing that the non-Gaussian excess of residuals is ∼ 3% of the stellar flux on average.

Figure 6 .
Figure 6.Residual fluxes scaled by their uncertainty (i.e.[normalized flux−model]/uncertainty) for all the stars in two different APOGEE pixels.The wavelength of each pixel, in the stellar rest-frame, is listed above the top panel in columns of figures.By highlighting a single pixel, we focus on the residual information from all of the spectra at a particular stellar rest-frame wavelength.Left: A well-modelled pixel, where the uncertainty-scaled flux residuals (blue in top panel) follow the expected unit normal (orange in top panel).In the middle panel, the residuals are binned by heliocentric velocity (number of spectra in each bin given by N * in the legend), and we see no differences between the velocity bins.The bottom two panels show the median and standard deviations of the velocity-binned histograms, again showing no noticeable trend in the residuals with velocity.Right: A pixel where the residuals do not follow a unit Gaussian, or any Gaussian for that matter.After binning by heliocentric velocity in the middle panel, and plotting the binned medians and standard deviations in the bottom two panels, we see a noticeable trend in the median and width of the residual distribution with velocity.

Figure 7 .
Figure 7. Smoothed residual spectra of all RC stars with |VHELIO| < 100 km/s, sorted by the star's VHELIO.Each column corresponds to a residual spectrum and each row corresponds to a pixel/wavelength, with bluest wavelengths at the top and reddest wavelengths at the bottom.This image focuses on 18 different features (wavelengths listed on the right edge of the figure)with each panel consisting of 21 pixels; these cutouts include 15 elemental features, one region we've identified as being mostly continuum, one region of high telluric contamination, and one region near the strongest known DIB in APOGEE (15272 Å feature).The residual spectra have been Gaussian-smoothed in the horizontal direction with a kernel width of 5 km/s.Where applicable, ionization levels are shown.For the C window, the dominant feature is from CO, and in the O window, the dominant feature is from OH.Some elemental regions (e.g. the DIB region, 3rd from the top) show a residual minimum that moves through the element window as a function of heliocentric velocity; this can be explained as an unmodelled residual feature that is not in the same reference frame as the star.

Figure 8 .
Figure 8. Binned heliocentric velocity in the Galactic disk for the RC stars in the APOGEE DR16 sample with |Z| < 0.75 kpc.The dashed circle of radius 8 kpc and the horizontal dashed line intersect at the approximate location of the sun.

Figure 9 .
Figure 9. Tracing the 15272 Å DIB feature location to estimate the velocity offset between DIB sources and RC stars.The smoothed residual data are the same as the 3rd-from-top panel of Figure 7.The white points showing the approximate location of the residual minimum at each pixel in stellar rest-frame, and the orange line shows the best fit line to these data.Using the known rest-frame wavelength of this DIB feature, the orange line is used to infer the VHELIO of the DIB source as a function of stellar VHELIO.

Figure 10 .
Figure10.Simulated comparison of the difference in heliocentric velocity between gas and stars when the intervening gas is at different fractions of the stellar LOS distance.The simulated stellar disk is rotating at 215 km/s while the simulated gaseous disk is rotating at 220 km/s, and measurements are generated for each RC star in Figure8.The orange line is the result of using the known DIB rest-frame wavelength with the the orange line in Figure9to estimate DIB velocity as a function of stellar velocity.The agreement between the orange line and the results of placing the DIB sources between 1% and 50% of the distance to each star suggests that the DIB absorption feature is being produced by intervening material, as we would expect.The tight relationship between ∆VHELIO and VHELIO explains why binning in stellar heliocentric velocity for the RC sample leads to increased DIB residual strength in Figure7.

Figure 11 .
Figure 11.Inverse-variance coadded residual spectra in the observer reference frame in different spectral SNR bins.The blue shaded regions highlight where the strongest telluric features (caused by CO2) occur in APOGEE.The top panel shows the raw coadded residual spectra, and each combined spectrum has a similar median SNR of ∼ 550 pixel −1 ; this panel demonstrates that the residual spikes of the lower SNR bins are larger than those of the high SNR bins.The bottom panel takes the coadded spectra from the top panel and subtracts off the combined spectrum of highest SNR bin (yellowest line in top panel) to better highlight the difference between the spectra of the SNR bins.

Figure 12 .
Figure 12.Fitting the 15272 Å DIB feature with an inverted Gaussian to the residuals of the spectrum in Figure 3.The blue data points show the residual data, the red line shows the best fit inverted Gaussian model with a 68% error envelope, and the dashed vertical black line marks the restframe wavelength of the 15272 Å DIB feature.

Figure 13 .
Figure 13.Equivalent width strength of the 15272 Å DIB feature versus K-band reddening for the spectra in our sample with well-measured 15272 Å DIB features.The different colors correspond to the 10 different EW bins used later in combining the residuals.There are approximately 2100 residual spectra per bin.

Figure 14 .
Figure 14.SNR-modeled residual spectra in the observer rest-frame.The colors of the spectra in this figure correspond to the SNR values in the colorbar above the top panel, and not the EW bins in Figure 13.The blue shaded regions highlight where the strongest telluric features (caused by CO2) occur in APOGEE.

Figure 15 .
Figure 15.Coadded residual spectra in the 15272 Å DIB rest-frame colored by the 15272 Å DIB strength as shown in Figure 13.The vertical red lines show the rest-frame wavelengths of the previously-known DIB features in Table1.There are many newly discovered features in emission and absorption that scale as a function of EW bin, making them new candidate DIB features.Thirty five of the possible DIBs we detect -which we will later measure to have significant correlation with the 15272 Å DIB strength -are shown by the blue points above the spectra with widths showing the region we use as a cutout when measuring equivalent widths.

Figure 16 .
Figure 16.Equivalent width strength of a previously-unknown DIB feature near 15706 Å for the different EW bins shown in Figure 13.The "Window" in the title defines the width of the wavelength cutout.The dashed vertical green line in the right panel shows the location of the local minimum identified in the smoothed highest DIB-strength residual spectrum, which is slightly redward of the feature minimum due to the feature's assymetric shape.The black horizontal line in the left panel and the corresponding grey shaded region are the population mean and width of the current DIB EW measurements (i.e.population fit to the y-axis values).This particular feature's EW shows a fairly strong correlation with the strength of the 15272 Å DIB.
Figure 17.Comparison of the population width distributions of EW measurements around the possible DIB absorption feature shown in Figure16when sorted and binned by different parameters.In each case, we measure the EW of the possible DIB feature from the combined residual spectra in different parameter bins, and then measure the population width of the EW measurements (e.g. the shaded grey region of Figure16).Left: Distributions on the EW population widths for different parameter bins; the black histogram is the resulting average of 5 realizations of random sorting and binning.Right: Distribution on the difference between the population widths in the left panel and the random population width.In this case, the largest population width is achieved after binning by DIB-strength, suggesting that the absorption feature in Figure16is best correlated with 15272 Å DIB-strength.
the 15272 Å DIB rest-frame, of the local extremum found by smoothing the highest DIB-strength spectrum in Figure 15 with a 5 pixel Gaussian kernel.b Hydrogen transition wavelengths that are within 2 Å of the feature's central wavelength.c Region, in the 15272 Å DIB rest-frame, that we measure the equivalent width from.The wavelengths used to estimate the local continuum are taken to be the 5 pixels nearest to the Wavelength Range values.d Equivalent width measured when binning by A K or the 15272 Å DIB strength, in the minimum and maximum bins.* Law, we can find the marginal posterior for p (σ r,i |data) as:p(σ r,i |data) ∝ p(σ r,i ) • σ1

Table 2 .
Pixel bitmasking of spectra a .
a This is a subset of Table1in Price-Jones & Bovy (

Table 3 .
Number of the 133 possible DIBs with best sorting parameter detection probabilities above different thresholds.

Table 4 .
Summary of the 84 DIB features (76 in absorption, 8 in emission) where the best sorting parameter is either AK or 15272 Å DIB-strength.Only features where Pr (σSORT > σRANDOM) > 0.5 are included.