The diffuse interstellar band around 8620 {\AA} I. Methods and application to the GIBS data set

We developed a set of procedures to automatically detect and measure the DIB around 8620 {\AA} (the Gaia DIB) for a wide range of temperatures. The DIB profile is fit with a Gaussian function. Specifically, the DIB feature is extracted from the spectra of late-type stars by subtracting the corresponding synthetic spectra. For early-type stars we applied a specific model based on the Gaussian process that needs no prior knowledge of the stellar parameters. The method was tested on $\sim$5000 spectra from the Giraffe Inner Bulge Survey (GIBS). After validation, we obtained 4194 reasonable fitting results from the GIBS database. An EW versus $E(J\,{-}\,K_{\rm S})$ relation is derived as $E(J\,{-}\,K_{\rm S})\,{=}\,1.875\,({\pm}\,0.152)\,{\times}\,{\rm EW}\,{-}\,0.011\,({\pm}\,0.048)$, according to $E(B\,{-}\,V)/{\rm EW}\,{=}\,2.721$, which is highly consistent with previous results toward similar sightlines. After a correction based on the VVV database for both EW and reddening, the coefficient derived from individual GIBS fields, $E(J\,{-}\,K_{\rm S})/{\rm EW}\,{=}\,1.884\,{\pm}\,0.225$, is also in perfect agreement with literature values. Based on a subsample of 1015 stars toward the Galactic center within $-3^{\circ}\,{<}\,b\,{<}\,3^{\circ}$ and $-6^{\circ}\,{<}\,l\,{<}\,3^{\circ}$, we determined a rest-frame wavelength of the Gaia DIB as 8620.55 {\AA}. A Gaussian profile is proved to be a proper and stable assumption for the Gaia DIB as no intrinsic asymmetry is found.


