Free Access
Issue
A&A
Volume 538, February 2012
Article Number A106
Number of page(s) 14
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/201116512
Published online 13 February 2012

© ESO, 2012

1. Introduction

Since the discovery of a triaxial structure in the central regions of the Galaxy from 24 micron observations (Blitz & Spergel 1991), COBE NIR data (Weiland et al. 1994; Dwek et al. 1995), or visible photometry (Stanek et al. 1994), numerous attempts have been made to characterize this structure and to investigate its origin. It is still unclear whether this structure, often called the outer bulge but sometimes the bar, had its origin in the early formation of the spheroid (like a typical bulge, similar to ellipsoidal galaxies) or was formed by a bar instability later in the disc. It is crucial to investigate this formation history as a benchmark in understanding formation of disc galaxies. Kinematical studies are very important for tracing the stellar motions and deducing a dynamical history of the populations in the inner regions. New surveys in the near future will help solve this question, particularly the Gaia mission of the European Space Agency.

In the meantime photometric surveys, especially in the infrared, can help to determine the shape and density laws of the stellar populations in place. The question of the shape, orientation and extent is still open, because one can find for the angle of the major axis direction with regard to the Sun-Galactic centre direction values between 10° and 45°, depending on the fields considered, the wavelength and the method used. Among all the attempts, the significant ones include Stanek et al. (1994), Dwek et al. (1995), Fux (1997), López-Corredoira et al. (2000), Freudenreich (1998), Picaud & Robin (2004), Babusiaux & Gilmore (2005), Benjamin et al. (2005), and Rattenbury et al. (2007). The difficulty in accurately measuring the orientation of the bulge might be aggravated if, actually, there is more than one population present in the region and if these populations have different origins.Martinez-Valpuesta & Gerhard (2011) also offer a clever explanation for the lack of consensus concerning the bar angle, which they think is partlyfrom the cone effect and partly from the presence of leading or trailing overdensities at the end of the bar, as also shown by Romero-Gomez et al. (2011).

In this context, metallicity distribution studies can be useful constraints for distinguishing populations. Zoccali et al. (2008) and Hill et al. (2011) highlight a bimodal metallicity distribution in Baade’s window and several fields on the bulge minor axis, thus identifying a mixture of populations. Concerning the kinematics, proper motions are difficult to obtain due to the large distance of the bulge, while radial velocities are expensive to process for large numbers of stars. Sumi et al. (2004) studied proper motions from the OGLEIII survey and compared them with a dynamical model. However, the OGLE fields avoid the Galactic plane and the northern part of the bulge, they are limited by extinction (because the observations were performed in the visible), and the mean proper motion error is too large to accurately determine velocities at the bulge distance.

Rich et al. (2007) have undertaken a wide survey (BRAVA, for bulge radial velocity assay) of radial velocities at longitudes between –10° and 10° and at latitudes –4° and –8°. Their first results from the latitude –4° data set show that the slope of the rotation curve flattens considerably at longitude |l| > 5° indicating a probable change in the dominant population at this point. Furthermore, the analysis of metallicity versus kinematics by Babusiaux et al (2010) for RGB stars in three bulge fields, at different latitudes along the minor axis, show that two populations with distinct chemo-dynamical signatures seem to be present in these fields. Since 2010, there have been several studies showing double clumps in a number of fields in the bulge (Nataf et al. 2010; McWilliam & Zoccali 2010; Saito et al. 2011). This is also a feature that should be taken into account in an overall view of the bulge region and that the models should explain.

Complete understanding of the bulge’s structure and formation would require studying the abundances and kinematics in various fields at different longitudes and latitudes,such as the APOGEE project (Majewski et al. 2007). In the mean time we propose to explore the bulge structure based on an analysis of existing photometric data. We intend to check whether a detailed modelling of the stellar populations and interstellar extinction present in the central region might lead to an answer. We use the Two Micron All Sky Survey, hereafter 2MASS (Skrutskie et al. 2006) photometry and a model-fitting procedure based on the Besançon Galaxy model (BGM) (Robin et al. 2003) and explore a wide region –20° < l < 20° and –10° < b < 10° to fit a multi-dimensional space of parameters defining the density laws of the bulge and disc populations. The analysis proceeds in several steps. First we attempt to fit the whole area with a Monte Carlo scheme to explore the parameter space using the differential star counts in bins of K and J − K in a set of selected windows. We compare the resulting model with integrated star counts in order to visually check the realism of the resulting model. In a second step we compare simulated CMDs in a sample of directions to check their agreement in more detail. At this stage we introduce a modulation of the density distribution of the bar structure to reproduce the double clump structure seen in medium latitude directions by Nataf et al. (2010), McWilliam & Zoccali (2010), and Saito et al. (2011). In a final step we compare model simulations of the metallicity distribution functions with data over the minor axis of the bar, placing constraints on the mean metallicity of the composite populations in the bulge region.

In Sect. 2 we describe the basis of the BGM and the free parameters which we adjust to the existing data. The data are described in Sect. 3 along with the method that corrects for the 3D extinction distribution towards the bulge. The fitting method and results are given in Sect. 4. Validation of our extinction method is provided in Sect. 5. In Sect. 6, we compare model predictions with CMDs in several directions and explain how the model can reproduce the double clump feature in the CMDs. In Sect. 7, we zoom in the internal bulge. In Sect. 8 we discuss the results, and the metallicity distribution function, and present the perspectives of this work.

2. Population synthesis model for the bulge

The stellar population model is based on a population synthesis scheme. Four distinct populations are assumed (a thin disc, a thick disc, a bulge, and a spheroid), each deserving specific treatment. In the central regions, only the bulge and thin disc are important. The thick disc is a minor component, and the halo completely negligible at the magnitudes and latitudes considered here. The thick disc contributes slightly to the star density but only when the bulge and the thin disc become less important, that is, at latitudes above about 8–10°.

2.1. The thin disc population

The thin disc contributes significantly to star counts in the Galactic central region, so its characteristics have to be fitted at the same time as the bulge parameters. The important parameters are the scale length and the size of the central hole. Other parameters and structures, such as outer radius, warp, and flare, are taken into account in the BGM and described in Reylé et al. (2009), but are not relevant here.

A standard evolution model is used to produce the disc population, based on a set of evolutionary tracks, a constant star formation rate (hereafter SFR) over 10 Gyr, and a two-slope initial mass function φ(M) = A × M − α with α = 1.6 for M < 1 M and α = 3.0 for M > 1 M. The preliminary tuning of the disc evolution parameters against relevant observational data was described in Haywood et al. (1997a,b) and later changes are explained in Robin et al. (2003).

With the evolution model we populate the thin disc dividing it into seven age components. For each subcomponent, the distribution in absolute magnitude and effective temperature is obtained, assuming a star formation history constant between 0 and 10 Gyr. The star counts are computed using the standard equation of stellar statistics from the distribution in MV and a spatial density distribution. It is therefor equivalent to assuming that the SFR is roughly the same over all the thin disc.

