Issue 
A&A
Volume 593, September 2016



Article Number  A108  
Number of page(s)  22  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201527535  
Published online  30 September 2016 
Models for the 3D axisymmetric gravitational potential of the Milky Way galaxy
A detailed modelling of the Galactic disk
^{1} Instituto de Astronomia, Geofísica e
Ciências Atmosféricas, Universidade de São Paulo, Cidade
Universitária, 05508090 São
Paulo, SP, Brasil
email: douglas.barros@iag.usp.br
^{2} UNIFEI, Instituto de Ciências Exatas,
Universidade Federal de Itajubá, Av. BPS 1303 Pinheirinho, 37500903
Itajubá, MG, Brasil
Received:
10
October
2015
Accepted:
7
July
2016
Aims. Galaxy mass models based on simple and analytical functions for density and potential pairs have been widely proposed in the literature. Disk models that are constrained solely by kinematic data only provide information about the global disk structure very near the Galactic plane. We attempt to circumvent this issue by constructing disk mass models whose threedimensional structures are constrained by a recent Galactic star counts model in the nearinfrared and also by observations of the hydrogen distribution in the disk. Our main aim is to provide models for the gravitational potential of the Galaxy that are fully analytical but also give a more realistic description of the density distribution in the disk component.
Methods. We produced fitted mass models from the disk model, which is directly based on the observations divided into thin and thick stellar disks and H I and H_{2} disks subcomponents, by combining three MiyamotoNagai disk profiles of any model order (1, 2, or 3) for each disk subcomponent. The MiyamotoNagai disks are combined with models for the bulge and dark halo components and the total set of parameters is adjusted by observational kinematic constraints. A model that includes a ring density structure in the disk, beyond the solar Galactic radius, is also investigated.
Results. The Galactic mass models return very good matches to the imposed observational constraints. In particular, the model with the ring density structure provides a greater contribution of the disk to the rotational support inside the solar circle. The gravitational potential models and their associated forcefields are described in analytically closed forms.
Conclusions. The simple and analytical models for the mass distribution in the Milky Way and their associated threedimensional gravitational potential are able to reproduce the observed kinematic constraints and, in addition, they are also compatible with our best knowledge of the stellar and gas distributions in the disk component. The gravitational potential models are suited for investigations of orbits in the Galactic disk.
Key words: Galaxy: disk / Galaxy: fundamental parameters / Galaxy: kinematics and dynamics / Galaxy: structure / methods: numerical / methods: statistical
© ESO, 2016
1. Introduction
Reliable models for the gravitational potential of the Galaxy are mandatory when studies of the structure and evolution of the Galactic mass components rely upon the characteristics of the orbits of their stellar content. In this sense, Galaxy mass models are regarded as the simplest way of assessing and understanding the global structure of the main Galactic components, and provide great insight into their mass distribution once a good agreement between the model predictions and the observations is obtained. A pioneer Galactic mass model was that of Schmidt (1956), which was contemporary with the early years of the development of radio astronomy and the first studies of the largescale structure of the Milky Way. With the subsequent improvement of the observational data, updated mass models have been undertaken by several authors (e.g. Bahcall & Soneira 1980; Caldwell & Ostriker 1981; Rohlfs & Kreitschmann 1988, among others), and with the advent of the Hipparcos mission and largescale surveys in the optical and nearinfrared, new observational constraints have been adopted in the more recent Galaxy mass models (e.g. Dehnen & Binney 1998; Lépine & Leroy 2000; Robin et al. 2003; Polido et al.2013).
In order to evaluate the capability of a given mass model in reproducing some observables, the force field associated with the resulting gravitational potential has to be compared with available dynamical constraints, such as the radial force in the plane given by the rotation curve and the force that is perpendicular to the plane of the disk along a given range of Galactic radii. Regarding the latter force, the associated masssurface density, to the best of our knowledge, is the integrated up to the height of 1.1 kpc of the Galactic midplane (Kuijken & Gilmore 1991; Bovy & Rix 2013), and as pointed out by Binney & Merrifield (1998), this constraint is not able to provide much information about the mass distribution some kiloparsecs above the plane. Owing to these shortcomings, a degeneracy in the set of best models is observed, which means that different mass models are able to reproduce the kinematic information of the observed data equally well. As stated by McMillan (2011), one possible way of circumventing such obstacles is by combining the kinematic data with star counts to improve the Galactic potential models and its force field above the plane.
Regarding the use of a Galactic potential model for the purpose of orbit calculations, the model of Allen & Santillan (1991) has been widely adopted. This model has the attractive characteristics that it is mathematically simple and completely analytical with closed forms for the potential and density, assuring both fast and accurate orbit calculations. Irrgang et al. (2013) recalibrated the Allen & Santillan (1991) model parameters using new and improved observational constraints.
The main goal of the present work is to provide a fully analytical, threedimensional description of the gravitational potential of the Galaxy, but with the novelty of expending considerable efforts in a detailed modelling of the disk component. The basic new aspects of the present Galactic mass model, all of which are related to the disk modelling process, can be summarized in the following way:

The structural parameters of the disk, i.e. scale length, scaleheight, and the radial scale of the central hole, are based on theGalactic star counts model in the nearinfrared developed byPolido et al. (2013, hereafter PJL)for the case of the stellar disk component; for the gaseous diskcounterpart, we adopt recent values returned by surveys of thedistribution of hydrogen, atomic H I, andmolecular H_{2}, in the disk.

The density and potential of the disk components are modelled by the commonly used MiyamotoNagai disk profiles (Eqs. (4) and (5) of Miyamoto & Nagai 1975), but here we also attempt to make use of the higher model orders 2 and 3 of MiyamotoNagai disks (Eqs. (6)–(9) of the abovereferred paper) to better fit some of the disk subcomponents. The approach followed to construct the MiyamotoNagai disks is based on the approach presented by Smith et al. (2015).

