A parametric description of the 3D structure of the Galactic bar/bulge using the VVV survey

Abstract

We study the structure of the inner Milky Way using the latest data release of the VISTA Variables in the Via Lactea (VVV) survey. The VVV is a deep near-infrared, multi-colour photometric survey with a coverage of 300 square degrees towards the bulge/bar. We use red clump (RC) stars to produce a high-resolution dust map of the VVV's field of view. From de-reddened colour–magnitude diagrams, we select red giant branch stars to investigate their 3D density distribution within the central 4 kpc. We demonstrate that our best-fitting parametric model of the bulge density provides a good description of the VVV data, with a median percentage residual of 5 per cent over the fitted region. The strongest of the otherwise low-level residuals are overdensities associated with a low-latitude structure as well as the so-called X-shape previously identified using the split RC. These additional components contribute only ∼5 per cent and ∼7 per cent respectively to the bulge mass budget. The best-fitting bulge is ‘boxy’ with an axial ratio of [1:0.44:0.31] and is rotated with respect to the Sun–Galactic Centre line by at least 20°. We provide an estimate of the total, full sky, mass of the bulge of |$M_\mathrm{bulge}^{\mathrm{Chabrier}} = 2.36 \times 10^{10} \,\,\mathrm{M}_{{\odot }}$| for a Chabrier initial mass function. We show that there exists a strong degeneracy between the viewing angle and the dispersion of the RC absolute magnitude distribution. The value of the latter is strongly dependent on the assumptions made about the intrinsic luminosity function of the bulge.

1 INTRODUCTION

Integrated light, star counts, gas kinematics and microlensing studies over the past 20 years have provided plenty of evidence that the heart of our Galaxy harbours a triaxial entity that is distinct from the disc (see e.g. Binney et al. 1991; Blitz & Spergel 1991; Dwek et al. 1995; Binney, Gerhard & Spergel 1997; Stanek et al. 1997; Freudenreich 1998; Alcock et al. 2000; Bissantz & Gerhard 2002; Evans & Belokurov 2002; Babusiaux & Gilmore 2005; McWilliam & Zoccali 2010; Nataf et al. 2010). Dominating the dynamics in the inner regions of the Milky Way (MW), this centrally concentrated and elongated structure resembles both a bar and a bulge, although not in the classical sense for the latter (see Kormendy & Kennicutt 2004, for details). However, as the intention of this paper is not to re-enact the decades-old battle between the bar and the bulge, we will refer to the central 4 kpc region of the Galaxy as the ‘bulge’, as have many authors before us (see e.g. Zoccali & Valenti 2016).

Near-infrared (NIR) bands, where the extinction is only around 10 per cent of that in the optical, are ideal for studying the dust-obscured Galactic Centre (GC). Therefore, infrared maps from surveys such as COBE/DIRBE (Diffuse Infrared Background Experiment) have long been used to construct 3D models of the bulge: Dwek et al. (1995) and Freudenreich (1998) have fitted parametric models while Binney et al. (1997) and Bissantz & Gerhard (2002) have obtained non-parametric models using the same data. The 2MASS (Skrutskie et al. 2006) catalogue added depth and resolution to the picture of the inner MW and helped to provide tight constraints on models of the GC (Robin et al. 2012). More recently, the VISTA Variables in the Via Lactea (VVV) project (Minniti et al. 2010) has staked a claim as the deepest wide-coverage NIR survey aimed at delivering an accurate 3D map of the Galactic bulge. The VVV survey has already been used to study the split red clump (RC) in the bulge (Saito et al. 2011), create a high-resolution extinction map of the region (Gonzalez et al. 2012), discover new Galactic star clusters (Borissova et al. 2011) and, most relevant to this study, probe the structure of the Galactic bar (Wegg & Gerhard 2013).

The work of Wegg & Gerhard (2013) is one of the latest in the long line of red giant branch (RGB) density modelling exercises that have helped to reveal the structural properties of the bulge region (e.g. Stanek et al. 1994, 1997; Nikolaev & Weinberg 1997; Unavane & Gilmore 1998; López-Corredoira et al. 2000; Benjamin et al. 2005; Nishiyama et al. 2005; Rattenbury et al. 2007). Many of the works above relied in particular on the RC population to disentangle the bulge from the disc, taking advantage of the well-defined location of the RC in colour–magnitude diagrams (CMD) and its superb performance as a standard candle (Paczyński & Stanek 1998). Importantly, the RC star counts analysis has so far produced the most consistent range of values for the viewing angle of the ellipsoid in the bulge: between 20° and 30°, while the same quantity estimated with other tracers appears to have a much larger uncertainty, with values as low as 10° and as high as 50° quoted (see e.g. Vanhollebeke, Groenewegen & Girardi 2009, for a summary). Of course, it is quite likely that this apparent dramatic variation in the basic property of the bar/bulge is not solely the sign of inconsistency between the methods but also an indication that the angle does change depending on the tracer used and the range of Galactic l and b explored. For example, there exists ample evidence for a long bar, clearly seen outside the bulge at larger |l|, with a viewing angle in the range 28°–45° (for a recent study, see Wegg, Gerhard & Portail 2015, and references therein). Similarly, the structure traced with the RR Lyrae might have very little triaxiality, if any at all, thus most likely representing the actual (classical) old bulge of the MW (see Dékány et al. 2013; Pietrukowicz et al. 2015; Kunder et al. 2016).

In this study, we use the most recent data release of the VVV and strive to find the most appropriate analytic function that best describes the full 3D bulge density distribution. Therefore, our adopted approach is different from that employed by e.g. Wegg & Gerhard (2013) who applied instead a symmetrized non-parametric modelling procedure to an earlier version of the same VVV data. There are clearly advantages and disadvantages associated with both parametric and non-parametric methods. On the one hand, describing the stellar distribution in the bulge with an analytic expression clearly gives a more portable solution that can be straightforwardly used in subsequent dynamical modelling of the influences of the bar/bulge. On the other hand, the actual density field in the central Galaxy might not be of a simple triaxial shape. This of course does indeed take place in the MW and is exemplified by the recently discovered cross-like component of the bulge as manifested by the split RC (see e.g. McWilliam & Zoccali 2010; Nataf et al. 2010).

The paper is organized as follows: the VVV data used are described in Section 2; the procedure used to build a high-resolution extinction map is detailed in Section 3; Section 4 contains a description of the parametric approach used to investigate the bulge stellar density law using a mock catalogue; in Section 5 we outline the fitting procedure; in Section 6 we explain the models used and provide the best-fitting results. Finally, in Sections 7 and 8, we discuss our results and draw the conclusions.

2 THE VVV SURVEY

In 2016, the VVV survey finished the 6 year long observing campaign and although many publications have already stemmed from it, it remains a largely unexplored resource. The VVV survey contains a stellar sample large enough to examine the structural properties of the bulge in unprecedented detail, in the J, H, K_s NIR bands and reaches up to 4 mag deeper than 2MASS. The depth of the survey allows us to select the RGB population visible across the whole survey footprint, apart from some highly reddened regions close to the Galactic plane, and thus build a 3D parametric model to describe their density distribution towards the bulge.

The VVV survey covers ∼315 square degrees of the inner MW, i.e. a region with Galactic longitude −10° < l < +10° and latitude −10° < b < 5°, as well as 220 square degrees of the Galactic plane, spanning between 295° < l < 350° and −2° < b < 2°. In this work, we exclusively use the Galactic bulge part of the survey. VISTA data are processed through the VISTA Data Flow System by the Cambridge Astronomical Survey Unit (CASU¹).

The photometric and astrometric calibration of the VVV survey is performed relative to the 2MASS point source catalogue ensuring that the calibration stars are measured at the same time as the target sources (see e.g. González-Fernández et al. 2017). For the J, H and K_s passbands employed here, typical overall calibration errors are around 1–2 per cent. However, the accuracy of these calibrations is field dependent and can be more uncertain in very crowded fields with high levels of spatially varying extinction, where the derived zero-points can vary from field to field by up to 0.1 mag. We therefore performed a full global quality check using overlap regions for every tile to facilitate correcting tiles with offset zero-points and/or unusual CMD.

3 EXTINCTION MAP

The biggest challenge to building an accurate map of the bulge stellar density arises from the rapid extinction variations across different lines of sight. Following a similar approach to the one described in Gonzalez et al. (2011), we use the RC giants in the VVV to build a high-resolution (1 arcmin × 1 arcmin sampling) reddening map (see Fig. 1) by comparing the mean colour of the RC stars observed in each field with that of the RC population in Baade's window, |$(J^{\prime }-K^{\prime }_{s})_{\mathrm{BW}} = 0.89$| (throughout the paper, we denote with J΄ and |$K^{\prime }_{s}$| the observed VVV magnitudes and with J and K_s the extinction-corrected magnitudes), under the assumption that the intrinsic colour of the RC remains, on average, the same across the bulge, i.e. J − K_s = 0.62. We adopt the Nishiyama et al. (2009) interstellar extinction law derived for the central bulge region since it was shown to reproduce most accurately the extinction of RC stars in the VVV (Gonzalez et al. 2012).

