Models for the 3D axisymmetric gravitational potential
of the Milky Way Galaxy
Key Words.:
Galaxy: fundamental parameters  Galaxy: kinematics and dynamics  Galaxy: structure  Methods: numericalAbstract
Context:
Aims:Galaxy mass models based on simple and analytical functions for the density and potential pairs have been widely proposed in the literature. Disk models constrained by kinematic data alone give information on the global disk structure only 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 with a more realistic description of the density distribution in the disk component.
Methods:From the disk model directly based on the observations (here divided into the thin and thick stellar disks and the H I and H disks subcomponents), we produce fitted mass models 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 which 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.
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, providing a 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), contemporary of 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, among others; Caldwell & Ostriker 1981, among others; 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 forcefield 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, as well as the force perpendicular to the plane of the disk along a given range of Galactic radii. Regarding the latter one, the associated masssurface density, to our better knowledge, is the one 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. Due 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 forcefield above the plane.
Regarding the use of a Galactic potential model for the purpose of orbit calculations, one which has been widely adopted is that of Allen & Santillan (1991). Such model has the attractive characteristics of being mathematically simple and completely analytical, with closed forms for the potential and density, assuring both fast and accurate orbit calculations. Irrgang et al. (2013) have recalibrated the Allen & Santillan (1991) model parameters using new and improved observational constraints.
The main goal of the present work is to provide a fullyanalytical, 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 related to the disk modelling process, can be summarized in the following way:

the structural parameters of the disk  scalelength, scaleheight, radial scale of the central ‘hole’  are based on the Galactic star counts model in the nearinfrared developed by Polido, Jablonski, & Lépine (2013, hereafter PJL) for the case of the stellar disk component; for the gaseous disk counterpart, we adopt recent values returned by surveys of the distribution of hydrogen, atomic H I and molecular H, in the disk;

