© The Authors 2023

1 Introduction

Diffuse interstellar bands (DIBs) are a set of ubiquitous interstellar absorption features that primarily exist in the optical and near-infrared wavelength range (about 0.4–2.4 µm) of the spectra of stars (Fan et al. 2019; Hamano et al. 2022; Ebenbichler et al. 2022), galaxies (e.g. Monreal-Ibero et al. 2018), and distant quasars (e.g. Monreal-Ibero et al. 2015). DIBs presumably originate from the electronic transitions of carbon-bearing molecules and are now well recognized as the signatures of complex molecules (Tielens 2014), although the exact species of their carriers remain largely unidentified. Based on highresolution spectrometry, DIBs can be used to probe the variation of the interstellar environments in clouds (Cordiner et al. 2013), and reveal the physical and chemical process of the interstellar medium (ISM; Welty 2014) and the formation and development of chemical complexity in space (Tielens 2014).

While in the past DIB studies were mainly concentrated on small dedicated sample sizes, large Galactic spectroscopic surveys during the last decade such as Gaia–ESO, APOGEE, SDSS, RAVE, and GALAH have given rise to numerous DIB detections in our Milky Way (e.g. Puspitarini et al. 2015; Zasowski et al. 2015; Elyajouri & Lallement 2019; Lan et al. 2015; Baron et al. 2015b; Kos et al. 2013; Vogrinčič et al. 2023), which has enabled the investigation of the kinematics (Zasowski et al. 2015; Zhao et al. 2021b) and three-dimensional (3D) distribution (Kos et al. 2014) of DIB carriers that trace the large-scale structures of our Milky Way. The latest Gaia data release 3 (DR3) contains the largest catalogue so far for the DIB at 862.1 nm in air (λ862.1). DIB λ862.1 was detected and measured in Gaia Radial Velocity Spectrometer (RVS; Seabroke et al. 2021) spectra of individual stars, by the General Stellar Parametrizer from spectroscopy (GSP-Spec) module (Recio-Blanco et al. 2023) of the Astrophysical parameters inference system (Apsis; Bailer-Jones et al. 2013; Creevey et al. 2023). Productive analysis and results of the DR3 DIB catalogue were performed and presented in Gaia Collaboration (2023), hereafter S23 in which we built a map of the median DIB λ862.1 strength, covering all the longitudinal directions, mainly within 3 kpc from the Sun and 1 kpc above and below the Galactic plane. The rest-frame wavelength of λ862.1 was determined as λ₀ = 8623.23 ± 0.019 Å in vacuum. An average scale height for the carrier of λ862.1 was estimated as ${98.60}_{-8.46}^{+11.10}$ pc assuming a simple exponential distribution of the carrier perpendicular to the Galactic plane. The longitudinal variations of the radial velocity of the λ862.1 carrier were clearly shown as well.

Nevertheless, due to the limitation of the signal-to-noise ratio (S/N) of the individual RVS spectra, DIB λ862.1 could be successfully measured in only ~10% of RVS objects. Restricting the DIB sample further to reliable high-quality measurement reduced the sample to only ~ 140 000 (see S23 Sect. 3 for the definition of the high-quality sample). Another major limitation is the usage of synthetic spectra in the process of measuring the DIB (Recio-Blanco et al. 2023), assuming that they represent the stellar components in the observed spectra perfectly. The complexity of the description of the full stellar physics in stellar atmosphere models as well as uncertainties in the atomic line list could easily lead to inappropriate modelling of stellar lines around the DIB signal, which would introduce further uncertainties in the fitting of the DIB profile.

To overcome the disadvantage of using synthetic spectra and the constraint of the S/N of individual RVS spectra for DIB measurement, we developed an Apsis module of the Gaia Data Processing and Analysis Consortium (DPAC), called Gaia DPAC/DIB-Spec, which processed 6.8 million RVS spectra and conducted new DIB measurements. To avoid using synthetic spectra, Kos et al. (2013) developed a data-driven method, called the BNM, to detect DIB signals in the spectra of late-type stars using artificial stellar templates constructed from real spectra observed at high latitudes that have a similar morphology to the spectrum of the target star but are likely to be free of the signature of ISM, in analogy to the distribution of interstellar extinction. Zhao et al. (2022), hereafter Z22 applied BNM to the publicly available RVS spectra within Gaia DR3 where they confirmed the presence of the weak DIB at 864.8 nm¹ in RVS spectra. On the other hand, stacking spectra in a given spatial volume would significantly increase the S/N of spectra and thus allows the detection of much weaker DIB signals (e.g. Kos et al. 2013; Baron et al. 2015b; Lan et al. 2015; Zhao et al. 2022, 2023). With these two techniques, DIB-Spec could detect and measure DIB signals in more distant zones compared to the results in DR3 and reveal the large-scale spatial distribution of the DIB carriers (stacking reduces the spatial resolution). Furthermore, an increased S/N of stacked spectra enables DIB-Spec to measure DIB λ864.8 as well. The Gaia Focused Product Release (FPR) contains the parameters of two DIBs, λ862.1 and λ864.8, fitted by DIB-Spec and the stacked ISM spectra (spectra only containing interstellar features) in each defined VOlume piXEL (voxel, or 3D display element). The aim of this paper is to introduce the DIB-Spec module and present a preliminary analysis of the DIB measurements.

The paper is outlined as follows: A brief description of the input Gaia RVS spectra is provided in Sect. 2. Section 3 explains the pipeline of the DIB-Spec module, including the construction of target and reference samples, deriving ISM spectra, the stacking of ISM spectra, and the DIB measurement. Section 4 describes and discusses in detail the outputs of DIB-Spec, the fitted DIB parameters, and their stacked ISM spectra. The performed validation of DIB-Spec outputs are presented in Sect. 5. We discuss in Sect. 6 detections of DIB λ862.1 inside the Local Bubble. We finish with some caveats about the usage of the DIB-Spec outputs in Sect. 7 and our main conclusions in Sect. 8.

Fig. 1 Schematic workflow of the DIB-Spec module.

2 Input Gaia RVS data

The input data for DIB-Spec are based on the Gaia RVS spectra that were processed by the Gaia DPAC Coordination Unit 6 (CU6). The processing includes removal of cosmic rays, wavelength calibration, normalization to the continuum, and resampling from 846 to 870 nm with a spacing of 0.01 nm (2400 wavelength bins, Sartoretti et al. 2018, 2023). The resolving power is R = λ/Δλ ~ 11500 (Cropper et al. 2018). The RVS spectra were then processed by GSP-Spec to estimate their stellar atmospheric parameters (effective temperature T_eff, surface gravity log ɡ, metallicity [M/H], and [α/Fe]) without taking into account any post-processing steps. The RVS spectra were renormalized and rebinned by GSP-Spec, from 2400 to 800 wavelength bins, sampled every 0.03 nm to increase their S/N. In total, there are 6 862 982 RVS spectra with S/N ⩾ 20 (S/N is provided by the CU6 analysis Seabroke et al. 2021). These RVS spectra, normalized by GSP-Spec, and their stellar atmospheric parameters, as well as three other parameters (parallax ϖ, stellar radial velocity V_rad, and the velocity uncertainty ${\sigma }_{V_{\text{rad}}}$), are used as the basic input for DIB-Spec. We want to stress that the stellar atmospheric parameters are only needed to speed up the procedure of deriving ISM spectra (see Sect. 3.2) but are not indispensable in the data-driven method. We also note that not all the 6.8 million spectra have published stellar parameters in Gaia DR3 (about 1.2 million were filtered out), but their parameter distributions are shown in this paper, like Fig. 5.

3 The pipeline of the DIB-Spec module

Figure 1 shows the general overview of the pipeline of DIB-Spec module, which was operated by CU8 at the Data Processing Center CNES (DPCC) in Toulouse, France.

DIB-Spec contains three main steps, which are explained in detail below. First, based on the quality of the RVS spectra and the results of GSP-Spec, DIB-Spec builds two samples: a sample of ‘target’ stars whose spectra are expected to contain DIB signals, and a sample of ‘reference’ stars at high latitudes whose spectra have presumably no DIB features. Second, the reference spectra are used to derive the ISM spectra for each target star. Finally, DIB-Spec stacks the ISM spectra of individual target stars in each voxel and fits the two DIBs in the stacked ISM spectra.

Fig. 2 Galactic spatial density distribution of the 6 143 681 RVS spectra used in DIB-Spec. This HEALPix (Górski et al. 2005) map has a level of 7, corresponding to a spatial resolution of 0.46°.

Fig. 3 Gaia G–band magnitude distribution of the 6 143 681 RVS objects used in DIB-Spec.

3.1 Target and reference samples

The RVS objects must fulfil a set of requirements in order to be retained in the DIB-Spec analysis. First, the RVS objects must have a measurement of V_rad derived from CU6 with ${\sigma }_{V_{\text{rad}}}⩽5$ km s⁻¹, because their ISM spectra have to be shifted from the stellar frame (as provided by CU6) back into the heliocentric frame before stacking. Then, the parallax must be above 0.1 mas, which corresponds to objects within 10 kpc if they have small parallax uncertainties. We expect to detect the DIB signals in very distant zones, but it should be noted that most of the RVS objects (91.5%) are within 4 kpc (see Fig. 5). These criteria result in a set of 6 143 681 objects to be used. Their number density distribution in Galactic coordinates and G–band magnitude distribution are shown in Figs. 2 and 3, respectively. Most of these RVS objects are located close to the Galactic plane and their G–band magnitudes are mainly within 8–14 mag.

The selected stars are separated into two sets according to their Galactic latitudes: 5 983 289 target stars with |b| < 65° and, 160 392 reference stars with |b| ⩾ 65°. The DIB strength has a strong dependence on the Galactic latitude since the stars at high latitudes usually contain very weak DIB signals in their spectra (see the maps of DIB strength in Lan et al. 2015, Baron et al. 2015a, and S23 for example). In Gaia DR3, we detected λ862.1 in 247 sightlines with |b| ⩾ 65° in the high-quality sample (~140 000 in total) of the DIB catalogue (see Sect. 3 in S23). Their mean depth is only 1.5% of the continuum with a very small mean A₀ (monochromatic extinction at 541.4 nm Fitzpatrick 1999) of 0.104 mag estimated by the Total Galactic extinction (TGE) map of Gaia DR3 (Delchambre et al. 2023). The mean A₀, on the other hand, for the reference stars is 0.103 mag (see Fig. 4). Therefore, 1.5% would be the maximum error we expect in the derived ISM spectra within the DIB region of target spectra introduced by the possible DIB signals at high latitudes when using the reference spectra. It should be noted that the real errors could be much lower because the TGE map would overestimate A₀ for nearby stars. Furthermore, as discussed in Kos et al. (2013), any stars with unusually high extinctions or strong DIBs will be averaged out because the stellar template is produced by averaging several reference spectra. In summary, the reference spectra could be treated as DIB-free spectra and used to model stellar components in the DIB region of the target spectra.

The target and reference stars have a similar distribution in T_eff, log ɡ, and [M/H] (Fig. 5). Therefore, most of the target spectra should have enough neighbour candidates selected by comparing their atmospheric parameters (see Sect. 3.2 for more details) to model their stellar components. The target and reference samples also have a similar distribution of spectral S/N, although the reference sample has a higher proportion of nearby stars (87.5% within 1 kpc) compared to the target sample (47.1%).

Fig. 4 Distribution of A0 from Gaia TGE map (Delchambre et al. 2023) for 160 392 reference stars.

3.2 Derivation of ISM spectra for individual targets with the BNM

To detect and measure the two DIBs, λ862.1 and λ864.8, in target spectra, the stellar components need to be subtracted by the BNM. We follow here the main principles in Kos et al. (2013). For each target, BNM selects a set of reference spectra with similar spectral morphology to that of the target spectrum and combines them by averaging their normalized flux weighted by the spectral S/N to build the stellar template for the target spectrum. In practice, BNM first constrains the stellar atmospheric parameters (from GSP-Spec) of the reference stars, to be around the values of the target, in the ranges T_eff ± 20%, log ɡ ± 0.6 dex, and [M/H] ± 0.4 dex. The ranges follow those in Kos et al. (2013). This step reduces the number of considered reference stars and consequently speeds up the procedure.

To compare the spectral morphology between the target and reference spectra, we first calculate the absolute difference and their normalized fluxes at each wavelength bin. Then we take the weighted mean of these flux differences, where the inverse of the weighted mean can be seen as the similarity of the spectral morphology. The weights are set to 0.0 in the range 860–868 nm where the two DIBs are located, and to 1.0 elsewhere. An exception to this is the central regions of the Ca II lines (849.43–851.03 nm and 853.73–855.73 nm) where the weights are set to 0.7 because the Ca II lines are dominating the RVS spectra of cool stars but do not affect the DIB measurements which rely more heavily on the modelling of the Fe I lines (Kos et al. 2013; Z22). Reference spectra with similarity values smaller than three times the S/N of the target spectrum are discarded. The remaining spectra are sorted by decreasing similarity, and the first 200 at most are selected as the best neighbours for the given target spectrum, as best neighbours do not significantly increase the accuracy of the stellar template (Z22). We note that the neighbourhood used in BNM refers to the parameter space rather than the physical space, that is the best neighbours and target should have similar stellar parameters and spectral morphology, but not necessarily similar sky positions. The best neighbours are averaged with their spectral S/N as weights in order to build a stellar template. The ISM spectrum is defined as the target spectrum divided by the stellar template. If the number of the best neighbours is less than 10, no stellar template or ISM spectrum is generated. DIB-Spec successfully generated stellar templates for 4 595 489 target spectra (76.8% of the target sample) and consequently derived their ISM spectra.

An illustration of BNM is shown in Fig. 6: the main steps are summarized on the left side and an example is presented on the right side. We made a test of BNM with the reference sample (see Appendix A for details) and found that the performance depends strongly on the spectral S/N. For S/N > 50, the average flux residuals between the RVS spectra and the derived stellar templates are mainly within 0.02 in the DIB window (861.2–866.0 nm). While the residuals also vary with the stellar atmospheric parameters, which is caused by the change in the number density of stars and the accuracy of GSP-Spec light for different types of stars. As BNM constrains the SED similarity by S/N and requires at least ten best neighbours for a given target spectrum, a target star that has an incorrect spectral type (e.g. an early-type star gains a very low T_eff) or rare reference stars in its vicinity will not get a generated stellar template.

Fig. 5 Distribution of Teff, log ɡ, [M/H], the distance of stars, and spectral S/N for both reference stars (in red) and target stars (in blue). These parameters have not been filtered out in the post-processing.

