Solid conﬁrmation of the broad DIB around 864.8 nm using stacked Gaia –RVS spectra

Context. Studies of the correlation between di ﬀ erent di ﬀ use interstellar bands (DIBs) are important for exploring their origins. However, the Gaia –RVS spectral window between 846 and 870nm contains few DIBs, the strong DIB at 862nm being the only convinc-ingly conﬁrmed one. Aims. Here we attempt to conﬁrm the existence of a broad DIB around 864.8nm and estimate its characteristics using the stacked Gaia –RVS spectra of a large number of stars. We study the correlations between the two DIBs at 862nm ( λ 862) and 864.8nm ( λ 864.8), as well as the interstellar extinction. Methods. We obtained spectra of the interstellar medium (ISM) absorption by subtracting the stellar components using templates constructed from real spectra at high Galactic latitudes with low extinctions. We then stacked the ISM spectra in Galactic coordinates ( (cid:96), b ) – pixelized by the HEALPix scheme – to measure the DIBs. The stacked spectrum is modeled by the proﬁles of the two DIBs, Gaussian for λ 862 and Lorentzian for λ 864.8, and a linear continuum. We report the ﬁtted central depth (CD), central wavelength, equivalent width (EW), and their uncertainties for the two DIBs. Results. We obtain 8458 stacked spectra in total, of which 1103 (13%) have reliable ﬁtting results after applying numerous conservative ﬁlters. This work is the ﬁrst of its kind to ﬁt and measure λ 862 and λ 864.8 simultaneously in cool-star spectra. Based on these measurements, we ﬁnd that the EWs and CDs of λ 862 and λ 864.8 are well correlated with each other, with Pearson coe ﬃ cients ( r p ) of 0.78 and 0.87, respectively. The full width at half maximum (FWHM) of λ 864.8 is estimated as 1 . 62 ± 0 . 33nm which compares to 0 . 55 ± 0 . 06nm for λ 862. We also measure the vacuum rest-frame wavelength of λ 864.8 to be λ 0 = 864 . 53 ± 0 . 14nm, smaller than previous estimates. Conclusions. We ﬁnd solid conﬁrmation of the existence of the DIB around 864.8nm based on an exploration of its correlation with λ 862 and estimation of its FWHM. The DIB λ 864.8 is very broad and shallow. That at λ 862 correlates better with E ( BP − RP ) than λ 864.8. The proﬁles of the two DIBs could strongly overlap with each other, which contributes to the skew of the λ 862 proﬁle.


Introduction
Diffuse interstellar bands (DIBs) are a set of absorption features with profiles that are much broader than those of the interstellar atomic lines (e.g., Na i lines).To date, over 600 DIBs have been identified at optical and near-infrared wavelengths (0.4-2.4 µm; Galazutdinov et al. 2017, Fan et al. 2019).However, our knowl-edge about their origins is still very limited.So far, only C + 60 has been confirmed as the carrier for five DIBs between 930 and 965 nm (see Linnartz et al. 2020 for a review).Gas-phase carbon-bearing molecules are thought to be the most probable DIB carriers (e.g., Snow 2014).
Despite the unknown nature of DIBs, absorption intensity maps have been built for several such features detected in large spectroscopic surveys (e.g., Kos et al. 2014;Zasowski et al. 2015;Baron et al. 2015b;Lan et al. 2015;Puspitarini et al. 2015) in order to explore the interstellar medium (ISM) and Galactic structure.Among them, the latest map was built for the DIB at 862 nm (DIB λ862) using data from the Gaia Radial Velocity Spectrometer (RVS; Gaia collaboration, Seabroke et al. 2022) in Data Release 3, which has the largest sky coverage to date (see Gaia collaboration, Schultheis et al. 2022 for more details).Their wavelengths and line widths are expressed in nanometers following the unit used in the Gaia analysis, while their strength is still expressed in Ångströms due to their small line depth ( 6%).The mentioned wavelengths in this paper are in the vacuum.
Studies of the correlation between different DIBs are important for identifying their origins (e.g., Friedman et al. 2011;Elyajouri et al. 2017).However, such studies are difficult for λ862 because the RVS spectral window (846-870 nm) has very few DIBs.In the review of possible DIBs in this region by Munari et al. (2008), the most promising one was thought to be located around 864.8 nm (λ864.8)with a very shallow profile.λ864.8 was first reported by Sanner et al. (1978) following positive support from Herbig & Leka (1991), Jenniskens &Desert (1994), andWallerstein et al. (2007).This DIB was clearly seen in spectra of early-type stars, but its analysis is complex because of the wing of the hydrogen Paschen 13 line and other stellar components.On the other hand, this band was not reported by Fan et al. (2019).Krełowski et al. (2019) argued that this feature was of stellar origin (He i lines) as it did not correlate with dust extinction.However, Baron et al. (2015a) confirmed weak correlations between dust extinction and many DIBs, which weakens this argument.Furthermore, a position correlation with dust extinction could be seen around 865 nm in Fig. 2 in Baron et al. (2015a), strengthening the DIB hypothesis for this feature.
In this work, we use the unique set of Gaia-RVS spectra of cool stars observed over the whole sky which has just been published as part of the third Gaia data release (DR3; Gaia Collaboration et al. 2016;Gaia Collaboration, Vallenari et al. 2022).We aim to confirm λ864.8 as a broad DIB rather than a stellar component by the use of stacked Gaia-RVS spectra to detect it and then analysis of its correlations with dust extinction and DIB λ862.