To build the extinction map, we combine the 2MASS–VVV colour transformations |$J^{\prime } = J^{\prime }_{\mathrm{2M}} - 0.065(J^{\prime }-K^{\prime }_{s})_{\mathrm{2M}}$| and |$K^{\prime }_{s} = K^{\prime }_{s\mathrm{2M}} + 0.010(J^{\prime }-K^{\prime }_{s})_{\mathrm{2M}}$| and obtain the Nishiyama et al. (2009) law for VVV photometry: |$A_{J} = 1.351E(J^{\prime }-K^{\prime }_{s})$| and |$A_{Ks} = 0.482E(J^{\prime }-K^{\prime }_{s})$|⁠. We then calculate the reddening |$E(J^{\prime } - K^{\prime }_{s})_{\mathrm{RC}}$| by comparing the observed RC colour in each field with its known extinction-free value in Baade's window, J − K_s = 0.62. We are able to use this technique in all bulge fields covered by the survey as VVV photometry reaches the RC even in fields close to the Galactic plane. Correction for extinction then proceeds on a star-by-star basis interpolating between the nearest four extinction map values.

The problem with using two-dimensional maps such as ours is the assumption that the extinction arises mainly from foreground material while in reality the sources are intermixed with dust at a range of distances, rendering the use of 3D reddening maps necessary especially at low latitudes (|b| < 3°). To compare the accuracy of 2D versus 3D dust maps, Schultheis et al. (2014) used a large spectroscopic sample of stars from the Apache Point Observatory Galactic Evolution Experiment (APOGEE) towards the MW bulge, to conclude that for sources at distances larger than ∼4 kpc from the Sun, 2D maps can be applied without significant systematic offset. Moreover, it was noticed that none of the existing 3D maps agrees with independent extinction calculations based the APOGEE spectra towards the MW bulge. Their study suggests that there is still room for improvement in the construction of 3D dust maps and that they do not outperform 2D maps for distances larger than 4 kpc from the Sun.

4 MODEL INGREDIENTS

In this section, we describe the main ingredients of our model that consists of synthetic thin and thick discs and a purely analytic bulge. The model takes into account the disc contribution to the observed luminosity function (LF) as a function of field of view and magnitude and therefore has an advantage over recent VVV studies that either fitted a polynomial background to the discs and RGB (from both bulge and discs) populations (Wegg & Gerhard 2013) or assumed a uniform disc contamination throughout the bulge area (Valenti et al. 2016).

The disc populations are built with galaxia² (Sharma et al. 2011), a C++ code able to generate synthetic surveys of the MW, implementing the Besançon³ Galaxy model (Robin et al. 2003).

4.1 Isochrones

We add to galaxia a recent set of PARSEC isochrones in the VISTA photometric system (Bressan et al. 2012), built using the default parameters of the online Padova CMD interface (we used PARSEC v1.2S + COLIBRI PR16;⁴ described in Marigo et al. 2017). For the table of isochrones, we choose 177 age bins and 39 metallicity bins, within the parameters range used in galaxia for a previous PARSEC release (Marigo et al. 2008).

VVV zero-point measurements (section 2.2 in Rubele et al. 2012) found a constant offset between the VISTA system calibration and the Vegamag system of the stellar models of 0.021 and 0.001 mag for the J and K_s bands, respectively. This offset was implemented in the PARSEC version we use.

4.2 The discs

The thin and thick disc populations are an important component of the central regions; they are implemented by averaging several repeated mock catalogue realizations with galaxia to minimize sampling noise. The disc models are shown in Fig. 2, in Galactic coordinates – these are mock stars within the same range of extinction-corrected magnitudes and colours, 12 < K_s < 14, 0.4 < J − K_s < 1.0, used for the selected VVV sample (see Section 5). Overlaid we mark two contours of constant reddening |$E(J^{\prime }-K^{\prime }_{s}) = 1.5$| (light grey) and |$E(J^{\prime }-K^{\prime }_{s}) = 1.0$| (thick grey), as in Fig. 1. The initial mass function (IMF), star formation rate, ages, metallicity and density laws of the discs implemented in galaxia are described in tables 1–3 in Robin et al. (2003). The thin and thick discs are assumed to have a warp and a flare that are modelled following the prescription of Robin et al. (2003).

4.3 The bulge LF

4.3.1 Analytic description of the bulge LF

We assume a 10 Gyr old (Bruzual et al. 1997) bulge population with solar metallicity [Fe/H] = 0 dex (e.g. Zoccali et al. 2003) and a metallicity dispersion of σ_[Fe/H] = 0.4 (Zoccali et al. 2008) with no metallicity gradients (default values of the Besançon Galaxy model).

With these specifications, we simulate the LF with galaxia using a Chabrier IMF, within the colour range 0.4 < J − K_s < 1.0, adopted to reduce contamination from non-giants in the VVV data. The Chabrier IMF is chosen over the galaxia default Salpeter IMF as it was shown to describe more accurately observations and model predictions in the Galactic bulge (see Chabrier 2003 and Portail et al. 2015, respectively). The resulting synthetic LF, ϕ_B(M_K; 10 Gyr), is shown in the third panel of Fig. 3 in red. We model it following the description in section 3.1 of Nataf et al. (2015), with three Gaussians to represent the asymptotic giant branch bump (AGBB), the RC and the red giant branch bump (RGBB) and an exponential function for the RGB stars. The RC was modelled with either a Gaussian or a skew-Gaussian distribution, with the latter providing a slightly better fit; the addition of the skewness parameter affects only the μ^RC and σ^RC values, as expected (see the results in Table 1). The only constraint was on the mean RGBB magnitude, set to μ^RGBB = −1.02, where the peak of the LF lies in the RGBB region. The fitted LF, with a skew-Gaussian model for the RC, is shown in green in the third panel of Fig. 3, overlapping the histogrammed LF simulated with galaxia (in red). We find μ^RC = −1.53, consistent with previous determinations of the mean RC magnitude in the K band (see table 1 in the RC review paper by Girardi 2016).

The first panel of the figure shows the number density distribution of the galaxia generated stars in a CMD with the specifications mentioned above, and the second panel the mean metallicity in each of those pixels. Notice that the theoretical models predict the AGBB, RC and RGBB mean magnitudes depend on metallicity; this variation is responsible for the magnitude dispersions in the mock LF. The dependence of the RC magnitude on the age and metallicity of the population has also been recorded empirically, e.g. the globular cluster 47 Tuc (age = 11  Gyr, [Fe/H] = −0.7 dex) has |$\mu _{K}^{\mathrm{RC}}=-1.28$| mag while the much younger open cluster NGC 2204 (age = 1.7  Gyr, [Fe/H] = −0.38 dex) has |$\mu _{K}^{\mathrm{RC}}=-1.67$| (see table 1 in Percival & Salaris 2003 and top panel of fig. 6 in Girardi 2016 that shows the RC in open clusters as a function of age and metallicity). It is therefore important, for building a reliable LF, to have accurate assumptions on the metallicity distribution function and the age of the bulge population as they influence both the mean magnitude of each Gaussian (AGBB/RC/RGBB) and the magnitude dispersion.

The photometric metallicity maps derived by Gonzalez et al. (2013) measure a metallicity gradient in the bulge of ∼0.04 dex deg⁻¹. The metal-rich stars ([Fe/H] ∼ 0) dominate the inner regions with little to no metallicity gradient within |b| < 5°, while more metal-poor stars prevail at distances further from the plane. The absolute majority of stars in the VVV survey are found at |b| < 5°, so our assumption of solar metallicity (adopted by Robin et al. 2003 and after, by Wegg & Gerhard 2013) with 0.4 dex dispersion should be a good approximation. Recent measurements have also shown that the metallicity distribution function of RC stars is bimodal (e.g. Zoccali et al. 2017) with two peaks, one slightly more metal rich and one less metal rich than the Sun; however, the two mean metallicities and their ratio vary across the 26 sparse fields observed by the GIRAFFE survey (see fig. 4 in Zoccali et al. 2017), making it difficult to implement these variations in a spatially varying LF over 315 square degrees.