The thin disc density distribution model follows the Einasto (1979) law: the distribution of each disc component (except for the very young one of age less than 150 million years, which is not important here) is described by an axisymmetric ellipsoid with an axis ratio depending on the age. The density law of the ellipsoid is described by the subtraction of two functions:

with , where:

  • R and Z are the cylindrical Galactocentric coordinates;

  • ϵ is the axis ratio of the ellipsoid. Values of ϵ as a function of age and local normalization are given in Robin et al. (2003);

  • Rd is the scale length of the disc and is around 2.2–2.5 kpc (Ruphy et al. 1996). This parameter will be fitted in the procedure described in this paper;

  • Rh is the scale length of the hole. It is also fitted here. The maximum density of the disc population is approximately at 2 kpc from the centre;

  • The normalization ρd0 is deduced from the local luminosity function (Jahreiß et al., priv. comm.), assuming that the Sun is located at R = 8 kpc and Z = 15 pc. Alternative values for the sun-Galactocentre distance R are also considered.

2.2. The outer bulge

Picaud & Robin (2004) undertook a detailed analysis of the outer bulge stellar density and luminosity function by fitting model parameters to a set of 94 windows in the outer bulge situated at –8° < l < 10° and –4° < b < 4°. The data were obtained by the DENIS survey team (Epchtein et al. 1997) using Ks magnitude and J − Ks colour distributions. Using a Monte Carlo method to explore an 11-dimensional space of bulge and disc density model parameters and a maximum likelihood test, they showed that the bulge follows a boxy exponential or a boxy sech2 profile like the one fitted to the DIRBE integrated flux in the near infrared (Freudenreich 1998). These are Ferrers ellipsoids (Ferrers 1877). Their resulting triaxial bulge has a major axis pointing towards the first quadrant with an angle of about 10° with respect to the Sun-Galactic centre direction. A full description of the parameter values of the bulge density law can be found in Picaud & Robin (2004).

Here we reconsider the density law fitted by Picaud & Robin (2004). The Freudenreich S function is:

with

(X,Y, Z) are the Cartesian coordinates in the referential of the triaxial structure (X being the major axis, Y the second axis and the Z the third). The parameters c ∥  and c ⊥  are important for exploring a wide range of shapes, from “disky” to “boxy”. This allows great flexibility: one can even have a “disky” shape in the plane, together with a boxy projection vertically.

In the following we also consider alternative functions for the overall shape, exponential (E), and Gaussian (G):

and

The density function is then multiplied by the cut-off function fc (distances given in kpc, and Rc is called the cut-off radius): Three angles define the orientation:

  • φ: orientation angle from the sun-centre direction of the projection on the Galactic plane of the bulge major axis;

  • β: pitch angle of the bulge major axis from the Galactic plane;

  • γ: roll angle around the bulge major axis.

In the Picaud & Robin (2004) study, the angle β was found very close to 0°, and the third angle γ was found to be ill defined because the minor axes Y and Z had similar scale lengths. In this study, we therefore adopt β = γ = 0°. Finally, the angle φ, the three scale lengths x0, y0, z0, the density at the centre ρ0, the cut-off radius Rc, and the two coefficients C and C are considered as free density parameters in the fitting process.

The bulge stellar population is taken from a single burst population of ages varying between 6 and 10 Gyr as tested by Picaud & Robin (2004). The favoured combinations between the age and evolution models are from Bruzual et al. (1997) with an age of 10 Gyr, or from Girardi et al. (2002) with an age of 7.94 Gyr (log(age) = 9.9). It appears that the Ks band counts are not very sensitive to the age, so a proper estimation of the age and age range of the populations in the bulge would require complementary data. This analysis thus takes the luminosity function from Girardi et al. (2002) with a log(age) of 9.9 as a reference.

2.3. One or two populations in the bulge?

In the bulge region there may be several populations, with several star formation histories, cohabiting. The coexistence of different populations in the bulge region can be tested using the distribution in colour-magnitude diagrams. In our fitting process we consider alternatively one or two bulge populations, and for each of them the S, E, or G functions.

3. Data description and extinction corrections

3.1. 2MASS data

The Picaud & Robin (2004) analysis was limited in longitude and latitude (–12°  < l < 8°, –4°  < b <  4°), used DENIS data, and did not contain many fields near the plane (less than 10% were situated at |l| < 1°). Here, we use the 2MASS point source catalogue data. We concentrate on the region –20°  < l <  20°, –10°  < b <  10°. Unlike the Picaud & Robin (2004) study, the selection goes down to latitude b = 0°. The data are used in three ways: (i) for determining the extinction using the method from Marshall et al. (2006); (ii) for fitting bulge model parameters on a selection of windows; (iii) for comparing the whole bulge region with the fitted model and checking the extinction distribution again.

3.2. Extinction

Since the extinction has a major effect on star counts close to the plane, a detailed analysis of the extinction field by field is needed to provide realistic model simulations. For this purpose, we used the method described in Marshall et al. (2006). This method adopts the standard BGM and determines the extinction by comparing the simulated colour-magnitude diagrams from this model and the observational data. Adjusting both the extinction and the bulge parameters in a single run is a difficult data analysis problem. Therefore we proceeded with iterations, first fitting the extinction with the preliminary bulge model (from Picaud & Robin 2004), then adjusting the bulge parameters assuming this extinction and finally verifying that the extinction derived from the new fitted bulge model is not significantly different from the original extinction distribution.

The goodness of fit used in the extinction determination is based on a chi-squared statistic for two binned distributions having different total number of elements (Press et al. 1992). As such, the extinction determination is sensitive to differences in colour and not to discrepancies between the number of observed and modelled stars.

In the present study, the Ks magnitude and the J − Ks colour were used to compare with simulations. Cuts were made in Ks and J at the completeness limit, defined from the histogram in magnitude at the bin before the maximum for limiting the risk of incompleteness that can occur even at the peak in magnitude. In the extinction determination and in the maps for visualizing the comparison between model and data (see Sect. 4.1 and Figs. 2–4) we further restrict the comparison at Ks < 12 because the extinction determination method is only based on the giants’ colour and that they can be separated from the dwarfs only up to this magnitude.

4. Fitting procedure

The first set of data we use is a selection of fields from the 2MASS point source catalogue. Two hundred windows have been chosen in order to cover the region defined by –20°  < l <  20° and –10°  < b <  10°. For each window, the extinction is determined using the method described in Sect. 3.2, as well as the completeness limit independently in each band. The BGM in its standard version is used to simulate the NIR data, taking the completeness into account. To illustrate, we show in Fig. 1 the 200 windows and the completeness limit in the Ks band. Then the following procedure is used to adjust the bulge model parameters.

thumbnail Fig. 1

Distribution in longitude and latitude of the 200 windows used in the fitting process. The colour coding corresponds to the limiting magnitude of the counts in Ks.