Data
The Gaia-RVS spectra have a wavelength coverage of [846 − 870] nm and a medium resolution of R ∼ 11 500 (Cropper et al. 2018).The stellar parameters and chemical abundances, including effective temperature (T eff ), surface gravity (log g), mean metallicity ([M/H]), and individual abundances of up to 13 chemical species, are derived by the General Stellar Parametrizer from spectroscopy (GSP-Spec) module of the Astrophysical parameters inference system (Apsis, Gaia collaboration, Creevey et al. 2022;Gaia collaboration, Recio-Blanco et al. 2022).We refer to Gaia collaboration, Recio-Blanco et al. (2022) for detailed descriptions of the parametrization and a validation of the measurement of λ862 in individual RVS spectra.
About one million Gaia-RVS spectra have been published in Gaia DR3.We use a sample of 648 944 spectra that meet the following criteria (given as an ADQL query example in Appendix A): 1. Stellar parameters (T eff , log g, [M/H]) and radial velocity (V rad ) are not NaN values.
3. Stellar distances derived from parallax measurements (1/ ) are within 6 kpc. 4. The signal-to-noise ratios (S/N) of the observed spectra are greater than 20. 5.The uncertainty of V rad is smaller than 5 km s −1 .This work processes stars with T eff 7500 K, defined as "cool stars".The whole sample is separated into two parts: 75 941 "reference spectra" at high latitudes (|b| 30 • ) and with low extinctions (E(BP − RP) < 0.02 mag) are used to construct the stellar templates for the cool stars; the other 573 003 are the "target spectra" in which we attempt to fit and measure the two DIBs.The reddening E(BP − RP) used in this work was derived by the GSP-Phot module using their BP/RP spectra (Gaia collaboration, Andrae et al. 2022).Our sample uses about one-tenth of the total RVS spectra processed by GSP-Spec (Gaia collaboration, Recio-Blanco et al. 2022).Kos et al. (2013) and Lan et al. (2015) made use of a similar number of spectra in their studies.Baron et al. (2015a,b)

