DIFFUSE INTERSTELLAR BAND AT 8620 Å IN RAVE: A NEW METHOD FOR DETECTING THE DIFFUSE INTERSTELLAR BAND IN SPECTRA OF COOL STARS

For all the analysis presented in this paper we used spectra of the data release 4 (hereafter DR4), which includes 482,430 spectra of 439,503 objects. Stellar parameters are calculated for 410,837 objects and photometric distances for 397,783 of these objects. DR4 (Kordopatis et al. 2013) is based on the pipeline presented in Kordopatis et al. (2011) and is more than five times bigger than data release 3. It includes updated parameters compared to the older data releases, which are now more accurate, since the spectral degeneracies and the 2MASS photometry are better taken into account. These improved parameters lead naturally to improved accuracy on the stellar distances.

Figure 1 shows the range of several properties of the stars and their spectra. The observed fields cover declinations below zero, but avoid the Galactic plane at −5° < b < 5° and most of the Galactic bulge in the region with −25° < b < 25° and 45° > l > 315°. This means that there are almost no regions with high interstellar extinction covered and all observed DIBs will be weak. In order to detect a weak DIB in a given spectrum, it must have a high S/N. To detect the DIB 8620 in spectra of RAVE stars directly, the reddening of the sightline must be above 0.2 for the vast majority of stars and close to 1 for an average star. There are only a few observed regions with such a high reddening.

Zoom In Zoom Out Reset image size

The procedure to model the stellar spectrum described in this paper requires a large number of observed stars with no or very low reddening. Such stars were selected from the fields with b < −65°, which include 7.5% of all observed stars, 36,308 in total. As the stars are less abundant with rising temperature, the method becomes unreliable for stars with T_eff > 7000 K, because there are not enough stars at b < −65° to construct a stellar template at this temperature. Because the shape of the spectrum is less dependent on log (g), the method works at all gravities and becomes unreliable only at very extreme values of log (g) < 0.5. The same holds for the metallicity.

Because the DIB feature is expected to be very weak in most of the spectra, a precise elimination of the stellar component is the most important step in the analysis. Usually synthetic spectra are used for the stellar component, but they have many disadvantages. They are calculated on a grid of parameters, most often with a step of 250–1000 K for cool and hot stars, respectively. Other parameters are sampled at similar relative steps. Such a grid proves to be too widely spaced, and some weak stellar lines suffer from inappropriate oscillator strengths. Divided spectra show significant artifacts if synthetic templates are used, mostly close to the line inflection points or line centers. The Fe i line at 8621 Å, the strongest line that overlaps with the DIB, is no exception. There are also some spectral lines missing in synthetic spectra, but not in the band around the DIB, so this does not affect the measurements of the DIB profile directly. For normal stars the divided spectrum shows artifacts as large as 4% of the continuum. DIB 8620 is expected to be shallower than this.

Around half a million observed stars make a database large enough to include many spectra of different stars with very similar parameters, much closer to each other than the grid size in the database of the synthetic spectra. Many stars have high Galactic latitudes, where only a negligible amount of DIB component is expected. We assume that their spectra are free of interstellar or any other but stellar component. The expected equivalent width of the 8620 DIB can be estimated from the color excess from Schlegel et al. (1998). Schlegel et al. (1998) give color excesses for objects outside the Galaxy, so the given reddening is only the upper limit for reddening for stars in our Galaxy. For b < −65° the mean E(B − V) is 0.016. This corresponds to an equivalent width of 5.8 mÅ or a depth of around 0.08% of the continuum, if an FWHM of 5.59 Å is assumed (Jenniskens & Desert 1994) and if the average relation between reddening and DIB's equivalent width from this paper holds. This is the maximum error we expect when using spectra at high Galactic latitudes for the most inconvenient case. Because the Schlegel et al. (1998) values are known to overestimate the reddening (Yasuda et al. 2007; Kohyama et al. 2013) and because most stars have lower reddening than the calculated limit, the real errors are much lower than that. Several spectra are averaged to produce the stellar spectrum with higher S/N, so any stars with possible peculiar high reddening are averaged out.

The procedure described in the following paragraphs is illustrated in Figure 2. The spectrum used in this illustration is one of the few high-S/N spectra that show the DIB and is therefore an unusual case, but very convenient for the illustration.

Zoom In Zoom Out Reset image size