For each window the star counts in the CMD diagram are binned in Ks and J − Ks in bins with equal numbers of stars, two bins in colour, and seven in magnitude, giving a total of 14 bins. In the model the counts of each population are counted separately to be able to change them as a function of the variable parameters. The parameter space is explored randomly from successive Monte Carlo drawings and a selection of the best models in the sense for the likelihood. The likelihood function is the same as in Picaud & Robin (2004). The process is stopped when the variation in the parameters from one iteration to the next is less than the required accuracy for each parameter. The process is done 20 to 50 times to ensure a robust solution.

In a second step, we compare the fitted model with 2MASS data in the whole region  −20°  < l < 20°and  −10°  < b < 10° and compute the relative difference between data and model on a grid of 15 × 15 arcmin fields. The aim is twofold. Firstly, it ensures that the 2MASS data are correctly fitted by the model all over the bulge area (not only in the 200 selected windows), without systematics in any region. Secondly, it verifies that the new extinction, recomputed from the new density model, has not varied significantly since the density fitting process, so we ensure that proceeding in two separate steps (fitting of the extinction, then fitting of the density model) is efficient and will converge.

We first attempted to fit a unique bulge population all over the considered region. The parameters that are fitted are the following: the bulge angle φ (angle between the major axis of the bulge and the Sun-Galactic centre direction), the bulge scale lengths (considering 3 axes, x0, y0, z0), the bulge mass, the extent in radial distance (cutoff Rc), two parameters describing the boxiness of the bulge c ∥  and c ⊥  (see above), the disc scale length, and the scale length describing the size of the hole inside the disc.

As we show below, the agreement was not good enough, and we had to consider alternative models. We attempted to fit a second structure, with the same kinds of parameters as for the first bulge structure. We also considered different formulae for the bulge density and for the second structure: sech2, Gaussian, or exponential.

We assume here a Sun to Galactocentric distance of 8 kpc (following recent results from Ghez et al. 2008), unlike Picaud & Robin (2004) who assume a value of 8.5 kpc. We investigate the effect of this value on the result in the discussion.

In the next subsections we describe the results for each kind of model fitted.

4.1. Fitting a single ellipsoid

The fit of a single ellipsoid on the large area considered here leads to significantly different results compared with the analysis done by Picaud & Robin (2004) where the data set considered was restricted to the region  −12°  < l < 8° and  −4°  < b < 4° with about 100 windows. The important thing here is probably not the number of windows, which is doubled, but rather the extent of the area.Moreover, unlike the Picaud & Robin (2004) study, we drop the selection of giants and consider a fit covering the whole colour range. Finally, we account for interstellar extinction in a more realistic manner than they did.

Table 1

Best fit models with 1 single bulge ellipsoid.

Results are given in Table 1. Model A indicates the parameters fitted on DENIS data by Picaud & Robin (2004) for comparison only. In the model A fitting process, the parameters c and c ⊥  were fixed to 2 and the Sun-galactocentric radius to 8.5 kpc. For each model and each parameter, the dispersion around the fitted values is indicated in the second line. When a model has a parameter fixed rather than fitted, it is in italics. The B model is the model with a free angle with the (default) sech2 shape. C is similar but with a fixed angle of 20°, D is similar but with a fixed angle of 25°. The G model is similar to the B model but with a Gaussian shape, the E model has an exponential shape (see below).

thumbnail Fig. 2