3.3 Stacking and fitting in each voxel

The ISM spectra of individual targets are stacked in each voxel according to their equatorial coordinates (α, δ) and the distances (1/ϖ) to the target stars to get a much higher spectral S/N than individuals’ for a reliable measurement of the two DIBs. As the DIB carriers should be located between the observer and the background star, we measure in each voxel the sum of the DIB materials from the observer up to the voxel. Considering the size of each voxel, the measured DIB feature would be averaged according to the spatial and distance resolution of the voxels. The pixelation in (α, δ) is done at level 5 in the HEALPix² scheme (Górski et al. 2005), corresponding to 12 288 pixels with an equivalent spatial resolution of 1.8°. For the binning in distance, we selected some adaptive steps instead of a uniform bin size because the distribution of the target stars is inhomogeneous in space. Furthermore, as the resolution of the sky (α, δ) is fixed as 1.8°, we would like to get as high a resolution in distance as possible, especially for nearby regions.

The S/N of the stacked spectra is a good criterion for determining the size of the distance bins, but the main restriction is that DIB-Spec runs with defined distance bins as input. Therefore, we cannot use the derived ISM spectra of individual stars, which were intermediate products of DIB-Spec, to find the best solution for the adaptive distance bins. Instead, $S\text{/}{N}^{\prime }=\sqrt{{\displaystyle \sum {\text{S/N}}_i^2}}$, where S/N_i is the S/N of the ith individual RVS spectra in a voxel, is used to characterize the stacked ISM spectra in each voxel. The size of the distance bins was determined by the requirement that in each bin the number of HEALPixels having S/N′ ⩾ 200 is larger than 85% of the total when distance is smaller than 4.5 kpc. And a constant size of 0.5 kpc was applied for bins between 4.5 and 10 kpc. In this way, there are finally 29 distance bins whose ranges are listed in Table 1, together with the number of HEALPixels in each bin with at least one target star. Certainly, S/N′ would be much higher than the true S/N of the stacked ISM spectra, as it does not account for the uncertainty caused by the subtraction of stellar components. In fact, only ~50% voxels within 4.5 kpc have S/N ⩾ 200 (here is the truly calculated S/N of the stacked ISM spectra, see its definition below). Table 1 lists the fraction of S/N ⩾ 200 in each distance bin. The fraction could exceed 85% when we consider S/N ⩾ 150, but it will also significantly reduce beyond 1.67 kpc.

The ISM spectra are stacked in each voxel by taking, in each wavelength bin, the median value of the fluxes in order to reduce the influence of the outliers. The flux uncertainty at each wavelength bin is taken as the mean of the individual flux uncertainties divided by $\sqrt{N_{\text{tar}}}$. The S/N of the stacked ISM spectra is calculated between 860.2 and 861.2 nm as mean(flux)/std(flux).

For the profiles of the two DIBs, λ862.1 was usually assumed to have a Gaussian profile in previous studies (e.g. Kos et al. 2013; Zhao et al. 2021a). Although some departures from a Gaussian profile were reported (e.g. Krełowski et al. 2019), the origin is more like the superposition of multiple DIB clouds, as no evidence supports an intrinsic asymmetry of the λ862.1 profile. Thus, we still assume a Gaussian profile for λ862.1 in this work and treat the possible departures as a source of uncertainty. Rare studies of λ864.8 make it harder to determine the shape of its profile. Zhao et al. (2022) chose a Lorentzian profile as it showed smaller residuals on the ISM spectra compared to a Gaussian profile. Additionally, the Lorentzian profile has been proved to be appropriate for the very broad DIB λ442.8 (Snow 2002) while λ864.8 has a broad profile as well. Therefore, DIB-Spec models the profiles of the two DIBs in stacked ISM spectra by a Gaussian function (Eq. (1)) for λ862.1, a Lorentzian function (Eq. (2)) for λ864.8, and a linear function for continuum (Eq. (3)): $$G\left(\lambda ;𝒟,{\lambda }_{\text{DIB}},{\sigma }_{\text{DIB}}\right)=-𝒟\times \exp \left(-\frac{\left(\lambda -{\lambda }_{\text{DIB}}\right)}{2{\sigma }_{\text{DIB}}^2}\right),$$(1) $$L\left(\lambda ;𝒟,{\lambda }_{\text{DIB}},{\sigma }_{\text{DIB}}\right)=\frac{-\left(𝒟{\sigma }_{\text{DIB}}^2\right)}{{\left(\lambda -{\lambda }_{\text{DIB}}\right)}^2+{\sigma }_{\text{DIB}}^2},$$(2) $$C\left(\lambda ;a_0,a_1\right)=a_0\times \lambda +a_1,$$(3)

where 𝒟 and σ_DIB are the depth and width of the DIB profile, λ_DIB is the measured central wavelength, a₀ and a₁ describe the linear continuum, and λ is the wavelength. Subscripts ‘862.1’ and ‘864.8’ are used to distinguish the profile parameters of the two DIBs. The full DIB model, M_Θ, is the sum of Eqs. (1)–(3), where Θ = {𝒟_862.1, λ_862.1, σ_862.1, 𝒟_864.8, λ_864.8, σ_864.8, a₀, a₁} are the adjusted model parameters. Given the stacked ISM spectrum {λ, y, σ_y) where y is the normalized flux and σ_y is the flux uncertainty, and the unnormalized posterior probability density function (PDF) is P(Θ|y) ∝ P(y|Θ)P(Θ). P(y|Θ) is the likelihood: $$P\left(y|\text{Θ}\right)=\frac{1}{\sqrt{2\pi }{\sigma }_y}\exp \left[\frac{1}{2{\sigma }_u^2}{\left(y-M_{\text{Θ}}\right)}^2\right],$$(4)

and P(Θ) represents the prior distributions of the parameters. Flat and independent priors were applied for the DIB parameters, that is 0–0.2 for 𝒟_862.1 and 𝒟_864.8, 861.2–863.0 nm for λ_862.1, 863.0–866.0 nm for λ_864.8, 0.01–0.5 nm for σ_862.1, and 0.1–1.5 nm for σ_864.8. No priors were used for a₀ and a₁.

The optimization of the eight parameters – three for each DIB plus two for the continuum – was done by sampling their full posterior distributions. We note that during the optimization, a masked region was applied between 866.0 and 866.8 nm for the residuals Ca II that were caused by downweighting the central regions of Ca II lines in the BNM (see Sect. 3.2). A Markov chain Monte Carlo (MCMC) procedure (Foreman-Mackey et al. 2013) was performed to implement the parameter estimates. The initial guess of 𝒟 and λ_DIB were determined by averaging flux near the lowest point (5 pixels) within 861.2–863.0 nm and 863.5–866.0 nm for λ862.1 and λ864.8, respectively. And initial σ_DIB was fixed as 0.12 nm for λ862.1 and 0.4 nm for λ864.8. a₀ and a₁ were initially set as 0 and 1 as well. One hundred walkers were progressed for 250 steps to complete the burn-in stage. The best fits derived by the last 200 steps were then used as the initial conditions to sample the posterior with 100 walkers and 2000 steps. The best estimates and their statistical uncertainties were taken in terms of the 50th, 16th, and 84th percentiles of the posterior PDF drawn from the last 1500 steps. According to Eqs. (1) and (2), the equivalent width (EW) for λ862.1 is calculated as $${\text{EW}}_{862.1}=\sqrt{2\pi }\times {𝒟}_{862.1}\times {\sigma }_{862.1},$$(5)

and for λ864.8 as $${\text{EW}}_{864.8}=\pi \times {𝒟}_{864.8}\times {\sigma }_{864.8}.$$(6)

The lower (16%) and upper (84%) confidence levels of EW were estimated by 𝒟 and σ_DIB drawn from the MCMC posterior samplings. Specifically, each {𝒟, σ_DIB} pair sampled during the MCMC fitting was used to calculate one value of EW, and finally we had an EW distribution from which the 16th, 50th, and 84th percentile values were calculated.

Figure 7 shows examples of the fits for five stacked ISM spectra in voxels towards the same direction, that is HEALPixel number of 10 450, corresponding to Galactic coordinates of the voxel of (ℓ_c, b_c) = (322.43°, −0.44°). With increasing voxel central heliocentric distance (d_c) from top to bottom, EW_862.1 and mean E(BP − RP) measured by GSP-Phot (Andrae et al. 2023) both increase and show a good correlation with each other (these values are indicated in Fig. 7). The profiles of λ862.1 are strong and prominent in all these voxels. On the other hand, the S/N of the stacked ISM spectra and the number of target spectra (N_tar) decrease with distance. Consequently, the fit to the profile of λ864.8 becomes worse in voxels with d_c = 2.40 and 3.23 kpc (bottom panels in Fig. 7) as the very broad and shallow profile of λ864.8 is more affected by the stellar residuals and uncertainties introduced by BNM and stacking than λ862.1.

The corner plots of these examples, presenting the one- and two-dimensional projections of the posterior distributions of the fitted parameters, are shown in Figs. B.1–B.5, respectively. The corner plots clearly show that the MCMC fitting is converged and the posterior PDF of all the parameters is Gaussian. 𝒟 and σ_DIB are correlated with each other during the fitting for both λ862.1 and λ864.8. The depth and central position of λ862.1 (𝒟_862.1 and λ_862.1) are not sensitive to the continuum placement (a₀ and a₁), while σ_862.1 presents a weak correlation with a₀ and a₁. The continuum placement has a heavier effect on the broad and shallow profile of λ864.8.

Fig. 6 Illustration of deriving ISM spectra for individual targets. Left: Schematic view of each main step. Detailed explanations are in Sect. 3.2. Right: Example for a target (Gaia ID: 4169294186397475072) having $T_{\text{eff}}={4445}_{-11}^{+15}$ K, log ɡ = 1.86 ± 0.03 dex, $\left[\text{M/H}\right]={0.12}_{-0.01}^{+0.02}$ dex, and spectral S/N of 110.3. From top to bottom panels: 1) Observed RVS spectrum of this target; 2) First 24 best neighbours for this target; 3) The black line is the target spectrum and the blue line is its stellar template built by averaging 200 best neighbours; 4) Derived ISM spectrum for this target. The positions of the two DIBs are also indicated. The orange shades indicate the regions of two Ca II lines where the weights are set to 0.7 when comparing the target spectrum and reference spectra. The cyan region (860–868 nm) indicates the spectral window used for fitting the two DIBs, with a masked region (green, 866.0–866.8 nm) during the fitting because of the residuals of the Ca II line.

Fig. 7 Examples of the fits to DIBs λ862.1 and λ864.8 in stacked ISM spectra in five voxels in the same direction, whose HEALPix number (Ipix = 10 450) and GC (ℓc, bc) = (322.43°, −0.44°) are marked in the top panel. The black and red lines are the ISM spectra and fitted DIB profiles, respectively, normalized by the fitted linear continuum. The error bars indicate the flux uncertainties at each pixel. Orange indicates the region that was masked during the fittings. The central heliocentric distance (dc), the number of target spectra (Ntar), mean E(BP − RP), EWs of two DIBs, and the S/N of the stacked ISM spectrum in each voxel are indicated as well.

4 DIB-Spec outputs

DIB-Specsuccessfully fitted the two DIBs in 235428 voxels. This is less than the total number (356 352) of a level 5 HEALPix binning × 29 distance bins, which is due to discarding voxels with no target stars (targets with failed generated stellar templates were not included as well). Each voxel on average contains 20 target stars, with a maximum number of 386. There are 23 875 voxels (10.1%) that only contain one target spectrum. The outputs of DIB-Spec are arranged into two tables in the Gaia Archive³: (1) ‘interstellar_medium_params’: the parameter table listing the fitted DIB parameters in each voxel; (2) ‘interstellar_medium_spectra’: the spectra table containing the stacked ISM spectra in each voxel. The column names, units, and descriptions of the two tables are given in Tables 2 and 3 respectively. In Table 2, the symbols of the full DIB parameters (Eqs. (1)–(3)) and EWs are indicated in the description of the corresponding columns.

There are some notes for the DIB-Spec outputs:

1. As the column names were pre-defined, the DIB at 862.1 nm was cited as ‘8620’ in related names, but we prefer to use ‘862.1’ in the descriptions and in the context of this Gaia product.

2. Second, the lower and upper confidence levels of the intercept of the linear continuum (‘dibcont_a1’) are not included in the parameter table because they were not defined at the time of the processing and archive table definition.

3. As in DIB-Spec the HEALPix binning was done in the equatorial system, following the Gaia convention⁴, the Galactic coordinates of the voxel centre (‘lc’ and ‘bc’ in the table) were converted from the equatorial coordinates of the centre of each HEALPixel.

4. The fitted DIB parameters result from the integration of their carriers from the voxels to the observer, like dust extinction, rather than from one voxel to the next.

5. About 5.4% (12 692) of stacked spectra have null flux uncertainties. This is due to the fact that the first pixels of the individual RVS spectra for stacking do have zero values. Their flux uncertainties are fixed as 0.01 for the MCMC fittings.

6. The spectra table contains 62 859 276 rows which equal 235 428 voxels × 267 wavelength bins of the stacked spectra, that is each row in the spectra table contains information of one wavelength bin. A simple Python script shown in Appendix C can be used to convert the spectra table to a fits file, in which each row stands for one stacked ISM spectra.

Below, we describe and discuss the fitted DIB parameters and their uncertainties.

4.1 S/N and DIB quality flag (QF)

The S/N of the stacked ISM spectra strongly affects the quality of the DIB fit. S/N is determined by the quality of individual RVS spectra, the number of target stars in each voxel, and the performance of the BNM on target spectra. For DIB-Spec results, 42% of the voxels have a stacked spectrum S/N > 200, but 59% of these voxels are within 1 kpc. The DIB signal will be more easily detected in spectra with higher S/N, while the DIB depth and strength are generally smaller in the solar neighbourhood than in distant zones. Thus, the goodness-of-fit for DIBs cannot be determined simply by S/N. On one hand, the distribution of S/N shown in Fig. 8 presents a dependence on both 𝒟 and σ_DIB. A main part of high S/N (≳200) is found with very small 𝒟, corresponding to nearby voxels containing weak or no DIB signals. On the other hand, large 𝒟 found in noisy spectra (low S/N) would indicate a needle-like spurious DIB signal caused by the random noise or the stellar residuals, especially in the regions with relatively small σ_DIB. Nevertheless, one can also find some regions with both high S/N and large 𝒟, indicating a pronounced DIB signal. Such regions are evident near σ_DIB ~ 0.22 nm for DIB λ862.1, and between 0.8 and 1.0 nm for DIB λ864.8, but not particularly obvious due to the fact that λ864.8 is much broader and shallower than λ862.1.