Introduction
Diffuse interstellar bands (DIBs) are a set of absorption features that were first discovered in 1919 (Heger 1922). These features originate in the interstellar medium (ISM; Merrill 1934Merrill , 1936 and usually contain broader widths than typical atomic lines (Herbig 1975;Hobbs et al. 2008). Herbig (1975) was the first to systematically discuss the behavior of 39 DIBs in the region of 4400-6850 Å. An extended search was made by Sanner et al. (1978) from 6500 to 8900 Å. Jenniskens & Desert (1994) made a systematic search for the DIBs on the spectra of four reddened early-type stars and presented a catalog containing 229 DIBs, of which 133 were newly detected. The total number of the DIBs increases with the quality and wavelength coverage of the spectra. The recently released Apache Point Observatory Catalog contains more than 500 DIBs covering optical and near-infrared (NIR) bands (Fan et al. 2019). More than 100 years have passed since the first discovery of DIBs, but we still know very little about their carriers. The correlation between the strength of the DIBs and interstellar extinction is a general property for many strong DIBs (Sanner et al. 1978;Lan et al. 2015). However, the lack of a linear polarization in strong DIBs (Cox et al. 2007) and their missing link to far-ultraviolet extinction (Desert et al. 1995;Xiang et al. 2017) result in the thought that large carbonaceous molecules in the gas phase rather than small dust grains are the carriers of the DIBs, for example, polycyclic aromatic hydrocarbons (PAHs, Leger & D'Hendecourt 1985;van der Zwet & Allamandola 1985;Salama et al. 1999;Cox & Spaans 2006) and fullerenes (Kroto 1988;Campbell et al. 2015). The Buckminsterfullerene, C + 60 , is the first and only identified DIB carrier for four DIBs λ9365, λ9428, λ9577, and λ9632, according to the match of the band wavelengths and the strength ratios between observational and laboratory data (Campbell et al. 2015;Campbell & Maier 2018;Lallement et al. 2018;Cordiner et al. 2019). Because DIBs are weak and easily blended with stellar lines (Kos et al. 2013), early works preferred high-quality early-type stars with only several to a few hundred observations at best. During the past ten years, the upcoming large spectroscopic surveys opened a new era in the DIB research, with a considerable number of spectra that allowed constructing a three-dimensional (3D) map of the DIBs and unveiling kinematic information and statistical properties of their carriers. Using the spectra from the Gaia-ESO Spectroscopic Survey (Gilmore et al. 2012), Chen et al. (2013) and Puspitarini et al. (2015) first detected the DIBs in the spectra of late-type stars with automated techniques by fitting the observed spectrum with a combination of a synthetic stellar spectrum, a synthetic telluric transmission, and empirical DIB profiles. In addition to the use of synthetic spectra, Kos et al. (2013) developed a method for detecting interstellar DIBs on cool-star spectra using artificial templates constructed from real spectra at high latitudes that are morphologically similar to the target spectrum. This method requires no prior knowledge of stellar parameters but can only be applied with large databases. Kos and collaborators applied the method to study the DIB around 8620 Å with ∼500,000 spectra from the Radial Velocity Experiment (RAVE; Steinmetz et al. 2006) and built a pseudo-3D map of the DIB strength covering about 3 kpc from the Sun with a spatial resolution between 0.075 and 0.8 kpc (Kos et al. 2014). Yuan & Liu (2012) also reported the detection of two optical DIBs λ5780 and λ6283 in about 2000 low-resolution spectra (R ∼ 2000) from the Sloan Digital Sky Survey (SDSS; Eisenstein et al. 2011). By stacking thousands of SDSS spectra of stars, galaxies, and quasars, Lan et al. (2015) successfully created an intensity map of 20 DIBs covering ∼5000 deg 2 and measured their correlations with various ISM tracers (atomic, molecule, and dust). The tight correlation between the strength of the DIBs and interstellar extinction was confirmed toward substantial sightlines. The strong DIB at λ = 1.527 µm (i.e., APOGEE DIB) was thoroughly studied using data from the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2016) by Zasowski et al. (2015), Elyajouri et al. (2016), and Elyajouri & Lallement (2019). In addition to the common correlation between its strength and extinction, various properties were investigated based on the large number of APOGEE spectra: Zasowski et al. (2015) derived the velocity curve of the DIB carrier and estimated the rest-frame wavelength of the APOGEE DIB; Elyajouri & Lallement (2019) revealed the depletion of the DIB carrier in dense clouds.
Based on large sky survey projects and new techniques, strong DIBs are identified to be a powerful tool for ISM tomography and consequently can probe the Galactic structure, although the carriers are unknown. The forthcoming third data release of the ESA Gaia mission that will contain the parameterization of several million spectra will be a leap forward in the sky coverage and spatial resolution of the DIB intensity map. These spectra are observed with the Gaia Radial Velocity Spectrometer (RVS; Recio-Blanco et al. 2016;Gaia Collaboration et al. 2018;Katz et al. 2019) for stars as faint as G ∼ 15.5 mag, with a spectral window from 847 to 871 nm at a resolution of ∼11,200. DIB λ8620 is also the strongest DIB covered by the Gaia-RVS spectra, known as the "Gaia DIB". It was first reported by Geary (1975) and has been widely studied for its correlation with interstellar extinction (Sanner et al. 1978;Munari 2000;Wallerstein et al. 2007;Munari et al. 2008;Kos et al. 2013;Damineli et al. 2016). Its carrier is not associated with dust grains (Cox et al. 2011) and is still not identified.
In this paper, we describe our automatic procedure for the detection and measurement of the Gaia DIB, which can be applied for large spectroscopic surveys such as the forthcoming Gaia DR3 release. We applied this method to nearly 5000 spectra from the Giraffe Inner Bulge Survey (GIBS; Zoccali et al. 2014) located in highly extincted regions. The full procedures of the DIB measurements, as well as the error analysis, are presented in Sect. 2. The GIBS data are introduced in Sect. 3. Section 4 shows the fitting results and the related discussions. Our main conclusions are summarized in Sect. 5.

Procedures of the DIB measurement
Most of the DIB studies have focused either on late-type stars (e.g., Kos et al. 2013) or on early-type stars (e.g., Munari et al. 2008) with a reasonable number of spectra to treat (several tens of thousands of stars). The challenge of this work is to implement a procedure that is valid for a wide temperature range and applicable to very large spectral surveys such as that of Gaia RVS. This requires a set of automatic procedures that has to be fast in terms of computing time and also reliable. Figure 1 shows the flowchart of our full procedures, and we describe our automatic procedures in detail below.

Inputs and spectral check
The global inputs of our procedure are the observed spectra corrected for their radial velocities together with their best-fit synthetic spectrum and the corresponding stellar parameters. We used these stellar parameters together with the corresponding synthetic spectra for stars with temperatures from 3500 K to 7000 K, which we call cool stars. We chose this limit in order to ensure that we did not encounter problems with the synthetic spectra at the border of their grid. For stars above 7000 K, which are called hot stars, we used a specific technique based on the Gaussian process that does not require synthetic spectra, as described in Kos (2017).
The input spectra and parameters (effective temperature T eff , signal-to-noise ratio S/N, radial velocity of the target star RV star , and its uncertainty σ(RV star )) were checked before further processes to eliminate invalid cases. Cool-star spectra should have nonzero flux for both observed and synthetic spectra, while for hot-star spectra, only a nonzero observed flux is required. Stars with T eff < 3500 K were discarded because those spectra are mainly dominated by molecular lines that cannot easily be reproduced well by synthetic spectra. In addition, in order to avoid fitting random-noise profiles instead of the true DIB profiles, we restricted our analysis to stars with S/N > 50. We describe in Sect. 2.6.1 the effect of the S/N on the error in the DIB measurement. RV star was used to convert the central wavelength measured in the stellar rest frame into the heliocentric frame. Targets with large radial velocity errors (σ(RV star ) > 5 km s −1 ) were discarded as well.

Interstellar spectra and renormalization
The interstellar spectra were derived by dividing the observed spectra by the corresponding synthetic spectra for cool stars. For hot stars, the Gaia DIB is usually not blended with stellar lines and can be directly measured on the observed spectra, while for cool stars, the stellar lines first need to be removed by using the  synthetic spectra. We refer here to the cool interstellar spectra and to the hot interstellar spectra as CIS and HIS, respectively. We analyzed and measured the Gaia DIB in a 35 Å wide region around its central wavelength (Jenniskens & Desert 1994;Galazutdinov et al. 2000;Munari et al. 2008), that is, 8605-8640 Å. Although the input spectra should be normalized, the interstellar spectra usually do not have uniform continua. Especially for hot-star spectra, heavily uneven continua can be found with the strong hydrogen Paschen 14 line (see Fig. 2 or the examples shown in Munari et al. 2008). Therefore, a specific renormalization technique was applied to the local spectra within the 8605-8640 Å spectral window. For HIS, the local spectrum was first fit by a second-order polynomial, where the differences of the flux of each pixel to the fitting curve were calculated, as well as their standard deviation. Pixels far away from the fitting curve were replaced by the corresponding points on the fitting curve. Specifically, for the pixels above the polynomial, they were replaced when their distances were larger than five times the standard deviation. When the pixel was below the fitting curve, the threshold was 0.5 times the standard deviation. Different rejected thresholds were set to ensure that the fitted continuum can access the real continuum and is not lowered by the stellar and/or DIB features. The remaining pixels, together with the points replacing outliers, were fit again by a second-order polynomial. After 20 iterations, the final fitted polynomial was used as the continuum to renormalize the original local spectrum. Figure 2 illustrates the local renormalization with five RAVE spectra of hot stars. The spectra and their atmospheric parameters were taken from RAVE-DR6 (Steinmetz et al. 2020a,b). The curvatures caused by the Paschen 14 line are alleviated, but the DIB and stellar features are kept.
The same technique was also applied to the cool-star spectrum, but using a linear form. The local renormalization and the derivation of CIS were made simultaneously following these steps: 1. Derive a rough interstellar spectrum, R rough = F λ /S λ , where F λ is the observed spectrum and S λ is the synthetic spectrum. F λ and S λ have the same spectral samplings. 2. Renormalize R rough and extract its continuum, F cont . 3. Renormalize the observed spectrum, F norm = F λ /F cont . 4. Derive the final interstellar spectrum: A renormalized CIS is shown in Fig. 3 with the corresponding fit of the DIB feature on it. The spectrum and stellar parameters come from RAVE-DR6 as well.

Preliminary detection
In order to process a large number of spectra, the fitting of the DIB profile was completely automated without any visual inspection. We therefore made a preliminary detection of the DIB profile to produce initial guesses for the fitting and eliminated cases whose noise is at the level of or exceeds the depth of the DIB feature. The detection was made within the wavelength range between 8614.3-8625.7 Å according to a radial velocity of ±200 km s −1 of the DIB carrier at the stellar frame. This is a reasonable assumption if the DIB carrier mainly traces the local ISM at several kiloparsecs from the Sun. When the largest depth of the spectrum in this region is larger than 3 × 1 S /N , we considered this DIB as a true detection, and the considered spectrum entered the main process of the DIB profile fitting (see Fig.  1), where the depth and its according position were used as the initial conditions of D and λ C in the DIB fitting. Otherwise, the case was discarded.

Main process: Fitting the DIB profile
The observed profile of the Gaia DIB along the line of sight could be the superposition of several features with different widths but at almost the same wavelength (Jenniskens & Desert 1994), which can be described by a Gaussian profile (Kos et al. 2013). We decided to fit the DIB feature on the spectra with a Gaussian profile because 1) previous studies revealed no intrinsic asymmetry of the Gaia DIB (Munari et al. 2008;Kos et al. 2013;Puspitarini et al. 2015), 2) the departures from a Gaussian profile caused by the multiple cloud superposition is smaller than other sources of uncertainty (Elyajouri et al. 2016), and 3) a Gaussian fit is easier, more stable, and faster in terms of computing time than the asymmetric Gaussian fit.