Star counts up to limiting magnitude Ks from 2MASS data (top left) compared with fitted 1 ellipsoid model B (top right), residuals (Nmod − Nobs)/Nobs (bottom left), and difference divided by the Poissonian counting error of the model (i.e. (Nobs − Nmod)/ (bottom right). In the residual map, contours are drawn at intervals of 20% model overestimate (black) and 20% model underestimate (white). The 1-ellipsoid model leaves significant X-shaped residuals. Near the Galactic centre the nuclear bar population is missing in the model. The residuals in the outer regions are not very significant due to the small number of stars in each bin (seen in the top left panel in dark violet and black).

The most significant difference between the new result (model B) and the old one (model A) is the new scale length of the bulge (major axis scale length x0), which was found to be 1.59 kpc in Picaud & Robin (2004) and is here about 4 kpc. The cutoff radius has also grown from 2.54 to about 6 kpc. It means that the limited region used in our first study has led to underestimating of the bulge’s extent. The second axis scale length is also greater than in Picaud & Robin (2004) (corresponding to the depth of the bulge along the line of sight if the angle is small), while the scale height (minor axis z0) has not changed (the latitude range covered by the previous paper was sufficient to cover the bulge). We also found that the scale length of the disc has slightly diminished and that the hole has diminished.

We plot in Fig. 2 star counts in Ks from 2MASS data (top left) compared with the fitted model with single population (top right), residuals as defined by (Nmod − Nobs) / Nobs (bottom left), and difference divided by the Poissonian counting error of the model (i.e. (Nobs − Nmod)/ (bottom right). When computing the extinction, the Ks band was limited to 12 or brighter. The other bands have no such limitation. This limit ensures lower uncertainty on the 2MASS PSC fluxes. This figure shows the B model described in Table 1 after the 3D extinction distribution has been refitted. The structure of the residuals shows that the shape of the fitted bulge model is not appropriate. In the in-plane region the model density agrees well with the data. But the fields at high latitudes and intermediate longitudes have X-shape residuals, where the model density is too low. It is probable that the resulting fit is a bad compromise between adjusting the densities in plane and at different latitudes. It is worth noting that in the top left-hand panel the 2MASS star count map indeed shows boxy iso-densities that are not reproduced by the model, even if the density laws were chosen explicitly to allow for a boxy shape.

Interestingly, the residual map shows a structure in the inner bulge at longitudes less than 2°, which we believe is the nuclear bar seen by Alard (2001). This population has not yet been included in the model, and is shown here by subtraction. We discuss the point in Sect. 7.

4.2. Influence of other parameters

The Sun-Galactocentric distance R is not adjusted here during the fitting process. It is fixed at either 7., 7.5, 8, or 8.5 kpc. Models with R of 8 kpc give the larger likelihood, but the data are not sensitive enough to distances along the line of sight to strongly constrain this parameter. Moreover, changing R has little impact on the other parameters, such as c and c, which stay of the order of 1 and 4 (resp).

We assume an age of about 8 Gyr for the bulge population (log(age) = 9.9), following the best result in the Picaud & Robin (2004) study. Our analysis is not very sensitive to the assumed age for an old population (age over 5 Gyr) in the red clump. The difference in the luminosity function would be noticeable if we had covered the red clump and the giant branch down to the turnoff of the population, which is not possible with 2MASS data. The turnoff position and the giant sequence are also influenced by the age. The turnoff would be slightly brighter in a younger population (but this is not visible in most 2MASS fields), and the giant sequence would be slightly bluer. But this effect is of second order and difficult to see in broad band photometry and even more in extincted fields. Complementary data would be necessary to constrain the age of the population at the level of a Gyr. Consequently, we consider an isochrone for the bulge with an age of about 8 Gyr from Girardi et al. (2002), which should be valid if the true age of the population is in the range 6–10 Gyr. Alternative isochrones within a reasonable range of age would not significantly change the conclusion on the bulge shape and the extinction.

4.3. About the bulge/bar angle

Many bulge studies find a significantly larger bulge angle than our result (from one population model B: 7.1°). Most of them are in the range 10° to 45°. We found 11.1° in model A. Since some of our parameters are slightly correlated, we attempted to fix the bulge angle at different values to see whether it is possible to find a good enough fit with a larger bulge angle and to see in which field this value is favoured. Models C and D in Table 1 have fixed bulges angle of 20° and 25°  respectively. These models have very low likelihood with regards to smaller angles.

We suspect that the bulge angle found in different studies might depend on the field selection, extent in longitude, distance to the Galactic plane, and depth of the survey. Martinez-Valpuesta & Gerhard (2011) have shown how the bar angle determination can be strongly biased by the cone volume effect, in the case of the long bar, but this effect is also present when determining the main bulge angle using the mean distance of the red clump.

To test these effects, we selected subsets of fields in longitude and latitude and reproduced the fitting procedure. It appears that for fields at |l| > 12° and close to the plane (|b| < 3°), the likelihood is higher with a large angle, while this large angle gives a bad fit at low longitude and high latitude (|b| > 3° and |l| < 5°). Fields at high longitudes and low latitudes are well fitted by large bulge angles, while fields at low longitudes and high latitudes favour small bulge angles. When fixing the angle of the bar at higher values (models C and D, 20 or 25°), the fields situated at large longitudes show similar likelihood, but fields at higher latitudes (|b| > 3°) get significantly worse. This is in good agreement with the analysis from Cabrera-Lavers et al. (2007, 2008), from 2MASS data as well, but limited to fewer fields. These findings might indicate that a single ellipsoid fit is not adequate. This is why in the following section we attempt to fit the sum of two ellipsoids on our data set.

4.4. Fitting two ellipsoids

We attempt to apply the procedure to fit two ellipsoids (in addition to the thin disc) having the same luminosity function, age, and metallicity, as these parameters do not influence the overall near infrared photometry very much (see the discussion below). The same parameter set describing each ellipsoid is fitted on both populations. We consider several combinations of shapes, S shape as before (S), exponential (E), Gaussian (G) as in Dwek et al. (1995).

Table 2 shows the parameters of the best fit with different shapes S, E, or G. Because it is unlikely that two populations have different angles from a dynamical point of view, and because the fit with different angles did not give significantly different angles, we have imposed the two ellipsoids to have the same orientation in the plane.

Table 2

Parameters and reduced likelihood Lr when two ellipsoids are fitted, including main sequence and giants, and forcing the two ellipsoids at the same orientation.

Among the best combinations obtained for the best model shapes, the parameters are very similar: the best fit is obtained with the two ellipsoids having an orientation of about 9° to 13°, the first ellipsoid with a scale length x0 of 1.6 and secondary scale length of 0.5 and 0.4 kpc. Values are slightly different when an exponential shape is assumed, because the exponential is significantly more like a disc than a sech2 or a Gaussian, and the different shape is compensated for by different scale lengths and normalization. Regarding the likelihood, the best solution is with a sech2 for the first component and an exponential for the second one.

The second ellipsoid is generally less massive than the first one but with much larger scale lengths (4./1.3/0.8). The disc scale length is about 2.2 to 2.3 kpc with a wide hole of about 1 to 1.2 kpc (but 0.5 kpc only when an E-shape is assumed for the bulge). As a result, the maximum of the disc density appears to be outside the bulge/bar maximum, as if the internal disc was swept out by the bar. The boxiness of the two ellipsoids is noticeable (c and c generally larger than 2) but their values are not very well constrained.

Maps of the star counts in the Ks band are shown in Fig. 3, to be compared with Fig. 2. Residuals are much smaller than in the 1-ellipsoid model, and the fit appears good in nearly all regions with residuals at the level of 10%. The population of the nuclear bar at |l| < 2° as in previous model appears by subtraction. To see where the differences are most significant, accounting for the Poisson noise, we have plotted the difference divided by the Poissonian counting error of the model (i.e. (Nobs − Nmod)/) in the bottom right-hand panel of Fig. 2.

thumbnail Fig. 3

Same as Fig. 2 but with two fitted ellipsoids.

4.5. Attempts with other shapes

As is the case for many disc galaxies, the Milky Way has a spiral structure even if it is not well defined at present. The spiral arms are mainly detected in the visible from young populations and from interstellar matter tracers such as HII regions. 2MASS star counts are dominated by K giants, which generally do not follow the spiral structure because of their age, but the counts could be slightly sensitive to them because in any case, a part of these K giants are indeed young objects. Another structure that has been well studied in CO is the molecular ring where a large proportion of the giant molecular clouds reside. It is not well established whether the stellar populations contribute to this ring. Sevenster & Kalnajs (2001) propose that a stellar ring significantly enhances the microlensing optical depth towards the bulge.

To check this point we investigate the model likelihood when a ring shape is included in the fit in place of the second ellipsoid. For this we use the Sevenster & Kalnajs (2001) elliptical ring density law, with free parameters: ring mean radius, ellipticity, and orientation in the plane, thickness in the plane, scale height, and normalization. We also tested cutoffs in longitude for this ring to allow for the existence of a portion of ring rather than a complete one. The result is that the likelihood is comparable to the 1 ellipsoid model, and the density of such a structure is found to be compatible with 0. It means that in the longitude and latitude zones considered here, the ring and the spiral arms do not significantly contribute to the 2MASS star counts.

thumbnail Fig. 4

Distribution of solutions for the main bulge population (first ellipsoid) with regards to its parameters. Columns are, from the left to the right: angle in degrees, scale length x0, y0, z0 in kpc, central mass, c and c . Rows are given in the same order from top to bottom. Red plus signs indicate all likelihood, squares hold for the 10 best likelihoods, the blue circle for the best solution.

thumbnail Fig. 5

Distribution of solutions for the second ellipsoid with regards to its parameters. Red plus signs indicate all solutions, squares hold for the 10 best likelihood solutions, blue circle for the best solution.

thumbnail Fig. 6

Distribution of solutions and correlation between parameters of the second ellipsoid with regard to the first ellipsoid. The first ellipsoid is in abscissa, second in ordinate. Red plus signs indicate all likelihoods, green squares hold for the 10 best likelihoods, blue circle for the best solution.

4.6. Uncertainties and correlations

Using our method we are able to estimate both the accuracy with which the parameters are determined and the correlations between them. A Monte-Carlo process was computed about 40 times in the case of the best model (at least 20 times in other cases) to compute the average best parameters and give estimates of the uncertainties and correlations between parameters. For the model S+E, correlations between parameters are shown in Fig. 4 for the first ellipsoid and in Fig. 5 for the second one. Cross-correlations between parameters of the two ellipsoids are shown in Fig. 6. In these figures for all trials, the best ten solutions (in the sense of the likelihood) and the best likelihood solution shown in the table. The correlations between parameters of a given population are rather weak in general, indicating that the method is robust. The only net correlation appears between scale lengths for the first ellipsoid (between x0 and y0 and z0) and with the normalization.

The uncertainties on the scale lengths is about 250 pc for the major axis, and only 50 pc on the secondary axis. The normalization is determined at the level of 15% for the main ellipsoid and the second ellipsoid. The boxyness parameters c and c ⊥  are nearly always found to be larger than two, being three in the mean, making both ellipsoids rather boxy, even if the secondary one is found with an exponential shape (which is more like a disc than the sech2).

The disc and hole-scale lengths are slightly anti-correlated, a longer scale length favouring a smaller hole. This is expected due to the uncertainty on the distance indicators in broad band photometry. A better fit of the disc scale length would need to include larger longitudes. But then the spiral structure should also be included, as we would reach arm tangents. The model also has a slight tendancy to produce too many stars at l <  −5° and too few stars at l > 5° close to the plane. We could interpret this by the ability of the hole to be ellipsoidal in the plane, as found by Binney et al. (1997). We attempted to model this ellipticity but it is poorly constrained, and the fit is only very slightly better than with a round hole. This case could be considered later in more detail with deeper data.

5. Extinction validation

The analysis was done assuming a 3D extinction distribution model obtained using the bulge model from Picaud & Robin (2004). We need to ensure that the extinction resulting from the new fit is consistent with the initial extinction we assumed. Figure 7 shows the comparison between the projected extinction map used originally with the map deduced from the fit with the new model. The differences are negligible. We also compared the distribution in distance of AK along some lines of sight and saw generally good agreement in the distances of the clouds, within the error bars of this distance (a few hundred parsecs). Figure 8 shows the variation in the distribution of extinction along four lines of sight for the standard model (old), the model with the bulge population simulated by one ellipsoid (F1), and the model with two ellipsoids (model 2, F2). The differences are generally negligible. In one case (longitude 0), we see that the cloud is found at a greater distance with the new model, but the total extinction in the bulge remains the same, so the extinction deduced using our method (Marshall et al. 2006) is sufficiently robust and not particularly model dependent. Conversely, we can emphasize that the 3D extinction model is reliable enough to avoid biasing our results concerning the stellar populations study.

thumbnail Fig. 7

Ak map deduced from Marshall et al. (2006) method. Top: based on the standard stellar population model from BGM. Bottom: based on the new stellar population model.

thumbnail Fig. 8

Distribution of extinction over four lines of sight (latitudes b =  −3° and longitudes between –3 and 0°) from the orignal model (Old), 1-ellipsoid model (F1), 2-ellipsoid model (F2).

6. Resulting colour–magnitude diagrams and double-clump feature

We present colour-magnitude diagrams for five directions in order to see the subtle differences between 1-ellipsoid and 2-ellipsoid models. To reproduce the natural width of the sequences in the CMD, we have applied a star-by-star variation of the extinction, about the mean extinction at the star distance given by the extinction map. The cosmic variation in the extinction is assumed to be 20% of the extinction, a value that is estimated by trial and error.

Figures 913 show these CMDs. The 2MASS data are in the first columns, 1-ellipsoid model in column 2, 2-ellipsoid model in Cols. 3, and histograms in Ks in Col. 5. In Col. 4 we also present a modified model that produces a double-clump feature (see below). One notices the incompleteness of the data starting probably at about K = 12.5 in fields at |b| ≤  −4°.

The 1-ellipsoid model gives star counts in the red giant branch that are too high. The agreement is generally much better with the 2-ellipsoid model, even if it is not perfect in some fields. There are still residuals in the shape of the giant branch that could probably be solved by varying the star by star dispersion in the extinction, especially in the case of the field with coordinates (l,b) = (5;   4.5). The important point here is to have a general model reproducing most fields, at the expense of some discrepancy in a limited number of fields, since we are interested in the large-scale structure.

In fields at intermediate latitudes (|b| > 5°) several studies have shown there are double clumps. These double clumps change as a function of longitude (McWilliam & Zoccali 2010). From our analysis limited to 2MASS differential star counts, and due to the lack of completeness of the data, our solution does not present double clumps, but just a modulation of the shape of the luminosity function due to the shape of the bar. The second structure has its clump peak at fainter magnitudes than the 2MASS completeness limit. If there is a structure creating the double clump it should be in the main bar.

Since our model reproduces the counts at the level of 10% in the whole bulge region, we can consider that what is creating the double clump feature is a modulation of the whole shape of the bar. Pfenniger (1985) did a numerical study of instability in the 3D family of periodic orbits in a barred galaxy model. He shows a family of resonance creating the growth of a box-shaped bulge. It shows how stars are diffused at high z. It is natural to hypothesize that the double clump can be produced by such an effect.

To test this hypothesis, we simulated a flaring bar by enhancing the scale height of the bar with a sinus fonction of the distance along the major axis of the bar. If x is the galactocentric distance along the major axis, and dz0 is the scale height along the z axis, the new vertical scale height is simulated by

The original value of dz0 (obtained by the fit above) was decreased slightly (from 390 to 330 pc) to account for the fact that we measured the mean scale height when we fixed this parameter in the fit. The parameters in the equation are a reasonable guess in order to fit the clump data at a medium latitude.

In Fig. 14 we show the star counts around the clump in various fields at latitude b =  −8° in order to compare with Fig. 3 of McWilliam & Zoccali (2010). It shows that a flaring bar as modelled here can qualitatively reproduce the double-clump feature. The proportion of stars in each clump varies as a function of the longitude due to the flare, and the significant differences between positive and negative longitudes. This figure can be compared with Fig. 3 of McWilliam & Zoccali (2010). It can also be seen on a zoom of the luminosity function in the field at l = 0, b =  −7° (Fig. 15) where we superimposed the observed counts with the model predictions without the flare and with the flare. We show that a simple modulation of our primary model with a slight bar flare reproduces the two-clump feature.

Since the 2MASS data are not complete, we cannot perform a quantitative fit of the flare parameters with these data. We only show here that a flaring bar can be a realistic explanation for the double-clump and that it is dynamically reasonable, as in the Pfenniger (1985) simulations. In a future paper we plan a more complete analysis of the shape of the flare using VVV data, for example, by comparing simulations with the density plots from Saito et al. (2011) in order to specify the characteristics of this flare more accurately as a function of the galacto-centric distance.

thumbnail Fig. 9

Colour–magnitude diagrams for the 2MASS field at l = 0, b =  −4. a) data; b) best-fit model with 1 ellipsoid; c) best-fit models with 2 ellipsoids; d) modified model with flared bar; e) histograms of data (red solid) and models (1 ellipsoid: green long dashed; 2 ellipsoids: blue dotted, flared bar: magenta short dashed).

