© ESO, 2013

1. Introduction

Interstellar extinction is a serious obstacle for the interpretation of stellar populations in the Galactic bulge. It shows an in-homogeneous clumpy distribution and is extremely high toward the Galactic center.

Several studies of cumulative extinction at the distance of the inner Galactic bulge have been made: Schultheis et al. (1999) used the DENIS near-infrared data set in combination with theoretical red giant branch (RGB)/asymptotic giant branch (AGB) isochrones from Bertelli et al. (1994). They showed that for an appropriate sampling area the observed sequence matched the isochrone with appropriate reddening well. The highest extinction that could be reliably derived from the J and Ks data available from DENIS was about A_V = 25^m. Schultheis et al. (1999) found that in some areas, presumably those with the highest extinctions, no J-band counterparts for sources were detected in Ks. For these regions, only lower limits to the extinctions could be obtained. Dutra et al. (2003) obtained similar results using the same technique with 2MASS data. Gonzalez et al. (2011) mapped the extinction along the minor axis of the bulge based on the Vista Variables in the Via Lactea (VVV) ESO public survey data (Saito et al. 2012) using the red clump giants. Gonzalez et al. (2012) presented the full extinction map over the entire VVV data set covering 315 sq. degrees using the red clump technique. This map extends beyond the inner Milky Way bulge also to higher Galactic latitudes and is the most complete study extending up to A_V ~ 35^m, which are much higher extinction values than were used in most previous studies. Recently, Nidever et al. (2012) determined high-resolution A_Ks maps using GLIMPSE-I, GLIMPSE-II, and GLIMPSE-3d data based on the Rayleigh-Jeans color excess (RJCE) method.

However, most of the studies are restricted to 2D extinction maps. Drimmel et al. (2003) built a theoretical large-scale 3D Galactic dust extinction map. Marshall et al. (2006) presented a 3D model of the extinction properties of the Galaxy by using the 2MASS data and the stellar population synthesis model of the Galaxy, the so-called Besançon model of the Galaxy (Robin et al. 2003). However, Marshall et al.’s study is limited by the confusion limit of 2MASS in the Galactic bulge region and by the limiting sensitivity of 2MASS in highly extincted regions (A_V > 30^m). Recently, Robin et al. (2012) improved the Besançon model significantly by adding a bar component which results in an excellent agreement for a comparsion of 2MASS color–magnitude diagrams in the metallicity as well as for the radial velocity distribution (see Babusiaux et al. 2010; Gonzalez et al. 2012; Uttenthaler et al. 2012). We used this updated model to build a new 3D extinction model using the GLIMPSE-II data, which map the Galactic plane within  ± 10 degrees in longitude from the Galactic center with a wavelength coverage between 3.6 to 8.0 μm, allowing one to map the most highly extincted regions, as shown in Schultheis et al. (2009). In addition, we used the recently published VVV data with a pixel size nearly ten times smaller and a sensitivity in Ks of 3–4 mag deeper than 2MASS to map the 3D extinction in the near infrared (IR).

The infrared color–color diagram J–Ks vs. K–[mid-IR] has been used by Jiang et al. (2006) and Indebetouw et al. (2005) to determine the extinction coefficients A_mid−ir along different lines of sight. Gao et al. (2009) traced the extinction coefficients based on the data from the GLIMPSE survey and found systematic variations of the extinction coefficients with Galactic longitude, which appears to correlate with the location of the spiral arms. Zasowski et al. (2009) found a Galactic radial dependence of the extinction law in the mid-IR using GLIMPSE data, while Nishiyama et al. (2009) determined the extinction law close to the Galactic center between 1.2 and 8.0 μm.

In this paper, we first introduce our data set in Sect. 2. We discuss in Sect. 3 the Besançon stellar population synthesis model together with the recent improvements. In Sect. 4 we describe the method for determining the extinction coefficients, and in Sect. 5 the method deriving the 3D-extinction. We finish in Sect. 6 with the determination of the extinction coefficient part and the conclusion (Sect. 7).

2. Data set

2.1. GLIMPSE-II survey

The Galactic Legacy Infrared Mid-Plane Survey Experiment (GLIMPSE) surveyed about 20 square degrees of the central region of the Galactic inner bulge using the Spitzer Space Telescope (Werner et al. 2004) equipped with the Infrared Array Camera (IRAC) (Fazio et al. 2004). It surveyed approximately 220 square degrees of the Galactic plane at four IRAC bands [3.6], [4.5], [5.8], and [8.0], centered at approximately 3.6, 4.5, 5.8, and 8.0 μm. The GLIMPSE-II data cover a range of longitudes from −10° to 10° and a range of latitudes from  ± 1 deg to  ± 2 deg depending on the longitude. The GLIMPSE-II coverage does not include the Galactic center region |l| ≤ 1° and |b| ≤ 0.75°, which was observed by the GALCEN program (Ramírez et al. 2008).

We used here the point-source archive (GLMIIA, Churchwell et al. 2009). The GLMIIA catalogs consist of point sources with a signal-to-noise ratio higher than 5 in at least one band and less stringent selection criteria than the catalog (see Churchwell et al. 2009 for more information). It contains the entire GLIMPSE-II survey region including the GALCEN data and GLIMPSE-I data at the boundary of the surveys (at l = 10° and l = −10°). Moreover, these catalogs were also merged with 2MASS, resulting in a seven-band catalog: three 2MASS bands and four IRAC bands. The photometric uncertainty of the GLIMPSE-II data is typically better than 0.2 mag and the astrometric accuracy is better than 0.3′′.

2.2. VVV-survey data