Method
This work is the first of its kind to fit and measure λ862 and λ864.8 simultaneously in cool-star spectra.To do so, first we derive the ISM spectrum for each target by subtracting their stellar components using the reference spectra.We then stack the ISM spectra of the targets according to their Galactic coordinates ( , b), per HEALPix.Finally we fit and measure the two DIBs, λ862 and λ864.8, in the stacked spectra.The first step follows the main principles of Kos et al. (2013).For a given target, we find a set of reference spectra with similar parameters to the GSP-spec values of the target that have T eff ± 20%, log g ± 0.6 dex, and [M/H] ± 0.4 dex.These ranges are arbitrary, but chosen to be larger than the parameter uncertainties (Kos et al. 2013).The constraints on the stellar parameters are not necessary but they can speed up the procedure.The similarity between the target spectrum and each reference spectrum is then calculated as the average absolute difference of their flux at all pixels, except a masked region of 860-868 nm where the two DIBs are located.When measuring the difference, the central regions of the Ca ii triplet are down-weighted to 30% because Ca ii lines dominate the whole Gaia-RVS spectra of cool stars whereas the similarity between Fe i lines (close to DIBs) is more important than between the Ca ii lines for constructing the templates (see examples in Contursi et al. 2021).In practice, the two Ca ii regions of concern are defined as 849.43-851.03nm and 853.73-855.73nm and are down-weighted.The Ca ii line near 866.5 nm is within the masked DIB region.Reference spectra with average differences of greater than the inverse of the square of the target S/N are discarded.For the rest, up to 500 (and at least 10) best-matching reference spectra are averaged and weighted by their S/N in order to build a stellar template.The target spectrum divided by the template then gives the ISM spectrum for this target.An illustration of this method can be found in Fig. 2 in Kos et al. (2013).Examples of RVS spectra and successful measurements of the DIB λ862 in individual RVS spectrum can be found in Gaia collaboration, Recio-Blanco et al. (2022).
In the second step, the ISM spectra are stacked according to their Galactic coordinates ( , b) to get a sufficiently high S/N for a reliable measurement of λ864.8.The pixelation of the sky is done by the HEALPix1 (Górski et al. 2005) scheme.We choose level 5 (N side = 32), corresponding to 12 288 pixels and a spatial resolution about 1.8 • .We note that 8458 HEALPix pixels (69% of the total) contain at least ten targets that have ISM spectra (N tar > 10).In each of these pixels, we stack all of the ISM spectra (54 on average with a maximum of 362) to generate one stacked spectrum.We stack the ISM spectra in each pixel by taking the median value of the fluxes in order to reduce the influence of outliers.The flux uncertainty of the stacked spectrum is the mean of the individual flux uncertainties divided by √ N tar , which is the standard error in the mean.The S/N of the stacked spectra is calculated between 860.2 and 861.2 nm as mean(flux)/std(flux).We have 8458 stacked spectra in total, one for each HEALPix pixel.
In the third and final step, a Markov chain Monte Carlo (MCMC) procedure (Foreman-Mackey et al. 2013) is used to fit each stacked spectrum between 860 and 868 nm (the DIB region) with a Gaussian function for the profile of λ862, a Lorentzian function (better than Gaussian when considering the goodness of fit to the line wings) for the profile of λ864.8, and a linear function for the continuum placement, masking 866-866.8nm for the Ca ii line residuals.The best estimates and statistical uncertainties are taken in terms of the 50th, 16th, and 84th percentiles of the posterior distribution.We note that the Ca ii region slightly changes after the shifting to the heliocentric frame, con-sidering that the radial velocities of most stars lie mainly within ±50 km s −1 , corresponding to |∆λ| 0.14 nm at 866.5 nm.Furthermore, the wings of the Ca ii lines are better modeled than their central parts (Gaia collaboration, Creevey et al. 2022;Gaia collaboration, Recio-Blanco et al. 2022).Thus, our masked region can properly prevent contamination by the Ca ii residuals and consequently the λ864.8profile can be fitted well.