A model with a ring density structure that is added to the disk density profile has been studied, with the ring feature placed somewhat beyond the solar Galactic radius. The inclusion of such ring structure is motivated by the attempt of modelling the local dip in the observed Galactic rotation curve that is also placed a little beyond the solar orbit radius. An explanation for the existence of such ring density structure is given by Barros et al. (2013, hereafter BLJ).
The organization of this paper is as follows: in Sect. 2, we present the details of the mass models of the Galactic disk and the steps for the construction of MiyamotoNagai disk versions of the “observed” disks. In Sect. 3, we give the expressions for the bulge and dark halo components as well as the functional form for the gravitational potential associated with the ring density structure. The group of observational constraints adopted for the fitting of the models are presented in Sect. 4, while the fitting scheme and estimation of uncertainties are presented in Sect. 5. In Sect. 6, we analyse the results of each mass model via a direct comparison with other models in the literature. Concluding remarks are drawn in the closing Sect. 7.
2. Mass models for the disk of the Galaxy
We model the Milky Way disk and separate it into the stellar (thin and thick disks) and gaseous (H I and H_{2} disks) components. In the following subsections, we present the observational basis taken as prior information to constrain the values of the parameters of the models and the steps taken to construct the mass and potential disk models. In this paper, the cylindrical coordinates (R, z) are used in the density and potential expressions. The solar Galactic radius is denoted as R_{0}.
2.1. The “observationbased” disk model
2.1.1. The stellar component
Our models for the density distribution in the thin and thick stellar disks of the Galaxy are based on the structural disk parameters presented by PJL. These authors performed a star counts model of the Galaxy using nearinfrared data of the 2MASS survey (Skrutskie et al. 2006) with lines of sight covering the entire sky and including the Galactic plane. The authors explored the parameter space and estimated its optimal values with statistical methods such as the Markov chain Monte Carlo (MCMC; Gilks et al. (1996)) and the nested sampling (NS) algorithm (Skilling 2004). Using a modified exponential law, PJL modelled the radial profile of the density of each subcomponent of the stellar disk based on the Galactic model of Lépine & Leroy (2000). Such profile is equivalent to the Freeman Type II disk brightness profile, which contains a depletion in the centre with respect to a pure exponential law (Freeman 1970; Kormendy 1977). The stellar surface densities Σ_{d★} for the thin and thick disks can then be written as (1)where Σ_{0d★} corresponds to the local disk stellar surface density (at R = R_{0}), R_{d} is the radial scale length, and R_{ch} is the radial length of the “central hole” in the density of each stellar disk i subcomponent (i = thin, thick). The hypothesis that the Galactic disk is hollow in its centre has been justified by some models that use observational data at infrared bands to describe the inner structure of the Galaxy (e.g. Freudenreich 1998; Lépine & Leroy 2000; LópezCorredoira et al. 2004; Picaud & Robin 2004). In the particular case of the PJL model, only the thin disk needs a density depression in its inner part; in contrast, the thick disk can be described by a simple radial exponential decay, i.e. R_{ch thick} = 0.
Structural parameters, local surface densities, and masses of the disk components taken as observational prior information to the Milky Way modelling.
For the stellar density variation perpendicular to the Galactic plane, PJL modelled the vertical profile of the thin and thick disks by exponential laws with scale height h_{z}. In that case, the authors introduced the variation of the scale height with the Galactic radius, h_{z} = h_{z}(R), which is known as the flare of the disk. Recently, Kalberla et al. (2014) compiled some published results in the literature and found compelling evidence for the increase of the scale heights with Galactocentric distance for different stellar distributions. In the present study, however, we do not attempt to model such a function for h_{z}(R), and we consider the scale height as a constant along the Galactic radius and with a value equal to the local scale height h_{z0} (at R_{0}) estimated by PJL, for each thin and thick disk. The reason for this approximation is justified by the fact that the introduction of the flaring of the disk requires a more careful analysis with respect to the form of the gravitational potential that would result by such a distribution of density. The volume density for both thin and thick disks is written in the form (2)The adopted values for the structural parameters of the thin and thick disks, i.e. the scale lengths R_{d}, radii of the central hole R_{ch} and scale heights h_{z}, are, as mentioned before, the bestfitting values reported in the PJL model; these values are listed in Table 1.
The local stellar surface densities for both thin and thick disks (Σ_{0d★} in Eq. (1)) are based on the model of Flynn et al. (2006, see also Holmberg & Flynn 2000, 2004. These authors discriminate between the contributions for the total Σ_{0d★} generated by different stellar components, namely mainsequence stars of different absolute magnitudes, red giants and supergiants, stellar remnants (white dwarfs, neutron stars, black holes), and brown dwarfs. The mainsequence stars and giants contribute with Σ_{°} = 28.3M_{⊙} pc^{2}, which can be compared with the recent determination by Bovy et al. (2012) of Σ_{°} = 30M_{⊙} pc^{2} using the SEGUE spectroscopic survey data. The stellar remnants and brown dwarfs in the Flynn et al. (2006) model contribute with Σ_{•} = 7.2M_{⊙} pc^{2}. The combination of the Bovy et al. value for Σ_{°} and the Flynn et al. value for Σ_{•}, as a constraint to the local stellar surface mass density, yields Σ_{0d★} = Σ_{°} + Σ_{•} = 37.2M_{⊙} pc^{2}, which is the same value adopted by Read (2014). Separating this last value between the thin and thick disks, we take Σ_{0d★, thick} = 7.0M_{⊙} pc^{2} for the thick disk, as in the Flynn et al. (2006) model, and Σ_{0d★, thin} = 30.2M_{⊙} pc^{2} for the thin disk, where we assigned the brown dwarfs and stellar remnants to the thin disk for practical purposes. These local surface densities along with the scale heights result in the local volume densities in the midplane of the Galaxy, i.e. ρ_{0d★, thick} = 0.0055M_{⊙} pc^{3} for the thick disk and ρ_{0d★, thin} = 0.0736M_{⊙} pc^{3} for the thin disk (or ρ_{0d★, thin} = 0.0561M_{⊙} pc^{3}, considering only mainsequence stars and giants). The thicktothin disk density ratio, ρ_{0d★, thick}/ρ_{0d★, thin} ~ 10% (neglecting the stellar remnants/brown dwarfs contribution), is close, given the errors, to the value measured by Jurić et al. (2008) of 12%. The values adopted as constraints for Σ_{0d★, thin} and Σ_{0d★, thick} are listed in Table 1. In Table 1, we also give the total masses M_{d} calculated for the thin and thick disks and the radial scale length R_{dexp} relative to the region of each disk subcomponent that presents the density exponential decay, and which is used in the modelling process described in Sect. 2.2.3. Since the thick disk is modelled by a single exponential, R_{d thick} = R_{dexp thick}.
2.1.2. The gaseous component
We consider the atomic H I and molecular H_{2} gaseous disks the major contributors to the density of the interstellar medium (ISM) with the proper correction factor for He. The radial profile of the surface density for these components is chosen as (3)with the exponents n = 3/2 for the H I and n = 1 for the H_{2} gas disk components, respectively. Similarly to Eq. (1), Σ_{0d} corresponds to the local (H I, H_{2}) surface density, R_{d} is the radial scale length, and R_{ch} is the radius of the central hole in the density distribution of H I and H_{2} disks. The form for the gas surface densities presented in Eq. (3), together with the values for the structural parameters R_{d} and R_{ch} of the H I and H_{2} disks listed in Table 1, are chosen to reproduce as closely as possible the density distributions observed in H I by Kalberla & Dedes (2008, cf. their Fig. 3), and in H_{2} by Nakanishi & Sofue (2006, cf. their Fig. 9).
The volume density for both H I and H_{2} disks is given in the form (4)with the scale heights z_{1/2} corresponding to the half width at half maximum of the density peaks of H I and H_{2} in the Galactic midplane. The Gaussian profile for the vertical distribution of hydrogen in the form presented in Eq. (4) has been widely adopted in the literature to represent the H I and H_{2} distributions (e.g. Sanders et al. 1984; Amôres & Lépine 2005; PJL). As pointed out by Amôres & Lépine (2005), the Gaussian vertical distribution for the gaseous disk correctly fits the observations, and is also an expected solution for hydrostatic equilibrium considerations.
The scale height of the H I distribution in the Milky Way has since long ago been observed to increase systematically with Galactic radius (e.g. Lozinskaya & Kardashev 1963; Burton 1976; Kalberla & Dedes 2008, among others). The flaring of the H I disk is naturally expected when one considers the fact that, for the case of hydrostatic equilibrium, the gravitational force perpendicular to the disk plane must balance the gas pressuregradient force, and since the vertical velocity dispersion σ_{z} of the H I gas is approximately constant with radius R (e.g. Spitzer 1968), the scale height of this component must increase with R. The Galactic distribution of molecular gas also shows a flaring that is consistent with that observed in H I (Kalberla et al. 2007, and references therein). Although the flaring of the gaseous disk component is a very important feature to be included in any realistic mass model of the Galaxy, we do not attempt to model this flaring for the same reason exposed in the case of the stellar disk. Therefore, here we consider the scale heights of the H I and H_{2} disk components as independent of radius and with values equal to the local values (at R_{0}) calculated from the model of Amôres & Lépine (2005, cf. their Eq. (7)). Table 1 lists the adopted prior values for the H I and H_{2} disk scale heights.
The model of Holmberg & Flynn (2000) for the ISM component, which is based on the original multiphase model of Bahcall et al. (1992), discriminates the contributions between the molecular gas, the warm (ionized) gas, and a split of the neutral H I into cold and hot components. The total local surface density in the gas form is Σ_{0dg} = 13M_{⊙} pc^{2} with an uncertainty of ~50%, according to this model. More recently, Hessman (2015) has warned about the underestimation of these traditional determinations of the neutral and molecular gas densities. According to this author, substantial amounts of “dark” gas in the form of optically thick cold neutral hydrogen and COdark molecular gas are known to contribute to the ISM density, which must raise the estimates of the local midplane gas densities by as much as ~60%. Moreover, if local density features such as the Local Bubble or the local spiral arms structure are taken into account, the corrections for a larger Σ_{0dg} are even higher. Therefore, as observational constraints to our gaseous disk model, we take the following values for the local surface gas densities: Σ_{0d HI} = 15M_{⊙} pc^{2}, as in the Hessman (2015) updated estimate, which already takes into account the correction factor of 1.36 for the mass in helium (He); and Σ_{0d H2} = 3M_{⊙} pc^{2}, as in the Holmberg & Flynn (2000) model. The warm gas disk component, which contributes locally with 2 M_{⊙} pc^{2} (Holmberg & Flynn 2000; Hessman 2015), is here incorporated to the H I disk for practical purposes, increasing the value of Σ_{0d HI} to 17 M_{⊙} pc^{2}. The adopted values for the local surface densities, total masses M_{d}, and radial exponential scale lengths R_{dexp} of the H I and H_{2} disks are listed in Table 1.
As can be seen from Table 1, our “observationbased” disk model comprises a total mass M_{d} = 4.47 × 10^{10}M_{⊙}, which is M_{d★} = 3.06 × 10^{10}M_{⊙} in stars and M_{dg} = 1.41 × 10^{10}M_{⊙} in the gaseous form. The total local disk masssurface density is Σ_{0d} = 57.2M_{⊙} pc^{2}, which can be compared with the updated density estimates of 58 M_{⊙} pc^{2} by Hessman (2015) and 54.2 ± 4.9M_{⊙} pc^{2} by Read (2014). The local disk midplane volume density is ρ_{0d} = 0.138M_{⊙} pc^{3}. This value is somewhat larger than the estimate of the local dynamical mass density of 0.102 ± 0.010M_{⊙} pc^{3} by Holmberg & Flynn (2000) of which 0.095 M_{⊙} pc^{3} is in visible matter. However, the density corrections discussed by Hessman (2015) should increase these traditional estimates to local densities as large as ~0.120 or even ~0.160M_{⊙} pc^{3}. Figure 1 shows the radial distribution of the surface density of our “observationbased” disk model, and the surface densities of each disk subcomponent. The surface densities of the thin and thick stellar disks drop to ~1M_{⊙} pc^{2} at the radius R ~ 15 kpc, while the H_{2} disk reaches this value at R ~ 10 kpc. The H I disk subcomponent extends to larger radii with a drop of its surface density to ~0.1M_{⊙} pc^{2} at R ~ 30 kpc, which is very similar to the mean surface density profile of H I derived by Kalberla & Dedes (2008).
The derivation of the gravitational potential directly from the density distributions in Eqs. (2) and (4), through the Poisson equation ∇^{2}Φ = 4πGρ, must involve the use of some mathematical tools and/or numerical techniques, such as the multipole expansion and numerical interpolation employed by Dehnen & Binney (1998) in their determination of the potential and force field associated with their Galactic mass models. Another way of handling such a task is to approximate the disk density distribution by an analytical functional form for which the associated potential is also known analytically, as is the case of the densitypotential pairs of MiyamotoNagai disks. It is known that a single MiyamotoNagai disk only provides a rough approximation to a real galactic disk density profile. But we show later in the next sections that the combination of three MiyamotoNagai disks is able to give good approximations to the observed Galactic disk mass distribution, and thereby for the gravitational potential for given ranges of Galactic radii and heights above the disk midplane. In the present paper, we adjust combinations of MiyamotoNagai disks to each component of our “observationbased” disk model. The disk potential then described in a complete analytical form provides quicker computations of its related force field, which is suitable for fast calculations of galactic orbits of large samples of stars or testparticles in a numerical simulation. This is one of the main advantages that we pursue for our future studies concerning stellar orbits in the Galactic disk. In the next subsection, we describe the method employed to the search for the best set of MiyamotoNagai disks (hereafter MNdisks) that reproduce the main features of the density distributions in Eqs. (2) and (4).
Fig. 1
Radial profile of the surface density of our “observationbased” model for the Milky Way’s disk. The curves indicate the surface densities of the thin disk (green), thick disk (brown), H I disk (blue), H_{2} disk (violet), and the total disk (red). 
2.2. Reproducing the disk mass model with MiyamotoNagai disks
2.2.1. The ToomreKuzmin disks
As a first step, we make use of the family of disk models of infinitesimal thickness (the razorthin disks) introduced by Toomre (1963) and Kuzmin (1956), whose densities are written in the form ρ(R,z) = Σ(R) δ(z). The “model 1” of this family, of what we refer to as ToomreKuzmin disks (TKdisks), is described by a surface density which is written as (5)where M_{TK1} is the total disk mass and a_{TK1} is related to the radial scale length of the disk. The choice of use of TKdisks, in the first step, is because all the information about the disk surface density can be recovered adjusting only two parameters, M_{TK} and a_{TK}. As shown later, the threedimensional structure of the disks is obtained after “inflating” the TKdisks with the introduction of the parameter b related to the scale height, in the same way as in the original method of Miyamoto & Nagai (1975).
From Eq. (5), it can be seen that the density Σ_{TK1} decreases with R^{3} at large radii. This is a slower decrease when compared to the observed exponential fall off of the brightness profiles in galaxy disks (Binney & Tremaine 2008). Therefore, a single TKdisk (and the correspondent MNdisk) poorly matches the surface density profile of a radially exponential disk (Smith et al. 2015). This motivated the use of a combination of three MNdisks (related to the “model 1” of TKdisks, Eq. (5)), which was first introduced by Flynn et al. (1996) and since then used for modelling the disk of the Milky Way (Smith et al. 2015, and references therein). In this procedure, each MNdisk has a different scale length a and mass M, with one of the masses with a negative value. The MNdisk with negative mass also has the largest scale length. This feature helps decrease the surface density at large radii (and improve the fit to an exponential disk), but also leads to the undesirable occurrence of regions with negative densities near the midplane of the disk and at large radii (Smith et al. 2015).
It is also known that the density Σ_{TK} in Eq. (5) is just the first of a family of possible forms for the potentialdensity pair that obey the Poisson equation. Toomre (1963) showed that, since the Poisson equation is linear in Φ and ρ, new potentialdensity pairs can be obtained through the derivation of Φ /an times with respect to a^{2} (Binney & Tremaine 2008). For example, the densities for “model 2” and “model 3” of the TKdisks are described by (6)and (7)respectively (Miyamoto & Nagai 1975). One can see that the densities Σ_{TK2} and Σ_{TK3} decrease with R^{5} and R^{7}, respectively, at large radii. In this way, a combination of models 2 or 3 of TKdisks could, in principle, not only result in an equally good fit to an exponential disk, but also circumvent the problem of the appearance of negative densities. In our present case, we also need to include a TKdisk with negative mass to the modelling of the central density depression that occurs in the thin stellar disk, the H I and H_{2} disk subcomponents. The difference here is that the scale length of this component with negative mass does not need to be the largest, thus avoiding the negative densities at large radii. Therefore, for each disk subcomponent (thin, thick, H I, H_{2}), we search for the combination of three TKdisks of each model (1, 2, 3) that better reproduces the radial profile of the surface density of each subcomponent. We take, for example, the case of the thin disk. We search for the set of parameters ^{1} of “model 1” of three TKdisks that generates the best fit to the radial distribution of surface density of the thin disk. These three TKdisks result in the total surface density . The same procedure is carried out with the three TKdisks in the forms of “model 2” and “model 3”, which result in the total densities and , respectively. We quantify ξ^{2} as the sum over a given radial range of the squares of the residuals between the surface density of the thin disk (obtained from Eq. (1) and parameters from Table 1) and the surface density resulted from the bestfitting combination of the three TKdisks, for each one of the models, (8)where k refers to the disk subcomponent (k = thin in the above example) and j denotes the model order of TKdisks; N_{R} is the number of radial bins over which the sum is calculated. The TKdisk model order that results in the lowest value for ξ^{2} is that which we consider the model that better represents the surface density Σ_{dk} of the disk subcomponent. The search for the best sets of parameters corresponding to the best TKdisk models for j was carried by applying the global optimization technique based on the crossentropy (CE) algorithm for parameters estimation. The CE algorithm (Rubinstein 1997, 1999; Kroese et al. 2006) provides a simple adaptive way of estimating the optimal set of reference parameters in the fitting process. Recent investigations have applied the CE technique to some astrophysical problems, e.g. Caproni et al. (2009), Monteiro et al. (2010), Monteiro & Dias (2011), Martins et al. (2014), Dias et al. (2014), Caetano et al. (2015), among others. Since we also employ the CE method on other occasions throughout this paper, in the next subsection we give a brief description of the main steps of the algorithm. Table 2 summarizes the model of three TKdisks found to better reproduce the surface density of each one of the disk subcomponents and the corresponding ξ^{2} (Eq. (8)) minimized within the CE algorithm.
Models of 3 TKdisks that provide the best match to the surface density of each disk subcomponent.
2.2.2. The crossentropy algorithm
For a detailed presentation and description of the CE method, see the papers by Monteiro et al. (2010), and Dias et al. (2014). Here, we give a brief overview of the method and how it works.
Supposing we have a set of data D with individual points d_{1}, d_{2}, ..., d_{ND} and we wish to study it in terms of an analytical model Θ with a vector of parameters θ_{1}, θ_{2}, ..., θ_{Np}. The main goal of the CE continuous multiextremal optimization method is to find a set of parameters for which the model provides the best description of the data, based on some statistical criterion. This is performed by randomly generating N independent sets of model parameters , where , under some chosen distribution, and the subsequent minimization of an objective function S(X) is used to transmit the quality of the fit during the run process. As the method converges to the “theoretical” exact solution, then S → 0, which means .
The CE method performs an iterative statistical coverage of the parameter space, where the following is carried out in each iteration (Monteiro et al. 2010):

(i)
random generation of the initial parameter sample, respectingsome underlying distributions and predefined criteria;

(ii)
selection of the best candidates based on some mathematical criterion, which will compose the elite sample array – samples with the lowest values for the objective function S(X);

(iii)
random generation of updated parameter samples from the previous best candidates, i.e. the elite sample, to be evaluated in the next iteration;

(iv)
optimization process repeats steps (ii) and (iii) until a prespecified stopping criterion is fulfilled.
In all the implementations of the CE algorithm carried out in this work, we followed the general step by step presented in Sect. 2.2 of Monteiro et al. (2010). The tuning parameters used in the run processes were N = 3 × 10^{3} sets of trial model parameters per iteration; a maximum number of 100 iterations; N_{elite} = 100, which is the number of sets of trial parameters that return the best solutions (the lowest values for the objective function) at a given iteration and that will be used to estimate the distribution parameters for the next iteration. For the smoothing parameters that reduce the convergence speed of the algorithm, preventing it from finding a nonglobal minimum solution, we have used α = 0.9, α′ = 0.7, and q = 7. For more details about these parameters and how they are implemented in the algorithm, see Monteiro et al. (2010).
2.2.3. The MiyamotoNagai disks – following the approach developed by Smith et al. (2015)
Once the models of TKdisks have been found, the next step is to “inflate” them vertically to find the correspondent MNdisks. As firstly performed by Miyamoto & Nagai (1975), this is carried out by replacing the term (a +  z ) in the potential function, which is associated with the density of a given model of TKdisk, with the term , where b expresses the vertical scale height of the disk. In this way, corresponding to the three models of TKdisks expressed by Eqs. (5)–(7), we have the following threedimensional density functions that compose the first three models of MNdisks: where . In the above expressions, for the sake of good readability, we avoided the excessive use of subscripts in the parameters M, a, b, and ζ that could distinguish between each one of the models. Here, we make it clear that for each disk subcomponent modelled with a combination of three TKdisks of a given model, we use the corresponding combination of three MNdisks of the equivalent model.
Parameters of the three MNdisk combinations for modelling each subcomponent of the Galactic disk.
Fig. 2
Lefthand panel: radial distribution of the surface mass density of the Galactic disk model; the “observationbased” disk model (blue curve) and the resultant from the sum of all the 3 MNdisk combinations fitted to each disk subcomponent (red curve). Righthand panel: vertical profile of the volume density of the Galactic disk model as a function of the height z from the midplane and at three different arbitrary radii: R = 2 kpc (solid lines); R = 8 kpc (R_{0}) (dashed lines); and R = 15 kpc (dotted lines). The blue curves are also related to the “observationbased” disk model and the red curves to the total 3 MNdisk combinations. 
The task now is to find the best values for the parameter b that reproduce the vertical density distribution of each subcomponent of the Galactic disk. For this purpose, we follow the approach recently developed by Smith et al. (2015). In that work, the authors create a procedure for modelling simple radially exponential disks from the combination of three MNdisks of “model 1”, allowing the construction of disks of any mass, scale length, and for an ample range of thicknesses, from infinitely thin disks to approximately spherical systems. The basic idea of the method is as follows: Starting from the model composed of the three TKdisks (for which b = 0) that better match the radial profile of the disk surface density, one chooses the corresponding model composed of three MNdisks with a small value for the b parameter (with the aim of remaining near the solution for the infinitesimal thickness disk), and whose parameters M and a of the MNdisks deviate from small amounts with respect to those of the TKdisks. With the continuous variation of b in small steps, the parameters M and a are searched for at intervals close to the solution found in the previous step, ensuring a smooth variation of both M and a with the variation of the disk thickness b. In fact, Smith et al. (2015) analyse the variation of the ratios M/M_{d} and a/R_{d} as a function of the variation of the ratio b/R_{d}, where M_{d} is the mass and R_{d} is the radial scale length of the exponential disk to be modelled. Therefore, for each small variation of the ratio b/R_{d}, one searches for the parameters M and a of the three MNdisks whose integral of the density ρ(z) in the vertical direction better reproduces the radial profile of the disk surface density Σ(R).
To relate the scale height b of the MNdisks with the scale height h_{z} of the vertically exponential disk, Smith et al. (2015) follow an approach that is analogous to that described above. In this case, the authors first sum up the fractional differences between the density of the three MNdisks combination and the exponential density, measured vertically from the midplane up to z = 5b. These authors vary the ratio b/R_{d} to minimize the sum, and then find the best match to the exponential distribution.
In this section, we only present the values of the parameters [M,a,b] calculated for each combination of three MNdisks that fit each disk subcomponent with their fixed values for [M_{d},R_{dexp},h_{z}] listed in Table 1. However, in Appendix A, we also present the variations of the ratios M/M_{d} and a/R_{dexp}, as a function of b/R_{dexp}, and the relations between b/R_{dexp} and h_{z}/R_{dexp} for each disk subcomponent, calculated in the same way as in Smith et al. (2015). Therefore, disks of any mass M_{d} and scales R_{dexp} and h_{z} can be built up using those relations, but remembering here that they are suited for models of disks with central density depressions. Since the majority of our disk subcomponents are no longer modelled by a single radially exponential law (cf. Eqs. (1) and (3)), coupled with the fact that we also use models of higher order (2 and 3) than model 1 of MNdisks, our calculated relations between the aforementioned parameters are somewhat different from those presented in Smith et al. (2015). But with equivalent utility, the relations in Appendix A allow the construction of models for disks with masses and sizes that are different from those modelled in this work.
In our search for the best values of the parameters [M,a,b] of the MNdisks, we also employed the CE technique described in the previous subsection, but now look for the minimization of the sum of the squares of the residuals between the surface density of the disk subcomponent and the vertically integrated volume density of the three MNdisk combinations (12)for the relations of M/M_{d} and a/R_{dexp} with b/R_{dexp}, and the minimization of the sum of the squares of the fractional differences between the volume densities of the three MNdisk combinations and the disk subcomponent over a given vertical range and at R = R_{0}(13)for the relations between b/R_{dexp} and h_{z}/R_{dexp}.
We return to the case of the thin stellar disk for exemplification. With the ratio h_{z}/R_{dexp} calculated from the values listed in Table 1 for the thin disk, and putting it into Eq. (A.1) with the coefficients for the thin disk given in Table A.1, both from Appendix A, one finds the corresponding value for the ratio b/R_{dexp}. Substituting this last value into Eq. (A.2) with the coefficients given in Table A.2, also from Appendix A, one obtains the values of the parameters M_{i} and a_{i} for the combination of three MNdisks of model 3 (Eq. (11)), which matches the density distribution ρ_{d★, thin} of the thin stellar disk better. We emphasize here that a single scale height b is used for all the three MNdisks of the combination, as originally carried out by Flynn et al. (1996) and followed by Smith et al. (2015). The entire procedure described above, which was illustrated with the thin disk, is repeated with the other subcomponents (thick disk, H I, and H_{2} disks) using their corresponding equations and tables in Appendix A. Table 3 lists the values of the parameters [M,a,b] of each combination of three MNdisks that fit each disk subcomponent and the correspondent models of MNdisks that had been previously found with modelling of TKdisks.
Fig. 3
Contours of isodensities log ρ for the disk models built in the present work. Lefthand panel: the “observationbased” disk model; righthand panel: the sum of all three MNdisk combinations fitted to each disk subcomponent. The values of the contours lie in the range log ρ = [−6; +1], with ρ in M_{⊙} pc^{3}. 
Figure 2 shows, in the lefthand panel, the radial profile of the surface mass density for the “observationbased” disk model (blue curve) and the total surface density resultant from all the three MNdisk combinations fitted to each disk subcomponent (red curve). A good agreement is seen between these two curves, which denotes the great effectiveness of the method. The fractional differences between these surface density distributions are <10% at 1 kpc ≲ R ≲ 21.5 kpc, reaching 50% at R ~ 26 kpc and 100% at R ~ 29 kpc, where the absolute values of the densities become lower than ~0.5M_{⊙} pc^{2}. These features show the tendency of the MNdisks to return larger densities at large radii. The local disk surface density returned by the MNdisk fit models is Σ_{0d} = 59.4M_{⊙} pc^{2}, which can be considered close to the prior value of 57.2 M_{⊙} pc^{2} of the “observationbased” disk model, if we take the uncertainty of ± 4.9M_{⊙} pc^{2} over this quantity as estimated by Read (2014). The local surface density integrated within 1.1 kpc of the disk midplane returned by the MNdisk fit models is Σ_{0d1.1 kpc} = 57.0M_{⊙} pc^{2}, in agreement with the values reported by Bienaymé et al. (2006). On the righthand panel of Fig. 2, we show the distribution of the volume density as a function of the height z from the disk midplane, taken at three arbitrary radii: R = 2 kpc (solid line), R = 8 kpc (dashed line), and R = 15 kpc (dotted line). The blue and red curves also correspond to the density distributions of the “observationbased” disk model and the MNdisk fit models, respectively. The fractional differences between the volume density distributions are lesser than 25% for  z  ≲ 2.3 kpc at R = 2 kpc, for  z  ≲ 1 kpc at R = 8 kpc, and for  z  ≲ 0.35 kpc at R = 15 kpc, reaching values that are greater than 25% beyond these heights. These features denote the systematic trend of the MNdisks to return higher densities at greater distances from the disk plane, while at relatively small z the density distributions show roughly good matches at radii up to at least R ~ 22 kpc. As pointed out by Smith et al. (2015), it is impossible for the three MNdisk models to perfectly reproduce the vertically exponential density profile or even the sech^{2}(z) type profile, since they are mathematically distinct. The same can be stated about the vertical Gaussian profile adopted for the H I and H_{2} disks in the present work; at heights lower than 1 kpc the three MNdisks start presenting great deviations from the Gaussian profile. This can be noticed by analysing the distributions of ρ(z) calculated at R = 15 kpc (dotted lines in the righthand panel of Fig. 2), where the density of the H I disk overtakes those of the other disk subcomponents (cf. Fig. 1).
Figure 3 shows the contours of isodensities in the form log ρ (ρ in M_{⊙} pc^{3}) in the meridional plane of the Galaxy: for the “observationbased” disk model (lefthand panel) and for the total three MNdisk combination fit models (righthand panel).
2.3. The gravitational potential of the disk
The gravitational potential expressions related through Poisson equation to the densities of the MiyamotoNagai disk models, given by Eqs. (9)–(11), are written as “model 1”, “model 2”, and “model 3”, respectively, in the form where again . Therefore, the gravitational potential of the thin stellar disk Φ_{d thin} is modelled with three components of the potential Φ_{MN3} (Eq. (16)), whose parameters M_{i}, a_{i}, and b (i = 1, 2, 3) are given in the first row of Table 3. Equivalently, the thick stellar disk potential Φ_{d thick} is modelled with three components of the potential Φ_{MN1} (Eq. (14)), the H I disk potential Φ_{d HI} is modelled with three components of Φ_{MN2} (Eq. (15)), and the H_{2} disk potential Φ_{d H2} with three components of the potential Φ_{MN3} (Eq. (16)); the parameters M_{i}, a_{i}, and b for these disks are given in the second, third, and fourth rows, respectively, of Table 3. The total gravitational potential of the disk Φ_{d} is then given by the sum of all the combinations of three MNdisks that fit each Galactic disk subcomponent, (17)Figure 4 shows the contours of equipotential curves in the plane R − z for the gravitational potential of the disk Φ_{d}, which is expressed in Eq. (17).
Fig. 4
Contours of equipotentials in the Galactic meridional plane for the gravitational potential resulting from the disk model composed of combinations of MNdisks, as described in Sect. 2.2.3. The contours of Φ_{d} are given in 10^{4} km^{2} s^{2}, as indicated in the colourbar. 
3. Galactic models
3.1. Bulge and dark halo components
For the mass density distribution in the spheroidal component of the Galaxy, here assumed as composed of the central bulge and the stellar halo as a smooth extension of the bulge itself, we adopt the profile proposed by Hernquist (1990), which reproduces the r^{1/4}de Vaucouleurs (1977) surface brightness law over a large range of galactic radius, (18)where M_{b} is the total mass and a_{b} is the scale radius of the bulge. The gravitational potential associated with the density distribution in Eq. (18) was shown by Hernquist (1990) to be in the form (19)In the PJL star counts model, the bulge is modelled as an oblate spheroid, which is described by a modified Hernquist profile with three free parameters: the oblateness parameter κ, the spheroidtodisk local stellar density ratio for normalization of the density profile, and the spheroid scale radius a_{H} = 0.4 ± 0.1 kpc. Considering the local stellar density of our “observationbased” disk model, the mass of the spheroid of the PJL model inside R = 3 kpc (a radius that encompasses ~90% of the total spheroid mass) results in M_{sph}(R< 3 kpc) ≈ 2.2 × 10^{10}M_{⊙}. We use this information to constrain the initial ranges of the parameters M_{b} and a_{b} of our bulge models to search for their best values in the fitting procedure described in Sect. 5. We adopt the initial guess interval M_{b} = 2.4–2.8 × 10^{10}M_{⊙} for the bulge mass; we use a_{b} = 0.3–0.5 kpc for the scale radius. These intervals comprise masses inside R = 3 kpc for our bulge models with possible values in the range 1.8–2.3 × 10^{10}M_{⊙}.
The fact that the observed rotation curve, which is adopted to constrain the Galactic models, is nearly flat in the interval R_{0} ≲ R ≲ 2R_{0} (see Sect. 4), and besides the fact that in our models (see Sect. 6) the sole disk plus bulge density distributions are not capable of giving such support to the rotation curve at these radii, leads us to take into consideration the contribution of an extra mass component to which we refer by the oftenused term “dark halo”. Since there is still some tension in the explanations for the behaviour of the rotation curves at large radii of some external spiral galaxies, as for the Milky Way, we do not reflect on the reality of such a dark halo component nor its material content, although the dark halo contributes an extra term of mass in our Galactic models. For this reason, we do not go deeper in the modelling of such a component and take simple forms for its density and potential functions.
We model the dark halo component with a logarithmic potential of the form (e.g. Binney & Tremaine 2008) (20)where r_{h} is the core radius; v_{h} is the circular velocity at large r, i.e. relative to the core radius; and q_{φ} is the axis ratio of the equipotentials. For simplicity, we consider a spherical dark halo, for which q_{φ} = 1. The density distribution corresponding to the potential Φ_{h} is given by (for q_{φ} = 1) (21)This profile yields a flat rotation curve at large radii. However, as we make clear in the following sections, we are interested in a description of the Galactic potential that can be suitable for the study of stellar orbits that do not reach great extensions of the outer disk, so the real behaviour of the rotation curve of the Galaxy at very large radii is not of concern in the present study.
The values of the parameters M_{b} and a_{b} that describe the bulge component, r_{h}, and v_{h} for the dark halo are found after fitting the contributions of such components to the Galactic rotation curve as well as to other observational constraints as described in Sect. 4.
3.2. The disk component
The model constructed for the mass distribution in the disk of the Galaxy, and its associated gravitational potential, is described in detail in the subsections of Sect. 2. In general, the overall disk is modelled by a superposition of distinct models of MiyamotoNagai disks; these models are comprised of combinations of three MNdisks for each one of the four subcomponents: thin, thick, H I, and H_{2} disks yielding 12 MNdisks in total. Each MNdisk has three free parameters (M, a, b), but since each combination of three MNdisks is modelled with a single value for b (see Sect. 2.2.3), we have a total of 28 parameters for the construction of the disk (all of these parameters are given in Table 3).
In the present work, in contrast to other studies, the values of the structural parameters of the disk (here represented by the length scales a and b of the MNdisks) are not constrained by the kinematic information from the rotation curve of the Galaxy. Here, we opt to constrain such parameters based only on the information brought by the starcounts model of PJL and the observed distribution of the gaseous component in the disk. However, some of the disk scale parameters in the Galactic model of PJL possess considerable uncertainties, as for instance the radial length of the central hole of the thin disk kpc, which the kinematic data would help to reduce. Although this can represent a drawback of our approach, in Sect. 5 we give an estimate of the uncertainties on the parameters a and b of the MNdisks based on the uncertainties on the scale parameters of the “observationbased” disk model. Our reason for leaving the scale parameters fixed at the values from Table 3 is justified owing to the way the MNdisks were constructed; the scale lengths a were found as functions of the scale heights b, and they should preserve the relationships during the fitting of the kinematic data. However, a proper fitting of the kinematic constraints should vary these scale lengths in an independent way to probe the correlations among the parameters. Moreover, it is well known in the literature that the kinematic data commonly used for the model fittings do not provide strong constraints on the vertical distribution of the Galaxy’s mass (e.g. Dehnen & Binney 1998; McMillan 2011). Therefore, different sets of disk scale heights would be found able to reproduce the kinematic data equally well. Another problematic issue in varying the scales a and b of all the MNdisks in the fitting of the models, taking into account the uncertainties on the values of R_{d}, R_{ch}, and h_{z} of all the disk subcomponents, is the long computing time involved in the whole process. We recognize this disadvantageous point of our approach in fitting MNdisks to the observed disk; other studies in the literature that do not make use of this approximation are able to vary the scale parameters of their disk potential models in a more straightforward way (e.g. McMillan 2011; Piffl et al. 2014a).
The total mass of the disk, otherwise, is left to be constrained not only by the rotation curve that provides information on the dynamical mass of the Galaxy, but also by prior information on the local disk surface density in visible matter. Therefore, for the fitting of the disk model to the observational constraints described in Sect. 4, we take the values of the scale lengths a_{i} and scale heights b of the MNdisks as fixed and equal to those given in Table 3. The total mass of the disk is adjusted by introducing the parameter f_{mass}, a scale factor by which the disk total mass is then multiplied. Once the best value for f_{mass} is found, then the masses M_{i} of all the MNdisks in Table 3 have to be multiplied by this factor. In this way, in the fitting process of the disk model based on the dynamical constraints, only the disk total mass is treated as a free parameter. Although the individual masses M_{i} of the MNdisks have already been found after the construction of the disk mass model in Sect. 2, in the next step we use their single values in Table 3 as first guesses for the optimal fit values to the dynamical constraints, but keeping their relative contributions to the disk total mass unaltered since all of them are rescaled by the same factor f_{mass}.
3.3. An alternative model: a disk with a ring density structure near the solar orbit radius
There are reasons to believe that the Galactic disk presents a local minimum of density associated with the corotation resonance radius of the Galactic spiral structure. Such a minimum density would be present in both the stellar and gaseous disks components. A structure like this has been observed in the hydrogen distribution of the disk, as discussed later in this section. A minimum in the density of Cepheids is seen at the same radius (BLJ, see their Fig. 14). However, as they are relatively young, the Cepheids possibly trace the gas distribution and therefore the recent past of star formation. There are theoretical argumentation and simulations given in BLJ that predict that the corotation resonance scatters the stars out from that radius, situated close to the solar radius, and produces a minimum in the stellar density. However, there is no direct strong evidence from photometric studies of the disk for such a ring of minimum density of stars. For this reason, our proposed ring density structure is still a speculative one, which we attempt to test in this work.
In Zhang (1996), based on Nbody simulation results of a spiral galaxy evolution and also on theoretical predictions about secular processes of energy and angular momentum transfer between stars and the spiral density wave, the author showed that a local minimum in the stellar density of the disk is formed, centred at the corotation radius. BLJ performed numerical integrations of testparticle orbits with a representative model for the potential due to the axisymmetric density distribution in the Galactic disk and models for the gravitational perturbation due to the spiral arms. The authors verified that a minimum of stellar density with relative amplitude ~30%–40% of the background density is formed at the corotation radius in a time interval ~3 billion years of the evolution of the system. Lépine et al. (2001), based on a simulation of gascloud dynamics in the spiral gravitational field model of the Galaxy, also verified a gap in the ISM distribution at the corotation radius. The authors compared their results (cf. their Fig. 4) with the observed H I radial distribution presented by Burton (1976, see his Fig. 6) and found close similarities between them. Such results are compatible with evidence for a ring that is void of gas, as observed by Amôres et al. (2009) as gaps in the H I density distributed in a ringlike structure with radius slightly outside the solar circle. Recent studies on the threedimensional distribution of the neutral and molecular hydrogen in the Galactic disk by Nakanishi & Sofue (2016) and Sofue & Nakanishi (2016) also corroborate the existence of the minimum gas density slightly beyond the solar Galactic radius. Since the corotation resonance is observed to be at a slightly larger but very close radius to the solar orbit radius (e.g. Dias & Lépine 2005), then the solar orbit would be placed close to the minimum of the disk surface density.
Lépine et al. (2001) associated the minimum in the stellar and gas density distributions with a local dip in the observed rotation curve of the Galaxy at R ~ 9 kpc. A common association between these two features was also modelled by Sofue et al. (2009). In the simulation results of BLJ, the minimum density at corotation is followed by a local maximum density at a slightly larger radius. BLJ modelled such a density feature as a wavy ring superposed on the surface density profile of the disk following Sofue et al. (2009).
Since we still have no clue of how the threedimensional distribution of this density feature might appear, we make the simplest assumption that it would be mainly detected in the density distribution very near the midplane of the disk. Hence we can take the zerothickness disk approximation for the surface density of the ring. We propose a tentative analytical function for the gravitational potential associated with the ring density (local minimum followed by a maximum density), which is written in the form (22)where A_{ring} is the amplitude of the potential related to the amplitude of the minimum and maximum densities of the ring, given in units of km^{2} s^{2}; β_{ring} is a parameter related to the ring width; R_{ring} is the ring node radius, where the inflection point between the minimum and maximum density features sits; and h_{zring} is the scale height of the ring potential. Solving the Poisson equation for the ring potential ∇^{2}Φ_{ring} = 4πGΣ_{ring}δz (within the zerothickness assumption), we can find the corresponding expression for the ring surface density Σ_{ring}.
In the model that includes the ring structure, the surface density Σ_{ring} is added to the disk surface density Σ_{d}, but since the ring is modelled by a minimum followed by a maximum density of similar amplitudes, which makes the net ring mass M_{ring} ~ 0, the disk total mass is then kept approximately unaltered. With the ring density structure, the Galactic models are adjusted with the inclusion of the ring free parameters: A_{ring}, β_{ring}, and R_{ring}; the scale height h_{zring} is chosen to be fixed at a predetermined value.
The inclusion of the ring density feature in the Galactic models is translated in a rescaling of the disk total mass towards larger values when compared to the disk models without the ring. This result is better discussed in Sect. 6.3.
4. Observational constraints
In this section we present the groups of observational data used to constrain the free parameters of the Galactic models introduced in Sect. 3. These groups basically comprise kinematic data from the rotation curve of the Galaxy, with tangent velocities at radii R<R_{0} and rotation velocities at R>R_{0}, the local angular rotation velocity Ω_{0}, and values for the estimated total local disk masssurface density Σ_{0d} already discussed in Sect. 2, as well as the surface density integrated within 1.1 kpc of the disk midplane. In the following, we discuss each group of constraints.
4.1. The local angular rotation velocity and the solar kinematics
In this work, we take the Galactocentric distance of the Sun R_{0} from the statistical analysis performed by Malkin (2013) on 53 R_{0} measurements published in the literature over the last 20 yrs, whose average value, for practical purposes, is recommended by the author as written as For the peculiar velocity of the Sun v_{⊙} relative to the local standard of rest (LSR), we take the reevaluation proposed by Schönrich et al. (2010, hereafter SBD), where we use a lefthanded system for (U, V, W), in which U is positive towards the Galactic anticentre, V is positive in the direction of Galactic rotation, and W is positive towards the direction of the North Galactic Pole; (u_{⊙},v_{⊙},w_{⊙}) are the Sun’s velocity components in such system. The abovequoted uncertainties are the systematic ones, since these dominate the total uncertainties as estimated by SBD.
To constrain the angular velocity of the Sun Ω_{⊙}, we consider the direct measurement of Sgr A^{∗} proper motion along the Galactic plane carried out by Reid & Brunthaler (2004), whose value is μ_{Sgr A∗} = 6.379 ± 0.024 mas yr^{1}. Thus we have Ω_{⊙} = μ_{Sgr A∗} = 30.24 ± 0.11 km s^{1} kpc^{1}. This last equality comes from the assumption that Sgr A^{∗} is at rest at the Galactic centre, so its apparent proper motion can be thought to be solely due to the sum of the Galactic rotation at the LSR and the solar peculiar motion with respect to the LSR in the same direction, i.e. Ω_{⊙} ≡ Ω_{0} + v_{⊙}/R_{0} (e.g. Honma et al. 2012). We then have (23)for the local angular rotation velocity Ω_{0}. With the abovequoted values for Ω_{⊙}, v_{⊙}, and R_{0}, one obtains Ω_{0} = 28.7 km s^{1} kpc^{1}. The uncertainty on Ω_{0} is calculated as , which returns the value σ_{Ω0} = 0.4 km s^{1} kpc^{1}; to the uncertainty σ_{v⊙} of 2 km s^{1} given by SBD, we added 1 km s^{1} to allow for a possible peculiar motion of Sgr A^{∗} at the Galactic centre (McMillan & Binney 2010). The corresponding rotation velocity at the LSR, V_{0} = Ω_{0}R_{0}, assumes the value V_{0} = 230 ± 8 km s^{1}.
4.2. The rotation curve
4.2.1. Tangent velocities
The tangent (or terminal) velocity V_{term} is usually assumed as the maximum velocity of the ISM gas along a given line of sight at Galactic coordinates b = 0° and −90° ≤ l ≤ 90°, which can be related to the circular velocity V_{c}(R) at the tangent point R = R_{0}sinl (or subcentral point), considering a circularly rotating gas. For the tangent velocity data, we use the COline data inside the solar circle from Table 2 of Clemens (1985)^{2}, which cover longitudes in the first Galactic quadrant. We also use the H I tangent point data from Table 2 of Fich et al. (1989), which cover both the first and fourth Galactic quadrants. We converted the LSR tangent velocities of these compiled data to heliocentric velocities and then back to LSR tangent velocities using the components of the peculiar solar motion adopted in this work. Then the Galactocentric distances R = R_{0}sinl and the rotation velocities V_{rot} = V_{term} + V_{0}sinl were calculated using the Galactic constants adopted in this work (R_{0}, V_{0}) = (8 kpc, 230 km s^{1}). We propagated the uncertainties on both R_{0} and V_{0} to the uncertainties on R and V_{rot} (σ_{R} and σ_{Vrot}), respectively.
It has often been argued that the central region of the Galaxy is strongly affected by nonaxisymmetric structures such as the bar, which can induce noncircular motions of the ISM. The true Galactic rotation curve would then be distorted by the measurement of nonuniform azimuthal velocities, making the tangentpoint method inappropriate at these regions. Recently, Chemin et al. (2015) attempted to quantify the asymmetries in the Galactic rotation curve derived by the tangentpoint method. Figure 3 of Chemin et al. (2015) shows that the largest discrepancies between the first and fourth quadrant rotation curves are in the interval 1 ≲ R ≲ 2 kpc. Perhaps the most delicate feature in the Milky Way inner rotation curve, as derived from the tangentpoint method, is the velocity peak at R ~ 300 pc. Based on a numerical simulation of a disk galaxy similar to the Milky Way, Chemin et al. (2015) concluded that the tangent velocities lead to an inner velocity profile with a peak that is not present in the true rotation curve, when the bar major axis is viewed with angles <45° with respect to the direction of the galaxy’s centre. However, we argue that the interpretation of this inner peak as due to the contribution of the bulge can still be defensible: at R ~ 300 pc, a centrally concentrated bulge might dominate the rotation curve, while the bar, as a more extended structure, might influence the rotation curve at radii larger than that.
Also based on the simulated galaxy, Chemin et al. (2015) show that the resulting velocity profile strongly deviates from the true rotation curve in the region R< 4 kpc; the tangentpoint method in the inner regions systematically select highvelocity gas along the bar and spiral arms, or lowvelocity gas in the less dense media. The authors state that the observed rotation curve of the Milky Way derived by the tangentpoint method is expected to be close to the true rotation curve only for radii 4 ≲ R ≤ 8 kpc. The unreliable rotation curve at 2 ≲ R ≲ 4 kpc in the Chemin et al. simulation is attributed to effects associated with the corotation of the bar, which have a pattern speed of the order of 59 km s^{1} kpc^{1} in their model. Other authors have obtained lower values for the bar pattern speed, for example 30–40 km s^{1} kpc^{1} (RodriguezFernandez & Combes 2008) or 33 km s^{1} kpc^{1} (Li et al. 2016). These lower patterns would put the bar corotation at radii larger than 4 kpc. The tangentpoint data from Fich et al. (1989) show a velocity difference on the two sides of the Galactic centre, which is of the order of 12 km s^{1} on average, in the range R = 2–4 kpc. We think that this amplitude of asymmetry, contrary to that observed in the region 1 <R< 2 kpc, does not influence our models. Given the considerations presented above, we decided to restrict the tangent velocity data to  sinl  ≥ 0.3 (e.g. Dehnen & Binney 1998), which in our case is equivalent to Galactic radii R ≥ 2.4 kpc. We end up with a data set that totalizes 280 rotation velocity measurements in the inner solar circle region (2.4 ≤ R ≤ 8.0 kpc).
4.2.2. Maser sources data
From very long baseline interferometry techniques, several maser sources associated with highmass starforming regions (HMSFRs) have been studied and their positions, parallaxes, and proper motions have been measured with high accuracy. Complementing these data with heliocentric radial velocities from Doppler shifts, we are able to access the full threedimensional location of each source in the Galaxy as well as their full space motion relative to the Sun. Since it is believed that the maser sources do not present large peculiar motions, we can use their velocity components in the direction of Galactic rotation as a proxy for the rotation curve of the Galaxy (Irrgang et al. 2013), taking for this the value of v_{⊙} estimated by SBD. The data for HMSFRs with maser emission were obtained from Table 1 of Reid et al. (2014), where the authors list parallaxes, proper motions, and LSR radial velocities of 103 regions measured with Very Long Baseline Interferometry (VLBI) techniques from different surveys and projects (the Bar and Spiral Structure Legacy (BeSSeL) Survey^{3} and the Japanese VLBI Exploration of Radio Astrometry (VERA)^{4}). We converted the tabulated LSR radial velocities to heliocentric radial velocities by adding back the components of the standard solar motion (Reid et al. 2009). With the coordinates, parallaxes, proper motions, heliocentric radial velocities, and their respective errors, we calculated the heliocentric U, V, W velocities and uncertainties σ_{U}, σ_{V}, and σ_{W} for each source, following the formalism described by Johnson & Soderblom (1987). Correcting for the SBD solar peculiar motion and for the LSR circular velocity V_{0}, we calculated the Galactocentric component V_{φ} of the space velocity of each maser source in the direction of Galactic rotation. The uncertainties on V_{φ} and σ_{Vφ} were obtained by propagation from the uncertainties on the parallaxes, proper motions, heliocentric radial velocities, and the uncertainties of the solar motion. The rotation velocities and uncertainties were then defined by V_{rot} = V_{φ} and σ_{Vrot} = σ_{Vφ}. The Galactic radii were obtained directly from the positions and parallaxes of the sources as well as their uncertainties. The distribution of Galactic radii comprises regions both inside and outside the solar circle. We discarded the sources with radii R< 4 kpc because most of them present values of V_{φ} with large deviations from the rotation curve that are traced by the tangent velocity data at these radii; a similar selection criterion was used by Reid et al. (2014). This selection reduces the number of maser sources used to probe the rotation curve to 94 objects.
Although the rotation curve for R>R_{0} that is traced by the maser sources is restricted to a few data points, these measurements are the best that we can think of, considering the precision in the distances, proper motions, and lineofsight velocities. The less accurate distances of H II regions, for instance, can produce false trends in the rotation curve traced by these objects, which is the main reason for not using them in the present study. Indeed, Binney & Dehnen (1997) pointed out that the apparent rising rotation curve traced by the H II regions outside the solar radius can be explained by the objects tending to be concentrated in a “ring” (which could be a segment of a spiral arm) with a mean radius larger than their estimated Galactic radii. We believe that these shortcomings disappear when we take the outer rotation curve traced by the maser sources data with welldetermined distances. A more quantitative analysis on the correlations between the velocity uncertainties and distance uncertainties of the maser sources is presented in Sect. 6.3.
4.3. Local surface densities
As presented in Sect. 2.1, our “observationbased” disk model is constructed to give a total local disk masssurface density of Σ_{0d} = 57.2M_{⊙} pc^{2}, which is based on recent determinations in the literature about the local surface densities in both stellar and gaseous components. This is also the value that we adopt as constraint to the local dynamical disk surface density, which can be compared to the estimate by Holmberg & Flynn (2004) of 56 ± 6M_{⊙} pc^{2} for this quantity. We adopt the same uncertainty on Σ_{0d} of σ_{Σ0d} = 6M_{⊙} pc^{2}. Holmberg & Flynn (2004) also estimated the local surface density integrated within 1.1 kpc of the disk midplane as being Σ_{01.1 kpc} = 74 ± 6M_{⊙} pc^{2}, which, according to the authors, takes into account both disk and dark halo contributions. We take such a value as an observational constraint on Σ_{01.1 kpc} and its uncertainty.
5. Fitting procedure
As exposed in Sect. 3, we attempt to construct two types of Galactic models: those that do not incorporate the ring density structure in the disk (Sect. 3.3) and that we refer to as model MI; and those that incorporate the ring density structure and are denoted by model MII. The free parameters of model MI are the bulge parameters, M_{b} and a_{b}; the dark halo parameters, r_{h} and v_{h}; and the disk mass scale factor, f_{mass}. The free parameters of model MII are the same as model MI, plus the ring parameters A_{ring}, β_{ring}, and R_{ring}; the scale height h_{zring} is chosen to be fixed at the value 0.65 kpc.
We search for the bestfit set of parameters for both models MI and MII using a χ^{2}minimization procedure, which is implemented through the crossentropy (CE) algorithm (Sect. 2.2.2). We first attempt initial guesses for the parameters sets {M_{b},a_{b},r_{h},v_{h},f_{mass}} of model MI and of model MII based on visual fits of the given model to the observed rotation curve. From these trial sets of parameters, we randomly generate N = 3 × 10^{3} independent sets of model parameters by considering initial uniform distributions centred on the given trial parameters and with halfwidths equal to the initial uncertainties chosen for each parameter. We select the best 100 sets candidates that compose the elite sample of sets based on the χ^{2}minimization criterion (the N_{elite} = 100 sets that return the lowest values for χ^{2}). From the mean μ and standard deviation σ of each ensemble of 100 parameters of this elite sample, we generate N other independent sets of trial parameters through Gaussian distributions N(μ,σ^{2}), which will be evaluated in the next iteration. This process is repeated by several iterations until we obtain a stable set of parameters for each model MI and MII.
The total gravitational potential Φ is calculated as the sum of the gravitational potentials of each individual Galactic component: bulge, disk, and dark halo, Φ(R,z) = Φ_{b}(R,z) + Φ_{d}(R,z) + Φ_{h}(R,z), in the case of model MI, and the addition of the ring potential in the case of model MII, Φ(R,z) = Φ_{b}(R,z) + Φ_{d}(R,z) + Φ_{h}(R,z) + Φ_{ring}(R,z). Considering the balance between the centrifugal force and gravity, the circular velocity V_{c}(R) of a given model, measured for instance in the Galactic midplane z = 0, is linked to the total gravitational potential by the form (24)Therefore, the gravitational potential Φ and the model rotation curve V_{c}(R) are totally defined by the set of parameters of each model MI and MII described above. The radial derivative of the potential dΦ(R,0)/dR in Eq. (24) is given by the sum of the radial derivatives of the potentials of each Galactic component, for which explicit expressions are presented in Appendix B. Each data point i of the rotation curve has a pair of values (R_{i}, V_{roti}). In the fitting process of the rotation curve, we search for the minimization of the residuals between the rotation velocities V_{roti} of the observational data and the circular velocities V_{c}(R_{i}) that resulted from the model at each radius R_{i}. We thus minimize the quantity (25)where we denote as the chisquared component relative to the rotation curve observational constraint. We separate the rotation curve data into four different groups: group 1 corresponds to the CO tangent velocity data; group 2 corresponds to the H I tangent velocities in the first quadrant (l> 0°); group 3 is relative to the H I tangent velocities in the fourth quadrant (l< 0°); and group 4 comprises the maser sources data. Then the in Eq. (25) is calculated for each group separately, where the subindex k denotes the number of the group.
As the data points of the rotation curve present errors in both R and V_{rot} (σ_{R}, σ_{Vrot}), a proper fitting procedure has to take into account the effects of both groups of errors on the model to be adjusted. We deal with this issue by following the solution adopted by Irrgang et al. (2013); the uncertainty σ_{R} is converted to an error in V_{rot} by estimating its effect on the model rotation curve according to the relation σ_{Vc} = (dV_{c}/ dR) σ_{R}, and then σ_{Vc} is added in quadrature to σ_{Vrot}. The total velocity uncertainty σ_{Vtot} in Eq. (25) is just given by .
Regarding the other observational constraints discussed in Sect. 4, namely, the local angular rotation velocity Ω_{0}, the local disk masssurface density Σ_{0d} and the surface density within  z  ≤ 1.1 kpc Σ_{01.1 kpc}, they contribute to the total χ^{2} in the form (26)where ψ_{j, obs} and σ_{ψj, obs} refer to each one of the three j abovementioned observables and their associated uncertainties, respectively; and ψ_{j, model} refers to the respective quantities resulted from the models. The total weighted χ^{2} is the sum (27)where N_{rck} is the number of data points in each group k of rotation curve data, and N_{other} = 3. Each contribution of the groups of observational constraints to the total χ^{2} is divided by the number of observational data actually used in the group. This is the same procedure adopted by Dehnen & Binney (1998) and Irrgang et al. (2013) to ensure that the fitting process is not dominated by the rotation curve owing to the larger number of individual data points in this group of observational constraints.
Uncertainty estimates.
We apply Monte Carlo techniques to obtain the uncertainties on the fitting parameters of the Galactic models. From the original observational data set compiled for the rotation curve, namely the Galactic radii R and rotation velocities V_{rot} of the sources, we create N_{run} = 100 new data sets by resampling the original one with replacement of the data to perform a bootstraplike procedure (e.g. Monteiro et al. 2010). We then run the CE algorithm that implements the fitting procedure described above N_{run} times to find the bestfitting parameter sets, each time with one of the “new rotation curves” resampled in the way explained above. We end up with 100 sets of parameters for each model, MI and MII, from which the uncertainties on the parameters are determined. It was found that the average values of the parameters obtained from these 100 parameter sets are very close to the bestfit values found after the fitting of the models using the original rotation curve data. However, we prefer to consider these latter values the best ones since they were directly obtained from the original data set.
The uncertainties on the scale lengths a, scale heights b, and masses M of the MiyamotoNagai disks depend directly on the uncertainties on the structural parameters and local surface densities adopted for the “observationbased” disk model. Therefore, a proper estimation of such uncertainties should take into consideration the combination of all errors of the parameters used to constrain the MNdisk models. A rough calculation, using Monte Carlo techniques, indicates that uncertainties on the a_{i} and b_{i} parameters (σ_{ai}, σ_{bi}) of ~10% of their values listed in Table 3 seem to be appropriate estimates. The uncertainties on the masses of the MNdisks are considered proportional to the uncertainties on the mass scale factor σ_{fmass}, i.e. σ_{Mi} = M_{i}σ_{fmass}, with M_{i} also given in Table 3.
Bestfitting values for the parameters of the Galactic models MI and MII.
Local properties resulting from the models MI and MII for the mass distribution and gravitational potential of the Galaxy.
Correlation matrix for the model parameters: the lowerleft triangle for model MI and the upperright triangle for model MII.
Fig. 5
Lefthand column, top panel: rotation curve resulting from the Galactic model MI. The circular velocities due to the three mass components are depicted as dotted lines for the bulge, dashed lines for the disk, and dashdotted lines for the dark halo. The total circular velocity is represented by the solid curve. The data for the observed rotation curve are presented as red points with orange error bars for the CO tangentpoint data (Clemens 1985) and the H I tangentpoint data (Fich et al. 1989); and blue points with light blue error bars for the maser sources data (Reid et al. 2014). Bottom panel: velocity residuals between the observed rotation velocities and the total circular model velocities calculated at each radius of the data points. Righthand column: as for the lefthand column but for the Galactic model MII. 
6. Results and discussions
The bestfitting values for the parameters of both Galactic models MI and MII are summarized in Table 4. Table 5 lists the values for the local dynamical properties resulting from the models, some of which can be compared to the values used as observational constraints and, given the uncertainties, a good overall agreement is observed. Figure 5 shows, in the top panels, the rotation curves resulting from each Galactic model (solid lines) and the curves relative to the contribution of each Galactic component. The data for the observed rotation curve are presented as red points with orange error bars for the CO tangentpoint data (Clemens 1985) and for the H I tangentpoint data (Fich et al. 1989) and blue points with light blue error bars for the maser sources data (Reid et al. 2014). The bottom panels present the residuals between the observed and the modelled rotation velocities at each radius of the data points. Figure 6 shows the tangentpoint data in the plane of observables V_{term} − sinl and the curves resultant from the models MI and MII calculated as V_{term} = V_{rot}(R) − V_{0}sinl, with sinl = R/R_{0}. A good agreement between the two curve models and the data can be observed.
In order to compare the goodness of fit between models MI and MII, we use the chisquared per degree of freedom statistics, χ^{2}/ d.o.f.. Model MI returns χ^{2}/ d.o.f. = 1.85, while model MII provides a fit with χ^{2}/ d.o.f. = 1.46. The relatively better fit provided by model MII can be explained, in part, by its better match to the dip in the observed rotation curve at radii in the interval ~8–10 kpc. As mentioned before, the correspondent dip in the model rotation curve is a consequence of the ring density structure added to the disk model density profile.
In Table 6 we give the correlation matrix of the fitted parameters for both models MI (the lowerleft triangle) and MII (the upperright triangle). In a given correlation matrix, the elements are distributed in the interval [−1; 1], where the value of 1 indicates a perfect correlation and −1 indicates a perfect anticorrelation, whereas 0 corresponds to no correlation. As can be seen from Table 6, the strongest correlations are between parameters associated with a single component, as is the case of the bulge parameters M_{b} and a_{b}, and also the halo parameters r_{h} and v_{h}. An anticorrelation between the asymptotic circular velocity of the dark halo v_{h} and the disk mass scale factor f_{mass} is observed from model MI; since the dark halo mass must be proportional to v_{h}, this anticorrelation means there is a competition between the disk and the dark halo for the contribution to the dynamical mass in the inner Galaxy. For model MII, anticorrelations are also seen between M_{b} and f_{mass}, and between f_{mass} and the ring parameter R_{ring}.
6.1. Comparison with studies in the literature
The masses attributed to the bulge in our models are relatively large (M_{b} ~ 2.6 × 10^{10}M_{⊙}) compared to other estimates in the literature. This is one cause for the small contribution of the disk component of model MI to the rotational support in the inner Galaxy and also, but to a lesser extent, in model MII (see next section). Most of the studies in the literature find bulge masses of the order of 1 × 10^{10}M_{⊙}. However, these studies also adopt simple exponential profiles for the disk component with masses in the central regions that are larger than those obtained with our disk models with central density depletion. As pointed out by McMillan (2011), his bestfitting Galactic model presents a stellar mass within the inner 3 kpc of ~2.4 × 10^{10}M_{⊙} (the same result obtained in the Flynn et al. 2006 model), which is in close agreement with the bulge mass found by Picaud & Robin (2004) assuming a disk model with a central hole. This is also the same value for the mass of the bulge in the model of Lépine & Leroy (2000), which is also based on a disk with a central density depletion. As described in Sect. 3.1, we based the initial guess values for the bulge parameters on the spheroidal component of the PJL model, which returns a mass inside 3 kpc of ≈2.2 × 10^{10}M_{⊙} with a scale radius parameter of 0.4 kpc. The values for the bulge masses M_{b} and scale radii a_{b} found in our models (see Table 4) reflect the consequence of this choice. The relatively short scale radii cause the rise of the peak in the rotation curves within the inner 1 kpc, as can be seen from Fig. 5. Therefore, according to the PJL star counts model, and also reflected in our models, the bulge is an important mass component in the central region of the Galaxy and its contribution to the model rotation curve is dominant at these radii. Our models return stellar masses enclosed within the inner 3 kpc of 2.7 × 10^{10}M_{⊙} (model MI, that is 1.98 × 10^{10}M_{⊙} owing to the bulge and 0.72 × 10^{10}M_{⊙} owing to the stellar thin and thick disks), and 3 × 10^{10}M_{⊙} (model MII, that is 1.98 × 10^{10}M_{⊙} owing to the bulge and 1.02 × 10^{10}M_{⊙} owing to the stellar thin and thick disks). We can therefore argue about the tendency of the disk models with central holes to redistribute the mass in the central regions from the disk to the bulge without significantly altering the total stellar mass in such regions.
Fig. 6
Tangentpoint data from Clemens (1985; green circles) and Fich et al.1989; orange circles). The tangent velocity curves resultant from model MI (blue curve) and model MII (red curve) are also plotted. 
The resulting disk stellar mass (thin disk + thick disk) is 3.2 × 10^{10}M_{⊙} from model MI and 4.74 × 10^{10}M_{⊙} from model MII. As a comparison, Flynn et al. (2006, and references therein) argue that models with disk scale lengths in the range 2–2.6 kpc, based on studies of the Galactic emission in the nearinfrared or studies of the local stellar distribution, correspond to disk stellar masses in the range 3.6–5.4 × 10^{10}M_{⊙}. The total stellar mass (bulge + thin disk + thick disk) of 5.81 × 10^{10}M_{⊙} from model MI is close to the interval of 4.85–5.5 × 10^{10}M_{⊙} estimated by Flynn et al. (2006), and the total stellar mass of 7.37 × 10^{10}M_{⊙} from model MII is relatively close to the bestfitting value of 6.61 × 10^{10}M_{⊙} in the model by McMillan (2011). The disk mass in the gaseous form, which is already corrected for the mass contribution of helium, is 1.42 × 10^{10}M_{⊙} from model MI and 2.1 × 10^{10}M_{⊙} from model MII. The total baryonic mass (stellar + gas) from the models are 7.23 × 10^{10}M_{⊙} (model MI) and 9.47 × 10^{10}M_{⊙} (model MII). These values are somewhat larger than the “back of the envelope” estimate by Flynn et al. (2006) of 6.1 ± 0.5 × 10^{10}M_{⊙}, but are compatible with the masses (bulge + disk) of 7–8.2 × 10^{10}M_{⊙} found in the models by Irrgang et al. (2013), for instance.
Table 5 also lists the values of the Oort’s constants A and B resulting from the models. Although these constants were not used as observational constraints in the modelling process, their returned values from model MI agree with those estimated in the literature (e.g. Feast & Whitelock 1997). The steeper gradient of the rotation curve at the solar Galactic radius R_{0} causes the deviation of the Oort’s constants of model MII from their common estimated values. The local escape velocities of 452 km s^{1} from model MI and 550 km s^{1} from model MII are close and within the 90% confidence interval of 492–587 km s^{1} as determined by Piffl et al. (2014b).
6.2. The disk support to the rotation curve
The disk mass scale factor f_{mass} of model MI implies no increase of the originally modelled disk mass. Although returning local surface densities that are compatible with the observed densities that are used as constraints, within the adopted uncertainties, the model MI returns a dynamical mass for the Galactic disk that moderately contributes to the total rotation curve in the inner Galaxy. Sackett (1997), based on a set of Galactic observational constraints, verified that the “maximal disk hypothesis” that is commonly applied to external spiral galaxies also gives a maximal disk when applied to the Milky Way. According to this definition, to be maximal, the exponential disk must provide 85% ± 10% of the total rotation velocity of the galaxy at the radius R = 2.2R_{d}, where the rotation curve of the disk presents a peak (R_{d} is the exponential disk scale length). Since our disk models are no longer simple exponentials, we just compare the peak rotation speed of the disk with the total rotation speed at the same radius. From Fig. 5, the peak in the rotation curve of the disk occurs at R ~ 7 kpc for both models MI and MII. We estimate that the disk of model MI is responsible for only 64% of the total circular velocity V_{c}(R = 7 kpc), which puts it in the condition of a “submaximal” disk, however.
On the other hand, the mass scale factor of model MII increases the disk mass by 49% of its original value. This increase in the disk mass is a consequence of the addition of the ring density structure to the density profile of the disk (see Sect. 6.3). Once the solar orbit radius is sitting inside the valley of the ring density feature (near the minimum density, see Fig. 8), to maintain the local disk surface density value within the range determined by the observations, the underlying surface densities at the other radii are proportionally increased. For instance, if it were not for the presence of the ring density giving a local disk surface density of Σ_{0d} = 61M_{⊙} pc^{2} (cf. Table 5), the disk mass of model MII would return Σ_{0d} = 90M_{⊙} pc^{2} instead. We estimate that the disk of model MII contributes with 78% of the total circular velocity at R = 7 kpc, putting it within the limit for a “maximal disk” condition. Therefore, compared with model MI, the disk of model MII contributes more expressively to the inner rotation curve of the Galaxy, giving rise to a less important dark halo component at such regions. It is still important to note that the conclusion by Sackett (1997) about the Galaxy supporting a maximal disk was based, among several observational constraints, on a local disk surface density with a value similar to that used in the present work. But Sackett adopted a local circular velocity of V_{0} = 210 ± 25 km s^{1}, while our adopted V_{0} of 230 km s^{1} sits closer to the high side of such interval. Considering now the total rotation supplied by the mass from both bulge and disk, ~82% and 93% of the total circular velocity at R = 7 kpc are provided by the baryonic matter from models MI and MII, respectively.
6.3. The ring structure
In Sect. 3.3, we give some arguments about our reasons for considering an alternative disk model that incorporates a ring density structure near the solar Galactic radius. In terms of the observational data treated in this work, we argue that the apparent velocity dip in the rotation curve of the maser sources, in the radial interval of ~7–11 kpc (compared to a flat rotation curve V_{c}(R) = V_{0}, for example), can be associated with a ringlike density structure of the type that we consider here. Before modelling such a ring structure, however, it is recommended that we check the reliability of the abovementioned velocity dip. Indeed, correlations between distance uncertainties and velocity uncertainties could create false trends in the plot of the rotation curve. Since the velocity dip of the rotation curve is traced by the maser sources, we restrict ourselves to these data in the following analysis. For each maser source, we decompose the uncertainty on the azimuthal velocity V_{φ}, σ_{Vφ}, into two parts: one component is formed by the combination of the uncertainty on the heliocentric radial velocity, the uncertainties on the two components of the proper motion, and the uncertainties on the solar motion; the second component is solely due to the uncertainty on the distance of the source σ_{d} that is derived from the uncertainty on its parallax. We express this second component as σ_{Vφ}(σ_{d}). To verify how much the correlations between the uncertainties on the velocities V_{φ} and the uncertainties on the distances d affect the velocity dip of the rotation curve, we compute, for each maser source, the difference in velocities (V_{φ} + σ_{Vφ}(σ_{d})) − (V_{φ} − σ_{Vφ}(σ_{d})) = 2σ_{Vφ}(σ_{d}), and compare this difference to the magnitude of the velocity dip V_{dip}. From the rotation curve in Fig. 5, we estimate V_{dip} ≈ 12 km s^{1}. The comparison between 2σ_{Vφ}(σ_{d}) and V_{dip} for the maser sources in the interval 7 kpc ≤ R ≤ 11 kpc is shown in Fig. 7. About 80% of the sources present 2σ_{Vφ}(σ_{d}) ≤ 5 km s^{1}, which is much smaller than V_{dip}. This result seems to indicate that the correlations between distance uncertainties and velocity uncertainties are not able to create a false trend in the rotation curve, which might appear as a dip in the velocity profile. In other words, the velocity dip in the radial interval 7−11 kpc seems to be a real feature of the observed rotation curve.
Fig. 7
Comparison of the rotation velocity difference 2σ_{Vφ}(σ_{d}) (black circles; see text for details) with the velocity dip V_{dip} (horizontal line) for the maser sources in the radial interval comprised by the dip of the rotation curve. 
Fig. 8
Surface density radial profile for the disks of model MI (dashed curve) and model MII (red solid curve). The black solid curve represents a disk with mass equivalent to that from model MII but without the ring density structure. The shaded area in light grey emphasizes the resulting difference in mass between the disks from both models. The blue point denotes the pair (R_{0}; Σ_{0d}). 
The ring density structure associated with the velocity dip in the rotation curve of model MII, whose parameters (A_{ring}, β_{ring}, R_{ring}) are given in Table 4, is depicted in Fig. 8 along with the radial profile for the surface density of the disk from the same model (red curve). The surface density of a disk with equivalent mass to that from model MII, but without the ring structure, is shown by the curve in the black solid line. The ring is formed by a minimum density at the radius of 8.3 kpc, a maximum at 9.8 kpc, and the node radius R_{ring} = 8.9 kpc, which is also the radius of the minimum point in the dip of the modelled rotation curve. The amplitude of the minimum is ≈−0.34 times the density at the same radius of the equivalent disk without the ring. The blue point in Fig. 8 indicates the approximately common value of Σ_{0d} returned from both models MI (indicated by the curve in dashed line) and model MII (red curve) at the solar radius R_{0}. As anticipated in Sects. 3.3 and 6.2, the ring structure causes the rescale of the disk total mass from model MII to a value as large as ~48% of the disk mass from model MI. This increase in mass can be checked out by observing the area highlighted in light grey between the Σ(R) curves of model MI and the equivalent disk of model MII in Fig. 8. The above result is dependent on the scale length measured for the disk of Model MII. The ring structure does not alter the global scale length of the disk, i.e. the global disk scale length of Model MII is the same of model MI. Then the abovequoted increase in the disk mass is satisfied. However, the ring tends to decrease the scale length of the Σ(R) distribution mainly in the radial range of its influence, that is 7 ≲ R ≲ 11 kpc. Once the disk scale length is estimated using this shorter interval of radius, the above considerations of mass increase might not be appropriate.
The density bump at R ~ 9.8 kpc induced by the ring (see Fig. 8) can also be a matter of debate. In fact, the amplitude of this local maximum of density may be too exaggerated, which is a consequence of the form chosen for the ring potential in Eq. (22). Regardless of the true amplitude of the density bump, the question about its plausibility in the observational point of view still needs to be answered. We argue here that since the ring structure was modelled over both stellar and gaseous components, part of such density bump must be associated with the hydrogen distribution of the disk. And indeed, recent works by Nakanishi & Sofue (2016) and Sofue & Nakanishi (2016) clearly show an increase in the hydrogen density of the disk at radii slightly beyond 10 kpc. A stellar counterpart of the density bump, if it existed, would be more difficult to notice for various reasons, such as the relative position of the Sun in the ring structure and the beginning of the flare of the stellar disk just after the solar radius. We believe that estimates of the stellar surface density in a great coverage of the Galactic disk exterior to the solar circle (not only in radius but also in azimuth and height from the plane) are still needed to confirm or refute the existence of the stellar ring structure. In a positive case, our ring model could serve as a starting point to the study of the properties of such a structure in the Galactic disk. These are some aspects that make model MII an “alternative model”, in the sense treated in this work.
7. Conclusions
We have developed models for the axisymmetric mass distribution of the Galaxy with the aim to derive fully analytical descriptions of its associated threedimensional gravitational potential. We followed an approach that intentionally expends more effort in a detailed modelling of the disk component. Based on photometric constraints for the stellar distribution in the disk given by the Galactic infrared star counts model of PJL and on the observed distribution of atomic and molecular hydrogen gas, we derived an “empirical basis” for the structural parameters (scale length, scale height, and radial scale of the disk central hole) of the thin and thick stellar disks and the H I and H_{2} disks subcomponents. With a priori values for the masses of each disk, based on the most recent determinations of the local stellar and gaseous disk surface densities, we constructed versions of MiyamotoNagai disks for each disk subcomponent mass model. The method follows the approach developed by Smith et al. (2015), but we allow ourselves to use the MiyamotoNagai disk models of higher orders 2 and 3 beside the commonly used model of order 1 (Miyamoto & Nagai 1975). Along with parametric models for the bulge and an extra unknown spherical mass component which for “convenience” we refer to by the oftenused term dark halo, we searched for the dynamical mass of each Galactic component by fitting the models to the kinematic constraints given by the observed rotation curve and some local Galactic measured properties.
We have shown that a disk model that includes a ring density pattern beyond but very close to the solar orbit radius is able to better reproduce an observed local dip in the Galactic rotation curve centred at R ~ 9.0 kpc; such dip is naturally explained by a ring density structure composed of a minimum followed by a maximum density of similar amplitudes. Furthermore, the model with the ring structure allows a more massive disk to increasingly contribute to the rotational support of the Galaxy inside the solar circle, thereby helping the Milky Way to satisfy the condition to be considered a “maximal disk” galaxy. The model with the ring structure also relies on a numerical study of the stellar disk, where a ring density structure develops as a consequence of interactions between the stars and the galactic spiral arms in resonance at the corotation radius (BLJ). Since we still have no information about the threedimensional structure of this ring density feature, we had to make crude approximations for the vertical profile of its associated gravitational potential. We believe that, with the forthcoming data of the GAIA mission (Perryman et al. 2001), the ring structure in the disk may become an object of examination and, in the case of confirmation of its existence, constraints on its threedimensional distribution in the global disk structure will help us to create more realistic disk models.
The method we applied for the construction of the disk mass model and its associated gravitational potential is very flexible in the sense that for any set of structural parameters, disks of different masses can be generated. We emphasize at this point that the models of Galactic gravitational potential presented in this work are aimed to be used in studies of orbits in the disk that do not extend too far away in Galactic radii and do not reach great heights above the disk midplane. These limitations are imposed by the spatial coverage of the observational data used to constrain the models in the sense that for radii greater than ~2R_{0} (the maximum radius of the data used for the rotation curve) and heights  z  ≳ 3 kpc, we do not guarantee that our models return confident representations of the Galactic potential, but they may provide reasonable representations at least. Also in this respect, no constraint on the mass at large radii is adopted here, as has been made by some studies that derive properties of the dark halo by requiring distant halo stars to be bound to the Milky Way potential. We included a density component associated with a spherical logarithmic potential to explain the observed rotation curve data at radii R_{0} ≲ R ≲ 2R_{0}. As pointed out by Dehnen & Binney (1998), however, instead of a distinct physical component, it is possible that we are measuring the dynamical effects of a disk and/or a bulge in which the masstolight ratio of their content strongly increases from the centre to the Galactic outskirts.
The models for the threedimensional gravitational potential of the Galaxy presented in this work, as they are fully analytical and easy to obtain the associated gravitational forcefield at any point, are suited for fast and accurate calculations of orbits of stellarlike objects belonging to the main populations of the Galactic disk. In a forthcoming study, we aim to present an application of these new potential models to a description of the distributions of orbital parameters for samples of Galactic open clusters and some expected correlations with their chemical abundance patterns.
The tangent velocity data in Table 2 of Clemens (1985) were corrected for 3 km s^{1} line width, and a correction of 7 km s^{1} for the LSR peculiar motion in the azimuthal direction was used in the calculation of the rotation velocities, as proposed by the author in the abovecited paper.
Acknowledgments
We acknowledge the referee whose comments and suggestions have significantly improved the present paper. D.A.B. received financial support for this work from the Brazilian research agency CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) through the PNPD postdoctoral programme.
References
 Allen, C., & Santillan, A. 1991, Rev. Mex. Astron. Astrofis., 22, 255 [NASA ADS] [Google Scholar]
 Amôres, E. B., & Lépine, J. R. D. 2005, AJ, 130, 659 [NASA ADS] [CrossRef] [Google Scholar]
 Amôres, E. B., Lépine, J. R. D., & Mishurov, Y. N. 2009, MNRAS, 400, 1768 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., & Soneira, R. M. 1980, ApJS, 44, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., Flynn, C., & Gould, A. 1992, ApJ, 389, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Barros, D. A., Lépine, J. R. D., & Junqueira, T. C. 2013, MNRAS, 435, 2299 (BLJ) [NASA ADS] [CrossRef] [Google Scholar]
 Bienaymé, O., Soubiran, C., Mishenina, T. V., Kovtyukh, V. V., & Siebert, A. 2006, A&A, 446, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J., & Dehnen, W. 1997, MNRAS, 287, L5 [NASA ADS] [Google Scholar]
 Binney, J., & Merrifield, M. 1998, Galactic Astronomy (Princeton University Press) [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton University Press) [Google Scholar]
 Bovy, J., & Rix, H.W. 2013, ApJ, 779, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., Rix, H.W., & Hogg, D. W. 2012, ApJ, 751, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Burton, W. B. 1976, ARA&A, 14, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Caetano, T. C., Dias, W. S., Lépine, J. R. D., et al. 2015, New Astron., 38, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Caldwell, J. A. R., & Ostriker, J. P. 1981, ApJ, 251, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Caproni, A., Monteiro, H., & Abraham, Z. 2009, MNRAS, 399, 1415 [NASA ADS] [CrossRef] [Google Scholar]
 Chemin, L., Renaud, F., & Soubiran, C. 2015, A&A, 578, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Clemens, D. P. 1985, ApJ, 295, 422 [NASA ADS] [CrossRef] [Google Scholar]
 de Vaucouleurs, G. 1977, AJ, 82, 456 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Dias, W. S., & Lépine, J. R. D. 2005, ApJ, 629, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Dias, W. S., Monteiro, H., Caetano, T. C., et al. 2014, A&A, 564, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feast, M., & Whitelock, P. 1997, MNRAS, 291, 683 [NASA ADS] [CrossRef] [Google Scholar]
 Fich, M., Blitz, L., & Stark, A. A. 1989, ApJ, 342, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Flynn, C., SommerLarsen, J., & Christensen, P. R. 1996, MNRAS, 281, 1027 [NASA ADS] [CrossRef] [Google Scholar]
 Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Freeman, K. C. 1970, ApJ, 160, 811 [NASA ADS] [CrossRef] [Google Scholar]
 Freudenreich, H. T. 1998, ApJ, 492, 495 [NASA ADS] [CrossRef] [Google Scholar]
 Gilks, W., Richardson, S., & Spiegelhalter, D. 1996, Markov chain Monte Carlo in practice (London: Chapman and Hall) [Google Scholar]
 Hernquist, L. 1990, ApJ, 356, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Hessman, F. V. 2015, A&A, 579, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Holmberg, J., & Flynn, C. 2000, MNRAS, 313, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Holmberg, J., & Flynn, C. 2004, MNRAS, 352, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Honma, M., Nagayama, T., Ando, K., et al. 2012, PASJ, 64, 136 [Google Scholar]
 Irrgang, A., Wilcox, B., Tucker, E., & Schiefelbein, L. 2013, A&A, 549, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864 [NASA ADS] [CrossRef] [Google Scholar]
 Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kalberla, P. M. W., & Dedes, L. 2008, A&A, 487, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalberla, P. M. W., Dedes, L., Kerp, J., & Haud, U. 2007, A&A, 469, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalberla, P. M. W., Kerp, J., Dedes, L., & Haud, U. 2014, ApJ, 794, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J. 1977, ApJ, 217, 406 [NASA ADS] [CrossRef] [Google Scholar]
 Kroese, D. P., Porotsky, S., & Rubinstein, R. Y. 2006, Methodology and Computing in Applied Probability, 8, 383 [Google Scholar]
 Kuijken, K., & Gilmore, G. 1991, ApJ, 367, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Kuzmin, C. G. 1956, Astron. Zh, 33, 27 [Google Scholar]
 Lépine, J. R. D., & Leroy, P. 2000, MNRAS, 313, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Lépine, J. R. D., Mishurov, Y. N., & Dedikov, S. Y. 2001, ApJ, 546, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Z., Gerhard, O., Shen, J., Portail, M., & Wegg, C. 2016, ApJ, 824, 13 [NASA ADS] [CrossRef] [Google Scholar]
 LópezCorredoira, M., CabreraLavers, A., Gerhard, O. E., & Garzón, F. 2004, A&A, 421, 953 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lozinskaya, T. A., & Kardashev, N. S. 1963, Soviet Ast., 7, 161 [NASA ADS] [Google Scholar]
 Malkin, Z. 2013, in IAU Symp. 289, ed. R. de Grijs, 406 [Google Scholar]
 Martins, L. P., Coelho, P., Caproni, A., & Vitoriano, R. 2014, MNRAS, 442, 1294 [Google Scholar]
 McMillan, P. J. 2011, MNRAS, 414, 2446 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J., & Binney, J. J. 2010, MNRAS, 402, 934 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
 Monteiro, H., & Dias, W. S. 2011, A&A, 530, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Monteiro, H., Dias, W. S., & Caetano, T. C. 2010, A&A, 516, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nakanishi, H., & Sofue, Y. 2006, PASJ, 58, 847 [NASA ADS] [Google Scholar]
 Nakanishi, H., & Sofue, Y. 2016, PASJ, 68, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Perryman, M. A. C., de Boer, K. S., Gilmore, G., et al. 2001, A&A, 369, 339 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Picaud, S., & Robin, A. C. 2004, A&A, 428, 891 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piffl, T., Binney, J., McMillan, P. J., et al. 2014a, MNRAS, 445, 3133 [NASA ADS] [CrossRef] [Google Scholar]
 Piffl, T., Scannapieco, C., Binney, J., et al. 2014b, A&A, 562, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Polido, P., Jablonski, F., & Lépine, J. R. D. 2013, ApJ, 778, 32 (PJL) [NASA ADS] [CrossRef] [Google Scholar]
 Read, J. I. 2014, J. Phys. G Nucl. Phys., 41, 063101 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Menten, K. M., Zheng, X. W., et al. 2009, ApJ, 700, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [Google Scholar]
 Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RodriguezFernandez, N. J., & Combes, F. 2008, A&A, 489, 115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rohlfs, K., & Kreitschmann, J. 1988, A&A, 201, 51 [NASA ADS] [Google Scholar]
 Rubinstein, R. Y. 1997, Eur. J. Oper. Res., 99, 89 [Google Scholar]
 Rubinstein, R. Y. 1999, Methodology and Computing in Applied Probability, 1, 127 [Google Scholar]
 Sackett, P. D. 1997, ApJ, 483, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, D. B., Solomon, P. M., & Scoville, N. Z. 1984, ApJ, 276, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, M. 1956, Bull. Astron. Inst. Netherlands, 13, 15 [Google Scholar]
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 (SBD) [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 2004, in AIP Conf. Ser. 735, eds. R. Fischer, R. Preuss, & U. V. Toussaint, 395 [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R., Flynn, C., Candlish, G. N., Fellhauer, M., & Gibson, B. K. 2015, MNRAS, 448, 2934 [NASA ADS] [CrossRef] [Google Scholar]
 Sofue, Y., & Nakanishi, H. 2016, PASJ, 68, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Sofue, Y., Honma, M., & Omodaka, T. 2009, PASJ, 61, 227 [NASA ADS] [Google Scholar]
 Spitzer, L. 1968, Diffuse matter in space (New York: Interscience Publication) [Google Scholar]
 Toomre, A. 1963, ApJ, 138, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, X. 1996, ApJ, 457, 125 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Parameters and relations for MiyamotoNagai disks models
In this Appendix, we present the variations of the ratios M/M_{d} and a/R_{dexp} as a function of b/R_{dexp} (M, a, and b are the free parameters of the MiyamotoNagai disks) and the relations between b/R_{dexp} and h_{z}/R_{dexp} for each disk subcomponent, calculated in the same way as in Smith et al. (2015). Here, for simplicity, we denote the exponential scale length R_{dexp} simply as R_{d}. As stated in the main text, these relations are intended to help in the construction of MiyamotoNagai disks models for disks of any mass M_{d}, scale length R_{d} and scale height h_{z}, which can be different from those modelled in the present work. The relations presented for the thin stellar disk, the H I and H_{2} disks, are for models of disks with density depressions in their central regions, while for the thick stellar disk the relations are correspondent to simple radially exponential disks. Regarding the vertical density distributions, the models for the thin and thick stellar disks consider exponential density decays with the increase of the height above the midplane and, for the H I and H_{2} disks, a Gaussian density profile is adopted for the models. We emphasize here that such relations are constructed for models of three MNdisk combinations, which are model 3 for the thin disk and the H_{2} disk (Eqs. (11) and (16)), model 1 for the thick disk (Eqs. (9) and (14)), and model 2 for the H I disk (Eqs. (10) and (15)).
Figure A.1 shows the variation of the thickness ratio b/R_{d} as a function of the ratio h_{z}/R_{d} for the thin and thick stellar disks and the H I and H_{2} disks subcomponents. Filled circles in the figure represent the calculated relations in the way described in Sect. 2.2.3 of the main text. The solid lines represent 4thorder polynomial fits to the distribution of points associated with each calculated relation and are given in the form (A.1)The values of the coefficients k_{0}–k_{4} of the above relation, and for each disk subcomponent, are given in Table A.1.
Figure A.2 shows the variations of the ratios M_{i}/M_{d} (lefthand panel) and a_{i}/R_{d} (righthand panel) as a function of the variation of the ratio b/R_{d} (filled circles), calculated for the modelling of the thin stellar disk. The subscript i = 1, 2, 3 denotes each one of the three MNdisks used in the combination. The points corresponding to the ratio b/R_{d} = 0
denote the values of M_{i} and a_{i}, which are the bestfitting solutions for the three ToomreKuzmin disks combination. The solid lines in each panel of Fig. A.2 represent 4thorder polynomial fits to the points that describe the variations of M/M_{d} and a/R_{d} as a function of b/R_{d}. Each parameter or can then be written as a function of b/R_{d} in the form (A.2)The values of the coefficients c_{0}–c_{4} for the three MNdisk fit models to the thin stellar disk for the range of thicknesses b/R_{d} from 0 to 1.5, are given in Table A.2.
Analogously, the relations between (M_{i}/M_{d}; a_{i}/R_{d}) and b/R_{d} for the thick stellar disk, the H I disk, and the H_{2} disk are shown in Figs. A.3–A.5, respectively. The solid lines in these figures also represent 4thorder polynomial fits to the distribution of points shown as filled circles, and are all also written in the form of Eq. (A.2). The coefficients c_{j} are given in Tables A.3–A.5, for the thick disk, the H I disk, and the H_{2} disk, respectively.
Fig. A.1
Variation of the thickness ratio b/R_{d} as a function of the ratio h_{z}/R_{d}, for the thin and thick stellar disks and the H I and H_{2} disks subcomponents. Filled circles represent the optimal solutions found after the application of the crossentropy algorithm, as described in Sect. 2.2.2 of the main text. Solid lines represent 4thorder polynomial fits to the distribution of points associated with each calculated relation. 
Fig. A.2
Lefthand panel: variation of the three mass parameters of the three MNdisk models as a function of the disk thickness ratio b/R_{d}, for the modelling of the thin stellar disk. Righthand panel: variation of the three scalelength parameters as a function of the ratio b/R_{d}. In the panels, the points represent the optimal solutions found after the application of the crossentropy algorithm. Solid lines are fourthorder polynomial fits to the distribution of points. 
Appendix B: Radial and vertical components of the gradients of the gravitational potential models for each Galactic component
Here we give the explicit forms for the radial and vertical components of the gradients of the gravitational potential expressions modelled for each Galactic mass component, which can be used in the construction of the gravitational forcefield model of the Galaxy.
Appendix B.1: The MiyamotoNagai disks
In the following, we write the components of the gradients for the three models of MiyamotoNagai gravitational potential, which are expressed by Eqs. (14)–(16) in the main text. In all expressions, we use the identity .
Appendix B.1.1: Model 1 (Eq. (14))
Appendix B.1.2: Model 2 (Eq. (15))
Appendix B.1.3: Model 3 (Eq. (16))
Appendix B.2: The bulge (Eq. (19))
Appendix B.3: The dark halo (Eq. (20))
Appendix B.4: The ring potential (Eq. (22))
With the functions ϕ_{R}(R) and ϕ_{z}(z) given in Eq. (22), we have (B.11)(B.12)with the argument in Eq. (B.11).
All Tables
Structural parameters, local surface densities, and masses of the disk components taken as observational prior information to the Milky Way modelling.
Models of 3 TKdisks that provide the best match to the surface density of each disk subcomponent.
Parameters of the three MNdisk combinations for modelling each subcomponent of the Galactic disk.
Local properties resulting from the models MI and MII for the mass distribution and gravitational potential of the Galaxy.
Correlation matrix for the model parameters: the lowerleft triangle for model MI and the upperright triangle for model MII.
All Figures
Fig. 1
Radial profile of the surface density of our “observationbased” model for the Milky Way’s disk. The curves indicate the surface densities of the thin disk (green), thick disk (brown), H I disk (blue), H_{2} disk (violet), and the total disk (red). 

In the text 
Fig. 2
Lefthand panel: radial distribution of the surface mass density of the Galactic disk model; the “observationbased” disk model (blue curve) and the resultant from the sum of all the 3 MNdisk combinations fitted to each disk subcomponent (red curve). Righthand panel: vertical profile of the volume density of the Galactic disk model as a function of the height z from the midplane and at three different arbitrary radii: R = 2 kpc (solid lines); R = 8 kpc (R_{0}) (dashed lines); and R = 15 kpc (dotted lines). The blue curves are also related to the “observationbased” disk model and the red curves to the total 3 MNdisk combinations. 

In the text 
Fig. 3
Contours of isodensities log ρ for the disk models built in the present work. Lefthand panel: the “observationbased” disk model; righthand panel: the sum of all three MNdisk combinations fitted to each disk subcomponent. The values of the contours lie in the range log ρ = [−6; +1], with ρ in M_{⊙} pc^{3}. 

In the text 
Fig. 4
Contours of equipotentials in the Galactic meridional plane for the gravitational potential resulting from the disk model composed of combinations of MNdisks, as described in Sect. 2.2.3. The contours of Φ_{d} are given in 10^{4} km^{2} s^{2}, as indicated in the colourbar. 

In the text 
Fig. 5
Lefthand column, top panel: rotation curve resulting from the Galactic model MI. The circular velocities due to the three mass components are depicted as dotted lines for the bulge, dashed lines for the disk, and dashdotted lines for the dark halo. The total circular velocity is represented by the solid curve. The data for the observed rotation curve are presented as red points with orange error bars for the CO tangentpoint data (Clemens 1985) and the H I tangentpoint data (Fich et al. 1989); and blue points with light blue error bars for the maser sources data (Reid et al. 2014). Bottom panel: velocity residuals between the observed rotation velocities and the total circular model velocities calculated at each radius of the data points. Righthand column: as for the lefthand column but for the Galactic model MII. 

In the text 
Fig. 6
Tangentpoint data from Clemens (1985; green circles) and Fich et al.1989; orange circles). The tangent velocity curves resultant from model MI (blue curve) and model MII (red curve) are also plotted. 

In the text 
Fig. 7
Comparison of the rotation velocity difference 2σ_{Vφ}(σ_{d}) (black circles; see text for details) with the velocity dip V_{dip} (horizontal line) for the maser sources in the radial interval comprised by the dip of the rotation curve. 

In the text 
Fig. 8
Surface density radial profile for the disks of model MI (dashed curve) and model MII (red solid curve). The black solid curve represents a disk with mass equivalent to that from model MII but without the ring density structure. The shaded area in light grey emphasizes the resulting difference in mass between the disks from both models. The blue point denotes the pair (R_{0}; Σ_{0d}). 

In the text 
Fig. A.1
Variation of the thickness ratio b/R_{d} as a function of the ratio h_{z}/R_{d}, for the thin and thick stellar disks and the H I and H_{2} disks subcomponents. Filled circles represent the optimal solutions found after the application of the crossentropy algorithm, as described in Sect. 2.2.2 of the main text. Solid lines represent 4thorder polynomial fits to the distribution of points associated with each calculated relation. 

In the text 
Fig. A.2
Lefthand panel: variation of the three mass parameters of the three MNdisk models as a function of the disk thickness ratio b/R_{d}, for the modelling of the thin stellar disk. Righthand panel: variation of the three scalelength parameters as a function of the ratio b/R_{d}. In the panels, the points represent the optimal solutions found after the application of the crossentropy algorithm. Solid lines are fourthorder polynomial fits to the distribution of points. 

In the text 
Fig. A.3
Same as Fig. A.2, but for the modelling of the thick stellar disk. 

In the text 
Fig. A.4
Same as Fig. A.2, but for the modelling of the H I disk. 

In the text 
Fig. A.5
Same as Fig. A.2, but for the modelling of the H_{2} disk. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.