thumbnail Fig. 10

Colour–magnitude diagrams for the 2MASS field at l =  + 3, b =  −3. Same coding as in Fig. 9.

thumbnail Fig. 11

Colour–magnitude diagrams for the 2MASS field at l = 5, b = 4.5. Same coding as in Fig. 9.

thumbnail Fig. 12

Colour–magnitude diagrams for the 2MASS field at l =  −17, b =  −4.25. Same coding as in Fig. 9.

thumbnail Fig. 13

Colour–magnitude diagrams for the 2MASS field at l = 0, b =  −8. Same coding as in Fig. 9. In this field we see clearly in the histogram the double clump feature reproduced by the modified model (panel d) and magenta line) and not by the originally fitted model (panel c) and blue line).

7. The inner bulge

Assuming that the model of stellar populations and 3D extinction now reproduces the populations of the outer bulge and the disc, we can investigate if any further structures are missing from our modelling. As a check we plot in figure 16 the higher resolution map of the difference between model predictions and observed star counts in the region  −4°  < l < 4° and  −3°  < b < 3°. The excess of stars in the data should represent the population of the inner bulge. It resembles the population seen by Alard (2001) and coincide with the position of the central molecular zone. Several authors (Morris & Serabyn 1996; An et al. 2011) have shown that this region contains young stellar populations. To investigate further this population, one would need to improve the 3D extinction model here. As we saw in Fig. 7, the extinction determined by our method saturates at about AK = 3 due to the limiting magnitude in 2MASS. Thus our extinction is most probably underestimated, implying that the density of the inner bulge population is overestimated. We do not investigate further the density of this population in this paper, as it could only be qualitative from 2MASS data alone. We envisage to apply our method to deeper data and at larger wavelengths, like Spitzer at 3.6 and 4.5 microns, in order to make a deeper 3D map of the extinction in this region and to study the inner bulge region, the probable nuclear bar or ring and its link with the central molecular zone.