Results
The aforementioned fit of the DIB was done for each of the 8458 stacked spectra (one per HEALPix pixel).To get reliable results, we only retain the stacked spectra that meet the following criteria, which leaves us with 1962 spectra: -S/N of the stacked ISM spectrum is greater than 100.
-The FWHM of λ862 and λ864.8 are greater than 0.2 nm and 0.5 nm, respectively.
Figure 2 shows the distributions of the FWHM of λ862 and λ864.8, as well as their joint distribution, for the 1962 stacked spectra.The FWHM of λ864.8 has a much wider distribution than that of λ862, and both slightly deviate from a Gaussian.We apply a Gaussian kernel density estimation (KDE) to their joint distribution with a bandwidth of 0.2826 nm (automatically determined by the Python package scipy).The median FWHM is to the upper left of the peak density estimated by the KDE (white star in Fig. 2), which is caused by the fits to some very shallow profiles that result in broader λ862 and narrower λ864.8.The outliers in low-density regions are due to the noise in stacked spectra.
Therefore, we discard the points that lie in the region with a density of less than 1.2 (red line in Fig. 2), which is about onethird of the peak density.After this final filtering, our sample comprises 1103 reliable measurements of the two DIBs, about 13% of the total stacked spectra.
The two-dimensional intensity map of λ862 and λ864.8 as well as E(BP − RP) are shown in Galactic coordinates in Fig. 3. Limited by our sample, the spatial distributions of EW 862 and EW 864.8 cannot be well described.However, it is clear that large EWs for both λ862 and λ864.8 concentrate in the Galactic plane.Furthermore, decreasing EW 862 and EW 864.8 with latitude (up to b = ±30 • ) can be found near Galactic center (GC, | | < 30 • ), which are similar to each other and that of E(BP − RP).We also apply the linear fits without any clipping.For the two DIBs, the coefficient from the no-clipping fit becomes 6% lower (1.548 ± 0.015) for their EWs and 5% lower (0.350 ± 0.006) for their CDs, compared to the 2σ-clipped results.The relative strength between DIB and dust is 13% above for λ862 (E(BP − RP)/EW 862 = 4.096 ± 0.038) and 21% above for λ864.8 (E(BP − RP)/EW 864.8 = 2.367 ± 0.034), respectively, when not clipping.This is not surprising, as significant scatter about the fit can be found around E(BP − RP) ∼ 1 mag and EW 862 ∼ 0.1 Å.Similar scatter exists for λ864.8 as well, but this latter is more significant.This scattering could be a result of the overestimation of E(BP − RP) as seen in Fig. 7 in Gaia collaboration, Schultheis et al. (2022).We prefer the 2σ-clipped results because the scattering is not only caused by the measurement uncertainties but may indicate a spatial disconnection between their carriers and dust grains.Thus, the fitted coefficients with a strong filtering would represent their mean relative strength in the regions where these ISM materials are well mixed.
Compared to EWs, the CDs of the two DIBs are better correlated with each other with r p = 0.87 and no systematic deviations.The depth of λ864.8 is about 37% of that of λ862, while its EW is over 1.6 times EW 862 .We note that EW 862 shows good correlation with E(BP − RP) with r p = 0.85 as expected, but their relative strength is about 20% below that derived from the DIB sample in Gaia-DR3 (4.507 ± 0.137, Gaia collaboration, Schultheis et al. 2022).The reason could be that different templates -from GSP-Spec and from reference spectra-would introduce differences in EW measurements (see e.g., Elyajouri & Lallement 2019).It should also be noted that E(BP − RP)/EW 862 fitted in this work with no clipping is closer to the result in Gaia collaboration, Schultheis et al. (2022, 9% difference) where the authors did not make any filtering either.Also, λ864.8 does not correlate as well as λ862 with E(BP − RP), showing a much smaller r p = 0.64 and a larger dispersion.
In a similar manner to Munari et al. (2008) and Zhao et al. (2021), in order to estimate the rest-frame wavelength (λ 0 ) of λ864.8, we assume that the median radial velocity in the Local Standard of Rest of the DIB carrier is zero toward the GC.There are 43 cases with | | 10 • , |b| 10 • , and small uncertainties in λ C (<0.1 nm for λ862 and <0.15 nm for λ864.8).For λ862, we derive a λ 0 of 862.319 ± 0.018 nm, which is highly consistent with the result of 862.323 ± 0.0019 nm in Gaia collaboration, Schultheis et al. (2022).For λ864.8, we derive λ 0 = 864.53± 0.14 nm in the vacuum, corresponding to 864.29 nm for air wavelength.This result is smaller than the commonly suggested air wavelength measured in the spectra of early-type stars, such as 865.0 nm (Sanner et al. 1978), 864.9 nm (Herbig & Leka 1991), and 864.83 nm (Jenniskens & Desert 1994).