the density and potential of the disk components are modelled by the commonly used MiyamotoNagai disk profiles (equations 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 (equations 6, 7, 8 and 9 of the abovereferred paper) in order to better fit some of the disk subcomponents. The approach followed for the construction of the MiyamotoNagai disks is based on the one presented by Smith et al. (2015);

a model with a ring density structure added to the disk density profile is 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 also placed a little beyond the solar orbit radius. An explanation for the existence of such ring density structure is given by Barros, Lépine, & Junqueira (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 through the construction of MiyamotoNagai disks versions of the ‘observed’ ones. 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 the estimation of uncertainties are presented in Sect. 5. In Sect. 6, we analyse the results of each mass model by 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’s disk separating it into the stellar (thin and thick disks) and gaseous (H I and H 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, as well as the steps for the construction of the mass and potential disk models. In this paper, we use the cylindrical coordinates (, , ) for the density and potential expressions. The solar Galactic radius is denoted as .
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 have 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 exploration of the parameter space and the estimation of its optimal values were done by the authors with the usage of statistical methods such as the Markov Chain Monte Carlo (MCMC) (Gilks et al. 1996) and the Nested Sampling (NS) algorithm (Skilling 2004). PJL have modelled the radial profile of the density of each subcomponent of the stellar disk by a modified exponential law, based on the Galactic model of Lépine & Leroy (2000). Such profile is equivalent to the Freeman’s Type II disk brightness profile, which contains a depletion in the center, with respect to a pure exponential law (Freeman 1970; Kormendy 1977). The stellar surface densities for the thin and thick disks can then be written as:
(1) 
where corresponds to the local disk stellar surface density (at ); is the radial scalelength; and 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 center 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; differently, the thick disk can be described by a simple radial exponential decay, i.e. .
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 scaleheight . In that case, the authors introduced the variation of the scaleheight with the Galactic radius, , 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 scaleheights with Galactocentric distance for different stellar distributions. In the present study, however, we do not attempt to model such function for , and we consider the scaleheight as a constant along the Galactic radius and with a value equal to the local scaleheight (at ) estimated by PJL, for each thin and thick disks. 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 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 scalelengths , radii of the central hole , and scaleheights , are, as mentioned before, the bestfitting values reported in the PJL model, which are listed in Table 1.
The local stellar surface densities for both thin and thick disks ( in Eq. 1) are based on the model of Flynn et al. (2006) (see also Holmberg & Flynn, 2000, 2004). These authors discriminate the contributions for the total 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 M pc, which can be compared with the recent determination by Bovy et al. (2012) of M pc using the SEGUE spectroscopic survey data. The stellar remnants and brown dwarfs in the Flynn et al. (2006) model contribute with M pc. Taking 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, we end up with M pc, the same value adopted by Read (2014). Separating this last value between the thin and thick disks, we take M pc for the thick disk, as in the Flynn et al. (2006) model, and M pc for the thin disk, where we have assigned the brown dwarfs and stellar remnants to the thin disk for practical purposes. These local surface densities along with the scaleheights result in the local volume densities in the midplane of the Galaxy: M pc for the thick disk; M pc for the thin disk (or M pc considering only mainsequence stars and giants). The thicktothin disk densityratio, (neglecting the stellar remnants/brown dwarfs contribution), is close, given the errors, to the value measured by Jurić et al. (2008) of . The values adopted as constraints for and are listed in Table 1. In Table 1, we also give the total masses calculated for the thin and thick disks, as well as the radial scalelength relative to the region of each disk subcomponent that presents the density exponential decay, and which will be used in the modelling process described in Sect. 2.2.3. Since the thick disk is modelled by a single exponential, .
Component  ^{}^{}For the H I and H disks, the scaleheights are the parameters expressed in Eq. 4.  ^{}^{} corresponds to the radial exponential scalelength fitted to the region of the disk subcomponent where the surface density profile is dominated by the exponential decay.  

(kpc)  (kpc)  (kpc)  (M pc)  ( M)  (kpc)  
thin disk  2.12  2.07  0.205  30.2  2.489  2.18 
thick disk  3.05  0.00  0.640  7.0  0.568  3.05 
H I disk  9.50  1.90  0.180  17.0  1.184  5.00 
H disk  1.48  4.20  0.100  3.0  0.227  1.52 
2.1.2 The gaseous component
We consider the atomic H I and molecular H gaseous disks as 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 for the H I and for the H gas disk components, respectively. Similarly to Eq. 1, corresponds to the local (H I, H) surface density; is the radial scalelength; and is the radius of the ‘central hole’ in the density distribution of H I and H disks. The form for the gas surface densities presented in Eq. 3, together with the values for the structural parameters and of the H I and H 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 Figure 3), and in H by Nakanishi & Sofue (2006, cf. their Figure 9).
The volume density for both H I and H disks is given in the form:
(4) 
with the scaleheights corresponding to the halfwidth at halfmaximum of the density peaks of H I and H 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 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, also being an expected solution for hydrostatic equilibrium considerations.
The scaleheight 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, among others; Burton 1976, among others; 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 of the H I gas is approximately constant with radius (e.g. Spitzer 1968), the scaleheight of this component must increase with . The Galactic distribution of molecular gas also shows a flaring consistent with that observed in H I (Kalberla et al. 2007, and references therein). Although being a very important feature to be included in any realistic mass model of the Galaxy, we do not attempt to model the flaring of the gaseous disk component for the same reason exposed in the case of the stellar disk. Therefore, here we consider the scaleheights of the H I and H disk components as independent of radius and with values equal to the local ones (at ) calculated from the model of Amôres & Lépine (2005, cf. their Equation 7). Table 1 lists the adopted prior values for the H I and H disks scaleheights.
The model of Holmberg & Flynn (2000) for the ISM component, which is based on the original multiphase model of Bahcall, Flynn, & Gould (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, according to this model, is M pc, with an uncertainty of . 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 both in the form of optically thick cold neutral hydrogen and COdark molecular gas is known to contribute to the ISM density, which must raise the estimates of the local midplane gas densities by as much as . 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 are even higher. Therefore, as observational constraints to our gaseous disk model, we take the following values for the local surface gas densities: M pc, as in the Hessman (2015) updated estimate, which already takes into account the correction factor of 1.36 for the mass in helium (He); M pc, as in the Holmberg & Flynn (2000) model. The warm gas disk component, which contributes locally with 2 M pc (Holmberg & Flynn 2000; Hessman 2015), is here incorporated to the H I disk for practical purposes, increasing the value of to 17 M pc. The adopted values for the local surface densities, total masses and radial exponential scalelengths of the H I and H disks are listed in Table 1.
As can be seen from Table 1, our ‘observationbased’ disk model comprises a total mass M, being M in stars and M in the gaseous form. The total local disk masssurface density is M pc, which can be compared with the updated density estimates of 58 M pc by Hessman (2015), and M pc by Read (2014). The local disk midplane volume density is M pc. This value is somewhat larger than the estimate of the local dynamical mass density of M pc by Holmberg & Flynn (2000), of which 0.095 M pc is in visible matter. However, the density corrections discussed by Hessman (2015) should increase these traditional estimates to local densities as large as or even M pc. Figure 1 shows the radial distribution of the surface density resulted from our ‘observationbased’ disk model, as well as the surface densities of each disk subcomponent. The surface densities of the thin and thick stellar disks drop to M pc at the radius kpc, while the H disk reaches this value at kpc. The H I disk subcomponent extends to larger radii, with a drop of its surface density to M pc at kpc, 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 , must involve the use of some mathematical tools and/or numerical techniques, like the multipole expansion and numerical interpolation employed by Dehnen & Binney (1998) in their determination of the potential and forcefield 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 so for the gravitational potential, for given ranges of Galactic radii and heights above the disk midplane. This is the approach that we adopt 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 forcefield, 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.
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 . The ‘model 1’ of this family, of what here we refer to as ToomreKuzmin disks (TKdisks), is described by a surface density which is written as:
(5) 
where is the total disk mass and is related to the radial scalelength of the disk. The choice of use of TKdisks, in the first step, is due to the fact that all the information about the disk surface density can be recovered adjusting the only two parameters and . As shown later, the threedimensional structure of the disks is obtained after ‘inflating’ the TKdisks with the introduction of the parameter related to the scaleheight, in the same way as in the original method of Miyamoto & Nagai (1975).
From Eq. 5, it can be seen that the density decreases with 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 is the motivation to the use of a combination of three MNdisks (related to the ‘model 1’ of TKdisks, Eq. 5), firstly 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 scalelength and mass , being one of the masses with a negative value. The MNdisk with negative mass is also the one with the largest scalelength, a feature that helps to decrease the surface density at large radii (and improve the fit to an exponential disk), but that 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 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, due to linearity between and in the Poisson equation, new potentialdensity pairs can be obtained through the derivation of times with respect to (Binney & Tremaine 2008). For example, the densities for ‘models 2’ and ‘3’ of the TKdisks are described by
(6) 
and
(7) 
respectively (Miyamoto & Nagai 1975). One can see that the densities and decrease with and , 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 the H disk subcomponents. The difference here is that the scalelength of this component with negative mass does not need to be the largest one, thus avoiding the negative densities at large radii. Therefore, for each disk subcomponent (thin, thick, H I, H), we search for the combination of 3 TKdisks, of each model (1, 2, 3), that better reproduces the radial profile of the surface density of each subcomponent. Let us take, for example, the case of the thin disk. We search for the set of parameters ^{2}^{2}2 , where denotes each one of the three TKdisk, and denotes the TKdisk model order. of ‘model 1’ of 3 TKdisks that generates the best fit to the radial distribution of surface density of the thin disk. These chosen 3 TKdisks result in the total surface density . The same procedure is done with the 3 TKdisks in the forms of ‘model 2’ and ‘model 3’, which result in the total densities and , respectively. We quantify 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 refers to the disk subcomponent ( = thin in the above example) and denotes the model order of TKdisks; is the number of radial bins over which the sum is calculated. The TKdisk model order that results in the lowest value for is the one which we consider as the model that better represents the surface density of the disk subcomponent. The search for the best sets of parameters corresponding to the ’s best TKdisk models were done 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. We refer the reader to recent papers that 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 in 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 3 TKdisks found to better reproduce the surface density of each one of the disk subcomponents, as well as the corresponding (Eq. 8) minimized within the CE algorithm.
Component  Model ^{}^{}The model order of 3 TKdisks that results in the lowest value for (Eq. 8), which also corresponds to the associated model of 3 MNdisks.  

thin disk  3 (Eq. 7)  4.35 
thick disk  1 (Eq. 5)  1.12 
H I disk  2 (Eq. 6)  0.25 
H disk  3 (Eq. 7)  0.74 
2.2.2 The crossentropy algorithm
For a detailed presentation and description of the CE method, we refer the reader to the papers by Monteiro, Dias, & Caetano (2010) and Dias et al. (2014). Here, we briefly show an overview of the method and how it works.
Let us suppose that we have a set of data with individual points , , …, and we wish to study it in terms of an analytical model with a vector of parameters , , …, . 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 independent sets of model parameters , where , under some chosen distribution, and the subsequent minimization of an objective function used to transmit the quality of the fit during the run process. As the method converges to the ‘theoretical’ exact solution, then , which means .
The CE method performs an iterative statistical coverage of the parameter space, where the following is done in each iteration (Monteiro et al. 2010):
 (i)

Random generation of the initial parameter sample, respecting some 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 lower values for the objective function ;
 (iii)

Random generation of updated parameter samples from the previous best candidates  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 done in this work, we followed the general step by step presented in section 2.2 of Monteiro et al. (2010). The tuning parameters used in the run processes were: sets of trial model parameters per iteration; a maximum number of 100 iterations; , 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: ; ; and . 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 done by replacing the term in the potential function, associated to the density of a given model of TKdisk, with the term , where expresses the vertical scaleheight of the disk. In this way, corresponding to the three models of TKdisks expressed by Eqs. 5, 6, and 7, we have the following threedimensional density functions that compose the first three models of MNdisks:
(9) 
(10) 
(11) 
where . In the above expressions, for the sake of good readability, we avoided the excessive use of subscripts in the parameters , , 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 3 TKdisks of a given model, now we use the corresponding combination of 3 MNdisks of the equivalent model.
The task now is to find the best values for the parameter 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, scalelength, and for an ample range of thicknesses, from infinitely thin disks to approximately spherical systems. The basic idea of the method consists in: starting from the model composed by the 3 TKdisks (for which ) that better match the radial profile of the disk surface density, one chooses the corresponding model composed by 3 MNdisks with a small value for the parameter (with the aim of remaining near the solution for the infinitesimal thickness disk), and whose parameters and of the MNdisks deviate from small amounts with respect to those of the TKdisks. With the continuous variation of in small steps, the parameters and are searched for at intervals close to the solution found in the previous step, ensuring a smooth variation of both and with the variation of the disk thickness . In fact, Smith et al. (2015) analyze the variation of the ratios and as a function of the variation of the ratio , being the mass and the radial scalelength of the exponential disk to be modelled. Therefore, for each small variation of the ratio , one searches for the parameters and of the 3 MNdisks whose integral of the density in the vertical direction better reproduces the radial profile of the disk surface density .
To relate the scaleheight of the MNdisks with the scaleheight of the vertically exponential disk, Smith et al. (2015) follow an approach analogous to that described above. In this case, the authors first sum up the fractional differences between the density of the 3 MNdisks combination and the exponential one, measured vertically from the midplane up to . They vary the ratio in order to minimize the sum, then finding the best match to the exponential distribution.
In this section, we only present the values of the parameters calculated for each combination of 3 MNdisks that fit each disk subcomponent with their fixed values for listed in Table 1. However, in Appendix A, we also present the variations of the ratios and as functions of , as well as the relations between and for each disk subcomponent, calculated in the same way as in Smith et al. (2015). Therefore, disks of any masses and scales and 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 the 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 interested user to construct models for disks with masses and sizes different from the ones modelled in this work.
In our search for the best values of the parameters of the MNdisks, we have also employed the CE technique described in the previous subsection, but now looking 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 3 MNdisks combination
(12) 
for the relations of and with , and the minimization of the sum of the squares of the fractional differences between the volume densities of the 3 MNdisks combination and the disk subcomponent over a given vertical range and at
(13) 
for the relations between and .
Let us take again the case of the thin stellar disk for exemplification. With the ratio calculated from the values listed in Table 1 for the thin disk, and putting it into Eq. 28 with the coefficients for the thin disk given in Table 7, both from Appendix A, one finds the corresponding value for the ratio . Substituting this last one into Eq. 29 with the coefficients given in Table 8, also from Appendix A, one obtains the values of the parameters and for the combination of 3 MNdisks of model 3 (Eq. 2.2.3) that better matches the density distribution of the thin stellar disk. We emphasize here that a single scaleheight is used for all the three MNdisks of the combination, as originally done by Flynn et al. (1996) and followed by Smith et al. (2015). All the procedure described above, which was illustrated with the thin disk, is repeated with the other subcomponents (thick disk, H I and H disks) using their corresponding equations and tables in Appendix A. Table 3 lists the values of the parameters of each combination of 3 MNdisks that fit each disk subcomponent, as well as the correspondent models of MNdisks that had been previously found with the TKdisks modelling.
Component  Model ^{}^{}The model of 3 MNdisks combination correspondent to the TKdisks model that results in the lowest value for the quantity in Eq. 8.  

( M)  (kpc)  ( M)  (kpc)  ( M)  (kpc)  (kpc)  
thin disk  2.106  3.859  2.162  9.052  1.704  3.107  0.243  3 (Eqs. 2.2.3 and 16) 
thick disk  0.056  0.993  3.766  6.555  3.250  7.651  0.776  1 (Eqs. 9 and 14) 
H I disk  2.046  9.021  2.169  9.143  3.049  7.758  0.168  2 (Eqs. 10 and 15) 
H disk  0.928  6.062  0.163  3.141  0.837  4.485  0.128  3 (Eqs. 2.2.3 and 16) 
Figure 2 shows, on the lefthand panel, the radial profile of the surface massdensity for the ‘observationbased’ disk model (blue curve), as well as the total surface density resultant from all the 3 MNdisks combinations fitted to each disk subcomponent (red curve). It can be seen a good agreement between these two curves, which denotes the great effectiveness of the method. The fractional differences between these surface density distributions are at , reaching 50% at kpc and 100% at kpc, where the absolute values of the densities become lower than M pc. These features show the tendency of the MNdisks in returning larger densities at large radii. The local disk surface density returned by the MNdisks fit model is M pc, which can be considered close to the prior value of 57.2 M pc of the ‘observationbased’ disk model, if we take the uncertainty of M pc over this quantity as estimated by Read (2014). The local surface density integrated within 1.1 kpc of the disk midplane returned by the MNdisks fit model is M pc, 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 from the disk midplane, taken at three arbitrary radii: kpc (solid line), kpc (dashed line), and kpc (dotted line). The blue and red curves also correspond to the density distributions of the ‘observationbased’ disk model and the MNdisks fit model, respectively. The fractional differences between the volume density distributions are lesser than for kpc at kpc, for kpc at kpc, and for kpc at kpc, reaching values greater than beyond these heights. These features denote the systematic trend of the MNdisks in returning higher densities at greater distances from the disk plane, while at relatively small the density distributions show roughly good matches at radii up to at least kpc. As pointed out by Smith et al. (2015), it is impossible for the 3 MNdisks model to perfectly reproduce the vertically exponential density profile, or even the type profile, since they are mathematically distinct. The same can be stated on the vertical Gaussian profile adopted for the H I and H disks in the present work; at heights lower than 1 kpc the 3 MNdisks start presenting great deviations from the Gaussian profile. This can be noticed analysing the distributions of calculated at 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 ( in M pc) in the meridional plane of the Galaxy: for the ‘observationbased’ disk model (lefthand panel), and for the total 3 MNdisks combinations fit model (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, 10 and 2.2.3, are written as ‘model 1’, ‘model 2’ and ‘model 3’, respectively, in the form:
(14) 
(15) 
(16) 
where again . Therefore, the gravitational potential of the thin stellar disk is modelled with 3 components of the potential (Eq. 16), and whose parameters , and () are given in the first row of Table 3. Equivalently, the thick stellar disk potential is modelled with 3 components of the potential (Eq. 14), the H I disk potential is modelled with 3 components of (Eq. 15), and the H disk potential with 3 components of the potential (Eq. 16); the parameters , and for these disks are given in the second, third and fourth rows, respectively, of Table 3. The total gravitational potential of the disk is then given by the sum of all the combinations of 3 MNdisks that fit each Galactic disk subcomponent:
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 by 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 de Vaucouleurs (1977) surface brightness law over a large range of galactic radius:
(18) 
where is the total mass and 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 spheroid to disk local stellar density ratio for normalization of the density profile, and the spheroid scale radius kpc. Considering the local stellar density of our ‘observationbased’ disk model, the mass of the spheroid of PJL’s model inside kpc (radius which encompasses % of the total spheroid mass) results in M. We use this information to constrain the initial ranges of the parameters and of our bulge models for the search for their best values in the fitting procedure described in Sect. 5. For the bulge mass we adopt the initial guess interval M; for the scale radius we use kpc. These intervals comprise masses inside kpc for our bulge models with possible values in the range M.
The fact that the observed rotation curve, which is adopted to constrain the Galactic models, is nearly flat in the interval (see Sect. 4), besides the fact that in our models (see Sect. 6) the sole disk+bulge density distributions are not capable to give such a support to the rotation curve at these radii, lead 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, here we do not discuss about the reality of such dark halo component, nor about its material content, although it contributes as an extra term of mass in our Galactic models. For this reason, we do not go deeper in the modelling of such 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 is the core radius; is the circular velocity at large (i.e., relative to the core radius); and is the axis ratio of the equipotentials. For simplicity, we consider a spherical dark halo, for which . The density distribution corresponding to the potential is given by (for ):
(21) 
This profile yields a flat rotation curve at large radii. However, as we make clear in the next sections, we are interested in a description of the Galactic potential which 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 our concern in the present study.
The values of the parameters and that describe the bulge component, and and 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 details in the subsections of Sect. 2. In general, the overall disk is modelled by a superposition of distinct models of MiyamotoNagai disks (combinations of 3 MNdisks for each one of the 4 subcomponents: thin, thick, H I and H disks 12 MNdisks in total). Each MNdisk has 3 free parameters (, , ), but since each combination of 3 MNdisks is modelled with a single value for (see Sect. 2.2.3), we have a total of 28 parameters for the construction of the disk (all of them are given in Table 3).
In the present work, differently from other studies, the values of the structural parameters of the disk (here represented by the length scales and 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. It must be said that, 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.1 we give an estimate of the uncertainties on the parameters and 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 by the fact that, due to the way the MNdisks were constructed, the scalelengths were found as functions of the scaleheights , 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, in order to probe the correlations among the parameters. Moreover, it is well known in the literature that the kinematic data commonly used for the models fitting 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 disks scaleheights would be found able to reproduce the kinematic data equally well. Another problematic issue in varying the scales and of all the MNdisks in the fitting of the models, taking into account the uncertainties on the values of , and of all the disk subcomponents, is the high 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 which gives information on the dynamical mass of the Galaxy, but also by prior informations 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 scalelengths and scaleheights of the MNdisks as fixed and equal to the ones given in Table 3. The total mass of the disk will be adjusted introducing the parameter , a scale factor by which the disk total mass is then multiplied. Once the best value for is found, then the masses 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 of the MNdisks have already been found after the construction of the disk mass model in Sect. 2, in the next step we will use their single values in Table 3 as first guesses for the optimal fit values to the dynamical constraints, but keeping unaltered their relative contributions to the disk total mass since all of them will be rescaled by the same factor .
3.2.1 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 one 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 Figure 14). However, being relatively young, the Cepheids possibly trace the gas distribution, and therefore the recent past of star formation. There is a theoretical argumentation and simulations given in BLJ, predicting that the corotation resonance scatters out the stars 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, centered at the corotation radius. BLJ performed numerical integrations of testparticles 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 of the background density is formed at the corotation radius, in a time interval billion years of the evolution of the system. Lépine, Mishurov, & Dedikov (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 Figure 4) with the observed H I radial distribution presented by Burton (1976, see his Figure 6) and found close similarities between them. Such results are compatible with evidences for a ring void of gas as observed by Amôres, Lépine, & Mishurov (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 one (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 kpc. A common association between these two features was also modelled by Sofue, Honma, & Omodaka (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 density feature as a wavy ring superposed on the surface density profile of the disk, in a similar way as done by Sofue et al. (2009).
Since we still have no clue of how could be the threedimensional distribution of this density feature, here we make the simplest assumption that it would be mainly detected in the density distribution very near the midplane of the disk, so we can take the zerothickness disk approximation for the surface density of the ring. Here 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 is the amplitude of the potential related to the amplitude of the minimum and maximum densities of the ring, given in units of km s; is a parameter related to the ring width; is the ring node radius  where sits the inflection point between the minimum and maximum density features; and is the scaleheight of the ring potential. Solving the Poisson equation for the ring potential (within the zerothickness assumption), we can find the corresponding expression for the ring surface density .
In the model that includes the ring structure, the surface density is added to the disk surface density , but since the ring is modelled by a minimum followed by a maximum density of similar amplitudes, which makes the net ring mass , 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: , and ; the scaleheight 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 and rotation velocities at , the local angular rotation velocity , and values for the estimated total local disk masssurface density 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 from the statistical analysis performed by Malkin (2013) on 53 measurements published in the literature over the last 20 years, which average value, for practical purposes, is recommended by the author as being
For the peculiar velocity of the Sun relative to the local standard of rest (LSR), we take the reevaluation proposed by Schönrich, Binney, & Dehnen (2010, hereafter SBD):
where we use a lefthanded system for (, , ), in which is positive towards the Galactic anticenter, is positive in the direction of Galactic rotation, and is positive towards the direction of the North Galactic Pole; 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 mas yr. Thus we have km s kpc. This last equality comes from the assumption that Sgr A is at rest at the Galactic center, 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. (e.g. Honma et al. 2012). We then have
(23) 
for the local angular rotation velocity . With the abovequoted values for , and , one obtains km s kpc. The uncertainty on is calculated as , which returns the value km s kpc; to the uncertainty of 2 km s given by SBD, we have added 1 km s to allow for a possible peculiar motion of Sgr A at the Galactic center (McMillan & Binney 2010). The corresponding rotation velocity at the LSR, , assumes the value km s.
4.2 The rotation curve
Tangent velocities:
The tangent (or terminal) velocity is usually assumed as the maximum velocity of the ISM gas along a given lineofsight at Galactic coordinates and , which can be related to the circular velocity at the tangent point (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)^{5}^{5}5The tangent velocity data on table 2 of Clemens (1985) were corrected for 3 km s line width, and a correction of 7 km s 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., which cover longitudes in the first Galactic quadrant. We also use the H I tangent point data from table 2 of Fich, Blitz, & Stark (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 and the rotation velocities were calculated using the Galactic constants adopted in this work (, ) = (8 kpc, 230 km s). We have propagated the uncertainties on both and to the uncertainties on and ( and ), respectively.
It has often been argued that the central region of the Galaxy is strongly affected by nonaxisymmetric structures like 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, Renaud, & Soubiran (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 kpc. Perhaps the most delicate feature in the Milky Way inner rotation curve, as derived from the tangentpoint method, is the velocity peak at 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 which is not present in the true rotation curve, when the bar major axis is viewed with angles with respect to the direction of the galactic center. However, we argue that the interpretation of this inner peak as being due to the contribution of the bulge can still be defensible: at 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 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 one only for radii kpc. The unreliable rotation curve at kpc in the Chemin et al.’s simulation is attributed to effects associated with the corotation of the bar, which have, in their model, a pattern speed of the order of 59 km s kpc. Other authors have obtained lower values for the bar pattern speed, e.g. km s kpc (RodriguezFernandez & Combes 2008), or 33 km s kpc (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 center which is of the order of 12 km s on average, in the range kpc. We think that this amplitude of asymmetry, contrary to the one observed in the region kpc, does not have influence on our models. Given the considerations presented above, we decided to restrict the tangent velocity data to (e.g. Dehnen & Binney 1998), which in our case is equivalent to Galactic radii kpc. We end up with a data set which totalizes 280 rotation velocity measurements in the inner solar circle region ( kpc).
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 can have access to 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 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^{6}^{6}6http://bessel.vlbiastrometry.org and the Japanese VLBI Exploration of Radio Astrometry (VERA)^{7}^{7}7http://veraserver.mtk.nao.ac.jp). 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 , , velocities and uncertainties , and for each source, following the formalism described by Johnson & Soderblom (1987). Correcting for the SBD’s solar peculiar motion and for the LSR circular velocity , we calculated the Galactocentric component of the space velocity of each maser source in the direction of Galactic rotation. The uncertainties on , , 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 are then defined by and . 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 have discarded the sources with radii kpc because most of them present values of with large deviations from the rotation curve 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 traced by the maser sources is restricted to few data points, these measurements are the best ones 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 M pc, 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 M pc for this quantity. We adopt the same uncertainty on of M pc. Holmberg & Flynn (2004) also estimated the local surface density integrated within 1.1 kpc of the disk midplane as being M pc, which, according to the authors, takes into account both disk and dark halo contributions. We take such value as observational constraint on and its uncertainty.
5 Fitting procedure
As exposed in Sect. 3, we attempt to construct two types of Galactic models: those which do not incorporate the ring density structure in the disk (Sect. 3.2.1), and which we refer to as model MI; and those which do incorporate the ring density structure and will be denoted by model MII. The free parameters of model MI are the bulge parameters and , the dark halo parameters and , and the disk mass scale factor . The free parameters of model MII are the same of model MI, plus the ring parameters , and ; the scaleheight 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 minimization procedure, which is implemented through the crossentropy (CE) algorithm (Sect. 2.2.2). We first try initial guesses for the parameters sets 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 independent sets of model parameters by considering initial uniform distributions centered 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 minimization criterion (the sets that return the lowest values for ). From the mean and standard deviation of each ensemble of 100 parameters of this elite sample, we generate other independent sets of trial parameters through Gaussian distributions 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, , in the case of model MI, and the addition of the ring potential in the case of model MII, . Considering the balance between the centrifugal force and gravity, the circular velocity of a given model, measured for instance in the Galactic midplane , is linked to the total gravitational potential by the form
(24) 
Therefore, the gravitational potential and the model rotation curve are totally defined by the set of parameters of each model MI and MII described above. The radial derivative of the potential 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 of the rotation curve has a pair of values (, ). In the fitting process of the rotation curve, we search for the minimization of the residuals between the rotation velocities of the observational data and the circular velocities resulted from the model at each rad