An age variation in the bulge stellar population is another factor that could shift the mean magnitude and increase the magnitude dispersion σ of the RC, AGBB and RGBB populations: for example, according to theoretical models, the RC of a 5 Gyr population is 0.1 mag brighter than that of a 10 Gyr population (see the violet LF in the right-hand panel of Fig. 3 and the last row of Table 1 that gives the model parameters with a skew-Gaussian RC distribution). In the literature, there is overall good consensus on the age (10 Gyr) of the bulge (see review by Zoccali & Valenti 2016 and references therein), with some studies postulating the existence of a younger bulge component (e.g. Ness et al. 2014; Catchpole et al. 2016; Haywood et al. 2016).

In the next subsection, we take into account other factors, such as the extinction residuals and the VVV photometric errors, that influence the magnitude dispersion of the RC, RGBB and AGBB (more often we just mention σ^RC, as the RC is the main population in our VVV sample), which are characteristic to the survey and our methods – rather than to the intrinsic bulge population. The LF is further discussed in Section 7.

4.3.2 VISTA photometric errors and residual extinction

The observed RC in the VVV CMD, shown in Fig. 7, will be more spread compared to the theoretical LF because of variations in distance, extinction, age and metallicity in the bulge population, coupled with survey-specific photometric errors. We calculate the total magnitude dispersion |$\sigma _{K_{s}}(l,b)$| caused by photometric errors and extinction corrections (see equation 6 in Wegg & Gerhard 2013) and show its spatial variation in Galactic coordinates in Fig. 4. To account for the effects of photometric errors and residual extinction correction, we convolve the theoretical LF (green curve in Fig. 3) with a Gaussian kernel with |$\sigma _{K_{s}}(l,b)$|⁠, for models E and S described in the next session. However, the effects of the convolution are negligible on the number of counts predicted by our model: convolving the intrinsic LF with a Gaussian with the maximum dispersion, max(⁠|$\sigma _{K_{s}}(l,b)) = 0.065$|⁠, found close to the exclusion boundary |$E(J^{\prime }-K^{\prime }_{s})$| = 0.1, causes a variation in the peak of the observed magnitude distribution N(K_s) smaller than 1.8 per cent. As the convolution slows down our fitting procedure and has negligible effects, it will not be performed for models with two components and the tests in the next sections.

4.4 The bulge density distribution

We do not attempt to build a parametric model of the bulge X-shape, which we discuss in detail in Section 7.

4.5 Other stellar populations in the inner MW

We assume that the halo population contributes a negligible fraction of RC stars to the whole VVV RC sample (see e.g. column 4 of table 3 in Zoccali et al. 2017). Our colour–magnitude selection box does not include stars with J − K_s < 0.4 mag, where we would expect to see the halo old metal-poor stars. We also assume that the ring, the spiral arms and the thin long bar (e.g. González-Fernández et al. 2014; Wegg et al. 2015) do not significantly contribute to the VVV star counts in the fields selected for the analysis [regions of low reddening, with |$E(J^{\prime }-K^{\prime }_{s}) <$| 1.0] as they are confined to the Galactic plane.

We therefore assume that in the regions probed by our analysis, the thin disc, thick disc and the bulge make up the bulk of the stellar density and the contribution of other populations is negligible.

5 THE FITTING PROCEDURE. TEST ON MOCK DATA

In this section, we outline the fitting procedure and demonstrate its robustness by testing it on a mock bulge data set generated with galaxia pre-defined parameters. These tests help assess the quality of our results when fitting the VVV data set in Section 6.

5.1 The mock data

The mock data (D^mock) is a catalogue of stars generated with galaxia (Sharma et al. 2011) using a specified set of parameters. It consists of a synthetic disc (⁠|$N_{{\rm d}}^{\mathrm{mock}}$|⁠) multiplied by a constant (Scale) and a mock sample of bulge/bar stars (⁠|$N_{{\rm B}}^{{\rm mock}}$|⁠) so that the number of stars in a given sightline at a given magnitude K_s is |$D^{\mathrm{mock}} = N_{{\rm d}}^{\mathrm{mock}}\times {\rm Scale} + N_{{\rm B}}^{\mathrm{mock}}$|⁠, where Scale = 0.4 (a value close to the observed value reported in Table 3). The parameter Scale, introduced to match the total number of disc stars, is a free parameter in our model and is needed to account for variations in the number density introduced by changing the isochrone set. We generate with galaxia a mock bulge catalogue of ∼13 million stars, to roughly match the number of VVV stars, following model G (equation 1) density distribution with parameters [x₀, y₀, z₀, |α|] = [1.6, 0.4, 0.4, 20°]. This distribution is shown in Fig. 5, projected on to the X–Y plane, where the GC is at [X, Y] = [0,0]. Only the latitude range −5° < b < 5° is plotted but we fit −10° < b < 5°, as in the data.

5.2 The model

In the fitting process, we calculate the predicted number of counts N_B(K_s) for 20 different values of K_s between 12 < K_s/mag < 14 so that the total number of counts is |$N_{\mathrm{B}}(12 < K_{s} < 14) = N_{\mathrm{B}} =\sum \limits _{j = 0}^{20} (N_{\mathrm{B}})_{j}(K_{s})$|⁠. The limits of the integral in equation (5) are set to 4 and 12 kpc as ρ_B approaches zero outside this range for 12 < K_s < 14.

5.3 The fitting procedure

5.3.1 Method

We use the python implemented limited-memory BFGS method (L-BFGS; Byrd, Nocedal & Zhu 1995) that calculates an estimate of the inverse Hessian matrix to facilitate search through the model parameter space. This method is particularly well suited for optimization problems with a large number of parameters, e.g. including the ones governing the shape of the bulge density, n, c_⊥, c_||. Moreover, it allows the user to constrain the parameters within reasonable physical values, equivalent to adopting a uniform Bayesian prior. We draw the initial (guess) values of the free parameters from uniform distributions, with the same limits as our grid search, and add three extra parameters, n, c_⊥, c_||, also drawn from random uniform distributions, with limits as in Fig. 6.

5.3.2 Optimal integration step estimation

Before fitting, we need to find the most appropriate integration steps (Δl, Δb and ΔD) for the bulge model (N_B), the only population we generate analytically and whose parameters we want to determine. Using the same density function and true values in both the synthetic bulge and the mock catalogue, we made several tests varying (Δl, Δb and ΔD) in our model to find the resolution that produces the best match to the mock catalogue. We found that small-sized fields (Δl, Δb) = (0|$_{.}^{\circ}$|1,0|$_{.}^{\circ}$|1) produced a reduced chi-squared |$\tilde{\chi }^{2}$| ≈ 1 while the distance step can be relatively large (e.g. 400 pc) given that we calculate N_B(K_s) for magnitude step sizes of dK_s = 0.1 mag. With this choice, the model gave a perfect match to the mock catalogue, with a reduced chi-squared of |$\tilde{\chi }^{2} = 1.02$|⁠. The number of sources varies rapidly from pencil beam to pencil beam in the inner Galaxy, so it is not surprising that in order to have an accurate model prediction, we need to calculate parameters over a very fine grid of Galactic coordinates.

To fit the VVV data, we increase the magnitude resolution to dK_s = 0.05 mag, for which we need to employ a finer distance integration step of ΔD = 50 pc.

5.4 Test results

In the test run, we compute chi-squared maps varying five parameters ϑ₀ = [x₀, y₀, z₀, |α|, S], over a wide range of values on a coarse grid [Δx₀, Δy₀, Δz₀, Δ|α|, ΔS] = [0.3 kpc, 0.1 kpc, 0.1 kpc, 2°, 0.2]. The high-reddening region where |$E(J^{\prime }-K^{\prime }_{s}) >1.0$| (grey mask in Fig. 4) is excluded from the test, to reproduce the data-fitting procedure. We compute N_B(12 < K_s < 14) for fields of 0|$_{.}^{\circ}$|1 × 0|$_{.}^{\circ}$|1, and in each pencil beam we calculate N_B for 20 values of magnitude K_s (equation 5). We group the model M (equation 4) in five magnitude bins to be fitted on five magnitude bins in the mock catalogue (VVV data in Section 6), in order to extract information about the depth distribution of sources in the bulge.

Fig. 6 shows the resulting slices of ln (χ²) maps as a function of the two varying parameters, while the other three are fixed to their true value. The absolute minima (violet bin in the figure) coincide with the true values used as an input of the model (yellow open circles in the figure), proving that the method works well. These maps also reveal the correlations between parameters, see e.g. α versus x₀.