The other data set comes from the VVV ESO VISTA near-IR public survey. Observations were carried out using the VIRCAM camera on the VISTA 4.1 m telescope (Emerson & Sutherland 2010) located at the ESO Cerro Paranal Observatory in Chile. The five-year-variability campaign of the survey, which started observations during 2010, will repeatedly observe in an area of  ~564 deg²Ks-band , including a region of the Galactic plane from 295° < l < 350° to −2° < b < 2°, and the bulge spanning from −10° < l < 10° and −10° < b <  +5°. During the first year of observations, the complete survey area was fully imaged in five bands: Z, Y, H, J, and Ks. In this work, we used the J, H, and Ks catalogs for the same field of the GLIMPSE-II data from the 314 sq.deg coverage of the Galactic bulge. For a complete description of the survey and the observation strategy we refer to Minniti et al. (2010). A detailed description of the data products used in this work, corresponding to the first VVV data release, can be found in Saito et al. (2012).

A description of the single-band photometric catalogs, produced by the Cambridge Astronomical Survey Unit (CASU), and the procedure to build the multi-band catalogs can be found in Gonzalez et al. (2011). These catalogs were used in Gonzalez et al. (2012) to construct the complete 2D extinction map of the bulge and are the same ones used for the present study. We provide here only a brief summary of the procedure. J, H, and Ks catalogs produced by the CASU are first matched using the STILTS (Taylor 2006) code based on the proximity of sky positions. Since VVV photometry is saturated for stars brighter than Ks = 12 mag, 2MASS photometry is used to complete the bright end of the catalogs. Therefore, VVV catalogs are calibrated into the 2MASS photometric system. As described in Gonzalez et al. (2011), this calibration zero point is obtained by matching VVV sources with those of 2MASS in the magnitude range of 12^m < Ks < 13^m. This magnitude range was selected to avoid the saturation limit of VVV and the faint 2MASS magnitudes with poor photometric quality. This last point is also supported by using only sources flagged with high photometric quality in the 2MASS catalogs. The calibration zero point is calculated and applied individually for each VVV tile catalog. The final J, H, and Ks magnitudes in the calibrated catalogs are in the usual 2MASS system.

2.3. Completeness limit

Owing to the different pixel sizes of the different surveys (GLIMPSE-II, 2MASS, and VVV), completeness limits change as a function of wavelength and depend on the field direction. We divided the GLIMPSE II fields into subfields of 15′  ×  15′ each to ensure a sufficient number of stars (at least 200) for our fitting process (see Sect. 4). The completeness limit for every subfield was then calculated in each filter. To estimate the completeness we used star counts as a function of magnitude where the peak of each histogram in each filter gives the corresponding completeness limit. The catalog becomes severely incomplete beyond the point where the star counts are highest. We are aware that this is a simple treatment of the completeness limit. However, we also clearly see here how strongly the completeness limit changes along different lines of sight as demonstrated (Saito et al. 2012) using artificial sources.

Figure 1 shows the star counts as a function of magnitude for the GLIMPSE-II and 2MASS catalogs for two subfields, one located at (l, b) = (0, +1) which is our reference field (the so-called C32 field) because of its very homogeneous and low extinction (Omont et al. 1999, 2003). The second field is located at (l, b) = (0, −1) selected here for its higher and more clumpy extinction (Ojha et al. 2003). Due to crowding and extinction, the completeness limit can change within one or two magnitudes in the infrared bands. The difference in the completeness limit at [3.6] and [4.5] clearly visible, compared to [5.8] and [8.0] for both fields, indicating the higher sensitivity of [3.6] and [4.5].

Figure 2 shows the same fields, this time with star count histograms from the VVV data. Note the “artificial” rise at Ks = 12^m, which is the upper brightness limit of VVV. We decided to use the VVV data only from Ks > 12^m in the following analysis. Figure 2 demonstrates a clear gain in sensitivity of VVV in comparison with 2MASS near-IR photometry. As shown by Saito et al. (2012), the VVV data confusion are limited. However, because spatial scales of the two catalogs are not fully matched and because the overlap between 2MASS and VVV is in the magnitude range where the completeness between the two catalogs is very different, we did not investigate this problem in more depth (this is beyond the scope of this paper), but it should be kept in mind when assessing the results.

We divided the VVV data into the same subfields as in the GLIMPSE-II data set.

3. The Besançon model

3.1. Main description

The Besançon Galaxy model is a stellar population synthesis model of the Galaxy based on the scenario of formation and evolution of the Milky Way (Robin et al. 2003). It simulates the stellar contents of the Galaxy by using the four distinct stellar populations: the thin disk, the thick disk, the bulge, and the spheroid. It also takes into account the dark halo and a diffuse component of the interstellar medium. For each population, a star formation-rate history and initial mass function are assumed, which allows one to generate stellar catalogs for any given direction, and returns the magnitude, color, and distance for each simulated star as well as the kinematics and other stellar parameters. The model in its present form can simulate observations in many combinations of different photometric systems from the ultraviolet (UV) to the IR. The photometry is computed using stellar atmosphere models (Basel 3.1) from Lejeune et al. (1997) and (1998). However, while the Basel 3.1 stellar library is reliable for temperatures higher than 4000 K, atmospheric models become much more complex through the appearance of molecules such as TiO, VO, and H₂O for cooler stars. Because we have a significant contribution of M giants in our study, we decided to adopt a more realistic T_eff vs. color relation that extends to cooler objects of about 2500 K. This new relation is described in Sect. 3.2.

Recently, a new bulge model has been proposed by Robin et al. (2012) as the sum of two ellipsoids: a main component that is the standard boxy bulge/bar, which dominates the counts up to latitudes of about 5°, and a second ellipsoid that has a thicker structure seen mainly at higher latitudes. The resulting 2MASS color–magnitude diagrams as well as star counts of this two component model clearly indicate the significant improvement (see Robin et al. 2012). Uttenthaler et al. (2012) also showed an excellent agreement in the kinematics based on radial velocity measurements.