Considering both S/N and the shape of the fitted DIB profile ({𝒟, σ_DIB}), we generated quality flags (QF) to describe the quality of the fits. This idea comes from Elyajouri et al. (2016) and was applied to DIB λ862.1 in Gaia DR3 (Recio-Blanco et al. 2023; Gaia Collaboration 2023) as well. In the present work, we follow the same scheme for λ862.1 as in DR3, but some borders of the DIB parameters are redefined according to the DIB-Spec results. The scheme for λ864.8 contains newly defined borders, as it has been little investigated to date.

Figure 9 shows the flowchart for the generation of the quality flag for the two DIBs, respectively, with QF = 0 indicating the best fit and QF = 5 the worst fit. First, we require the measured λ_DIB to be consistent with a Doppler wavelength shift within about ±200 km s⁻¹ for the radial velocity of the two DIBs, which should be a reasonable velocity range for the ISM within 3 kpc from the Sun where most of the DIBs were detected. The DIB radial velocity is calculated from the estimated λ_DIB and the rest-frame wavelength λ₀ of the two DIBs reported in S23 and Z22. We note that the upper limit of λ_DIB for λ864.8 is a bit larger than a wavelength corresponding to 200 km s⁻¹ because in the literature this DIB was suggested to be located around 8648 Å in air. Secondly, based on S23 and Z22, 𝒟 is required to be smaller than 0.15 for λ862.1 and 0.10 for λ864.8, respectively. 𝒟 is then compared to R_C, which is defined as the standard deviation of the difference between observed and modelled flux of the ISM spectra, std(data – model), in a range from λ_DIB − 3σ_DIB to λ_DIB + 3σ_DIB to represent the noise level around the position of the DIB feature. We note that R_C is not published in the DIB-Spec tables, but one can recalculate R_C from the DIB parameters and the stacked ISM spectra. Last, we need to define some acceptable values of σ_DIB in order to evaluate the QFs, that is a “best range” of σ_DIB for QF = 0,1 and a “secondary range” for QF = 2, 3, 4. DIBs with σ_DIB out of these two ranges will be marked as QF = 5. The best range of σ_862.1 is set as 0.18–0.26 nm, a width of 0.04 nm to the median σ_862.1 (0.22 nm) that is about 2.5 times the average error of σ_862.1 for detections with QF = 0. The secondary range of σ_862.1 is 0.1–0.34 nm, threefold the span of the best range. Determining the ranges for σ_864.8 is much harder, as very few previous studies can serve as a reference. We determined a median σ_864.8 = 0.81 nm in Z22, but the median σ_864.8 in DIB-Spec is close to 1.0 nm when we impose S/N > 200. Thus, we select a loose range, 0.8–1.2 nm, for the best range of σ_864.8, and 0.6–1.4 nm as the secondary range. The distribution of σ_DIB is discussed in detail in Sect. 4.2. This suggests that careful analysis of the DIB-Spec results is needed to refresh our knowledge about λ864.8. It should be noted that when 𝒟 < R_C (case (g), Fig. 9), σ_DIB is only compared with the range of 0.1–0.18 nm for λ862.1 and 0.6–0.8 nm for λ864.8 (case (k), Fig. 9). This choice follows the idea that when 𝒟 is smaller than the noise level, a large σ_DIB is more likely a spurious feature caused by the fit to the successive correlated noise. Figure 10 shows the distribution of the QF for DIBs λ862.1 and λ864.8. Nearly half of λ862.1 have QF = 5, and this fraction reaches over 66% for DIB λ864.8, which dramatically reduces the DIB sample size for further analysis, especially for λ864.8. The much lower proportion of QF = 0, 2 for λ864.8 than QF = 1, 3 is due to its very small 𝒟, since only 4.3% of 𝒟_864.8 is larger than 3R_C. Thus, only 16.1% of DIB λ864.8 have QF ⩽ 2, the recommended high quality in Gaia DR3 (S23). DIB λ862.1 has a larger proportion (26.6%) because λ864.8 is much broader than λ862.1 and therefore under a heavier effect of the noise. The QF provides us with an evaluation of the fitted profile of the DIB feature and consequently the goodness-of-fit for the ISM spectra. But the present QF still contains some shortcomings, for example, the QF has discrete values and has hard and artificial borders for σ_DIB. Therefore, the QF should be used with caution for DIB λ864.8, which is not well-studied.

Figure 11 presents the distribution on the sky of QFs for the two DIBs in Galactic coordinates. Each level of QF is shown in one HEALPix map with level 5. The colour scale represents the count of DIB detections in each HEALPixel at different distances. Because QFs 0–5 are classified by {𝒟, σ_DIB, R_C}, the QF sky distribution would be related to DIB properties. At high latitudes, the stacked ISM spectra generally have small S/N due to the decreasing density of RVS objects (S/N is also affected by distance), which results in a large R_C. Moreover, 𝒟 is also expected to decrease with |b|. Thus, we can find numerous detections with QF = 1 and 3 at |b| ≳ 20° but rare with QF = 0 and 2. The detections with QF = 0 and 2 mainly occupy the Galactic middle plane where one expects more abundant ISM and stronger DIB signals (more RVS objects as well) and only extend to higher latitudes (|b| ~ 30°) towards the inner disc (|ℓ| < 30°) and the Galactic anti-centre (ℓ ~ 180°) where numerous molecular clouds exist. The two DIBs show similar QF distribution, but λ862.1 has a wider latitude distribution than λ864.8 for QF= 0 and 2 because 𝒟_864.8 is only ~35% of 𝒟_862.1 (see Sect. 5.4).

On the other hand, the QF sky distribution is similar between QF = 0 and QF = 2, although they contain different ranges of σ_DIB (the same for QF = 1 and QF = 3). Moreover, we can find more detections of λ862.1 with QF = 0 than QF = 2 and less with QF = 1 than QF = 3, which indicates that the fit of σ_DIB is heavily affected by 𝒟/R_C which can be treated as a measure of the S/N of DIB signal.

The detected DIB signals with QF = 4 are weak and noisy with a main distribution out of the middle plane. DIB detections with QF = 5 are complicated, containing the cases with λ_DIB out of the reasonable range, too small or too big σ_DIB, and too low 𝒟. QF = 5 is distributed in almost the full sky with |b| < 65° as it takes half of the DIB sample. It is interesting that the empty or low-density regions in QF = 4 and 5 are complementary to those in QF = 0.

QFs of 0 and 2 are recommended as the best level, and 1 and 3 are the secondary level. While 4 and 5 are the worst and are suggested to be better used in a statistical way.

Fig. 8 Distribution of the S/N of the stacked ISM spectra in the DIB depth-width plane (𝒟 vs. σDIB) for DIBs λ862.1 and λ864.8, respectively. The colour represents the mean S/N in each 0.001 × 0.003 nm bin for λ862.1 and 0.001 × 0.01 nm bin for λ864.8.

Fig. 9 Flowchart of the criteria to generate quality flags (QFs) for DIBs λ862.1 (upper panel) and λ864.8 (lower panel), respectively. The flag numbers and the corresponding arrow paths for classification are listed in the bottom box in each panel.

Fig. 10 Distribution of DIB quality flags (QFs) at each level (0–5, 0 is the best and 5 is the worst) for DIBs λ862.1 (upper panel) and λ864.8 (lower panel). The number of DIB detections and the fraction are indicated at the top of each bar.

4.2 DIB width

The width of the DIB profile contains vast information about the properties of the intervening ISM clouds and the DIB carriers, such as the profile broadening that is related to the gas kinetic and rotational temperature in different physical environments (e.g. Lai et al. 2020; Krelowski et al. 2021). The fitted σ_DIB in DIB-Spec is a measure of the average width under different ISM environments and may be affected by Doppler splitting or broadening, especially for distant voxels. An investigation of the Doppler effect and the decomposition of multiple velocity components with the DIB-Spec results will be done in a forthcoming work.

The distribution of σ_DIB of the two DIBs with the full DIB-Spec results (235 428 voxels) is shown in the left panels in Fig. 12 (entire sample). The histogram of the width of λ862.1 contains three peaks (upper left panel). The first peak at σ_862.1 = 0.02 nm indicates ISM spectra with low S/N or extremely weak DIB signals. The second peak at σ_862.1 = 0.12 nm, the initial guess of σ_862.1, is also caused by the fit to low-S/N spectra, where the fitted parameters are just around their initial values. The third peak at σ_862.1 = 0.22 nm, as we have seen in Fig. 8, marks the best fits of DIB λ862.1. For λ864.8 (lower left panel), the first (σ_864.8 = 0.1 nm) and second (σ_864.8 = 0.4 nm) peaks can also be found. But the width distribution is then very flat for larger σ_864.8, and the peak of σ_864.8 that indicates the best λ864.8 profiles is smoothed by the noisy detections, which could only be seen in the control sample (see below).

To determine the “best” and “secondary” ranges of σ_DIB for QF evaluation (Sect. 4.1), we apply a quality control of S /N > 100 and 𝒟 > 3R_C, and get 56 816 fit results for λ862.1 and 9927 for λ864.8, whose width distributions are shown in the right panels in Fig. 12. The peaks at the initial guess of σ_DIB disappear, and a quasi-Gaussian distribution of σ_DIB can now be found for both DIBs with a lower cut that excludes very narrow spurious features. The selected ranges of σ_DIB for QF determination are marked in Fig. 12 and discussed in Sect. 4.1. We note that the upper limit of the acceptable range for σ_864.8 in QF evaluation is enlarged because of the long tail of the distribution of σ_864.8 showing more detections with σ_864.8 ≳ 1.1 nm than with σ_864.8 ≲ 0.8 nm. The cause is unknown and DIB λ864.8 with a very broad profile needs further exploration.

The full width at half maximum (FWHM) of the two DIBs is calculated as ${\text{FWHM}}_{862.1}=\sqrt{8\ln \left(2\right)}\times {\sigma }_{862.1}$ for the Gaussian profile (Eq. (1)) of λ862.1 and FWHM_864.8 = 2 × σ_864.8 for the Lorentzian profile (Eq. (2)) of λ864.8. The joint distribution of FWHM_862.1 and FWHM_864.8 for 9778 detections with S/N > 100 and 𝒟 > 3R_C (for the two DIBs) is shown in Fig. 13. The distribution is Gaussian with a long tail for FWHM_864.8. The median FWHM of DIB λ862.1 is 0.52 ± 0.05 nm, which is consistent with Zhao et al. 2022, (0.55 ± 0.06 nm) but larger than previous results based on early-type stars, that is 4.3 Å of Herbig & Leka (1991), 4.38 Å of Jenniskens & Desert (1994), and 4.69 Å of Maíz Apellániz (2015). As discussed in Z22, the increase in FWHM_862.1 of our result could be explained by a Doppler broadening caused by our stacking strategy. Other effects, such as the observational instrument and the stellar residuals, may be proposed as well. With a machine-learning approach on published RVS spectra, Saydjari et al. (2023) got a σ_862.1 = 1.9 Å, corresponding to a FWHM of 4.47 Å. This consistency indicates that Gaia has no significant instrumental effects on the DIB measurement. Furthermore, a data-driven method could significantly reduce the influence of the residuals of specific stellar lines (e.g., Fe I for λ862.1) on the DIB width. On the other hand, Hobbs et al. (2009) and Fan et al. (2019) reported a FWHM of 3.56 Å and 3.98 Å for λ862.1, respectively, which are smaller than mentioned results above. Puspitarini et al. (2015) mentioned that many Gaia–ESO spectra of the λ862.1 region are contaminated by sky emission residuals which fall within the red wing of λ862.1. In principle, any emission residuals could make the DIB appear narrower than it actually is, even in the case of hot stars. Therefore, we propose that FWHM ~ 4.7 Å (σ_862.1 around 2 Å) could be proper for DIB λ862.1. But we should keep in mind that the DIB width could vary under different ISM environments. We note that DIB λ862.1 could contain multiple components (Jenniskens & Desert 1994) and present an asymmetric profile (e.g. Krełowski et al. 2019). But considering the resolution of RVS spectra, our study cannot reveal an accurate shape of λ862.1. Nevertheless, the high S/N of the stacked spectra would allow us to resolve different DIB velocity components in further analysis.

The median FWHM of DIB λ864.8 is 1.91 ± 0.44 nm, larger than 1.62 ± 0.33 nm in Z22, 1.4 nm in Herbig & Leka (1991) and 0.42 nm in Jenniskens & Desert (1994). Doppler broadening should have much less effect on the FWHM of λ864.8 due to its very large breadth. Z22 also measured the median σ_864.8 by stacking RVS spectra but in a much smaller sample (1103 detections) with weaker quality control. Figure 12 shows that involving noisy detections would reduce and smooth the peak of σ_864.8. The measurements in early studies could be questionable due to the broad span of λ864.8 and its superposition with some blended stellar lines. The median σ_DIB reported in this work is from a quasi-Gaussian distribution.

Fig. 11 Galactic spatial distribution of each level of DIB QF (0–5) for λ862.1 (left panels) and λ864.8 (right panels). The colour represents the counts of DIB detections in each HEALPixel at various selected distances. This HEALPix map has a level of 5, according to a spatial resolution of 1.8°.

Fig. 12 Distribution of σDIB for DIBs λ862.1 and λ864.8 (top and bottom, respectively). The dashed red and green lines correspond to the “best range” and the “secondary range” of σDIB defined in Sect. 4.1 for the DIB QF, respectively. The distribution of the full DIB-Spec results (235 428 voxels) is shown in the left panels, while the right panels show a quality-controlled sample with S/N > 100 and 𝒟 > 3RC (see Sect. 4.2). The latter criterion implies QF values of 0 or 2. The number of detected DIBs and the percentage after the quality control is indicated as well.

Fig. 13 Joint distribution of the FWHM of DIBs λ862.1 and λ864.8, measured in 9778 stacked ISM spectra, with S/N > 100 and 𝒟 > 3RC (see Sect. 4.2), generated by a Gaussian kernel density estimation (KDE). The white star indicates the peak density, and the red line in the central panel indicates the contour of the 2σ level. The orange lines indicate the median FWHM of the two DIBs.

4.3 Scaling factor for the lower and upper limits of DIB EW