For each spectrum from which we want to extract the ISM spectrum, we find several spectra at b < −65° that most closely resemble the original spectrum. These spectra are called nearest neighbors, as they lie closest to the original spectrum in the parameter space.

The first step is to limit the number of spectra based on their stellar parameters. Only spectra with a derived effective temperature in the range of ±20%, log (g) in the range of ±0.6 dex, metallicity in the range of ±0.4 dex, and vsin (i) in the range of ±40% around the calculated values are retained. These limits are taken very liberally, as the errors of the calculated stellar parameters are always much smaller. In this way, the selection of spectra is limited to around 500 stars. The next step is to compare the target spectrum to these 500 stars and to find the most similar spectra. Spectra are shifted into the same velocity frame and interpolated to have the same sampling. The estimator for the similarity of the spectra is the sum of the positive differences of all pixels. However, several areas must be weighted differently. The area around the DIB (between 8612 Å and 8628 Å) is excluded from the difference calculations, and the weights on the central part of all three Ca ii lines are put to 30%. Ca ii lines might show chromospheric activity (Žerjal et al. 2013). Such stars are most common among cool dwarfs (Matijevič et al. 2012). The differences between the spectra are highly dependent on Ca ii lines, because they are the strongest lines in the observed spectral region. But also the chromospherically active stars can be used to model the stellar component of the spectrum. Therefore, the influence of the Ca lines is intentionally reduced. The inclusion of active stars with Ca emission shows itself in the final spectrum as artifacts at the wavelengths corresponding to Ca ii lines, but does not influence the band around the DIB. The whole observed spectral range is compared to find the nearest neighbors. Apart from the Ca ii lines, the spectra are dominated by Fe lines and a good match of the whole spectrum also means a good match of the DIB's neighborhood, at it is also dominated by Fe lines.

When all 500 spectra are compared with the original spectrum, up to 25 best matching spectra are averaged (with weights corresponding linearly to their S/N) into one stellar spectrum. For stars with more unusual parameters the number of spectra is smaller than 500. The smallest number of spectra that we ever used was 25. In this case, only the eight best matches are averaged into the final stellar spectrum. The target spectrum is then divided by this spectrum. We will call the divided spectrum the ISM spectrum. The generated stellar spectrum is calculated in the same way, regardless of the S/N of the target spectrum. Because we are aiming for the stellar spectrum to have high S/N, the target spectrum has almost always lower S/N than its generated stellar spectrum.

The mode of the S/N distribution is around 25. With a few exceptions, this is too low to detect the DIB in an individual star. Spectra of stars in an arbitrary volume can be averaged in order to achieve a better S/N. Before the averaging, the spectra are shifted back into the local standard of rest. If there are any spectra left that do not have the stellar component properly eliminated, this process smooths out such deviations. The measured DIB equivalent width is then proportional to the averaged gas density of the DIB's carrier toward each star in the region. For example, in a volume with average E(B − V) = 0.2, a spectrum with S/N around 300 is required for the DIB to have a depth equal to 3 standard deviations of the noise. With an S/N distribution with mode at 25, ∼100 spectra, occupying a certain volume, have to be averaged to achieve such a high S/N. In reality, the number is somewhat lower, because the weighted average is calculated with weights linearly depending on the S/N of the ISM spectrum.

To obtain the results presented in this paper, an average of 48 spectra was used. This is also the limiting factor in the spatial resolution due to the limited star density of the survey. The star density in the observed volume varies significantly, and so does the limiting spatial resolution.

The FWHM of DIB 8620 is around 5 Å, which is more than the resolution of the RAVE spectra (around 1.15 Å at the wavelength of 8620). This range is covered by around 12 pixels. Jenniskens & Desert (1994) actually report two features, described by two Gaussians with FWHMs of 1.86 and 5.59 Å at almost the same position, the narrower one being about three times weaker than the wider one. However, because the blend of both features is not far off from a Gaussian shape and because fitting two Gaussians to low-resolution RAVE data would be hard to automate and do it precisely, we used only one Gaussian to represent the DIB profile.

The equivalent width of the DIB is measured in two ways. The first one is by using the definition (integration)

$$\begin{equation} W_e=\int _{\lambda _1}^{\lambda _2}\left(\frac{F_c-F_\lambda }{F_c}\right)d\lambda , \end{equation} \tag{ 1 }$$

