© The Authors 2024

1 Introduction

Diffuse interstellar bands (DIBs) are a set of broad interstellar absorption features that are ubiquitously observed from the optical to the near-infrared (NIR) wavelength range (Fan et al. 2019; Ebenbichler et al. 2022; Hamano et al. 2022; Vogrinčič et al. 2023; Castellanos et al. 2024). Since the first discovery of the DIBs λ5780 and λ15797¹ by Heger (1922), many efforts have been dedicated to identifying the carriers of DIBs (e.g., Iglesias-Groth 2007; Krelowski et al. 2011; Zanolli et al. 2023), but with the exception of ${\text{C}}_{60}^{+}$ (Campbell et al. 2015), which has been found to be responsible for five DIBs in the NIR wavelength range (see a review by Linnartz et al. 2020), the carriers of DIBs remain unknown.

In attempts to identify the carriers of DIBs and guide laboratory searches for them, attention has been turned to the characterization of DIBs based on astronomical observations, such as their intrinsic profiles in high-resolution spectra and their correlations with other tracers of the interstellar medium (ISM). For example, MacIsaac et al. (2022) constrained the sizes of possible carrier molecules by comparing the peak-to-peak separation in the DIBs λ5797, λ6379, and λ6614 using the extensive spectra from the ESO diffuse interstellar bands large exploration survey (EDIBLES). To better understand the behavior of DIBs in different interstellar environments, Fan et al. (2017) analyzed how the pattern of the equivalent widths (EWs) of DIBs – normalized by the corresponding extinction – changes with the molecular fraction f_H2.

The first observed DIBs, and those seen to have the strongest signal (Merrill 1930), the DIBs λ5780, λ5797, and λ6614, as well as their interrelations, have been extensively studied and discussed in the literature. Owing to the proximity of the DIBs λ5780 and λ5797, they are usually investigated as a pair. The well-known “σ − ζ” effect was found from the varying strength ratio of these two prominent DIBs (e.g., Krelowski & Westerlund 1988; Vos et al. 2011; Kos & Zwitter 2013). These works indicated that the EWs of the DIB λ5780 in the “ζ–type” sight lines tend to be weaker than those in the “σ–type” sight lines, while the corresponding EWs of the DIB λ5797 are similar. Upon discovering that the relationship between DIBs and reddening is linear, but with considerable dispersion (including the DIBs λ5780, λ5797, and λ6614; Raimond et al. 2012; Puspitarini et al. 2013), and that there exists a “skin effect” (Snow & Cohen 1974; Herbig 1995), the correlations with other ISM tracers, such as OH, CH, CO, H I and H₂, have been explored (Weselak et al. 2004; Friedman et al. 2011; Lan et al. 2015; Fan et al. 2017; Weselak 2019). It is now well established from those studies that the DIB λ5797 is more closely related to the column density of H₂ than DIB λ5780, while this latter shows a stronger correlation with that of H I. In addition, through studies of the correlations between DIBs, an almost perfectly correlated DIB pair has been identified, namely DIB λ6196-λ6614 (Cami et al. 1997; Moutou et al. 1999; Galazutdinov et al. 2002; McCall et al. 2010; Bailey et al. 2016), as have some relatively well correlated DIB pairs, such as DIB λ5780-λ6614 and DIB λ5797-λ5850 (Friedman et al. 2011; Smith et al. 2021; Fan et al. 2022), and even some anticorrelated DIB pairs, such as DIB λ4984-λ7559 and λ5418-λ7562 (Fan et al. 2022), which can help to group the DIBs into specific DIB families. Furthermore, in high-resolution spectra, the longward blue wing in the DIB λ5780 profile (Snell & vanden Bout 1981; Danks & Lambert 1976), the asymmetry profile in the DIB λ5797 (MacIsaac et al. 2022), and the triple-peak fine structure in the DIB λ6614 profile (Herbig & Soderblom 1982; Galazutdinov et al. 2008), all suggest that their carriers are complex macromolecules.

However, the majority of the aforementioned studies are mainly based on the high-resolution spectra of early-type stars without contamination from the stellar lines. It is challenging to carry out larger-scale and more statistically significant analyses as a result of the limited sight lines toward the early-type stars. Kos et al. (2013) introduced a best neighbor matching method (BNM) to extract the DIB λ8621 from approximately 500000 spectra of late-type stars in the Radial Velocity Experiment survey (RAVE; Steinmetz et al. 2006) and provided a statistical ratio between DIB λ8621 and extinction. Although the objects observed in the Sloan Digital Sky Survey (SDSS; Yanny et al. 2009) do not concentrate on the Galactic plane where the reddening is high, Lan et al. (2015) still managed to measure over 20 DIBs from the spectra of cool stars by stacking the spectra with similar extinction, and mapped out the distribution of DIBs at high latitudes in the Milky Way. Later on, using the cool stars from the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017), Zasowski et al. (2015) presented a projected Galactic map of the single DIBs at 1.5273 µm and were the first to estimate its precise central wavelength based on kinematics information. Similarly, Gaia Collaboration (2023b) and Gaia Collaboration (2023a) further elaborated the properties of the DIBs λ8621 and λ8648 in the Gaia radial velocity spectrometer (RVS) spectra (Cropper et al. 2018; Seabroke et al. 2021) of both cool and hot stars. Furthermore, Cox et al. (2024) constructed a 3D Galactic map tracing the DIB λ8621 based on the measurements of Gaia Collaboration (2023b). Using the information provided by the GALactic Archaeology with HERMES (GALAH; Buder et al. 2021) for approximately 872 000 cool stars, brand new DIBs were found in the residuals of these stellar spectra (Vogrinčič et al. 2023). Moreover, Saydjari et al. (2023) and Zhao et al. (2024) recently employed the methods based on machine learning to improve measurements of DIBs from the spectra of cool stars. To measure DIBs from cool stellar spectra, the stellar lines must be removed from spectra in order to extract the ISM composition. A critical distinction exists between the methods used by authors to obtain the stellar lines. On the one hand, studies such as those of Chen et al. (2013), Puspitarini et al. (2015), Zhao et al. (2021a), and Gaia Collaboration (2023b) use theoretical methods to synthesize the stellar lines. On the other hand, observed spectra with low or no extinction are used to pairwise match the stellar components (e.g., Cordiner et al. 2008a,b; Kos & Zwitter 2013; Kos et al. 2014; Zasowski et al. 2015 and Gaia Collaboration 2023a), which is also the approach adopted in the present study.

In addition to the works of Yuan & Liu (2012) and Lan et al. (2015), which are based on low-resolution spectra (R ≈ 2000) obtained by SDSS, certain teams have also made efforts to glean hidden information from low-resolution spectra. Using spectra with a resolving power of R = 3300, Cordiner et al. (2008a) investigated the DIBs λ5780, λ5797, λ6203, λ6283, and λ6614 in the ISM of M31. These authors then used the spectra with similar resolution to extend their work to the DIBs in the ISM of M33 (Cordiner et al. 2008b), and carried out the first survey of DIBs observed in the spectra of B-type supergiants in M31 (Cordiner et al. 2011). By means of the DIBs λ5780 and λ5797 measured from spectra (Farhang et al. 2015a) with R = ~ 2000 in the northern hemisphere, Farhang et al. (2015b) probed the Local Bubble and its surroundings. Later, combining the DIB measurements derived from the spectra with a resolution of 5500 in the southern observations, they presented a 3D map of the hot Local Bubble (Farhang et al. 2019). Although the quantity of low-resolution spectra used in these studies has risen tremendously compared to high-resolution spectra, it is still insufficient to support a broader perspective for investigating the properties of DIBs, such as their spatial distribution and kinematics in the Milky Way.

In the present work, we take full advantage of the LAMOST low-resolution survey (LAMOST LRS; Zhao et al. 2012; Deng et al. 2012; Luo et al. 2015). The survey almost covers the entire northern sky from the declination of -10° to nearly 90°, with particularly dense coverage of the ISM-rich Galactic plane. So far, a few million spectra have been obtained for cool stars across the optical wavelength range, which allow us to detect several DIBs simultaneously. We selected the DIBs λ5780, λ5797, and λ6614 as our subjects, because their higher intensities compared to the other DIBs (Fan et al. 2019; Vogrinčič et al. 2023) mean that we are more likely to be able to measure them in the individual spectra of cool stars. Although the DIBs λ4430 and λ6283 are also strong, they are not included in this work because of the complexity of the spectral region in which they are found. For the DIB λ4430, it is hard to determine the continuum level due to the adjacent strong hydrogen Balmer lines and its own broad profile.

For the DIB λ6283, the O₂ telluric band wraps around it (Herbig 1975; Jenniskens & Desert 1994) and is not cleanly removed in the LAMOST spectra. As a consequence, we focus on the DIBs λ5780, λ5797, and λ6614, build the largest sample of measurements of these three DIBs to date, and conduct a comprehensive and statistically robust analysis of their properties.

The remainder of this paper is organized as follows. In Sect. 2, we describe the data used in this work. The pipeline for measuring the DIBs λ5780, λ5797, and λ6614 in the spectra of cool stars is detailed in Sect. 3. In Sect. 4, we control the quality of our measurements and carry out a series of validation tests to ensure the reliability of the measurements. We present and discuss our results in Sect. 5. Finally, we summarize our findings in Sect. 6.

2 Data

LAMOST (the Large Sky Area Multi-Object-fiber Spectroscopic Telescope; Zhao et al. 2012; Deng et al. 2012; Luo et al. 2015) is a reflecting Schmidt telescope with a 5° field-of-view and an effective aperture of 4 meters. Thanks to 4000 fibers averagely distributed among 16 spectrometers (i.e., each one accepts 250 fibers), LAMOST is able to simultaneously observe 4000 objects in theory. There are two resolving modes for LAMOST: one is the low resolution (R = 1800 at 5500 Å) and the other is the medium resolution (R = 7500 at 5163 Å and 6593 Å). When switching to the low-resolution mode employed here, the designed optical band covers the range from 3700 Å to 9000 Å.