The EW of the DIBs λ862.1 and λ864.8 (‘ew8620’ and ‘ew8648’ in Table 2) was calculated by Eqs. (5) and (6), respectively. Nevertheless, when estimating the lower and upper confidence levels of EW by the {𝒟, σ_DIB} pairs from MCMC samplings (see Sect. 3.3), the EW values of each pair were calculated by numerical integration of Eqs. (1) and (2) from 𝒟 − 3σ_DIB to 𝒟 + 3σ_DIB for the two DIBs, respectively, instead of using Eqs. (5) and (6). This approximation is suitable for a Gaussian function, as the integration equals to $\text{erf}\left(3\text{/}\sqrt{2}\right)\times \sqrt{2\pi }\times {𝒟}_{862.1}\times {\sigma }_{862.1}\approx 0.997{\text{EW}}_{862.1}$. But it is problematic for the Lorentzian profile of λ864.8 for the integration equals to 2arctan(3) × 𝒟_864.8 × σ_864.8 ≈ 0.795 EW_864.8. Therefore, we propose a scaling factor of 1.258 for ‘ew8648_lower’ and ‘ew8648_upper’ in Table 2, that is the table values times 1.258 are the correct lower and upper levels for EW_864.8. A scaling factor of 1.003 can also be used for ‘ew8620_lower’ and ‘ew8620_upper’.

Fig. 14 Number density of DIB detections in DIB-Spec as a function of λDIB and σDIB, without any cuts, for λ862.1 (upper panel) and λ864.8 (lower panel). The dashed blue lines indicate the range of permitted λDIB in the QF evaluation (see Sect. 4.1), and the dashed red and green lines correspond to the “best range” and the “secondary range” of σDIB, respectively.

4.4 Effect of stellar residuals

Because most of the target stars processed by DIB-Spec are late-type stars whose spectra contain abundant stellar components, the stellar residuals from the unresolved stellar features or the not well-modelled stellar lines would affect the DIB detection and measurement, such as the shift of λ_DIB or the broadening of σ_DIB. In extreme cases, the detected signal is an artefact that comes from the stellar residuals. As pointed out by Saydjari et al. (2023), for example, the DIB detections in the Gaia DR3 catalogue would be contaminated by not perfectly modelled lines, such as Fe I, when σ_DIB is very small.

Figure 14 shows the distribution of DIB detections as a function of λ_DIB (in stellar frame) and σ_DIB, which can straightforwardly identify the effect of stellar residuals. A vertical stripe at the initial guess of σ_DIB can be seen for both the two DIBs that corresponds to the noisy cases, mostly with 𝒟 < R_C. λ_862.1 is uniform and widespread when σ_862.1 < 0.04 nm, where DIB detection is noise dominated. For larger σ_862.1, especially 0.16–0.28 nm, λ_862.1 is concentrated between 862.2 and 862.5 nm, representing the reliable measurements. Our detection of λ862.1 is not or at most very weakly affected by the stellar residuals, as no clustering near the stellar lines is found in the σ_862.1 − λ_862.1 plane. This is because the stacking increases the S/N of the ISM spectra and averages out the large residuals in individual spectra.

The very broad profile of λ864.8 makes it difficult for accurate measurements of λ_864.8 and σ_864.8, resulting in a very scattered λ_864.8 versus σ_864.8 distribution. However, the large σ_864.8 can hardly be affected by the stellar residuals.

4.5 Uncertainties of DIB parameters

Since 𝒟 and σ_DIB are correlated with each other during the MCMC fitting, the distribution of the uncertainty of the DIB parameters presents a dependency on both 𝒟 and σ_DIB. Figure 15 shows the distribution of the fractional uncertainties of 𝒟, σ_DIB, and EW in the 𝒟 − σ_DIB plane for λ862.1 (left panels) and λ864.8 (right panels). The uncertainty of 𝒟, σ_DIB, and EW are taken as the mean difference between their lower, median, and upper values.

Large uncertainties can generally be found in regions with small 𝒟 or σ_DIB. DIB detections with very small 𝒟 will get large QF (low reliability) by comparing with R_C, and those with small σ_DIB are ruled out by the width border defined in the QF evaluation. It is noted that detections with small uncertainties can be found very close to the upper limit of σ_DIB in priors, 0.5 nm for λ862.1 and 1.5 nm for λ864.8. These small uncertainties are because of the narrow sampling range in the MCMC fitting and do not represent good measurements. Figure 12 also shows that the number of detections with extremely large σ_DIB will significantly decrease when applying strict constraints. For λ862.1, detections with low uncertainty (<10%) are located in similar fields for 𝒟_862.1, σ_862.1, and EW_862.1 and correspond to the “best range” of σ_862.1 defined in Sect. 4.1. On the other hand, low-uncertainty regions for λ864.8 do not present a well-constrained range for σ_864.8.

The uncertainty of EW in some regions, mainly the bottom in each panel with small σ_DIB and large 𝒟, is much smaller than that of σ_DIB, although EW was calculated by 𝒟 and σ_DIB. The reason is that in each MCMC chain, in spite of the width with large uncertainty, the depth could be stable with small uncertainty. Consequently, the calculated EW also concentrates and has smaller relative uncertainty. Besides the method used in DIB-Spec, the EW uncertainty is also proposed to be calculated by the span of the profile (3×FWHM), the pixel resolution (δλ = 0.03 nm/pixel), and the noise level of the line centre (R_C = std(data − model)), as ${\sigma }_{\text{EW}}=\sqrt{6\text{\hspace{0.17em}FWHM}\delta \lambda }\times R_C$ (Z22). Similar formulas were also given by Vos et al. (2011) and Vollmann & Eversberg (2006). Figure 16 shows the difference between these two kinds of EW uncertainty. The difference, in general, is small (over 90% within 0.02 Å for λ862.1 and within 0.07 Å for λ864.8), but our method systematically gives out smaller σ_EW for λ862.1 and larger σ_EW for λ864.8.

5 Validation tests

We present in this section a number of validation tests for the DIB results from DIB-Spec, including the comparison with Gaia DR3 DIB catalogue for the measurements of EW_862.1, the EW maps of the two DIBs in a Galactic view at different distances, and the correlation between two DIBs and dust reddening.

5.1 Comparison between fitted and integrated DIB EW

Differing from the direct measurement of DIB EW by integrating the ISM spectrum (e.g. Hobbs et al. 2008; Fan et al. 2019), the DIB EW in this work was calculated by the analytic function. The integrated EW is not affected by the asymmetry of the DIB profile, while the fitted EW is less affected by the noise and the overlapping of different DIB profiles. Specifically, the EW of λ862.1 cannot be directed integrated because the profiles of λ862.1 and λ864.8 could be overlapped. Therefore, we first subtracted the fitted profile of λ864.8 from the ISM spectrum and then integrated the rest part of the ISM spectrum within λ_862.1 ± 3σ_862.1 to calculate the integrated EW of λ862.1. Figure 17 shows the comparison between fitted and integrated EW_862.1 for 6240 cases with S / N ⩾ 100, QF_862.1 = 0, and QF_864.8 = 0. We did not do such a comparison for λ864.8 as the spectra suffer there the residuals of the CaII line. The fitted and integrated EW_862.1 are highly consistent with each other with a mean difference of only 0.001 Å and a standard deviation of 0.020 Å. This proves that the possible asymmetry of the DIB profile has little effect on EW_862.1. After the check of the ISM spectra, we find that the outliers with much larger fitted EW_862.1 than integrated EW_862.1 are caused by the wrongly fitted central position of DIB λ862.1. The λ_862.1 of these cases are close to 863 nm which is the upper limit of λ_862.1 during the MCMC fitting. A larger upper limit could improve the results of these fittings.

Fig. 15 Distributions of the fractional uncertainties of 𝒟, σDIB, and DIB EW as a function of 𝒟 and σDIB for λ862.1 (left panels) and λ864.8 (right panels). Colour represents the mean fractional errors in each 0.001 × 0.003 nm bin for λ862.1 and 0.001 × 0.01 nm bin for λ864.8.

Fig. 16 Histogram of the difference between the EW uncertainty estimated in this work ${\sigma }_{{\text{EW}}_1}$ and the one calculated by the formula ${\sigma }_{{\text{EW}}_2}=\sqrt{6\text{\hspace{0.17em}FWHM}\delta \lambda }\times R_C$ (see Sect. 4.5) for λ862.1 (upper panel) and λ864.8 (lower panel), respectively. The mean (Δ) and standard deviation (σ) of the differences are also indicated.

Fig. 17 Upper panel: comparison between the fitted and integrated EW862.1. The colour represents the number density (estimated by a Gaussian KDE) of the data points in a linear scale. The grey colour bars show the uncertainty of fitted EW862.1. The dashed red line traces the one-to-one correspondence. A zoom-in panel shows the distribution of the EW difference (ΔEW = EWfit − EWint). The mean (Δ) and standard deviation (σ) of the EW difference are indicated. Lower panel: the distribution EW difference as a function of fitted EW862.1.

Fig. 18 Upper panel: comparison between the EW862.1 measured by DIB-Spec and the mean EW862.1 in each voxel taken from the DIB catalogue in DR3 for 8963 voxels. The colour represents the number density (estimated by a Gaussian KDE) of the data points in a linear scale. The grey colour bars show the EW uncertainty in DIB-Spec and the standard deviation of EW in DR3. The dashed red line traces the one-to-one correspondence. A zoom-in panel shows the distribution of the EW difference (ΔEW = EWFPR − EWDR3). The mean (Δ) and standard deviation (σ) of the EW difference are indicated. Lower panel: the variation of the mean EW difference in each EW bin with a step of 0.01 Å with EW862.1 in FPR. The orange shades the range of 1σ.

5.2 Comparison with DIB measurements in Gaia DR3

To make a direct comparison of the measurements of EW_862.1 between FPR and DR3, the mean EW_862.1 of the full DIB catalogue in DR3 is calculated in each voxel defined in DIB-Spec. Figure 18 presents the comparison between EW_862.1 from FPR and mean EW_862.1 from DR3 for 8963 voxels (3.81% of the amount in FPR) that contain at least ten DIB detections in DR3. The difference of EW_862.1 (FPR-DR3) presents a quasi-Gaussian distribution with an extremely small mean value of 0.002 Å and a standard deviation of 0.025 Å. The EW difference is very close to zero between 0.06 and 0.2 Å of EW_862.1 from FPR. For smaller realms, the DIB results in DR3 present larger values of EW than those in FPR, and their mean difference increases with the decreasing EW. This is because in DR3, we only successfully detected relatively strong DIB signals limited by the S/N of individual RVS spectra, whereas the EW measured in FPR represents an average DIB strength in each voxel by stacking a set of ISM spectra. This difference, between the strong signals captured in DR3 and the mean strength measured in FPR, also implies a variation in DIB strength in the solar neighbourhoods. On the contrary, FPR gets a systematically larger EW_862.1 than DR3 when EW_862.1 ≳ 0.2 Å, because large EW generally comes from distant voxels (or dense clouds) where detections in FPR and DR3 may trace different ISM environments. The DIB detections made in individual RVS spectra are mainly located in diffuse and intermediate regions, while DIB signals in denser regions can be measured in stacked spectra with higher S/N. Moreover, λ864.8 may also contribute to the difference as in DR3 λ864.8 was not considered for fittings. Thus, if the profile of λ864.8 is treated as the continuum placement for normalization (the two DIBs are close to each other, see Fig. 7 for example), EW_862.1 would be underestimated. This effect will be investigated in detail in follow-up work.

To compare the spatial distribution of EW_862.1 between FPR and DR3, Fig. 19 shows the EW maps in the Galactic (XY), meridian (XZ), and rotational (YZ) planes where the Sun is located at the origin with the GC as the primary direction. The two EW maps (FPR on the left and DR3 on the right) were constructed in different ways. For DR3, we make use of the high-quality sample (see its definition in S23) with an additional constraint of 0.18 ⩽ σ_DIB ⩽0.26 nm (the best range of σ_862.1 determined in DIB-Spec), resulting in 52 180 DIB detections. Then the median EW_862.1 was taken in 0.1 × 0.1 kpc bins for the Galactic, meridian, and rotational planes, respectively, in the Cartesian system using detections within ±50 pc above and below each corresponding plane. For FPR, EW_862.1 in Galactic, meridian, and rotational planes were taken from the voxels that are crossed with each corresponding plane, and the crossed sections were painted by the EW862.1 in that voxel. To have a clean map with reliable detections, we require S / N ⩾ 100 and QF_862.1 ⩽ 2.

Similar large-scale structures can be found in each plane between FPR and DR3 results, while λ862.1 in FPR can be detected in more distant regions, such as between the Perseus Arm and the Outer Arm and beyond the Scutum-Centaurus Arm. Two sightlines, ℓ ~ −70° and ℓ ~ −117°, have significant low EW_862.1 reaching over 4 kpc, indicating two void regions with less abundant ISM between the Galactic main arms. On the other hand, much fewer DIB signals can be detected towards ~80° − ~90° and − ~ 100° − ~ − 90° for both FPR and DR3 because in these two ranges, we are looking to directions that are parallel to the Local Arm with less intervening DIB clouds.

5.3 Equivalent width (EW) map: In a Galactic view

Figure 20 shows the Galactic distribution of EW_862.1 (left panels) and EW_864.8 (right panels) at four distances, d_c = 1.05, 1.53, 2.12, and 3.23 kpc from the DIB results in FPR. As DIB-Spec stacked spectra in (α, δ, d) voxels, the distances indicated in the figure are the centre of the voxels from the Sun. S / N ⩾ 100 and QF ⩽ 3 (for λ862.1 and λ864.8, respectively) are used to control the quality of the DIB detections. The number of selected DIB detection decreases with distance for both λ862.1 and λ864.8 (from d_c = 1.05 kpc to d_c = 1.53 kpc for λ864.8 is an exception). And at 2.12 kpc, there are 4316 detections of λ862.1 that cover 35% of the full sky. The amount of λ864.8 is around 70% of λ862.1 at each d_c.