where F_c and F_λ are fluxes of the continuum and spectrum, respectively. The error is estimated from the noise in the regions near the DIB. The second method consists of fitting a Gaussian profile to the DIB using the Levenberg–Marquardt method. This method gives position, width, and equivalent width of the fitted profile, as well as the errors from the covariance matrix. The resolution and S/N are too low to fit more than just a simple Gaussian profile.

The limited S/N is not the only source of error on the equivalent width. A large error arises from the uncertainty of the continuum determination. Even though the used spectra were normalized, the ISM spectrum did not have a uniform continuum. The deviations from a constant continuum were of the order of 1%. Legendre polynomial continuums of different orders (up to 35th order) were fitted to the whole ISM spectrum (around 800 data points), producing a series of spectra with different continuum normalization. The above procedure of parameter determination was performed on each of the spectra in the series, and the scatter in each parameter was added to the measured errors. The average error of the W_e after this operation is around 10%, which is normal for such spectra and spectral lines. Not normalizing the ISM spectrum prior to the profile fitting would require an additional free parameter (the level of the continuum) to be introduced and would result in less stable calculation of the profile fit.

The fitting of the DIB profile is completely automated; hence, there may be the concern that in the low-S/N spectra with weak DIB the method might overestimate the equivalent width of the DIB. Several constraints are given for the DIB profile, one of them being that the DIB is an absorption feature. The DIB only covers around 12 pixels, so it is likely that in some cases the noise will resemble the DIB. Because the position of the DIB is also fitted, the procedure will search for such a dip in the area where DIB is expected to be. A simple integration over the same area is less sensitive to this effect, because coincidental rises in the continuum are equally probable and average out the dips. Figure 3 shows the extent of this problem. A typical overestimation at W_e = 0 is 3.5 mÅ and gradually drops to zero at larger equivalent widths. Data described in Section 4 were used to make this comparison.

Zoom In Zoom Out Reset image size

A nice test of the nearest neighbor method is to compare the stellar parameters and to classify the individual stars and their nearest neighbors.

The method of classification is described in Matijevič et al. (2012) and was used to classify most of the stars in DR4, used in this paper. Most of the stars (82%) are classified as normal stars, around 7% of the stars are unclassified for different reasons or have nearest neighbors with different classification (these stars will be named unclassified in this paper), 4% of stars are classified as stars with emission in the Ca ii triplet, 3% are spectra with problems with continuum fitting or wavelength calibration, and the remaining 4% are classified as binary stars, TiO band stars, hot (T > 7000 K) stars, cool (T < 3500 K) stars, peculiar stars, and carbon stars (in order from the most frequent to the least frequent classes).

Whether or not the individual stars share the same classification with their nearest neighbors is a powerful test of our method. For normal stars on average 98.7% of the neighbors are normal stars as well. The remaining neighbors have mostly unclassified spectra. The opposite holds for the unclassified spectra, where most of the neighbors are normal stars. Ca ii emission spectra have nearest neighbors that are in the same class only in 55% of the cases. Because most of these stars have weak emission and because different weights are set for the Ca ii triplet, the normal stars (representing the other 45% of the neighbors) are good matches to the spectrum with weak emission. Spectra with TiO band have nearest neighbors that belong in the same class in 17% of the cases. However, these stars are so rare that there are not enough stars of the same type in the DR4, and if we force the method to find a certain number of neighbors, it will usually include some of the normal stars. For other classes the comparison is difficult due to an even lower number of stars in each class.

This test gave positive results, because the stars have nearest neighbors of the same class much more often than pure coincidence would imply.

The second test was to compare the derived stellar parameters of individual stars and their nearest neighbors. Temperature, gravity, metallicity, and projected rotational velocities were compared. Figure 4 sums up the results.

Zoom In Zoom Out Reset image size

The width of the distribution is a combined effect of the final number of spectra from which the neighbors are selected and the accuracy of the stellar parameters that is achievable at a certain S/N. Temperature is the leading parameter that drives the shape of the spectrum. Therefore, it is calculated most precisely. A typical relative error is around 5%. In this case, the width of the distribution is mostly affected by the final number of spectra in the database. In order to construct a good stellar spectrum template, several spectra must be averaged. If there are not enough spectra with exactly the same temperature (within the error), spectra with similar temperature are taken. Gravity is calculated with lesser accuracy, and the width is mostly governed by noise. This is even more so in the case of metallicity, where the width of the distribution matches the accuracy of the calculated metallicity. We also compared the projected rotational velocity. Special care is needed in the interpretation of the results, as the rotational velocity at these noise levels is almost indistinguishable from the effect of the changing resolution of the spectrograph. Because we are only interested in finding the best-matching spectra, it is not important which of the two parameters the v_rsin (i) really traces.