The data used in this work are from the data release 10 of LAMOST low-resolution spectroscopic survey (LAMOST LRS DR10)². A total of 11817430 spectra are published in this release, which contains 11 473 644 stellar spectra, 263 444 galaxy spectra, and 80342 quasar spectra. All the published spectra are reduced by wavelength calibration, sky background subtraction, telluric correction, and relative flux calibration (Zhao et al. 2012; Luo et al. 2015). Note that the telluric absorption is negligible in the regions of the DIBs λ5780, λ5797, and λ6614 investigated in this work (see Fig. 16 in Matheson et al. 2000 and Fig. 2 in Kimeswenger et al. 2015 for more details). All the LAMOST targets are cross-matched with Gaia DR3 sources (Gaia Collaboration 2023c) within a radius of 3 arcsec, and the Gaia IDs, namely gaia_source_id, are included in the official LAMOST DR10 catalog³. The stellar spectra we focus on are corresponding to ~2 000 000 unique stars and the majority of them concentrate on the Galactic plane. For the cool stars, that is, the F, G, K, and late A-type stars, the LAMOST stellar parameter pipeline (LASP) provides the effective temperature (T_eff), surface gravity (log ɡ), metallicity abundance ([Fe/H]), and radial velocity (RV) with the precision of 100 K, 0.19 dex, 0.14 dex, and 6 km s⁻¹, respectively⁴. In addition, the continuum-normalized flux are calculated for these spectra with stellar parameters (Song et al. 2012).

Benefiting from the estimations of geometric distances from parallaxes in Gaia DR3 (Bailer-Jones et al. 2021), the relatively reliable distance of an object in LAMOST can be simply derived. Combined with the three-dimensional dust reddening map (version name “bayestar2019”, Green et al. 2019), once the celestial coordinate and astrometric distance of an object is given, the fine-grained extinction E(B – V)⁵ can be calculated at the sight line of each object.

3 Method

Figure 1 shows an overview of the pipeline for measuring the DIBs λ5780, λ5797, and λ6614 of cool stars in the LAMOST LRS DR10. Four swimlanes in the figure represent the four modules of the pipeline, which are responsible for building the target and reference dataset, deriving the ISM residual spectra of the target dataset, measuring the DIBs in the residual spectra, and quality control of the final measurements. The details of each module are described in the following subsections.

3.1 Building the dataset

For the cool stars, the regions nearby DIBs are usually blended with the stellar absorption lines, which need to be thoroughly removed as much as possible. Inspired by the cases for the success in detecting the DIBs in the spectra of cool stars (Kos et al. 2013; Yuan & Liu 2012; Lan et al. 2015; Vogrinčič et al. 2023), we build two datasets, that is, the target dataset in which the DIBs are likely to be embedded in the spectra, and the reference dataset in which the spectra are expected to contain the pure stellar features within the region of the DIB to be measured.

At first, both the target and reference datasets are filtered by the basic data quality control for the better stellar parameter matching in Sect. 3.2. They are required to have valid stellar parameters within the error of 100 K, 0.2 dex, 0.1 dex, and 10 km s⁻¹ for T_eff, log ɡ, [Fe/H], and radial velocity (RV), respectively.

In general, given that the absorption depths of the DIBs λ5780, λ5797, and λ6614 are commonly larger than one percent of continuum (Vogrincic et al. 2023), we conservatively select the spectra with the signal-to-noise ratio in r band (the keyword snrr in the LAMOST catalog, hereafter S/N, which is derived from the average signal-to-noise ratio of the pixels within the r band, i.e., from 5600 to 6800 Å) greater than 100 as our target dataset. In that the reference dataset are used to subtract the components that do not contribute to the DIB in the target dataset, they are supposed to meet the following conditions to ensure that the interstellar clouds in the foreground are as thin as possible: (i) S/N > 50; (ii) E(B – V) < 0.03 mag; (iii) Galactic latitude, |b| > 60°. There are 2 188 240 spectra in the target dataset and 424 164 spectra in the reference dataset after the above filtering. The distribution of three atmospheric parameters for the two datasets are shown in Fig. 2. The JS divergences (Menéndez et al. 1997), which range from 0 to 1 and are often used to measure the similarity between two distributions, of T_eff, log ɡ, and [Fe/H] between the target and reference dataset are 0.0589, 0.0410, and 0.0278, respectively. These values closer to zero indicate it is of high similarity between the two datasets, which makes the reference dataset matchable for the subtraction of the stellar features in the target dataset.

3.2 Getting the ISM residual spectra

An ISM residual spectrum is the integral signals of DIB at the line-of-sight of a target object, so that an observed spectrum in the target dataset can be viewed as the product of the ISM residual and the stellar components. Although it is impossible to resolve the authentic stellar features out of a target spectrum, one of the workarounds is stacking the neighboring spectra of the target in the reference dataset to approximate the stellar ingredients we call the template spectrum (Kos et al. 2013). Then, the ISM residual spectrum can be obtained by the ratio of the target spectrum to the template spectrum from the stacked closest spectra, as can be seen in Fig. 3.

Fig. 1 Schematic flowchart of the pipeline for measuring the DIBs λ5780, λ5797, and λ6614 of cool stars in the LAMOST LRS DR10. Four colors encode the four parts of the pipeline, and the output of each modular function is used for the input of the subsequent one.

Fig. 2 Distribution of the three atmospheric parameters, namely, Teff (left panel), log ɡ (middle panel), and [Fe/H] (right panel), for the target and reference datasets. The JS divergences of the three parameters between the two datasets are annotated in each panel.

3.2.1 Matching neighbors

The quality of the best neighbor matching determines whether the stellar components in the target spectrum can be cleanly deducted. It is essential to ensure that these neighbors are as close to the target possible. We impose constraints on the selection of the best neighboring samples based on the atmospheric parameter space, match number, S/N, and similarity in spectral morphology.

We make a preliminary selection within a certain range of atmospheric parameters. For the samples in the reference dataset that satisfy ΔT_eff < 100 K, Δlog ɡ < 0.5 dex, and Δ[Fe/H] < 0.2 dex with respect to the target, if the number of those samples is no more than 100, we drop them out. Then, to balance the computation efficiency and match accuracy, the top 1000 samples (if enough) with the highest S/N in r band as the neighboring candidates of the target are retained.

However, the proximity of atmospheric parameter space cannot guarantee the similarity of spectral morphology due to the measurement error of parameters. We further adopt the reciprocal of Euclidean distance between the target and the neighbor to represent the similarity in spectral morphology. Each continuum-normalized spectrum is processed for the alignment of similarity computation as follows: (i) mask the bad flux pixels according to the ormask flag; (ii) shift the spectrum to the rest frame; (iii) convert the wavelength from vacuum to air; (iv) truncate the spectrum to [4500, 7500] Å to avoid the bias of Balmer series at the blue end and the low response efficiency at the red end, and then rebin the spectrum with a step of 1 Å; (v) mask the regions that do not contribute to the metal lines, including Hα, Hβ, strong DIBs λ5780, λ5797, λ6283, and λ6614 (see Table for their detailed wavelength ranges). Ultimately, the top 25% similar candidate samples, namely those with no fewer than 25 and at most 250, are selected as the closest neighbors prepared for the stacked template spectrum in Sect. 3.2.3. For instance, the first 30 closest neighbors of a target spectrum of obsid = 837603202 can be seen in the second panel of Fig. 3.

Fig. 3 Example measurement of the DIBs λ5780 and λ5797 from a target spectrum. From top to bottom, (i) the first panel shows the target spectrum to be measured, with the red text indicating its obsid (the unique ID of the spectrum in the LAMOST catalog) and atmospheric parameters. (ii) The second panel provides the first 30 best neighboring spectra of the target spectrum, which are used to build the template spectrum displayed in the third panel. The purple-shaded areas in the first and second panels are the masked regions that do not participate in the similarity computation. (iii) The third panel illustrates the target and template spectrum after local normalization. The most prominent stellar absorptions that may affect the DIB extraction are annotated in the panel. (iv) The fourth panel presents the ISM residual spectrum (in black) of the target and the fit (in red) of the DIBs λ5780 and λ5797 by MCMC. The shallow DIBs λ5762 and λ5815 in the ISM residual spectrum are also annotated but not fitted in the panel.

3.2.2 Preprocessing

Similar to Sect. 3.2.1, a series of preprocessing steps are required to align the target and its neighbors before extracting the ISM residual spectrum of the target. These preprocessing steps still include removing the bad flux pixels, shifting the spectrum to the rest frame, and converting the wavelength from vacuum to air. Unlike Sect. 3.2.1, the renormalized flux within the local region instead of the continuum-normalized flux is utilized. Although the continuum-normalized flux provided by LASP meets the requirement of spectral morphology matching in Sect. 3.2.1, it is too coarse to measure the weak DIBs when anchoring the local region of the DIBs. Therefore, a more careful local normalization for the flux within the local DIB region needs to be performed.

Initially, the wavelength range of the local DIB region is determined. For the DIB λ5780 and λ5797, they are merged into a single local region, namely [5755, 5855] Å owing to their close central wavelengths. For the DIB λ6614, the blue end keeps away from the Hα line to alleviate the influence of the quite strong Hα feature on the blue wing of the DIB 26614. In addition, the red end of the local region is moderately extended to compensate for the accuracy of the local continuum normalization. The final selection for the DIB 26614 is [6590, 6690] Å. Subsequently, the flux in these two local regions is linearly interpolated and rebinned to an interval of 0.8 Å, resulting in a total of 125 pixels. We refer to the method by Zhao et al. (2021a) as the local normalization of the flux within the selected local DIB region. The only difference is that we use a fifth-order polynomial curve instead of a second-order one to fit the local spectrum because of the cool stars in this work. A 20-iteration continuum fitting is performed as follows: (i) fit the local spectrum by a fifth-order polynomial curve; (ii) compute the differences of the local spectrum to the fitting curve, as well as their standard deviation; (iii) replace the flux values with the fitted polynomial ones if the differences between them are larger than five times of the standard deviation when the pixels are located above the polynomial curve, or 0.5 times when the pixels are below the polynomial curve; (iv) use the remaining and replaced pixels as the new local spectrum to repeat the above steps. The final normalized spectrum is obtained by dividing the last local spectrum by the last fitted continuum. Figure 4 provides two examples of the local renormalization for the target spectra with obsid = 894515122 and obsid = 794616199. It is clear that the local renormalization can effectively correct the continuum level of some absorption features.