8. Discussion

We have shown in this analysis that the bulge region can be simulated using the sum of two ellipsoids: a main component that is a quite standard boxy bulge/bar, which dominates the counts up to latitudes of about 5°, and a second ellipsoid that is a thicker structure seen mainly at higher latitudes. This thick component seems to be about as thick as the local thick disc. It could be either the classical bulge, flattened by the effect of the barred potential, or the inner counter part of the thick disc population. To distinguish how these structures can be differentiated and to understand which population they belong to, it is useful to look at the spectroscopic surveys showing metallicity distributions at different latitudes.

8.1. Metallicity as a function of latitude

Among studies measuring the metallicity distribution function (MDF) in different fields of the bulge region several have considered the characteristics of the populations along the minor axis at different latitudes. Zoccali et al. (2008) studied the MDF and radial velocity distribution from the RGB, Gonzalez et al. (2011) completed the data with the abundances of alpha abundances, Babusiaux et al (2010) studied the correlation between metallicity and kinematics and Hill et al. (2011) computed the MDF from the red clump giants for Baade’s window. The overall conclusion of these studies is that the MDF looks bimodal and the contribution of the components changes with distance from the plane. At a low latitude (–4°), the dominant population has a metallicity slightly above-solar and a small contribution comes from a second component of metallicity around –0.5. Babusiaux et al (2010) show that the high-metallicity component has a high vertex deviation compared to the low-metallicity component that does not. They conclude that the vertical metallicity gradient is the effect of mixing of populations. On the other hand, Ness et al. (2011) decompose the bulge populations from an analysis of the ARGOS survey in up to five components to reproduce the complete MDF. They find a smaller gradient (about 0.4 dex/kpc) compared with Zoccali et al. (2008), who found a mean gradient of 0.6 dex/kpc.

Can we explain these MDFs by our two-ellipsoid model? To answer this question, we simulated the observed MDF from Zoccali et al. (2008) in three latitude ranges, on the minor axis and deduced what would be the necessary metallicity of each of our components to explain the observations.

thumbnail Fig. 14

Luminosity functions for the red clump region in fields at latitude –8°. Coordinates in longitude and latitude are indicated in each panel. The double-clump feature is simulated by adding a slight flare to the original bar model (see text).

thumbnail Fig. 15

Luminosity functions for the red clump region in a field at latitude –7°and longitude 0° from 2MASS (red dots with error bars) compared with the model without flare (in green solid line) and with the flaring bar (in blue dotted line).

thumbnail Fig. 16

Difference map for the inner bulge region (Nobs − Nmod) / Nobs. The population in excess in the data is most likely the nuclear bar.

We assume that the boxy bar has a mean solar metallicity and that the thick component has a mean metallicity of –0.35 dex. We did not fit the density of each population, which is the one obtained from the 2MASS star-count fitting. We applied the selection function explained in Zoccali et al. (2008) for the RGB and in Hill et al. (2011) for the red clump sample. Figure 17 shows the comparison between the observed metallicity distribution and the simulated one for RGB giants at l = 0° and b =  −4°, –6°, and –12°, and for red clump giants at b =  −4°. At the bottom we also show the decomposition of the distribution in the three major components: the thin disc, the “boxy-bar”, and the “thick-bulge”. We see good agreement beween our model predictions and the data in all four cases, and even clearer at high latitudes, within the small Poisson statistics of the samples. The error bars are only the Poisson noise on the counts in the selected sample. It does not take observational errors on the metallicity into account or the noise due to the selection function. The only significant disagreement comes from a lower number of low-metallicity stars in the model for the field at b =  −6° compared to the data. This feature will be explored in more detail as soon as a larger sample will be available. It will then be possible to determine whether a metallicity gradient should be added to one of the two populations or whether the relative density of the two populations is represented well by this model or should be revised given the errors.

thumbnail Fig. 17

Metallicity distribution for RGB samples from Zoccali et al. (2008) in fields at l = 0° and b =  −4°, –6°, and –12°, and for the red clump samples from Hill et al. (2011) in Baade’s window (columns from left to right). Top: the observed data (green dots with Poisson noise) are compared with model simulations, applying the same selection function as in the observations (red solid lines). Bottom: decomposition in populations of the model simulation. Red solid line: total; green long dash: disc, short dash blue: boxy bar; dotted magenta: thick bulge. The error bars take only the Poisson noise into account and not the observational errors or the metallicity or the error in the selection function estimation.

8.2. Conclusion and future plans