Models for the cool and hot Interstellar spectra
The DIB profiles on CIS and HIS were fit by different techniques that are described below. CIS was modeled using a Gaussian function that describes the DIB profile and a constant that accounts for the continuum, where D and σ are the depth and width of the DIB profile, λ C is the measured central wavelength, C is the constant continuum, and x is the spectral wavelength. However, a simple Gaussian model is not suitable for HIS because they are usually distorted by the strong Paschen 13 and 14 lines and sometimes contain a strong N I line around 8629 Å (see, e.g., HD 150850 and HD 167745 in Fig. 2). To fit the DIB profile together with the distorted continuum and possible stellar lines, we applied a similar method as in Kos (2017) using the Gaussian process (GP) described in detail below.
The GP is defined as a collection of random variables, any finite number of which have a joint multivariate Gaussian distribution (Schulz et al. 2018). Formally, let the input space be X, and f denotes a function mapping the input space to reals: f : X → R. Then, f is a GP if for any vector of in- . . , f (x n )] T is Gaussian distributed. GP is specified by a mean function m(x) reflecting the expected function value at input x, and a kernel (also called covariance function) k(x, x ) models the dependence between the output values at different input points (Schulz et al. 2018). GP can be used as a supervised learning technique for classification and regression.
Gaussian process regression (GPR) is a nonparametric Bayesian approach to regression problems (Gershman & Blei 2012). The output y of a function f at input x can be written as where ∼ N(0, σ 2 ) represents the observational error. f (x) ∼ GP(m(x), k(x, x )) is distributed as a GP (Schulz et al. 2018). GPR can capture many different relations between inputs and outputs by using a theoretically infinite number of parameters (Williams 1998).
For HIS, our goal is to apply GPR to fit the DIB profile and the remaining spectrum simultaneously. The prior mean function is often set to m(x) = 0 in order to avoid expensive posterior computations. Because we wish to extract the information of the DIB feature, however, a Gaussian mean function (Eq. 1) is applied with C ≡ 1. For the kernels, we followed the strategy of Kos (2017): the exponential-squared kernel models the stellar absorption lines, and a Matérn 3/2 kernel models the correlated noise, where a scales the kernels, and l is the characteristic width of each kernel. In principle, the fitting technique based on GP can be applied to spectra of both hot and cool stars. Because it is computationally expensive, however, we only applied it to hot-star spectra, which take only a small fraction of the substantial spectra in large spectroscopic surveys such as Gaia RVS. Nevertheless, as an illustration of this method, we applied it to GIBS data (see Sect. 3).

Parameter optimization and MCMC fit
Maximum likelihood estimation was used to optimize the parameters in the Gaussian model for CIS, that is, Θ = {D, λ C , σ, C}. Given the spectrum {X, y, σ 2 y }, where X is the wavelength, y is the flux, and σ 2 y is the observational uncertainties (if σ 2 y is not accessible, it was fixed to 0.001), the log marginal likelihood is where r = y − f Θ (X) is the residual vector and f Θ is the Gaussian model. N is the pixel size of the spectrum. K is the covariance matrix, To implement GPR for HIS, we optimized five parameters, three for the DIB profile (D, λ C , σ), and two for the kernels (l se and l m3/2 ). The scaling factor a of the kernel can be estimated as the variance of the noise and does not need to be fit. We used the square of the inverse of the S/N to approximate it. The optimal parameters were estimated by maximizing the type II maximum likelihood (Rasmussen & Williams 2006). Its log marginal likelihood is almost the same as Eq. 5, but the covariance matrix becomes nondiagonal, where σ i is the observational error, δ i j is the Kronecker delta, and k(x i , x j ) is the element of the specified kernel. A Markov chain Monte Carlo (MCMC) procedure (Foreman-Mackey et al. 2013) was performed to implement the parameter estimates for the Gaussian fit and GPR. The initial conditions are perturbed by a normal distribution around the initial guess with a standard deviation of 0.01. Different walkers of the MCMC can therefore start with different conditions. One hundred walkers were progressed for 50 steps to complete the burn-in stage. The best fits were then used as the initial conditions to sample the posterior with 100 walkers and 200 steps. The best estimate and its statistical uncertainty were taken in terms of the 50th, 16th, and 84th percentiles of the posterior distribution.
The initial conditions of D and λ C were measured by the preliminary detection (see Sect. 2.3). The initial values of l se and l m3/2 were set to 0.3 and 0.15 for all the cases. C has an initial guess of 1.0 assuming a well-normalized CIS. The initial guess of σ is hard to determine. Based on the GIBS results (see Sect. 4), σ 0 = 1.2 is a proper guess. Strong DIB profiles are not sensitive to the initial guess, while weak profiles in general show a good fitting behavior with this value. Examples of the DIB fittings for CIS and HIS are shown in Fig. 3 and 4, respectively. We indicate the first initial fit and a second fit for the error analysis (see Sect. 2.6.2). The final selected fits are marked as solid lines.

Priors
As a Bayesian approach, priors can be used to prevent unphysical or unreasonable fittings. For the Gaussian fit applied to the CIS, we adopted flat priors, More rigorous priors are needed for HIS to avoid treating the DIB profile as the correlated noise and fit by the kernels of GP. The priors of l se and l m3/2 provided by Kos (2017) were taken. 0.22 Å was assumed as the boundary of the characteristic widths of stellar feature and the random noise, and therefore it is the lower limit of l se and the upper limit of l m3/2 . l m3/2 has a lower limit of 0.08 Å. l se is flat at high values and gradually decreases to −∞ at the value of the DIB width. D has a flat prior the same as for CIS. The priors of λ C and σ are the Gaussian priors centered at their initial conditions with a width of 0.5 Å, which are simpler than those of Kos (2017) because they lack the preliminary fit of the DIB profile. The prior of λ C is stricter than that of σ because its initial guess can be determined by the preliminary detection. Some examples of the priors of l se , l m3/2 , λ C , and σ are presented in Fig. 5.