The mock catalogue can also be used to test the accuracy of the L-BFGS method employed in Section 6 to find the parameters of the real bulge density distribution. The 10 filled red circles in Fig. 6 are the initial parameters of 10 L-BFGS independent runs with initial values drawn from uniform distributions. The identical final solutions of each run (red open circles) where the algorithm converges coincide with the fiducial value (yellow circles) ϑ = [1.629, 0.4026, 0.3996, 20.01, 3.88, 1.925, 1.993, 0.3865] ± [0.004, 0.0009, 0.0009, 0.02, 0.015, 0.007, 0.004, 0.0009] and are very close to the true value ϑ₀ = [1.6, 0.4, 0.4, 20, 4, 2, 2, 0.4], included in Table 2. The statistical errors in the parameters, obtained from the covariance matrix, are extremely small, and we consider them to be negligible in comparison to the systematic errors (see Section 7). However, the test on the mock data also provides an indication of how well the fitting method recovers the true parameters of the bulge density distribution (see the first two lines in Table 2, which shows the true values and the best-fitting values found using the L-BFGS minimization method). The agreement here leads us to conclude that the random error component is well captured by the method.

An estimate of the systematic errors can be obtained by adopting different density distributions in the models as compared to the mock data. In Table 2, we list the results obtained from fitting a hyperbolic secant model and an exponential form of the density distribution, with r the Galactocentric radius (last two rows of the table) and from comparing the results of the three fits to the true values ϑ₀, we notice that the best-fitting parameters ϑ of the correct density law are close to their true value ϑ₀(δϑ = [x, y, z] − [x₀, y₀, z₀] = [0.029, 0.003, 0.000]) while when fitting the ‘wrong’ density laws, the results diverge significantly from the true value (δϑ = [0.234, 0.056, 0.029] and δϑ = [0.361, 0.068, 0.054] for the two ‘wrong’ models). In conclusion, the choice of the density distribution alone introduces systematic errors much larger than the standard errors.

6 FITTING THE VVV DATA

In this section, we fit a triaxial density model to the observed VVV magnitude distribution in the RC region of the CMD. We focus on describing the inner MW large-scale structure by finding a model capable of reproducing the distributions seen in the majority of the fields at the expense of discrepancy in a small number of fields, which include fields close to the Galactic plane/Centre and those at b < −5°, where a double RC is observed (e.g. McWilliam & Zoccali 2010; Saito et al. 2012).

6.1 Data sample

The photometry of every source in the single band tile catalogues produced by CASU is labelled with a morphological classification flag. In the following sections, we use detections classified as ‘stellar’, ‘borderline stellar’ or ‘saturated stars’ in the K_s magnitude band, which has the best seeing. We avoid applying a strict cleaning procedure in order to work with a sample of bulge stars as complete as possible.

We choose sources with photometric errors σ_J, |$\sigma _{K_{s}} < 0.2$| mag and colours (J − K_s) > 0.4, to minimize contamination from disc main-sequence stars and (J − K_s) < 1.0 to separate red giant stars from spurious objects and highly reddened K and M dwarfs. In our analysis, we apply a faint-end magnitude cut of K_s = 14 mag, which we consider adequate to minimize the effects of incompleteness given the substantial extinction over the field of view, and a bright-end limit of K_s = 12 mag to avoid saturation and non-linearity (Gonzalez et al. 2013), which sets in around |$K^{\prime }_{s} < 11.5$|⁠. These cuts form our colour–magnitude selection for the VVV working sample and are marked with a red box in Fig. 7. The two panels show the distribution of all the VVV sources with |$E(J^{\prime }-K^{\prime }_{s}) < 1$|⁠, totalling 96 million objects, in a CMD before (left-hand panel) and after (right-hand panel) extinction correction. The observed |$K_{s}^{^{\prime }}$|-band magnitudes for each star are corrected for the effects of extinction, |$K_{s} = K_{s}^{^{\prime }} - A_{K_{s}} (l, b)$|⁠, using the reddening maps |$E(J^{\prime }-K^{\prime }_{s})$| constructed from the VVV data (see Section 3). In the left-hand panel, most of the RC stars are outside our selection box due to extinction and reddening effects that cause the stars to appear fainter and cooler. The number density distribution of all the stars within the red selection box is then shown in Galactic coordinates in Fig. 8. Again, the contours delineate the highly reddened regions |$E(J^{\prime }-K^{\prime }_{s}) = 1.5$| (light grey) and |$E(J^{\prime }-K^{\prime }_{s}) = 1.0$| (thick grey).

We illustrate the completeness of the survey by plotting the variation of the observed magnitude distributions with latitude in Fig. A1 of Appendix A, both before (top row) and after (bottom row) extinction correction. The stars were selected in 1|$_{.}^{\circ}$|0 × 0|$_{.}^{\circ}$|4 fields at constant longitude l = 0° and decreasing latitude b. The fields within the |$E(J^{\prime }-K^{\prime }_{s})=1$| mag contour are incomplete with the high extinction causing a depletion of stars in the 13–15 |$K^{\prime }_{s}$| magnitude range and a rapid drop in the star counts for |$K^{\prime }_{s} > 15$| mag (see the orange dotted line in the figure). However, further from the plane, the incompleteness limit increases and at b ∼ −3° the distribution is complete up to |$K^{\prime }_{s} \sim 15.2$| or K_s ∼ 15, making our working magnitude range (red vertical lines in the figure) a fairly conservative cut.

6.2 An example data fit

The functions involved in the data-fitting procedure are illustrated in Fig. 9. The apparent magnitude distribution N_B of the bulge (thick blue, equation 5) is the integral of the absolute magnitude LF (green, labelled ‘LF+14.5’) with a distance modulus of 14.5 added to the absolute magnitudes, |$M_{K_{s}}$|⁠, times the density law ρ_B(l, b, D) (in the figure we show ρ_B(l, b, μ) in dotted line, equation 2). The total apparent magnitude distribution M = 〈 N_d〉 × Scale + N_B (green triangles, labelled ‘MODEL E’) is given by the sum of the bulge N_B (labelled BULGE E) and the discs 〈 N_d〉 (labelled ‘DISCS’). To obtain N_B (dark green continuous line, labelled BULGE E), equation (5) was calculated on a grid of 20 × 10 (l, b)_i fields of view between 0° < l < 2° and 3° < b < 4° for 50 K_s apparent magnitude bins between 12 < K_s < 14 mag.

We do not fit directly for the Sun–GC distance, which is a fixed parameter, R_⊙ = 8 kpc (see table 1 in Vallée 2017), but we allow the mean magnitude of the RC to vary by Δμ^RC. Here Δμ^RC can translate to a new value for R_⊙, as it is effectively a distance modulus shift, or a new value for the mean magnitude of the RC. For all models, we found Δμ^RC = −0.1 or |$\mu _{10}^{\mathrm{RC}} = -1.63$| at R_⊙ = 8 kpc, similar to Laney, Joner & Pietrzyński (2012) and Hawkins et al. (2017) who found |$M_{K}^{\mathrm{RC}} = -1.61$| using nearby Hipparcos and Gaia stars, respectively, and Wegg et al. (2015) who found |$M_{K}^{\mathrm{RC}} = -1.67$| at R_⊙ = 8 kpc (note that their LFs were constructed interpolating the BASTI isochrones in the 2MASS photometric system; Pietrinferni et al. 2004).

The model plotted in Fig. 9 is a result of the full global fit of the whole VVV data set, discussed in the next subsection. We performed the fits for two cut-off radii, R_c = 2.5 and 4.5 kpc, and found that R_c = 4.5 kpc provides a better fit; therefore, in the remainder of this work, it is kept constant.

6.3 The full bulge sample fit

We compute the model on a continuous grid of 6 arcmin × 6 arcmin fields covering the whole area of the VVV survey amounting to 200 × 150 (l, b)_i lines of sight, for 40 magnitude bins K_s (green triangles in Fig. 9) to be fitted to the data (red circles in the same figure). We first consider a single bulge population, and next we add a second component to improve the agreement in the central regions at low latitudes where the residuals highlight an overdensity.

The three models we consider are as follows:

• model E: the density model is described by an exponential-type model (equation 2) and has a 10 Gyr old input LF (see Fig. 3). The model M has 10 free parameters, seven describing the density law, a scaling for the thin discs and one for the thick disc and the RC peak magnitude shift Δμ^RC;

• model S: this density law is described by a hyperbolic secant density distribution (equation 3) and has a 10 Gyr old input LF, the same as model E. The model has nine free parameters, six describing the density law, a scaling factor for the thin discs and one for the thick disc and Δμ^RC;

• model S + E: Robin et al. (2012) using 2MASS data found that a combination of an exponential-type (E) and a sech² (S) density law provides the best description of the bulge (see their table 2), and we use the same description for the VVV data. In total, the model has 16 free parameters, seven for model E, six for model S, two scaling factors for the discs and one for the ratio between the two bulge components. The magnitude shift is not a free parameter; we set to |$\Delta _{\mu }^{\mathrm{RC}} = -0.1$| in accordance with the results found in the single-component fit.

