Gaia Focused Product Release: Spatial distribution of two diffuse interstellar bands

Di ff use interstellar bands (DIBs) are absorption features seen in optical and infrared spectra of stars and extragalactic objects that are probably caused by large and complex molecules in the galactic interstellar medium (ISM). Here we investigate the Galactic distribution and properties of two DIBs identified in almost six million stellar spectra collected by the Gaia Radial Velocity Spectrometer. These measurements constitute a part of the Gaia Focused Product Release to be made public between the Gaia DR3 and DR4 data releases. In order to isolate the DIB signal from the stellar features in each individual spectrum, we identified a set of 160 000 spectra at high Galactic latitudes ( | b | ⩾ 65 ◦ ) covering a range of stellar parameters which we consider to be the DIB-free reference sample. Matching each target spectrum to its closest reference spectra in stellar parameter space allowed us to remove the stellar spectrum empirically, without reference to stellar models, leaving a set of six million ISM spectra. Using the star’s parallax and sky coordinates, we then allocated each ISM spectrum to a voxel (VOlume piXEL) on a contiguous three-dimensional grid with an angular size of 1 . 8 ◦ (level 5 HEALPix) and 29 unequally sized distance bins. Identifying the two DIBs at 862.1nm ( λ 862.1) and 864.8nm ( λ 864.8) in the stacked spectra, we modelled their shapes and report the depth, central wavelength, width, and equivalent width (EW) for each, along with confidence bounds on these measurements. We then explored the properties and distributions of these quantities and compared them with similar measurements from other surveys. Our main results are as follows: (1) the strength and spatial distribution of the DIB λ 862.1 are very consistent with what was found in Gaia DR3, but for this work we attained a higher signal-to-noise ratio in the stacked spectra to larger distances, which allowed us to trace DIBs in the outer spiral arm and beyond the Scutum–Centaurus spiral arm; (2) we produced an all-sky map below ± 65 ◦ of Galactic latitude to ∼ 4000pc of both DIB features and their correlations; (3) we detected the signals of DIB λ 862.1 inside the Local Bubble ( ≲ 200pc); and (4) there is a reasonable correlation with the dust reddening found from stellar absorption and EWs of both DIBs with a correlation coe ffi cient of 0.90 for λ 862.1 and 0.


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 A38, page 2 of 33 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 λ 0 = 8623.23± 0.019 Å in vacuum.An average scale height for the carrier of λ862.1 was estimated as 98.60 +11.10  −8.46 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 nm1 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. 2022Zhao et al. , 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.

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. 2018Sartoretti et al. , 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 g, 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 σ V 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.

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.

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 σ V rad ⩽ 5 km s −1 , 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, andS23 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 0 (monochromatic extinction at 541.4 nm Fitzpatrick 1999) of 0.104 mag estimated by the Total Galactic The mean A 0 , 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 0 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 g, 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%).

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 g ± 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 A38, page 4 of 33 Fig. 5. Distribution of T eff , log g, [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.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.0nm).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.

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 HEALPix2 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 −11 K, log g = 1.86 ± 0.03 dex, [M/H] = 0.12 +0.02 −0.01 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.8nm) during the fitting because of the residuals of the Ca II line.
find the best solution for the adaptive distance bins.Instead, S /N ′ = S/N 2 i , 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 √ N 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)): where D and σ DIB are the depth and width of the DIB profile, λ DIB is the measured central wavelength, a 0 and a 1 describe the linear continuum, and λ is the wavelength.Subscripts '862.1' and '864.8' are used to distinguish the profile parameters of the 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 D 862.1 and D 864.8 , 861.2-863.0nm for λ 862.1 , 863.0-866.0nm 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 0 and a 1 .
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. was fixed as 0.12 nm for λ862.1 and 0.4 nm for λ864.8. a 0 and a 1 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 and for λ864.8 as The lower (16%) and upper (84%) confidence levels of EW were estimated by D and σ DIB drawn from the MCMC posterior samplings.Specifically, each {D, σ 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.-Intercept of the linear continuum fitted to the stacked ISM spectrum (a 1 ) 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.D 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 (D 862.1 and λ 862.1 ) are not sensitive to the continuum placement (a 0 and a 1 ), while σ 862.1 presents a weak correlation with a 0 and a 1 .The continuum placement has a heavier effect on the broad and shallow profile of λ864.8.

DIB-Spec outputs
DIB-Specsuccessfully fitted the two DIBs in 235 428 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 Archive3 : (1) 'interstellar_medium_params': the parameter table listing the fitted DIB parameters in each voxel; (2) 'interstellar_medium_spectra': the spectra table containing the  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 convention4 , 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.

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 D and σ DIB .A main part of high S/N (≳200) is found with very small D, corresponding to nearby voxels containing weak or no DIB signals.On the other hand, large D 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 D, 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 ({D, σ DIB }), we generated quality flags (QF) to describe the quality of the fits.This idea comes from Elyajouri et al. (2016) A38, page 9 of 33 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 −1 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 restframe wavelength λ 0 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 −1 because in the literature this DIB was suggested to be located around 8648 Å A38, page 10 of 33 in air.Secondly, based on S23 and Z22, D is required to be smaller than 0.15 for λ862.1 and 0.10 for λ864.8,respectively.D 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.26nm, 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.34nm, 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.2nm, for the best range of σ 864.8 , and 0.6-1.4nm 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 D < R C (case (g), Fig. 9), σ DIB is only compared with the range of 0.1-0.18nm for λ862.1 and 0.6-0.8nm for λ864.8 (case (k), Fig. 9).This choice follows the idea that when D 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 D, since only 4.3% of D 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 {D, σ 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, D 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 D 864.8 is only ∼35% of D 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 D/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 D. 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.

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;Krełowski 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   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 D > 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 FWHM 862.1 = √ 8ln(2) × σ 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 D > 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.

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 {D, σ 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 D − 3σ DIB to D + 3σ DIB A38, page 13 of 33 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 erf(3/ √ 2) × √ 2π × D 862.1 × σ 862.1 ≈ 0.997 EW 862.1 .But it is problematic for the Lorentzian profile of λ864.8 for the integration equals to 2arctan(3) × D 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'.

Effect of stellar residuals
Because most of the target stars processed by DIB-Spec are latetype 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 D < 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.28nm, λ 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.

Uncertainties of DIB parameters
Since D 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 D and σ DIB .Figure 15 shows the distribution of the fractional uncertainties of D, σ DIB , and EW in the D − σ DIB plane for λ862.1 (left panels) and λ864.8 (right panels).The uncertainty of D, σ 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 D or σ DIB .DIB detections with very small D 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 D 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, lowuncertainty 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 D, is much smaller than that of σ DIB , although EW was calculated by D 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 σ EW = √ 6 FWHM δλ × 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.

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.

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 A38, page 14 of 33  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 A38, page 15 of 33 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.

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 EW 862.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.

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     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.

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),  EW 864.8 ,between D 862.1 and D 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 A38, page 18 of 33 Number (a) healpix lc bc dc p08620 p18620 p28620 ew8620 flags8620 Notes. (a) The number of voxels matches those marked in Fig. 22.The rest column names correspond to those in Table 2.
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 D 862.1 and D 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% 5 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 D 862.1 and D 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 sigmaclipping 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 5 The differences between seven studies are listed in Table 3 in Gaia Collaboration (2023) except Wallerstein et al. (2007).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.

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.
A38, page 19 of 33 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.

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 A38, page 20 of 33 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 D ('p08620' and 'p08648' in Table 2) can be found for λ862.1 (166) and λ864.8 (198), although the prior of D is always positive.These cases are all badly fitted due to the low-S/N spectra and/or the weak DIB signals.
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.8nm) 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.

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 235 428 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.

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 (160 392, 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 The flux residuals in the spectra present some patterns, especially in low-S/N spectra, instead of being uniform or noisedominated.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 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.
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.0nm), 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).
A38, page 25 of 33   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.

