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

M. Schultheis, H. Zhao, T. Zwitter, C. Bailer-Jones, R. Carballo, R. Sordo, R. Drimmel, C. Ordenovic, F. Pailler, M. Fouesneau, et al.

► To cite this version:

M. Schultheis, H. Zhao, T. Zwitter, C. Bailer-Jones, R. Carballo, et al. Gaia Focused Product Release: Spatial distribution of two diffuse interstellar bands. Astronomy and Astrophysics - A&A, 2023, 680, pp.A38. 10.1051/0004-6361/202347103. insu-04353955

HAL Id: insu-04353955

https://insu.hal.science/insu-04353955

Submitted on 20 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gaia Focused Product Release:
## Spatial distribution of two diffuse interstellar bands

Gaia Collaboration: M. Schultheis¹,*, H. Zhao¹,², T. Zwitter³, C. A. L. Bailer-Jones⁴, R. Carballo⁵, R. Sordo⁶, R. Drimmel⁷, C. Ordenovic¹, F. Pailler⁸, M. Fouesneau⁴, O. L. Creevey¹, U. Heiter⁹, A. Recio-Blanco¹, G. Kordopatis¹, P. de Laverny¹, D. J. Marshall¹⁰, T. E. Dharmawardena⁴, A. G. A. Brown¹¹, A. Vallenari⁶, T. Prusti¹², J. H. J. de Bruijne¹², F. Arenou¹³, C. Babusiaux¹⁴, A. Barbier⁸, M. Biermann¹⁵, C. Ducourant¹⁶, D. W. Evans¹⁷, L. Eyer¹⁸, R. Guerra¹⁹, A. Hutton²⁰, C. Jordi²¹,²²,²³, S. A. Klioner²⁴, U. Lammers¹⁹, L. Lindegren²⁵, X. Luri²¹,²²,²³, F. Mignard¹, S. Randich²⁶, P. Sartoretti¹³, R. Smiljanic²⁷, P. Tanga¹, N. A. Walton¹⁷, U. Bastian¹⁵, M. Cropper²⁸, D. Katz¹³, C. Soubiran¹⁶, F. van Leeuwen¹⁷, R. Andrae⁴, M. Audard¹⁸,²⁹, J. Bakker¹⁹, R. Blomme³⁰, J. Castañeda³¹,²¹,²³, F. De Angeli¹⁷, C. Fabricius²³,²¹,²², Y. Frémat³⁰, L. Galluccio¹, A. Guerrier⁸, E. Masana²³,²¹,²², R. Messineo³², C. Nicolas⁸, K. Nienartowicz³³,²⁹, P. Panuzzo¹³, F. Riclet⁸, W. Roux⁸, G. M. Seabroke²⁸, F. Thévenin¹, G. Gracia-Abril³⁴,¹⁵, J. Portell²¹,²²,²³, D. Teyssier³⁵, M. Altmann¹⁵,³⁶, K. Benson²⁸, J. Berthier³⁷, P. W. Burgess¹⁷, D. Busonero⁷, G. Busso¹⁷, H. Cánovas³⁵, B. Carry¹, N. Cheek³⁸, G. Clementini³⁹, Y. Damerdji⁴⁰,⁴¹, M. Davidson⁴², P. de Teodoro¹⁹, L. Delchambre⁴⁰, A. Dell'Oro²⁶, E. Fraile Garcia⁴³, D. Garabato⁴⁴, P. García-Lario¹⁹, N. Garralda Torres⁴⁵, P. Gavras⁴³, R. Haigron¹³, N. C. Hambly⁴², D. L. Harrison¹⁷,⁴⁶, D. Hatzidimitriou⁴⁷, J. Hernández¹⁹, S. T. Hodgkin¹⁷, B. Holl¹⁸,²⁹, S. Jamal⁴, S. Jordan¹⁵, A. Krone-Martins⁴⁸,⁴⁹, A. C. Lanzafame⁵⁰,⁵¹, W. Löffler¹⁵, A. Lorca²⁰, O. Marchal⁵², P. M. Marrese⁵³,⁵⁴, A. Moitinho⁴⁹, K. Muinonen⁵⁵,⁵⁶, M. Nuñez Campos²⁰, I. Oreshina-Slezak¹, P. Osborne¹⁷, E. Pancino²⁶,⁵⁴, T. Pauwels³⁰, M. Riello¹⁷, L. Rimoldini²⁹, A. C. Robin⁵⁷, T. Roegiers⁵⁸, L. M. Sarro⁵⁹, C. Siopis⁶⁰, M. Smith²⁸, A. Sozzetti⁷, E. Utrilla²⁰, M. van Leeuwen¹⁷, K. Weingrill⁶¹, U. Abbas⁷, P. Ábrahám⁶²,⁶³, A. Abreu Aramburu⁴⁵, C. Aerts⁶⁴,⁶⁵,⁴, G. Altavilla⁵³,⁵⁴, M. A. Álvarez⁴⁴, J. Alves⁶⁶, F. Anders²¹,²²,²³, R. I. Anderson⁶⁷, T. Antoja²¹,²²,²³, D. Baines⁶⁸, S. G. Baker²⁸, Z. Balog¹⁵,⁴, C. Barache³⁶, D. Barbato¹⁸,⁷, M. Barros⁶⁹, M. A. Barstow⁷⁰, S. Bartolomé²³,²¹,²², D. Bashi⁷¹,⁷², N. Bauchet¹³, N. Baureau⁷³, U. Becciani⁵⁰, L. R. Bedin⁶, I. Bellas-Velidis⁷⁴, M. Bellazzini³⁹, W. Beordo⁷,⁷⁵, A. Berihuete⁷⁶, M. Bernet²¹,²²,²³, C. Bertolotto³², S. Bertone⁷, L. Bianchi⁷⁷, A. Binnenfeld⁷⁸, A. Blazere⁷⁹, T. Boch⁵², A. Bombrun⁸⁰, S. Bouquillon³⁶,⁸¹, A. Bragaglia³⁹, J. Braine¹⁶, L. Bramante³², E. Breedt¹⁷, A. Bressan⁸², N. Brouillet¹⁶, E. Brugaletta⁵⁰, B. Bucciarelli⁷,⁷⁵, A. G. Butkevich⁷, R. Buzzi⁷, E. Caffau¹³, R. Cancelliere⁸³, S. Cannizzo⁸⁴, T. Carlucci³⁶, M. I. Carnerero⁷, J. M. Carrasco²³,²¹,²², J. Carretero⁸⁵,⁸⁶, S. Carton⁸⁴, L. Casamiquela¹⁶,¹³, M. Castellani⁵³, A. Castro-Ginard¹¹, V. Cesare⁵⁰, P. Charlot¹⁶, L. Chemin⁷⁸, V. Chiaramida³², A. Chiavassa¹, N. Chornay¹⁷,²⁹, R. Collins⁴², G. Contursi¹, W. J. Cooper⁸⁸,⁷, T. Cornez⁸⁴, M. Crosta⁷,⁸⁹, C. Crowley⁸⁰, C. Dafonte⁴⁴, F. De Luise⁹⁰, R. De March³², R. de Souza⁹¹, A. de Torres⁸⁰, E. F. del Peloso¹⁵, M. Delbo¹, A. Delgado⁴³, S. Diakite⁹², C. Diener¹⁷, E. Distefano⁵⁰, C. Dolding²⁸, K. Dsilva⁶⁰, J. Durán⁴³, H. Enke⁶¹, P. Esquej⁴³, C. Fabre⁷⁹, M. Fabrizio⁵³,⁵⁴, S. Faigler⁷¹, M. Fatović⁹³, G. Fedorets⁵⁵,⁹⁴, J. Fernández-Hernández⁴³, P. Fernique⁵², F. Figueras²¹,²²,²³, Y. Fournier⁶¹, C. Fouron⁷³, M. Gai⁷, M. Galinier¹, A. Garcia-Gutierrez²³,²¹,²², M. García-Torres⁹⁵, A. Garofalo³⁹, E. Gerlach²⁴, R. Geyer²⁴, P. Giacobbe⁷, G. Gilmore¹⁷,⁹⁶, S. Girona⁹⁷, G. Giuffrida⁵³, R. Gomel⁷¹, A. Gomez⁴⁴, J. González-Núñez⁹⁸, I. González-Santamaría⁴⁴, E. Gosset⁴⁰,⁹⁹, M. Granvik⁵⁵,¹⁰⁰, V. Gregori Barrera²³,²¹,²², R. Gutiérrez-Sánchez³⁵, M. Haywood¹³, A. Helmer⁸⁴, A. Helmi¹⁰¹, K. Henares⁶⁸, S. L. Hidalgo¹⁰²,¹⁰³, T. Hilger²⁴, D. Hobbs²⁵, C. Hottier¹³, H. E. Huckle²⁸, M. Jabłońska¹⁰⁴,¹⁰⁵, F. Jansen¹⁰⁶, Ó. Jiménez-Arranz²¹,²²,²³, J. Juaristi Campillo¹⁵, S. Khanna⁷,¹⁰¹, A. J. Korn⁹, Á. Kóspál⁶²,⁴,⁶³, Z. Kostrzewa-Rutkowska¹¹, M. Kun⁶², S. Lambert³⁶, A. F. Lanza⁵⁰, J.-F. Le Campion¹⁶, Y. Lebreton¹⁰⁷,¹⁰⁸, T. Lebzelter⁶⁶, S. Leccia¹⁰⁹, I. Lecoeur-Taibi²⁹, G. Lecoutre⁵⁷, S. Liao¹¹⁰,⁷,¹¹¹, L. Liberato¹,¹¹², E. Licata⁷, H. E. P. Lindstrøm⁷,¹¹³,¹¹⁴, T. A. Lister¹¹⁵, E. Livanou⁴⁷, A. Lobel³⁰, C. Loup⁵², L. Mahy³⁰, R. G. Mann⁴², M. Manteiga¹¹⁶, J. M. Marchant¹¹⁷, M. Marconi¹⁰⁹, D. Marín Pina²¹,²²,²³, S. Marinoni⁵³,⁵⁴, J. Martín Lozano³⁸, J. M. Martín-Fleitas²⁰, G. Marton⁶², N. Mary⁸⁴, A. Masip²³,²¹,²², D. Massari³⁹, A. Mastrobuono-Battisti¹³, T. Mazeh⁷¹, P. J. McMillan²⁵, J. Meichsner²⁴, S. Messina⁵⁰, D. Michalik¹²,