We perform several runs of the L-BFGS minimization starting from random initial points to test whether the final result is dependent on the initialization position and to avoid being trapped in a local minimum. We consider our best solution to be the one that maximizes the log-likelihood (equation 6) across all trials; note that several runs may actually converge to the same best-fitting solution. We do not have tight constraints on the parameters, limiting them to an initial set of random values that would have a physical meaning as described in the previous section. The model (S + E) has 16 free parameters so to avoid getting not meaningful results we choose the starting point for the main component S to be the best-fitting solution of the one population fit (model S), while the starting point for the E component is drawn from a uniform random distribution.

6.4 Final results and MCMC checks

The best-fitting parameters describing the bulge density laws are listed in Table 3 for each model considered together with the maximum log-likelihood ln(L_P) and the total number of stars in each component. The statistical errors for all model parameters are negligible compared to the complex systematic errors; thus, they are not reported.

In Section 5, we tested our fitting procedure on a mock catalogue and proved that it successfully recovered the input values. We thus expect, within the limitations of the model and of the LF assumptions, to have found the true values of the bulge density distribution. None the less, we execute an additional check of the results and of their statistical errors with a Markov chain Monte Carlo (MCMC) in a Bayesian framework. In Fig. B2 in Appendix B, we compare the results (green dots) obtained using the L-BFGS minimization method and an MCMC, for model S (label ‘σ^RC free’ in Table 3, see the discussion in Section 7.1).

We have used the affine-invariant MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013) to explore the full posterior distribution with 30 walkers. The contours of the two-dimensional projection of the posterior probability distribution of the parameters of model S are nearly Gaussian (Fig. B2), but there is a clear correlation between the scalelengths of the bulge density law, and σ^RC and the viewing angle α. The marginalized distribution for each parameter is shown along the diagonal, and the 0.5 quantiles with the upper and lower errors (16 per cent and 84 per cent quantiles) are given above each 1D histogram. For the initial state of the Markov chain, we chose the best-fitting results reported in Table 3 for model S, label ‘σ^RC free’: with 30 000 iterations (20 per cent of the chain was discarded as burn-in) and an acceptance fraction of 43 per cent, the final results and errors on the parameters are in complete agreement with the values found with the L-BFGS method. The statistical errors are equally small and, as can be seen from Fig. B1 showing a sample of 10 chains for 8 parameters, they converge to a value close to the initial one (the green line in Figs B1 and B2). We note that Wegg et al. (2015) also report small statistical errors: ‘We have performed an MCMC to estimate errors on both the model parameters (...) The resultant statistical errors are extremely small, significantly smaller than the difference between fits with different models. We therefore consider the statistical errors to be negligible in comparison to systematic errors’.

There may exist various causes for the systematic errors: the extinction-related dispersion and photometric errors (see Fig. 4); confusion of point sources in overcrowded regions; classification and calibration errors; and the fact that the model is not an exact description of the data. In addition, at the end of Section 5, we have shown that the choice of the bulge density distribution alone can introduce systematic errors larger than the standard errors.

The next subsections look into more detail at the models reported in Table 3 and their suitability in describing the data.

6.5 Data versus model comparison

We here discuss Figs 10–15, which show the data and the best-fitting models (labelled ‘main’ in Table 3) in different projections.

The log-likelihood values in Table 3 suggest that model E performs slightly better than model S; however, there are many similarities between the two sets of results: for both, we obtain α ≃ −19|$_{.}^{\circ}$|5 and axial ratios [1:0.44:0.31]. As we will see in the residual maps, a two-component model (S + E) provides a better description of the data especially in the inner regions, but at the expense of a larger number of free parameters.

Using the first data release of the VVV, Wegg & Gerhard (2013) fitted exponentials to the density profiles of the major, intermediate and minor axes shown (see their fig. 15) and found scalelengths of 0.70, 0.44 and 0.18 kpc, respectively, corresponding to axial ratios [1:0.63:0.25]. In the discussions section, we compare our results with other literature estimates for the shape and viewing angle.

In Fig. 10, we show the number density distribution of VVV stars in five magnitude bins (left column) and the best-fitting models E and (S + E) in Galactic coordinates (middle and right column, respectively). Model S is very similar to model E and thus not shown. A magnitude slice is almost equivalent to selecting stars within a certain distance range, where the brightest magnitude slice 12 < K_s < 12.4 selects the nearest stars, in front of the bulge, and the faintest magnitude slice 13.6 < K_s < 14.0 selects the furthest stars, behind the bulge.

The middle and right-hand panels of Fig. 11 show the difference between the data and the E and S + E models, divided by the Poisson noise in the model, (N_i−M_i)/|$\sqrt{M_{i}}$|⁠, for each magnitude bin. The addition of the second component E (left column) improves the fit in the central region, especially in the bright magnitude slice (top-right panels of Fig. 11). Model S + E assumes a young population (5 Gyr), described by model E and an old one (10 Gyr), described by model S, which is the main component. The younger component has a 0.1 mag brighter peak than the 10 Gyr old population (see Fig. 3 where |$\mu _{5}^{\mathrm{RC}}=$| −1.64, hence |$\mu _{5}^{\mathrm{RC}}=$| −1.74 after applying the Δμ^RC = −0.1 shift), which helps to provide a better fit in the central regions of the brighter magnitude bins. The LFs for both populations are shown in the right-hand panel of Fig. 3, where the integrated area under each curve gives the total number of stars expected for each component for the best-fitting parameters. We provide several interpretations for the additional component in Section 7.

To summarize our results, we show the residuals for all three models in one magnitude slice in Fig. 12. As previously observed, models E and S do not provide a good match for the data in the area close to the GC, |b| < 3° and |l| < 2°, though in general they perform very well. The model (S + E) matches the observed distribution remarkably well (with a median percentage residual of 5 per cent) apart from some discrepancies within the central tiles b306, b320, b347. It is not clear at this stage whether this is entirely due to data processing or calibration issues or possibly some real feature present in the central data.

The best-fitting density laws projected along the Z-axis (bottom panels) and X-axis (top panels) for the models E, S and S + E (first three columns of the figure) are shown in Fig. 13. In the figure, we mark the cut-off radius R_c, the viewing angle and the limiting longitudes of the VVV survey. Overlaid on model E, we add a density law that is a complex parametrization of the bulge X-shaped component, provided by López-Corredoira (2017) in their equation 4. This model predicts a double RC at intermediate latitudes. We only show a slice of the density law within −0.7 < Z/kpc < −2.0 (below the red line in the Y–Z plane) corresponding to b < −5°, projected on to the X–Y plane. The density contours of the foreground RC (X < 0) are marked in violet and those of the background RC (X > 0) in yellow. The same colour scheme is adopted for the Y–Z plot projection where we do not make any Z cuts, to show that the nearer RC will be dominant at positive longitudes while the further RC will be dominant at negative longitudes, symmetrically above and below the plane. In addition, we can expect (1) the contribution to the magnitude distribution of the background RC to diminish and the longitude to increase towards positive values, as the contribution of the foreground RC increases and (2) the separation between the two RCs in the magnitude distribution to increase as we look away from l = 0°. This separation will depend on the angles formed by the X-shape with the Galactic plane as well as on the Galactic coordinates or distance from the plane.

The Y–Z distributions display a ‘boxy’ bulge with a flattened distribution along the Z direction. The difference between the E and S models (fourth column) illustrates that the E density law is more centrally concentrated. The E − (S + E) residuals in the X–Y plane (fifth column) are irregular because the two components in the S + E model have different viewing angles.

In Fig. 14, we show eight examples of magnitude distributions in 2° × 1° fields where the data points we fit are the red circles and the best-fitting models are represented with a black dotted line (model S + E), green triangles (model E) and a thin blue line (model S). The disc populations are shown in grey (the updated discs, discussed in Section 7.2, in pink). The centres of the fields are marked with letters in the middle column of Fig. 11, where each letter represents a column of Fig. 14. The panels are organized as follows.

• Column (a) shows the magnitude distributions in two symmetric fields centred at (l, b) = (1°, 3|$_{.}^{\circ}$|5) and (1°, −3|$_{.}^{\circ}$|5): the star counts are symmetric with respect to the Galactic plane, and there is good agreement between the data and the best-fitting model.

• Column (b) shows the magnitude distributions in two tiles (b306 and b347) that have pronounced residuals in Fig. 12 with centres at (l, b) = (−1°, 1|$_{.}^{\circ}$|5) and (1°, −1|$_{.}^{\circ}$|5), though interestingly these are symmetric with respect to the GC. The fit is poorest in the first two brighter magnitude slices, which produce the overdensities observed in the top panels of Fig. 11. As we have shown, adding an extra component partially solves the problem: model S + E gives the best fit to the data in central regions, making the overall fit for the bright magnitude slices better, as can be seen in the top-right panels of Fig. 12.