Fig. 2 .
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. 6 .
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 eff = 4445 +15−11 K, log g = 1.86 ± 0.03 dex, [M/H] = 0.12 +0.02 −0.01 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.8nm) during the fitting because of the residuals of the Ca II line.

Fig. 7 .
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 (I pix = 10 450) and GC (ℓ c , b c ) = (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 (d c ), the number of target spectra (N tar ), mean E(BP − RP), EWs of two DIBs, and the S/N of the stacked ISM spectrum in each voxel are indicated as well.

Fig. 8 .
Fig. 8. Distribution of the S/N of the stacked ISM spectra in the DIB depth-width plane (D 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 .
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 .
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.
(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 A38, page 11 of 33

Fig. 11 .Fig. 12 .
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. 13 .
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 D > 3R C (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.

Fig. 14 .
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.

Fig. 16 .
Fig. 16.Histogram of the difference between the EW uncertainty estimated in this work (σ EW 1 ) and the one calculated by the formula σ EW 2 = √ 6 FWHM δλ × 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 .
Fig. 17.Upper panel: comparison between the fitted and integrated EW 862.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 EW 862.1 .The dashed red line traces the one-to-one correspondence.A zoom-in panel shows the distribution of the EW difference (∆EW = EW fit − EW int ).The mean (∆) and standard deviation (σ) of the EW difference are indicated.Lower panel: the distribution EW difference as a function of fitted EW 862.1 .

Fig. 18 .
Fig. 18.Upper panel: comparison between the EW 862.1 measured by DIB-Spec and the mean EW 862.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 oneto-one correspondence.A zoom-in panel shows the distribution of the EW difference (∆EW = EW FPR − EW DR3 ).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 EW 862.1 in FPR.The orange shades the range of 1σ.
Figure20shows 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.6in S23 for a median EW map at HEALPix level 5), but the distance-sliced maps in Fig.20contain 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

Fig. 19 .
Fig. 19.Distribution of EW 862.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 .
Fig. 20.Galactic distributions of EW 862.1 (left panels) and EW 864.8 (right panels) from DIB results in FPR in Mollweide projection at HEALPix level 5, at four distances.The distance (d c ) and the number of voxels (N) are indicated in each subpanel.
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.

Fig. 21 .
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 (r p ) and the number of DIB detections (N).Dashed black lines are previous results from Z22, and the dashed blue line is from S23.

Fig. 22 .
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).
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.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 10 (MAR) = −0.16− 0.99 × log 10 (S/N) (the dashed black line in Fig. A.

Fig
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.

Fig. A. 2 .
Fig. A.2. Variation of the mean absolute residual (MAR) between observed and modelled RVS spectra, calculated within the DIB window (861.2-866.0nm), 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.

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

import
numpy a s np import p a n d a s a s pd from tqdm import tqdm from a s t r o p y .i o import f i t s d i b s = pd .r e a d _ c s v ( s o m e p a t h + ' p a r a m e t e r _ t a b l e .c s v ' ) t a b e = pd .r e a d _ c s v ( s o m e p a t h + ' s p e c t r a l _ t a b l e .c s v ' ) s p e c _ l c = t a b e .l c .v a l u e s s p e c _ b c = t a b e .bc .v a l u e s s p e c _ d c = t a b e .dc .v a l u e s c a t a _ l c = d i b s .l c .v a l u e s c a t a _ b c = d i b s .bc .v a l u e s c a t a _ d c = d i b s .dc .v a l u e s wave = np .u n i q u e ( t a b e [ ' lambda ' ] .v a l u e s ) f l u x = np .z e r o s ( ( d i b s .s h a p e [ 0 ] , wave .s h a p e [ 0 ] ) , d t y p e = ' f l o a t ' ) f e r r = np .z e r o s _ l i k e ( f l u x ) f o r i i n tqdm ( range ( d i b s .s h a p e [ 0 ] ) ) : t = ( s p e c _ l c == c a t a _ l c [ i ] ) & ( s p e c _ b c == c a t a _ b c [ i ] ) & ( s p e c _ d c == c a t a _ d c [ i ] ) f l u x [ i ] = d a t a [ ' f l u x ' ] [ t ] f e r r [ i ] = d a t a [ ' f l u x _ u n c e r t a i n t y ' ] [ t ] h d r 0 = f i t s .H e a d e r ( ) h d r 0 [ 'COMMENT1' ] = ' s t a c k e d ISM s p e c t r a ' h d r 0 [ 'EXTNAME' ] = ( ' f l u x ' , ' n o r m a l i z e d f l u x ' ) h d r 1 = f i t s .H e a d e r ( ) h d r 1 [ 'EXTNAME' ] = ( ' f e r r ' , ' f l u x u n c e r t a i n t y ' ) h d r 2 = f i t s .H e a d e r ( ) h d r 2 [ 'EXTNAME' ] = ( ' wave ' , ' w a v e l e n g t h b i n s i n vacuum ' ) hdu0 = f i t s .PrimaryHDU ( f l u x , h e a d e r = h d r 0 ) hdu1 = f i t s .ImageHDU ( f e r r , h e a d e r = h d r 1 ) hdu2 = f i t s .ImageHDU ( wave , h e a d e r = h d r 2 ) hdu = f i t s .HDUList ( [ hdu0 , hdu1 , hdu2 ] ) hdu .w r i t e t o ( s o m e p a t h + ' D I B S p e c _ I S M _ s p e c t r a .f i t s ' , o v e r w r i t e = T r u e ) A38, page 29 of 33

Table 1 .
Distance bins defined for the stacking of ISM spectra.
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:

Table 2 .
Column names and their descriptions of the parameter table of DIB-Spec (interstellar_medium_params).

Table 3 .
Column names and their descriptions of the spectra table of DIB-Spec (interstellar_medium_spectra).

Table 4 .
DIB parameters of λ862.1 and filed information of 12 voxels within the Local Bubble (see Sect. 6).