The flux error is estimated while performing the local normalization because the flux error after local normalization is not only a key input parameter for the subsequent MCMC fitting in Sect. 3.3, but also an important indicator for evaluating the quality of the ISM residual spectra in Sect. 4. Although the LAMOST provides the error of the flux in a format of the inverse variance (1/error²), many works, such as Xiang et al. (2015) and Ho et al. (2017) have reported this error is overestimated. Ho et al. (2017) and Zhang et al. (2020) have given the same approximations of the inverse variance of the flux after normalization, which can be expressed as Eq. (1), $${\text{ivar}}_{\text{norm }}={\text{F}}_{\text{cont }}^2\times {\text{ivar}}_{\text{rel}},$$(1)

where F_cont is the fitted continuum and ivar_rel is the inverse variance of the flux. The final renormalized flux error equals the inverse of the square root of ivar_norm.

Fig. 4 Examples of the local renormalization for two target spectra. The upper panel is the target spectrum within the local region covering the DIBs λ5780 and λ5797, and the lower panel is the target spectrum within the local region including the DIB λ6614. In each panel, the black line represents the original spectrum, and the orange line is the renormalized spectrum. The obsid of the target spectrum and its atmospheric parameters are also annotated in each panel.

3.2.3 Isolating the ISM residual spectra