• Column (c) shows the magnitude distributions in symmetric fields with respect to l = 0°, centred at (l, b) = (7°, 2|$_{.}^{\circ}$|5) and (−7°, 2|$_{.}^{\circ}$|5). The peak of the RC is shifted between positive (bright RC) and negative (faint RC) longitudes, a consequence of the orientation of the bulge with regard to the Sun–GC direction.

• Column (d) shows magnitude distributions in two fields at intermediate latitudes, centred at (l, b) = (1°, −6|$_{.}^{\circ}$|5) and (−1°, −6|$_{.}^{\circ}$|5). In these fields, the double RC is clearly visible in the VVV data as already shown by previous studies (e.g. Saito et al. 2012). Because the fields are symmetrical to l = 0°, the two RCs contribute similarly to the magnitude distribution, even though the foreground RC will contribute slightly more to the distribution at l = 1° compared to l = −1°; for the background RC, the opposite is true. The X-bulge density law proposed by López-Corredoira (2017) describes well this phenomenon: by looking at Fig. 13, we can clearly see how the l = 1° line of sight will pass closer to foreground RC than to the background RC, but thanks to the viewing angle under which we observe the bar, both RCs will be visible in the magnitude distribution, even though in different contributions. The double-peaked magnitude distribution shown in yellow is the López-Corredoira (2017) model calculated using equation (4) where ρ_B is the X-shaped density distribution that we have shown in Fig. 13. The model is successful at predicting the double RC, as expected.

We now explore the signal of the double RC at intermediate longitudes l = ±5° and varying distance below the GP by comparing the data with our models and the model based on the X-shaped density law (which should be accurate for b < −5°) in Fig. 15. The centres of the 2° × 1° fields are marked with ‘x’ in Fig. 11. The top and bottom rows show the apparent magnitude variations at constant longitudes, l = 5° and −5°, respectively, with the latitude decreasing from left to right between −1|$_{.}^{\circ}$|5 and −7|$_{.}^{\circ}$|5.

The fields close to the Galactic plane (first column, −2° < b < −1°) are largely incomplete because we mask the high-extinction regions; thus, the amplitude of the RC is not well fitted (please note that this is the result of the global fit). At −4° < b < −3° (second column), the models perform well but at lower latitudes −5° < b < −6° (third column), where the double RC becomes visible, we can clearly see that the amplitude of the RC peak in the data is not well matched by the model. The excess of stars supposedly belongs to the peanut shape that is not well represented by a single simple triaxial shape. Even further from the Galactic plane, where the number counts are much lower (fourth column, −7° < b < −7°), the discrepancy is more obvious at all magnitudes.

At b < −5°, the model using the X-shaped density law (López-Corredoira 2017) is doing a better job at predicting the magnitude distribution. However, while the model predicts a double RC, the data only show one. This could be explained if the second RC contributed less than predicted by the López-Corredoira (2017) model at intermediate longitudes (this could easily be the case, e.g. imagine an l = 5° line crossing the X–Y plane in Fig. 13). From the magnitude variation of the peaks at constant longitude (left to right), it is clear that with increasing distance from the plane the foreground RC gets brighter (positive longitude, top row) and the background RC gets fainter (negative longitude, bottom row), behaviour specific to an X-shaped distribution. To roughly quantify what proportion of the VVV stars are found in the X-shaped residuals relative to the total number of bulge stars, we consider the fields at intermediate latitudes, outside |b| < 4°. Here the feature is most clearly visible in the VVV data, even though it may continue to lower latitudes. The model tends to overfit the region between the X-shaped arms (see blue depression at l = 0° in Fig. 12) and underfit the X-shaped arms themselves; hence, we conjecture that the model has to be reduced by some fraction <1 such that only ∼5 per cent of the residuals are negative, and the positive excess that remains would theoretically belong to the arms. With this assumption and mirroring the region at b < −5° to positive latitudes, b > 5°, we find that ∼7 per cent of the stars belong to the X-shaped arms relative to the entire modelled bulge, within −10° < l < 10°, −10° < b < 10°.

7 DISCUSSION

7.1 The luminosity function

7.1.1 Variations of the LF

The intrinsic LF of the bulge is sensitive to both age and metallicity variations but in the absence of accurate and consistent measurements of these quantities across the bulge, we do not attempt to build a spatially varying LF.

The main effect of not having a simple bulge population as the one we assumed in Section 4.3 will be to blur the intrinsic LF. We therefore convolve the galaxia generated LF with a Gaussian kernel of σ as a free parameter (therefore σ_free^RC|$=\sqrt{\sigma ^{2}+(\sigma ^{\mathrm{RC}})^{2}}$| becomes a free parameter, where σ^RC = 0.06 mag) and re-fit the VVV data with models E and S. The best-fitting results are listed in Table 3 for both models and are labelled ‘σ^RC free’.

Compared to the approach presented in this paper, Wegg & Gerhard (2013) followed a different route and attempted to constrain the bulge properties non-parametrically. Instead of estimating the parameters of a bulge density law, they deconvolve the observed RC magnitude distributions to extract line-of-sight RC densities from which they build a 3D map assuming eightfold mirror symmetry. The result they obtained, α = ( − 26|$_{.}^{\circ}$|5 ± 2°), falls in the middle of our extreme values for σ^RC = 0.06 mag, α = −19|$_{.}^{\circ}$|4 (red dot) and σ_free^RC = 0.26 mag, α = −34|$_{.}^{\circ}$|4 (green dot).

In the absence of a fully reliable calibration of the bulge RC dispersion in the literature (the latest measurement in the local neighbourhood with Gaia finds σ^RC = 0.17 mag – see Hawkins et al. 2017), we investigated the effect of the σ^RC variation on the density law parameters by convolving the galaxia generated LF with a Gaussian of increasing σ in order to obtain σ^RC = 0.09, 0.12, 0.15, 0.18 and 0.21 mag for a 10 Gyr population. The results shown in Fig. 16 and listed in Table 4 confirm the expectation that a higher σ^RC is responsible for a higher viewing angle α. This was also observed by Wegg & Gerhard (2013) who found that |α| reduces from 26|$_{.}^{\circ}$|5 to 25|$_{.}^{\circ}$|5 (triangles in Fig. 16) when σ^RC decreases from 0.18 to 0.15 mag and explain ‘this is because increasing the line-of-sight geometric dispersion has given the illusion that the bar is closer to end on’. We find a stronger dependence of the viewing angle on σ^RC, with α decreasing from 24|$_{.}^{\circ}$|69 to 22|$_{.}^{\circ}$|64 for the same σ^RC variation. Our work is probably more sensitive to the viewing angle as we do not need to assume eightfold mirror symmetry, thanks to continuous coverage in all quadrants of the newer VVV data release.

While the viewing angle is strongly dependent on σ^RC, the axial ratios do not vary significantly as seen in the middle and right-hand panels of Fig. 16. In fact, the parameters controlling the shape of the density distribution are not significantly affected by the choice of σ^RC even though a trend can be observed in the middle and right-hand panels of the figure: in the middle panel, as the σ^RC increases (the values are colour-coded as in the left-hand panel), the bar gets puffier in both the vertical and Y-axis directions but as it reaches σ^RC = 0.15 mag it starts to get thinner again along the Y-axis while it continues to get thicker in the vertical directions. The parameters controlling the shape of the distribution are also correlated, with n decreasing as σ^RC increases while the c_∥/c_⊥ ratio also increases. Notice that these variations are small and do not change the shape of the distribution significantly, the viewing angle α being the most affected parameter.

7.2 The discs

7.2.1 The thin disc

In galaxia, the thin disc density distribution has a scalelength of R_d = 2.53 kpc and a disc hole scalelength of R_h = 1.32 kpc. More recently, Robin et al. (2012) used 2MASS photometry to fit simultaneously for the parameters of the bulge population and the shape of the inner disc. The updated parameters of the thin disc in the best bulge model (model S + E, see their table 2) are R_d = 2.17 kpc and R_h = 1.33 kpc: the disc hole scalelength is almost identical and the disc scalelength is shorter by 400 pc compared to the previous values.

7.2.2 The thick disc

Robin et al. (2014) used photometric data at high and intermediate latitudes from Sloan Digital Sky Survey (SDSS) and 2MASS surveys to re-analyse the shape of the thick disc. Their updates include a new IMF best fit, dN/dm ∝ m^−0.22, a slightly smaller scaleheight of h_Z = (535.2 ± 4.6) pc and a shorter scalelength of h_R = ( 2362 ± 25) pc plus a new age best fit of 12 Gyr.

We have performed all the fits described in the paper with the updated thin and thick disc versions, and the results are reported in the last rows of Table 3. They are also shown in pink in Figs 14 and 15. Finally, we have decided to use the default (Robin et al. 2003) values throughout the paper, for the following reasons.