The sky coverage of λ862.1 with reliable detections in FPR is similar to that in DR3 (see Fig. 6 in S23 for a median EW map at HEALPix level 5), but the distance-sliced maps in Fig. 20 contain more details. At d_c = 1.05 kpc and 1.53 kpc, some clumpy regions with large EW are consistent with nearby molecular clouds, such as Phoenix (ℓ ~ 30°–40°) and Cygnus complex (ℓ ~ 80°) in the middle plane and Ophiuchus (ℓ ~ 0° and b ~ 15°) at high latitude, although the resolution of the present DIB EW map is pretty low (1.8°). Nevertheless, some clouds seem to disappear in the DIB EW map, like the Cepheus and Polaris Flare (ℓ ~ 110°–120°), which can be clearly seen in the Gaia Total Galactic Extinction map (see Fig. 24 in Delchambre et al. 2023) and the Planck dust map (see Fig. 3 in Planck Collaboration Int. XLVIII 2016). Cepheus and Polaris Flare are nearby clouds (<400 pc) with strong CO emission, but no strong DIB signals were detected in this region in either this FPR or DR3, in spite of many RVS observations there. With the increase of distance, at 3.23 kpc, the high-EW regions are linked up into a bright stripe along the galactic plane.

At low galactic latitudes, λ862.1 was detected in most longitudinal directions, but a gap can be found at ℓ = − 120°, where the CO emission is also weak (see Dame et al. 2001). It seems that λ864.8 is more concentrated around the regions with larger EW, which should be a bias due to the difficulty in measuring weak λ864.8 with the present spectral S/N level. On the other hand, λ862.1 and λ864.8 seem to occupy similar latitude ranges, that they are mainly distributed within ±30° and this range decreases with distance. In principle, the DIB signals detected in nearby voxels should also be seen in distant voxels along the same direction as EW is an integrated variable of the abundance of DIB carriers. But in fact, the number of DIB detections about ±30° of latitude decreases with the distance for the control sample shown in Fig. 20. The reason could be that some of the DIB fittings in distant voxels were filtered out by the criteria of the control sample due to the low quality of their ISM spectra. Furthermore, the measured DIB strength in distant voxels would be a biased mean, as faint stars among or behind dense clouds are hard to be observed by Gaia. Therefore, it is possible to see the decrease of DIB EW with distance even if the voxels all contain reliable DIB detections. This observational bias could also happen in nearby regions behind dense clouds.

Fig. 19 Distribution of EW862.1 in the Galactic plane (XY), meridian plane (XZ), and rotational plane (YZ) for the FPR results (left panel) and the DR3 results (right panel), respectively, plotted over the Milky Way sketch created by Robert Hurt and Robert Benjamin (Churchwell et al. 2009). Some log-periodic spiral arms described in Reid et al. (2019) are also presented by coloured lines: Scutum-Centaurus Arm in orange; Sagittarius-Carina Arm in purple; Local Arm in black; Perseus Arm in green; Outer Arm in cyan; and the spur between the Local and Sagittarius-Carina arms in blue. The Galactic centre is located at (X, Y, Z) = (8.15, 0, 0).

Fig. 20 Galactic distributions of EW862.1 (left panels) and EW864.8 (right panels) from DIB results in FPR in Mollweide projection at HEALPix level 5, at four distances. The distance (dc) and the number of voxels (N) are indicated in each subpanel.

5.4 Correlations between DIBs and dust reddening

The linear correlation between DIB strength and dust reddening is a general property for many strong DIBs (e.g. Friedman et al. 2011; Lan et al. 2015) and could be treated as a validation test of the DIB measurement. For this comparison, we only select DIB detections with QF_862.1 = 0 and QF_864.8 = 0, resulting in 6278 cases (2.67% of total) to increase the reliability and to focus on the Galactic middle plane (see QF distribution in Fig. 10) where the interstellar materials are generally well mixed with each other, so a tight linear correlation can be expected between λ862.1, λ864.8, and dust. There are 3 677 773 (60%) RVS objects that have E(BP − RP) from GSP-Phot (Andrae et al. 2023). The mean E(BP − RP) in a voxel was calculated with all the stars with E(BP − RP) in that voxel. Thus, the number of used stars in a voxel for mean E(BP - RP) and DIB fitting could be different. There are 6038 (2.56%) voxels containing at least ten stars having E(BP − RP) and QF = 0 for the fittings of both λ862.1 and λ864.8, which were used to compare DIB EW and E(BP − RP) (see the third and fourth panels from top to bottom in Fig. 21). On the other hand, for the comparison of EW and depth between λ862.1 and λ864.8 (see the first and second panels from top to bottom in Fig. 21), we only required QF = 0 and got 6278 (2.67%) voxels. The requirement of the highest level of QF is the main constraint to the sample size. The comparison including other QFs and discussions will be made in a follow-up work. The linear fit to each correlation was achieved by an ordinary least squares regression using the Python package statsmodels, with some upper limits on the variables (indicated in Fig. 21) to exclude the nonlinear realms. The fitted slope (α), intercept (t), and Pearson coefficient (r_p) are marked in each panel in Fig. 21 as well.

Tight linearity can be found between EW_862.1 and EW_864.8, between 𝒟_862.1 and 𝒟_864.8, and between EW_862.1 and E(BP − RP), and a weaker one for EW_864.8 and E(BP − RP). In the correlation of EW for the two DIBs, some cases present greater EW_864.8 than expected when EW_862.1 ≳ 0.5 Å. Similar deviation from the linearity could also be found for their depth. This deviation could indicate the departure of the carriers of the two DIBs (if we assume they have different origins) as the trend is gradual and continuous with EW and the ISM spectra of these cases show prominent DIB features. No apparent deviation is seen between EW and E(BP − RP), but the linear fit needs to be limited to EW_864.8 < 1 Å to get a small intercept between EW_864.8 and E(BP − RP).

Linear fits done in S23 and Z22 are also shown for comparison. We note that Z22 fixed the intercept as zero (except 𝒟_862.1 and 𝒟_864.8) and applied 2σ-clipping for the linear fits due to its small sample size (1103 DIB detections). The E(BP − RP)/EW_862.1 ratios from FPR and S23 are consistent with each other with a difference of 9.20%, but the degree of dispersion is much lower in FPR because of stronger quality control and the higher S/N in general after stacking RVS spectra. As a reference, the difference of E(B − V) /EW_862.1 between literature studies varies from 4% to 41% with a mean of 20%⁵ The comparison between the DIB results in FPR and Z22 is more meaningful because they both use BNM to build the stellar templates and fit DIB profiles in stacked RVS spectra with the same model, only with differences in the sample size and stacking strategy. By an internal comparison, a consistent tendency can be found for each correlation between this FPR and Z22, with a similar degree of dispersion of the scatter plots, but the DIB results in FPR contain much stronger DIB signals. For example, EW_862.1 is mainly within 0.15 Å in Z22 and 0.4 Å in this FPR. Both the slope and intercept are highly consistent with each other for 𝒟_862.1 and 𝒟_864.8 in this FPR and Z22, with a difference of only 6.50%. The difference in the correlation between EW_862.1 and EW_864.8 is slightly larger (~ 10%), and the fit in FPR has a positive and non-negligible intercept. The cause is unclear. Strict QF control in FPR, which leads to the lack of λ864.8 with EW_864.8 < 0.1 Å, could have an impact but cannot fully explain the offset as the tight linear correlation keeps until EW_862.1 ~ 0.5 Å. Including all the other QFs (1–5) could reduce the intercept from 0.084 to 0.076, but the degree of dispersion will significantly increase. The biggest difference between FPR and Z22 occurs for E(BP − RP)/EW_862.1 (26.93%), which is caused by the sigma-clipping in Z22. We refit E(BP − RP)/EW_862.1 using the data in Z22 without sigma-clipping and get a ratio of 4.468 ± 0.081. The difference then becomes much smaller (9.50%). Additionally, the 2σ-clipping caused the linear fit in Z22 to be dominated by DIBs with EW_862.1 between 0.04 and 0.1 Å, where much fewer cases were seen in this FPR after QF control. Thus, the E(BP − RP)/EW_862.1 ratio would also vary in different EW_862.1 ranges. The sigma-clipping, on the other hand, has a very light effect on the correlation between E(BP − RP) and EW_864.8. The difference is only 3.77%, although in Z22 the intercept was fixed as zero and in FPR we got a positive one of 0.014. Besides the sample size and sigma-clipping, many other factors can affect the fitting results as well, such as the QF control and the source of dust reddening. More detailed discussions are beyond the scope of this work. And the deviation and variation in the correlation between different ISM, like DIB and dust, need to be understood by considering the interstellar environment as well.

Fig. 21 Diverse correlations between λ862.1, λ864.8, and mean E(BP − RP) from GSP-Phot (Andrae et al. 2023) in each voxel. The data points are coloured by their number densities estimated by a Gaussian KDE. The red lines are the linear fit to the data points. The dashed green lines indicate the upper limits on the variables for the linear fits. And the fitted slope (α) and intercept (t) are marked in each panel, together with the Pearson correlation coefficient (rp) and the number of DIB detections (N). Dashed black lines are previous results from Z22, and the dashed blue line is from S23.

6 The Local Bubble

It is widely known that the Local Bubble has a much lower density than the average of the ISM in the solar neighbourhood because of its harsh environment (high temperature and low density; Welsh et al. 2010; Lallement et al. 2014), but Farhang et al. (2019) reported the detections of DIBs λ578.0 and λ579.7 in the Local Bubble in the spectra of 359 early-type stars and mapped the 3D density distribution of their carriers. S23 also found some relatively strong signals of λ862.1 that are very close to the Sun and generally suggested that λ862.1 could also be detected in the Local Bubble. However, by reanalysing the public RVS spectra (about one million) in DR3, Saydjari et al. (2023) claimed no detection of λ862.1 with high confidence levels in their analysis within the Local Bubble. As S23 only presented the median EW_862.1 distribution in the Galactic plane and Saydjari et al. (2023) only analysed a small part of the RVS spectra, we plan to make a thorough investigation to see if we can reliably find λ862.1 in the Local Bubble.

As a preliminary investigation, we focus in this work only on 145 DIB detections with a high level of S/N and QF, that is S/N ⩾ 300, QF_862.1 ⩽ 2, QF_864.8 ⩽ 2, and d_c < 300 pc. Then their ISM spectra were further visually inspected. As DIB signals within the Local Bubble would be very weak, a highly reliable DIB detection needs to satisfy that the flux uncertainty of its ISM spectrum within the DIB profile is smaller than the depth of each wavelength bin. Finally, we found four DIB detections whose voxels are mostly inside the Local Bubble and 12 detections whose voxels are crossed with the surface of the Local Bubble determined by Pelgrims et al. (2020) using the 3D dust map of Lallement et al. (2019). The DIB parameters of λ862.1 and field information of these 12 voxels are listed in Table 4. And the projection of these 12 voxels and the surface of the Local Bubble in the Galactic (XY), meridian (XZ), and rotational (YZ) planes are shown in the upper panels in Fig. 22, and their ISM spectra and DIB fittings are presented in the lower panel. The profile of the DIB λ862.1 is conspicuous in the noisy spectra, with all QF_862.1 ⩽ 2, although the DIB fitting is not good, especially for λ864.8 (whose profile is too shallow for the S/N of the spectra). The Local Bubble is known to contain assembling molecular clouds on its wall so that the dust abundance is significantly different inside and on the surface of the Local Bubble. Besides the selection bias (weaker signals cannot be detected by present data), the DIB profile would be broadened by the heavy noise, or the continuum placement was affected by the poorly fitted λ864.8 profile, which leads to an overestimation of EW. The four internal voxels mostly inside the Local Bubble, containing apparent DIB profiles, seem to indicate the possible detection of λ862.1 in the Local Bubble, but further investigation is necessary, especially taking into account all available ISM spectra.

Fig. 22 Upper panels: spatial distribution of 12 reliable DIB detections by visual inspection. Their voxels are projected into the Galactic (XY), meridian (XZ), and rotational (YZ) planes. The reds are inside the Local Bubble and the blues are crossed with the surface of the Local Bubble. The grey marks the surface of the Local Bubble determined by Pelgrims et al. (2020). Lower panel: black lines are stacked ISM spectra, and the red/blue lines are the DIB fitting results of the corresponding DIB profiles. The orange indicates the masked spectral region between 866 and 866.8 nm in the fitting. The vertical dashed magenta line indicates the rest-frame wavelength of DIB λ862.1 of 862.323 nm determined in Gaia Collaboration (2023).

7 Caveats and known issues

We list the caveats and the known issues of the first results of DIB-Spec as presented in this FPR. Future developments for DR4 aim to tackle and possibly remove these issues.

1. The inverse of parallax (1/ϖ) of the background stars was directly used in DIB-Spec for stacking spectra in different voxels. But 1 / ϖ would not be the true distance of the stars in distant zones (several kpc, Bailer-Jones 2015). The detected DIB signal accounts for an integration of the DIB carriers between the background stars and us. If the distance between the DIB carriers and the background stars is much larger than the distance uncertainty of 1 / ϖ, the DIB measurement would be safe. Otherwise, additional uncertainty of the DIB detection will be introduced. Distant voxels or voxels with few targets will get heavier influence.

2. About 5% (12 692) of stacked ISM spectra have zero flux uncertainty at each wavelength bin. Some of them are due to the zero flux error in observed RVS spectra used for stacking. But the cause for others is still unknown. The flux uncertainties of these spectra were fixed as 0.01 during the DIB fitting.

3. The flux uncertainty at each wavelength bin was taken as the standard error in the mean (SEM) from the individual RVS spectra (see Sect. 3.3). Nevertheless, as the median of individual fluxes was taken for stacking, the flux uncertainty should be 1.253 × SEM. Presently used values in the spectra table underestimate the uncertainties of the stacked ISM spectra.

4. Negative 𝒟 (‘p08620’ and ‘p08648’ in Table 2) can be found for λ862.1 (166) and λ864.8 (198), although the prior of 𝒟 is always positive. These cases are all badly fitted due to the low-S/N spectra and/or the weak DIB signals.

5. A scaling factor is suggested to be used to correct the lower and upper confidence levels of EW for λ862.1 and λ864.8 (‘ew8620_lower’, ‘ew8620_upper’, ‘ew8648_lower’, ‘ew8648_upper’ in Table 2), respectively (see Sect. 4.3 for details).

6. There are a set of detections reported in the output table ‘interstellar_medium_params’ with inconsistencies in their EW and the lower or upper confidence levels of EW, even after the correction by the scaling factors (see Sect. 4.3). Specifically, there are 421 cases with EW_862.1 < EW_862.1,lower, 351 cases with EW_862.1 < EW_862.1,upper, 2791 cases with EW_864.8 < EW_864.8,lower, and 440 cases with EW_864.8 < EW_864.8,upper. The cause of this inconsistency could be some problems in recording the lower and upper confidence levels when producing the DIB results. We note that this inconsistency does not mean a bad DIB fitting, but the EW confidence levels would be problematic for such cases.