Quality flags
To select reliable DIB profiles, Elyajouri et al. (2016) applied a series of tests to the fit parameters (D, λ C , σ) and generated different quality flags (QF). We follow their main principles. The flowchart of QF is schematically shown in Fig. 6, and below we describe in detail our procedure to determine QF ranging from QF = 5 (highest quality) to QF = 0 (lowest quality). Cases with negative QF were not fit, that is, QF = −1 was rejected by the preliminary detection, and QF = −2 means invalid spectra.
1. Global test: The first test gives the upper limit of the depth D and the realistic range of the measured central wavelength λ C . Here λ C is converted from the stellar frame into the heliocentric frame using RV star . The cases with D > 0.15 (the deepest absorption detected on GIBS spectra) was eliminated (arrow (a) in Fig. 6). These are spurious features generated mainly by the mismatch between the observed and synthetic spectra. We also eliminated the fittings with λ C outside the range 8614.3-8625.7 Å, which is the same interval as we applied in the preliminary detection. 2. Test of the DIB depth: When the first test was passed successfully (arrow (b)), we compared their depth with the standard deviation of the fitting residuals, R = std (data-model). R was calculated for two regions: R A is for the global spectrum [8605-8640] Å, and R B is in a region close to the DIB feature [λ C − 3σ, λ C + 3σ] Å. When D was larger than the maximum of R A and R B (arrow (d)), then the test on the width was directly applied. When the interstellar spectrum is globally too noisy to detect the DIB or the DIB is too shallow (arrow (c)), we only compared D with the local standard deviation R B . This step allowed us to recover the DIB on the spectra that are noisy in some regions far from the DIB, but have good quality near the DIB. The cases passing this test (arrow (e)) were subjected to the same test on width as the previous ones (arrow (d)). Failed cases (arrow (g)) were examined differently. 3. Test of the DIB width: In the final test, we defined some limits to select DIBs with reasonable widths. The profiles that exceed the global or local noise level (arrow (d) and (e)) gain high QF with 1.2 σ 3.2 Å or low QF with 0.6 σ 1.2 Å. The shallow DIBs (arrow (g)) were directly tested within the range of 0.6-1.2 Å. Any case with σ < 0.6 Å was discarded (arrow (h)) and marked as QF = 0. Both the lower (0.6 Å) and upper (3.2 Å) limits were derived from the GIBS results (see Sect. 4.2.3). Profiles with σ < 0.6 Å were likely to come from random noise. Extremely broad profiles (σ > 3.2 Å) are due to unphysical features originating from the data processing or the fitting process. We set 1.2 Å as the boundary of the two ranges of σ to 1) select narrow DIBs through arrow (f) and 2) eliminate the flat and elongated features of uncertain origin (Elyajouri et al. 2016) for shallow DIBs (arrow (g)).
We created six QFs based on the tests to evaluate the fit quality. QF = 5 represents the best fits and the detected DIBs with proper parameters {D, λ C , σ}. Recovered DIBs with QF = 4 are locally detected. Narrow DIBs with 0.6 σ 1.2 Å are flagged  as 3 or 2 because they exceed the global or local noise level. QF = 1 corresponds to spectra with very low S/N or shallow DIBs. A failed detection is marked as QF = 0.

Equivalent width and error analysis
The equivalent width (EW) is proportional to the column density of the DIB carriers and reflects the relative oscillator strength (Jenniskens & Desert 1994). With the Gaussian profile, the EW is calculated by the depth D and width σ, where I 0 and I λ are fluxes of the continuum and spectrum, respectively. For CIS, the calculated EW was further scaled by C because the fit C is usually not unit. There are two main sources of the EW errors, σ EW , one associated with the random noise (σ noise ), and the other (σ spect ) contributed by the continuum (HIS) or the mismatches between observed and synthetic spectra (CIS). The total error is considered as We estimated σ noise for different DIB profiles by a randomnoise simulation. σ spect was accessed through a second fit for each interstellar spectrum. In the following, the estimation of these two errors is explained in detail.

Random-noise simulation
The random noise discussed in this section mainly refers to the observational uncertainty, which is Gaussian and independent. The random-noise level of a spectrum is usually characterized by its S/N. The uncertainty introduced during the data reduction and interstellar spectra derivation is discussed and estimated in Sect. 2.6.2. Although the random noise is assumed to be Gaussian, the local noise might still distort the DIB profile and affect the fit parameters and consequently the physical quantities such as EW and radial velocity. To account for the error of EW contributed by the random noise, Elyajouri et al. (2017) made a conservative estimation: σ noise = 2 √ 2 σ δ depth , where σ is the fit width of the DIB profile, and δ depth is the uncertainty of the DIB depth. Puspitarini et al. (2015) applied a similar formula with a scaling factor 1 √ N , where N is the number of pixels covering the DIB width. Their formulas were derived from a series of simulations with varying Gaussian noise and can quickly approximate σ noise . The use of σ in their formulas might lead to a strong overestimation for large DIB profiles, however.
For a more comprehensive study and more accurate estimate of the effect of the random noise on the DIB fitting, we performed a series of random-noise simulation in the wavelength range between 8605 Å to 8640 Å with a pixel size of 0.1 Å, containing different Gaussian DIB profiles ({D 0 , µ 0 , σ 0 }) and constant continua (C 0 ≡ 1). Then for every spectrum, a Gaussian noise ( ) was added according to an assigned S/N, that is, ∼ N(0, (S/N) −2 ). The parameter grids were constructed as D 0 ranges from 0.01 to 0.20 with a step of 0.01, µ 0 ≡ 8620.0 Å, and σ 0 ranges from 0.05 to 5.0 Å with a step of 0.05 Å. For 20 S/N 100, the step size is 1. For an S/N within 100-300, the step size 5. Higher S/N (300-1000) were assigned a step size of 50. Finally, the sample contains 270,000 pseudo-spectra in total with different DIB profiles and S/N. These spectra were fit by a Gaussian model with the Levenberg-Marquardt method, and the fit parameters {D f , µ f , σ f , C f } were used to study the effect of the random noise. The true (EW 0 ) and fit (EW f ) EWs were calculated with Eq. 9. Figure 7 shows the distribution of the fractional error between EW 0 and EW f (|EW 0 −EW f |/EW 0 ) in the D 0 −σ 0 plane for some specific S/N, overlapped with some contours of EW 0 calculated by the according D 0 and σ 0 . The change in fractional errors is rough because we only fit each spectrum once. The shown σ 0 is limited to 3.2, the same as the largest valid width detected on GIBS spectra. Generally, the fractional error decreases with the increase in S/N and EW 0 . For S/N > 200, the fractional errors are smaller than 10% (white regions in subpanels in Fig. 7) for most of the fitting results, but for S/N = 100, EW 0 has to be as large as 0.4 Å to ensure that most of the fractional errors are within 10%. Nevertheless, some shallow profiles could still gain large errors up to 20%. If the spectra have S/N ≈ 50, the random noise can cause fractional errors as large as 20% even for EW 0 > 0.5 Å, which is stronger than most of the Gaia DIBs detected in previous works (Sanner et al. 1978;Munari et al. 2008;Puspitarini et al. 2015). Therefore we only regard sources with S/N higher than 50. Moreover, for a given EW 0 , shallow DIB profiles tend to gain larger errors than narrow ones because the shallow profiles cover more pixels. If σ 0 is approximate to the pixel size of the spectra, the fractional error maintains a high level and does not significantly decrease with S/N. This also occurs when D 0 is approximate to 1 S/N . This implies the detection limits for the width and depth of DIB. The effects of the random noise on D 0 and σ 0 are similar to that on EW, while the error of µ 0 is more sensitive to D 0 than σ 0 . The random noise has almost no effect on the continuum C 0 for well-normalized spectra.
Based on the random-noise simulation, the effect of random noise on the DIB fitting was studied in detail. Estimating of σ noise is not straightforward, however, because in practice we can only access the fit parameters and not their true values. That is to say, we have to use {D, µ, σ} instead of {D 0 , µ 0 , σ 0 } to estimate the error of EW f contributed by the random noise, that is, σ noise . We try to build a model based on the random forest regression, which is an ensemble machine-learning method combining a large number of decision trees (Breiman 2001). The model returns σ noise when given {D, σ, S/N}. µ and C were not used because in our simulation µ 0 and C 0 were fixed. A quarter of the simulation results that were uniformly selected with EW 0 constituted the training set, and the test set consisted of the remaining part. The regression was completed by the Python scikit-learn package (Pedregosa et al. 2011). We used 100 trees in the forest (n_estimators=100) and followed the default values of other main parameters. The differences between the true and estimated σ noise are mainly within 0.05 Å for the training set and 0.1 Å for the test set, and they do not significantly change with EW 0 . The uncertainty for S/N < 50 could be up to 0.2 Å and higher. The performance of the model is limited by the fact that the features we used {D, σ, S/N} in the algorithm are not enough to fully access the true σ noise . The estimate of σ noise is accurate for large DIB or high-quality spectra, but it is less reliable for small DIB or low-S/N spectra.