• In the fitting procedure, we exclude the area close to the Galactic plane where the reddening |$E(J^{\prime }-K^{\prime }_{s})$| is higher than 1 mag and where the majority of the thin disc population lies. The |$E(J^{\prime }-K^{\prime }_{s}) = 1.0$| contour extends, on average, to latitudes −1|$_{.}^{\circ}$|5 < b < 1|$_{.}^{\circ}$|5 (see Figs 1 and 2), which is equivalent to a height perpendicular to the Galactic plane of −0.2 < Z/kpc < 0.2 at 8 kpc. This region is masked out in the fitting procedure described in Section 6. The scaleheight of the oldest thin disc component is ∼300 pc; therefore, a large part of the thin disc population is excluded from the fit.

• The fitting results confirmed that the bulge density law parameters remain almost identical, regardless of the discs (original or updated version) assumed.

• To avoid multiple modifications to the galaxia code.

7.3 Stellar mass of the bulge

galaxia relies on the isochrone table as input to the code to assign a mass to each star produced as part of the mock catalogue for a given IMF. Pietrinferni et al. (2004) undertook a detailed analysis of the outer bulge stellar density and LF, assuming a Salpeter (1955) IMF, by fitting model parameters to a set of 94 windows in the outer bulge situated at −8° < l < 10° and |b| < 4° in the Deep Near-Infrared Survey (DENIS). The same IMF was used in the Besançon model and implemented in galaxia. The bulge stars lie at distances between 4 and 12 kpc, corresponding to absolute magnitudes M_K in the range −3.5 to 1 mag in the 12–14 apparent magnitude range K_s. At these bright magnitudes, the intrinsic LF is not influenced by the choice of IMF; therefore, using a different IMF would not change the results in Table 3.

The stellar mass in the bulge is however sensitive to the IMF choice. We calculate the mass of the bulge assuming a bottom-light Chabrier IMF, similar to a Kroupa IMF, which we have used throughout the paper. The slope of the latter becomes significantly shallower than the Salpeter IMF below 1 M_⊙ and describes more accurately observations in the Galactic bulge (see also Section 4.3). The majority of stars in the bulge have masses <1 M_⊙ so we would expect to greatly overestimate the total mass of the bulge using a Salpeter IMF. In addition, Sharma et al. (2016) test the theoretical predictions of galaxia with asteroseismic information for 13 000 red giant stars observed by Kepler and find that the distribution of the masses predicted by the galaxia Galactic model overestimates the number of low-mass stars. We therefore report the bulge mass estimates using the Chabrier IMF that is similar to the Kroupa IMF, excluding the Salpeter IMF.

Using the best-fitting model (E) in galaxia, we can predict the fraction of bulge stellar mass in the RGB and the RGB luminosity that we should observe within the VVV coverage, with the colour–magnitude cuts 0.4 < J − K_s < 1.0 and 12 < K_s < 14, and hence use the observed population (and the measured fraction) to constrain the total stellar mass of the bulge. For a given IMF, the PARSEC isochrones enable us to calculate the remnant population and therefore the total mass of the original population. For the total MW bulge mass, we obtain M_bulge^Chabrier = 2.36 × 10¹⁰ M_⊙, including the mass of the remnants (49 per cent of the total mass of the bulge) using a Chabrier IMF. For the two-component S + E best fit, we obtain 2.21 × 10¹⁰ M_⊙ (S component) and 1.2 × 10⁹ M_⊙ (E component), for Chabrier IMFs. If the low-latitude residuals are caused by the presence of an extra component (component E), then this component would contribute ∼5 per cent to the total bulge mass budget.

Other recent measurements of the bulge stellar mass include the following.

• Robin et al. (2012) who estimate the total mass of the bulge by fitting a population synthesis model to 2MASS and for the two-component best fit, find M_bulge = 6.36 × 10⁹ M_⊙ where the mass of the S component is 6.1 × 10⁹ M_⊙ and the mass of the smaller component E is 2.6 × 10⁸ M_⊙.

• Portail et al. (2015) construct dynamical models of the MW bulge using RC giants and kinematic data from the Bulge Radial Velocity Assay survey. They compute the total mass of the Galactic box/peanut bulge and find a stellar mass M_bulge = 1.25–1.6 × 10¹⁰ M_⊙, depending on the total amount of dark matter in the bulge. Wegg et al. (2015) use a Kroupa IMF at 10 Gyr and report a bulge mass estimate of 1.81 × 10¹⁰ M_⊙ that includes, in addition to a bulge population, a thin and a superthin bar.

• Valenti et al. (2016) scale the total magnitude distribution obtained from all VVV fields (320 square degrees) to the one from Zoccali et al. (2003, an 8 arcmin × 8 arcmin field at (l, b) = (0°, −6°)), to obtain a stellar mass of M_bulge = 2 × 10¹⁰ M_⊙. They find that a ±0.3 dex variation in their assumed IMF slope α = −1.33 ± 0.07 (Zoccali et al. 2000) would change the mass estimate by less than 12 per cent; moreover, the authors assume constant disc contamination. In contrast, we included in our model the disc contamination that is calculated as a function of field and magnitude.

7.4 Structural comparison with previous results

Fig. 17 compares the bulge parameters measured in this work to those obtained in other studies with a variety of tracers, e.g. RR Lyrae, Mira variables and RC stars. Note that the results of modelling the Galactic bulge with an exponential and sech² density laws as reported in Table 3 are indistinguishable and are marked in red in Fig. 17. The red error bars give the maximum variation seen in the parameters as a result of varying σ^RC: only the viewing angle seems to be significantly affected (see Fig. 16).

We take advantage of a large catalogue of bulge studies up to 2009 in Vanhollebeke et al. (2009) to show the range of the viewing angle values and the axial ratios with black dots, adding the error bars when present. Additionally, four more recent results are marked with colours: one from a parametric model describing the 2MASS data (cyan; Robin et al. 2012), one from a non-parametric study of VVV RC stars (magenta; Wegg & Gerhard 2013), and two from the OGLE survey, with RR Lyrae (yellow; Pietrukowicz et al. 2015) and RC stars (blue; Cao et al. 2013). Robin et al. (2012) obtain an orientation of about 13° for the main structure, in their two-component fit. Cao et al. (2013), who use number counts of RC giants from the OGLE III survey to fit triaxial parametric models for the bulge, obtain axial ratios within the parameter space spanned by our results when testing different σ^RC; we mark their result in both Figs 16 and 17.

7.5 Central region

As in previous studies (e.g. Robin et al. 2012), we find that one-component bulge models do not provide a good match to the data in the central Galactic regions. We notice a clear overdensity at |b| < 3° and |l| < 2° (left and middle panels in Fig. 12). The match between the model and the data can be improved by considering a combination of two populations (right-hand panel in the same figure), where the central one has an RC peak brighter than the RC peak of the main population.

The presence of the overdensity can be explained using several theories.

• The presence of a ‘discy’ pseudo-bulge – in addition to a ‘boxy’ bulge (Dwek et al. 1995), our Galaxy might host a ‘discy’ pseudo-bulge, observed in many early- and late-type b/p bulge galaxies (Bureau et al 2006). These types of bulges grow ‘secularly’ out of disc gas that was transported inwards along the bars (e.g. Kormendy & Kennicutt 2004). The centrally concentrated residuals with short vertical scaleheight we observe might be indicative of this type of disc-like, high-density pseudo-bulge, also seen in N-body models by Gerhard & Martinez-Valpuesta (2012). In addition, a change of slope at |l| ∼ 4° can be seen in the VVV RC star count maxima longitude profiles close to the Galactic plane (see Gerhard & Martinez-Valpuesta 2012) while further from the Galactic plane, the change in slope disappears. Gerhard & Martinez-Valpuesta (2012) argue that a pseudo-bulge does not need the presence of a secondary nuclear bar to explain the change of slope in the longitude profile.

• Spiral arms/thin bar overdensity – the bar and buckling instabilities in the stellar disc lead to a boxy bulge that extends to a longer in-plane bar (Athanassoula 2005). The bar couples with the spiral arms in the disc, giving rise alternatively to leading, straight or trailing bar ends (Martinez-Valpuesta, Shlosman & Heller 2006). It is possible that a density enhancement produced by spiral arms along the Sun–GC line of sight might explain the residuals observed in Fig. 11 (the spiral arms are confined to the Galactic plane but a height of 100 pc would subtend an angle of ∼1|$_{.}^{\circ}$|5 at 4 kpc, while the same vertical height would only subtend ∼0|$_{.}^{\circ}$|7 at 8 kpc).