<footnote>* Corresponding author: M. Schultheis, e-mail: mathias.schultheis@oca.eu</footnote>
N. R. Millar¹⁷, A. Mints⁶¹, D. Molina²²,²¹,²³, R. Molinaro¹⁰⁹, L. Molnár⁶²,¹¹⁸,⁶³, G. Monari⁵²,
M. Monguiò²¹,²²,²³, P. Montegriffo³⁹, A. Montero³⁸, R. Mor¹¹⁹,²²,²³, A. Mora²⁰, R. Morbidelli⁷, T. Morel⁴⁰,
D. Morris⁴², N. Mowlavi¹⁸, D. Munoz⁸⁴, T. Muraveva³⁹, C. P. Murphy¹⁹, I. Musella¹⁰⁹, Z. Nagy⁶², S. Nieto⁴³,
L. Noval⁸⁴, A. Ogden¹⁷, C. Pagani¹²⁰, I. Pagano⁵⁰, L. Palaversa⁹³, P. A. Palicio¹, L. Pallas-Quintela⁴⁴,
A. Panahi⁷¹, C. Panem⁸, S. Payne-Wardenaar¹⁵, L. Pegoraro⁸, A. Penttilä⁵⁵, P. Pesciullesi⁴³, A. M. Piersimoni⁹⁰,
M. Pinamonti⁷, F.-X. Pineau⁵², E. Plachy⁶²,¹¹⁸,⁶³, G. Plum¹³, E. Poggio¹,⁷, D. Pourbaix†,⁶⁰,⁹⁹, A. Prša¹²¹,
L. Pulone⁵³, E. Racero³⁸,¹²², M. Rainer²⁶,¹²³, C. M. Raiteri⁷, P. Ramos¹²⁴,²¹,²³, M. Ramos-Lerate³⁵,
M. Ratajczak¹⁰⁴, P. Re Fiorentin⁷, S. Regibo⁶⁴, C. Reylé⁵⁷, V. Ripepi¹⁰⁹, A. Riva⁷, H.-W. Rix⁴,
G. Rixon¹⁷, N. Robichon¹³, C. Robin⁸⁴, M. Romero-Gómez²¹,²²,²³, N. Rowell⁴², F. Royer¹³,
D. Ruz Mieres¹⁷, K. A. Rybicki¹²⁵, G. Sadowski⁶⁰, A. Sáez Núñez²³,²¹,²², A. Sagristà Sellés¹⁵,
J. Sahlmann⁴³, V. Sanchez Gimenez²³,²¹,²², N. Sanna²⁶, R. Santoveña⁴⁴, M. Sarasso⁷, C. Sarrate Riera³¹,²¹,²³,
E. Sciacca⁵⁰, J. C. Segovia³⁸, D. Ségransan¹⁸, S. Shahaf¹²⁵, A. Siebert⁵²,¹²⁶, L. Siltala⁵⁵, E. Slezak¹,
R. L. Smart⁷,⁸⁸, O. N. Snaith¹³,¹²⁷, E. Solano¹²⁸, F. Solitro³², D. Souami¹⁰⁷,¹²⁹, J. Souchay³⁶, L. Spina⁶,
E. Spitoni¹,¹³⁰, F. Spoto¹³¹, L. A. Squillante³², I. A. Steele¹¹⁷, H. Steidelmüller²⁴, J. Surdej⁴⁰, L. Szabados⁶²,
F. Taris³⁶, M. B. Taylor¹³², R. Teixeira⁹¹, K. Tisanić⁹³, L. Tolomei³², F. Torra³¹,²¹,²³,
G. Torralba Elipe⁴⁴,¹³³,¹³⁴, M. Trabucchi¹³⁵,¹⁸, M. Tsantaki²⁶, A. Ulla¹³⁶,¹³⁷, N. Unger¹⁸, O. Vanel¹³,
A. Vecchiato⁷, D. Vicente⁹⁷, S. Voutsinas⁴², M. Weiler²³,²¹,²², L. Wyrzykowski¹⁰⁴, J. Zorec¹³⁸,
L. Balaguer-Núñez²³,²¹,²², N. Leclerc¹³, S. Morgenthaler¹³⁹, G. Robert⁸⁴, and S. Zucker⁷⁸

(Affiliations can be found after the references)

Received 6 June 2023 / Accepted 5 September 2023

ABSTRACT

Diffuse 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.1 nm (λ862.1) and 864.8 nm (λ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 ~4000 pc of both DIB features and their correlations; (3) we detected the signals of DIB λ862.1 inside the Local Bubble (≤200 pc); and (4) there is a reasonable correlation with the dust reddening found from stellar absorption and EWs of both DIBs with a correlation coefficient of 0.90 for λ862.1 and 0.77 for λ864.8.

**Key words.** ISM: lines and bands – dust, extinction

# 1. Introduction

Diffuse interstellar bands (DIBs) are a set of ubiquitous inter-
stellar absorption features that primarily exist in the optical and
near-infrared wavelength range (about 0.4–2.4 µm) of the spec-
tra 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 presum-
ably 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 high-
resolution 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, here-
after 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 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<fnref n="1" /> 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.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic workflow of the DIB-Spec module.</figcaption></figure>

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.

## 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 = \lambda/\Delta\lambda \sim 11500$ (Cropper et al. 2018). The RVS spectra were then processed by GSP-Spec to estimate their stellar atmospheric parameters (effective temperature $T_{\text{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 \ge 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 $\omega$, stellar radial velocity $V_{\text{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.

<footnote n="1">¹ The accurate rest-frame wavelength of the DIB λ864.8 has not been determined, and we therefore name it λ864.8 following previous suggestions.</footnote>
<figure><img src="image_3.png" /><figcaption>**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°.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Gaia G-band magnitude distribution of the 6 143 681 RVS objects used in DIB-Spec.</figcaption></figure>

## 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_{rad}} \le 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^\circ$ and, 160 392 reference stars with $|b| \ge 65^\circ$. 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 $\lambda862.1$ in 247 sightlines with $|b| \ge 65^\circ$ 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

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Distribution of $A_0$ from Gaia TGE map (Delchambre et al. 2023) for 160 392 reference stars.</figcaption></figure>

extinction (TGE) map of Gaia DR3 (Delchambre et al. 2023). 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%).

## 3.2. Derivation of ISM spectra for individual targets with the BNM

To detect and measure the two DIBs, $\lambda862.1$ and $\lambda864.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} \pm 20\%$, $\log g \pm 0.6$ dex, and [M/H] $\pm 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
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Distribution of $T_{\text{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.</figcaption></figure>

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 4595–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_{\text{eff}}$) or rare reference stars in its vicinity will not get a generated stellar template.

### 3.3. Stacking and fitting in each voxel

The ISM spectra of individual targets are stacked in each voxel according to their equatorial coordinates $(\alpha, \delta)$ and the distances $(1/\sigma)$ 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 $(\alpha, \delta)$ is done at level 5 in the HEALPix<fnref n="2" /> scheme (Górski et al. 2005), corresponding to 12 288 pixels with an equivalent spatial resolution of $1.8^\circ$. 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 $(\alpha, \delta)$ is fixed as $1.8^\circ$, 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

<footnote n="2">² https://healpix.sourceforge.io</footnote>
<figure><img src="image_7.png" /><figcaption>**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^{+15}_{-11}$ K, $\log g = 1.86 \pm 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.8 nm) during the fitting because of the residuals of the Ca II line.</figcaption></figure>

find the best solution for the adaptive distance bins. Instead, $S/N' = \sqrt{\sum 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' \ge 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' \ge 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' \ge 200$ in each distance bin. The fraction could exceed 85% when we consider $S/N' \ge 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, $\lambda862.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 $\lambda862.1$ profile. Thus, we still assume a Gaussian profile for $\lambda862.1$ in this work and treat the possible departures as a source of uncertainty. Rare studies of $\lambda864.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 $\lambda442.8$ (Snow 2002) while $\lambda864.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 $\lambda862.1$, a Lorentzian function (Eq. (2)) for $\lambda864.8$, and a linear function for continuum (Eq. (3)):

$$G(\lambda; \mathcal{D}, \lambda_{\text{DIB}}, \sigma_{\text{DIB}}) = -\mathcal{D} \times \exp\left(-\frac{(\lambda - \lambda_{\text{DIB}})^2}{2\sigma_{\text{DIB}}^2}\right), \quad (1)$$