For more detailed information about the Besançon model we refer to Robin et al. (2003) and the most recent update in Robin et al. (2012).

3.2. Implementing IRAC colors and a new temperature–color relation

Reliable temperature–color relations are of extreme importance for comparing the synthetic color−magnitude and color–color diagrams with the observed ones. We used the reference field at (l, b) = (0, +1) to check the new synthetic colors in the [IRAC] filters. Figure 3 shows a kink in the Ks–[4.5] vs. [4.5] color−magnitude diagram at  [Ks]−[4.5]  ~  +0.1, [4.5]  ~ 10, which we clearly identify also in the corresponding temperature vs. color relation (see Fig. 4c). While the Basel 3.1 library predicts decreasing J–Ks values for T_eff < 3000 K, our new adopted relation shows a steady increase in J–Ks until T_eff = 2500   K. We also noticed that the full parameter space in the log T_eff vs. log g plane was not covered in the previous model, in particular, there is a gap in the temperature range between 2500–4000 K for 0.5 < log    g < 3.

We updated the temperature–color relation and filled in the gap, using the isochrones from Girardi et al. (2010). These isochrones include the TP-AGB phase (Marigo et al. 2008) including mass-loss. Starting from a group of isochrones with ages ranging between 5–10 Gyr, we derived a new temperature vs. color relation for metallicities between [Fe/H] = −2 and [Fe/H] = 0.0. These ages and metallicities are characteristic of the bulge populations (e.g. Zoccali et al. 2003, 2008; Bensby et al. 2011). The comparison between the color−temperature relations from the Besançon model (so-called old relation) and the relation derived by us (so-called new relation) is shown in Figs. 4b and c. For comparison we show the J–Ks vs. T_eff relation obtained by Houdashelt et al. (2000) for M giants. The clear difference between the new and old temperature vs. color relations for T_eff < 4000 K is obviously responsible for the bent in the color−magnitude diagram (Fig. 3), which disappears after adopting our new relation. We also see that using the new relation the J–Ks colors will be systematically redder by about 0.05 mag.

Figure 5 shows the comparison of a color–magnitude diagram in a low-extinction field, the so-called Baade window. We used here the catalog of Uttenthaler et al. (2010). One clearly notices that the predicted color in the old model is too blue compared to the observations and we see in addition a much higher percentage of foreground objects (Ks− [4.5]  < 0.5), which we do not find in the data.

4. Method