Discussion
Until now, the FWHM of the λ864.8feature has not been well determined.Sanner et al. (1978) noted that this DIB profile extends from 863 to 866 nm and its strength did not correlate with spectral type, but they did not measure its FWHM.We found two measurements of FWHM λ864.8 in the literature that are dramatically different from each other: 1.4 nm by Herbig & Leka (1991) and 0.42 nm by Jenniskens & Desert (1994).The difficulty is that this band is shallow and superposed across several blended stellar lines, such as the He i line at 864.83 nm in early B-type supergiants (Herbig & Leka 1991).This may be the reason why a stellar origin has been claimed in the literature (Krełowski et al. 2019).In this work, we made use of the spectra of cool stars, and  subtracted the stellar lines (mainly Fe i and Ca ii) using reference spectra.This makes our measurements more accurate than those previous works based on hot-star spectra.We find λ864.8 to be a very broad and shallow DIB with a FWHM of 1.62 ± 0.33 nm, which is not significantly less than the strongest DIB at 442.8 nm (e.g., 1.7 nm; Galazutdinov et al. 2020).This provides strong evidence for the DIB interpretation of this line because no known stellar components or Doppler broadening (a velocity dispersion over 600 km s −1 ) can explain such a large width.On the other hand, the derived FWHM of λ862 in this work is 0.55 ± 0.06 nm, which is slightly larger than previously reported values such as 0.43 nm (Herbig & Leka 1991), 0.47 nm (Maíz Apellániz 2015), and 0.37 nm (Fan et al. 2019).The increase in FWHM could be explained by a Doppler broadening caused by our stacking strategy, with a velocity dispersion about 30-50 km s −1 .This value is consistent with the average magnitude of the V rad dispersion in each pixel (about 36 km s −1 ) which can be treated as an upper limit to the velocity dispersion of the DIB clouds, because closer objects should have a smaller radial velocity than the background stars.
Although the DIB profile could be broadened due to the superposition effect, our stacking has only a very weak influence on EW measurements, because DIBs are weak spectral features and we are therefore in a linear regime where the EW produced by multiple DIB clouds can be simply added up (Munari & Zwitter 1997).The linear correlation between EW 862 and EW 864.8 provides good evidence for a DIB interpretation of this line, in spite of the superposition effect, because the stellar components or their residuals cannot have such a behavior with interstellar features.Another important point is that the profiles of the two DIBs could be overlapped with each other.Therefore, it is better to measure the two DIBs simultaneously; although this requires a very high S/N.
Limited by our sample and the conservative approach to filtering it, further analysis of the two DIBs -such as the kinematics of their carriers and the nature of the Lorentzian profile of λ864.8isnot feasible.Nevertheless, it is significant to confirm a second DIB close to the strong one λ862 in such a big spectroscopic survey, which enables us to use them as tracers of the ISM environment if they have different origins, or facilitates the study of their origin if they share a common carrier (with two DIBs, their relative strength and central wavelengths could provide clues as to their carrier).Dedicated research into these two DIBs will be carried out using all available RVS spectra and the results will be published in next Gaia data release.

Conclusions
Based on the measurements in 1103 stacked Gaia-RVS spectra, we provide solid confirmation of the DIB around 864.8 nm through its correlation with λ862 and its clear and broad profile in the stacked spectra.λ864.8 is a very broad and shallow DIB.The FWHM of λ864.8 is estimated to be 1.62 ± 0.33 nm.Using 43 high-quality measurements toward GC (| | 10 • , |b| 10 • ), the rest-frame wavelength of λ864.8 is determined as λ 0 = 864.53± 0.14 nm in the vacuum, which is smaller than previous reported measurements.EW 862 correlates better with E(BP − RP) than EW 864.8 .
Our work shows the power of using a large set of cool-star spectra to study the DIBs in the ISM.The extremely small depth of these lines (CD 864.8 3%) and the ability to assess their properties is a clear demonstration of the quality of Gaia-RVS spectra as published in Gaia DR3.

Fig. 1 .
Fig. 1.Four examples of the fits to DIBs λ862 and λ864.8 in stacked Gaia-RVS spectra.The black lines are the ISM spectra normalized by the fitted linear continuum.The error bars indicate the flux uncertainties at each pixel.The red lines are the fitted curves that are normalized by the continuum as well.Orange marks the Ca ii line residuals.The median E(BP − RP) and EWs of the two DIBs in each HEALPix are also indicated.

Fig. 2 .
Fig. 2. Distributions of the FWHM of λ862 (red histogram) and λ864.8 (blue histogram), as well as their joint distribution (middle colored map), measured in 1962 stacked spectra after applying the general filtering.The colors represent the densities of the joint FWHM distribution, estimated by a Gaussian KDE.The white star indicates the peak density.The red line in the central panel indicates a contour with a probability density of 1.2, about one-third of the peak density.The orange lines indicate the median FWHM of λ862 and λ864.8,respectively.

Fig. 3 .
Fig. 3. Galactic distributions of the 1103 reliable measurements of the strength of λ862 (top), λ864.8 (middle), and E(BP − RP) (bottom), in Mollweide projection, with the Galactic center in the middle and Galactic longitude increasing to the left.

Fig. 4 .
Fig. 4. Intensity correlations between λ862, λ864.8, and E(BP − RP).Different correlations are shown in each panel with the Pearson coefficient (r p ) and the slope of a 2σ-clipped linear fit (α) with zero intercept shown in red.Data points with error bars are colored by their number density estimated by a Gaussian KDE.The black dashed line shown in the second panel from the top is the linear relation derived in Gaia collaboration, Schultheis et al. (2022).
Article number, page 4 of 7 H.Zhao (赵赫) et al.: Solid confirmation of the broad DIB around 864.8 nm using stacked Gaia-RVS spectra