We have presented here an attempt to fit differential star counts, CMDs, and metallicity distribution functions in a wide area over the Milky Way bulge region. The model accounts for the metallicity distribution on the minor axis of the bulge, the dissymmetry of the star counts as a function of longitude, and the existence of double clumps at medium latitudes. This model has two components. The main one is modelled by boxy, S-shape triaxial Ferrer ellipsoids, having scale lengths of 1.46/0.49/0.39 kpc, an angle of the main axis with regards to the Sun-Galactic centre direction of 13°, and a total mass of 6.1 × 109 M. The age is assumed to be 8 Gyr or more, and the mean metallicity is approximately solar. The second component is modelled by a triaxial exponential Ferrer ellipsoid less massive (2.6 × 108M, with scale lengths of 4.44/1.31/0.80 kpc, also old with a poorer metallicity of –0.35 dex. A slight flare of the bar component is able to explain the double clumps qualitativeley at latitudes higher than 6°.

This model with two components explains the observed metallicity distribution function and its gradient along the minor axis by a transition between the two components dominating the counts, and also better reproduces the detailed distribution in the CMDs. A preliminary comparison with Brava data Rich et al. (2007) shows that the main component is a fast rotator, so is probably a boxy bar (the object of a paper in preparation). The second structure (the thicker bulge) has a smaller velocity dispersion, thus could be either a classical bulge flattened by the potential of the bar or formed by early mergers, as in the Bournaud et al. (2007) scenario. Testing of these scenarios are on-going, using our model and comparing predictions with kinematical and metallicity data in the whole bulge region. The forthcoming surveys, APOGEE project Majewski et al. (2007), Gaia-ESO spectroscopic survey, among others, and later the ESA Gaia mission (http://www.rssd.esa.int/index.php?project=GAIA) are all expected to offer new insights into the kinematics and abundances of the central region of the Galaxy, giving strong constraints on any scenario for the formation of the bulge and bar stellar populations. We emphasize that our model can be used for preparing these future surveys and will be a useful tool for their interpretation.

Acknowledgments

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. The CDSClient package was used for the remote querying of the 2MASS dataset. This material is based upon work supported in part by the National Science Foundation under Grant No. 1066293 and the hospitality of the Aspen Center for Physics. Simulations were executed on computers from the Utinam Institute of the Université de Franche-Comté, supported by the Région de Franche-Comté and Institut des Sciences de l’Univers (INSU). We acknowledge the French spatial agency, the Centre National d’Étude Spatiale, which has provided funding for D. J. Marshall, and the support of the french Agence Nationale de la Recherche under contract ANR-2010-BLAN-0508-01OTP. We thank Lia Athanassoula for fruitful discussions and Carine Babusiaux and Vanessa Hill for providing their data and advices.

References

  1. Alard, C. 2001, A&A, 379, L44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. An, D., Ramírez, S. V., Sellgren, K., et al. 2011, ApJ, 736, 133 [NASA ADS] [CrossRef] [Google Scholar]
  3. Babusiaux, C., & Gilmore, G. 2005, MNRAS, 358, 1309 [NASA ADS] [CrossRef] [Google Scholar]
  4. Babusiaux, C., Gomez, A., Hill, V., et al. 2010, A&A, 519, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Benjamin, R. A., Churchwell, E., Babler, B. L., et al. 2005, ApJ, 630, L149 [NASA ADS] [CrossRef] [Google Scholar]
  6. Binney, J., Gerhard, O., & Spergel, D. 1997, MNRAS, 288, 365 [NASA ADS] [CrossRef] [Google Scholar]
  7. Blitz, L., & Spergel, D. N. 1991, ApJ, 379, 631 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237 [Google Scholar]
  9. Bruzual, G., Barbuy, B., Ortolani, S., et al. 1997, AJ, 114, 1531 [NASA ADS] [CrossRef] [Google Scholar]
  10. Cabrera-Lavers, A., González-Fernández, C., Garzón, F., Hammersley, P. L., & López-Corredoira, M. 2008, A&A, 491, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Cabrera-Lavers, A., Hammersley, P. L., González-Fernández, C., et al. 2007, A&A, 465, 825 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Einasto, J. 1979, in The Large Scale Characteristics of the Galaxy, ed. W. B. Burton, IAU Symp. 84 [Google Scholar]
  13. Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716 [CrossRef] [Google Scholar]
  14. Epchtein, N., de Batz, B., Capoani, L., et al. 1997, The Messenger, 87, 27 [NASA ADS] [Google Scholar]
  15. Ferrers, N. M. 1877, Quart. J. Pure Appl. Math., 14, 1 [Google Scholar]
  16. Freudenreich, H. T. 1998, ApJ, 492, 495 [NASA ADS] [CrossRef] [Google Scholar]
  17. Fux, R. 1997, A&A, 327, 983 [NASA ADS] [Google Scholar]
  18. Fux, R. 1999, A&A, 345, 787 [NASA ADS] [Google Scholar]
  19. Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044 [NASA ADS] [CrossRef] [Google Scholar]
  20. Girardi, L., Bertelli, G., Bressan, A., et al. 2002, A&A, 391, 195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2011, A&A, 530, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Haywood, M., Robin, A. C., & Crézé, M. 1997a, A&A, 320, 428 [NASA ADS] [Google Scholar]
  23. Haywood, M., Robin, A. C., & Crézé, M. 1997b, A&A, 320, 440 [NASA ADS] [Google Scholar]
  24. Hill, V., Lecureur, A., Gómez, A., et al. 2011, A&A, 534, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Howard, C. D., Rich, R. M., Clarkson, W., et al. 2009, ApJ, 702, L153 [NASA ADS] [CrossRef] [Google Scholar]
  26. López-Corredoira, M., Hammersley, P. L., Garzón, F., Simonneau, E., & Mahoney, T. J. 2000, MNRAS, 313, 392 [NASA ADS] [CrossRef] [Google Scholar]
  27. López-Corredoira, M., Hammersley, P. L., Garzón, F., et al. 2001, A&A, 373, 139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. López-Corredoira, M., Cabrera-Lavers, A., & Gerhard, O. E. 2005, A&A, 439, 107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Majewski, S. R., Skrutskie, M. F., Schiavon, R. P., et al. 2007, BAAS, 39, 962 [NASA ADS] [Google Scholar]
  30. Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Martinez-Valpuesta, I., & Gerhard, O. 2011, ApJ, 734, L20 [NASA ADS] [CrossRef] [Google Scholar]
  32. McWilliam, A., & Zoccali, M. 2010, ApJ, 724, 1491 [NASA ADS] [CrossRef] [Google Scholar]
  33. Morris, M., & Serabyn, E. 1996, ARA&A, 34, 645 [NASA ADS] [CrossRef] [Google Scholar]
  34. Nataf, D. M., Udalski, A., Gould, A., Fouque, P., & Stanek, K. Z. 2010, ApJ, 721, L28 [NASA ADS] [CrossRef] [Google Scholar]
  35. Ness, M., et al. 2011, In Assembling the puzzle of the Milky Way, Proc. Conference held in Grand-Bornand, France, on 18th to 22th april, ed. C. Reylé, et al. [Google Scholar]
  36. Nishiyama, S., Nagata, T., Baba, D., et al. 2005, ApJ, 621, L105 [NASA ADS] [CrossRef] [Google Scholar]
  37. Pfenniger, D. 1985, A&A, 150, 112 [NASA ADS] [Google Scholar]
  38. Picaud, S., & Robin, A. C. 2004, A&A, 428, 891 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992 (Cambridge: University Press), 2nd ed. [Google Scholar]
  40. Rattenbury, N. J., Mao, S., Sumi, T., & Smith, M. C. 2007, MNRAS, 378, 1064 [NASA ADS] [CrossRef] [Google Scholar]
  41. Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Rich, R. M. 1988, AJ, 95, 828 [NASA ADS] [CrossRef] [Google Scholar]
  43. Rich, R. M., & Origlia, L. 2005, ApJ, 634, 1293 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  44. Rich, R. M., Reitzel, D. B., Howard, C. D., & Zhao, H. 2007, ApJ, 658, L29 [NASA ADS] [CrossRef] [Google Scholar]
  45. Robin, A.C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Rodríguez-Fernández, N. J. 2011, Mem. Soc. Astron. Ital. Suppl., 18, 195 [NASA ADS] [Google Scholar]
  47. Romero-Gomez, M., Athanassoula, E., Antoja, T., & Figueras, F. 2011, MNRAS, 418, 1176 [NASA ADS] [CrossRef] [Google Scholar]
  48. Saito, R. K., Zoccali, M., McWilliam, A., et al. 2011, AJ, 142, 76 [NASA ADS] [CrossRef] [Google Scholar]
  49. Sevenster, M. N., & Kalnajs, A. J. 2001, AJ, 122, 885 [NASA ADS] [CrossRef] [Google Scholar]
  50. Shen, J., Rich, R. M., Kormendy, J., et al. 2010, ApJ, 720, L72 [NASA ADS] [CrossRef] [Google Scholar]
  51. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
  52. Stanek, K. Z., Mateo, M., Udalski, A., et al. 1994, ApJ, 429, L73 [NASA ADS] [CrossRef] [Google Scholar]
  53. Sumi, T., Wu, X., Udalski, A., et al. 2004, MNRAS, 348, 1439 [NASA ADS] [CrossRef] [Google Scholar]
  54. Weiland, J. L., Arendt, R. G., Berriman, G. B., et al. 1994, ApJ, 425, L81 [NASA ADS] [CrossRef] [Google Scholar]
  55. Zoccali, M., Hill, V., Lecureur, A., et al. 2008, A&A, 486, 177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1

Best fit models with 1 single bulge ellipsoid.

Table 2

Parameters and reduced likelihood Lr when two ellipsoids are fitted, including main sequence and giants, and forcing the two ellipsoids at the same orientation.

All Figures

thumbnail Fig. 1

Distribution in longitude and latitude of the 200 windows used in the fitting process. The colour coding corresponds to the limiting magnitude of the counts in Ks.

In the text
thumbnail Fig. 2

Star counts up to limiting magnitude Ks from 2MASS data (top left) compared with fitted 1 ellipsoid model B (top right), residuals (Nmod − Nobs)/Nobs (bottom left), and difference divided by the Poissonian counting error of the model (i.e. (Nobs − Nmod)/ (bottom right). In the residual map, contours are drawn at intervals of 20% model overestimate (black) and 20% model underestimate (white). The 1-ellipsoid model leaves significant X-shaped residuals. Near the Galactic centre the nuclear bar population is missing in the model. The residuals in the outer regions are not very significant due to the small number of stars in each bin (seen in the top left panel in dark violet and black).

In the text
thumbnail Fig. 3

Same as Fig. 2 but with two fitted ellipsoids.

In the text
thumbnail Fig. 4

Distribution of solutions for the main bulge population (first ellipsoid) with regards to its parameters. Columns are, from the left to the right: angle in degrees, scale length x0, y0, z0 in kpc, central mass, c and c . Rows are given in the same order from top to bottom. Red plus signs indicate all likelihood, squares hold for the 10 best likelihoods, the blue circle for the best solution.

In the text
thumbnail Fig. 5

Distribution of solutions for the second ellipsoid with regards to its parameters. Red plus signs indicate all solutions, squares hold for the 10 best likelihood solutions, blue circle for the best solution.

In the text
thumbnail Fig. 6

Distribution of solutions and correlation between parameters of the second ellipsoid with regard to the first ellipsoid. The first ellipsoid is in abscissa, second in ordinate. Red plus signs indicate all likelihoods, green squares hold for the 10 best likelihoods, blue circle for the best solution.

In the text
thumbnail Fig. 7

Ak map deduced from Marshall et al. (2006) method. Top: based on the standard stellar population model from BGM. Bottom: based on the new stellar population model.

In the text
thumbnail Fig. 8

Distribution of extinction over four lines of sight (latitudes b =  −3° and longitudes between –3 and 0°) from the orignal model (Old), 1-ellipsoid model (F1), 2-ellipsoid model (F2).

In the text
thumbnail Fig. 9

Colour–magnitude diagrams for the 2MASS field at l = 0, b =  −4. a) data; b) best-fit model with 1 ellipsoid; c) best-fit models with 2 ellipsoids; d) modified model with flared bar; e) histograms of data (red solid) and models (1 ellipsoid: green long dashed; 2 ellipsoids: blue dotted, flared bar: magenta short dashed).