On average all the distributions are centered at 0. However, there are non-negligible deviations in the distribution at a varying value of the parameter. Because the spectra have low S/N, a small perturbation in the spectrum can be attributed to a perturbation in more than one parameter. The correlation between the parameters is plotted in Figure 5. This does not mean that the nearest neighbors are not good matches; it just shows that there is a certain degree of uncertainty in spectra with mode of the S/N distribution around 25 (Kordopatis et al. 2011).

Zoom In Zoom Out Reset image size

We conclude that the nearest spectra have good matching parameters and classification to their target stars.

We want to check if the measured equivalent width of the DIB 8620 is independent of the temperature of the stars used to calculate the ISM spectrum. We took different regions on the sky where the DIB is expected to be strongest and calculated its equivalent width only with stellar spectra within a temperature range of 100 K. The regions were all located along 230° < l < 315° and in two bands between 5° < b < 10° and −10° < b < −5°. All regions covered 5° × 5° in l and b and a 0.5 kpc range in distance. We used regions with the nearest distance of 1 kpc, 1.5 kpc, and 2 kpc. Together this makes 102 different regions. Only stars with a temperature range of 100 K were used each time the equivalent width was calculated. Equivalent widths were then normalized to the average equivalent width in each region, so regions with different density of DIB carrier can be compared.

Figure 6 shows the result of the analysis. The average value of the normalized equivalent width in all regions is plotted against the temperature of the stars. The result is consistent with 1 for all temperatures, so the calculated equivalent widths are independent from the temperature of the stars used in the procedure. For a precise measurement of the DIB profile, we had to use regions as close to the Galactic plane as possible. Because the stars were divided even further into the temperature classes, only the regions with highest star density could be used, as the spectra of several stars must be combined to gain an S/N high enough to detect a DIB. This leads to a small number of regions considered in the calculation and to a wide distribution of the measurements along the W = 1 line.

Zoom In Zoom Out Reset image size

In order to have another correlation with which to compare the results of the new method, we repeated the measurements done in Munari et al. (2008). We detected the DIB in spectra of 144 stars with an effective temperature larger than 8000 K. We fitted a Gaussian profile to the DIB, together with a second-order polynomial, which represents a local continuum. Only the region between 8605 Å and 8660 Å, an area between the 14th Paschen hydrogen line and a blend of Ca ii and 13th Paschen hydrogen line, was taken into account when fitting the Gaussian profile and the local continuum. A narrow N i line at 8630 Å was also omitted.

The resolving power of the RAVE spectra is around 7500. This is high enough to resolve a broad DIB feature, but low enough to discard any substructure or peculiarities in the DIB profile when fitting a simple Gaussian curve to it. Spectra of hot stars are smooth enough, with no distinctive features around the DIB to make the elimination of the stellar component unnecessary. There is also no contamination by telluric lines. Measurements of the DIB were compared with the reddening from the literature (Neckel et al. 1980—18 stars; Savage et al. 1985—6 stars; Guarinos 1992—32 stars; Bailer-Jones 2011—82 stars; Gudennavar et al. 2012—3 stars). Some stars are present in more than one catalog. In addition, we used stars from Munari et al. (2008) and used reddenings from the same paper. However, we were only able to detect the DIB in 31 spectra out of 68. In other spectra from Munari et al. (2008), the measured DIB's equivalent width was consistent with zero on a one standard deviation level due to high noise. The measurements are shown in Figure 7. A linear relation was fitted to the measurements, taking the measurement errors into account.

Zoom In Zoom Out Reset image size

All the coefficients a for the linear relation E(B − V) = a · W_e between reddening and the DIB's equivalent width are listed in Table 1. Measurements based on a sample of hot stars have all been carried out with different detectors and/or with different methods to measure the equivalent width. Sanner et al. (1978) observed with a Reticon detector and integrated the spectrum to get the equivalent width, Munari (2000) observed with a thick CCD detector and measured the equivalent width by integration, and Munari et al. (2008) used RAVE observations, made with a thinned CCD detector, the same as in this paper, except that we fitted the continuum differently and also fitted the DIB, instead of integrating it. Results from Wallerstein et al. (2007) are inconclusive as discussed in Munari et al. (2008).