Spectral contribution
To obtain σ spect , each spectrum was fit twice. The first fit is detailed in Section 2.4. The second fit considered the effect of the observed-synthetic mismatch for CIS and continuum for HIS. The difference of EWs between these two fits was used to estimate σ spect .
For CIS, the second fit is still a Gaussian fit, but with five masked regions centered at 8611.8 Å, 8616.3 Å, 8621.6 Å, 8626.2 Å, and 8634.1 Å, with a width of 1 Å for each of them. These regions correspond to some strong stellar lines, for instance, the Fe I lines at 8610.602 Å, 8611.804 Å, 8613.935 Å, 8616.276 Å, 8621.601 Å, and 8632.412 Å; the Ca I line at 8633.933 Å; and the Ti I line at 8618.425 Å (the strength of these stellar lines would vary with the stellar types and metallicities), that may be poorly modeled by the synthetic spectra. Figure 3 shows that the mismatches in the masked regions are higher than average. Consequently, a large σ spect is obtained. Although this method is incomplete because we cannot mask all of the abundant stellar lines in the DIB analysis interval, it is still a good estimate of σ spect for strong DIBs, as discussed by Puspitarini et al. (2015).
Although the local renormalization corrects for the curved continuum of HIS, some curvatures could still remain and lead to an underestimated EW. After extracting the fitted DIB profile, we therefore applied polynomials from first to sixth order to fit the remaining HIS to approximate the possible curved continuum. Then the best fit was used to renormalize HIS again, and we refit the DIB profile with the Gaussian model (Eq. 1).
To ensure that the second fit is reasonable, the new DIB profile was preferred only if it was stronger and deeper than the first. Otherwise, we retained the results from the first fit. For example, the second fit was accepted for star HD 149349 (Fig. 3), while it was rejected for star HD 166167 (Fig. 4). Furthermore, we obtained an EW of 0.234 ± 0.039 Å for star HD 166167, which is consistent with the value of 0.217 Å reported by Munari et al. (2008), who applied a sixth-order polynomial to fit the continuum and then calculated the EW by integration. Additionally, the cases with large differences between the measured central wavelength of the two fits were eliminated, that is, ∆λ C > 0.5 Å (the sixth step in Fig. 1). These cases were also marked QF = 0.

Outputs and summary
The final output of each fitting was included in the fit parameters ( {D, λ C , σ, C} for CIS, {D, λ C , σ, l se , l m3/2 } for HIS), calculated EW and its errors (σ EW , σ noise , σ spect ), and QF. For discarded cases, all the parameters were set as -1.

GIBS spectra
The GIBS is a survey of red clump (RC) stars selected from the Vista Variables in the Via Lactea (VVV) catalogs (Minniti et al. 2010) in the Milky Way bulge (Zoccali et al. 2014). We used 4797 low-resolution spectra from 20 observational fields in the GIBS survey (see Fig. 8). The GIBS spectra analyzed here are from the GIRAFFE LR8 setup at the resolution R = 6500 with the spectral coverage of 8206 Å < λ < 9400 Å. For the analysis, we selected a smaller range of 8450-8950 Å because beyond the shorter interval adopted for the analysis, many skylines affect the spectra and we are interested in the region around the calcium triplet lines where the Gaia DIB is located. We calculated the S/N of each spectra between 8850 Å and 8858.5 Å, where no strong stellar lines are present. We also used the templates of RCs to subtract stellar components from the DIB measurement. The synthetic spectra used in this work were generated by the Turbospectrum code (Alvarez & Plez 1998), the MARCS atmosphere model (Gustafsson et al. 2008), and the line list of the Gaia-ESO survey (Heiter 2020, in prep) with T eff = 4500 K, log g = 2.5 and the particular metallicity of each star. As the residuals in the interstellar spectra are too high because of the mismatch between observed and synthetic spectra, we did not apply the CIS method to the GIBS data set.