In the text
thumbnail Fig. 10

Colour–magnitude diagrams for the 2MASS field at l =  + 3, b =  −3. Same coding as in Fig. 9.

In the text
thumbnail Fig. 11

Colour–magnitude diagrams for the 2MASS field at l = 5, b = 4.5. Same coding as in Fig. 9.

In the text
thumbnail Fig. 12

Colour–magnitude diagrams for the 2MASS field at l =  −17, b =  −4.25. Same coding as in Fig. 9.

In the text
thumbnail Fig. 13

Colour–magnitude diagrams for the 2MASS field at l = 0, b =  −8. Same coding as in Fig. 9. In this field we see clearly in the histogram the double clump feature reproduced by the modified model (panel d) and magenta line) and not by the originally fitted model (panel c) and blue line).

In the text
thumbnail Fig. 14

Luminosity functions for the red clump region in fields at latitude –8°. Coordinates in longitude and latitude are indicated in each panel. The double-clump feature is simulated by adding a slight flare to the original bar model (see text).

In the text
thumbnail Fig. 15

Luminosity functions for the red clump region in a field at latitude –7°and longitude 0° from 2MASS (red dots with error bars) compared with the model without flare (in green solid line) and with the flaring bar (in blue dotted line).

In the text
thumbnail Fig. 16

Difference map for the inner bulge region (Nobs − Nmod) / Nobs. The population in excess in the data is most likely the nuclear bar.

In the text
thumbnail Fig. 17

Metallicity distribution for RGB samples from Zoccali et al. (2008) in fields at l = 0° and b =  −4°, –6°, and –12°, and for the red clump samples from Hill et al. (2011) in Baade’s window (columns from left to right). Top: the observed data (green dots with Poisson noise) are compared with model simulations, applying the same selection function as in the observations (red solid lines). Bottom: decomposition in populations of the model simulation. Red solid line: total; green long dash: disc, short dash blue: boxy bar; dotted magenta: thick bulge. The error bars take only the Poisson noise into account and not the observational errors or the metallicity or the error in the selection function estimation.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.