• The discs – we have used a 2D reddening map to correct for the effects of extinction, assuming that all the sources are behind the dust. Disc objects (e.g. nearby dwarfs) or bulge stars that are not behind all the dust (and are situated between the Sun and the GC) will be overcorrected by various amounts, which could possibly lead to an apparent overdensity in the central regions of the bulge in the bright magnitude slices that contain the sources closest to the Sun. We have tested this hypothesis by adding 3D extinction to a sample of galaxia generated mock thin disc stars. We have then corrected for extinction using a 2D map, therefore considering only the highest value of the 3D map. If this method was accurate, the number density of the 2D extinction-corrected sample should coincide with that of the initial mock sample. Instead, we found only a 4 per cent increase in the number of counts in the 12.0 < K_s < 12.8 magnitude range after correction for the stars near the plane but for stars with |b| > 1|$_{.}^{\circ}$|5 no increase was found. We thus believe that using a 2D reddening map does not produce the central overdensity observed.

• Younger central population – altering the disc parameters does not improve the agreement between data and model in the magnitude range considered (12 < K_s < 14). We can only account for it by adding an extra component with a smaller size ([ x₀, y₀, z₀] = [1.52, 0.24, 0.27] kpc) and a brighter RC, with ΔM = 0.1 between the main component and the secondary one. A younger population, of 5 Gyr and solar metallicity (as the main population), would produce such a shift but it could also be achieved considering other combinations of ages and metallicities between the two populations. The presence of a young population inside the bulge is not inconsistent with simulations (e.g. Ness et al. 2014) and Mira variable (e.g. Catchpole et al. 2016) observations.

7.6 The X-shape of the bulge: testing the model on a peanut bulge N-body simulation

In this work, we have not implemented a parametric description of the X-shaped component that is apparent in the RC population (e.g. Ness & Lang 2016). To test if our model is capable of robustly recovering basic parameters such as the orientation angle in realistic mock bulges, we test its performance on an N-body model (Shen et al. 2010) that has been shown to predict the split RC (Li & Shen 2012). Li & Shen (2012) found that 7 per cent of the light in the whole boxy bulge region resided in the X-shaped structure. To use the same observables (the number density distribution in different K_s magnitude ranges) as in the VVV, we convolve the N-body model with the bulge LF, ϕ_B.

The N-body model has ∼1 million particles but within the VVV footprint only ∼400 000 remain. Because of the low number of particles, we have performed two fits, one using 6 arcmin × 6 arcmin bins as in the data and one with double the bin size 12 arcmin × 12 arcmin. The simulation, the best-fitting model and the residuals are shown in Fig. 18, for the larger bin size needed to satisfy our assumption that the measurements are distributed according to a Gaussian function. We also mask the high-reddening regions to mimic the data-fitting procedure, but exclude the discs as they do not contribute significantly to the simulation. The model therefore contains only the analytic form of the bulge. Using the larger bins 12 arcmin × 12 arcmin, we find a viewing angle of 15°, on the lower end of the expected value of |$20^{\circ } _{-5^{\circ }}{^{+5^{\circ }}}$| for the simulation while with the smaller bins we find α = 20°. Nataf et al. (2015) also test their model on this N-body simulation, assuming four viewing angles α = 0°, 15°, 30°, 45° and find that α = 30°, 45° provide the best fit, with α = 15° being nearly as consistent.

The fit in the central regions of the brightest slice 12.0 < K_s < 12.4 is good, suggesting that the overdensity observed in the data does not belong to a structure specific to the peanut bulge that we could not fit due to our simple assumption of a triaxial bulge. In the second slice, 12.4 < K_s < 12.8, there is an overdensity in the simulation running diagonally from the centre towards the left, both at positive and negative latitudes; the same feature can be observed in the data residuals (see Fig. 11), but mainly at negative latitudes and positive longitudes (see also the residuals for the WISE survey, in fig. 3 in Ness & Lang 2016). Due to the orientation of the bulge, the arms of the X-shape appear asymmetric, with the arms at positive latitudes more visible in the residuals while the arms further from the Sun, behind the GC, almost invisible. This is only a projection effect caused by the orientation of the bulge as discussed in the previous section. At fainter magnitudes, the model overpredicts the number of counts at intermediate longitudes on each side of the GC; however, in the faintest two slices, the underdensity is observed mainly at negative latitudes, −2° < l < −6° and −3° < b < 3°, in the same region as in the data but slightly closer to the GC. The overdensity in the model is created by the redistribution of the total number of stars present in the peanut bulge into a boxy shape that the model is constrained to have.

Other tracers such as F dwarfs (<5 Gyr old; López-Corredoira 2016), RR Lyrae (>10 Gyr old; e.g. Gran et al. 2016; Kunder et al. 2016) and short-period Mira variables (old stars; see Catchpole et al. 2016; López-Corredoira 2017) do not trace the X-shaped component, thus emphasizing that the mechanism and the time of formation of the Galactic bulge are not 100 per cent clear.

8 CONCLUSIONS

We have used a parametric approach to find the best analytic function that describes the 3D bulge density observed by the VVV survey in the NIR. We test our modelling methodology on a range of mock data sets, including those generated using self-consistent simulations of the MW Galaxy. We thus prove that the technique delivers robust results and take care to quantify random and systematic uncertainties associated with the method. Our best-fitting model produces a faithful representation of the stellar density field in the inner Galaxy, with reassuringly low levels of residuals across the entire VVV field of view.

We find the bulge/bar to be ‘boxy’, with an axial ratio of [1:0.44:0.31] and an angle between the major axis and the Sun–GC line of sight of ∼20° (top row, Table 3). We show that there exists a strong degeneracy between the viewing angle and the dispersion of the RC absolute magnitude distribution. To reproduce the stellar line-of-sight distribution at a range of locations on the sky, models with larger intrinsic RC dispersions prefer larger viewing angles. We demonstrate that as a result of the degeneracy, assuming σ^RC = 0.18 leads to inferring viewing angles of ∼25° (closer to the value reported in Wegg & Gerhard 2013, α = (−26|$_{.}^{\circ}$|5 ± 2°)). However, in the analysis presented here, we also found evidence for a much larger intrinsic dispersion of the RC absolute magnitudes in the bulge area, namely σ_free^RC ∼ 0.26, larger than the latest measurement in the solar neighbourhood (σ^RC = 0.17 mag; Hawkins et al. 2017) where the younger thin disc population is dominant. We envisage that with new next-generation bulge surveys it should be possible to measure the true σ^RC across the bulge and thus, at last, settle the viewing angle argument.

The most pronounced residuals of the one-component bulge model appear to cluster tightly in the central regions of the Galaxy. One interpretation of this residual density distribution is that the core portions of the bulge are better described by a combination of two triaxial components: the main old one and a younger smaller one (last two rows, Table 3). Apart from this additional small-scale ellipsoidal component indicated by the low-latitude counts of the VVV RC giants, the only other noticeable outlier is the pattern induced by the split RC. We choose not to augment our bulge model with an ingredient responsible for the cross-shaped pattern, but provide an estimate of the total stellar mass locked in this population. We show that ∼7 per cent of the total bulge/bar could be residing in this non-ellipsoidal constituent, in agreement with some of the numerical simulations of the Galactic bar (e.g. Shen et al. 2010). However, this is probably only a lower estimate of the total fraction of stars contributing to the X-shape as it does not include the stars on orbits that could be fitted within the triaxial component.

Finally, the mass of the bulge/bar is a key ingredient in the recipe for the total density distribution of the MW. The bulge mass has to satisfy microlensing constraints (see e.g. Sumi et al. 2013) and ought to affect the dark matter density distribution as the halo reacts to the presence of baryons in the inner Galaxy (see e.g. Cole & Binney 2017). In light of the most recent measurements of the bulge mass function (e.g. Calamida et al. 2015), we argue for an MW bulge with a total mass of ∼2.3 × 10¹⁰ M_⊙, some 50 per cent of which is contributed by dark stellar remnants.

Acknowledgements

This work relies on DR2 of the VVV survey taken with the VISTA telescope under the ESO Public Survey programme ID 179.B-2002, pipeline processed, calibrated and provided by the Cambridge Astronomical Survey Unit. The authors are grateful for the enlightening feedback from all Cambridge Streams enthusiasts and the anonymous referee. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) through the Gaia Research for European Astronomy Training (GREAT-ITN) Marie Curie Network Grant Agreement no. 264895 and the ERC Grant Agreement no. 308024. SK thanks the United Kingdom Science and Technology Council (STFC) for the award of Ernest Rutherford fellowship (grant number ST/N004493/1). JS acknowledges support from the Newton Advanced Fellowship awarded by the Royal Society and the Newton Fund. ZYL is sponsored by Shanghai Yangfan Research Grant (no. 14YF1407700).

REFERENCES

APPENDIX A: COMPLETENESS OF THE VVV FIELDS AS A FUNCTION OF LATITUDE

APPENDIX B: MCMC RESULTS