Extinction
Based on the VVV survey, Gonzalez et al. (2011Gonzalez et al. ( , 2012 built the first complete bulge extinction map (G12 henceforth) with a differential method. They first derived the mean (J−K S ) color of the RC stars in 1835 subfields and then compared it to the color of RCs in a referred region with known extinction, that is, Baade's Window. The extinction calculated with the BEAM calculator 1 is shown in Fig. 8. The resolution varies from 2 to 6 . With a newly developed JK S photometry catalog of the VVV survey  Surot et al. (2020) calculated the mean (J −K S ) of RC+RGB (red giant branch) stars in finer bins. A calibration was then made by comparison with G12 in areas with |b| > 3 • to derive the absolute extinction values. The improved extinction map (S20 henceforth) has a higher resolution that can reach subarcmin in lowlatitude regions.
We derived the extinctions of the GIBS targets from S20 according to their spatial positions. Neither S20 nor G12 were able to resolve the extinction of individual GIBS targets for |b| > 3 • because the spatial resolution of these maps decrease with increasing latitudes.

Specific model for GIBS spectra
We used the GIBS spectra to validate our procedures of the DIB detection and measurement presented in Sect. 2 and also test the fitting technique based on GPR. We applied the GPR method to fit the DIB profile on the derived interstellar spectra from GIBS because 1) the data size is small, therefore the computational time is acceptable, and 2) the templates with constant T eff and log g cause considerable mismatch and correlated noise that prevent fitting the simple Gaussian model. The simple Gaussian model has been widely used to fit various DIB profiles on a substantial number of spectra, for example, Kos et al. (2013), Lan et al. (2015), Zasowski et al. (2015), and Elyajouri et al. (2017). We therefore did not select a specific sample to test the Gaussian fit used for cool-star spectra. An illustration is shown in Fig. 3.
Furthermore, we applied both Gaussian and asymmetric Gaussian models to the GIBS spectra to study the amplitude and the effect of the asymmetry caused by the velocity dispersion of DIB carriers along different sightlines. We chose a simple method presented in Kos (2017) to implement the asymmetry, which was introduced by making the width a function of the wavelength, where asym is the amplitude of the asymmetry, and it is limited within ±0.75 in the fittings.

Results and discussions
In this section, we study the correlation between EW and reddening, as well as the properties of the Gaia DIB based on the fitting results. In general, the Gaussian fit and the asymmetric Gaussian fit yielded similar profiles (see Section 4.2.1). We therefore base our results and discussions on the Gaussian fit alone. Because of the applied synthetic spectra, the local renormalization for the GIBS spectra could lead to nonunit continua and cause an overestimation of the EW. We therefore performed a calibration to the EW according to the residual spectra (datamodel). The fit profiles were shifted according to the mean flux of the residual spectra where the EW were recalculated as well. Most of the mean flux is within 0.95-1.0, implying an overestimation of the EW. For 378 cases the local continua are above 1.0. After visual inspection, they were eliminated because the fit profiles were not physical or too noisy. The cases that did not pass the preliminary detection were also discarded. Finally, we obtained 4194 valid GIBS spectra for our DIB analysis. As the average S/N of GIBS spectra is about 80 per pixel, the fit parameters (D, λ C , σ) and EW are noise dominated, especially for small DIBs.

Linear correlation between EW and reddening
Several studies have revealed a linear correlation between EW and reddening in the optical bands for the Gaia DIB, usually taking the form of E(B − V) = a × EW. Some early estimates of a are based on several dozen hot stars, for example, 2.85 (Sanner et al. 1978) (the coefficient was calculated by Kos et al. (2013)), 2.69 (Munari 2000), 4.61 (Wallerstein et al. 2007), and 2.72 (Munari et al. 2008). The measurement of E(B − V) in Wallerstein et al. (2007) is doubtful, as discussed in Munari et al. (2008). By merging several thousand RAVE cool-star spectra, Kos et al. (2013) derived a = 2.49 with an offset of 0.028. Kos et al. (2013) also studied 114 RAVE hot stars (including 31 objects from Munari et al. (2008)), which yielded a highly consistent value of 2.48. Puspitarini et al. (2015) also derived a linear relation between EW of the Gaia DIB and A 0 (extinction at λ = 5500 Å) toward the Galactic anticenter based on 64 cool stars from Gaia-ESO, but no coefficient was given. From Fig. 7 in Puspitarini et al. (2015), we estimate a coefficient of A 0 /EW 8620 = 2.3/0.35 = 6.57. Applying the CCM89 model (Cardelli et al. 1989), we have E(B − V)/EW = 2.12 (R V = 3.1), a value significantly lower than others. While Damineli et al. (2016) reported a quadratic relation between A K S and EW based on ∼100 hot stars in and around the stellar cluster Westerlund 1. Their relation is close to Munari et al. (2008) for EW < 0.5 Å. As discussed below, we derived a linear correlation between EW and E(J − K S ) for our GIBS sample by averaging over different reddening bins (Sect. 4.1.1) or individual fields (Sect. 4.1.2).

EW versus E(J − K S ) relation
The correlation between the EW of the Gaia DIB and the reddening E(J − K S ) from S20 for individual GIBS targets is shown in Fig. 9, overlapped with the measurements from Munari et al. (2008). We derived the linear relation over the GIBS targets by taking the median values from different reddening bins, ranging from E(J − K S ) = 0.2 to 0.9 with a bin size of 0.1 mag (the white and red dots in Fig. 9). The linear fitting for these median points yields E(J − K S ) = 1.875 (± 0.152) × EW − 0.011 (± 0.048). The coefficients and their standard errors were derived with the Python package statsmodels (Seabold & Perktold 2010).
The median points in bins with E(J − K S ) < 0.2 mag were not used because they deviate significantly from the linear relation. This is because a) spectra containing very small DIBs (EW < 0.05 Å before calibration) did not pass the preliminary detection because of their low S/N. Therefore the median value of EW is higher than expected. b) Most of the small DIBs are noise dominated and might be overestimated by the local increase of noise. It is not possible to select and discard the overestimated cases even by visual inspection. Reliable fittings with EW ≈ 0.1-0.2 Å can also be found in low-extinction regions, however.
In the absence of optical photometry and in order to compare our results with previous works, a conversion between E(J − K S ) and E(B − V) is needed. For the NIR bands, we applied the extinction law derived by Nishiyama et al. (2009) toward the Galactic center: 528, which is widely used for Galactic Bulge studies, while for the optical-NIR conversion, we still used the CCM89 model with R V = 3.1 and the corresponding ratio Gonzalez et al. (2011Gonzalez et al. ( , 2012 used the same ratio to calculate E(J − K S ) for RC stars, which has been used by Surot et al. (2020) to calibrate their extinction map. Therefore E(J − K S ) used in this work already implies a specific ratio between A K S and E(B − V). 2) Although many works have pointed out that the extinction law toward the inner Milky Way deviates from the CCM89 model with R V = 3.1 (e.g., Indebetouw et al. 2005;Nishiyama et al. 2006Nishiyama et al. , 2009Nataf et al. 2016;Damineli et al. 2016), the optical-NIR relation is not studied as well as for NIR bands, and the ratio of  Fig. 9). The comparison with other works is shown in Fig. 10, and the coefficients for the linear relation in the NIR bands are listed in Table 1.
The choice of extinction law significantly affects the comparison. Assuming the CCM89 model in NIR bands (R V = 3.1), With the extinction laws derived by Damineli et al. (2016) from optics to NIR, EW = 2.833. Although without the optical photometry for GIBS targets, the high consistency between optical and NIR relations against the Gaia DIB still needs to be confirmed by more studies, our comparison and the consistency found in this work and Damineli et al. (2016) imply that the correlations of EW with optical and NIR extinctions are at least not very far away from each other. However, the relations between EW 8620 and extinction in different bands and the extinction law have been studied very little. On the other hand, the correlation between EW and extinction is also related to the dust properties along the line of sight. Ramírez-Tannus et al. (2018) reported a linear correlation between extinction-normalized EW, EW/A V , and R −1 V for 14 DIBs in M17. They derived a relation between EW and E(B − V). Li et al. (2019) used E(B − V) to normalize the EW and reported no relation with R −1 V . They suggested that the hydrogen column density, N H , is a more appropriate normalization than extinction and reddening. Theoretically, small dust grains would present a steep extinction curve with small R V , and very large grains experience a flat curve with R V → ∞ (Draine 2003). This means that different R V values should indicate different cor-  Sanner et al. (1978) 1.964 0.11 Munari (2000) 1.853 0.03 Wallerstein et al. (2007) 3.176 0.56 Munari et al. (2008) 1.874 0.03 Kos et al. (2013) 1.716 0.23 Puspitarini et al. (2015) 1.461 -