7. When we compare the target and reference spectra to build the stellar template, we have to reduce the weights of the Ca II lines to have a better model of weaker lines such as Fe I. Thus, we cannot model Ca II lines very well and mask the specific region (866.0–866.8 nm) during the fitting. Such a region falls within the profile wing of DIB λ864.8, introducing an uncertainty on the determination of the precise boundary of the DIB. We will try other data-driven methods in the future to model all the stellar lines without downweighting the Ca II lines.

8 Summary and conclusions

We summarize here the processing and validation of two published tables produced by the DIB-Spec module in the Gaia FPR, one for fitted DIB parameters (Table 2, ‘interstel-lar_medium_params’) and the other for stacked ISM spectra (Table 3, ‘interstellar_medium_spectra’). DIB-Spec derived ISM spectra for 5 983 289 RVS objects with |b| < 65° as targets using the other 160 392 RVS spectra (|b| ⩾ 65°) as references and the BNM. The individual ISM spectra were stacked to increase the spectral S/N in defined 3D voxels with a resolution of 1.8° (level 5 HEALPix binning) on the celestial sphere and 0.07–1 kpc in distance, and on average 20 spectra in each voxel. DIB-Spec fitted and measured the two DIBs in 235428 voxels, with a DIB model applying a Gaussian profile for λ862.1 and a Lorentzian profile for λ864.8. A median FWHM was determined as 0.52 ± 0.05 nm for λ862.1 and 1.91 ± 0.44 nm for λ864.8, which can be referred to as a typical value considering a large space coverage. Users are encouraged to use QF to control the quality of the DIB fittings for specific investigations and the related discussions in Sect. 4 are also useful.

Taking advantage of the stacking procedure, DIB-Spec extends the DIB detection in distance compared to the DIB catalogue in Gaia DR3 (Gaia Collaboration 2023). The DIB strength of λ862.1 is highly consistent between DR3 and this FPR (a mean difference of 0.002 Å with a standard deviation of 0.025 Å), and some systematic differences along with EW_862.1 would be caused by the selection bias between the two samples. We provide several other validation tests as well: The DIB EW map in a Galactic view at different distances shows the integration of DIB carriers along the sightlines, revealing some clumpy dense regions corresponding to notable molecular clouds, such as Phoenix, Cygnus, and Ophiuchus, while Cepheus and Polaris Flare were not seen in either of the DIB EW maps. Based on a high-quality subsample with QF = 0 for both λ862.1 and λ864.8 (6278 detections, 2.67% of total), linear correlations between λ862.1, λ864.8, and dust reddening (E(BP − RP)) were found to be consistent with those in Gaia Collaboration (2023) and Zhao et al. (2022), with the smallest difference of 3.77% and biggest one of 26.93%. We also found some detections of λ862.1 inside and around the surface of the Local Bubble with prominent profiles in the derived ISM spectra.

The DIB work in this FPR is not only complementary to Gaia DR3, but a pathfinder for future releases. The acquired experience and caveats will benefit the development of the DIB-Spec module for future Gaia releases. DIB results in this FPR have already shown the power of using numerous RVS spectra to map both the intermediate and strong DIB λ862.1 and the broad and shallow DIB λ864.8 in the solar neighbourhood and reaching over 4 kpc. In particular, we note that broad DIBs such as λ864.8 were never measured in a sample of hundreds of thousands of spectra before.

Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/Gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/Gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. The full acknowledgements are available in Appendix D.

Appendix A Error analysis of the BNM

In this section, we estimate the magnitudes of the errors introduced by the BNM, which was used in DIB-Spec to build the stellar templates for target spectra (Sect. 3.2). For each spectrum in the reference sample (160392, see Sect. 3.1), we apply the BNM to generate its stellar template using all the other reference spectra. In total, 159 591 (99.5%) spectra have enough best neighbours (>10) to build their stellar templates. The distribution of the flux residuals between the observed RVS spectra and the generated stellar templates (observed – modelled) is shown on the left-hand side in Fig. A.1, and four generated stellar templates are shown as examples in the right panels together with the corresponding observed RVS spectra.

Fig. A.1 Left panel: Distribution of flux residuals (observed – modelled) with spectral wavelength for 159 591 reference spectra. Each row presents one spectrum and the spectra are sorted by their S/N. The dashed blue and red lines indicate the position of S/N equalling 50 and 100, respectively. Some typical stellar lines within the RVS spectral region determined by Contursi et al. (2021) are indicated as well. The lower panel is a zoom-in plot of the upper one to show the distribution of the residuals in the DIB window (860–868 nm). Right panel: Four examples of the reference spectra (black lines with observed flux errors) and their derived stellar templates (red lines). The orange lines are the flux residuals (observed – modelled) with dashed black lines indicating ±5% of the continuum.

The flux residuals in the spectra present some patterns, especially in low-S/N spectra, instead of being uniform or noise-dominated. Positive residuals (red region) can be found near the Ca II triplets because of their reduced weights and some other stellar lines like Si I at 8658.4 Å, Fe I at 8623.97 Å and 8677.13 Å, indicating an overestimation of the depth of these stellar lines. On the other hand, negative residuals (blue region) mainly appear in spectra with S/N < 50, which could be due to a failed modelling of the lines by BNM in low-S/N spectra or to the improper normalization. The positions of some strong stellar lines are indicated in Fig. A.1. These lines are determined by Contursi et al. (2021) for the RVS spectra.

To quantitatively characterize the uncertainty introduced by BNM into the ISM spectra and consequently to the DIB measurement, we calculate the mean absolute residuals (MAR = mean(|observed – modelled|)) of normalized flux between 861.2 and 866.0 nm (region of the two DIBs in priors). The relationship between MAR and spectral S/N is shown in Fig. A.2. MAR has a strong dependence on spectral S/N. We applied a linear fit to them in logarithmic scale with 3σ clipping, and got a relation of log₁₀(MAR) = −0.16 − 0.99 × log₁₀(S/N) (the dashed black line in Fig. A.2) that can be used to describe the detections in the main branch, that is MAR linearly increases with 1/(S/N) for most of the stars, although MAR would be slightly larger than expected for very large S/N. Other 11.1% of stars in a secondary branch have larger MAR, which could be caused by the bad normalization of RVS spectra (most with S/N < 50) and/or the low number density in the vicinity of the queried stars. The reason for the dependence of MAR on the spectral S/N could be that high-S/N spectra are less affected by the random noise and have best neighbours with higher S/N, as BNM rejects the reference spectra with morphological differences larger than 3/(S/N) of the queried spectrum.

Fig. A.2 Variation of the mean absolute residual (MAR) between observed and modelled RVS spectra, calculated within the DIB window (861.2−866.0 nm), with the spectral S/N. The dashed green lines indicate S/N = 50 and MAR = 0.02, respectively. The dashed black line is fitted to MAR and S/N on a logarithmic scale.

From our test, we conclude that the uncertainty introduced by BNM strongly depends on the S/N of the queried spectra. When S/N > 50, MAR, the average magnitude of the flux residuals in the DIB vicinity (861.2−866.0 nm), is smaller than 0.02 (2% of the continuum) for 96.8% of spectra in the reference sample, and smaller than 0.01 for 61.7% of spectra. We expect BNM to have a similar performance on the RVS target sample because the target and reference samples have similar S/N distribution (see the right panel in Fig. 5).

Appendix B Corner plots of the DIB fittings

Figures B.1–B.5 show the corner plots of the DIB fitting in the voxels according to the examples shown in Fig. 7. The histograms and scatter plots show the one- and two-dimensional projections of the posterior distributions of the fitted parameters, with red squares and lines indicating the best estimates. Because the intermediate quantities of DIB-Spec were dropped, the posterior distributions of the parameters were drawn by refitting the stacked ISM spectra in the same way as DIB-Spec, some tiny differences can thus be found between the best estimates (red) and the fitting results in the output table of DIB-Spec (blue).

Fig. B.1 Corner plot of the DIB fitting in the voxel with Ipix = 10450 and dc = 1.05 kpc (the first panel in Fig. 7 from top to bottom). The histograms and scatter plots show the one- and two-dimensional projections of the posterior distributions of the fitted parameters. The red squares and lines indicate the best-fit estimates for each parameter in the reproduced fitting. And the dashed blue lines mark the fitted parameters in the output table of DIB-Spec.

Fig. B.2 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 1.27 kpc (the second panel in Fig. 7 from top to bottom)

Fig. B.3 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 2.12 kpc (the third panel in Fig. 7 from top to bottom)

Fig. B.4 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 2.40 kpc (the forth panel in Fig. 7 from top to bottom)

Fig. B.5 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 3.23 kpc (the fifth panel in Fig. 7 from top to bottom)

Appendix C Python script for converting the spectra table

The python script below shows a simple method to convert the spectra table to a fits file, in which each row stands for one stacked ISM spectra.

import numpy as np

import pandas as pd

from tqdm import tqdm

from astropy.io import fits

dibs = pd . read_csv (somepath+ ’parameter_table.csv’)
tabe = pd . read_csv (somepath+ ’spectral_table.csv’)

spec_lc = tabe.lc.values

spec_bc = tabe.bc.values

spec_dc = tabe.dc.values

cata_lc = dibs.lc.values

cata_bc = dibs.bc.values

cata_dc = dibs.dc.values

wave = np.unique (tabe [’lambda’].values)

flux = np.zeros((dibs.shape[0], wave.shape[0]),dtype=‘float’)

ferr = np.zeros_like (flux)

for i in tqdm(range(dibs.shape [0])):

   t = (spec_lc == cata_lc [i]) & (s pec_bc == cata_bc [i]) & (spec_dc == cata_dc [i])

   flux [i] = data [’flux’] [t]

   ferr[i] = data [’flux_uncertainty’] [t]

hdr0 = fits.Header ()

hdr0 [‘COMMENT1’] = ’stacked⊔SM⊔spectra’

hdr0 [‘EXTNAME’] = (’flux’, ‘normalized⊔flux’)

hdr1 = fits.Header()

hdr1 [‘EXTNAME’] = (’ferr’, ’flux⊔uncertainty’)

hdr2 = fits.Header()

hdr2 [‘EXTNAME’] = (’wave’, ’wavelength⊔bins⊔in⊔vacuum’)

hdu0 = fits.PrimaryHDU(flux, header = hdr0)

hdu1 = fit s.ImageHDU(ferr, header=hdr1)

hdu2 = fits.ImageHDU(wave , header=hdr2)

hdu = fits.HDUList ([hdu0 , hdu1, hdu2])

hdu.write to (somepath+ ‘DIBS pec_ISM_spectra.fits’, overwrite =True)

Appendix D Acknowledgements

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

The Gaia mission and data processing have financially been supported by, in alphabetical order by country:

• the Algerian Centre de Recherche en Astronomie, Astrophysique et Géophysique of Bouzareah Observatory;

• the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF) Hertha Firnberg Programme through grants T359, P20046, and P23737;

• the BELgian federal Science Policy Office (BEL-SPO) through various PROgramme de Développement d’Expériences scientifiques (PRODEX) grants of the European Space Agency (ESA), the Research Foundation Flanders (Fonds Wetenschappelijk Onderzoek) through grant VS.091.16N, the Fonds de la Recherche Scientifique (FNRS), and the Research Council of Katholieke Universiteit (KU) Leuven through grant C16/18/005 (Pushing AsteRoseismology to the next level with TESS, GaiA, and the Sloan DIgital Sky SurvEy – PARADISE);

• the Brazil-France exchange programmes Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Coordenação de Aperfeicoamento de Pessoal de Nível Superior (CAPES) - Comité Français d’Evaluation de la Coopération Universitaire et Scientifique avec le Brésil (COFECUB);

• the Chilean Agencia Nacional de Investigación y Desarrollo (ANID) through Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) Regular Project 1210992 (L. Chemin);

• the National Natural Science Foundation of China (NSFC) through grants 11573054, 11703065, and 12173069, the China Scholarship Council through grant 201806040200, and the Natural Science Foundation of Shanghai through grant 21ZR1474100;

• the Tenure Track Pilot Programme of the Croatian Science Foundation and the École Polytechnique Fédérale de Lausanne and the project TTP-2018-07-1171 ‘Mining the Variable Sky’, with the funds of the Croatian-Swiss Research Programme;

• the Czech-Republic Ministry of Education, Youth, and Sports through grant LG 15010 and INTER-EXCELLENCE grant LTAUSA18093, and the Czech Space Office through ESA PECS contract 98058;

• the Danish Ministry of Science;

• the Estonian Ministry of Education and Research through grant IUT40-1;

• the European Commission’s Sixth Framework Programme through the European Leadership in Space Astrometry (ELSA) Marie Curie Research Training Network (MRTN-CT-2006-033481), through Marie Curie project PIOF-GA-2009-255267 (Space AsteroSeismology & RR Lyrae stars, SAS-RRL), and through a Marie Curie Transfer-of-Knowledge (ToK) fellowship (MTKD-CT-2004-014188); the European Commission’s Seventh Framework Programme through grant FP7-606740 (FP7-SPACE-2013-1) for the Gaia European Network for Improved data User Services (GENIUS) and through grant 264895 for the Gaia

Research for European Astronomy Training (GREAT-ITN) network;

• the European Cooperation in Science and Technology (COST) through COST Action CA18104 ‘Revealing the Milky Way with Gaia (MW-Gaia)’;

• the European Research Council (ERC) through grants 320360, 647208, and 834148 and through the European Union’s Horizon 2020 research and innovation and excellent science programmes through Marie Sklodowska-Curie grants 687378 (Small Bodies: Near and Far), 682115 (Using the Magellanic Clouds to Understand the Interaction of Galaxies), 695099 (A sub-percent distance scale from binaries and Cepheids – CepBin), 716155 (Structured ACCREtion Disks – SACCRED), 745617 (Our Galaxy at full HD –Gal-HD), 895174 (The build-up and fate of self-gravitating systems in the Universe), 951549 (Sub-percent calibration of the extragalactic distance scale in the era of big surveys – UniverScale), 101004214 (Innovative Scientific Data Exploration and Exploitation Applications for Space Sciences – EXPLORE), 101004719 (OPTICON-RadioNET Pilot), 101055318 (The 3D motion of the Interstellar Medium with ESO and ESA telescopes – ISM-FLOW), and 101063193 (Evolutionary Mechanisms in the Milky waY; the Gaia Data Release 3 revolution – EMMY);