We obtained the extinction by using the color excess between the intrinsic color (derived from the model defined as Co_ins) and the observational apparent color (derived from the observation, here defined as Co_obs). Assuming an extinction law, we obtain: ${\mathit{A}}_{\mathit{\lambda }}\mathrm{=}{\mathit{c}}_{\mathit{\lambda }}\mathrm{\times }\mathrm{(}\overline{\mathit{C}{\mathit{o}}_{\mathrm{obs}}}\mathrm{-}\overline{\mathit{C}{\mathit{o}}_{\mathrm{ins}}}\mathrm{)}\mathit{,}$(1)where c_λ is related to the extinction coefficient. As we trace our fields close to the Galactic center, we used the extinction coefficients from Nishiyama et al. (2009). We ignored here the effect of the broadband extinctions and the reddenings on the spectral energy distribution (SED) of the individual sources because our method does not allow for a more sophisticated treatment, which could imply systematic errors. This results in $\begin{array}{ccc}& & {\mathit{A}}_{\mathit{K}\mathrm{s}\mathit{,J}\mathrm{-}\mathit{K}\mathrm{s}}\mathrm{=}\mathrm{0.528}\mathrm{\times }\mathrm{\left(}\mathrm{\right(}\mathit{J}\mathrm{-}\mathit{K}\mathrm{s}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{\left(}\mathit{J}\mathrm{-}\mathit{K}\mathrm{s}{\mathrm{)}}_{\mathrm{ins}}\mathrm{\right)}\\ & & {\mathit{A}}_{\mathit{K}\mathrm{s}\mathit{,H}\mathrm{-}\mathit{K}\mathrm{s}}\mathrm{=}\mathrm{1.61}\mathrm{\times }\mathrm{\left(}\mathrm{\right(}\mathit{H}\mathrm{-}\mathit{K}\mathrm{s}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{\left(}\mathit{H}\mathrm{-}\mathit{K}\mathrm{s}{\mathrm{)}}_{\mathrm{ins}}\mathrm{\right)}\\ & & {\mathit{A}}_{\mathrm{\left[}\mathrm{3.6}\mathrm{\right]}\mathit{,K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{3.6}\mathrm{\right]}}\mathrm{=}\mathrm{1.005}\mathrm{\times }\mathrm{\left(}\mathrm{\right(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{3.6}\mathrm{\right]}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{\left(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\right[}\mathrm{3.6}\mathrm{\left]}{\mathrm{)}}_{\mathrm{ins}}\mathrm{\right)}\\ & & {\mathit{A}}_{\mathrm{\left[}\mathrm{4.5}\mathrm{\right]}\mathit{,K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{4.5}\mathrm{\right]}}\mathrm{=}\mathrm{0.640}\mathrm{\times }\mathrm{\left(}\mathrm{\right(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{4.5}\mathrm{\right]}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{\left(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\right[}\mathrm{4.5}\mathrm{\left]}{\mathrm{)}}_{\mathrm{ins}}\mathrm{\right)}\\ & & {\mathit{A}}_{\mathrm{\left[}\mathrm{5.8}\mathrm{\right]}\mathit{,K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{5.8}\mathrm{\right]}}\mathrm{=}\mathrm{0.562}\mathrm{\times }\mathrm{\left(}\mathrm{\right(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{5.8}\mathrm{\right]}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{\left(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\right[}\mathrm{5.8}\mathrm{\left]}{\mathrm{)}}_{\mathrm{ins}}\mathrm{\right)}\\ & & {\mathit{A}}_{\mathrm{\left[}\mathrm{8.0}\mathrm{\right]}\mathit{,K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{8.0}\mathrm{\right]}}\mathrm{=}\mathrm{0.748}\mathrm{\times }\mathrm{\left(}\mathrm{\right(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\left[}\mathrm{8.0}\mathrm{\right]}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{\left(}\mathit{K}\mathrm{s}\mathrm{-}\mathrm{\right[}\mathrm{8.0}\mathrm{\left]}{\mathrm{)}}_{\mathrm{ins}}\mathrm{\right)}\mathit{.}\end{array}$ Below we apply the same method as Marshall et al. (2006) assuming that the distance grows as the apparent color/extinction increases. Figure 6b and c show an example of this relation for our reference field. In contrast to Marshall et al. (2006), we used the full information of the stellar populations in the color−magnitude diagram, that is, we also included the M-dwarf population. While the local dwarf population (d < 1–2 kpc) can be easily identified (see Marshall et al. 2006) by simple photometric color criteria (e.g. J–Ks), this is not the case for dwarfs located at larger distances where they cannot be separated from giants. This is especially true for the VVV data where we are able to detect dwarf stars up to larger distances. However, the model simulations show that the contribution of dwarf stars within the VVV completeness limit used in our study is rather small (~1%) but can increase up to 20–30% going to fainter magnitudes beyond the VVV completeness limit.

The method can be described in following steps:

• 1.

  First we add photometric errors and a diffuse extinction to theintrinsic colors Co_ins from the model for each subfield. The assumption of a diffuse extinction is necessary to ensure that the colors of the model increase with distance. We assume a diffuse extinction of 0.7 mag/kpc in the V band. We use exponential errors for the GLIMPSE-II data in all seven filters. Figure 6a shows an example of the photometric errors as a function of the [3.6] mag for our reference field (l, b) = (0.00, +1.00). The error is much smaller for the VVV data and we adopted a constant error of 0.05 mag (Saito et al. 2012). This procedure provides a first set of simulated colors (Co_sim,ini)

• 2.

  Colors are then sorted for observed and simulated data (Co_obs and Co_sim,ini), and normalized by the corresponding number of colors in each bin given by n_bin = floor(min([N_obs,N_bes] )/100). N_obs and N_bes are the total number of observed stars (in each subfield) and of stars in the Besançon model. This ensures that we have at least 100 stars per bin.

• 3.

  In each corresponding color bin, the Co_ins of the stars from the model data and the Co_obs from the observational data are used to calculate the extinction using Eq. (1), as well as the median distance from the model. A first distance and extinction estimate can be thus obtained (A_1st). Due to the saturation limit of VVV (Ks = 12), we do not get enough stars in the first distance bins (0–3 kpc). Therefore our 3D map starts at 4 kpc for VVV, and at 3 kpc for GLIMPSE-II data.

• 4.

  The extinction and the distance estimated from the third step is then directly applied to the intrinsic magnitudes in the model, which gives us a new simulated color (Co_sim,1 = Co_ins + A_1st).

• 5.

  We construct histograms for each color of the new simulated data (Co_sim,1) and the observational data (Co_ins) using the same binsize (0.05 mag). The χ² statistics (Press et al. 1992; Marshall et al. 2006) are then used to evaluate the similarity of the two histograms, given by $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{j}}(\sqrt{{\mathit{N}}_{\mathrm{obs}}\mathit{/}{\mathit{N}}_{\mathrm{sim}}}{\mathit{n}}_{{\mathrm{sim}}_{\mathit{j}}}\\ & & \mathrm{-}\sqrt{{\mathit{N}}_{\mathrm{sim}}\mathit{/}{\mathit{N}}_{\mathrm{obs}}}{\mathit{n}}_{{\mathrm{obs}}_{\mathit{j}}}{)}^{\mathrm{2}}>􏼮\mathrm{(}{\mathit{n}}_{{\mathrm{sim}}_{\mathit{j}}}\mathrm{+}{\mathit{n}}_{{\mathrm{obs}}_{\mathit{j}}}\mathrm{)}\mathit{,}\end{array}$(2)where N_obs and N_sim are the total number of observed and simulated stars in each subfield, while n_obsj and n_simj show the number of stars for the jth color bin of the observations and simulated data. With the first extinction estimate applied we now have a new set of simulated data. For this new set of data, some of the stars may be fainter and fall outside the completeness limit while others may be brighter and come inside the limit. We return to the second step with the new Co_sim,1 that would be used iteratively to refine our results. In total we perform 20 iterations for each subfield and each color. The minimum χ² is then our final result. We have chosen 20 loops to ensure that we always find a minimum value. In most of the cases the χ² converges after one or two loops to a minimum and becomes stable in the following loops.

Figure 6c shows an example of the iteration process using the Ks–[3.6] color of the GLIMPSE-II data in our reference field (l, b) = (0.00, 1.00). We clearly see the difference between the first extinction estimate and our final extinction determination after 20 loops. We adapted the same completeness limit as in the observations for the model simulations (see Sect. 2.3).

To test our method, we used artificial data by manually adding extinction values to the intrinsic magnitude of the model. We used different lines of sight with a different number of clouds. Figure 7 shows two examples of several clouds. For field (l, b) at (0.00, 1.00), one located at a distance of 2 kpc with A_V = 4.5 mag, the second one at 6 kpc with A_V = 2.5 mag. For field (l, b) = (5.00, 0.00), we assumed three clouds: A_V = 2.20 mag at 2.0 kpc, 2.50 mag at 5.0 kpc and 2.70 mag at 7.0 kpc. We applied our method to our simulated data. We found 328 000 sources for the first field and 182 000 for the second field. We artificially cut the data at Ks = 18 mag, which goes much deeper than the observational data. Figure 7 shows that our method successfully derives the extinction along the distance and that we are able to trace this multi-cloud structure. Modifying the input parameters of the Galactic model such as Galactic scale length, the size of the hole of the Galactic disk or changing the bulge luminosity function could lead to systematic differences in our derived extinction map. But this systematics, as shown in Marshall et al. (2006), is on the order of the random uncertainty of the method.

5. Results and discussion

The full results are listed in Tables 1 and 2, both available at the CDS. Each row of Tables 1 and 2 contains the information for one line of sight: Galactic coordinates along with the measured quantities for each distance bin E(J–Ks), E(H–Ks) and respective uncertainties (VVV, Table 1) as well as E(Ks–[3.6]), E(Ks–[4.5]), E(Ks–[5.8]), E(Ks–[8.0]) and their uncertainties (GLIMPSE-II, Table 2). These results will also be added to the BEAM calculator¹ webpage (Gonzalez et al. 2012) where the 2D VVV extinction map is already available to the community. Users of the BEAM calculator can choose to retrieve the extinction calculation adopting a specific reddening law and distance.

5.1. Extinction along different lines of sight

As an example of our results we discuss in more detail two specific lines of sight: the reference field (l, b) = (0.00, 1.00) and the field at (l, b) = (0.00, −1.00).

Figures 8 and 9 show the corresponding color–magnitude diagrams and color distribution for the observational data and the simulated data with the best-fitted final extinction in each filter. In addition to the general good agreement in the shape of the color–magnitude diagram and in the color distribution, we also emphasize that the number densities in the color–magnitude diagrams between the observations and the model agree very well. However, we notice some interesting features by comparing the color–magnitude diagrams of the model and the observations:

• We clearly see in the [5.8] vs. Ks–[5.8] and especially in the [8.0] vs. Ks–[8.0] GLIMPSE-II diagram the AGB phase. The tip of the AGB is approximately slightly brighter than 8.0 mag at  [8.0] . This phase is not reproduced by the model and needs to be improved.

• In contrast to GLIMPSE II, the model predicts too many stars for VVV in all filters at the fainter end of the color–magnitude diagram. We suspect that the completeness limit of VVV is underestimated. Saito et al. (2012) analyzed for two VVV fields the completeness limit by using artificial sources. They found that the source detection efficiency reaches 50% for stars with 16.4 < Ks < 16.9 mag for tile b314 located at (l, b) = (352.21, −0.92). Our derived completeness limit is Ks = 14.9, which reaches the 75% completeness level according to Saito et al. (2012). A detailed analysis of the completeness limits of the VVV survey using artificial sources would be needed.

• We notice that the model has produced a second branch in the J–Ks vs. Ks diagram at about Ks = 12.5 mag, which we do not find in the VVV data (see Fig. 9). This branch is related mainly to the Galactic disk K-giant population, which is clearly too prominent in the model. In contrast to the J–Ks and H–Ks color, the J–H color shows a large scatter in the temperature vs. color relation. This scatter is mainly due to the strong influence of J–H on metallicity and gravity (see also Bessell et al. 1989). We therefore decided to exclude the J–H color from our analysis.

Figure 10 to Fig. 13 show the change in interstellar extinction with distance in the different bands for four subfields. Different symbols stand for different data set: squares for VVV data and triangles for GLIMPSE-II data, while different colors denote the various filters. The errors of the derived distances and the extinction were calculated using the bootstrap method (Wall & Jenkins 2003). We re-sampled the original datasets of the observational and the model data 1000 times. The rms scatter gives the uncertainty of our results. Figure 14 shows a typical example. It shows the distributions of our results obtained for the 1000 samples obtained for one color bin (Step 3 in Sect. 4) for the GLIMPSE-II data in [3.6], left for the distance and right for the extinction. A Gaussian fit was made to the distributions, where the full-width-half maximum of the distribution gives the uncertainty in distance and extinction. This process is applied for each bin and each filter. The extinction changes along with the distance. As expected, it is higher for increasing distance. As seen in Fig. 10 to Fig. 13, the extinction becomes flatter beyond 7 kpc. As shown by Marshall et al. (2006), this is probably because inside the molecular ring the density of interstellar matter is lower than outside. It has also been shown that this region contains the dust lanes of the bar (Marshall et al. 2009), which can explain that at positive longitudes one can see an increase of the absorbing matter after this plateau, while it appears to be farther away at negative longitudes.

As already mentioned in Sect. 2.3, the VVV data used here are only for Ks fainter than 12 mag. For this reason, the first bin of the VVV data is 4 kpc. The VVV data on the other hand go much deeper (d ~ 10−12 kpc). But 2MASS starts to be incomplete at 12 mag depending on the line of sight. Therefore, there is in general a very small overlap between 2MASS and VVV, which is also seen in the color–magnitude diagrams (see Figs. 8 and 9) where we miss the bright foreground branch of K giants in the VVV data that is clearly visible in GLIMPSE II. For the four [IRAC] bands we see very similar extinction vs. distance relations.

5.2. Three-dimensional extinction map

To visualize the distribution of extinction in 3D, we divided the extinction between subsequent bins by the distance between them, similar as in Marshall et al. (2006). The distance intervals are interpolated in bins of 1 kpc. All 3D maps using the six different color excesses i.e. E(J − Ks), E(H − Ks), E(K −  [3.6] ), E(K −  [4.5] ), E(K −  [5.8] ), and E(K −  [8.0] ), are presented in Appendix A. The units are in kpc^-1. In addition we also show the error maps for each pair of filters. These errors are calculated using the bootstrap method. On average the errors are lowest between 4–8 kpc. We see a very similar picture of the dust distribution very concentrated to the Galactic plane in all maps. As expected, extinction becomes higher as the distance from the Sun increases.

We converted the seven color excesses to A_Ks using the extinction law derived in Sect. 6 separately. The extinction map of the weighted mean values is shown in Fig. 15. The maps are presented as extinction in Ks band and are also in units of kpc^-1. To better visualize the dust extinction along the line of sight, Fig. 16 shows a view of our results (upper panel) together with that from Marshall et al. (2006) (lower panel) from the north Galactic pole toward the Galactic plane at |b| < 0.25 similar to Fig. 9 of Marshall et al. (2006). Our results show generally similar structures compared to Marshall et al. (2006). Except for l = 0⁰ and l = −7⁰ we notice a very smooth and featureless dust extinction beyond 6 kpc. Clearly visible in this comparison is also the much lower extinction until 6 kpc compared to Marshall et al. (2006), which we discuss in Sect. 5.3. The highest concentration of dust is around the central molecular zone, a giant molecular cloud complex with an asymmetric distribution of molecular gas (see e.g. Morris & Serabyn 1996). The dust extinction map shows a non-axisymmetric structure very similar to that observed in OH and CO. However, the spatial resolution of our map is not high enough to trace features like the 100 pc elliptical and twisted ring of cold dust discovered by Molinari et al. (2011).

The highest extincted regions are best traced by E(H − Ks) using the VVV data and E(Ks −  [3.6] ) going up to A_Ks = 3.5 mag. Thanks to the high sensitivity of VVV, we are able to trace highly extincted regions in 3D (such as star-forming regions, molecular clouds, etc.) up to distances of 10 kpc.

Concerning the E(Ks −  [8.0] ) map, we emphasize that at 8 μm PAHs become an important component of emission that weakens the stellar flux (see e.g. Cotera et al. 2006). We find systematically higher extinction due to the PAHs emission.

Figure 17 shows the comparison in the two-dimensional space between the color excess using E(J − Ks) and the E(H − Ks), E(Ks −  [3.6] ), E(Ks −  [4.5] ), E(Ks −  [5.8] ), and E(Ks −  [8.0] ) to check whether we introduced any systematic offsets between the different filter pairs. Points with bootstrapping errors larger than 1σ were excluded. In general we see a very good relation with a difference not larger than 0.1 mag. The dispersion is higher for the GLIMPSE II data, which is due to the larger photometric errors.

5.3. Comparison with other 3D maps

We compared our results with those from Marshall et al. (2006) and Drimmel et al. (2003) for the extinction in the Ks band. Figure 18 shows the extinction v.s. distance diagram for two different fields, one located at (l, b) = (0.00, 1.00) and the other at extremely high extinction located at the Galactic center (l = 0.00, b = 0.00) compared to the dust model of Drimmel et al. (2003) and to Marshall et al. (2006).

We can see that they show the same trend as mentioned in Sect. 5.1. The black solid line is from Drimmel et al. (2003) and the dotted line from Marshall et al. (2006). To compare our results with those of Marshall et al. (2006), we converted their A_Ks values assuming the extinction coefficients of Nishiyama et al. (2009). Figure 18 shows that we obtain a similar A_Ks vs. distance relation for d > 6 kpc for our reference field (l, b) = (0, +1) while Marshall et al. predict a rather steep slope for d < 6 kpc. For the high-extinction field at (l, b) = (0, 0) there is quite a significant difference. We obtain a steeper slope and systematically higher extinction values for d > 4 kpc. We also added for comparison a data point of the 2D map of Gonzalez et al. (2012) that is marked as star symbol on the diagram. The distance was calculated by taking the difference of the dereddened mean magnitude of the red clump and the absolute magnitude of the red clump, which we assume to be M_Ks = −1.55 mag. Evidently the estimated distances from Gonzalez et al. (2012) agree remarkably well with our distances. Schödel et al. (2010) used adaptive optics observations of the central parsec of the Galactic center in the Ks band to map the extinction. They obtained a mean extinction value around Sagittarius A^∗ of A_Ks = 2.46^m with a distance to the Galactic center of R₀ = 8.02 kpc, while Fritz et al. (2011) derived A_Ks = 2.42^m. Our derived value of A_Ks = 2.34^m is within the errors of the cited papers.

For low-extinction fields, our results agree with Marshall et al. (2006). Only for larger distances do we obtain slightly higher A_Ks values than Marshall et al. (2006). However, for the highly extincted fields such as the Galactic center (right panel of Fig. 18), we notice quite significant differences: While the extinction seems to flatten at a distance of about 6 kpc for Marshall et al. (2006), we notice a steady increase in extinction, reaching values of A_Ks = 2.5^m compared to only A_Ks = 1.5^m for Marshall et al. (2006). This low limit is due to the large incompleteness of 2MASS in this highly extincted region compared to the VVV. We emphasize that compared to Marshall et al. (2006) our distance steps shown here are much smaller and we therefore trace the derived distance-color relations better.

Thanks to its deeper photometry and higher spatial resolution, the VVV survey catalog has stars at larger distances than 2MASS. For a red clump star with M_Ks = −1.65^m for example, a brightness at Ks = 12^m, which is the reliable measurement by 2MASS toward crowded regions, means a distance of about 3.5 kpc with A_Ks = 1 mag, and 5.5 kpc if no extinction is applied. Since the Marshall results are derived based on the 2MASS data, there are difficulties in probing far and highly extincted regions. In contrast, at Ks = 15 mag the red clump star in the VVV survey can be detected at a distance of 13.4 kpc if A_Ks = 1 mag and 5.4 kpc if A_Ks = 3 mag. Considering the RGB and AGB stars that are mostly brighter than red clump stars, the VVV survey detected many stars at a distance farther than 5 kpc. This depth of the VVV makes the results at larger distance more solid. This also explains the systematically higher extinction we obtained for the Galactic center field at l = 0.0 and b = 0.0 (see Fig. 18) compared to Marshall et al. (2006).

Our results are clearly very different from those of Drimmel et al. (2003), in particular for the Galactic center region. The Drimmel model predicts much higher extinction than ours or Marshall’s. The dust model of Drimmel et al. (2003) strongly depends on the dust temperature. However, in the Galactic center region the gas pressure and temperature is higher than in the Galactic disk (see Serabyn & Morris 1996), which explains the large difference in the extinction values compared to Drimmel et al. (2003). Our map is not sensitive to dust temperature like theirs, but to the modeled K/M giant stellar population. Because our map is restricted to a spatial resolution of 15′  ×  15′, small-scale variations such as seen, e.g., by Gosling et al. (2008) cannot be resolved.

5.4. Two-dimensional extinction maps

We compared our extinction maps with the 2D-maps of Schultheis et al. (1999) and Gonzalez et al. (2012). While Schultheis et al. (1999) mainly used the RGB/AGB star population, Gonzalez et al. (2012) restricted themselves to the red clump stars to trace the extinction. To compare this with our results, we integrated our extinction map up to a distance to 8 kpc. The resulting 2D-maps of the combined A_Ks are shown in Fig. 19. The individual maps of the different color excesses are presented in Appendix B. They clearly show dust features very similar to those described in Schultheis et al. (1999) and Gonzalez et al. (2012). We clearly see the small-scale variation of the interstellar dust clouds concentrated toward the Galactic plane, which is not symmetric. In addition to the clear dust features of the central molecular zone,we also see a clear correlation between the high extinction and the number of detected infrared dark clouds based on molecular line observations (Jackson et al. 2008), especially at the region around (l, b) = (353.0, 0.0) with a strong concentration of young stellar objects.

For a more quantitative comparison we smoothed the Gonzalez et al. (2012) maps to the same spatial resolution as ours (15′  ×  15′) and corrected for differences in the extinction law. Figure 20 shows the comparison. We see an excellent agreement between these maps with no systematic offset. For A_Ks > 1.5 mag the dispersion between our map and that of Schultheis et al. (1999) increases which is caused mainly by the limited sensitivity of the DENIS survey in the J- and Ks bands. In contrast, the dispersion compared to the VVV extinction is quite low, even for higher A_K. We are therefore confident that our extinction determination is highly reliable.

Figure 21 shows that the intrinsic dispersion grows with increasing extinction. Lada et al. (1994) demonstrated that the form of the observed σ_disp versus A_V relation can be used to place constraints on the nature of the spatial distribution of extinction. However, photometric uncertainties as well as uncertainties in the extinction law dominate this relation. The dispersion in the VVV is much lower than for 2MASS because of its much smaller photometric errors. A higher spatial resolution would be needed to study this feature in more detail. The 2D maps of Gonzalez et al. (2012) have a 2′  ×  2′ resolution resolution but their dispersion is dominated by the distribution of dust along the distance in that line of sight.

Figure 22 shows the total number of observed stars for each VVV subfield and the corresponding A_KS. We clearly see that the number of stars gives an indication of the extinction, i.e., lower numbers mean higher extinction.

6. Extinction coefficients

Because we determined the extinction maps using E(λ − Ks), we can calculate the extinction coefficients. For the extinction coefficients A_λ/A_Ks, we used here the values from Nishiyama et al. (2009) as initial values. If λ is the wavelength we are going to derive and α is the assumed wavelength, we can write ${\mathit{A}}_{\mathit{\lambda }}\mathit{/}{\mathit{A}}_{\mathit{K}\mathrm{s}}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{A}}_{\mathit{\alpha }}\mathit{/}{\mathit{A}}_{\mathit{K}\mathrm{s}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{\times }\frac{\mathit{E}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathit{K}\mathrm{s}\mathrm{)}}{\mathit{E}\mathrm{(}\mathit{\alpha }\mathrm{-}\mathit{K}\mathrm{s}\mathrm{)}}\mathit{,}$(3)where A_α/A_Ks is the initial value for the assumed wavelength from Nishiyama et al. (2009). E(λ − Ks)/(E(α − Ks)) is the color excess derived from the iteration of χ² test (see Sect. 4). Using Eq. (3), we derive the assumed wavelengths as J, H, [3.6],[4.5],[5.8],[8.0] to calculate the extinction coefficients for each subfield. The values we adopted are taken from Nishiyama et al. (2009) because the line of sight study is taken toward the Galactic center. However, considering the variation of the extinction law that is suggested by several studies, using a unique value of A_α/A_Ks could introduce some errors.

6.1. Mean extinction coefficients

Different choices of the assumed A_α/A_Ks would lead to different results for A_λ/A_Ks. We list our results in Table 3, varying the assumed extinction coefficient A_α/A_Ks. The mean extinction coefficients of all our fields are also listed in Table 3 and are compared to Gao et al. (2009), Indebetouw et al. (2005), and Nishiyama et al. (2009). The assumed values are noted in bold font in Table 3. Figure 23 shows the derived extinction curve. The black circles indicate our mean value (see Table 3). The red line indicates a simple power law A_λ ∝ λ^-1.75 from Draine (1989), the blue line the power law of Cardelli et al. (1989) with A_λ ∝ λ^-1.61, and the green line denotes the power law of Fitzpatrick & Massa (2009) with A_λ ∝ λ^-1.84. Clearly these simple power laws cannot explain the flattening of the extinction curve between 3–8 μm. As shown by Gao et al. (2009) (see their Fig. 7), a larger mean dust grain size (e.g. 0.3 μm) would give a steeper extinction law (i.e., lower A_λ/A_Ks ratios) than that of the diffuse ISM where the mean dust size is 0.1 μm. Our derived extinction curve is very similar to that of Nishiyama et al. (2009). Interestingly, we find a lower A _[5.8] /A_Ks value than Nishiyama et al. (2009). However, our value is within the errors close to the value of Nishiyama et al. (2009). Fritz et al. (2011) studied the near-IR and mid-IR extinction law toward Sgr   A^∗ using hydrogen emission lines. They found a power-law slope of alpha = −2.11 ± 0.06 shortward of 2.8 μm. This agree with our derived value of α = −2.175 and those found in the literature. Figure 23 shows that our derived extinction curve agrees remarkably well with that of Fritz et al. (2011) (indicated by the dot-dashed line) and confirms our low extinction value at [5.8]. As pointed out by Fritz et al. (2011), classical grain models (Li & Draine 2001; Weingartner & Draine 2001), using mainly silicate and graphite grains with different size distributions fails in reproducing the observed extinction law and especially the change of slope in the extinction law. Water ice features could be responsible for the slope change in the mid-IR. However, additional modeling is necessary.

The comparison shows that our mean values in the IRAC bands are comparable to those found in the literature. However, it depends on the assumed A_α/A_Ks. It can be seen from Table 3 that A_λ/A_Ks has lower values when A_α/A_Ks is bigger, which is a natural result from Eq. (3) because (A_α/A_Ks − 1) is negative for any wavelength longer than Ks. The values of A_J/A_Ks assumed or derived have a minimum of 2.56, a maximum of 2.89, and a mean of 2.70. In comparison, Gao et al. (2009) used A_J/A_Ks = 2.52.

6.2. Variation of the extinction coefficients

The extinction coefficient may differ along different lines of sight. Gao et al. (2009) showed an interesting distribution of extinction ratios as a function of Galactic longitude with 131 GLIMPSE fields. In this work, we assumed the extinction coefficients A_λ/A_Ks for all subfields in our field. We took the median value of the A_λ/A_Ks in each bin of longitude or latitude. To see the variation clearly, we present here the relative value ($\frac{\mathit{x}\mathrm{-}\mathit{x̅}}{\mathit{x̅}}$). Figure 24b shows a 2D map of A _[3.6] /A_Ks assuming A_J/A_Ks = 2.86, and Fig. 24a, c show the variation of the extinction coefficients (i.e., A _[3.6] /A_Ks,A _[4.5] /A_Ks,A _[5.8] /A_Ks, and A _[8.0] /A_Ks) as a function of Galactic longitude (Fig. 24a) and latitude (Fig. 24c). Overplotted in Fig. 24a are also the values from Gao et al. (2009). These authors used binsizes larger (1 degree) than ours. We notice a similar behavior in all [IRAC] bands in the variation of the extinction coefficient along the Galactic longitude, which follows the extinction curve variation of Gao et al. (2009). While this variation is quite small, we see a surprisingly strong peak visible toward the Galactic plane, indicating a strong variation with Galactic latitude (see Fig. 24b). This variation seems to be stronger for longer wavelengths except for 8.0 μm, where the amplitude becomes small again. This might be related to the 9.7 μm silicate absorption feature (Gao et al. 2009). Figure 24 indicates that the extinction curve becomes flatter in the mid-IR, approaching the Galactic center region. A detailed comparison with interstellar dust models is needed to explain this feature, but this is beyond the scope of this paper.

7. Conclusion

Using an improved version of the Besançon model, we presented here 3D extinction maps in the J, H, Ks, [3.6], [4.5], [5.8], and [8.0] bands using GLIMPSE II and VVV data. All extinction maps are available online at CDS (Tables 1 and 2), or through the BEAMER webpage². We derived new temperature−color relations for M-giants, which match the observed color–magnitude diagrams better.

We presented 3D maps from 1.25 μm until 8 μm. Owing to the high sensitivity of the VVV data, we were able to trace high extinction until 10 kpc. Our maps integrated along the line of sight up to 8 kpc show an excellent agreement with the 2D-maps from Schultheis et al. (2009) and Gonzalez et al. (2012). These maps show the same dust features and consistent A_Ks values.

Using the initial value from Nishiyama et al. (2009), we derived the mean extinction coefficient of our field in the seven bands A_J/A_Ks = 2.70 ± 0.13, A_H/A_Ks = 1.63 ± 0.03, A _[3.6] /A_Ks = 0.53 ± 0.02, A _[4.5] /A_Ks = 0.40 ± 0.02, A _[5.8] /A_Ks = 0.32 ± 0.04, and A _[8.0] /A_Ks = 0.42 ± 0.05. The variation of the coefficients indicates that the wavelength dependence of interstellar extinction in the mid-IR varies for different lines of sight, which means that there is no universal IR extinction law. This also agree with the study of Gao et al. (2009).

Online material

Appendix A: 3D reddening maps

Appendix B: 2D reddening maps

---

Acknowledgments

We want to thank the referee for his/her fruitful comments. We gratefully acknowledge use of data from the ESO Public Survey program ID 179.B-2002 taken with the VISTA telescope, data products from the Cambridge Astronomical Survey Unit and funding from the FONDAP center for Astrophysics 15010003, the BASAL CATA center for Astrophysics and Associated Technologies PFB-06, the MILENIO Milky Way Millennium Nucleus from the Ministry of Economy ICM grant P07-021-F, and the FONDECYT the Proyecto FONDECYT Regular No. 1090213. We are grateful to the GLIMPSE-II project for providing access to the data. B.Q.C. was supported by a scholarship of the China Scholarship Council (CSC).

References

All Tables

All Figures