$$L(\lambda; \mathcal{D}, \lambda_{\text{DIB}}, \sigma_{\text{DIB}}) = \frac{-(\mathcal{D}\sigma_{\text{DIB}}^2)}{(\lambda - \lambda_{\text{DIB}})^2 + \sigma_{\text{DIB}}^2}, \quad (2)$$

$$C(\lambda; a_0, a_1) = a_0 \times \lambda + a_1, \quad (3)$$

where $\mathcal{D}$ and $\sigma_{\text{DIB}}$ are the depth and width of the DIB profile, $\lambda_{\text{DIB}}$ is the measured central wavelength, $a_0$ and $a_1$ describe the linear continuum, and $\lambda$ is the wavelength. Subscripts ‘862.1’ and ‘864.8’ are used to distinguish the profile parameters of the
<table><caption>**Table 1.** Distance bins defined for the stacking of ISM spectra.</caption><thead><tr><th>Distance bin (pc)</th><th>Pixel number<sup>(a)</sup></th><th>S/N &ge; 200 (%)</th></tr></thead><tbody><tr><td>0–130</td><td>11 179</td><td>48.7%</td></tr><tr><td>130–200</td><td>11 202</td><td>67.3%</td></tr><tr><td>200–270</td><td>11 204</td><td>68.0%</td></tr><tr><td>270–350</td><td>11 207</td><td>66.6%</td></tr><tr><td>350–430</td><td>11 203</td><td>52.9%</td></tr><tr><td>430–530</td><td>11 204</td><td>50.2%</td></tr><tr><td>530–650</td><td>11 200</td><td>48.4%</td></tr><tr><td>650–780</td><td>11 189</td><td>51.3%</td></tr><tr><td>780–950</td><td>11 167</td><td>64.8%</td></tr><tr><td>950–1150</td><td>11 146</td><td>66.1%</td></tr><tr><td>1150–1400</td><td>11 135</td><td>63.4%</td></tr><tr><td>1400–1670</td><td>11 093</td><td>53.0%</td></tr><tr><td>1670–2000</td><td>11 008</td><td>46.3%</td></tr><tr><td>2000–2250</td><td>10 656</td><td>30.9%</td></tr><tr><td>2250–2550</td><td>10 416</td><td>27.3%</td></tr><tr><td>2550–2950</td><td>10 050</td><td>27.3%</td></tr><tr><td>2950–3500</td><td>9464</td><td>26.1%</td></tr><tr><td>3500–4500</td><td>9060</td><td>25.0%</td></tr><tr><td>4500–5000</td><td>6659</td><td>7.9%</td></tr><tr><td>5000–5500</td><td>5943</td><td>5.2%</td></tr><tr><td>5500–6000</td><td>5164</td><td>3.0%</td></tr><tr><td>6000–6500</td><td>4528</td><td>1.2%</td></tr><tr><td>6500–7000</td><td>3876</td><td>0.5%</td></tr><tr><td>7000–7500</td><td>3392</td><td>0.3%</td></tr><tr><td>7500–8000</td><td>2959</td><td>0.1%</td></tr><tr><td>8000–8500</td><td>2541</td><td>0.0%</td></tr><tr><td>8500–9000</td><td>2147</td><td>0.1%</td></tr><tr><td>9000–9500</td><td>1838</td><td>0.0%</td></tr><tr><td>9500–10000</td><td>1598</td><td>0.1%</td></tr></tbody></table>

**Notes.** (a) Number of HEALPixels with $N_{\text{tar}} \ge 1$ in each distance bin.

two DIBs. The full DIB model, $M_Θ$, is the sum of Eqs. (1)–(3), where $Θ = \{D_{862.1}, λ_{862.1}, σ_{862.1}, D_{864.8}, λ_{864.8}, σ_{864.8}, a_0, a_1\}$ 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) \propto P(y|Θ)P(Θ)$. $P(y|Θ)$ is the likelihood:

$$P(y|Θ) = \frac{1}{\sqrt{2\pi}\sigma_y} \exp \left[ -\frac{1}{2\sigma_y^2} (y - M_Θ)^2 \right], \quad (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 $D_{862.1}$ and $D_{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_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. 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 $D$ and $λ_{\text{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 $σ_{\text{DIB}}$

<figure><img src="image_8.png" /></figure>

**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_{\text{pix}} = 10450$) and GC ($l_c, b_c = (322.43^\circ, -0.44^\circ)$) 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_{\text{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.

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

$$\text{EW}_{862.1} = \sqrt{2\pi} \times D_{862.1} \times \sigma_{862.1}, \quad (5)$$

and for $λ_{864.8}$ as

$$\text{EW}_{864.8} = \pi \times D_{864.8} \times \sigma_{864.8}. \quad (6)$$

The lower (16%) and upper (84%) confidence levels of EW were estimated by $D$ and $σ_{\text{DIB}}$ drawn from the MCMC posterior samplings. Specifically, each $\{D, σ_{\text{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.
<table><caption>**Table 2.** Column names and their descriptions of the parameter table of DIB-Spec (interstellar_medium_params).</caption><thead><tr><th>Column</th><th>Name</th><th>Unit</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>solution_id</td><td>–</td><td>Solution identifier</td></tr><tr><td>2</td><td>healpix</td><td>–</td><td>HEALPix identification</td></tr><tr><td>3</td><td>lc</td><td>deg</td><td>Central Galactic longitude of each voxel</td></tr><tr><td>4</td><td>bc</td><td>deg</td><td>Central Galactic latitude of each voxel</td></tr><tr><td>5</td><td>dc</td><td>kpc</td><td>Central heliocentric distance of each voxel</td></tr><tr><td>6</td><td>n_targets</td><td>–</td><td>Number of target stars in each voxel</td></tr><tr><td>7</td><td>snr</td><td>–</td><td>S/N of the stacked ISM spectrum</td></tr><tr><td>8</td><td>ew8620</td><td>nm</td><td>Equivalent width of DIB at 862.1 nm (EW<sub>862.1</sub>)</td></tr><tr><td>9</td><td>ew8620_lower</td><td>nm</td><td>Lower confidence level (16%) of equivalent width of DIB λ862.1</td></tr><tr><td>10</td><td>ew8620_lower</td><td>nm</td><td>Upper confidence level (84%) of equivalent width of DIB λ862.1</td></tr><tr><td>11</td><td>flags8620</td><td>–</td><td>Quality flag of the parameters of DIB at 862.1 nm</td></tr><tr><td>12</td><td>p08620</td><td>–</td><td>Depth of DIB at 862.1 nm (D<sub>862.1</sub>)</td></tr><tr><td>13</td><td>p08620_lower</td><td>–</td><td>Lower confidence level (16%) of depth of DIB λ862.1</td></tr><tr><td>14</td><td>p08620_lower</td><td>–</td><td>Upper confidence level (84%) of depth of DIB λ862.1</td></tr><tr><td>15</td><td>p18620</td><td>nm</td><td>Central wavelength of DIB at 862.1 nm (λ<sub>862.1</sub>)</td></tr><tr><td>16</td><td>p18620_lower</td><td>nm</td><td>Lower confidence level (16%) of central wavelength of DIB λ862.1</td></tr><tr><td>17</td><td>p18620_lower</td><td>nm</td><td>Upper confidence level (84%) of central wavelength of DIB λ862.1</td></tr><tr><td>18</td><td>p28620</td><td>nm</td><td>Gaussian width of DIB at 862.1 nm (σ<sub>862.1</sub>)</td></tr><tr><td>19</td><td>p28620_lower</td><td>nm</td><td>Lower confidence level (16%) of Gaussian width of DIB λ862.1</td></tr><tr><td>20</td><td>p28620_lower</td><td>nm</td><td>Upper confidence level (84%) of Gaussian width of DIB λ862.1</td></tr><tr><td>21</td><td>ew8648</td><td>nm</td><td>Equivalent width of DIB at 864.8 nm (EW<sub>864.8</sub>)</td></tr><tr><td>22</td><td>ew8648_lower</td><td>nm</td><td>Lower confidence level (16%) of equivalent width of DIB λ864.8</td></tr><tr><td>23</td><td>ew8648_lower</td><td>nm</td><td>Upper confidence level (84%) of equivalent width of DIB λ864.8</td></tr><tr><td>24</td><td>flags8648</td><td>–</td><td>Quality flag of the parameters of DIB at 864.8 nm</td></tr><tr><td>25</td><td>p08648</td><td>–</td><td>Depth of DIB at 864.8 nm (D<sub>864.8</sub>)</td></tr><tr><td>26</td><td>p08648_lower</td><td>–</td><td>Lower confidence level (16%) of depth of DIB λ864.8</td></tr><tr><td>27</td><td>p08648_lower</td><td>–</td><td>Upper confidence level (84%) of depth of DIB λ864.8</td></tr><tr><td>28</td><td>p18648</td><td>nm</td><td>Central wavelength of DIB at 864.8 nm (λ<sub>864.8</sub>)</td></tr><tr><td>29</td><td>p18648_lower</td><td>nm</td><td>Lower confidence level (16%) of central wavelength of DIB λ864.8</td></tr><tr><td>30</td><td>p18648_lower</td><td>nm</td><td>Upper confidence level (84%) of central wavelength of DIB λ864.8</td></tr><tr><td>31</td><td>p28648</td><td>nm</td><td>Lorentzian width of DIB at 864.8 nm (σ<sub>864.8</sub>)</td></tr><tr><td>32</td><td>p28648_lower</td><td>nm</td><td>Lower confidence level (16%) of Lorentzian width of DIB λ864.8</td></tr><tr><td>33</td><td>p28648_lower</td><td>nm</td><td>Upper confidence level (84%) of Lorentzian width of DIB λ864.8</td></tr><tr><td>34</td><td>dibcont_a0</td><td>–</td><td>Slope of the linear continuum fitted to the stacked ISM spectrum (a<sub>0</sub>)</td></tr><tr><td>35</td><td>dibcont_a0_lower</td><td>–</td><td>Lower confidence level (16%) of the slope of continuum</td></tr><tr><td>36</td><td>dibcont_a0_lower</td><td>–</td><td>Upper confidence level (84%) of the slope of continuum</td></tr><tr><td>37</td><td>dibcont_a1</td><td>–</td><td>Intercept of the linear continuum fitted to the stacked ISM spectrum (a<sub>1</sub>)</td></tr></tbody></table>

Figure 7 shows examples of the fits for five stacked ISM spectra in voxels towards the same direction, that is HEALPix number of 10 450, corresponding to Galactic coordinates of the voxel of $(l_c, b_c) = (322.43^\circ, -0.44^\circ)$. With increasing voxel central heliocentric distance ($d_c$) from top to bottom, EW<sub>862.1</sub> 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 $\sigma_{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 $\lambda_{862.1}$) are not sensitive to the continuum placement ($a_0$ and $a_1$), while $\sigma_{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.

## 4. DIB-Spec outputs

DIB-Spec successfully 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 Archive<fnref n="3" />: (1) 'interstellar_medium_params': the parameter table listing the fitted DIB parameters in each voxel; (2) 'interstellar_mediumspectra': the spectra table containing the

<footnote n="3">³ https://gea.esac.esa.int/archive/</footnote>
<table><caption>**Table 3.** Column names and their descriptions of the spectra table of DIB-Spec (interstellar_medium_spectra).</caption><thead><tr><th>Column</th><th>Name</th><th>Unit</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>solution_id</td><td>–</td><td>Solution identifier</td></tr><tr><td>2</td><td>healpix</td><td>–</td><td>HEALPix identification</td></tr><tr><td>3</td><td>lc</td><td>deg</td><td>Central Galactic longitude of each voxel</td></tr><tr><td>4</td><td>bc</td><td>deg</td><td>Central Galactic latitude of each voxel</td></tr><tr><td>5</td><td>dc</td><td>kpc</td><td>Central heliocentric distance of each voxel</td></tr><tr><td>6</td><td>lambda</td><td>Å</td><td>Wavelength in vacuum</td></tr><tr><td>7</td><td>flux</td><td>–</td><td>Normalized flux</td></tr><tr><td>8</td><td>flux_uncertainty</td><td>–</td><td>Uncertainty of the normalized flux</td></tr></tbody></table>

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_al’) 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<fnref n="4" />, 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% (12692) 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 62859276 rows which equal 235428 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

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Distribution of the S/N of the stacked ISM spectra in the DIB depth–width plane ($D$ vs. $\sigma_{DIB}$) for DIBs $\lambda862.1$ and $\lambda864.8$, respectively. The colour represents the mean S/N in each $0.001 \times 0.003$ nm bin for $\lambda862.1$ and $0.001 \times 0.01$ nm bin for $\lambda864.8$.</figcaption></figure>

determined simply by S/N. On one hand, the distribution of S/N shown in Fig. 8 presents a dependence on both $D$ and $\sigma_{DIB}$. A main part of high S/N ($\ge 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 $\sigma_{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 $\sigma_{DIB} \sim 0.22$ nm for DIB $\lambda862.1$, and between 0.8 and 1.0 nm for DIB $\lambda864.8$, but not particularly obvious due to the fact that $\lambda864.8$ is much broader and shallower than $\lambda862.1$.

Considering both S/N and the shape of the fitted DIB profile ($\{D, \sigma_{DIB}\}$), we generated quality flags (QF) to describe the quality of the fits. This idea comes from Elyajouri et al. (2016)

<footnote n="4">⁴ See Bastian & Portell (2020): Source Identifiers – Assignment and Usage throughout DPAC (GAIA-C3-TN-ARI-BAS-020), which can be accessed through https://www.cosmos.esa.int/web/gaia/public-dpac-documents</footnote>
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

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

$\lambda_{\text{DIB}}$ to be consistent with a Doppler wavelength shift within about $\pm 200 \text{ 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 $\lambda_{\text{DIB}}$ and the rest-frame wavelength $\lambda_0$ of the two DIBs reported in S23 and Z22. We note that the upper limit of $\lambda_{\text{DIB}}$ for λ864.8 is a bit larger than a wavelength corresponding to $200 \text{ km s}^{-1}$ because in the literature this DIB was suggested to be located around 8648 Å
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

in air. Secondly, based on S23 and Z22, $\mathcal{D}$ is required to be smaller than 0.15 for λ862.1 and 0.10 for λ864.8, respectively. $\mathcal{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.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 $\mathcal{D} < 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 $\mathcal{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 $\mathcal{D}$, since only 4.3% of $\mathcal{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 $\{\mathcal{D}, \sigma_{\text{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, $\mathcal{D}$ is also expected to decrease with $|b|$. Thus, we can find numerous detections with QF = 1 and 3 at $|b| \gtrsim 20^\circ$ 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| \sim 30^\circ$) towards the inner disc ($|\ell| < 30^\circ$) and the Galactic anti-centre ($\ell \sim 180^\circ$) 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 $\mathcal{D}_{864.8}$ is only ~35% of $\mathcal{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 $\mathcal{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 $\mathcal{D}$. QF = 5 is distributed in almost the full sky with $|b| < 65^\circ$ 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.

## 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; Krawł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 (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
<figure><img src="image_12.png" /><figcaption>**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 HEALPix at various selected distances. This HEALPix map has a level of 5, according to a spatial resolution of 1.8°.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

**Fig. 12.** Distribution of $\sigma_{\text{DIB}}$ for DIBs $\lambda862.1$ and $\lambda864.8$ (top and bottom, respectively). The dashed red and green lines correspond to the “best range” and the “secondary range” of $\sigma_{\text{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 $D > 3R_c$ (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.

S/N or extremely weak DIB signals. The second peak at $\sigma_{862.1} = 0.12$ nm, the initial guess of $\sigma_{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 $\sigma_{862.1} = 0.22$ nm, as we have seen in Fig. 8, marks the best fits of DIB $\lambda862.1$. For $\lambda864.8$ (lower left panel), the first ($\sigma_{864.8} = 0.1$ nm) and second ($\sigma_{864.8} = 0.4$ nm) peaks can also be found. But the width distribution is then very flat for larger $\sigma_{864.8}$, and the peak of $\sigma_{864.8}$ that indicates the best $\lambda864.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 $\sigma_{\text{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 $\lambda862.1$ and 9927 for $\lambda864.8$, whose width distributions are shown in the right panels in Fig. 12. The peaks at the initial guess of $\sigma_{\text{DIB}}$ disappear, and a quasi-Gaussian distribution of $\sigma_{\text{DIB}}$ can now be found for both DIBs with a lower cut that excludes very narrow spurious features. The selected ranges of $\sigma_{\text{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 $\sigma_{864.8}$ in QF evaluation is enlarged because of the long tail of the distribution of $\sigma_{864.8}$ showing more detections with $\sigma_{864.8} \gtrsim 1.1$ nm than with $\sigma_{864.8} \lesssim 0.8$ nm. The cause is unknown and DIB $\lambda864.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(2)} \times \sigma_{862.1}$ for the Gaussian profile (Eq. (1)) of $\lambda862.1$ and $\text{FWHM}_{864.8} = 2 \times \sigma_{864.8}$ for the Lorentzian profile (Eq. (2)) of $\lambda864.8$. The joint distribution of $\text{FWHM}_{862.1}$ and $\text{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 $\text{FWHM}_{864.8}$. The median FWHM of DIB $\lambda862.1$ is $0.52 \pm 0.05$ nm, which is consistent with Zhao et al. (2022, $0.55 \pm 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 $\text{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

<figure><img src="image_14.png" /></figure>

**Fig. 13.** Joint distribution of the FWHM of DIBs $\lambda862.1$ and $\lambda864.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\sigma$ level. The orange lines indicate the median FWHM of the two DIBs.

RVS spectra, Saydjari et al. (2023) got a $\sigma_{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 $\lambda862.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 $\lambda862.1$, respectively, which are smaller than mentioned results above. Puspitarini et al. (2015) mentioned that many Gaia-ESO spectra of the $\lambda862.1$ region are contaminated by sky emission residuals which fall within the red wing of $\lambda862.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 $\sim 4.7$ Å ($\sigma_{862.1}$ around $2$ Å) could be proper for DIB $\lambda862.1$. But we should keep in mind that the DIB width could vary under different ISM environments. We note that DIB $\lambda862.1$ could contain multiple components (Jenniskens & Desert 1994) and present an asymmetric profile (e.g. Krewlowski et al. 2019). But considering the resolution of RVS spectra, our study cannot reveal an accurate shape of $\lambda862.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 $\lambda864.8$ is $1.91 \pm 0.44$ nm, larger than $1.62 \pm 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 $\lambda864.8$ due to its very large breadth. Z22 also measured the median $\sigma_{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 $\sigma_{864.8}$. The measurements in early studies could be questionable due to the broad span of $\lambda864.8$ and its superposition with some blended stellar lines. The median $\sigma_{\text{DIB}}$ reported in this work is from a quasi-Gaussian distribution.

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

The EW of the DIBs $\lambda862.1$ and $\lambda864.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, \sigma_{\text{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\sigma_{\text{DIB}}$ to $D + 3\sigma_{\text{DIB}}$.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Number density of DIB detections in DIB-Spec as a function of $\lambda_{\text{DIB}}$ and $\sigma_{\text{DIB}}$, without any cuts, for $\lambda_{862.1}$ (upper panel) and $\lambda_{864.8}$ (lower panel). The dashed blue lines indicate the range of permitted $\lambda_{\text{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 $\sigma_{\text{DIB}}$, respectively.</figcaption></figure>

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}(3/\sqrt{2}) \times \sqrt{2\pi} \times D_{862.1} \times \sigma_{862.1} \approx 0.997 \text{EW}_{862.1}$. But it is problematic for the Lorentzian profile of $\lambda_{864.8}$ for the integration equals to $2\arctan(3) \times D_{864.8} \times \sigma_{864.8} \approx 0.795 \text{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’.

### 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 $\lambda_{\text{DIB}}$ or the broadening of $\sigma_{\text{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 $\sigma_{\text{DIB}}$ is very small.

Figure 14 shows the distribution of DIB detections as a function of $\lambda_{\text{DIB}}$ (in stellar frame) and $\sigma_{\text{DIB}}$, which can straightforwardly identify the effect of stellar residuals. A vertical stripe at the initial guess of $\sigma_{\text{DIB}}$ can be seen for both the two DIBs that corresponds to the noisy cases, mostly with $D<R_C$. $\lambda_{862.1}$ is uniform and widespread when $\sigma_{862.1} < 0.04$ nm, where DIB detection is noise dominated. For larger $\sigma_{862.1}$, especially 0.16–0.28 nm, $\lambda_{862.1}$ is concentrated between 862.2 and 862.5 nm, representing the reliable measurements. Our detection of $\lambda_{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 $\sigma_{862.1} - \lambda_{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 $\lambda_{864.8}$ makes it difficult for accurate measurements of $\lambda_{864.8}$ and $\sigma_{864.8}$, resulting in a very scattered $\lambda_{864.8}$ versus $\sigma_{864.8}$ distribution. However, the large $\sigma_{864.8}$ can hardly be affected by the stellar residuals.

### 4.5. Uncertainties of DIB parameters

Since $D$ and $\sigma_{\text{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 $\sigma_{\text{DIB}}$. Figure 15 shows the distribution of the fractional uncertainties of $D$, $\sigma_{\text{DIB}}$, and EW in the $D - \sigma_{\text{DIB}}$ plane for $\lambda_{862.1}$ (left panels) and $\lambda_{864.8}$ (right panels). The uncertainty of $D$, $\sigma_{\text{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 $\sigma_{\text{DIB}}$. DIB detections with very small $D$ will get large QF (low reliability) by comparing with $R_C$, and those with small $\sigma_{\text{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 $\sigma_{\text{DIB}}$ in priors, 0.5 nm for $\lambda_{862.1}$ and 1.5 nm for $\lambda_{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 $\sigma_{\text{DIB}}$ will significantly decrease when applying strict constraints. For $\lambda_{862.1}$, detections with low uncertainty (<10%) are located in similar fields for $D_{862.1}$, $\sigma_{862.1}$, and EW$_{862.1}$ and correspond to the “best range” of $\sigma_{862.1}$ defined in Sect. 4.1. On the other hand, low-uncertainty regions for $\lambda_{864.8}$ do not present a well-constrained range for $\sigma_{864.8}$.

The uncertainty of EW in some regions, mainly the bottom in each panel with small $\sigma_{\text{DIB}}$ and large $D$, is much smaller than that of $\sigma_{\text{DIB}}$, although EW was calculated by $D$ and $\sigma_{\text{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 ($\delta\lambda = 0.03$ nm/pixel), and the noise level of the line centre ($R_C = \text{std}(\text{data} - \text{model}))$, as $\sigma_{\text{EW}} = \sqrt{6 \text{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 $\lambda_{862.1}$ and within 0.07 Å for $\lambda_{864.8}$), but our method systematically gives out smaller $\sigma_{\text{EW}}$ for $\lambda_{862.1}$ and larger $\sigma_{\text{EW}}$ for $\lambda_{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
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Distributions of the fractional uncertainties of $\mathcal{D}$, $\sigma_{\text{DIB}}$, and DIB EW as a function of $\mathcal{D}$ and $\sigma_{\text{DIB}}$ for $\lambda_{862.1}$ (left panels) and $\lambda_{864.8}$ (right panels). Colour represents the mean fractional errors in each $0.001 \times 0.003$ nm bin for $\lambda_{862.1}$ and $0.001 \times 0.01$ nm bin for $\lambda_{864.8}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**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{FWHM}} \delta\lambda \times R_C$ (see Sect. 4.5) for $\lambda_{862.1}$ (upper panel) and $\lambda_{864.8}$ (lower panel), respectively. The mean ($\Delta$) and standard deviation ($\sigma$) of the differences are also indicated.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**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 ($\Delta$EW = EW$_{\text{fit}}$ - EW$_{\text{int}}$). The mean ($\Delta$) and standard deviation ($\sigma$) of the EW difference are indicated. Lower panel: the distribution EW difference as a function of fitted EW$_{862.1}$.</figcaption></figure>

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 $\lambda_{862.1}$ cannot be directed integrated because the profiles of $\lambda_{862.1}$ and $\lambda_{864.8}$ could be overlapped. Therefore, we first subtracted the fitted profile of $\lambda_{864.8}$ from the ISM spectrum and then integrated the rest part of the ISM

spectrum within $\lambda_{862.1} \pm 3\sigma_{862.1}$ to calculate the integrated EW of $\lambda_{862.1}$. Figure 17 shows the comparison between fitted and integrated EW$_{862.1}$ for 6240 cases with $S/N \ge 100$, $\text{QF}_{862.1} = 0$, and $\text{QF}_{864.8} = 0$. We did not do such a comparison for $\lambda_{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
<figure><img src="image_19.png" /></figure>

**Fig. 18.** *Upper panel:* comparison between the EW<sub>862.1</sub> measured by DIB-Spec and the mean EW<sub>862.1</sub> 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 = EW<sub>FPR</sub> − EW<sub>DR3</sub>). 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<sub>862.1</sub> in FPR. The orange shades the range of 1σ.

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<sub>862.1</sub>. After the check of the ISM spectra, we find that the outliers with much larger fitted EW<sub>862.1</sub> than integrated EW<sub>862.1</sub> 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.

## 5.2. Comparison with DIB measurements in Gaia DR3

To make a direct comparison of the measurements of EW<sub>862.1</sub> between FPR and DR3, the mean EW<sub>862.1</sub> of the full DIB catalogue in DR3 is calculated in each voxel defined in DIB-Spec. Figure 18 presents the comparison between EW<sub>862.1</sub> from FPR and mean EW<sub>862.1</sub> 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<sub>862.1</sub> (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<sub>862.1</sub> 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<sub>862.1</sub> than DR3 when EW<sub>862.1</sub> ≥ 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<sub>862.1</sub> would be underestimated. This effect will be investigated in detail in follow-up work.

To compare the spatial distribution of EW<sub>862.1</sub> 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 ≤ σ<sub>DIB</sub> ≤ 0.26 nm (the best range of σ<sub>862.1</sub> determined in DIB-Spec), resulting in 52 180 DIB detections. Then the median EW<sub>862.1</sub> 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<sub>862.1</sub> 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<sub>862.1</sub> in that voxel. To have a clean map with reliable detections, we require S/N ≥ 100 and QF<sub>862.1</sub> ≤ 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, $\ell \sim -70^\circ$ and $\ell \sim -117^\circ$, have significant low EW<sub>862.1</sub> 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 $\sim 80^\circ - \sim 90^\circ$ and $\sim -100^\circ - \sim -90^\circ$ 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<sub>862.1</sub> (left panels) and EW<sub>864.8</sub> (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 $(\alpha, \delta, 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 ($\ell \sim 30^\circ-40^\circ$) and Cygnus complex ($\ell \sim 80^\circ$) in the middle plane and Ophiuchus ($\ell \sim 0^\circ$ and $b \sim 15^\circ$) at high latitude, although the resolution of the present DIB EW
<figure><img src="image_20.png" /><figcaption>**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).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**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.</figcaption></figure>
map is pretty low (1.8°). Nevertheless, some clouds seem to dis-
appear 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 longi-
tudinal 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 sam-
ple 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. Fur-
thermore, 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 con-
tain reliable DIB detections. This observational bias could also
happen in nearby regions behind dense clouds.

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 $\lambda862.1$, $\lambda864.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 $\lambda862.1$ and $\lambda864.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 $\lambda862.1$ and $\lambda864.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 ($\alpha$), intercept ($t$),

<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Diverse correlations between $\lambda862.1$, $\lambda864.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 ($\alpha$) 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.</figcaption></figure>

and Pearson coefficient ($r_p$) are marked in each panel in Fig. 21
as well.

Tight linearity can be found between EW<sub>862.1</sub> and EW<sub>864.8</sub>, between D<sub>862.1</sub> and D<sub>864.8</sub>, and between EW<sub>862.1</sub> and E(BP-RP), and a weaker one for EW<sub>864.8</sub> and E(BP-RP). In the correlation of EW for the two DIBs, some cases present greater EW<sub>864.8</sub> than expected when EW<sub>862.1</sub> ≥ 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
<table><caption>**Table 4.** DIB parameters of λ862.1 and filed information of 12 voxels within the Local Bubble (see Sect. 6).</caption><thead><tr><th>Number<sup>(a)</sup></th><th>healpix</th><th>lc<br>(deg)</th><th>bc<br>(deg)</th><th>dc<br>(kpc)</th><th>p08620</th><th>p18620<br>(nm)</th><th>p28620<br>(nm)</th><th>ew8620<br>(nm)</th><th>flags8620</th></tr></thead><tbody><tr><td>1</td><td>6212</td><td>304.28</td><td>36.91</td><td>0.165</td><td>0.0082</td><td>862.20</td><td>0.19</td><td>0.0038</td><td>1</td></tr><tr><td>2</td><td>6348</td><td>284.14</td><td>47.74</td><td>0.165</td><td>0.0075</td><td>862.28</td><td>0.25</td><td>0.0047</td><td>1</td></tr><tr><td>3</td><td>7685</td><td>4.39</td><td>19.37</td><td>0.065</td><td>0.0087</td><td>862.33</td><td>0.22</td><td>0.0049</td><td>1</td></tr><tr><td>4</td><td>4793</td><td>59.97</td><td>-33.73</td><td>0.165</td><td>0.0081</td><td>862.25</td><td>0.34</td><td>0.0068</td><td>2</td></tr><tr><td>5</td><td>2312</td><td>40.31</td><td>43.28</td><td>0.165</td><td>0.0055</td><td>862.24</td><td>0.24</td><td>0.0033</td><td>1</td></tr><tr><td>6</td><td>8047</td><td>51.13</td><td>16.46</td><td>0.165</td><td>0.0112</td><td>862.35</td><td>0.13</td><td>0.0036</td><td>2</td></tr><tr><td>7</td><td>9695</td><td>272.04</td><td>29.70</td><td>0.235</td><td>0.0086</td><td>862.28</td><td>0.18</td><td>0.0039</td><td>1</td></tr><tr><td>8</td><td>10 008</td><td>257.39</td><td>11.41</td><td>0.235</td><td>0.0069</td><td>862.27</td><td>0.25</td><td>0.0043</td><td>1</td></tr><tr><td>9</td><td>5853</td><td>188.17</td><td>-18.54</td><td>0.165</td><td>0.0087</td><td>862.35</td><td>0.24</td><td>0.0052</td><td>1</td></tr><tr><td>10</td><td>5235</td><td>226.30</td><td>-8.49</td><td>0.235</td><td>0.0099</td><td>862.37</td><td>0.24</td><td>0.0059</td><td>1</td></tr><tr><td>11</td><td>10 041</td><td>249.64</td><td>12.91</td><td>0.235</td><td>0.0115</td><td>862.32</td><td>0.23</td><td>0.0067</td><td>0</td></tr><tr><td>12</td><td>7347</td><td>7.37</td><td>13.09</td><td>0.165</td><td>0.0151</td><td>862.25</td><td>0.22</td><td>0.0082</td><td>0</td></tr></tbody></table>

**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\sigma$-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 $\lambda864.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\sigma$-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.

## 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 $\lambda578.0$ and $\lambda579.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 $\lambda862.1$ that are very close to the Sun and generally suggested that $\lambda862.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 $\lambda862.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 $\lambda862.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 $\ge$ 300, $QF_{862.1} \le$ 2, $QF_{864.8} \le$ 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 $\lambda862.1$ and field information of these 12 voxels are listed in Table 4.

<footnote>$^5$ The differences between seven studies are listed in Table 3 in Gaia Collaboration (2023) except Wallerstein et al. (2007).</footnote>
<figure><img src="image_23.png" /><figcaption>**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).</figcaption></figure>

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} \le 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.

# 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/\sigma$) of the background stars was directly used in DIB-Spec for stacking spectra in different voxels. But $1/\sigma$ 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/\sigma$, 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 *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.

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<sub>862.1</sub> < EW<sub>862.1,lower</sub>, 351 cases with EW<sub>862.1</sub> < EW<sub>862.1,upper</sub>, 2791 cases with EW<sub>864.8</sub> < EW<sub>864.8,lower</sub>, and 440 cases with EW<sub>864.8</sub> < EW<sub>864.8,upper</sub>. 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 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<sub>862.1</sub> 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.

**References**

Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21

Albareti, F. D., Allende Prieto, C., Almeida, A., et al. 2017, ApJS, 233, 25

Andrae, R., Fouesneau, M., Sordo, R., et al. 2023, A&A, 674, A27

Astropy Collaboration (Price-Whelan, A., et al.) 2018, AJ, 156, 123

Bailer-Jones, C. A. L. 2015, PASP, 127, 994

Bailer-Jones, C. A. L., Andrae, R., Arcay, B., et al. 2013, A&A, 559, A74

Baron, D., Poznanski, D., Watson, D., et al. 2015a, MNRAS, 451, 332

Baron, D., Poznanski, D., Watson, D., Yao, Y., & Prochaska, J. X. 2015b, MNRAS, 447, 545

Boch, T., & Fernique, P. 2014, in ASP Conf. Ser., 485, Astronomical Data Analysis Software and Systems XXIII, eds. N. Manset, & P. Forshay, 277

Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, A&AS, 143, 33

Breddels, M. A., & Veljanoski, J. 2018, A&A, 618, A13

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, ArXiv e-prints [arXiv:1612.05560]

Churchwell, E., Babler, B. L., Meade, M. R., et al. 2009, PASP, 121, 213

Contursi, G., de Laverny, P., Recio-Blanco, A., & Palicio, P. A. 2021, A&A, 654, A130

Cordiner, M. A., Fossey, S. J., Smith, A. M., & Sarre, P. J. 2013, ApJ, 764, L10

Creevey, O. L., Sordo, R., Pailler, F., et al. 2023, A&A, 674, A26

Cropper, M., Katz, D., Sartoretti, P., et al. 2018, A&A, 616, A5

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792

Delchambre, L., Bailer-Jones, C. A. L., Bellas-Velidis, I., et al. 2023, A&A, 674, A31

Ebenbichler, A., Postel, A., Przybilla, N., et al. 2022, A&A, 662, A81

Elyajouri, M., & Lallement, R. 2019, A&A, 628, A67

Elyajouri, M., Monreal-Ibero, A., Remy, Q., & Lallement, R. 2016, ApJS, 225, 19

Fabricius, C., Høg, E., Makarov, V. V., et al. 2002, A&A, 384, 180

Fan, H., Hobbs, L. M., Dahlstrom, J. A., et al. 2019, ApJ, 878, 151

Farhang, A., van Loon, J. T., Khosroshahi, H. G., Javadi, A., & Bailey, M. 2019, Nat. Astron., 3, 922

Fitzpatrick, E. L. 1999, PASP, 111, 63
Flewelling, H. A., Magnier, E. A., Chambers, K. C., et al. 2020, ApJS, 251, 7

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Friedman, S. D., York, D. G., McCall, B. J., et al. 2011, ApJ, 727, 33

Gaia Collaboration (Schultheis, M., et al.) 2023, A&A, 674, A40

Gilmore, G., Randich, S., Worley, C. C., et al. 2023, A&A, 666, A120

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hamano, S., Kobayashi, N., Kawakita, H., et al. 2022, ApJS, 262, 2

Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalogue: II/336

Herbig, G. H., & Leka, K. D. 1991, ApJ, 382, 193

Hobbs, L. M., York, D. G., Snow, T. P., et al. 2008, ApJ, 680, 1256

Hobbs, L. M., York, D. G., Thorburn, J. A., et al. 2009, ApJ, 705, 32

Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

Huber, D., Bryson, S. T., Haas, M. R., et al. 2016, ApJS, 224, 2

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jenniskens, P., & Desert, F. X. 1994, A&AS, 106, 39

Kos, J., Zwitter, T., Grebel, E. K., et al. 2013, ApJ, 778, 86

Kos, J., Zwitter, T., Wyse, R., et al. 2014, Science, 345, 791

Krełowski, J., Galazutdinov, G. A., Godunova, V., & Bondar, A. 2019, Acta Astron., 69, 159

Krełowski, J., Galazutdinov, G. A., Gnaciński, P., et al. 2021, MNRAS, 508, 4241

Lai, T. S. Y., Witt, A. N., Alvarez, C., & Cami, J. 2020, MNRAS, 492, 5853

Lallement, R., Vergely, J. L., Valette, B., et al. 2014, A&A, 561, A91

Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135

Lan, T.-W., Ménard, B., & Zhu, G. 2015, MNRAS, 452, 3629

Lasker, B. M., Lattanzi, M. G., McLean, B. J., et al. 2008, AJ, 136, 735

Luo, A. L., Zhao, Y.-H., Zhao, G., et al. 2015, Res. Astron. Astrophys., 15, 1095

Magnier, E. A., Chambers, K. C., Flewelling, H. A., et al. 2020a, ApJS, 251, 3

Magnier, E. A., Schlafly, E. F., Finkbeiner, D. P., et al. 2020b, ApJS, 251, 6

Magnier, E. A., Sweeney, W. E., Chambers, K. C., et al. 2020c, ApJS, 251, 5

Maíz Apellániz, J. 2015, Mem. Soc. Astron. Italiana, 86, 553

Monreal-Ibero, A., Weilbacher, P. M., Wendt, M., et al. 2015, A&A, 576, A3

Monreal-Ibero, A., Weilbacher, P. M., & Wendt, M. 2018, A&A, 615, A33

Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23

Onken, C. A., Wolf, C., Bessell, M. S., et al. 2019, PASA, 36, e033

Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., & Montier, L. 2020, A&A, 636, A17

Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21

Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109

Puspitarini, L., Lallement, R., Babusiaux, C., et al. 2015, A&A, 573, A35

Randich, S., Gilmore, G., Magrini, L., et al. 2023, A&A, 666, A121

R Core Team 2013, R: A Language and Environment for Statistical Computing (Vienna, Austria: R Foundation for Statistical Computing)

Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2023, A&A, 674, A29

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131

Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440

Sartoretti, P., Katz, D., Cropper, M., et al. 2018, A&A, 616, A6

Sartoretti, P., Marchal, O., Babusiaux, C., et al. 2023, A&A, 674, A6

Saydjari, A. K., M. Uzsoy, A. S., Zucker, C., Peek, J. E. G., & Finkbeiner, D. P. 2023, ApJ, 954, 141

Seabroke, G. M., Fabricius, C., Teyssier, D., et al. 2021, A&A, 653, A160

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Snow, T. P. 2002, ApJ, 567, 407

Steinmetz, M., Guiglion, G., McMillan, P. J., et al. 2020a, AJ, 160, 83

Steinmetz, M., Matijevič, G., Enke, H., et al. 2020b, AJ, 160, 82

Tielens, A. G. G. M. 2014, in The Diffuse Interstellar Bands, 297, eds. J. Cami, & N. L. J. Cox, 399

Taylor, M. B. 2005, in ASP Conf. Ser., 347, Astronomical Data Analysis Software and Systems XIV, eds. P. Shopbell, M. Britton, & R. Ebert, 29

Taylor, M. B. 2006, in ASP Conf. Ser., 351, Astronomical Data Analysis Software and Systems XV, eds. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 666 van Leeuwen, F. 2007, A&A, 474, 653

Vogrinčič, R., Kos, J., Zwitter, T., et al. 2023, MNRAS, 521, 3727

Vollmann, K., & Eversberg, T. 2006, Astron. Nachr., 327, 862

Vos, D. A. I., Cox, N. L. J., Kaper, L., Spaans, M., & Ehrenfreund, P. 2011, A&A, 533, A129

Wallerstein, G., Sandstrom, K., & Gredel, R. 2007, PASP, 119, 1268

Waters, C. Z., Magnier, E. A., Price, P. A., et al. 2020, ApJS, 251, 4

Welsh, B. Y., Lallement, R., Vergely, J. L., & Raimond, S. 2010, A&A, 510, A54

Welty, D. E. 2014, in The Diffuse Interstellar Bands, 207, eds. J. Cami, & N. L. J. Cox, 153

Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9

Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44

Zacharias, N., Finch, C., Subasavage, J., et al. 2015, AJ, 150, 101

Zasowski, G., Ménard, B., Bizyaev, D., et al. 2015, ApJ, 798, 35

Zhao, H., Schultheis, M., Recio-Blanco, A., et al. 2021a, A&A, 645, A14

Zhao, H., Schultheis, M., Rojas-Arriagada, A., et al. 2021b, A&A, 654, A116

Zhao, H., Schultheis, M., Zwitter, T., et al. 2022, A&A, 666, A12

Zhao, H., Schultheis, M., Arentsen, A., et al. 2023, MNRAS, 519, 754

1 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France

2 Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023, PR China

3 Faculty of Mathematics and Physics, University of Ljubljana, Jadranska ulica 19, 1000 Ljubljana, Slovenia

4 Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

5 Dpto. de Matemática Aplicada y Ciencias de la Computación, Univ. de Cantabria, ETS Ingenieros de Caminos, Canales y Puertos, Avda. de los Castros s/n, 39005 Santander, Spain

6 INAF – Osservatorio astronomico di Padova, Vicolo Osservatorio 5, 35122 Padova, Italy

7 INAF – Osservatorio Astrofísico di Torino, via Osservatorio 20, 10025 Pino Torinese (TO), Italy

8 CNES Centre Spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

9 Observational Astrophysics, Division of Astronomy and Space Physics, Department of Physics and Astronomy, Uppsala University, Box 516, 751 20 Uppsala, Sweden

10 IRAP, Université de Toulouse, CNRS, UPS, CNES, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France

11 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

12 European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201AZ, Noordwijk, The Netherlands

13 GEPI, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France

14 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

15 Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr. 12-14, 69120 Heidelberg, Germany

16 Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France

17 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

18 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland

19 European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

20 Aurora Technology for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

21 Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), Martí i Franquès 1, 08028 Barcelona, Spain

22 Departament de Física Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), c. Martí i Franquès 1, 08028 Barcelona, Spain

23 Institut d'Estudis Espacials de Catalunya (IEEC), c. Gran Capità 2-4, 08034 Barcelona, Spain

24 Lohrmann Observatory, Technische Universität Dresden, Mommsenstraße 13, 01062 Dresden, Germany

25 Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, 22100 Lund, Sweden

26 INAF – Osservatorio Astrofísico di Arcetri, Largo Enrico Fermi 5, 50125 Firenze, Italy

27 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

28 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK
29 Department of Astronomy, University of Geneva, Chemin d'Ecogia 16, 1290 Versoix, Switzerland

30 Royal Observatory of Belgium, Ringlaan 3, 1180 Brussels, Belgium

31 DAPCOM Data Services, c. dels Vilabella 5-7, 80500 Vic, Barcelona, Spain

32 ALTEC S.p.a., Corso Marche 79, 10146 Torino, Italy

33 Sednai Sàrl, Geneva, Switzerland

34 Gaia DPAC Project Office, ESAC, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

35 Telespazio UK S.L. for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

36 SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, LNE, 61 avenue de l'Observatoire, 75014 Paris, France

37 IMCCE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Lille, 77 av. Denfert-Rochereau, 75014 Paris, France

38 Serco Gestión de Negocios for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

39 INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Piero Gobetti 93/3, 40129 Bologna, Italy

40 Institut d'Astrophysique et de Géophysique, Université de Liège, 19c Allée du 6 Août, 4000 Liège, Belgium

41 CRAAG – Centre de Recherche en Astronomie, Astrophysique et Géophysique, Route de l'Observatoire, Bp 63, Bouzareah 16340 Algiers, Algeria

42 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

43 RHEA for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

44 CIGUS CITIC – Department of Computer Science and Information Technologies, University of A Coruña, Campus de Elviña s/n, A Coruña, 15071, Spain

45 ATG Europe for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

46 Kavli Institute for Cosmology Cambridge, Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK

47 Department of Astrophysics, Astronomy and Mechanics, National and Kapodistrian University of Athens, Panepistimiopolis, Zografos, 15783 Athens, Greece

48 Donald Bren School of Information and Computer Sciences, University of California, Irvine, CA 92697, USA

49 CENTRA, Faculdade de Ciências, Universidade de Lisboa, Edif. C8, Campo Grande, 1749-016 Lisboa, Portugal

50 INAF – Osservatorio Astrofisico di Catania, via S. Sofia 78, 95123 Catania, Italy

51 Dipartimento di Fisica e Astronomia “Ettore Majorana”, Università di Catania, Via S. Sofia 64, 95123 Catania, Italy

52 Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 11 rue de l'Université, 67000 Strasbourg, France

53 INAF – Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monte Porzio Catone (Roma), Italy

54 Space Science Data Center – ASI, Via del Politecnico SNC, 00133 Roma, Italy

55 Department of Physics, University of Helsinki, PO Box 64, 00014 Helsinki, Finland

56 Finnish Geospatial Research Institute FGI, Vuorimiehentie 5, 02150 Espoo, Finland

57 Institut UTINAM CNRS UMR6213, Université de Franche-Comté, OSU THETA Franche-Comté Bourgogne, Observatoire de Besançon, BP1615, 25010 Besançon Cedex, France

58 HE Space Operations BV for European Space Agency (ESA), Keplerlaan 1, 2201AZ, Noordwijk, The Netherlands

59 Dpto. de Inteligencia Artificial, UNED, c/ Juan del Rosal 16, 28040 Madrid, Spain

60 Institut d'Astronomie et d'Astrophysique, Université Libre de Bruxelles CP 226, Boulevard du Triomphe, 1050 Brussels, Belgium

61 Leibniz Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany

62 Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), MTA Centre of Excellence, Konkoly Thege Miklós út 15–17, 1121 Budapest, Hungary

63 ELTE Eötvös Loránd University, Institute of Physics, 1117, Pázmány Péter sétány 1A, Budapest, Hungary

64 Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium

65 Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands

66 University of Vienna, Department of Astrophysics, Türkenschanzstraße 17, A1180 Vienna, Austria

67 Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

68 Quasar Science Resources for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

69 LASIGE, Faculdade de Ciências, Universidade de Lisboa, Edif. C6, Campo Grande, 1749-016 Lisboa, Portugal

70 School of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK

71 School of Physics and Astronomy, Tel Aviv University, Tel Aviv 6997801, Israel

72 Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK

73 Telespazio for CNES Centre Spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

74 National Observatory of Athens, I. Metaxa and Vas. Pavlou, Palaia Penteli, 15236 Athens, Greece

75 University of Turin, Department of Physics, Via Pietro Giuria 1, 10125 Torino, Italy

76 Depto. Estadística e Investigación Operativa. Universidad de Cádiz, Avda. República Saharaui s/n, 11510 Puerto Real, Cádiz, Spain

77 EURIX S.r.l., Corso Vittorio Emanuele II 61, 10128, Torino, Italy

78 Porter School of the Environment and Earth Sciences, Tel Aviv University, Tel Aviv 6997801, Israel

79 ATOS for CNES Centre Spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

80 HE Space Operations BV for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain

81 LFCA/DAS, Universidad de Chile, CNRS, Casilla 36-D, Santiago, Chile

82 SISSA – Scuola Internazionale Superiore di Studi Avanzati, via Bonomea 265, 34136 Trieste, Italy

83 University of Turin, Department of Computer Sciences, Corso Svizzera 185, 10149 Torino, Italy

84 Thales Services for CNES Centre Spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

85 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain

86 Port d'Informació Científica (PIC), Campus UAB, C. Albareda s/n, 08193 Bellaterra (Barcelona), Spain

87 Instituto de Astrofísica, Universidad Andres Bello, Fernandez Concha 700, Las Condes, Santiago RM, Chile

88 Centre for Astrophysics Research, University of Hertfordshire, College Lane, AL10 9AB, Hatfield, UK

89 University of Turin, Mathematical Department “G. Peano”, Via Carlo Alberto 10, 10123 Torino, Italy

90 INAF – Osservatorio Astronomico d’Abruzzo, Via Mentore Maggini, 64100 Teramo, Italy
91 Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, Rua do Matão, 1226, Cidade Universitaria, 05508-900 São Paulo, SP, Brazil

92 Mésocentre de calcul de Franche-Comté, Université de Franche-Comté, 16 route de Gray, 25030 Besançon Cedex, France

93 Ruđer Bošković Institute, Bijenička cesta 54, 10000 Zagreb, Croatia

94 Astrophysics Research Centre, School of Mathematics and Physics, Queen's University Belfast, Belfast BT7 1NN, UK

95 Data Science and Big Data Lab, Pablo de Olavide University, 41013, Seville, Spain

96 Institute of Astrophysics, FORTH, Crete, Greece

97 Barcelona Supercomputing Center (BSC), Plaça Eusebi Güell 1-3, 08034-Barcelona, Spain

98 ETSE Telecomunicación, Universidade de Vigo, Campus Lagoas-Marcosende, 36310 Vigo, Galicia, Spain

99 F.R.S.-FNRS, Rue d'Égmont 5, 1000 Brussels, Belgium

100 Asteroid Engineering Laboratory, Luleå University of Technology, Box 848, 981 28 Kiruna, Sweden

101 Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

102 IAC – Instituto de Astrofísica de Canarias, Via Láctea s/n, 38200 La Laguna S.C., Tenerife, Spain

103 Department of Astrophysics, University of La Laguna, Via Láctea s/n, 38200 La Laguna S.C., Tenerife, Spain

104 Astronomical Observatory, University of Warsaw, Al. Ujazdowskie 4, 00-478 Warszawa, Poland

105 Research School of Astronomy & Astrophysics, Australian National University, Cotter Road, Weston, ACT 2611, Australia

106 European Space Agency (ESA, retired), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201AZ, Noordwijk, The Netherlands

107 LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 Place Jules Janssen, 92190 Meudon, France

108 Université Rennes, CNRS, IPR (Institut de Physique de Rennes) – UMR 6251, 35000 Rennes, France

109 INAF – Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131, Napoli, Italy

110 Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, PR China

111 University of Chinese Academy of Sciences, No.19(A) Yuquan Road, Shijingshan District, Beijing 100049, PR China

112 São Paulo State University, Grupo de Dinâmica Orbital e Planetologia, CEP 12516-410, Guaratinguetá, SP, Brazil

113 Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, 2100 Copenhagen Ø, Denmark

114 DXC Technology, Retortvej 8, 2500 Valby, Denmark

115 Las Cumbres Observatory, 6740 Cortona Drive Suite 102, Goleta, CA 93117, USA

116 CIGUS CITIC, Department of Nautical Sciences and Marine Engineering, University of A Coruña, Paseo de Ronda 51, 15071 A Coruña, Spain

117 Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK

118 MTA CSFK Lendület Near-Field Cosmology Research Group, Konkoly Observatory, MTA Research Centre for Astronomy and Earth Sciences, Konkoly Thege Miklós út 15–17, 1121 Budapest, Hungary

119 Pervasive Technologies s.l., c. Saragossa 118, 08006 Barcelona, Spain

120 School of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK

121 Villanova University, Department of Astrophysics and Planetary Science, 800 E Lancaster Avenue, Villanova PA 19085, USA

122 Departamento de Física de la Tierra y Astrofísica, Universidad Complutense de Madrid, 28040 Madrid, Spain

123 INAF – Osservatorio Astronomico di Brera, via E. Bianchi 46, 23807 Merate (LC), Italy

124 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

125 Department of Particle Physics and Astrophysics, Weizmann Institute of Science, Rehovot 7610001, Israel

126 Centre de Données Astronomique de Strasbourg, Strasbourg, France

127 University of Exeter, School of Physics and Astronomy, Stocker Road, Exeter, EX2 7SJ, UK

128 Departamento de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESA-ESAC, Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

129 naXys, Department of Mathematics, University of Namur, Rue de Bruxelles 61, 5000 Namur, Belgium

130 INAF – Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, 34131, Trieste, Italy

131 Harvard-Smithsonian Center for Astrophysics, 60 Garden St., MS 15, Cambridge, MA 02138, USA

132 H H Wills Physics Laboratory, University of Bristol, Tyndall Avenue, Bristol BS8 1TL, UK

133 Escuela de Arquitectura y Politécnica – Universidad Europea de Valencia, Spain

134 Escuela Superior de Ingeniería y Tecnología – Universidad Internacional de la Rioja, Spain

135 Department of Physics and Astronomy G. Galilei, University of Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

136 Applied Physics Department, Universidade de Vigo, 36310 Vigo, Spain

137 Instituto de Física e Ciencias Aeroespaciais (IFCAE), Universidade de Vigo, Á Campus de As Lagoas, 32004 Ourense, Spain

138 Sorbonne Université, CNRS, UMR7095, Institut d'Astrophysique de Paris, 98bis bd. Arago, 75014 Paris, France

139 Institute of Mathematics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland
# 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 Fig. A.1, and four generated stellar templates are shown as examples in the right panels together with the corresponding observed RVS spectra.

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\sigma$ clipping, and got a relation of $\log_{10}(MAR) = -0.16 - 0.99 \times \log_{10}(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.

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).
<figure><img src="image_24.png" /></figure>

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

<figure><img src="image_25.png" /></figure>

**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.
## 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).

<figure><img src="image_26.png" /></figure>

**Fig. B.1.** Corner plot of the DIB fitting in the voxel with $I_{\text{pix}} = 10450$ and $d_c = 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.

<figure><img src="image_27.png" /></figure>

**Fig. B.2.** Same as Fig. B.1, but for the voxel with $I_{\text{pix}} = 10450$ and $d_c = 1.27$ kpc (the second panel in Fig. 7 from top to bottom)

<figure><img src="image_28.png" /></figure>

**Fig. B.3.** Same as Fig. B.1, but for the voxel with $I_{\text{pix}} = 10450$ and $d_c = 2.12$ kpc (the third panel in Fig. 7 from top to bottom)
<figure><img src="image_29.png" /><figcaption>**Fig. B.4.** Same as Fig. B.1, but for the voxel with $I_{pix} = 10450$ and $d_c = 2.40$ kpc (the fourth panel in Fig. 7 from top to bottom)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. B.5.** Same as Fig. B.1, but for the voxel with $I_{pix} = 10450$ and $d_c = 3.23$ kpc (the fifth panel in Fig. 7 from top to bottom)</figcaption></figure>
## 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.

```python
import numpy as np
import pandas as pd
from tqdm import tqdm
from astropy.io import fits
```

```python
dibs = pd.read_csv(somepath+'parameter_table.csv')
tabe = pd.read_csv(somepath+'spectral_table.csv')
```

```python
specLC = tabe.lc.values
specbc = tabe.bc.values
specdc = tabe/dc.values
```

```python
cataLC = dibs.lc.values
cataBC = dibs.bc.values
cataDC = dibs/dc.values
```

```python
wave = np.unique(tabe['lambda'], values)
```

```python
flux = np(zeros((dibs.shape[0], wave.shape[0]), dtype='float'))
```

```python
ferr = np(zeros_like(flux))
```

```python
for i in tqdm(range(dibs.shape[0])):
    t = (specLC == cataLC[i]) & (specBC == cataBC[i]) & (specDC == cataDC[i])
    flux[i] = data['flux'][t]
    ferr[i] = data['flux_uncertainty'][t]
```

```python
hdr0 = fits Header()
```

```python
hdr0['COMMENT1'] = 'stacked_ISM_spectra'
```

```python
hdr0['EXTNAME'] = ('flux', 'normalized_flux')
```

```python
hdr1 = fits Header()
```

```python
hdr1['EXTNAME'] = ('ferr', 'flux_uncertainty')
```

```python
hdr2 = fits Header()
```

```python
hdr2['EXTNAME'] = ('wave', 'wavelength_bins_in_vacuum')
```

```python
hdu0 = fits PrimaryHDU(flux, header=hdr0)
```

```python
hdu1 = fits.ImageHDU(ferr, header=hdr1)
```

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

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

```python
hdu.writeTo(somepath+'DIBSpec_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, Astronomie 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 (BELSPO) 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 Aperfeiçoamento 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 PIOFGA-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 Skłodowska-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 Astrofisica (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 'M'ining 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 HAR-MONIA 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 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 SSpace 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 Astrofisico 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 Astrofisica 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 Digital 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 Astrofisica, 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