• the European Science Foundation (ESF), in the framework of the Gaia Research for European Astronomy Training Research Network Programme (GREAT-ESF);

• the European Space Agency (ESA) in the framework of the Gaia project, through the Plan for European Cooperating States (PECS) programme through contracts C98090 and 4000106398/12/NL/KML for Hungary, through contract 4000115263/15/NL/IB for Germany, through PROgramme de Développement d’Expériences scientifiques (PRODEX) grants 4000132054 for Hungary and through contract 4000132226/20/ES/CM;

• the Academy of Finland through grants 299543, 307157, 325805, 328654, 336546, and 345115 and the Magnus Ehrnrooth Foundation;

• the French Centre National d’Études Spatiales (CNES), the Agence Nationale de la Recherche (ANR) through grant ANR-10-IDEX-0001-02 for the ‘Investissements d’avenir’ programme, through grant ANR-15-CE31-0007 for project ‘Modelling the Milky Way in the Gaia era’ (MOD4Gaia), through grant ANR-14-CE33-0014-01 for project ‘The Milky Way disc formation in the Gaia era’ (ARCHEOGAL), through grant ANR-15-CE31-0012-01 for project ‘Unlocking the potential of Cepheids as primary distance calibrators’ (UnlockCepheids), through grant ANR-19-CE31-0017 for project ‘Secular evolution of galaxies’ (SEGAL), and through grant ANR-18-CE31-0006 for project ‘Galactic Dark Matter’ (GaDaMa), the Centre National de la Recherche Scientifique (CNRS) and its SNO Gaia of the Institut des Sciences de l’Univers (INSU), its Programmes Nationaux: Cosmologie et Galaxies (PNCG), Gravitation Références Astronomie Métrologie (PNGRAM), Planétologie (PNP), Physique et Chimie du Milieu Interstellaire (PCMI), and Physique Stellaire (PNPS), supported by INSU along with the Institut National de Physique (INP) and the Institut National de Physique nucléaire et de Physique des Particules (IN2P3), and co-funded by CNES; the ‘Action Fédératrice Gaia’ of the Observatoire de Paris, and the Région de Franche-Comté;

• the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt e.V., DLR) through grants 50QG0501, 50QG0601, 50QG0602, 50QG0701, 50QG0901, 50QG1001, 50QG1101, 50QG1401, 50QG1402, 50QG1403, 50QG1404, 50QG1904, 50QG2101, 50QG2102, and 50QG2202, and the Centre for Information Services and High Performance Computing (ZIH) at the Technische Universität Dresden for generous allocations of computer time;

• the Hungarian Academy of Sciences through the János Bolyai Research Scholarship (G. Marton and Z. Nagy), the Lendület Programme grants LP2014-17 and LP2018-7 and the Hungarian National Research, Development, and Innovation Office (NKFIH) through grant KKP-137523 (‘SeismoLab’);

• the Science Foundation Ireland (SFI) through a Royal Society - SFI University Research Fellowship (M. Fraser);

• the Israel Ministry of Science and Technology through grant 3-18143 and the Israel Science Foundation (ISF) through grant 1404/22;

• the Agenzia Spaziale Italiana (ASI) through contracts I/037/08/0, I/058/10/0, 2014-025-R.0, 2014-025-R.1.2015, and 2018-24-HH.0 and its addendum 2018-24-HH.1-2022 to the Italian Istituto Nazionale di Astrofísica (INAF), contract 2014-049-R.0/1/2, 2022-14-HH.0 to INAF for the Space Science Data Centre (SSDC, formerly known as the ASI Science Data Center, ASDC), contracts I/008/10/0, 2013/030/I.0, 2013-030-I.0.1-2015, and 2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC S.p.A.), INAF, and the Italian Ministry of Education, University, and Research (Ministero dell’Istruzione, dell’Università e della Ricerca) through the Premiale project ‘MIning The Cosmos Big Data and Innovative Italian Technology for Frontier Astrophysics and Cosmology’ (MITiC);

• the Netherlands Organisation for Scientific Research (NWO) through grant NWO-M-614.061.414, through a VICI grant (A. Helmi), and through a Spinoza prize (A. Helmi), and the Netherlands Research School for Astronomy (NOVA);

• the Polish National Science Centre through HARMONIA grant 2018/30/M/ST9/00311 and DAINA grant 2017/27/L/ST9/03221 and the Ministry of Science and Higher Education (MNiSW) through grant DIR/WK/2018/12;

• the Portuguese Fundação para a Ciência e a Tecnologia (FCT) through national funds, grants 2022.06962.PTDC and 2022.03993.PTDC, and work contract DL 57/2016/CP1364/CT0006, grants UIDB/04434/2020 and UIDP/04434/2020 for the Instituto de Astrofísica e Ciências do Espaço (IA), grants UIDB/00408/2020 and UIDP/00408/2020 for LASIGE, and grants UIDB/00099/2020 and UIDP/00099/2020 for the Centro de Astrofísica e Gravitação (CENTRA);

• the Slovenian Research Agency through grant P1-0188;

• the Spanish Ministry of Economy (MINECO/FEDER, UE), the Spanish Ministry of Science and Innovation (MCIN), the Spanish Ministry of Education, Culture, and Sports, and the Spanish Government through grants BES-2016-078499, BES-2017-083126, BES-C-2017-0085, ESP2016-80079-C2-1-R, FPU16/03827, RTI2018-095076-B-C22, PID2021-122842OB-C22, PDC2021-121059-C22, and TIN2015-65316-P (‘Computación de Altas Prestaciones VII’), the Juan de la Cierva Incorporación Programme (FJCI-2015-2671 and IJC2019-04862-I for F. Anders), the Severo Ochoa Centre of Excellence Programme (SEV2015-0493) and MCIN/AEI/10.13039/501100011033/EU FEDER and Next Generation EU/PRTR (PRTR-C17.I1); the European Union through European Regional Development Fund ‘A way of making Europe’ through grants PID2021-122842OB-C21, PID2021-125451NA-I00, CNS2022-13523 and RTI2018-095076-B-C21, the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M, the University of Barcelona’s official doctoral programme for the development of an R+D+i project through an Ajuts de Personal Investigador en Formació (APIF) grant, the <mon>Spanish Virtual Observatory project funded by MCIN/AEI/10.13039/501100011033/through grant PID2020-112949GB-I00; the Centro de Investigación en Tecnologías de la Información y las Comunicaciones (CITIC), funded by the Xunta de Galicia through the collaboration agreement to reinforce CIGUS research centers, research consolidation grant ED431B 2021/36 and scholarships from Xunta de Galicia and the EU - ESF ED481A-2019/155 and ED481A 2021/296; the Red Española de Supercomputación (RES) computer resources at MareNostrum, the Barcelona Supercomputing Centre - Centro Nacional de Supercomputación (BSC-CNS) through activities AECT-2017-2-0002, AECT-2017-3-0006, AECT-2018-1-0017, AECT-2018-2-0013, AECT-2018-3-0011, AECT-2019-1-0010, AECT-2019-2-0014, AECT-2019-3-0003, AECT-2020-1-0004, and DATA-2020-1-0010, the Departament d’Innovació, Universitats i Empresa de la Generalitat de Catalunya through grant 2014-SGR-1051 for project ‘Models de Programació i Entorns d’Execució Parallels’ (MPEXPAR), and Ramon y Cajal Fellowships RYC2018-025968-I, RYC2021-031683-I and RYC2021-033762-I, funded by MICIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR and the European Science Foundation (‘Investing in your future’); the Port d’Informació Científica (PIC), through a collaboration between the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT) and the Institut de Física d’Altes Energies (IFAE), supported by the call for grants for Scientific and Technical Equipment 2021 of the State Program for Knowledge Generation and Scientific and Technological Strengthening of the R+D+i System, financed by MCIN/AEI/ 10.13039/501100011033 and the EU NextGeneration/PRTR (Hadoop Cluster for the comprehensive management of massive scientific data, reference EQC2021-007479-P);

• the Swedish National Space Agency (SNSA/Rymdstyrelsen);

• the Swiss State Secretariat for Education, Research, and Innovation through the Swiss Activités Nationales Complémentaires and the Swiss National Science Foundation through an Eccellenza Professorial Fellowship (award PCEFP2_194638 for R. Anderson);

• the United Kingdom Particle Physics and Astronomy Research Council (PPARC), the United Kingdom Science and Technology Facilities Council (STFC), and the United Kingdom Space Agency (UKSA) through the following grants to the University of Bristol, Brunel University London, the Open University, the University of Cambridge, the University of Edinburgh, the University of Leicester, the Mullard Space Sciences Laboratory of University College London, and the United Kingdom Rutherford Appleton Laboratory (RAL): PP/D006503/1, PP/D006511/1, PP/D006546/1, PP/D006570/1, PP/D006791/1, ST/I000852/1, ST/J005045/1, ST/K00056X/1, ST/K000209/1, ST/K000756/1, ST/K000578/1, ST/L002388/1, ST/L006553/1, ST/L006561/1, ST/N000595/1, ST/N000641/1, ST/N000978/1, ST/N001117/1, ST/S000089/1, ST/S000976/1, ST/S000984/1, ST/S001123/1, ST/S001948/1, ST/S001980/1, ST/S002103/1, ST/V000969/1, ST/W002469/1, ST/W002493/1, ST/W002671/1, ST/W002809/1, EP/V520342/1, ST/X00158X/1, ST/X001601/1, ST/X001636/1, ST/X001687/1, ST/X002667/1, ST/X002683/1 and ST/X002969/1.

The Gaia project and data processing have made use of:

• the Set of Identifications, Measurements, and Bibliography for Astronomical Data (SIMBAD, Wenger et al. 2000), the ‘Aladin sky atlas’ (Bonnarel et al. 2000; Boch & Fernique 2014), and the VizieR catalogue access tool (Ochsenbein et al. 2000), all operated at the Centre de Données astronomiques de Strasbourg (CDS);

• the National Aeronautics and Space Administration (NASA) Astrophysics Data System (ADS);

• the SPace ENVironment Information System (SPENVIS), initiated by the Space Environment and Effects Section (TEC-EES) of ESA and developed by the Belgian Institute for Space Aeronomy (BIRA-IASB) under ESA contract through ESA’s General Support Technologies Programme (GSTP), administered by the BELgian federal Science Policy Office (BELSPO);

• the software products TOPCAT, STIL, and STILTS (Taylor 2005, 2006);

• Matplotlib (Hunter 2007);

• IPython (Pérez & Granger 2007);

• Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018);

• R (R Core Team 2013);

• the HEALPix package (Górski et al. 2005, http://healpix.sourceforge.net/);

• Vaex (Breddels & Veljanoski 2018);

• the HIPPARCOS-2 catalogue (van Leeuwen 2007). The HIPPARCOS and Tycho catalogues were constructed under the responsibility of large scientific teams collaborating with ESA. The Consortia Leaders were Lennart Lindegren (Lund, Sweden: NDAC) and Jean Kovalevsky (Grasse, France: FAST), together responsible for the HIPPARCOS Catalogue; Erik Høg (Copenhagen, Denmark: TDAC) responsible for the Tycho Catalogue; and Catherine Turon (Meudon, France: INCA) responsible for the HIPPARCOS Input Catalogue (HIC);

• the Tycho-2 catalogue (Høg et al. 2000), the construction of which was supported by the Velux Foundation of 1981 and the Danish Space Board;

• The Tycho double star catalogue (TDSC, Fabricius et al. 2002), based on observations made with the ESA HIPPARCOS astrometry satellite, as supported by the Danish Space Board and the United States Naval Observatory through their double-star programme;

• data products from the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006), which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center (IPAC) / California Institute of Technology, funded by the National Aeronautics and Space Administration (NASA) and the National Science Foundation (NSF) of the USA;

• the ninth data release of the AAVSO Photometric All-Sky Survey (APASS, Henden et al. 2016), funded by the Robert Martin Ayers Sciences Fund;

• the first data release of the Pan-STARRS survey (Chambers et al. 2016; Magnier et al. 2020a; Waters et al. 2020; Magnier et al. 2020c,b; Flewelling et al. 2020). The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration (NASA) through grant NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation through grant AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation;

• the second release of the Guide Star Catalogue (GSC2.3, Lasker et al. 2008). The Guide Star Catalogue II is a joint project of the Space Telescope Science Institute (STScI) and the Osservatorio Astrofísico di Torino (OATo). STScI is operated by the Association of Universities for Research in Astronomy (AURA), for the National Aeronautics and Space Administration (NASA) under contract NAS5-26555. OATo is operated by the Italian National Institute for Astrophysics (INAF). Additional support was provided by the European Southern Observatory (ESO), the Space Telescope European Coordinating Facility (STECF), the International GEMINI project, and the European Space Agency (ESA) Astrophysics Division (nowadays SCI-S);

• the eXtended, Large (XL) version of the catalogue of Positions and Proper Motions (PPM-XL, Roeser et al. 2010);

• data products from the Wide-field Infrared Survey Explorer (WISE), which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration (NASA);

• the first data release of the United States Naval Observatory (USNO) Robotic Astrometric Telescope (URAT-1, Zacharias et al. 2015);

• the fourth data release of the United States Naval Observatory (USNO) CCD Astrograph Catalogue (UCAC-4, Zacharias et al. 2013);

• the sixth and final data release of the Radial Velocity Experiment (RAVE DR6, Steinmetz et al. 2020a,b). Funding for RAVE has been provided by the Leibniz Institute for Astrophysics Potsdam (AIP), the Australian Astronomical Observatory, the Australian National University, the Australian Research Council, the French National Research Agency, the German Research Foundation (SPP 1177 and SFB 881), the European Research Council (ERC-StG 240271 Galactica), the Istituto Nazionale di Astrofísica at Padova, the Johns Hopkins University, the National Science Foundation of the USA (AST-0908326), the W.M. Keck foundation, the Macquarie University, the Netherlands Research School for Astronomy, the Natural Sciences and Engineering Research Council of Canada, the Slovenian Research Agency, the Swiss National Science Foundation, the Science & Technology Facilities Council of the UK, Opticon, Strasbourg Observatory, and the Universities of Basel, Groningen, Heidelberg, and Sydney. The RAVE website is at https://www.rave-survey.org/;

• the first data release of the Large sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST DR1, Luo et al. 2015);

• the K2 Ecliptic Plane Input Catalogue (EPIC, Huber et al. 2016);

• the ninth data release of the Sloan Digitial Sky Survey (SDSS DR9, Ahn et al. 2012). Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the United States Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University;