For correlations calculated with cool stars we were not able to use individual stars. High-S/N (S/N of more than 300 is preferred) cool stars that have high reddening are not abundant in the RAVE survey, and in addition the published reddenings in the literature matched only a few stars from DR4.

Instead, we combined several spectra in a small region to get one averaged high-S/N spectrum where we can measure the DIB profile with high precision, even if the reddening is small. For the densest regions we used 2° × 2° wide areas in l and b that were 0.1 kpc deep in distance. For less populated regions we used areas 5° × 5° wide and 0.2 kpc deep. We did this for −35° < b < 35°. The reddening of the region was calculated from the reddenings for all stars in the region, where we set the weights for averaging in the same way as the weights of the ISM spectra.

The reddening for individual stars is calculated as a by-product of the distance calculations. A Bayesian distance finding algorithm is presented in Burnett & Binney (2010). To include the effects of the interstellar dust, several modifications were made (Binney et al. 2013). The input parameters are J, H, and K magnitudes from the 2MASS catalog (Skrutskie et al. 2006), stellar parameters from RAVE, and color excesses from Schlegel et al. (1998), corrected according to Arce & Goodman (1999). In addition, several assumptions for the shape of the Galaxy (Jurić et al. 2008), metallicity and age distribution (Haywood 2001; Aumer & Binney 2009; Carollo et al. 2010), and isochrones (Bertelli et al. 2008) are made.

The correlation between the reddening E(B − V) and equivalent width of the DIB 8620 is represented in Figures 8 (correlation with Schlegel et al. 1998) and 9 (correlation with Bayesian reddening). The parameters of the linear relation are given in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Schlegel et al. (1998) do not give the reddening for individual objects. The product of Schlegel et al. (1998) is a map of extinction for the objects outside the Galaxy. For every object, or in our case every area on the sky, the given reddening is just the upper limit of the reddening toward the object. In addition, the reddening from Schlegel et al. (1998) is known to be overestimated (Arce & Goodman 1999; Yasuda et al. 2007), so the reddenings for the areas used in this study were corrected according to Yasuda et al. (2007). With just the upper limit, the linear relation between the equivalent width and the reddening cannot be acquired in the usual way, but the envelope must be fitted.

The difference between the Schlegel et al. (1998) and Bayesian reddening is evident from Figures 8 and 9. If we compared the two reddenings for every data point, the Schlegel et al. (1998) reddening would still represent the upper limit and would always be larger (within the error) than the Bayesian reddening.

A condition used to fit a linear relation to Schlegel et al. (1998) data was to find the flattest line, where 32% of the points lying below the line have 1σ errorbars that do not extend above the line. Error lines were defined as linear relations, where the number of points under the line changes by 32%. No offset in the linear relation was considered.

When comparing the Bayesian reddening with measured equivalent widths, a strong offset of (0.028 ± 0.002) mag is evident. The x-intercept is at −0.011 Å, which is significantly more than −0.0058 Å, which would be the x-intercept if the maximum error for the reddening of high Galactic latitude stars is assumed. No such offset is seen in Figure 8. There are also no points with reddening close to 0, which is unusual, as there are regions (like the ones in the Local Bubble) included where the reddening is expected to be extremely low (Jones et al. 2011). Whether or not there is a systematic overestimation of Bayesian reddening can be tested by plotting reddening as it depends on distance. We used 20 regions, 2° wide in b and stretching from l = 230° to l = 315°, and calculated the average reddening for stars at different distances. If we assume an exponential model for the dust density in the z-direction, the function

$$\begin{equation} E(B-V)=A[ 1-\exp {(-d/d_0)} ]+B \end{equation} \tag{ 2 }$$

can be fitted to any region. d is the distance, d₀ is the distance where the dust density drops e times (scale height divided by sin (b)), and A and B are free parameters depending on the Galactic latitude b and overestimation of the reddening, respectively. Figure 10 shows four examples of such an analysis. The average offset in all 20 regions was (0.026 ± 0.009) mag, which matches perfectly the offset of the linear fit in Figure 9. We conclude that the Bayesian reddening is overestimated by (0.026 ± 0.009) mag and that the systematic error for the equivalent width of the DIB arising from our assumption of unreddened stars at high Galactic latitudes is negligible.

Zoom In Zoom Out Reset image size