Notes.
(a) reddening bins (b) individual fields (c) after correction relations between EW and extinction, especially in the ultraviolet and optical bands. Although the variation in R V in a single region is always smaller than the uncertainty in EW and extinction, we could investigate it from different sightlines: the linear coefficient between EW and E(J − K S ) toward the Galactic Bulge (this work and Munari 2000; Munari et al. 2008) is apparently larger than the value toward the Galactic anticenter (Puspitarini et al. 2015), although this difference might also be caused by the dependence of the extinction laws on different lines of sight. Kos et al. (2013) derived a mediate value from a substantial number of sightlines. Nevertheless, our studies are also affected by the method of extinction calculation and DIB measurement, which undermines the credibility of the variation of EWextinction correlation for different sightlines. On the other hand, the environmental dependence of the DIB carriers complicates this question as well. For example, the deviation from linear relation in high-extinction regions (A V ≈ 10 mag) may be caused by the carrier depletion in dense cores (Elyajouri & Lallement 2019). The forthcoming Gaia-RVS spectra are expected to bring new insights into the DIB properties. Its large spatial coverage and uniform DIB and extinction measurements will give us an opportunity to unveil the relation between the DIB strength and the corresponding dust properties.

EW versus E(J − K S ) relation from individual GIBS fields
We decreased the EW dispersion by averaging the measurements by calculating the median values of the EW and E(J − K S ) for each GIBS field individually. The correlation between EW and E(J − K S ) derived from the 20 GIBS fields is shown in Fig.  11 (left panel). After discarding fields with E(J − K S ) < 0.1 (indicated by the dot-dashed green line), we gain a linear relation of E(J − K S ) = 1.802 (± 0.258) × EW + 0.004 (± 0.065), corresponding to E(B − V) EW = 2.615, which is slightly smaller than the value derived before (2.721). The EW dispersion in each field does not notably decrease compared to that in the reddening bins. A possible reason is that the low spatial resolution of E(J − K S ) obscures the environmental variation in each GIBS field, which is traced by the Gaia DIB, leading to large dispersion of EW but small dispersion of E(J − K S ) (the dispersion of E(J − K S ) is usually smaller than its uncertainty of individual targets). This cannot account for the large dispersion in the fields with large EW, however, especially for the field (l, b) = (8, −2) (indicated by "F1" in Fig. 11). The problem may come from the contamination of stars in the GIBS target selection that are not RC stars. The RC stars in the GIBS survey are selected based on the J versus (J − K S ) color-magnitude diagram (CMD) with a limit in the J magnitude and a lower cut of J − K S , while J − K S is not stringently constrained at the red end (see Fig. 3 in Zoccali et al. 2014). Therefore the RC sample might be contaminated by highly reddened dwarfs and/or RGB stars. The spectra of the contaminators that might be very different from the RC template may give rise to pseudo-features on the interstellar spectra, causing incorrect fittings and calculations of EW. Although E(J − K S ) we used for the targets come from S20, which is not sensitive to stellar type, their map was calibrated by Gonzalez et al. (2012), which was based on RC stars. This means that contaminators far away from the peak color (J − K S ) will also gain incorrect E(J − K S ) values.
We therefore performed the correction for both reddening and EW based on VVV-DR2 catalog 2 (Minniti et al. 2017). We constructed a purer sample of RC stars by applying an additional color cut. For each field, RC candidates were first selected from the VVV catalog in a circular region located at the field center with a radius of 0.5 deg. Then we fit the (J − K S ) colors within the range of the J magnitudes given by the GIBS targets (see the dashed orange lines in Fig. A.1) with a Gaussian function to obtain the peak color as well as the 1σ width. This criterion ensures that our RC sample is as pure as possible, with the disadvantage that we loose stars. The percentage of the rejected stars differs for different fields. In total, we obtained 2437 targets in the purer sample, compared to 4194 in original sample. This means that about 42% stars are discarded. Assuming an intrinsic color of the RC stars of (J − K S ) 0 = 0.674 (Gonzalez et al. 2011), we obtain an average E(J − K S ) for each GIBS field where the standard deviation of J − K S is treated as the uncertainty of E(J − K S ), ∆E(J − K S ), including not only the error of E(J − K S ) calculated by RC stars, but also the dispersion of E(J − K S ) in each field. For low-reddening fields, ∆E(J − K S ) are similar to the mean errors of E(J − K S ) given by S20. The average value of ∆E(J − K S ), 0.066, is also close to the resultant error of RC stars when a photometric error of 0.03 for the J and K S bands and a spread of (J − K S ) 0 of 0.03 is assumed. The increase in ∆E(J − K S ) in highly extincted fields is caused by the J-magnitude range applied for GIBS targets. A wider Jmagnitude range covers a wider range of extinction and results in a larger dispersion of E(J − K S ). ∆E(J − K S ) is consequently dominated by the dispersion and can reach values above 0.1 mag.
Furthermore, the median EW of each field was recalculated by only considering the targets with J − K S within the 1σ region (dashed green lines in Fig. A.1). The correlation between EW and E(J − K S ) for individual fields after correction is presented in the right panel in Fig. 11. The dispersion of EW in each field markedly decreases and now is comparable to the uncertainty of E(J − K S ). A tighter linearity of the correlation was also derived by considering fields with E(J − K S ) > 0.25 with a Pearson correlation coefficient of 0.95, compared to the value 0.88 before the correction. The coefficient of the linear relation is E(J − K S ) = 1.884 (± 0.225) × EW − 0.012 (± 0.072), corresponding to E(B − V) EW = 2.734, which is highly consistent with Munari et al. (2008) as well. Using our new purer sample (see Fig. A.2), we obtain a similar relation (as discussed in detail in Appendix A). The relative strength (EW/A V ≈ 0.1 Å mag −1 ) and the tight correlation with reddening confirm that Gaia DIB is a powerful tracer of ISM species, independent of the foreground extinction, as suggested by Kos et al. (2013).
Four fields were selected to demonstrate the correction effect: F1 (l, b) = (8, −2), F2 (l, b) = (0, −1), F3 (l, b) = (−5, −2), and F4 (l, b) = (−4, −6). They are indicated in Fig. 11. The applied correction performs very well for highly extincted regions. However, it intensifies the deviation of the fields with low extinctions because a) the S/N of the GIBS target limits the detection of small DIBS (see also  extinction is lower, J − K S is not very sensitive to small-scale variation in the reddening. Optical data such as B − V would be more sensitive.

Asymmetry of the DIB profile
Different DIBs have diverse profiles, from single profiles (Sarre et al. 1995) to resolved substructures (Galazutdinov et al. 2002). Asymmetric shapes originate from the distortion of unresolved substructures or blended DIBs (Kos 2017). The Gaia DIB is suggested to be blended by two DIBs (Jenniskens & Desert 1994), while no intrinsic asymmetry has been unveiled by previous works. For most of the GIBS results, the difference between Gaussian and asymmetric Gaussian EW is comparable to their uncertainties (Fig. 12). As shown in Fig. 13, large asymmetries only occur in small DIBs, which probably originate in noise. No signature of intrinsic asymmetry of the Gaia DIB is revealed by Figs. 12 and 13 for large or small EWs.

Rest-frame wavelength
The rest-frame wavelength of the Gaia DIB is reported as 8620.8 Å by Galazutdinov et al. (2000) from one single star, and 8621.2 Å by Jenniskens & Desert (1994) from four hot stars, and 8620.4 Å by Munari (2000) and Munari et al. (2008) from dozens of RAVE hot stars. The determination of Munari et al. (2008) was based on the assumption that the average velocity of their carriers, which are close to the Galactic center, is essentially zero, after adopting the ISM radial velocity map of Brand & Blitz (1993). Therefore the average central wavelength λ C represents the rest-frame wavelength. We followed this method and selected spectra from the GIBS fields with −3 • < b < 3 • and −6 • < l < 3 • . Finally, we obtained 1015 spectra. The observed λ C is in the stellar frame, and it has to be converted into the heliocentric frame based on the stellar radial velocity. The converted λ C presents a Gaussian distribution (Fig. 14,  but the large uncertainty makes it still not sufficiently definite. A more accurate method is investigating measurements toward the Galactic anticenter, as illustrated in Zasowski et al. (2015), which we intend to apply to the Gaia-RVS spectra. Figure 14 shows the distributions of the measured DIB parameters from the GIBS spectra: depth D versus EW, line center λ C , and width σ. Unreasonable results with very small or large measurements were eliminated. Small EWs increase with depth, while the relation deviates from linearity when EW > 0.2 Å. Profiles with D > 0.15 come from spectra without proper normalization and were discarded. The limit of D in the QF test (Sect. 2.5) was therefore set as 0.15. The peak value of λ C for all the GIBS targets is the same as that of the subsample for deriving the rest-frame wavelength. The distribution of the subsample can be well fit by a Gaussian function. However, the λ C distribution of all the GIBS targets apparently deviates from the Gaussian profile and contains a bump around 8624 Å. The origin of this second bump is not clear and will be studied in a forthcoming paper.

Parameter distributions
The valid widths are within 0.6-3.2 Å, with a peak value of 1.74 Å. This value is close to the peak value of the APOGEE DIB λ1.5273 derived by Zasowski et al. (2015). These two DIBs also have a similar relative strength. An abrupt decrease in σ occurs at ∼1.4 Å; it is larger than that of the APOGEE DIB but consistent with APOGEE measurements for sightlines with l < 40 deg. As the inner Galaxy generally contains broader features than the outer Galaxy (Zasowski et al. 2015), we set a lower value σ = 1.2 for the QF test to distinguish narrow and shallow DIBs. We did not find any DIB width larger than 3.2 Å but a steep decrease from 2 to 3. However, Zasowski et al. (2015) revealed that the APOGEE DIB can be as broad as 6 Å toward the Galactic center. The upper limit of σ in the QF test could change if we were to find physical DIB profiles with larger width in further studies.

Conclusions
The main goal of this work was to develop a procedure for the automatic detection and measurement of the Gaia DIB (λ ≈ 8620 Å). A preliminary detection was applied to exclude low-S/N spectra and/or DIBs below the detection limit. The DIB feature was extracted from the cool-star spectra using synthetic spectra, while for hot stars, we applied a specific model based on GP (Kos 2017) to directly measure the DIB feature on the observed spectra without any stellar templates. The DIB profile was fit by a Gaussian function, and the EW was also calculated. A simulation based on pseudo-spectra with different S/N illustrated the effect of the random noise on the DIB fitting, as well as the EW calculation. Based on these simulations, a minimum S/N of 50 is required to detect DIBs. The error contributed by the synthetic spectra for cool stars and local continua for hot stars was also considered through a second fit. Furthermore, some tests on fitted parameters, {D, λ C , σ}, similar to Elyajouri et al. (2016), were used to assess their qualities.
These procedures and techniques were applied on a sample of 4979 GIBS spectra. The main results are summarized below.
By taking the median values from different reddening bins, we derived a linear relation between the EW of the Gaia DIB and the reddening: E(J − K S ) = 1.875 (± 0.152) × EW − 0.011 (± 0.048).
Applying the CCM89 model and the NIR extinction laws toward Galactic center from Nishiyama et al. (2009), we find E(B − V)/EW = 2.721, which is highly consistent with the results of Munari (2000) (2.69) and Munari et al. (2008) (2.72). Additionally, the difference of the coefficient between our result and other studies with different sightlines implies a possible variation in the relation for different ISM conditions. The median measurements from individual GIBS fields presented a relation with a coefficient of E(J − K S )/EW = 1.802 ± 0.258, with a relatively large dispersion of the EW in each field due to the contamination of non-RC stars in the GIBS sample. We eliminated them by using an additional color cut. This led to a smaller dispersion and improved the linearity of the EW-E(J − K S ) correlation for individual fields. The corrected relation, E(J − K S ) = 1.884 (± 0.225) × EW − 0.012 (± 0.072), also compares well with other results.
Assuming that the average radial velocity of the DIB carrier is zero when they are distributed close to the Galactic center (Brand & Blitz 1993), we determined the rest-frame wavelength of the Gaia DIB as λ 0 = 8620.55 ± 0.55 Å.
We also fit the GIBS spectra with an asymmetric Gaussian model. The results are in general consistent with those from the Gaussian model. No intrinsic asymmetry is found. This shows that the Gaussian profile is a proper assumption for the Gaia DIB and can be applied to the spectra from other spectroscopic surveys.