• the thirteenth release of the Sloan Digital Sky Survey (SDSS DR13, Albareti et al. 2017). Funding for SDSS-IV has been provided by the Alfred P. Sloan Foundation, the United States Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is https://www.sdss.org/. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University;

• the second release of the SkyMapper catalogue (SkyMapper DR2, Onken et al. 2019, Digital Object Identifier 10.25914/5ce60d31ce759). The national facility capability for SkyMapper has been funded through grant LE130100104 from the Australian Research Council (ARC) Linkage Infrastructure, Equipment, and Facilities (LIEF) programme, awarded to the University of Sydney, the Australian National University, Swinburne University of Technology, the University of Queensland, the University of Western Australia, the University of Melbourne, Curtin University of Technology, Monash University, and the Australian Astronomical Observatory. SkyMapper is owned and operated by The Australian National University’s Research School of Astronomy and Astrophysics. The survey data were processed and provided by the SkyMapper Team at the Australian National University. The SkyMapper node of the All-Sky Virtual Observatory (ASVO) is hosted at the National Computational Infrastructure (NCI). Development and support the SkyMapper node of the ASVO has been funded in part by Astronomy Australia Limited (AAL) and the Australian Government through the Commonwealth’s Education Investment Fund (EIF) and National Collaborative Research Infrastructure Strategy (NCRIS), particularly the National eResearch Collaboration Tools and Resources (NeCTAR) and the Australian National Data Service Projects (ANDS);

• the Gaia-ESO Public Spectroscopic Survey (GES, Gilmore et al. 2023; Randich et al. 2023). The Gaia-ESO Survey is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 188.B-3002. Public data releases are available through the ESO Science Portal. The project has received funding from the Leverhulme Trust (project RPG-2012-541), the European Research Council (project ERC-2012-AdG 320360-Gaia-ESO-MW), and the Istituto Nazionale di Astrofísica, INAF (2012: CRA 1.05.01.09.16; 2013: CRA 1.05.06.02.07).

The GBOT programme (GBOT) uses observations collected at (i) the European Organisation for Astronomical Research in the Southern Hemisphere (ESO) with the VLT Survey Telescope (VST), under ESO programmes 092.B-0165, 093.B-0236, 094.B-0181, 095.B-0046, 096.B-0162, 097.B-0304, 098.B-0030, 099.B-0034, 0100.B-0131, 0101.B-0156, 0102.B-0174, 0103.B-0165, 0104.B-0081, 0106.20ZA.001 (OmegaCam), 0106.20ZA.002 (FORS2), 0108.21YF; and under INAF programs 110.256C, 112.266Q; and (ii) the Liverpool Telescope, which is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias with financial support from the United Kingdom Science and Technology Facilities Council, and (iii) telescopes of the Las Cumbres Observatory Global Telescope Network.

In addition to the currently active DPAC (and ESA science) authors of the peer-reviewed papers accompanying the data release, there are large numbers of former DPAC members who made significant contributions to the (preparations of the) data processing. In addition to the DPAC consortium, past and present, there are numerous people, mostly in ESA and in industry, who have made or continue to make essential contributions to Gaia, for instance those employed in science and mission operations or in the design, manufacturing, integration, and testing of the spacecraft and its modules, subsystems, and units. Many of those will remain unnamed yet spent countless hours, occasionally during nights, weekends, and public holidays, in cold offices and dark clean rooms. They are acknowledged in the Gaia Documentation

References

All Tables

All Figures

	Fig. 1 Schematic workflow of the DIB-Spec module.
In the text

	Fig. 2 Galactic spatial density distribution of the 6 143 681 RVS spectra used in DIB-Spec. This HEALPix (Górski et al. 2005) map has a level of 7, corresponding to a spatial resolution of 0.46°.
In the text

	Fig. 3 Gaia G–band magnitude distribution of the 6 143 681 RVS objects used in DIB-Spec.
In the text

	Fig. 4 Distribution of A0 from Gaia TGE map (Delchambre et al. 2023) for 160 392 reference stars.
In the text

	Fig. 5 Distribution of Teff, log ɡ, [M/H], the distance of stars, and spectral S/N for both reference stars (in red) and target stars (in blue). These parameters have not been filtered out in the post-processing.
In the text

Fig. 6 Illustration of deriving ISM spectra for individual targets. Left: Schematic view of each main step. Detailed explanations are in Sect. 3.2. Right: Example for a target (Gaia ID: 4169294186397475072) having $T_{\text{eff}}={4445}_{-11}^{+15}$ K, log ɡ = 1.86 ± 0.03 dex, $\left[\text{M/H}\right]={0.12}_{-0.01}^{+0.02}$ dex, and spectral S/N of 110.3. From top to bottom panels: 1) Observed RVS spectrum of this target; 2) First 24 best neighbours for this target; 3) The black line is the target spectrum and the blue line is its stellar template built by averaging 200 best neighbours; 4) Derived ISM spectrum for this target. The positions of the two DIBs are also indicated. The orange shades indicate the regions of two Ca II lines where the weights are set to 0.7 when comparing the target spectrum and reference spectra. The cyan region (860–868 nm) indicates the spectral window used for fitting the two DIBs, with a masked region (green, 866.0–866.8 nm) during the fitting because of the residuals of the Ca II line.

In the text

Fig. 7 Examples of the fits to DIBs λ862.1 and λ864.8 in stacked ISM spectra in five voxels in the same direction, whose HEALPix number (Ipix = 10 450) and GC (ℓc, bc) = (322.43°, −0.44°) are marked in the top panel. The black and red lines are the ISM spectra and fitted DIB profiles, respectively, normalized by the fitted linear continuum. The error bars indicate the flux uncertainties at each pixel. Orange indicates the region that was masked during the fittings. The central heliocentric distance (dc), the number of target spectra (Ntar), mean E(BP − RP), EWs of two DIBs, and the S/N of the stacked ISM spectrum in each voxel are indicated as well.

In the text

	Fig. 8 Distribution of the S/N of the stacked ISM spectra in the DIB depth-width plane (𝒟 vs. σDIB) for DIBs λ862.1 and λ864.8, respectively. The colour represents the mean S/N in each 0.001 × 0.003 nm bin for λ862.1 and 0.001 × 0.01 nm bin for λ864.8.
In the text

	Fig. 9 Flowchart of the criteria to generate quality flags (QFs) for DIBs λ862.1 (upper panel) and λ864.8 (lower panel), respectively. The flag numbers and the corresponding arrow paths for classification are listed in the bottom box in each panel.
In the text

	Fig. 10 Distribution of DIB quality flags (QFs) at each level (0–5, 0 is the best and 5 is the worst) for DIBs λ862.1 (upper panel) and λ864.8 (lower panel). The number of DIB detections and the fraction are indicated at the top of each bar.
In the text

	Fig. 11 Galactic spatial distribution of each level of DIB QF (0–5) for λ862.1 (left panels) and λ864.8 (right panels). The colour represents the counts of DIB detections in each HEALPixel at various selected distances. This HEALPix map has a level of 5, according to a spatial resolution of 1.8°.
In the text

Fig. 12 Distribution of σDIB for DIBs λ862.1 and λ864.8 (top and bottom, respectively). The dashed red and green lines correspond to the “best range” and the “secondary range” of σDIB defined in Sect. 4.1 for the DIB QF, respectively. The distribution of the full DIB-Spec results (235 428 voxels) is shown in the left panels, while the right panels show a quality-controlled sample with S/N > 100 and 𝒟 > 3RC (see Sect. 4.2). The latter criterion implies QF values of 0 or 2. The number of detected DIBs and the percentage after the quality control is indicated as well.

In the text

	Fig. 13 Joint distribution of the FWHM of DIBs λ862.1 and λ864.8, measured in 9778 stacked ISM spectra, with S/N > 100 and 𝒟 > 3RC (see Sect. 4.2), generated by a Gaussian kernel density estimation (KDE). The white star indicates the peak density, and the red line in the central panel indicates the contour of the 2σ level. The orange lines indicate the median FWHM of the two DIBs.
In the text

	Fig. 14 Number density of DIB detections in DIB-Spec as a function of λDIB and σDIB, without any cuts, for λ862.1 (upper panel) and λ864.8 (lower panel). The dashed blue lines indicate the range of permitted λDIB in the QF evaluation (see Sect. 4.1), and the dashed red and green lines correspond to the “best range” and the “secondary range” of σDIB, respectively.
In the text

	Fig. 15 Distributions of the fractional uncertainties of 𝒟, σDIB, and DIB EW as a function of 𝒟 and σDIB for λ862.1 (left panels) and λ864.8 (right panels). Colour represents the mean fractional errors in each 0.001 × 0.003 nm bin for λ862.1 and 0.001 × 0.01 nm bin for λ864.8.
In the text

	Fig. 16 Histogram of the difference between the EW uncertainty estimated in this work ${\sigma }_{{\text{EW}}_1}$ and the one calculated by the formula ${\sigma }_{{\text{EW}}_2}=\sqrt{6\text{\hspace{0.17em}FWHM}\delta \lambda }\times R_C$ (see Sect. 4.5) for λ862.1 (upper panel) and λ864.8 (lower panel), respectively. The mean (Δ) and standard deviation (σ) of the differences are also indicated.
In the text

Fig. 17 Upper panel: comparison between the fitted and integrated EW862.1. The colour represents the number density (estimated by a Gaussian KDE) of the data points in a linear scale. The grey colour bars show the uncertainty of fitted EW862.1. The dashed red line traces the one-to-one correspondence. A zoom-in panel shows the distribution of the EW difference (ΔEW = EWfit − EWint). The mean (Δ) and standard deviation (σ) of the EW difference are indicated. Lower panel: the distribution EW difference as a function of fitted EW862.1.

In the text

Fig. 18 Upper panel: comparison between the EW862.1 measured by DIB-Spec and the mean EW862.1 in each voxel taken from the DIB catalogue in DR3 for 8963 voxels. The colour represents the number density (estimated by a Gaussian KDE) of the data points in a linear scale. The grey colour bars show the EW uncertainty in DIB-Spec and the standard deviation of EW in DR3. The dashed red line traces the one-to-one correspondence. A zoom-in panel shows the distribution of the EW difference (ΔEW = EWFPR − EWDR3). The mean (Δ) and standard deviation (σ) of the EW difference are indicated. Lower panel: the variation of the mean EW difference in each EW bin with a step of 0.01 Å with EW862.1 in FPR. The orange shades the range of 1σ.

In the text

Fig. 19 Distribution of EW862.1 in the Galactic plane (XY), meridian plane (XZ), and rotational plane (YZ) for the FPR results (left panel) and the DR3 results (right panel), respectively, plotted over the Milky Way sketch created by Robert Hurt and Robert Benjamin (Churchwell et al. 2009). Some log-periodic spiral arms described in Reid et al. (2019) are also presented by coloured lines: Scutum-Centaurus Arm in orange; Sagittarius-Carina Arm in purple; Local Arm in black; Perseus Arm in green; Outer Arm in cyan; and the spur between the Local and Sagittarius-Carina arms in blue. The Galactic centre is located at (X, Y, Z) = (8.15, 0, 0).

In the text

	Fig. 20 Galactic distributions of EW862.1 (left panels) and EW864.8 (right panels) from DIB results in FPR in Mollweide projection at HEALPix level 5, at four distances. The distance (dc) and the number of voxels (N) are indicated in each subpanel.
In the text

Fig. 21 Diverse correlations between λ862.1, λ864.8, and mean E(BP − RP) from GSP-Phot (Andrae et al. 2023) in each voxel. The data points are coloured by their number densities estimated by a Gaussian KDE. The red lines are the linear fit to the data points. The dashed green lines indicate the upper limits on the variables for the linear fits. And the fitted slope (α) and intercept (t) are marked in each panel, together with the Pearson correlation coefficient (rp) and the number of DIB detections (N). Dashed black lines are previous results from Z22, and the dashed blue line is from S23.

In the text

Fig. 22 Upper panels: spatial distribution of 12 reliable DIB detections by visual inspection. Their voxels are projected into the Galactic (XY), meridian (XZ), and rotational (YZ) planes. The reds are inside the Local Bubble and the blues are crossed with the surface of the Local Bubble. The grey marks the surface of the Local Bubble determined by Pelgrims et al. (2020). Lower panel: black lines are stacked ISM spectra, and the red/blue lines are the DIB fitting results of the corresponding DIB profiles. The orange indicates the masked spectral region between 866 and 866.8 nm in the fitting. The vertical dashed magenta line indicates the rest-frame wavelength of DIB λ862.1 of 862.323 nm determined in Gaia Collaboration (2023).

In the text

Fig. A.1 Left panel: Distribution of flux residuals (observed – modelled) with spectral wavelength for 159 591 reference spectra. Each row presents one spectrum and the spectra are sorted by their S/N. The dashed blue and red lines indicate the position of S/N equalling 50 and 100, respectively. Some typical stellar lines within the RVS spectral region determined by Contursi et al. (2021) are indicated as well. The lower panel is a zoom-in plot of the upper one to show the distribution of the residuals in the DIB window (860–868 nm). Right panel: Four examples of the reference spectra (black lines with observed flux errors) and their derived stellar templates (red lines). The orange lines are the flux residuals (observed – modelled) with dashed black lines indicating ±5% of the continuum.

In the text

	Fig. A.2 Variation of the mean absolute residual (MAR) between observed and modelled RVS spectra, calculated within the DIB window (861.2−866.0 nm), with the spectral S/N. The dashed green lines indicate S/N = 50 and MAR = 0.02, respectively. The dashed black line is fitted to MAR and S/N on a logarithmic scale.
In the text

Fig. B.1 Corner plot of the DIB fitting in the voxel with Ipix = 10450 and dc = 1.05 kpc (the first panel in Fig. 7 from top to bottom). The histograms and scatter plots show the one- and two-dimensional projections of the posterior distributions of the fitted parameters. The red squares and lines indicate the best-fit estimates for each parameter in the reproduced fitting. And the dashed blue lines mark the fitted parameters in the output table of DIB-Spec.

In the text

	Fig. B.2 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 1.27 kpc (the second panel in Fig. 7 from top to bottom)
In the text

	Fig. B.3 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 2.12 kpc (the third panel in Fig. 7 from top to bottom)
In the text

	Fig. B.4 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 2.40 kpc (the forth panel in Fig. 7 from top to bottom)
In the text

	Fig. B.5 Same as Fig. B.1, but for the voxel with Ipix = 10450 and dc = 3.23 kpc (the fifth panel in Fig. 7 from top to bottom)
In the text