Once the locally normalized spectra of the target and its best neighbors are ready, the high-quality template spectrum for the target can be obtained by averaging those spectra of the best neighbors weighted by their S/N in r band, which characterizes the stellar features in the target spectrum. Then, the ISM residual spectrum of the target is isolated by dividing the target spectrum by the template spectrum, as shown in the third and fourth panels of Fig. 3. By comparing the template spectrum to the spectral standard star Arcturus (a K1.5 III star with the luminosity of 170 L_⊙) using an interactive database (Lobel 2011) called SpectroWeb (http://spectra.freeshell.org/spectroweb.html), we also list the most prominent stellar lines that may contaminate the DIBs in the template spectrum in the third panel of Fig. 3.

The error of the ISM spectrum is deduced by the error propagation as Eq. (2), $${\text{E}}_{\text{res}}=\frac{{\text{F}}_{\text{tar}}}{ {\text{F}}_{\text{temp}}}\times \sqrt{{\left(\frac{{\text{E}}_{\text{tar}}}{ {\text{F}}_{\text{tar}}}\right)}^2+{\left(\frac{{\text{E}}_{\text{temp}}}{ {\text{F}}_{\text{temp}}}\right)}^2},$$(2)

where F_tar and F_temp are the normalized flux of the target and template spectra, respectively, and E_tar and E_temp are their corresponding errors. For the E_temp, it is calculated by the same operation as that of F_temp, that is, the weighted average of the flux errors of the best neighbors.

3.3 Measuring diffuse interstellar bands

Most of the DIB profiles are fitted by the symmetric Gaussian function, such as Lan et al. (2015) and Zhao et al. (2021a), but some studies based on the early-type stars in high resolution have found that the profiles of some DIBs are asymmetric and multiple components. For instance, Galazutdinov et al. (2008) reported the DIB λ6614 displays a triple-peak fine structure. Considering the low resolution of LAMOST, the instrumental broadening effect (see Appendix A) dominates the shape of the DIBs λ5780, λ5797, and λ6614 and their fine structures are blurred. In such a resolution, Yuan & Liu (2012), Lan et al. (2015), and Farhang et al. (2019) have proven the Gaussian function is still a good approximation for the profile of those DIBs. Thus, we adopt the Gaussian function to fit the DIBs in the ISM residual spectra. The process of measuring the DIBs consists of three steps: (i) pre-detect the central wavelength and depth of the DIB to determine whether there is a DIB signal; (ii) use the parameters obtained in the first step and the fixed Gaussian widths (3.85, 2, and 2 Å for the DIBs λ5780, λ5797, and λ6614, respectively) as the initial values of the least-square curve fitting to get the three parameters of the Gaussian profile, namely, the central wavelength, depth, and Gaussian width; (iii) further optimize the parameters obtained in the second step using the MCMC method to get the final fitting results. An MCMC fitting demo of the DIBs λ5780 and λ5797 is presented in the fourth panel of Fig. 3.

3.3.1 Pre-detection

In the context of big data mining, pre-detection not only allows the rapid identification of samples that probably contain DIBs, but also provides reasonable initial values for the curve fitting in Sect. 3.3.2. We take the rest-frame central wavelength and FWHM of the DIBs reported by Fan et al. (2019) into account, and set the wavelength range of the pre-detection segment to be about ±300 km s¹ around the central wavelength of the DIB. For the DIBs λ5780, λ5797, and λ6614, their pre-detection segments are [5775, 5785], [5793, 5801], and [6608, 6620] Å, respectively. Such a wide wavelength range is chosen to deal with those potential DIBs with large RVs or broad profiles. If the maximum depth of the ISM residual spectrum within the pre-detection segment exceeds $\frac{1}{ \text{S}/\text{N}}$ of the target spectrum, where S/N is the average signal-to-noise ratio within the r band that covers the DIBs λ5780, λ5797, and λ6614, we proceed to the next step of curve fitting. The maximum depth and the wavelength corresponding to it are taken as the initial depth and central wavelength of the DIB profile. Otherwise, we drop out the target.

3.3.2 Curve fitting

Both curve fitting and MCMC fitting in Sect. 3.3.3 are done in a narrower region called the “focus region”, which is dynamically determined according to the central wavelength of the DIB pre-detected in Sect. 3.3.1. The flux values beyond the focus region are set to the continuum level to prevent the other potential features, such as certain weaker DIBs (see the bottom panel of Fig. 3), from disturbing the fitting. The wavelength interval of the focus region is designed based on the 2σ principle, where σ refers to the Gaussian width. For instance, if a DIB has a σ of 2 Å, its entire profile will span approximately 2 × 3σ, namely 12 Å. Therefore, a wavelength interval of at least 2 × 2σ = 8 Å is necessary for a reasonable fitting. Specifically, for the DIBs λ5797 and λ6614, the width of the focus region is set to 6 Å on both sides of the pre-detected central wavelength. For the DIB λ5780, the width is set to 10 Å because the typical FWHM of the DIB λ5780 is larger than those of the DIBs λ5797 and λ6614 (Hobbs et al. 2008; Fan et al. 2019). On top of the above cases, if the DIBs λ5780 and λ5797 are both pre-detected, the blue end of the focus region is set to 6 Å on the left side of the pre-detected central wavelength of the DIB λ5780, and the red end is set to 6 Å on the right side of the pre-detected central wavelength of the DIB λ5797. Similar to the wavelength interval of the predetection in Sect. 3.3.1, such a focus region is wide enough to handle the broadened DIBs, which may be caused by multiple DIB components or the instrumental broadening effect.

We employ the method of least-square curve fitting (the function curve_fit in the Python module scipy.optimize) to obtain the central wavelength, depth and FWHM of the single Gaussian profile of the DIB as Eq. (3). $$f_{\theta }(x;\text{D},\mu ,\sigma )=\text{D}\times \exp \left(-\frac{{(x-\mu )}^2}{2{\sigma }^2}\right)+1,$$(3)

where x means the wavelength, and D, µ, σ denote the depth, central wavelength, and Gaussian width, respectively. In the case where the DIBs λ5780 and λ5797 are both detected, a double Gaussian profile is adopted to simultaneously fit the two DIBs, which can overcome the blend of the two close DIBs when the DIBs are greatly broadened. The initial values of the depth and central wavelength are taken from the pre-detection results in Sect. 3.3.1, and the initial values of the Gaussian width are set to 3.85, 2, and 2 Å for the DIBs λ5780, λ5797, and λ6614, respectively. Such initial values are larger than the typical values of the DIBs owing to the instrumental broadening effect in the low resolution (Yuan & Liu 2012; Lan et al. 2015). Actually, there is little influence on the final fitting results while using the canonical Gaussian width, for instance, the FWHM from Fan et al. (2019), as the initial value.

3.3.3 MCMC fitting

The prior in the Bayesian inference can rule out some cases that obviously violate the physics, for example, the central wavelength of DIB is far away from the well-known value. And the MCMC method is able to efficiently sample the multi-parameter distribution from the posterior constructed by the product of the likelihood and prior based on the Bayesian theory (Eq. (4)). $$\text{P}(\theta \mid \text{Y})\propto \text{P}(\text{Y}\mid \theta )\times \text{P}(\theta ),$$(4)

where θ represents the parameters of the DIB profile, that is, θ = (D,µ, σ), and Y denotes the observed data. Specifically, the likelihood and prior are expressed as Eqs. (5) and (6) in logarithmic form, respectively. $$\ln \text{P}(\text{Y}\mid \theta )=-\frac{1}{2}{\displaystyle \sum \left[\frac{{\left(y-f_{\theta }(x)\right)}^2}{y_{\text{err}}^2}+\ln \left(y_{\text{err}}^2\right)\right]},$$(5)

where x, y, and y_err as the observed data Y are the wavelength, flux, and flux error of the ISM residual spectrum, respectively, and f_θ(x) is the model of DIB profile which is set out in Eq. (3). $$\ln \text{P}(\theta )=\ln \text{P}(\text{D})+\ln \text{P}(\mu )+\ln \text{P}(\sigma ),$$(6)

where P(D), P(µ), and P(σ) are the priors of the depth, central wavelength, and Gaussian width, respectively. The prior of D satisfies the uniform distribution, and the priors of µ and σ are both set to be the Gaussian distributions with the means of the curve fitting results in Sect. 3.3.2 and the standard deviations of 0.5 Å. Table 2 summarizes the detailed upper and lower limits of the priors of D, µ, and σ for the DIBs λ5780, λ5797, and λ6614.

The Python package emcee is called to carry out the MCMC fitting. We set 100 walkers and 250 steps for each walker, where the first 50 steps are used as the burn-in stage and the last 200 steps are used to sample the posterior distribution. The initial guesses of the walkers are randomly seeded around the curve fitting results in Sect. 3.3.2 with a standard deviation of 0.01, which can greatly reduce the length of the chain as a result of the fast convergence of the MCMC fitting. We take the 50th per-centile and half of the difference between the 16th and the 84th percentiles of the posterior distribution sampled by the MCMC as our best estimate and statistical uncertainty. The equivalent width (EW) of the DIB is calculated by $\text{EW}=\sqrt{2\pi }\times \text{D}\times \sigma$. The EW error is estimated by the error propagation from the errors of D and σ.

4 Quality control and validation

4.1 Quality control

4.1.1 Defining the high-quality samples

A total of 2 188 240 target spectra are processed by the pipeline described in Sect. 3, and the statistics of all the valid measurements for the DIBs λ5780, λ5797, and λ6614 are summarized in Table 3. To analyze the results based on the reliable DIB measurements, we define the high-quality (HQ) samples. For each HQ sample, it must meet the following criteria:

1. The depth fitted by the MCMC ${\text{D}}_{\text{mc}}>\frac{3}{ \text{S}/\text{N}}$. The threshold of $\frac{3}{ \text{S}/\text{N}}$ is set to ensure the DIB signal is significantly above the noise level.

2. The coefficient of variation (CV) of the depths {D_pr, D_cf, D_mc} < 10%, where the CV is defined as the standard deviation divided by the mean value, and the subscripts pr, cf and mc represent the pre-detection, curve fitting and MCMC fitting in Sect. 3, respectively.

3. The standard deviation of the central wavelengths {µ_pr, µ_cf, µ_mc} < 0.5 Å.

4. |µ_mc − λ₀| < 3 Å, where λ₀ is the typical central wavelength of DIB (Herbig 1975; Fan et al. 2019; Vogrinčič et al. 2023), that is, λ₀ = 5780.6, 5797.1, 6613.6 Å for the DIBs λ5780, λ5797, and λ6614, respectively. This criterion is mainly used to confine the RV of the DIB carrier to within around ±200 km s⁻¹, which is a reasonable assumption if the DIB carrier primarily tracks the ISM at a distance of several kiloparsecs from the Sun (Zhao et al. 2021a).

5. The Gaussian width σ_mc > 1 Å and σ_mc < 3 Å, where the lower limit of 1 Å is set due to the instrumental broadening effect (see Appendix A), and the upper limit of 3 Å is set based on the typical FWHM of DIB (Vogrinčič et al. 2023), which will be elaborated in Sect. 4.2.1.

6. EW_DIB/EW_cont > 33% for the DIB λ5780 and >25% for the DIBs λ5797 and λ6614, where EW_DIB is the equivalent width of the DIB, and EW_cont is the integral equivalent width of the continuum outside the DIB region, namely [µ_mc − 3σ_mc, µ_mc + 3σ_mc]. The constraints can ensure the DIBs rather than the stellar residuals or noises significantly contribute to the EWs.

Besides the above criteria, we also tighten the constraints on the CV and central wavelength for the M- and K-type (T_eff < 5200 K) stars to get rid of the fake DIB signals caused by the stronger metal absorptions in cooler stars. Specifically, the CV must be less than 5% and the standard deviation of {µ_pr, µ_cf, µ_mc} must be less than 0.25 Å. After applying the above criteria, we obtain the HQ samples, and the summary of the HQ samples are also listed in Table 3. There are a total of 7681 spectra simultaneously contain the HQ measurements of the DIBs λ5780, λ5797, and λ6614. Besides the public access to three catalogs of the HQ samples for the DIBs λ5780, λ5797, and λ6614, here we also release the results of these 7681 HQ measurements⁶ via https://nadc.china-vo.org/res/r101404/.

4.1.2 Defining the marginal-quality samples

Aside from the HQ samples, the marginal-quality (MQ) samples are also available online. Although the MQ samples cannot provide DIB measurements as reliable as the HQ samples, these low or null values can offer important boundary constraints for other studies, for example, 3D mapping. For each MQ DIB measurement, it is required to meet the following two criteria.

1. The EW fitted by the MCMC, EW_DIB ≤ 0.03 Å, where the upper limit is set to 0.03 Å according to the typical uncertainty of 0.04 Å for the DIB measurements validated in Sect. 4.2.2, and the lowest EW of around 0.03 Å measured in the HQ samples.

2. EW_cont ≤ 0.5 Å, where EW_cont is the integral equivalent width of the continuum outside the DIB region, namely [µ_mc − 3σ_mc, µ_mc + 3σ_mc]. The upper limit is set to 0.5 Å to ensure that the stellar residuals are cleaned up to the continuum level.

After applying the above criteria, there are a total of 27 598, 42735, and 21635 MQ measurements of the DIBs λ5780, λ5797, and λ6614, respectively, corresponding to 23 920, 36 825, and 18 145 unique stars. We note that we do not set a specific threshold to further distinguish between marginal and null values, as the uncertainties estimated by MCMC are provided in the catalogs. Therefore, users can determine the threshold based on their needs. However, we suggest using a ratio of err(EW_DIB)/EW_DIB of 1 as the dividing line: values below 1 can be considered marginal, while values above 1 can be classified as null.

4.2 Validation

4.2.1 Central wavelength versus Gaussian width

The λ–σ 2D histogram has been a powerful diagnostic tool for visually distinguishing real DIBs from false positives in the context of automatic measurement of DIBs from a large spectroscopic survey (Saydjari et al. 2023; Zhao et al. 2024). The false-positive DIBs, which may be stellar absorption, are usually far away from the typical central wavelength of the DIBs, and their Gaussian widths are narrower than the typical values of the real DIBs as the FWHM of a DIB is generally broader than that of metal absorption line. In the 2D histogram, the true DIB measurements are located around the rest-frame central wavelength and typical Gaussian width of the DIB, but with a certain scatter, whereas the stellar features exhibit smaller Gaussian widths and are highly concentrated in a region.

Figure 5 illustrates the λ–σ distribution of the DIBs λ5780, λ5797, and λ6614 for the entire samples and the HQ samples. It is clear that for either the entire samples or the HQ samples, the DIBs λ5780 and λ6614 are mainly distributed nearby the crossed red dashed lines that represent the statistical central wavelength fitted by the method in Sect. 5.3 and the broadened Gaussian width (see Appendix A) of the DIB measured by Vogrinčič et al. (2023). However, the distribution of the DIB λ5797 for the entire samples displays an abnormal high-density region centered at about (σ, λ) = (0.9, 5797.8) Å, which deviates from the expected position marked by the red dashed lines. Herbig (1975) has reported a Si I line at 5797.859 Å can obliterate the DIB λ5797 in the spectra of the type F and later stars, which explains such an abnormal high-density region. Nevertheless, after the strict quality control, particularly with the criterion of 1 < σ_mc < 3 Å, this anomalous area is almost completely eliminated, and the distribution of the DIB λ5797 for the HQ samples is consistent with the expected position.

Fig. 5 Two-dimensional histograms of the measured DIBs λ5780, λ5797, and λ6614 as a function of the Gaussian width σ and the central wavelength λ in the rest frame. (a) The three upper panels present the density of all the samples that may have DIBs, (b) and the other three lower panels show the density of the HQ samples for which we have reliable DIB measurements. All of the histograms have the same bin size of 100. The red horizontal dashed lines denote the central wavelengths of the DIBs, which are fitted using the method described in Sect. 5.3. The red vertical dashed lines indicate the typical Gaussian widths of the DIBs considering the instrumental broadening of LAMOST LRS (see Appendix A), which are taken from Vogrinčič et al. (2023).

4.2.2 Integral EW versus fitted EW

To quantify the reliability of the DIB measurements, we compute the integral EW (EW_int) of the three DIBs in the HQ samples within the same wavelength range, that is, µ_mc ± 3σ_mc as the fitted EW (EW_fit). As shown in Fig. 6, the EW_int is highly consistent with the EW_fit, and the difference between them (∆EW) is mainly concentrated in a tiny region that is less than 0.04 Å for over 99% of the samples among three DIBs. For all three DIBs, however, there are certain samples whose EWint are far away from their EW_fit, which stem from the fact that the continuum level is not well determined (Farhang et al. 2015b; Farhang et al. 2019). For the DIB λ5780, the blend with a very broad DIB λ5778 is a confounding factor in the determination of the continuum level (Herbig 1975; Krelowski & Walker 1987; Herbig 1993). Similar with the DIB λ5780, the blue wing of the DIB λ5797 is also overlapped with a quite broad DIB λ5795 (Krelowski & Walker 1987; Friedman et al. 2011). Nevertheless, Vos et al. (2011) has reported that the broader DIBs λ5778 and λ5795 are largely eliminated from the ISM residual spectra in the local normalization, and they are too weak to significantly contaminate the EW measurements of the DIBs λ5780 and λ5797, even for high extinction. Such broader DIB blend does not exist for the DIB λ6614 (Herbig & Soderblom 1982; Galazutdinov et al. 2008), which gives rise to better agreement (see Fig. 6(c)) between the EW_int and the EW_fit than the other two DIBs.

A sequential eye check for certain ISM residual spectra with large ∆EW reveals that the blend with the broader DIB does destroy the continuum level and makes the fitted Gaussian width wider than the specified DIB. Assuming the continuum level is well defined, a multi-Gaussian profile may be a better choice to fit the DIB with a DIB blend. For example, Lan et al. (2015) has adopted a double Gaussian profile to simultaneously fit the DIB λ5780 and λ5778. In addition, the extra structural features are found among those ISM residual spectra, which are likely to be the stellar residuals (Zhao et al. 2024). One probable reason is that the stellar features are not completely erased by the subtraction of the stellar template from the target spectrum, which can be further attributed to the mismatch between the stellar template and the target spectrum, or the incorrect alignment caused by the false radial velocity. Nonetheless, the overall agreement between the EW_int and the EW_fit can support the reliability of the DIB measurements and reinforce the following analysis in statistical properties.

Fig. 6 Comparison between the integral EW and the fitted EW of the DIBs λ15780, λ5797, and λ16614 for their HQ samples, from the subfigures (a) to (c). For each subfigure, the upper panel shows the Gaussian KDE number density plot of the integral EW as a function of the fitted EW, with the red dashed line representing the one-to-one relation. The lower panel displays the Gaussian KDE number density plot of the difference between the integral EW and the fitted EW, namely ∆EW = EWint − EWfit as a function of the fitted EW, with the green dashed line representing the zero benchmark. The distribution of ∆EW is also overplotted in the upper panel with its mean and standard deviation annotated.

4.2.3 Multiple epochs

Benefiting from the high efficiency in the spectrum acquisition, LAMOST has accumulated numerous sources with multiple epochs. These sources that have been observed many times can provide a good opportunity to validate the variation of DIB measurement. We select the sources with over two epochs from the HQ samples by grouping the sources with the same uid (unique source identifier in LAMOST) and present the distribution of the number of epochs in the panels in the first column of Fig. 7. There are 15 839 891, and 9939 sources with over two epochs for the DIBs λ5780, λ5797, and λ16614, respectively. These sources have four epochs on average and the maximum number of epochs is up to 33. The panels in the second column of Fig. 7 illustrate the cumulative distribution of the standard deviation of the EWs (σ_EW) for multiple-epoch sources. Almost 95% of the sources have σ_EW < 0.05 Å for three DIBs, which is comparable to the typical value, namely 0.04 Å of ∆EW investigated in Sect. 4.2.2. Such similarity between σ_EW and ∆EW implies that the variation of the DIB measurements among multiple epochs of the same source is resulted from the measurement error rather than the intrinsic property of the DIB. Without the influence of blend with a broader DIB on the DIB λ6614 (see Sect. 4.2.2), it is not surprising that the cumulative distribution of the DIB λ6614 accumulates toward 1 more rapidly than that of the DIBs λ5780 and λ5797.

In order to present the relative variation of the EWs, we also supply the cumulative distribution of the coefficient of variation (CV = σ_Ew/EW_mean) in the panels in the third column of Fig. 7. Well over three quarters of the sources have CV < 0.2 for the DIB λ5797, and this proportion is more than 90% for the DIBs λ5780 and λ6614. The large CV of the multiple-epoch sources generally comes from the small EW, which appears to be more sensitive to the measurement error, as can be seen in the panels in the fourth column of Fig. 7. It is noted that both the cumulative distribution of the σ_EW and the CV demonstrate the long tail. Besides the measurement error mentioned above, such an imbalance is driven by the differences among the atmospheric parameters of multiple spectra for the same source. Specifically, although these spectra originate from the same source, the differences in observation conditions, processing procedures, instrument running status, etc., result in variations in atmospheric parameter measurements among these spectra. Consequently, their matched closest neighboring template spectra also differ, which ultimately leads to a noticeable discrepancy in the ISM residual spectra. In the visual inspection, we found that the maximum difference in T_eff among multiple-epoch spectra of the same source can reach 500 K, far exceeding the match threshold of 100 K set in Sect. 3.2.1. Zhao et al. (2024) and Saydjari et al. (2023) have circumvented the issues of atmospheric parameters and directly extracted the DIB λ8621 signal from the parameter-free spectra using machine learning methods. The applicability of this method, however, requires further validation, which is also one of our future directions.

4.2.4 Residuals

In addition to validating the measurements of DIBs, the quality of ISM residual spectra is also verified. Ideally, after subtracting the stellar features using the template spectrum, the flux of the ISM residual spectrum outside the DIB region which is from µ_mc − 3σ_mc to µ_mc + 3σ_mc should be at the continuum level with no stellar components. Thus, we calculate the standard deviation of the residuals between the flux of ISM residual spectrum outside the DIB region and that of the continuum, and then provide the trend of STD with S/N and T_eff in Fig. 8. As can be seen in the figure, the STD decreases as S/N increases, which indicates that higher S/N leads to smaller residuals in the ISM residual spectrum or cleaner subtraction of stellar components. The STD also decreases with increasing T_eff as a consequence of the weakening of metal lines in hotter stars. However, this decreasing trend of T_eff is more gradual compared to that of S/N, reflecting the robustness of measurement, which is less disturbed by absorption lines in cool stars. Furthermore, the magnitude of STD is comparable to the level of normalized flux error estimated in Sect. 3.2.2, which indicates that the pipeline does not introduce additional errors and that the ISM residual spectrum is quite clean, with no significant residual stellar components. It is noteworthy that although there probably exist the DIB components other than the interesting ones within the range of the ISM spectrum, they are usually weak absorptions, such as the DIBs λ5762 and λ5815 in Fig. 3, documented in the previous DIB catalogs (Fan et al. 2019; Vogrinčič et al. 2023). Moreover, after instrumental broadening by the LAMOST LRS, the depths of those weak DIBs and the flux errors share the close values, which cannot cause any significant bias in the final results.

Fig. 7 Three distinct colors encode the three DIBs λ5780, λ5797, and λ6614 for the HQ samples with over two epochs. The panels in the first column present the distribution of the number of epochs. The panels in the second column depict the cumulative distribution of the standard deviation of the EWs (σEW). The panels in the third column display the cumulative distribution of the coefficient of variation (σEW/EWmean). The panels in the fourth column show the scatter of the σEW/EWmean as a function of the EWmean.

5 Results and discussion

In this section, we review and explore the spatial distribution, correlation, central wavelength and kinematics of the DIBs λ5780, λ5797, and λ6614 in the context of statistics based on the HQ samples. Due to the multiple-epoch sources mentioned in Sect. 4.2.3, we combine the measurements of DIB from the same source and take the mean value as its final result. After combination, the number of the HQ sources for the DIBs λ5780, λ5797, and λ6614 is 142 074, 11 480, and 85 301 respectively.

5.1 Spatial distribution

Figure 9 presents the number density of the background stars and the sky distribution of the median EW of the DIBs λ5780 and λ6614. The absence of the DIB λ5797 in the figure is due to the scarcity of the HQ sources caused by the low detection rate (see Table 3). In order to present a more generalized sky distribution of DIB in the lower panels of Fig. 9, we remove the pixels in which the number density of the background stars is less than 3. Among the removed pixels, some exhibit significantly larger EWs than the values of DIB in the surrounding areas. Upon inspecting the spectra of the DIBs with large EW, we find that their background stars are located at considerable distances, which can lead to the large extinctions and EWs. However, there are some large EWs of DIB at the sight line of the short-distance background stars as well, which are likely to be the fake DIB signals investigated in Sect. 4 or the incomplete correlation between the DIB and extinction discussed in Sect. 5.2.1.

The DIBs λ5780 and λ6614 share a similar sky distribution, that is, the Galactic plane is the prominent region with the strong DIBs and the strength of the DIBs decreases with the increasing Galactic latitude. The lower number of samples for DIB λ6614 compared to DIB λ5780 at high latitudes is caused by the weaker intensity of DIB λ6614 per unit extinction (Vogrinčič et al. 2023), making it less detectable in the regions of low extinction at high galactic latitudes. The high-EW regions of the DIBs λ5780 and λ6614, under the Galactic view, are both spatially correlated with the well-known molecular clouds, such as Taurus and Orion. This behavior is similar to that of the DIB λ8621 investigated in Gaia RVS spectra (Gaia Collaboration 2023b; Saydjari et al. 2023). In addition, there is a notable extent to high latitudes (b ≈ −45°) in the directions of Galactic anti-center, which does not entirely align with the distribution of molecular clouds. This stretch from the plane to high latitudes is similar to that of the DIB λ8621 depicted by Gaia Collaboration (2023b), Saydjari et al. (2023), and Zhao et al. (2024). The causes of such a phenomenon remain unclear, but it is likely attributable to the selection effect of the observation, as can be seen in the sky distribution of the number density in the upper panels of Fig. 9. Another possible explanation is that the resolution of map is not high enough to resolve the sky distribution of the DIBs, which makes it appear as though the region is interconnected. It is anticipated that with the inclusion of DIB measurements from hot stars of LAMOST in the future, more detailed and accurate northern sky distribution of the DIBs will be obtained. Additionally, a three-dimensional spatial distribution of the DIBs, such as Kos et al. (2014) and Farhang et al. (2019), will be constructed, which can provide a more comprehensive understanding of the DIBs and their positional and physical relationship with the ISM.

Fig. 8 Trend of the standard deviation of the residuals between the flux of the ISM residual spectrum beyond the DIB region and that of the continuum as a function of S/N in the upper panels, and Teff in the lower panels. Three distinct colors encode the three DIBs λ5780, λ5797, and λ6614. In each panel, a violin plot is used to present the distribution of the standard deviation of the residuals for the binned HQ samples. The median of each bin is annotated with a white dot, and the 2.5 and 97.5 percentiles are marked with black dots. The red dashed line connects the median of each bin to show the trend of the standard deviation of the residuals.

5.2 Correlation

There have been numerous studies on the correlation of DIB with DIB and DIB with other ISM tracers, such as extinction, CO gas, atomic hydrogen, and molecular hydrogen (e.g., Herbig 1975; Weselak et al. 2004; Friedman et al. 2011; Lan et al. 2015). However, these works are more or less limited by the number of sight lines or combined the spectra along the same sight line for high S/N, which sacrifice the statistical generalization. For the sake of statistical and homogeneous comparison, a total of 7681 common samples for the DIBs λ5780, λ5797, and λ6614 from the HQ samples are picked out to investigate the correlation of DIB V.S. extinction and DIB V.S. DIB.

5.2.1 Correlation between DIBs and E(B–V)

To avoid zero or inaccurate values of E(B − V) at close distances, the heliocentric distance d is limited to greater than 500 pc and a total of 6514 common HQ sources are selected. Figure 10 compares the EW of DIB with the extinction E(B − V) derived from Green et al. (2019) by the linear fit to the median EW in each ∆E(B − V) = 0.05 mag bin. All three DIBs agree with a linear positive correlation with E(B − V) in Table 4, but not perfect, as known for a few decades (Herbig 1995). Nonetheless, the areas with high number density (in red) show smaller scatters and have better linear relationship (red solid lines). Even in sparsely sampled regions, such as those near low extinction of 0.15 mag and high extinction of over 0.6 mag, the medians do not deviate from the fitted curve. We do not include the samples with the higher extinction (>1 mag) in the fitting process due to the well-known “skin effect” (Herbig 1995), namely the strength of DIB cannot indefinitely increase with the extinction as expected. Besides, we also attempt a linear fit forced through the origin (green dashed lines) to facilitate the comparison with the previous works, for instance, Herbig (1975). It is clear that the green dashed lines are not suitable to represent the correlation of the DIBs with E(B − V), as they have larger fitting errors and deviate from the red solid lines.

As one of the strongest and earliest discovered DIBs, the correlation of DIB λ5780 with E(B − V) is the most studied. Our slope of 0.606 Å mag⁻¹ is close to the slope of 0.64 Å mag⁻¹ reported by Vos et al. (2011)⁷. The slopes of 0.549 Å mag⁻¹ computed by Herbig (1975)⁸ , 0.505 Å mag⁻¹ fitted by Friedman et al. (2011)⁹, 0.525 Å mag⁻¹ reported in Raimond et al. (2012), and 0.43 Å mag⁻¹ derived by Lan et al. (2015)¹⁰ are all smaller than our slope. Some other works have even obtained the slopes lower than those found in the above works, such as the slope of 0.401 Å mag⁻¹ from Puspitarini et al. (2013) and the slope of 0.419 Å mag⁻¹ from Vos et al. (2011)¹¹. There are three aspects that can cause the discrepancy of these slopes. One is the measurement error of the DIBs and the extinction. Apart from Lan et al. (2015), the other studies focused on their data within an absolutely narrow range of EW or E(B − V). For instance, the mean EW of the DIB λ5780 in Puspitarini et al. (2013) is about 60 mÅ, with the measurement errors exceeding the EWs themselves (see Fig. 6 in their work). Similarly, the samples with large EW in Friedman et al. (2011) hardly influence the fitting and significantly deviate from the fitted curve. Consequently, the fitting results are predominantly driven by the samples with E(B − V) < 0.2 mag (refer to the enlarged region in Fig. 4 of their work). Another reason is the different profile definition of the DIB λ5780 (see Sect. 4.2.1 for more details). Lan et al. (2015) separated the broader DIB λ5778 from the DIB λ5780 and simultaneously fitted them with a double Gaussian profile, whereas the other works, including ours, treated them as a single DIB. Despite the discovery of the DIB λ5778 as early as 1975 (Herbig 1975), the subsequent research on the λ5780 often overlooked this blended broad DIB. If the contribution of the DIB λ5778 by Lan et al. (2015) is attributed to the DIB λ5780, the co-added result will align with the slope obtained in this work. The last key factor is the selection of the sight lines. Vos et al. (2011) elaborated the influence of σ- and ζ-type sight lines on DIB carrier according to the stars from the upper Scorpius. They concluded that the DIB λ5780 in the σ-type sight line is stronger than that in the ζ-type (Krelowski & Westerlund 1988; Kos & Zwitter 2013), which may be caused by the different ultraviolet (UV) radiation fields, that is, the external ionized regions for the σ-type and the UV-shielded could cores for the ζ-type. It is also worth noting that, compared to the DIBs λ5797 and λ6614, the dispersion of EWs in each extinction bin for the DIB λ5780 increases with E(B − V), which further validates the σ − ζ effect on the DIB λ5780. Such a σ − ζ effect on our slope is beyond the scope of this work, but our correlations of the DIB λ5780 with E(B − V) offer more generalizable conclusions, and circumvent the sample selection effect (the limited number of sight lines) and the measurements that are sensitive to their own errors. Furthermore, the large number of samples in this work is likely to encompass more σ- and ζ-type sight lines, which will help us to understand the DIB λ5780 along the different sight lines in a more comprehensive way.

For the DIBs λ5797 and λ6614, their correlations with E(B − V) resembles those of previous works. Our slope of 0.2 Å mag⁻¹ for the DIB λ5797 is consistent with 0.219 Å mag⁻¹ derived by Herbig (1975), 0.174 Å mag⁻¹ computed by Friedman et al. (2011), and 0.18 Å mag⁻¹ fitted by Lan et al. (2015). By similar token, the slopes of 0.262 Å mag⁻¹ from Herbig (1975), 0.216 Å mag⁻¹ reported in Friedman et al. (2011) and 0.22 Å mag⁻¹ documented by Lan et al. (2015) for the DIB λ6614 are in good agreement with our slope of 0.274 Åmag⁻¹. Such consistency of the slopes among different works are unsurprising, as Kos & Zwitter (2013) indicated that the DIBs λ5797 and λ6614 are unaffected by the σ − ζ effect. Additionally, although the DIB λ5797 is influenced by the nearby DIB λ5795 (Krelowski & Walker 1987; Friedman et al. 2011), this impact is confined to the blue wing of the DIB λ5797 and not as significant as the complete overlap of the DIB λ5778 with the broad DIB λ5780. As to the DIB λ6614, it is not blended with any other DIB (Herbig & Soderblom 1982; Galazutdinov et al. 2008), which makes the correlation with E(B − V) more straightforward. The probable cause of the slight difference in slope still stems from the measurement error of the DIBs and the extinction considering weaker EW per unit extinction of the DIBs λ5797 and λ6614 compared to the DIB λ5780.

We also take the same data used to linear fittings to calculate the Pearson correlation coefficient between EW and E(B − V) for the three DIBs, which are 0.54, 0.51, and 0.64 for the DIBs λ5780, λ5797, and λ6614, respectively. None of these DIBs reaches a significant level of correlation with extinction, consistent with the conclusions of many recent studies (Friedman et al. 2011; Fan et al. 2017), which is caused by the variable interstellar environments as mentioned before. Furthermore, the differences of the number of data points that are used to compute correlation coefficients can also lead to the discrepancy of the interpretation of correlation. To generalize our conclusions further, besides the 6514 common HQ sources, we perform the same fitting procedure on all HQ sources at heliocentric distances greater than 500 pc. The results are also summarized in Table 4, which are comparable to those of the common HQ sources.

Fig. 9 Galactic distribution of the DIBs λ5780 and λ6614 for the HQ samples. The upper panels show the number density of the background stars. The lower panels display the sky distribution of the strength of the DIBs with the color bar representing the median EW in each pixel. The black contours of the corresponding median E(B − V) at the level of 0.35 and 0.7 magnitude are overplotted in the median EW map. The resolution of the map is about 0.92° (Nside = 64) calculated by the HEALPix package.

Fig. 10 Gaussian KDE number density plot of the EWs of the DIBs λ5780, λ5797, and λ6614 (from top to bottom) as a function of the extinction E(B − V) derived from Green et al. (2019). For each panel, the red solid line represents a linear fit to the median EW (which is marked with a white dot, and its standard deviation marked is with a black error bar) in each ΔE(B − V) = 0.05 mag bin. The green dashed line is the same as the red solid line, but the intercept is forced to be zero. The coefficients of the linear fit and the Pearson correlation coefficient between the EW and E(B − V) are also annotated in each panel.

5.2.2 Correlations between diffuse interstellar bands

The correlation between DIBs can be used to infer whether their carriers come from the same DIB family or even share a common origin (Moutou et al. 1999; Friedman et al. 2011; Smith et al. 2021). Figure 11 illustrates the mutual relation of DIB pairs and their Pearson correlation coefficients.

DIB λ5780−λ5797: This pair of DIBs has been widely studied in the literature since the first discovery of the DIBs λ5780 and λ5797 (Heger 1922). The general ratio of EW₅₇₉₇/EW₅₇₈₀ is 0.263, which is lower than the value of 0.384 reported by Friedman et al. (2011). The high scatter of the data points used for fitting in Friedman et al. (2011) is the main cause of this discrepancy. In Fig. 10 of Friedman et al. (2011), namely the figure showing the fitting results of EW₅₇₉₇ and EW₅₇₈₀ , the fitting is predominantly driven by their samples with EW₅₇₈₀ of less than 140 mÅ (see the inset of their Fig. 10), but approximately one-third of the data points lie below the fitted line. In particular, there is a strong deviation trend for the samples with EW₅₇₈₀ between 150 and 300 mÅ. Additionally, our fitting is performed on the medians within each bin, rather than on all data points together as done by Friedman et al. (2011). Results of the later approach is often more susceptible to outliers. The famous σ − ζ effect (see the discussion in Sect. 5.2.1) can be a probable cause of such a high scatter in Friedman et al. (2011), as seen, for instance, in the sight line towards Herschel 36 and 6 Cas annotated in their Fig. 10. It is also worth noting that while Friedman et al. (2011) used the OB stars within the Galactic plane, our sample includes the sources at higher Galactic latitudes, which can also contribute to the statistical differences between the two studies. The ratio of EW₅₇₉₇ /EW₅₇₈₀, however, are more complex than hitherto thought, and cannot be adequately represented by a simple linear relationship. In addition to the σ − ζ effect, Farhang et al. (2019) proposed that the carrier of the DIB λ5780 can resist X-ray photodissociation and sputtering by fast ions, whereas the carrier of the DIB λ5797 breaks down. A more sophisticated method, such as a 3D projection map (Alves et al. 2020; Saydjari et al. 2023), is anticipated to depict it accurately.

DIB λ5780−λ6614: Among all three DIB pairs, this pair shows the most significant correlation ranked by the Pearson correlation coefficient of 0.788. Interestingly, and somewhat paradoxically, many studies have provided conflicting conclusions regarding their classification into specific families or groups according to the different perspectives. For example, based on their intensity behavior through different types of clouds (σ- and ζ-types), Kos & Zwitter (2013) contended that the DIB λ6614 belongs to type I (where the linear relations with extinction are similar) and the DIB λ5780 is type II (where they changed substantially). In contrast, Fan et al. (2022) categorized the DIBs λ5780 and λ6614 into the same σ-DIB group using a hierarchical agglomerative clustering method based on EW of DIB normalized to extinction. Additionally, in the studies conducted by Moutou et al. (1999), Friedman et al. (2011), and Smith et al. (2021), including ours, the correlation of this pair ranked among the top across all DIB pairs. One reasonable explanation for the conflict is that the carriers of the DIBs λ5780 and λ6614 display different properties in various physical environments, as Fan et al. (2022) argued that the DIBs form a rather continuous sequence. It is worth emulating a multidimensional approach by Lan et al. (2015) to decouple different effects on DIB. Through combining extra information such as the sky distribution and kinematics with the ISM tracers, a more exhaustive understanding of DIB will be achieved.

DIB λ6614−λ5797: There is no evidence that supports the close relationship between the DIBs λ6614 and λ5797 (Moutou et al. 1999; Friedman et al. 2011; Smith et al. 2021), as the Pearson correlation coefficient is 0.673, ranking the lowest among the three DIB pairs. Although the DIBs λ5780 and λ5797 unaffected by the σ − ζ effect (Vos et al. 2011), Fan et al. (2022) suggested that they belong to distinct DIB families.

Fig. 11 As in Fig. 10, but for the Gaussian KDE number density plot of EW V.S. EW. From the top panel to the bottom panel, the DIBs λ5780, λ5797, and λ6614 are compared with each other. For each panel, the red solid line represents a linear fit to the median EW in each ΔEW = 0.05 Å bin. Except for the bottom panel, the bin size of the EW is ΔEW = 0.02 Å. The coefficients of the linear fit and the Pearson correlation coefficient between the EWs are also annotated in each panel.

5.3 Central wavelength in the rest frame

The central wavelength of DIB in the rest frame¹² is one of the most important properties, which can identify the DIB carrier through the match of the band wavelengths and the strength ratios between the observational and laboratory measurements (Campbell et al. 2015; Campbell & Maier 2018; Lallement et al. 2018; Cordiner et al. 2019). Table 5 provides the central wavelength of the DIBs λ5780, λ5797, and λ6614 from five canonical DIB catalogs, which are measured using the Doppler shift of the well-identified interstellar atomic or molecular lines, such as Na I D doublet, K I line, and CH line (see Table 2 in Jenniskens & Desert 1994). This traditional method is based on the strong assumption that the DIB carriers co-move with the atoms or molecules so that they share the same Doppler shift. However, all the five catalogs focused on a single or a dozen of sight lines, and there are probably multiple ISM components along the sight line, which can lead to the confusion of the ISM Doppler shift (Herbig & Soderblom 1982). Another method to measure the central wavelength of DIB is based on an empirical assumption that, in the Local Standard of Rest (LSR), the matter moving towards the Galactic center (GC) or the Galactic anticenter (GAC) tends to be at rest, meaning the statistical radial velocity of the matter approaches zero (Bovy et al. 2012). Without the interstellar reference lines, Zasowski et al. (2015) were the first to successfully measure the central wavelength of the DIB at 1.5273 µm using this statistical method, so did Zhao et al. (2021a), Gaia Collaboration (2023b), and Saydjari et al. (2023) for the DIB λ8621.

Although the spectral wavelength range of LAMOST LRS includes the interstellar absorption lines such as the Ca ii K line at 3933.66 Å and Na i D doublet at 5889.95 and 5895.92 Å, these lines are prone to blending with stellar absorption lines at the same wavelengths due to the cool stars. Additionally, owing to the increasing modernization around the observation site, the light emitted by high-pressure sodium lamps used for street lighting has contaminated the Na I D doublet in the spectra (private communication, Chao Liu, 2023). Poznanski et al. (2011) also concluded that Na i D absorption at low resolution is a bad proxy for extinction. As a consequence, the statistical method instead of the ISM Doppler shift method is employed in this work, but we plan to compare the differences in the central wavelength of DIB between the two methods when the measurements of DIB from hot stars are available in the future. The samples with |b| < 2°, 170° < ℓ < 190°, d ≤ 4 kpc, err(µ_mc) < 0.5 Å, and err(RV) < 5 km s⁻¹ are selected to fit the central wavelength of the DIBs. For the DIB λ5797, we loosen the constraint of |b| < 2° to |b| < 10° because of its limited number of the HQ sources. There are 4934, 1542, and 3186 HQ sources in the GAC for the DIBs λ5780, λ5797, and λ6614, respectively. Figure 12 shows the observed central wavelength of DIB as a function of the angular distance from the GAC. By the linear fit to the medians of each bin (∆ℓ = 1°) using the least squares method, the observed central wavelengths of the DIBs λ5780, λ5797, and λ6614 at ℓ = 180° are 5780.68 ± 0.01, 5797.14 ± 0.02, and 6613.87 ± 0.01 Å, respectively. We further take the effect of the solar motion into account, and the final central wavelength λ0 can be expressed as λ_obs × c/(c + U_⊙), where λ_obs is the observed central wavelength, c is the speed of light, and U_⊙ = 10.6 km s⁻¹ is the radial solar motion (Reid et al. 2019). As listed in Table 5, the final central wavelengths λ₀ of the DIBs λ5780, λ5797, and λ6614 are 5780.48 ± 0.01, 5796.94 ± 0.02, and 6613.64 ± 0.01 Å, respectively, which are all consistent with the central wavelengths from the cross-reference. Nonetheless, the central wavelength of the DIB λ5797 is slightly lower than those from the cross-reference. The obtained slopes are also kind of smaller than the value of 23 ± 3.4 mÅ deg⁻¹ derived by Gaia Collaboration (2023b) and the value of 19 mÅ deg⁻¹ estimated by Saydjari et al. (2023) for the DIB λ8621 in the GAC, but are quite lower than the value of 57 ± 8 mÅ deg⁻¹ reported by Zasowski et al. (2015) for the DIB at 1.5273 µm.

Fig. 12 Observed central wavelength (λobs) of the DIBs λ5780, λ5797, and λ6614 as a function of the angular distance (∆ℓ) from the GAC, from top to bottom. The gray dots and their error bars represent the measurements of the central wavelength of DIB and the corresponding uncertainties. The white dots and their black error bars denote the medians and the standard deviations in each ∆ℓ = 1° bin. The red lines are the linear fits to the medians of each bin using the least squares method. The fitting results, that is, the slopes k, the intercepts b, and the final central wavelengths λ0 considering the solar motion, are annotated in each panel.

5.4 Kinematics

Despite the unknown DIB carriers, many works have attempted to use DIB as a tracer to study the kinematics of the ISM since the DIB carriers are proven to be interstellar. Zasowski et al. (2015) provided what 3D distribution the carriers of the DIB at 1.5273 µm have, and were the first to elaborately depict how these DIB carriers move with the Milky Way by the mean Galactic velocity curve. Following this work, Zhao et al. (2021b) and Gaia Collaboration (2023b) displayed the kinematics of the DIB λ8621 in the same way, while the average galactocentric distance of the DIB λ8621 estimated by the velocity curve is about 7.5 kpc, smaller than 9 kpc derived from Zasowski et al. (2015) due to the differences of the distances of the observed stars. In addition, a comparison of the LSR velocity of the DIB 48621 with that of CO gas (Dame et al. 2001) shows they own the same kinematics pattern, which was quantified by the linear correlation of a slope of 0.95 using the algorithm called “peak-finding” (Saydjari et al. 2023). Inspired by these works, we investigate the LSR velocity of the DIBs 45780 and 46614, $V_{\text{DIB}}^{\text{LSR}}$, as a function of the Galactic longitude ℓ, as presented in Fig. 13. The comparison of $V_{\text{DIB}}^{\text{LSR}}$ with the LSR velocity of the CO gas and the other tracers of the ISM will be explored in the forthcoming work.

We select the HQ sources with |b| < 12°, err(µ_mc) < 0.5 Å, and err(RV) < 5 km s⁻¹ to study the kinematics of DIB and there are a total of 71 061 and 47 011 HQ sources for the DIBs 45780 and 46614. To convert the observed central wavelength of DIB to the LSR velocity, we take the effect of the solar motion V_⊙ = (U_⊙, V_⊙, W_⊙) = (10.6, 10.7, 7.6) km s⁻¹ (Reid et al. 2019) relative to the LSR into account, and $V_{\text{DIB}}^{\text{LSR}}$ can be expressed as Eq. (7). $$V_{\text{DIB}}^{\text{LSR}}=c\times \left({\lambda }_{\text{obs}}/{\lambda }_0-1\right)+V_{\odot }/\widehat{r},$$(7)

where c is the speed of light, λ_obs is the observed central wavelength of DIB, λ₀ is the final central wavelength of DIB fitted by the statistical method in Sect. 5.3, and $\widehat{r}$ is the unit direction vector from the Sun to the source. The Galactic velocity curves with diverse heliocentric distances d are computed using the Model A5 of Reid et al. (2019). It can be seen from Fig. 13 that the fluctuation of $V_{\text{DIB}}^{\text{LSR}}$ for the whole selected samples reflects the motion of the DIB carrier with the Galactic rotation. The medians of $V_{\text{DIB}}^{\text{LSR}}$ in each bin (∆ℓ = 10°) are both consistent with the velocity curve with the heliocentric distance d = 0.5 kpc at the 30° < ℓ < 180° and the d = 1.5 kpc at the ℓ > 180°. The results agree with the distance of the DIB 48621 estimated Gaia Collaboration (2023b), but may be related to the fact that most of the background stars are located within 3 kpc.

Fig. 13 Tow-dimensional histogram of the LSR velocity of the DIBs 45780 (upper panel) and 46614 (lower panel) as a function of the Galactic longitude. The medians and the standard deviations in each ∆ℓ = 10° bin are denoted by the white dots and their black error bars. The Galactic velocity curves with different heliocentric distances d calculated by the Model A5 of Reid et al. (2019) are overlaid in the plots.

6 Summary and conclusions

We successfully extracted 2 188 240 ISM residual spectra from the late-type stellar spectra in the LAMOST LRS using the method of subtracting the stellar templates. By fitting the ISM residual spectrum with the Gaussian profile, we obtained the largest sample of reliable measurements of the DIBs 45780, 45797, and 46614 in the northern sky to date. The selected 176 831, 13 473, and 110152 HQ samples of the DIBs 45780, 45797, and 46614, respectively, corresponding to 142 074, 11 480, and 85 301 unique sources, not only allow for a more comprehensive review of DIB correlations from a statistical perspective, but also provide the first characterization of the sky distribution and kinematic properties of these three DIBs. In addition, we release a total of 27 598, 42 735, and 21 635 MQ measurements of the DIBs 45780, 45797, and 46614, respectively, corresponding to 23 920, 36 825, and 18 145 unique stars. These marginal or null measurements allow further studies, such as 3D mapping or the investigation of additional correlations. Our primary conclusions are as follows:

• (i)

  We release to the community the largest sample of HQ measurements of the DIBs λ5780 and λ6614 to date, and by virtue of these data, we also present their Galactic distribution for the first time, at a resolution of 0.92°/pixel.

• (ii)

  In our statistical analysis with a sample size on the order of tens of thousands, we provide more consistent and robust EWs of DIBs per unit E(B – V). For the DIB λ5780, we determine a value of 0.565 Å mag⁻¹, which is larger than those found in previous works, while for the DIBs λ5797 and λ6614, we obtain values of 0.176 and 0.256 Å mag⁻¹, respectively, which agree with previous findings.

• (iii)

  Without assuming the co-motion of the DIB carriers and the ISM, we use a kinematic fitting method to obtain precise measurements of the central wavelengths of the DIBs 45780, λ5797, and λ6614, which are 5780.48 ± 0.01, 5796.94 ± 0.02, and 6613.64 ± 0.01 Å, respectively.

• (iv)

  For our HQ sources of the DIBs λ5780 and λ6614, the kinematic distances of their carriers, characterized by their velocity curves, are both mainly between 0.5 and 1.5 kpc from the Sun.

Data availability

The final catalog of the DIB measurements of the 7681 common HQ samples, and the catalogs of the HQ and the MQ samples for the DIBs 45780, 45797, and 46614 are publicly available at https://nadc.china-vo.org/res/r101404/ in FITS tables. The codes used to perform the pipeline of measuring DIB are also provided via https://github.com/iScottMark/LAMOST_DIB.

The other data used in this work and the codes that generate the figures and support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements

This work is supported by National Key R&D Program of China (Grant NO. 2019YFA0405102) and National Natural Science Foundation of China (grant Nos. 12261141689, 11973060 and 12090044). Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the Project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. He Zhao acknowledges the National Natural Science Foundation of China (grant No. 12203099), the China Postdoctoral Science Foundation (No. 2022M723373), and the Jiangsu Funding Program for Excellent Postdoctoral Talent. Software: Astropy (Astropy Collaboration 2013; Astropy Collaboration 2018; Astropy Collaboration 2022), TOPCAT (Taylor 2005), VizieR (Ochsenbein et al. 2000), Simbad (Wenger et al. 2000), emcee (https://emcee.readthedocs.io/en/stable/), uncertainties (https://uncertainties.readthedocs.io/en/latest/).

Appendix A Instrumental broadening in LAMOST LRS

Compared with the high-resolution spectrographs, the low- resolution spectrographs have a more significant instrumental broadening effect. For the LAMOST LRS, the instrumental broadening effect is mainly caused by the slit width, ∆λ. Although Cui et al. (2012) and Luo et al. (2012) have reported the design of LAMOST including ∆λ, it was an experimental parameter at that time. The actual resolution of LAMOST LRS is about 1800 (private communication, Yonghui Hou, 2023) at 5500 Å in red end (from 3300 Å to 5900 Å) and at 7500 Å in blue end (from 5700 Å to 9000 Å). The resolution R is defined as λ/∆λ, where λ is the corresponding central wavelength. Hence, ∆λ, namely the instrumental FWHM, at 5500 Å and 7500 Å are roughly 3.06 Å and 4.17 Å.

We refer to the canonical FWHM of DIB (see the column FWHM_ref in Table A.1) from Vogrinčič et al. (2023) as our benchmark to estimate the practical Gaussian width of DIB against the LAMOST LRS. The broadened Gaussian width of DIB can be approximated by Eq. (A.1). $${\sigma }_{\text{DIB}}=\sqrt{{\sigma }_{\text{ins}}^2+{\sigma }_{\text{ref}}^2},$$(A.1)

where $\sigma =\text{FWHM}/2\sqrt{2\ln 2}$ ln 2, and the subscripts ins and ref stand for the instrumental broadening and the reference. Table A.1 provides the practical Gaussian widths of DIBs λ5780, λ5797, and λ6614 against the LAMOST LRS, which are annotated in Fig. 5.

References

All Tables

All Figures

	Fig. 1 Schematic flowchart of the pipeline for measuring the DIBs λ5780, λ5797, and λ6614 of cool stars in the LAMOST LRS DR10. Four colors encode the four parts of the pipeline, and the output of each modular function is used for the input of the subsequent one.
In the text

	Fig. 2 Distribution of the three atmospheric parameters, namely, Teff (left panel), log ɡ (middle panel), and [Fe/H] (right panel), for the target and reference datasets. The JS divergences of the three parameters between the two datasets are annotated in each panel.
In the text

Fig. 3 Example measurement of the DIBs λ5780 and λ5797 from a target spectrum. From top to bottom, (i) the first panel shows the target spectrum to be measured, with the red text indicating its obsid (the unique ID of the spectrum in the LAMOST catalog) and atmospheric parameters. (ii) The second panel provides the first 30 best neighboring spectra of the target spectrum, which are used to build the template spectrum displayed in the third panel. The purple-shaded areas in the first and second panels are the masked regions that do not participate in the similarity computation. (iii) The third panel illustrates the target and template spectrum after local normalization. The most prominent stellar absorptions that may affect the DIB extraction are annotated in the panel. (iv) The fourth panel presents the ISM residual spectrum (in black) of the target and the fit (in red) of the DIBs λ5780 and λ5797 by MCMC. The shallow DIBs λ5762 and λ5815 in the ISM residual spectrum are also annotated but not fitted in the panel.

In the text

Fig. 4 Examples of the local renormalization for two target spectra. The upper panel is the target spectrum within the local region covering the DIBs λ5780 and λ5797, and the lower panel is the target spectrum within the local region including the DIB λ6614. In each panel, the black line represents the original spectrum, and the orange line is the renormalized spectrum. The obsid of the target spectrum and its atmospheric parameters are also annotated in each panel.

In the text

Fig. 5 Two-dimensional histograms of the measured DIBs λ5780, λ5797, and λ6614 as a function of the Gaussian width σ and the central wavelength λ in the rest frame. (a) The three upper panels present the density of all the samples that may have DIBs, (b) and the other three lower panels show the density of the HQ samples for which we have reliable DIB measurements. All of the histograms have the same bin size of 100. The red horizontal dashed lines denote the central wavelengths of the DIBs, which are fitted using the method described in Sect. 5.3. The red vertical dashed lines indicate the typical Gaussian widths of the DIBs considering the instrumental broadening of LAMOST LRS (see Appendix A), which are taken from Vogrinčič et al. (2023).

In the text

Fig. 6 Comparison between the integral EW and the fitted EW of the DIBs λ15780, λ5797, and λ16614 for their HQ samples, from the subfigures (a) to (c). For each subfigure, the upper panel shows the Gaussian KDE number density plot of the integral EW as a function of the fitted EW, with the red dashed line representing the one-to-one relation. The lower panel displays the Gaussian KDE number density plot of the difference between the integral EW and the fitted EW, namely ∆EW = EWint − EWfit as a function of the fitted EW, with the green dashed line representing the zero benchmark. The distribution of ∆EW is also overplotted in the upper panel with its mean and standard deviation annotated.

In the text

Fig. 7 Three distinct colors encode the three DIBs λ5780, λ5797, and λ6614 for the HQ samples with over two epochs. The panels in the first column present the distribution of the number of epochs. The panels in the second column depict the cumulative distribution of the standard deviation of the EWs (σEW). The panels in the third column display the cumulative distribution of the coefficient of variation (σEW/EWmean). The panels in the fourth column show the scatter of the σEW/EWmean as a function of the EWmean.

In the text

Fig. 8 Trend of the standard deviation of the residuals between the flux of the ISM residual spectrum beyond the DIB region and that of the continuum as a function of S/N in the upper panels, and Teff in the lower panels. Three distinct colors encode the three DIBs λ5780, λ5797, and λ6614. In each panel, a violin plot is used to present the distribution of the standard deviation of the residuals for the binned HQ samples. The median of each bin is annotated with a white dot, and the 2.5 and 97.5 percentiles are marked with black dots. The red dashed line connects the median of each bin to show the trend of the standard deviation of the residuals.

In the text

Fig. 9 Galactic distribution of the DIBs λ5780 and λ6614 for the HQ samples. The upper panels show the number density of the background stars. The lower panels display the sky distribution of the strength of the DIBs with the color bar representing the median EW in each pixel. The black contours of the corresponding median E(B − V) at the level of 0.35 and 0.7 magnitude are overplotted in the median EW map. The resolution of the map is about 0.92° (Nside = 64) calculated by the HEALPix package.

In the text

Fig. 10 Gaussian KDE number density plot of the EWs of the DIBs λ5780, λ5797, and λ6614 (from top to bottom) as a function of the extinction E(B − V) derived from Green et al. (2019). For each panel, the red solid line represents a linear fit to the median EW (which is marked with a white dot, and its standard deviation marked is with a black error bar) in each ΔE(B − V) = 0.05 mag bin. The green dashed line is the same as the red solid line, but the intercept is forced to be zero. The coefficients of the linear fit and the Pearson correlation coefficient between the EW and E(B − V) are also annotated in each panel.

In the text

Fig. 11 As in Fig. 10, but for the Gaussian KDE number density plot of EW V.S. EW. From the top panel to the bottom panel, the DIBs λ5780, λ5797, and λ6614 are compared with each other. For each panel, the red solid line represents a linear fit to the median EW in each ΔEW = 0.05 Å bin. Except for the bottom panel, the bin size of the EW is ΔEW = 0.02 Å. The coefficients of the linear fit and the Pearson correlation coefficient between the EWs are also annotated in each panel.

In the text

Fig. 12 Observed central wavelength (λobs) of the DIBs λ5780, λ5797, and λ6614 as a function of the angular distance (∆ℓ) from the GAC, from top to bottom. The gray dots and their error bars represent the measurements of the central wavelength of DIB and the corresponding uncertainties. The white dots and their black error bars denote the medians and the standard deviations in each ∆ℓ = 1° bin. The red lines are the linear fits to the medians of each bin using the least squares method. The fitting results, that is, the slopes k, the intercepts b, and the final central wavelengths λ0 considering the solar motion, are annotated in each panel.

In the text

Fig. 13 Tow-dimensional histogram of the LSR velocity of the DIBs 45780 (upper panel) and 46614 (lower panel) as a function of the Galactic longitude. The medians and the standard deviations in each ∆ℓ = 10° bin are denoted by the white dots and their black error bars. The Galactic velocity curves with different heliocentric distances d calculated by the Model A5 of Reid et al. (2019) are overlaid in the plots.

In the text