# A Parametric Galactic Model toward the Galactic Bulge Based on Gaia and
Microlensing Data

Naoki Koshimoto  Laboratory for Exoplanets and Stellar Astrophysics,
NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, MD 20742, USA  Junichi Baba
National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan
David P. Bennett  Laboratory for Exoplanets and Stellar Astrophysics,
NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, MD 20742, USA

###### Abstract

We developed a parametric Galactic model toward the Galactic bulge by fitting
to spatial distributions of the Gaia DR2 disk velocity, VVV proper motion,
BRAVA radial velocity, OGLE-III red clump star count, and OGLE-IV star count
and microlens rate, optimized for use in microlensing studies. We include the
asymmetric drift of Galactic disk stars and the dependence of velocity
dispersion on Galactic location in the kinematic model, which has been ignored
in most previous models used for microlensing studies. We show that our model
predicts a microlensing parameter distribution significantly different from
those typically used in previous studies. We estimate various fundamental
model parameters for our Galaxy through our modeling, including the initial
mass function (IMF) in the inner Galaxy. Combined constraints from star counts
and the microlensing event timescale distribution from the OGLE-IV survey, in
addition to a prior on the bulge stellar mass, enable us to successfully
measure IMF slopes using a broken power-law form over a broad mass range,
$\alpha_{\rm bd}=0.22^{+0.20}_{-0.55}$ for $M<0.08\,M_{\odot}$, $\alpha_{\rm
ms}=1.16^{+0.08}_{-0.15}$ for $0.08\,M_{\odot}\leq\,M<M_{\rm br}$, and
$\alpha_{\rm hm}=2.32^{+0.14}_{-0.10}$ for $M\geq\,M_{\rm br}$, as well as a
break mass at $M_{\rm br}=0.90^{+0.05}_{-0.14}\,M_{\odot}$. This is
significantly different from the Kroupa IMF for local stars, but similar to
the Zoccali IMF measured from a bulge luminosity function. We also estimate
the dark matter mass fraction in the bulge region of $28\pm 7$% which could be
larger than a previous estimate. Because our model is purely parametric, it
can be universally applied using the parameters provided in this paper.111 A
tool for microlensing simulation using our Galactic model has been published
(Koshimoto & Ranc, 2021), and can be downloaded at
https://github.com/nkoshimoto/genulens (catalog
https://github.com/nkoshimoto/genulens).

Milky Way Galaxy (1054), Galactic bulge (2041), Initial mass function (796),
Gravitational microlensing (672)

##  1 Introduction

Gravitational microlensing is a unique tool that can study the population of
planetary to black hole mass objects in our Galaxy (Gaudi et al., 2019). It
has been applied to study the population of exoplanets since the proposal by
Mao & Paczynski (1991), along with more than 100 exoplanets, has been
discovered as of 2020222https://exoplanetarchive.ipac.caltech.edu/ (catalog
NASA Exoplanet Archive) (Akeson et al., 2013). Suzuki et al. (2016) used a
statistical sample of 30 microlensing planets to reveal a likely peak at
$q\sim 10^{-4}$ in the mass-ratio function of planets beyond the H2O snow line
for the first time, challenging the core accretion theory (Suzuki et al.,
2018). Microlensing also enables us to study the population of free-floating
planets (Sumi et al., 2011; Mróz et al., 2017). Further, some candidates of
isolated black holes have also been discovered (Bennett et al., 2002;
Poindexter et al., 2005; Wyrzykowski et al., 2016; Wyrzykowski & Mandel,
2020).

Difficulty in measuring the lens mass and distance makes microlensing studies
complex. Although four physical quantities, namely the lens mass $M_{\rm L}$,
distance from lens $D_{\rm L}$, distance from source $D_{\rm S}$, and lens-
source relative proper motion $\mu_{\rm rel}$, are involved in each
microlensing event, the Einstein radius crossing time

$\displaystyle t_{\rm E}=\frac{\theta_{\rm E}}{\mu_{\rm rel}},$ (1)

is the only quantity that is measurable for all events, where $\theta_{\rm E}$
is the angular Einstein radius given by

$\theta_{\rm E}=\sqrt{\kappa M_{\rm L}\pi_{\rm rel}}$ (2)

with $\kappa=8.144~{}{\rm mas}~{}M_{\odot}^{-1}$ and $\pi_{\rm rel}={\rm
1\,AU}(D_{\rm L}^{-1}-D_{\rm S}^{-1})$. For microlensing events toward the
Galactic bulge, the source star is assumed to be a bulge star, i.e., $D_{\rm
S}\sim 8~{}$kpc, but the three remaining quantities still degenerate in a
$t_{\rm E}$ value.

There are three observable quantities each of which provides a mass-distance
relation: the angular Einstein radius, $\theta_{\rm E}$, microlens parallax,
$\pi_{\rm E}$, and lens star flux. The degeneracy in $t_{\rm E}$ could be
disentangled when any two of the three are observed. The angular Einstein
radius is given by Eq. (2) while the microlens parallax is given by

$\pi_{\rm E}=\frac{\pi_{\rm rel}}{\theta_{\rm E}},$ (3)

and these two quantities can be measured through higher order effects in the
light curve, finite source effect (Yoo et al., 2004) and annual parallax
effect caused by the Earth’s orbital motion (Alcock et al., 1995; An et al.,
2002), respectively. However, these effects are rarely measured, particularly
for single-lens events that account for $\sim$90 % of all microlensing events.
We can measure the lens flux with high-angular resolution follow-up imaging
using adaptive optics (AO) or the Hubble Space Telescope (HST) (Batista et
al., 2015; Bennett et al., 2015; Bhattacharya et al., 2018). However, such
measurements typically require us to wait for $\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}5$ years until the lens-source
separation can be measured, depending on the relative proper motion and
contrast between lens and source stars. These limitations associated with the
lens mass and distance measurements have resulted in many microlensing studies
suffering from a chronic lack of information.

From this perspective, a Galactic model, which refers to a combination of a
stellar mass function and stellar density and velocity distributions in our
Galaxy, has played a crucial role in providing prior probability distributions
for microlensing events toward the Galactic bulge. For example, a large number
of studies on individual event analysis use a Galactic model to calculate a
posterior distribution of the lens mass and distance (e.g., Alcock et al.,
1995; Beaulieu et al., 2006; Koshimoto et al., 2014; Bennett et al., 2014).
Further, it has also been used for statistical studies to estimate fundamental
parameters in our Galaxy, such as the slope of initial mass function (IMF)
(Sumi et al., 2011; Mróz et al., 2017; Wegg et al., 2017) or dark matter
fraction in the bulge region (Wegg et al., 2016). Galactic models are
occasionally used to diagnose whether measurements of microlensing parameters
are contaminated by systematic errors (Penny et al., 2016; Koshimoto &
Bennett, 2020).

However, there are concerns regarding the use of Galactic models because the
results depend on the choice of models, with most Galactic models used in
microlensing studies employing several simplified features (Yang et al.,
2020). This is a concern specifically for statistical studies in which small
effects for many individual events can be combined. For example, Koshimoto &
Bennett (2020) compared three different Galactic models by Sumi et al. (2011),
Bennett et al. (2014), and Zhu et al. (2017) with 50 microlens parallax
measurements by the 2015 Spitzer microlensing campaign (Zhu et al., 2017) and
concluded that the parallax measurements are significantly contaminated due to
systematic errors based on diagnoses using the Anderson-Darling (AD) tests.
Here, all of the three models use simplified models for disk kinematics with
flat rotation curves, constant velocity dispersions, and no asymmetric drift,
which contradict the observed velocity distributions (Gaia Collaboration et
al., 2018). Although the conclusions seem robust because they are based on the
discovered correlations between $p$-values from the AD tests and some
characteristics of events, which are likely related to vulnerability to
systematic errors (e.g., source brightness or peak coverage), the absolute
$p$-values might be misestimated due to the simplified features in the
Galactic models.333In Section 6.1, we show that our improved Galactic model
does not have a significant effect on the results of Koshimoto & Bennett
(2020) that were calculated with simplified Galactic models.

Another case motivating this concern is the IMF in brown dwarf mass range of
$M<0.08\,M_{\odot}$, parameterized by a slope, $\alpha_{\rm bd}$, with
$dN/dM\propto M^{-\alpha_{\rm bd}}$, and inferred by the $t_{\rm E}$
distribution from the OGLE-IV survey (Mróz et al., 2017). Although Mróz et al.
(2017) shows the likelihood peak at $\alpha_{\rm bd}\sim 0.8$ using the
Galactic model based on Han & Gould (1995, 2003) with simplified disk
kinematics, Specht et al. (2020) shows the likelihood peak at $\alpha_{\rm
bd}=-0.1$, and 0.8 (-0.8 in their definition) is significantly disfavored
according to their Figure 4. Specht et al. (2020) used the Besançon model
(Robin et al., 2003, 2012, 2017) for their analysis using a more realistic
model for disk kinematics. However, the Besançon model also has some
unrealistic features, as discussed by Penny et al. (2019), such as the small
Galactic bar angle of 13∘ compared to $25^{\circ}$–$30^{\circ}$ implied by
other observations (see a review by Bland-Hawthorn & Gerhard, 2016). Further,
the IMF for bulge stars in the public version has a minimum mass at
$0.15\,M_{\odot}$, a too massive cut-off for applying microlensing analysis,
which is sensitive down to planetary-mass objects.

Several non-parametric dynamical models of our Galaxy are developed with the
aid of $N$-body simulations (Wegg & Gerhard, 2013; Portail et al., 2015,
2017). These models are better than parametric models with respect to
consistency with the dynamics, but they have limited resolution. Numerical
models use particles of thousands of solar masses and it is not easy to employ
these models for simulations of microlensing events along a specific line of
sight. Also, they are not open to communities outside the group, and thus not
useful in terms of accessibility or difficulty of reproduction. A recent
version of the Besançon model uses results from an $N$-body simulation
(Gardner et al., 2014) for the kinematics of bulge stars, and it is difficult
for other people to reproduce the model for their own use.

In this study, we develop a parametric Galactic model by fitting to spatial
distributions of the median velocity and velocity dispersion from Gaia DR2
(Gaia Collaboration et al., 2018), OGLE-III red clump (RC) star count (Nataf
et al., 2013), VIRAC proper motion data (Smith et al., 2018; Clarke et al.,
2019), BRAVA radial velocity data (Rich et al., 2007; Kunder et al., 2012),
and the OGLE-IV star count and microlensing event data (Mróz et al., 2017,
2019). We parameterize the model with 40–48 parameters depending on the
selected options for the bulge density profile, and we present the best-fit
values for all the parameters to enable ease of reproduction or make further
updates. Our model aims to fulfill the gap in the market between the
simplified parametric models used for microlensing studies and dynamical
models based on $N$-body simulations that are difficult to use for
microlensing work.

This model enables, for the first time, simultaneous measurements of the three
slopes over entire mass range in a broken power law IMF of $dN/dM\propto
M^{-\alpha}$, $\alpha_{\rm bd}=0.22^{+0.20}_{-0.55}$ in $M<0.08\,M_{\odot}$,
$\alpha_{\rm ms}=1.16^{+0.08}_{-0.15}$ in $0.08\,M_{\odot}\leq M<M_{\rm br}$,
and $\alpha_{\rm hm}=2.32^{+0.14}_{-0.10}$ in $M\geq M_{\rm br}$, as well as a
break mass at $M_{\rm br}=0.90^{+0.05}_{-0.14}\,M_{\odot}$. The IMF
corresponds to the stellar mass-to-light ratio in $K$-band of
$\Upsilon_{K}=0.72^{+0.05}_{-0.02}\,M_{\odot}/L_{K_{\odot}}$, which is
significantly different from $\Upsilon_{K}=1.04\,M_{\odot}/L_{K_{\odot}}$ with
the Kroupa (2001) IMF for local stars but consistent with
$\Upsilon_{K}=0.75\,M_{\odot}/L_{K_{\odot}}$ with the Zoccali et al. (2000)
IMF for bulge stars. We estimated the dark matter mass inside the VVV bulge
box ($\pm 2.2\times\pm 1.4\times\pm 1.2$ kpc, Wegg & Gerhard, 2013) of
$M_{\cal DM}=0.45^{+0.02}_{-0.06}\times 10^{10}\,M_{\odot}$ by comparing our
$\Upsilon_{K}$ value with those of the five versions of dynamical Galactic
models in Portail et al. (2015). This is consistent with another independent
estimate of $M_{\cal DM}=0.51^{+0.12}_{-0.11}\times 10^{10}\,M_{\odot}$
derived by combining our estimate on the stellar mass inside the VVV bulge box
of $1.34^{+0.10}_{-0.11}\times 10^{10}\,M_{\odot}$ and the dynamical mass in
the same region of $(1.85\pm 0.05)\times 10^{10}\,M_{\odot}$ constrained by
Portail et al. (2017).

Using the developed models, we repeated a part of the analysis by Koshimoto &
Bennett (2020) and confirmed their finding of a significant discrepancy
between the predicted and observed $\pi_{\rm E}$ distributions with the new
model, although the predicted $\pi_{\rm E}$ distribution significantly
differed from the previous one. Further, we applied the model to calculate the
lens-source proper motion $\mu_{\rm rel}$ prior for OGLE-2011-BLG-0950, the
only ambiguous event suffering from degeneracy between the planetary and
binary solutions out of the 29 events in the Suzuki et al. (2016) combined
sample. Our prior probability distribution for $\mu_{\rm rel}$ indicates that
the binary solution is more probable than the planetary solution.

This paper is organized as follows. In Section 2, we describe our parametric
Galactic model and introduce all of the 40–48 fit parameters. Fitting to the
Gaia DR2 data is conducted to determine the 10 parameters for the disk
velocity model in Section 3. Fitting to the other data, which is toward the
bulge sky, is separately conducted in Section 4 to determine the 26–34
parameters for the bulge density and velocity models, as well as the 4
parameters for the IMF model, where the determined disk velocity model is used
in the fit. In Section 5, we discuss the determined parameters by comparing
them with previous studies. We apply the developed Galactic models to
microlensing analysis in Section 6. Section 7 presents the summary and
conclusion.

##  2 Model Parameterization

In this section, we describe our parametric Galactic model with a barred bulge
and a multi-component of the stellar thin and thick disks. We disregard a
stellar halo in this study because it rarely contributes to the considered
data, which are mostly toward the Galactic bulge. We refer to a bar with a
rigid-body rotation as the ‘bulge’ while an inner part of our Galaxy as the
‘bulge region,’ which means stars in the bulge region refer to both the bulge
and disk stars present there. These Galactic structures consist of stars with
mass functions (Section 2.1) and have the density and velocity distributions
as functions of the Galactocentric coordinate ($x,y,z$) or Galactocentric
cylindrical coordinate ($R,\phi,z$). We use different density and velocity
distributions for each component, where we use $\rho_{\rm d}$ and $v_{\rm d}$
to refer to the disk density and velocity, respectively (Section 2.2), and
$\rho_{\rm B}$ and $v_{\rm B}$ to refer to the bulge density and velocity,
respectively (Section 2.3).

Here we introduce 4 fit parameters to model the stellar mass function, 10 fit
parameters to model the $v_{\rm d}$ distribution, 7–15 fit parameters to model
the $\rho_{\rm B}$ distribution, and 19 fit parameters to model the $v_{\rm
B}$ distribution, respectively. Whereas the $\rho_{\rm d}$ model does not
contain any fit parameters, two options of flat-scale and linear-scale height
models are introduced for the thin disk. The 10 parameters for the $v_{\rm d}$
distribution are determined by fitting to the Gaia DR2 data (Gaia
Collaboration et al., 2018) in Section 3, while the other 30–38 parameters are
determined in Section 4 by fitting to the other data toward the bulge sky,
namely the OGLE-III red clump star count (Nataf et al., 2013), VIRAC proper
motion data (Smith et al., 2018; Clarke et al., 2019), BRAVA radial velocity
data (Rich et al., 2007; Kunder et al., 2012), and OGLE-IV star count and
microlensing event data (Mróz et al., 2017, 2019).

The OGLE-III red clump star count data from Nataf et al. (2013) is sensitive
to the Sun’s Galactic position and the bar angle. In fact, Cao et al. (2013)
used the same dataset to measure the distance to the Galactic center of
8.1–8.2 kpc and the bar angle of 27–32 deg. However, because the number of fit
parameters is already large, 40 to 48, we fix these two parameters to mitigate
degeneracies among fit parameters. For the solar position in our Galaxy, we
use $R_{\odot}=8160~{}{\rm pc}$, which is consistent with a recent precise
measurement of $8178\pm 13_{\rm stat.}\pm 22_{\rm sys.}~{}{\rm pc}$ on the
distance to Sgr A* (Gravity Collaboration et al., 2019), and
$z_{\odot}=25~{}{\rm pc}$ (Bland-Hawthorn & Gerhard, 2016). We use
$\alpha_{\rm bar}=27^{\circ}$ for the bar angle (Wegg & Gerhard, 2013; Bland-
Hawthorn & Gerhard, 2016). The $R_{\odot}$ and $\alpha_{\rm bar}$ values are
consistent with the OGLE-III red clump star count data because these two
values are both within the above respective ranges derived by Cao et al.
(2013). Thus, this choice generates no tension.

We also fix the solar velocity to be
$(v_{\odot,x},v_{\odot,y},v_{\odot,z})=(-10,243,7)$ km/s (Bland-Hawthorn &
Gerhard, 2016). The tangential velocity $v_{\odot,y}=243$ km/s is marginally
lower than $248\pm 3$ km/s suggested by Bland-Hawthorn & Gerhard (2016) (a
recent thorough review paper on our Galaxy) because we found a better
agreement to the VIRAC proper motion data with 243 km/s. This is consistent
with Clarke et al. (2019), who found that the dynamical model of Portail et
al. (2017), with a slower tangential velocity of $v_{\odot,y}=245$ km/s, gave
an improved match to the VIRAC proper motion data, compared to the original
value of $v_{\odot,y}=250.24$ km/s used by Portail et al. (2017).

###  2.1 Stellar Mass Function

The stellar mass function refers to a present-day mass function, which can be
calculated by a combination of IMF, star formation rate, and initial–final
mass relationships for remnants, i.e., for white dwarfs, neutron stars, and
black holes. Here, we introduce the IMF, star formation rate, and
initial–final mass relationships used in this paper.

We use a broken-power law form IMF given by

$\frac{dN}{dM}\propto\begin{cases}M^{-\alpha_{\rm hm}}&\text{ when $M_{\rm
br}<M<120\,M_{\odot}$}\\\ M^{-\alpha_{\rm ms}}&\text{ when
$0.08\,M_{\odot}\leq M<M_{\rm br}$}\\\ M^{-\alpha_{\rm bd}}&\text{ when
$10^{-3}\,M_{\odot}\leq M<0.08\,M_{\odot}$}\ .\end{cases}$ (4)

We use four parameters, the three slopes $(\alpha_{\rm hm},\alpha_{\rm
ms},\alpha_{\rm bd})$ and a break mass $M_{\rm br}$, as fit parameters in
Section 4, and fit them to the data including the OGLE-IV microlensing $t_{\rm
E}$ distribution, which has a sensitivity to objects with brown dwarf mass to
black hole mass. By contrast, we use the local IMF by Kroupa (2001),
$(\alpha_{\rm hm},\alpha_{\rm ms},\alpha_{\rm bd})$ = (2.30, 1.30, 0.30) with
$M_{\rm br}=0.50\,M_{\odot}$, in the fit to the Gaia DR2 data, in Section 3,
to derive the values of the local density for red giants, which are listed in
Table 1. We selected different IMF values because the Gaia data are dominated
by local stars, whereas the data used in Section 4 are dominated by bulge
stars.

Further, a star formation rate $\propto\exp[\frac{T}{7{\rm Gyr}}]$ (Bovy,
2017) is applied to the thin disk star, where $T$ is stellar age, while we
assume a mono-age thick disk at $T=12~{}$Gyr (Robin et al., 2014) and the
Gaussian distribution of $9\pm 1~{}$Gyr (Koshimoto et al., 2020) for the bulge
stars. The PARSEC isochrone models (Bressan et al., 2012; Chen et al., 2014;
Tang et al., 2014) are used for the stellar lifetime of a given initial mass.

For initial–final mass relationships of remnants, we apply the model used by
Lam et al. (2020), where a linear relation of Kalirai et al. (2008), $M_{\rm
WD}=0.109M_{\rm ini}+0.394M_{\odot}$, is used for white dwarfs and a modified
probabilistic relation based on Raithel et al. (2018) is used for neutron
stars and black holes. See Appendix C of Lam et al. (2020) and references
therein for more details. Further, we add a birth kick velocity in a random
direction of 350 km/s or 100 km/s to the velocity of the progenitor star for a
neutron star or a black hole, respectively (Lam et al., 2020).

###  2.2 Galactic Disk Model

In this subsection, we describe the density and velocity models of a thin and
a thick disk, where the thin disk is further divided into seven components
depending on ages following the classification by Robin et al. (2003). These
seven plus one components have different scale height values and velocity
dispersions. The velocity dispersions differ from each other depending on the
10 fit parameters introduced in the $v_{\rm d}$ model and determined later in
Section 3. The disk model contributes to both fits in Sections 3 and 4.

####  2.2.1 Disk density ($\rho_{\rm d}$ model)

A thin and a thick disk are considered with density profiles of

$\displaystyle\rho_{\rm d}^{\rm thin}(R,z)$
$\displaystyle=\begin{cases}\rho^{\rm thin}_{\rm d,\odot}\,\frac{z_{\rm
d,\odot}^{\rm thin}}{z_{\rm d}^{\rm
thin}(R)}\,\exp{\left(-\frac{R-R_{\odot}}{R_{\rm d}^{\rm thin}}\right)}\,{\rm
sech}^{2}{\left(-\frac{|z|}{z_{\rm d}^{\rm thin}(R)}\right)}&\text{ if
$R>R_{\rm d,break}$}\\\ \rho^{\rm thin}_{\rm d,\odot}\,\frac{z_{\rm
d,\odot}^{\rm thin}}{z_{\rm d}^{\rm thin}(R)}\,\exp{\left(-\frac{R_{\rm
d,break}-R_{\odot}}{R_{\rm d}^{\rm thin}}\right)}\,{\rm
sech}^{2}{\left(-\frac{|z|}{z_{\rm d}^{\rm thin}(R)}\right)}&\text{ if $R\leq
R_{\rm d,break}$}\end{cases}$ (5) $\displaystyle\rho_{\rm d}^{\rm thick}(R,z)$
$\displaystyle=\begin{cases}\rho^{\rm thick}_{\rm
d,\odot}\exp{\left(-\frac{R-R_{\odot}}{R_{\rm d}^{\rm thick}}\right)}\,{\rm
exp}{\left(-\frac{|z|}{z_{\rm d,\odot}^{\rm thick}}\right)}&\text{ if
$R>R_{\rm d,break}$}\\\ \rho^{\rm thick}_{\rm d,\odot}\exp{\left(-\frac{R_{\rm
d,break}-R_{\odot}}{R_{\rm d}^{\rm thick}}\right)}\,{\rm
exp}{\left(-\frac{|z|}{z_{\rm d,\odot}^{\rm thick}}\right)}&\text{ if $R\leq
R_{\rm d,break}$,}\end{cases}$ (6)

where $\rho_{\rm d,\odot}$ is the stellar density in the solar neighborhood,
hereafter referred to as the local (stellar) density, $R_{\rm d}$ is the disk
scale length, and $z_{\rm d}$ is the disk scale height.

The disk profile in the inner Galaxy region is still uncertain. A constant
surface density model is adopted within $R<R_{\rm d,break}$ following the
findings of Portail et al. (2017) who developed an $N$-body dynamical model
reproducing extensive photometric and kinematic data across our Galaxy. In
their model, the disk has an exponential feature in the outer part of our
Galaxy ($\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}5$ kpc) while it
has a flat surface density feature in the inner part ($\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}5$ kpc). We use $R_{\rm
d,break}=5.3~{}{\rm kpc}$ for our model.

For the thin disk scale height $z_{\rm d}^{\rm thin}$, two options are
considered: a flat- and a linear-scale height model. We adopt

$\displaystyle z_{\rm d}^{\rm thin}(R)$ $\displaystyle=\begin{cases}z_{\rm
d,\odot}^{\rm thin}-(z_{\rm d,\odot}^{\rm thin}-z_{\rm d,4.5}^{\rm
thin})\frac{R_{\odot}-R}{R_{\odot}-4.5{\rm kpc}}&\text{ if $R>4.5$ kpc}\\\
z_{\rm d,4.5}^{\rm thin}&\text{ if $R\leq 4.5$ kpc}\end{cases}$ (7)

for the linear-scale height model with $z_{\rm d,4.5}^{\rm thin}=0.6\,z_{\rm
d,\odot}^{\rm thin}$ (Wegg et al., 2016) while for the flat-scale height model
$z_{\rm d}^{\rm thin}(R)=z_{\rm d,\odot}^{\rm thin}$ is used. Two options have
been considered as we were motivated by a similar attempt of Wegg et al.
(2016) who considered the linear-scale height model to smoothly connect the
disk to the inner thin long bar with a scale height of $\sim 200$ pc suggested
by Wegg et al. (2015). As described in Section 2.3.1, adding the long bar
component to our bulge model did not result in a significant improvement.
However, whether we used the flat- or linear-scale height model in the disk
model did affect $\chi^{2}$ or $\tilde{\chi}^{2}$ values in the fit.

Table 1: Scale lengths, scale heights and local densities for thin and thick
disks.

| Age $T$ | $R_{\rm d}$ | $z_{\rm d,\odot}$ |  $z_{\rm d,4.5}$ aaFor the linear scale height model. |  $\rho_{\rm d,\odot}^{\rm MS}$ bbLocal densities are given for main sequence stars (MS, including brown dwarfs), white dwarfs (WD), and red giants (RG). In the fit in Section 4, the values are updated with a given IMF. |  $\rho_{\rm d,\odot}^{\rm WD}$ bbLocal densities are given for main sequence stars (MS, including brown dwarfs), white dwarfs (WD), and red giants (RG). In the fit in Section 4, the values are updated with a given IMF. |  $n_{\rm d,\odot}^{\rm RG}$ bbLocal densities are given for main sequence stars (MS, including brown dwarfs), white dwarfs (WD), and red giants (RG). In the fit in Section 4, the values are updated with a given IMF.
---|---|---|---|---|---|---|---
| [Gyr] | [pc] | [pc] | [pc] | [$M_{\odot}$ pc-3] | [$M_{\odot}$ pc-3] | [pc-3]
Thin disk | 0 - 0.15 | 5000 | 61 | 36 | $5.1\times 10^{-3}$ | $5.5\times 10^{-5}$ | $6.9\times 10^{-6}$
| 0.15 - 1 | 2600 | 141 | 85 | $5.0\times 10^{-3}$ | $2.2\times 10^{-4}$ | $3.3\times 10^{-5}$
| 1 - 2 | 2600 | 224 | 134 | $3.8\times 10^{-3}$ | $2.9\times 10^{-4}$ | $4.2\times 10^{-5}$
| 2 - 3 | 2600 | 292 | 175 | $3.2\times 10^{-3}$ | $3.3\times 10^{-4}$ | $2.1\times 10^{-5}$
| 3 - 5 | 2600 | 372 | 223 | $5.9\times 10^{-3}$ | $7.8\times 10^{-4}$ | $6.5\times 10^{-5}$
| 5 - 7 | 2600 | 440 | 264 | $6.3\times 10^{-3}$ | $1.0\times 10^{-3}$ | $6.1\times 10^{-5}$
| 7 - 10 | 2600 | 445 | 267 | $1.3\times 10^{-2}$ | $2.4\times 10^{-3}$ | $1.3\times 10^{-4}$
Sum/Mean |  |  | 329 | 197 | $4.2\times 10^{-2}$ | $5.1\times 10^{-3}$ | $3.6\times 10^{-4}$
Thick disk | 12 | 2200 | 903 |  | $1.7\times 10^{-3}$ | $4.4\times 10^{-4}$ | $9.1\times 10^{-6}$

Note. — $R_{\rm d}$ refers to $R_{\rm d}^{\rm thin}$ in the lines for thin
disk, and $R_{\rm d}^{\rm thick}$ in the line for thick disk. The same is true
for $z_{\rm d,\odot}$, $z_{\rm d,4.5}$, $\rho_{\rm d,\odot}$, and $n_{\rm
d,\odot}$.

The thin disk is divided into 7 components depending on the age ($T$)
following the Besançon model (Robin et al., 2003, 2012, 2017). Table 1 lists
the scale length and height values for each component of the thin and thick
disks. The scale height for each component is calculated using the age-scale
height relation for the axis ratio $\epsilon$, that is $\epsilon={\rm
Min}[0.0791,0.104\sqrt{(T/{\rm Gyr}+0.1)/10.1}]$ (Sharma et al., 2014).
Because $\epsilon$ is designed as the ratio of the scale height to the scale
length in a disk with Einasto laws (Einasto, 1979), we used the surface-to-
volume density ratio in the Einasto disk calculated with $\epsilon$ to derive
the $z_{\rm d,\odot}$ value for each component.

Table 1 also lists the local stellar density values for main sequence stars,
white dwarfs, and red giants for each disk component. The total thin disk
local mass density is normalized by $4.2\times 10^{-2}\,M_{\odot}$ pc-3 for
main sequence stars (Bovy, 2017), while the thick disk local mass density is
normalized so that its local contribution to the thin disk becomes $4\%$
(Bland-Hawthorn & Gerhard, 2016).

The local density value for stars at each evolutional stage in each disk
component is calculated using the local IMF by Kroupa (2001) combined with the
star formation rate and the initial–final mass relationship described in
Section 2.1. The red giants are selected based on their absolute magnitude
${\cal M}_{G}<3.9$ mag and intrinsic color $(G_{\rm BP}-G_{\rm RP})_{0}>0.95$
in the Gaia bands using the PARSEC isochrones (Bressan et al., 2012; Chen et
al., 2014; Tang et al., 2014). This criteria for red giants is same as that
used for the selection of the giant sample by Gaia Collaboration et al.
(2018), and the $n_{\rm d,\odot}^{\rm RG}$ values listed in the table are used
in the fit to the Gaia DR2 data in Section 3.

An important note here is that the local density values are updated with a
given IMF and are different in the fit in Section 4. This apparent
inconsistent treatment can be partially justified by considering that our IMF
presented in Section 4 is for the disk stars located in the inner Galaxy
region ($R\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}4$ kpc)
contributing to the data used in Section 4, and the Kroupa (2001) IMF is
appropriate for the nearby stars ($R\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}4$ kpc) contributing to the Gaia
data used in Section 3.

The local surface density with the Kroupa (2001) IMF is $30.01~{}M_{\odot}$
pc-2 for main sequence stars, which is compared with the measurements of
$23.0\pm 1.5~{}M_{\odot}$ pc-2 (Bovy, 2017) and $28.2\pm 2.7~{}M_{\odot}$ pc-2
(McKee et al., 2015). The local surface density of white dwarfs is
$4.87~{}M_{\odot}$ pc-2, and this is consistent with measurements of $5\pm
1~{}M_{\odot}$ pc-2 (Bovy, 2017) and $4.9\pm 0.6~{}M_{\odot}$ pc-2 (McKee et
al., 2015). The local number density of red giants is $3.6\times 10^{-4}$
pc-3, and this fairly well agrees with the measurement by Bovy (2017) of
$(3.9\pm 0.1)\times 10^{-4}$ pc-3.

####  2.2.2 Disk kinematics ($v_{\rm d}$ model)

Our work is highly motivated by Gaia DR2 (Gaia Collaboration et al., 2018) in
which skewed distributions for the azimuthal velocity $v_{\phi}$ and clear
dependencies of the velocity dispersion and mean azimuthal velocity of disk
stars on their location are shown. Such disk kinematic structures have not
been included in most of the models used for microlensing analysis like the
three models (Sumi et al., 2011; Bennett et al., 2014; Zhu et al., 2017) used
in Koshimoto & Bennett (2020). To include those dependencies as a function of
the Galactocentric cylindrical coordinate ($R,\phi,z$) in addition to a skewed
$v_{\phi}$ distribution, we follow the parameterization of a disk velocity
model by Sharma et al. (2014).

We assume that the Galaxy is in a dynamical equilibrium, and use a modified
Shu distribution function (DF) model developed by Schönrich & Binney (2012)
and Sharma & Bland-Hawthorn (2013) to represent the distribution of disk
azimuthal velocity $v_{\rm d,\phi}$. Gaussian velocity models are used for
$v_{\rm d,\phi}$ distributions in other Galactic models for microlens analysis
(Sumi et al., 2011; Bennett et al., 2014; Zhu et al., 2017; Jung et al., 2018)
and in the Besançon model (Robin et al., 2003, 2012, 2017). However, a real
$v_{\phi}$ distribution is highly skewed to low $v_{\phi}$ (e.g., Nordström et
al., 2004; Gaia Collaboration et al., 2018), and the Shu DF (Shu, 1969)
provides a much better approximation for it (Binney & Tremaine, 2008; Sharma
et al., 2014; Bland-Hawthorn & Gerhard, 2016).

We introduce the guiding-center radius $R_{g}$ as the radius of a circular
orbit with specific angular momentum $L_{z}$, i.e., $R_{g}=L_{z}/v_{c}$, where
$v_{c}$ is the circular velocity (Binney & Tremaine, 2008). The modified Shu
DF model provides a joint distribution of the Galactocentric radius $R$ and
$R_{g}$,

$\displaystyle
P(R,R_{g})=\frac{(2\pi)^{2}\Sigma(R_{g})}{g\left(\frac{1}{2a^{2}}\right)}\exp{\left[\frac{2\ln(R_{g}/R)+1-R_{g}^{2}/R^{2}}{2a^{2}}\right]},$
(8)

where $a=\sigma_{v_{{\rm d},R}}(R_{g})/v_{c}$,
$g(c)=\frac{e^{c}\Gamma(c-1/2)}{2c^{c-1/2}}$ with the velocity dispersion
along radial direction $\sigma_{v_{{\rm
d},R}}(R)=\sigma_{R,\odot}\exp[{-\frac{R-R_{\odot}}{R_{\sigma_{R}}}}]$ and the
Gamma function $\Gamma(x)$. $R_{\sigma_{R}}$ is the scale length for the
$\sigma_{v_{{\rm d},R}}(R)$ distribution given by Eqs. (11)-(12) below.
$\Sigma(R_{g})$ is a function that controls disk surface density, and we use
an empirical formula proposed by Sharma & Bland-Hawthorn (2013),

$\displaystyle\Sigma(R_{g})=\frac{e^{-R_{g}/R_{\rm d}}}{2\pi R_{\rm
d}^{2}}-\frac{c_{3}a_{0}^{c_{4}}}{R_{\rm d}^{2}}\times
s\left(\frac{R_{g}}{c_{1}R_{\rm d}(1+q/c_{2})}\right),$ (9)

where $s(x)=31.53e^{-x/0.2743}((x/0.6719)^{2}-1)$ and $q=R_{\rm
d}/R_{\sigma_{R}}$. We select
$(c_{1},c_{2},c_{3},c_{4})=(3.822,0.524,0.00567,2.13)$ from Table 1 of Sharma
& Bland-Hawthorn (2013) that is for the rising rotation curve of
$v_{c}\propto(R/R_{\rm d})^{0.2}$. This is because we use a similar rising
rotation curve of $v_{c}(R)$ from Bland-Hawthorn & Gerhard (2016), which comes
from the $N$-body dynamical model of Portail et al. (2015) for $R_{\rm d}=2.6$
kpc.

A conditional probability for $R_{g}$ given $R$, $P(R_{g}|R)$, which is
calculated using Eq. (8), is used to model the disk azimuthal velocity $v_{\rm
d,\phi}$ distribution through the relation between $v_{\rm d,\phi}$ and
$R_{g}$,

$\displaystyle v_{\rm d,\phi}(R,z)$ $\displaystyle=v_{c}(R_{g},z)R_{g}/R$
$\displaystyle=\frac{v_{c}(R_{g})R_{g}/R}{1+0.0374|z/{\rm kpc}|^{1.34}},$ (10)

where we apply $v_{c}(R,z)=v_{c}(R)(1+0.0374|z/{\rm kpc}|^{1.34})^{-1}$
(Sharma et al., 2014) for the vertical dependency of $v_{c}$. Again, the
rotation curve from Bland-Hawthorn & Gerhard (2016) is used for $v_{c}(R)$.

For the disk velocity along radial ($v_{{\rm d},R}$) and vertical ($v_{{\rm
d},z}$) directions, we use the Gaussian distribution with mean velocity of 0
(i.e. dynamical equilibrium) with the velocity dispersion given by

$\displaystyle\sigma_{v_{{\rm d},i}}^{\rm thin}(R)=\sigma_{i,\odot}^{\rm
thin}\left(\frac{T+T_{\rm min}}{T_{\rm max}+T_{\rm
min}}\right)^{\beta_{i}}\exp\left[-\frac{R-R_{\odot}}{R_{\sigma_{i}}^{\rm
thin}}\right]$ (11)

for the thin disk and

$\displaystyle\sigma_{v_{{\rm d},i}}^{\rm thick}(R)=\sigma_{i,\odot}^{\rm
thick}\exp\left[-\frac{R-R_{\odot}}{R_{\sigma_{i}}^{\rm thick}}\right]$ (12)

for the thick disk, where $i$ takes $R$ or $z$, and we use $T_{\rm min}=0.01$
Gyr and $T_{\rm max}=10$ Gyr in this study. Further, we introduce the
dependence on stellar age $T$ for the thin disk to include the age-velocity
dispersion relation owing to secular heating in the disk. With this formula,
the local velocity dispersion value calculated for thin disk,
$\sigma_{i,\odot}^{\rm thin}$, is for stars with $T=10$ Gyr.

In Section 3.3, we investigate acceptable combinations of the following 10 fit
parameters by comparing with the data from the giant sample by Gaia
Collaboration et al. (2018); the local velocity dispersion values,
$\sigma_{i,\odot}^{\rm thin}$ and $\sigma_{i,\odot}^{\rm thick}$, slope of
age-velocity dispersion relation, $\beta_{i}$, and scale lengths of the
velocity dispersion distribution, $R_{\sigma_{i}}^{\rm thin}$ and
$R_{\sigma_{i}}^{\rm thick}$, where $i$ takes $R$ or $z$. Note that additional
parameters are not needed to represent the distribution of velocity dispersion
along the azimuthal direction, such as the above parameters with $i=\phi$,
because Eq. (8) naturally relates the $v_{{\rm d},\phi}$ distribution to the
$v_{{\rm d},R}$ distribution.

###  2.3 Barred Bulge Model

Compared to the disk model, an analytical approximated expression for the
bulge dynamical model is less developed because of its difficulty in the
treatment of a non-axisymmetric property of the bar. An $N$-body model is
dynamically correct; however, fitting the model to the observational data is
difficult. A probable optimal technique is a made-to-measure method (Syer &
Tremaine, 1996), where weights of particles are updated such that observables
of the model match a given dataset during the simulation. Portail et al.
(2017) developed an $N$-body dynamical model matching extensive photometric
and kinematic data across our Galaxy using the made-to-measure method.
However, such a dynamical simulation is beyond the scope of this study because
we aim to develop a parametric Galactic model, which can be easily
implemented, reproduced, and updated by anybody. Although some studies
developed parametric models for the bar (Dwek et al., 1995; Rattenbury et al.,
2007b; Robin et al., 2012; Cao et al., 2013), they lack constraints from some
recent data such as the one by Mróz et al. (2019), who performed the largest
statistical study for single-lens microlensing events through the OGLE-IV
Galactic bulge survey, which is especially important for microlensing studies.

In this subsection, we describe our parameterization for the bulge density
($\rho_{\rm B}$) and velocity ($v_{\rm B}$) models. We consider a total of
four different shapes for the $\rho_{\rm B}$ model: two ‘one-component’ models
and two ‘two-components’ models, in which 7 and 15 fit parameters are
introduced, respectively. The one-component model is designed following the
parameterization used in the previous studies (Dwek et al., 1995; Rattenbury
et al., 2007b; Robin et al., 2012; Cao et al., 2013), while the two-components
model is designed to express the X-shape structure (Nataf et al., 2010), which
is not considered in the previous parametric models. We consider a bar’s
rigid-body rotation and a streaming motion in the $v_{\rm B}$ model and
introduce 19 fit parameters to model it. Note that the bulge model contributes
to the fits in Section 4 rather than those in Section 3.

####  2.3.1 Bulge density ($\rho_{\rm B}$ model)

Our one-component bulge model follows the parameterization by Robin et al.
(2012), in which we consider each of E (exponential) and G (Gaussian) models
given by

$\displaystyle\rho_{\rm B}^{(i)}=\rho_{\rm
0,B}\,f_{i}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})\,{\rm
Cut}\left[\frac{R-R_{\rm c}}{0.5\,{\rm kpc}}\right]\hskip 17.07182pt(i={\rm
E,G}),$ (13)

where ${\rm Cut}(x)$ is a cut-off function given by

$\displaystyle{\rm Cut}(x)=\begin{cases}\exp(-x^{2})&\text{ if $x>0$}\\\
1&\text{ if $x\leq 0$}\end{cases}$ (14)

and $R_{\rm c}$ is the cut-off radius. The two functions, $f_{\rm
E}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})$ and $f_{\rm
G}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})$, are defined as

$\displaystyle f_{\rm E}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})$
$\displaystyle=\exp[-r_{s}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})],$
$\displaystyle f_{\rm G}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})$
$\displaystyle=\exp[-0.5r_{s}^{2}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s}})]$
(15)

with

$\displaystyle
r_{s}(x^{\prime},y^{\prime},z^{\prime};{\bm{p}_{r_{s}}})=\left\\{\left[\left(\frac{x^{\prime}}{x_{0}}\right)^{C_{\perp}}+\left(\frac{y^{\prime}}{y_{0}}\right)^{C_{\perp}}\right]^{C_{\parallel}/C_{\perp}}+\left(\frac{z^{\prime}}{z_{0}}\right)^{C_{\parallel}}\right\\}^{1/C_{\parallel}}$
(16)

and $\bm{p}_{r_{s}}=(x_{0},y_{0},z_{0},C_{\perp},C_{\parallel})$. We use
$(x^{\prime},y^{\prime},z^{\prime})$ to refer to a Galactocentric coordinate
system rotated around the $z$-axis by an angle $\alpha_{\rm bar}$ such that
the $x^{\prime}$ axis is aligned with the major axis of the Galactic bar,
where $\alpha_{\rm bar}=27^{\circ}$ is applied as the bar angle. The
parameters $(x_{0},y_{0},z_{0})$ are the scale lengths along
$(x^{\prime},y^{\prime},z^{\prime})$ axes, and $C_{\perp}$ and $C_{\parallel}$
allow the bar to take various shapes (Robin et al., 2012).

Motivated by the X-shape structure confirmed both observationally and
dynamically (McWilliam & Zoccali, 2010; Nataf et al., 2010; Wegg & Gerhard,
2013), we consider a two-components model, $\rho_{\rm B}=\rho_{\rm
B}^{(i)}+\rho_{\rm X}^{(j)}$, where $\rho_{\rm B}^{(i)}$ is given by Eq. (13)
and

$\displaystyle\rho_{\rm X}^{(j)}=\rho_{\rm 0,X}\left[f_{j}(x^{\prime}-b_{\rm
X}z^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s},{\rm
X}})+f_{j}(x^{\prime}+b_{\rm
X}z^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s},{\rm X}})\right]{\rm
Cut}\left[\frac{R-R_{\rm c,X}}{0.5\,{\rm kpc}}\right]\ \hskip
17.07182pt(j={\rm E,G}).$ (17)

The parameter $b_{\rm X}$ controls the slope of an X-shape, and we use another
parameter set, $\bm{p}_{r_{s},{\rm X}}=(x_{\rm 0,X},y_{\rm 0,X},z_{\rm
0,X},C_{\rm\perp,X},C_{\rm\parallel,X})$, which is different from the
$\bm{p}_{r_{s}}$ for the first component given by Eq. (13). Note that the
X-shape structure with this expression is centered on the Galactic center,
although Portail et al. (2015) found a slightly off-centered X-shape structure
in their $N$-body dynamical model. We considered all four combinations of
$i={\rm E,G}$ and $j={\rm E,G}$, and found no significant difference among
those combinations with respect to agreement with the fitted data. Hence,
herein, we present results of two combinations of $(i,j)$ = (E, E) and (G, G)
among the four. Hereafter, we refer to these models as the E+EX and G+GX
models, respectively.

Furthermore, other structures known in the bulge region are locally
significant but not captured by either of the four $\rho_{\rm B}$ models. For
example, a nuclear stellar disk exists in a central sub-kpc region (Launhardt
et al., 2002; Nishiyama et al., 2013; Portail et al., 2017). The scale height
and the outer edge of the nuclear stellar disk are $\sim 45~{}{\rm pc}$
(Nishiyama et al., 2013) and $\sim 230~{}{\rm pc}$ (Bland-Hawthorn & Gerhard,
2016), respectively, which indicates no influence on the used data in this
study ranged in $|b|\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}2^{\circ}$. A long bar component
was found to be distributed in the outer bulge region along the major axis by
Wegg et al. (2015); hence, it was added in our model to observe its effect,
but no significant improvement was observed regarding the
${\tilde{\chi}^{2}_{\rm sum}}$ value defined in Section 4.3. This is probably
because the used data lacks the sky region in $|l|>10^{\circ}$, where the long
bar component becomes prominent. Therefore, we consider each of the E, G, E+EX
and G+GX models with no additional components.

The 7 fit parameters for the E and G models are:
$x_{0},y_{0},z_{0},C_{\perp},C_{\parallel},R_{c}$, and $\rho_{\rm 0,B}$. In
the E+EX and G+GX models, the additional 8 fit parameters are: $x_{\rm
0,X},y_{\rm 0,X},z_{\rm 0,X},C_{\rm\perp,X},C_{\rm\parallel,X},R_{\rm c,X}$,
$b_{\rm X}$, and $f_{\rm 0,X}$, where $f_{\rm 0,X}\equiv\rho_{\rm
0,X}/\rho_{\rm 0,B}$.

####  2.3.2 Bulge kinematics ($v_{\rm B}$ model)

For the velocity distribution of the bulge star, we use the Gaussian
distribution with a mean velocity and velocity dispersion, both varying as a
function of $(x^{\prime},y^{\prime},z^{\prime})$. The mean velocity is
calculated by combining the rigid-body rotation of the bar and streaming
motion along the bar. We denote the angular velocity of the rigid-body
rotation or the bar pattern speed by $\Omega_{\rm p}$ and consider the
streaming motion along the major axis as

$\displaystyle v_{x^{\prime}}^{\rm str}(y^{\prime})=v_{0}^{\rm
str}(1-\exp[-(y^{\prime}/y_{0}^{\rm str})^{2}]),$ (18)

where $v_{0}^{\rm str}$ is the streaming velocity at $y^{\prime}\gg y_{0}^{\rm
str}$, and $y_{0}^{\rm str}$ is the scale length along $y^{\prime}$ axis. This
form of distribution for the streaming motion is motivated by the bottom
panels in Figure 14 of Sanders et al. (2019a) that show the
$v_{x^{\prime}}^{\rm str}$ of their dynamical model, which represents our
Galaxy, increasing from $y^{\prime}=0$ to $|y^{\prime}|>0$ along $y^{\prime}$
axis.

For the velocity dispersion, we use

$\displaystyle\sigma_{v_{\rm
B},i}(x^{\prime},y^{\prime},z^{\prime})=\sigma_{v_{\rm B},i,0}+\sigma_{v_{\rm
B},i,1}\,f_{\rm
E}(x^{\prime},y^{\prime},z^{\prime};\bm{p}_{r_{s},\sigma_{i}})\hskip
22.76228pt(i=x^{\prime},y^{\prime},z^{\prime}),$ (19)

where we apply
$\bm{p}_{r_{s},\sigma_{x^{\prime}}}=\bm{p}_{r_{s},\sigma_{y^{\prime}}}\neq\bm{p}_{r_{s},\sigma_{z^{\prime}}}$
and denote the parameter set for $\sigma_{x^{\prime}}$ and
$\sigma_{y^{\prime}}$ by $\bm{p}_{r_{s},\sigma_{R}}$, i.e.,
$\bm{p}_{r_{s},\sigma_{R}}=\bm{p}_{r_{s},\sigma_{x^{\prime}}}=\bm{p}_{r_{s},\sigma_{y^{\prime}}}=(x_{\rm
0,\sigma_{R}},y_{\rm 0,\sigma_{R}},z_{\rm
0,\sigma_{R}},C_{\rm\perp,\sigma_{R}},C_{\rm\parallel,\sigma_{R}})$. The
constant $\sigma_{v_{\rm B},i,0}$ provides a minimum value of $\sigma_{v_{\rm
B},i}$ while $\sigma_{v_{\rm B},i,0}+\sigma_{v_{\rm B},i,1}$ provides the
maximum $\sigma_{v_{\rm B},i}$ value at the Galactic center.

This model is not dynamically consistent with the density model described in
Section 2.3.1; however, the profile of Eq. (19) peaking at the Galactic center
and gradually decreasing as it goes around is motivated by the velocity
dispersion field of the dynamical model of our Galaxy by Sanders et al.
(2019a) illustrated in their Figure 15.

There are a total of 19 fit parameters for the velocity model; $\Omega_{\rm
p},v_{0}^{\rm str},y_{0}^{\rm str},\sigma_{v_{\rm B},i,0},\sigma_{v_{\rm
B},i,1}(i=x^{\prime},y^{\prime},z^{\prime})$, $x_{\rm 0,\sigma_{R}},y_{\rm
0,\sigma_{R}},z_{\rm
0,\sigma_{R}},C_{\rm\perp,\sigma_{R}},C_{\rm\parallel,\sigma_{R}}$, $x_{\rm
0,\sigma_{z^{\prime}}},y_{\rm 0,\sigma_{z^{\prime}}},z_{\rm
0,\sigma_{z^{\prime}}},C_{\rm\perp,\sigma_{z^{\prime}}},$ and
$C_{\rm\parallel,\sigma_{z^{\prime}}}$.

##  3 Fitting for the Disk Velocity Parameters

In this section, we determine the 10 fit parameters for the $v_{\rm d}$ model
by fitting our disk model to the spatial distribution of the median velocity
and the velocity dispersion for giant stars in the Gaia DR2 (Gaia
Collaboration et al., 2018), where the data are given in grids of 200 pc by
200 pc in $(R,z)$. Bulge stars rarely contribute to the Gaia data and the
barred bulge model is not used in the fit in this section. The fit is
conducted through a grid search. The disk model with the determined parameters
in this section is used for the other fits in Section 4.

###  3.1 Gaia DR2 Velocity Data

We use the median velocity and velocity dispersion distributions of the red
giant sample consisting of 3,153,160 sources by Gaia Collaboration et al.
(2018) as a function of the Galactocentric radius $R$ and the height from the
Galactic plane $z$. The stars in this sample are selected based on their
absolute magnitude ${\cal M}_{G}<3.9$ mag and intrinsic color $(G_{\rm
BP}-G_{\rm RP})_{0}>0.95$ in the Gaia bands. This is the same data plotted in
Figure 11 of Gaia Collaboration et al. (2018), and we obtained the data
through private communication with the lead author, D. Katz. The medians and
dispersions are given in grids of 200 pc by 200 pc in $(R,z)$ over
$3340~{}{\rm pc}<R<13340~{}{\rm pc}$, $-3400~{}{\rm pc}<z<3400~{}{\rm pc}$.
The center of $(i_{R},i_{z})$th grid is
$(R,z)=(200\,i_{R}+3440,200\,i_{z}-3300)~{}{\rm pc}$. The total grids in the
range of $R$ and $z$ are $50\times 34=1700$; however, we do not use the grids
with $<30$ stars contributing to the statistics, which results in 1207 grids
being available.

The medians of the $v_{R}$ and $v_{z}$ distribution are assumed as 0; hence,
only the data of the median of the azimuthal velocity $\overline{v_{\phi}}$
and velocity dispersions along the three axes, $\sigma_{R},\sigma_{\phi}$, and
$\sigma_{z}$, are used. The data distribution for these four parameters is
plotted in the far-left panels in Fig. 1.

###  3.2 Definition of Goodness of Fit

The far-right column of Table 1 lists the local number density of red giants
in our model, which is calculated using the same criteria as for the Gaia DR2
giant sample. We use these $n_{\rm d,\odot}^{\rm RG}$ values instead of
$\rho_{\rm d,\odot}$ in Eqs. (5)-(6) to calculate model values for the median
$\overline{v_{\phi}}$ and the velocity dispersions $\sigma_{R},\sigma_{\phi}$,
and $\sigma_{z}$. Thereafter, Monte Carlo random sampling is used to calculate
the model values for each grid of the data with a particular combination of
the 10 fit parameters, i.e., $\sigma_{i,\odot}^{\rm thin}$,
$\sigma_{i,\odot}^{\rm thick}$, $\beta_{i}$, $R_{\sigma_{i}}^{\rm thin}$, and
$R_{\sigma_{i}}^{\rm thick}$ ($i=R,z$). We assume that data incompleteness
does not affect the kinematic statistics in each grid and do not consider the
completeness correction for the comparison between the data and model values.

Figure 1: Far left: The Gaia DR2 giant sample’s $R$–$z$ distributions for
median azimuthal velocity ($\overline{v_{\phi}}$) and velocity dispersions
along the radial ($\sigma_{R}$), azimuthal ($\sigma_{R}$), and vertical
($\sigma_{z}$) directions from top to bottom (reproduction of Figure 11 of
Gaia Collaboration et al., 2018). The right four columns: The same
distributions from each indicated model. Only the grids inside the purple
boxes are used in the calculation of $\tilde{\chi}^{2}$ for each model. Figure
2: Residuals (model $-$ data) corresponding to the four models plotted in Fig.
1.

To evaluate an agreement with the data at $(i_{R},i_{z})$th grid, we use

$\displaystyle\tilde{\chi}^{2}_{i_{R},i_{z}}=\sum_{p=\overline{v_{\phi}},\sigma_{R},\sigma_{\phi},\sigma_{z}}w_{i_{R},i_{z}}(p_{i_{R},i_{z}}^{\rm
simu}-p_{i_{R},i_{z}}^{\rm obs})^{2},$ (20)

where

$\displaystyle w_{i_{R},i_{z}}=\frac{\left(1/N_{i_{R},i_{z}}^{\rm
obs}+1/N_{i_{R},i_{z}}^{\rm
simu}\right)^{-1}}{\sum_{j_{R},j_{z}}\left(1/N_{j_{R},j_{z}}^{\rm
obs}+1/N_{j_{R},j_{z}}^{\rm simu}\right)^{-1}},$ (21)

and $N_{i_{R},i_{z}}^{\rm obs}$ is the number of observed stars in a grid,
while $N_{i_{R},i_{z}}^{\rm simu}$ is the number of simulated stars in the
same grid. $N_{i_{R},i_{z}}^{\rm obs}$ takes 30–34243 depending on grid
position, where the median is $\sim 250$ and 67 out of 1207 grids have
$N_{i_{R},i_{z}}^{\rm obs}>10^{4}$. The weight for each grid,
$w_{i_{R},i_{z}}$, disregards each measurement uncertainty because we were not
provided the uncertainties of measurements of the median velocity and velocity
dispersion through the private communication. This corresponds to an
assumption that the average measurement uncertainty in each grid is the same
regardless of the grid; however, this is not true because velocity
measurements in a grid further from the Sun tend to have larger uncertainties
than those in a closer grid due to its relative faintness, and hence a closer
grid should have a larger weight considering this effect. Nevertheless, the
relative weights among different grids are, qualitatively, correctly set by
the current form of $w_{i_{R},i_{z}}$ because a closer grid tends to have a
larger $N_{i_{R},i_{z}}^{\rm obs}$ value, which makes its relative weight
larger. Therefore, although there should be an underestimation of relative
weight for a grid closer to the Sun, our result is not expected to change
significantly with the missed effect.

Another concern with the $w_{i_{R},i_{z}}$ expression by Eq. (21) is its
dependency on $N_{i_{R},i_{z}}^{\rm simu}$. To reduce the dependency on the
simulated number of stars but not increase the computation time significantly,
we adopt $N_{i_{R},i_{z}}^{\rm simu}={\rm min}(10^{4},4\,N_{i_{R},i_{z}}^{\rm
obs})$ for each grid. The maximum simulation number of $10^{4}$ indicates that
the $N_{i_{R},i_{z}}^{\rm simu}$ term becomes dominant in $w_{i_{R},i_{z}}$
when $N_{i_{R},i_{z}}^{\rm obs}>10^{4}$. We set this maximum primarily to
reduce the computation time, but also to avoid placing too much weight for
grids with large $N_{i_{R},i_{z}}^{\rm obs}$ values, because we want to have a
model that matches a wide range of grid points instead of being optimized to a
small fraction of the grid points that is highly weighted.

A particular set of the 10 fit parameters is evaluated by
$\tilde{\chi}^{2}=\sum_{i_{R},i_{z}}\tilde{\chi}^{2}_{i_{R},i_{z}}$. The
summation for $\tilde{\chi}^{2}$ only run over 651 grids with $R\leq
8440~{}{\rm pc}$ out of a total of 1207 grids because our primary science
interest is in microlensing events toward the Galactic bulge and because the
outer disk has a warp and/or a flare feature, which is not considered in our
model. The cut at $8440~{}{\rm pc}$ is motivated by the starting
Galactocentric radius of the warp used in the Besançon model (Robin et al.,
2003) of $R=8400~{}{\rm pc}$. We refer to a model fitted with this data
selection as an all-$z$ model. Additionally, we consider another option, which
is highly optimized for microlensing studies toward the Galactic bulge, where
the summation for $\tilde{\chi}^{2}$ run over 234 grids with $R\leq
8440~{}{\rm pc}$ and $|z|\leq 900~{}{\rm pc}$. With this option, we also
ignore the agreement in $\sigma_{R}$, i.e., the
$\tilde{\chi}^{2}_{i_{R},i_{z}}$ in Eq. (20) just run over
$p=\overline{v_{\phi}},\sigma_{\phi},\sigma_{z}$ and not over $p=\sigma_{R}$.
This is because $\sigma_{R}$ is irrelevant to observables of microlensing
events in fields near the Galactic center. We refer to a model from this data
selection as a low-$z$ model.

###  3.3 Grid Search

Because the Monte Carlo simulation to calculate $\tilde{\chi}^{2}$ is
computationally expensive, we conduct a grid search to find a combination of
parameters with a better agreement to the data. We conduct a sparse initial
grid search of the 10 fit parameters with a large interval between adjacent
grids at first, and then repeat it in narrower parameter space with a smaller
interval. The initial grid search run over the following ranges;
$\sigma_{R,\odot}^{\rm thin}=30$–$50~{}{\rm km/s}$, $\sigma_{z,\odot}^{\rm
thin}=15$–$30~{}{\rm km/s}$, $\sigma_{R,\odot}^{\rm thick}=40$–$70~{}{\rm
km/s}$, $\sigma_{z,\odot}^{\rm thick}=35$–$65~{}{\rm km/s}$,
$\beta_{R}=0.1$–0.4, $\beta_{z}=0.2$–0.8, $R_{\sigma_{R}}^{\rm
thin}=8$–$30~{}{\rm kpc}$, $R_{\sigma_{z}}^{\rm thin}=8$–$30~{}{\rm kpc}$,
$R_{\sigma_{R}}^{\rm thick}=8$–$30~{}{\rm kpc}$, and $R_{\sigma_{z}}^{\rm
thick}=8$–$30~{}{\rm kpc}$. Once a best-fit combination was found at an edge
of the parameter space searched, we expanded the parameter space in the next
iteration. This procedure was repeated until no significant $\tilde{\chi}^{2}$
improvement was found with a smaller interval or an expansion of the parameter
space.

We conduct the above search for each of the four options, i.e., all
combinations of the two options for the data selection (all-$z$ or low-$z$)
and thin disk scale height (flat or linear), respectively. Table 2 shows the
best-fit parameters given by the calculations. Figs. 1 and 2 show color maps
of each model values and residuals (model value $-$ data value), respectively.
In each panel of these figures, we indicate the selected grid region with
magenta dashed boxes. Fig. 2 shows that our models fail to reproduce the outer
disk distributions, in particular, the $\overline{v_{\phi}}$, but this is
expected because our density model is designed for the inner disk and does not
consider the warp or flare structure seen in the outer disk. Similarly, the
$\sigma_{z}$ distribution of the low-$z$ model significantly overestimates the
values in $z>900$ pc which is outside the magenta box. This is because the
best-fit $\sigma_{z,\odot}^{\rm thick}$ value for a low-$z$ model (61.4 or
59.0 km/s) is much higher than that for an all-$z$ model (49.2 or 47.8 km/s),
as shown in Table 2. However, this difference has little effect in the
microlensing region toward the Galactic bulge because the thick disk stars are
relatively rare in this region. Focusing on the region related to the
microlensing study (i.e., inside the magenta box of a low-$z$ model), all the
four models show moderate agreements with the Gaia data.

The $\tilde{\chi}^{2}$ value is smallest for the low-$z$ $+$ linear $z_{\rm
d}^{\rm thin}$ model; however, comparison of the $\tilde{\chi}^{2}$ values
between low-$z$ and all-$z$ models is unreasonable because the grids
contributing to the $\tilde{\chi}^{2}$ are different. The $\tilde{\chi}^{2}$
values, defined by Eq (20), is the weighted root mean square of deviation of
the model from the data, that is, $\sqrt{\tilde{\chi}^{2}}$ provides a
weighted average of the deviations in km/s. Notably, the linear scale height
models are preferred for both the all-$z$ and low-$z$ models. Although this
might indicate that the scale height is not constant inside the solar radius,
the flat scale height models are favored in the fits to the data toward bulge
regions as described in Section 4.3.2.

We keep all the four models as options for the disk model, and use each of
them combined with a bulge model in the fits conducted in Section 4, and the
four models are compared in Section 4.3.2 with respect to the best-fit
$\chi^{2}$ values to the bulge data. We discuss the determined values in Table
2 by comparing with the previous studies in Section 5.1.

Table 2: Best-fit parameters for the disk kinematic model.  Model | $N_{\rm grid}$ | $\sigma_{R,\odot}^{\rm thin}$ | $\sigma_{z,\odot}^{\rm thin}$ | $\sigma_{R,\odot}^{\rm thick}$ | $\sigma_{z,\odot}^{\rm thick}$ | $\beta_{R}$ | $\beta_{z}$ | $R_{\sigma_{R}}^{\rm thin}$ | $R_{\sigma_{z}}^{\rm thin}$ | $R_{\sigma_{R}}^{\rm thick}$ | $R_{\sigma_{z}}^{\rm thick}$ | $\tilde{\chi}^{2}_{\overline{v_{\phi}}}$ | $\tilde{\chi}^{2}_{\sigma_{R}}$ | $\tilde{\chi}^{2}_{\sigma_{\phi}}$ | $\tilde{\chi}^{2}_{\sigma_{z}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  | [km/s] | [km/s] | [km/s] | [km/s] |  |  | [kpc] | [kpc] | [kpc] | [kpc] |  |  |  |
all-$z$ $+$ flat $z_{\rm d}^{\rm thin}$ | 651 | 42.0 | 24.4 | 75 | 49.2 | 0.32 | 0.77 | 14.3 | 5.9 | 180 | 9.4 | 19.5 | 16.3 | 15.9 | 9.6
all-$z$ $+$ linear $z_{\rm d}^{\rm thin}$ | 651 | 44.0 | 25.4 | 68 | 47.8 | 0.34 | 0.81 | 21.4 | 8.1 | 57.6 | 15.6 | 14.9 | 13.7 | 14.1 | 6.7
low-$z$ $+$ flat $z_{\rm d}^{\rm thin}$ | 234 | 35.2 | 22.2 | 75 | 61.4 | 0.22 | 0.77 | 9.5 | 10.4 | 90.0 | 6.9 | 12.6 | 7.2 | – | 4.2
low-$z$ $+$ linear $z_{\rm d}^{\rm thin}$ | 234 | 37.6 | 23.4 | 68 | 59.0 | 0.30 | 0.82 | 11.1 | 7.8 | 47.0 | 52.0 | 8.7 | 5.8 | – | 3.5

##  4 Fitting for the Bulge and IMF parameters

In this section, we determine the 4 fit parameters for the IMF model, 7–15 fit
parameters for the $\rho_{\rm B}$ model, and 19 fit parameters for the $v_{\rm
B}$ model, through the Markov Chain Monte Carlo (MCMC) fitting to the observed
distributions toward the bulge sky of the OGLE-III RC star count (Nataf et
al., 2013), VIRAC proper motion (Smith et al., 2018; Clarke et al., 2019),
BRAVA radial velocity (Rich et al., 2007; Kunder et al., 2012), and star and
microlensing event count by OGLE-IV (Mróz et al., 2017, 2019).

We use the bulge model combined with the disk model in the fit, where we fix
and use the 10 fit parameters for the $v_{\rm d}$ model determined in Section
3. However, the local density values are recalculated in the fit for every
given set of the 4 fit parameters in the IMF.

###  4.1 Data and Corresponding $\chi^{2}$ Values

This section describes details of each dataset used to constrain the 30–38 fit
parameters. For each dataset, we also describe model values compared with the
observed values and introduce the corresponding $\chi^{2}$ values. All the
used data are plotted in Figs. 4, 6, and 8.

####  4.1.1 OGLE-III red clump star count

Nataf et al. (2013) divided the OGLE-III Galactic bulge fields over 90.25 deg2
in $-10^{\circ}<l<10^{\circ}$ and $2^{\circ}<|b|<7^{\circ}$ into 9019 areas
and measured the red clump (RC) star count in each area. They used a
luminosity function model that primarily consists of a Gaussian RC component
on an exponential red giant branch continuum for fit to the observed $I$-mag
distribution in each area. As a result of the fits, they provided a catalog of
9019 sets of the number count $N_{\rm RC}$, mean distance modulus ${\rm
DM}_{\rm RC}$, variance of distance modulus $\sigma_{\rm DM}^{2}$, and error
matrix for them, as well as the extinctions $A_{I}$ and $A_{V}$, where $N_{\rm
RC}$ and ${\rm DM}_{\rm RC}$ are equivalent to the area and the peak of the
Gaussian, respectively. The variance of distance modulus, $\sigma_{\rm
DM}^{2}$, was determined by subtracting the sum of variances of the extinction
in the area ($=\sigma_{A_{I}}^{2}$) and of the intrinsic brightness of RC
($=\sigma_{I,{\rm RC},0}^{2}$) from the variance of the Gaussian
($=\sigma_{I,{\rm RC}}^{2}$).

Nataf et al. (2013) used $\sigma_{I,{\rm RC},0}=0.09$ with no uncertainty to
derive the $\sigma_{\rm DM}^{2}$ values. To be conservative, we apply
$\sigma_{I,{\rm RC},0}=0.15\pm 0.06$, i.e., we modify the data by $\sigma_{\rm
DM}^{2}=\sigma_{\rm DM,org}^{2}+0.09^{2}-0.15^{2}$ and increase the
uncertainty accordingly. This choice for $\sigma_{I,{\rm RC},0}$ is motivated
by Hawkins et al. (2017) who measured the RC mag dispersions of $0.20\pm 0.02$
in both $G$\- and $J$-bands. Because the wavelength of $I$-band is between
these bands, we conservatively take $0.15\pm 0.06$ so that its 1-$\sigma$
range includes both 0.09 used by Nataf et al. (2013) and 0.20 measured for
$G$\- and $J$-bands.

Moreover, because the measurements of $N_{\rm RC},{\rm DM}_{\rm RC}$, and
$\sigma_{\rm DM}^{2}$ are all from the Gaussian fit, it could
overestimate/underestimate the RC population when it is not distributed
following a Gaussian shape. In particular, because our disk model has a
constant surface density at $R<5.3$ kpc (see Eqs. 5–6), which is continuously
distributed in the bulge region mildly, most or part of the disk RC population
is expected to be absorbed into the exponential red giant branch continuum
component in the fit by Nataf et al. (2013), although the fraction of
absorption probably depends on the line of sight. To account for this
uncertainty, we considered only the bulge population in the model observables
given below in Eqs. (22)–(24), in addition to increasing the uncertainties of
$N_{\rm RC}$ by 10% and adding 0.04 mag error to the uncertainties of ${\rm
DM}_{\rm RC}$ in quadrature.

Cao et al. (2013) modeled a bulge density distribution by fitting to the Nataf
et al. (2013) data, and we similarly follow their parameterization but with a
modification on the integration range. For a particular $i$th line of sight
toward $(l_{i},b_{i})$, the model RC number count is expressed as

$\displaystyle N_{{\rm RC},i}^{\rm mod}=\Omega_{i}N_{{\rm
B},i}\times\frac{\sum_{j}N_{{\rm RC},j}^{\rm obs}}{\sum_{j}\Omega_{j}N_{{\rm
B},j}},$ (22)

where $\Omega_{i}$ is the sky area of the $i$th field, and $N_{{\rm B},i}$ is
the model number of bulge stars integrated along $i$th line of sight, which is
defined as

$\displaystyle N_{{\rm B},i}\equiv\int_{s_{{\rm min},i}}^{s_{{\rm
max},i}}n_{\rm B}(l_{i},b_{i},s)s^{2}ds,$

and $n_{\rm B}(l_{i},b_{i},s)$ is the number density of bulge stars at the
distance $s$ from the Sun toward $(l_{i},b_{i})$. The summation for $j$ in the
second factor in Eq. (22) runs over all the 9019 lines of sight, and the
factor is for a normalization to let the total $N_{{\rm RC},i}^{\rm mod}$ be
same as the observed one, i.e., to make $\sum_{i}N_{{\rm RC},i}^{\rm
mod}=\sum_{i}N_{{\rm RC},i}^{\rm obs}$.

The mean and variance of distance modulus are expressed as

$\displaystyle{\rm DM}_{{\rm RC},i}^{\rm mod}$ $\displaystyle=\frac{1}{N_{{\rm
B},i}}\int_{s_{{\rm min},i}}^{s_{{\rm max},i}}n_{\rm
B}(l_{i},b_{i},s)s^{2}{\rm DM}(s)ds$ (23)

and

$\displaystyle(\sigma_{{\rm DM},i}^{2})^{\rm mod}$
$\displaystyle=\frac{1}{N_{{\rm B},i}}\int_{s_{{\rm min},i}}^{s_{{\rm
max},i}}n_{\rm B}(l_{i},b_{i},s)s^{2}[{\rm DM}(s)]^{2}ds-({\rm DM}_{{\rm
RC},i}^{\rm mod})^{2},$ (24)

respectively, where ${\rm DM}(s)$ denotes the distance modulus at distance
$s$, which is given by ${\rm DM}(s)=5\log_{10}[s/(10{\rm pc})]$.

For the integration range, we use a pair of $s_{{\rm min},i}$ and $s_{{\rm
max},i}$ that satisfies

$\displaystyle{\rm DM}(s_{{\rm min},i})$ $\displaystyle={\rm DM}_{{\rm
RC},i}^{\rm obs}-{\rm min}(3\sigma_{I,{\rm RC},i}^{\rm obs}\,,1.5)$ (25)
$\displaystyle{\rm DM}(s_{{\rm max},i})$ $\displaystyle={\rm DM}_{{\rm
RC},i}^{\rm obs}+{\rm min}(3\sigma_{I,{\rm RC},i}^{\rm obs}\,,1.5),$ (26)

where $\sigma_{I,{\rm RC},i}^{\rm obs}$ is the standard deviation of the
Gaussian fit for the RC component in the $I$-band luminosity function of $i$th
field. We chose these values rather than $s_{{\rm min},i}=3~{}{\rm kpc}$ and
$s_{{\rm max},i}=13~{}{\rm kpc}$ used by Cao et al. (2013) because the RC
stars outside of $\pm 3$-$\sigma$ from the mean have a negligible contribution
to the measurements of the three observables, which are equivalent to the
area, peak position, and variance of the Gaussian distribution. The value of
1.5 for the maximum range of integration originates from the limit on the
fitting range of magnitude, $-1.5<I-I_{\rm RC}<1.5$, set by Nataf et al.
(2013).

Following Cao et al. (2013), we calculate the $\chi^{2}$ for this dataset by

$\displaystyle\chi^{2}_{\rm RC}=\sum_{i}({\bm{X}}_{i}^{\rm
mod}-{\bm{X}}_{i}^{\rm obs})^{T}{\bm{S}}^{-1}_{i}({\bm{X}}_{i}^{\rm
mod}-{\bm{X}}_{i}^{\rm obs}),$ (27)

where ${\bm{X}}_{i}=(N_{{\rm RC},i},{\rm DM}_{{\rm RC},i},\sigma_{{\rm
DM},i}^{2})$ and ${\bm{S}}^{-1}_{i}$ is the covariance matrix of the
uncertainties. We further define the following four $\chi^{2}$ values to
quantify a contribution to $\chi^{2}_{\rm RC}$ from each observable;

$\chi^{2}_{p}=\sum_{i}\left(\frac{p_{i}^{\rm mod}-p_{i}^{\rm obs}}{p_{i}^{\rm
err}}\right)^{2}\ \ \ \ (p=N_{\rm RC},{\rm DM}_{\rm RC},\sigma_{{\rm
DM}}^{2}),$

and $\chi^{2}_{\rm cov}=\chi^{2}_{\rm RC}-\chi^{2}_{N_{\rm RC}}-\chi^{2}_{{\rm
DM}_{\rm RC}}-\chi^{2}_{\sigma_{{\rm DM}}^{2}}$, where $p_{i}^{\rm err}$
denotes the error-bar of the $i$th data of parameter $p$.

####  4.1.2 VIRAC red giants’ proper motions

Smith et al. (2018) provided the VVV infrared astrometric catalogue (VIRAC),
which is a near-infrared proper motion catalog of the five years VISTA
Variables in the Via Lactea (VVV) survey (Minniti et al., 2010), which
includes 312,587,642 sources over 560 deg2 of the bulge and southern disk.
Clarke et al. (2019) calibrated the VIRAC proper motions by comparing the
VIRAC values to the corresponding Gaia values. They carefully selected red
giants with $11.8<K_{s,0}<13.6$ from the calibrated catalog, where $K_{s,0}$
indicates the extinction-corrected $K_{s}$ magnitude. Then they split each of
the VVV tiles located in the bulge region into 4 sub-tiles and presented means
$\left\langle\mu_{i}\right\rangle$ and dispersions $\sigma_{\mu_{i}}$
($i=l,b$) of proper motions of the selected giant stars in each sub-tile. We
use the data for $\left\langle\mu_{l}\right\rangle$, $\sigma_{\mu_{l}}$, and
$\sigma_{\mu_{b}}$ in 676 sub-tiles distributed roughly over
$-9^{\circ}<l<9.4^{\circ}$ and $-9.4^{\circ}<b<4.2^{\circ}$. Out of the 676
sub-tiles, we do not use 90 sub-tiles with $|b|\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1^{\circ}$ because of high
extinction values. Thus, we use the data in 586 sub-tiles for the fit. We do
not use the $\left\langle\mu_{b}\right\rangle$ data because they are little
sensitive to the fit parameters and not very useful to constrain them.

Figure 3: Comparison of proper motion dispersions (left $\sigma_{\mu_{l}}$,
right $\sigma_{\mu_{b}}$) between VVV and HST. Differences of measured values
are plotted as a function of extinction for 35 fields from Kozłowski et al.
(2006) and 4 fields from Table 2 of Terry et al. (2020). 22 fields with $\leq
40$ secs exposures from Kozłowski et al. (2006) are classified as the short
exp. fields while the other 17 fields are classified as the long exp. fields.
An outlier pointed by Terry et al. (2020) is indicated.

Seeing-limited observations toward the Galactic bulge often suffer from
systematic errors due to blended stars. Because the data are derived from
bright stars with $11.8<K_{s,0}<13.6$, they are not expected to be affected by
such systematics compared to fainter stars. To verify this, we compared the
proper motion dispersion measurements from VIRAC with those from the Hubble
Space Telescope (HST) observations. Fig. 3 shows the comparison as a function
of extinction $A_{I}$. We used 35 fields observed by Kozłowski et al. (2006)
and 4 fields summarized in Table 2 of Terry et al. (2020) for the comparison,
and the HST values are taken from the two papers. The 4 HST values from Terry
et al. (2020) were originally measured by Calamida et al. (2014), Kuijken &
Rich (2002), and Terry et al. (2020). If there are no systematic errors, the
mean value of $\sigma_{\mu_{i},{\rm VVV}}-\sigma_{\mu_{i},{\rm HST}}$ should
be consistent with 0. However, we found the mean and standard deviation values
of $\left\langle\sigma_{\mu_{l},{\rm VVV}}-\sigma_{\mu_{l},{\rm
HST}}\right\rangle=0.124\pm 0.135$ mas/yr and
$\left\langle\sigma_{\mu_{b},{\rm VVV}}-\sigma_{\mu_{b},{\rm
HST}}\right\rangle=0.153\pm 0.144$ mas/yr when we used the 39 fields without
the outlier indicated in the figure that is specified by Terry et al. (2020).
This implies the existence of systematic offset between the VVV and HST
measurements. We, at first, checked whether there was any correlation between
the offset and the extinction, but Fig. 3 shows no such a correlation except
for the outlier.

Thereafter, we divide the sample of 39 fields into two subsamples depending on
the exposure time of the HST observations because certain fields have
significantly shorter exposure times than others. One of the two subsamples
consists of 22 short exposure fields from Kozłowski et al. (2006) with each
exposure less than 40 s, and the other consists of 17 long exposure fields
with each exposure longer than 100 s. Note that Kozłowski et al. (2006) used
1st epoch data taken by the WFPC2/PC camera, which is much less sensitive than
the ACS camera used to take their 2nd epoch data. The mean and standard
deviation values are $\left\langle\sigma_{\mu_{l},{\rm
VVV}}-\sigma_{\mu_{l},{\rm HST}}\right\rangle=0.205\pm 0.116$ mas/yr and
$\left\langle\sigma_{\mu_{b},{\rm VVV}}-\sigma_{\mu_{b},{\rm
HST}}\right\rangle=0.246\pm 0.108$ mas/yr for the 22 short exposure fields and
$\left\langle\sigma_{\mu_{l},{\rm VVV}}-\sigma_{\mu_{l},{\rm
HST}}\right\rangle=0.023\pm 0.096$ mas/yr and
$\left\langle\sigma_{\mu_{b},{\rm VVV}}-\sigma_{\mu_{b},{\rm
HST}}\right\rangle=0.034\pm 0.093$ mas/yr for the 17 long exposure fields
without the outlier. The long exposure subsample shows no clear difference
between VVV and HST measurements while the short exposure subsample does.

Thus, we conclude that the implied offset is due to the systematic error in
the HST measurements for the 22 short exposure fields. Given the standard
deviation of $\sim 0.1$ mas/yr for the long exposure subsample, we use 0.1
mas/yr for the uncertainty of $\sigma_{\mu_{l}}$ and $\sigma_{\mu_{b}}$ of
each sub-tile because formal statistic errors are comparatively small. For the
uncertainty of $\left\langle\mu_{l}\right\rangle$, we use 0.14 mas/yr by
additionally considering the error of calibration to the Gaia scale of $\sim
0.1$ mas/yr (Clarke et al., 2019).

For a particular $i$th subtile’s coordinate of $(l_{i},b_{i})$, the model
$\left\langle\mu_{l}\right\rangle$ value is calculated by

$\displaystyle\left\langle\mu_{l}\right\rangle_{i}^{\rm mod}$
$\displaystyle=\frac{\int_{3{\rm kpc}}^{16{\rm kpc}}\left[n_{\rm d}^{\rm
RG}(l_{i},b_{i},s)\left\langle\mu_{l,{\rm
d}}(l_{i},b_{i},s)\right\rangle+n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)\left\langle\mu_{l,{\rm
B}}(l_{i},b_{i},s)\right\rangle\right]w(s)ds}{\int_{3{\rm kpc}}^{16{\rm
kpc}}\left[n_{\rm d}^{\rm RG}(l_{i},b_{i},s)+n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)\right]w(s)ds},$ (28)

where $n_{\rm d}^{\rm RG}(l_{i},b_{i},s)$ and $n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)$ are the number densities of red giants of the disk and
bulge components at $(l_{i},b_{i},s)$, respectively. We used the same
definition for red giants in Section 2.2.1, i.e., stars with ${\cal
M}_{G}<3.9$ mag and $(G_{\rm BP}-G_{\rm RP})_{0}>0.95$.
$\left\langle\mu_{l,{\rm d}}(l_{i},b_{i},s)\right\rangle$ and
$\left\langle\mu_{l,{\rm B}}(l_{i},b_{i},s)\right\rangle$ are mean proper
motions calculated using our disk and bulge velocity model at
$(l_{i},b_{i},s)$, respectively. Calculations of proper motion require the
solar velocity, and again we use
$(v_{\odot,x},v_{\odot,y},v_{\odot,z})=(-10,243,7)$ km/s in this study. The
weight $w(s)$ is given by

$\displaystyle w(s)=s^{2}\int_{11.8}^{13.6}L_{{\cal
M}_{K_{s}}}\left(K_{s,0}-{\rm DM}(s)\right)dK_{s,0},$ (29)

where $L_{{\cal M}_{K_{s}}}({\cal M}_{K_{s}})$ is luminosity function for the
red giant’s absolute magnitude ${\cal M}_{K_{s}}$, given by Eqs. (2)–(5) of
Clarke et al. (2019). Further, the model $\sigma_{\mu_{j}}$ $(j=l,b)$ values
are expressed as

$\displaystyle\sigma_{\mu_{j},i}^{\rm mod}$
$\displaystyle=\sqrt{\left\langle\mu_{j}^{2}\right\rangle_{i}^{\rm
mod}-\left(\left\langle\mu_{j}\right\rangle_{i}^{\rm mod}\right)^{2}},$ (30)

where

$\left\langle\mu_{j}^{2}\right\rangle_{i}^{\rm mod}=\frac{\int_{3{\rm
kpc}}^{16{\rm kpc}}\left[n_{\rm d}^{\rm
RG}(l_{i},b_{i},s)\left\langle\mu_{j,{\rm
d}}^{2}(l_{i},b_{i},s)\right\rangle+n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)\left\langle\mu_{j,{\rm
B}}^{2}(l_{i},b_{i},s)\right\rangle\right]w(s)ds}{\int_{3{\rm kpc}}^{16{\rm
kpc}}\left[n_{\rm d}^{\rm RG}(l_{i},b_{i},s)+n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)\right]w(s)ds}.$

$\chi^{2}$ for this dataset is denoted by

$\chi^{2}_{\mu}=\chi^{2}_{\mu_{l}}+\chi^{2}_{\sigma_{\mu_{l}}}+\chi^{2}_{\sigma_{\mu_{b}}},$

where

$\chi^{2}_{p}=\sum_{i}\left(\frac{p_{i}^{\rm mod}-p_{i}^{\rm obs}}{p_{i}^{\rm
err}}\right)^{2}\ \ \ \ (p=\mu_{l},\sigma_{\mu_{l}},\sigma_{\mu_{b}}).$

####  4.1.3 BRAVA radial velocity data

The Bulge Radial Velocity Assay (BRAVA, Rich et al., 2007; Kunder et al.,
2012) is a large spectroscopic survey of M giant stars in the Galactic bulge
to constrain the bulge dynamics by measuring their radial velocities (RVs). We
use the mean and dispersion values of the RV measurements in 82 fields where
80 fields are located at $-10^{\circ}\leq l\leq 10^{\circ}$ and
$3^{\circ}\leq|b|\leq 8^{\circ}$, while the other 2 are located at
$(l,b)=(0^{\circ},-1^{\circ})$ and $(l,b)=(0^{\circ},-2^{\circ})$. Although RV
is not very relevant to our interest in microlensing observables, we use these
data as a constraint because RV provides direct information of velocity
compared to proper motion, which is a combination of distance and velocity.

The BRAVA mean RV values are given in the Galactocentric frame by a conversion
from the originally observed values in the heliocentric frame. The conversion
was done using the Sun velocity of $(-9.0,231.9,7.0)$ km/s along $(x,y,z)$
axes (Howard et al., 2008), which is slightly different from our values of
(-10, 243, 7) km/s. Thus, we reconverted the mean RV values using our values
and then use them as data for our fits.

Similar to $\left\langle\mu_{l}\right\rangle_{i}^{\rm mod}$, the model value
for mean RV is given by

$\displaystyle\left\langle\rm RV\right\rangle_{i}^{\rm mod}$
$\displaystyle=\frac{\int_{3{\rm kpc}}^{16{\rm kpc}}\left[n_{\rm d}^{\rm
RG}(l_{i},b_{i},s)\left\langle{\rm RV}_{\rm
d}(l_{i},b_{i},s)\right\rangle+n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)\left\langle{\rm RV}_{\rm
B}(l_{i},b_{i},s)\right\rangle\right]s^{0.6}ds}{\int_{3{\rm kpc}}^{16{\rm
kpc}}\left[n_{\rm d}^{\rm RG}(l_{i},b_{i},s)+n_{\rm B}^{\rm
RG}(l_{i},b_{i},s)\right]s^{0.6}ds},$ (31)

where the weight $s^{0.6}$ follows Portail et al. (2015) who combined the
volume effect ($\propto s^{2}$) with an approximate luminosity function of
giant stars of $\propto 10^{0.28{\cal M}_{K}}$ (Wegg & Gerhard, 2013), which
resulted in a dependency on the distance of $\propto s^{-1.4}$. Further, the
model RV dispersion value is calculated using the following

$\displaystyle\sigma_{{\rm RV},i}^{\rm mod}$
$\displaystyle=\sqrt{\left\langle\rm RV^{2}\right\rangle_{i}^{\rm
mod}-\left(\left\langle\rm RV\right\rangle_{i}^{\rm mod}\right)^{2}}.$ (32)

$\chi^{2}$ for this dataset is denoted by

$\chi^{2}_{\rm BRA}=\chi^{2}_{\rm RV}+\chi^{2}_{\sigma_{\rm RV}},$

where

$\chi^{2}_{p}=\sum_{i}\left(\frac{p_{i}^{\rm mod}-p_{i}^{\rm obs}}{p_{i}^{\rm
err}}\right)^{2}\ \ \ \ (p=\left\langle\rm RV\right\rangle,\sigma_{\rm RV}).$

####  4.1.4 OGLE-IV star count data

Mróz et al. (2019) analyzed long-term photometric observations of the Galactic
bulge by the OGLE group using their fourth-generation wide-field camera OGLE-
IV. The analyzed fields consist of 9 high-cadence fields and 112 low-cadence
fields, with each field further divided into 32 subfields having an individual
area of $\sim 0.044~{}{\rm deg}^{2}$. They measured the number of stars with
$I<21~{}{\rm mag}$ and $I<18~{}{\rm mag}$, denoted by $N_{I<21}$ and
$N_{I<18}$ respectively, for each subfield in addition to the microlens-
related observables described in Section 4.1.5. Note that $N_{I<21}$ is used
as the number of candidate source stars in the microlensing events. Although
an extinction value $A_{I}$ is needed for each line of sight to model
$N_{I<I_{\rm c}}$ ($I_{\rm c}=18,21$), there is no extinction catalog in
$I$-band covering all the OGLE-IV fields. Thus, we only use 1456 subfields
covered by the $A_{I}$ map of Nataf et al. (2013) for our analysis.

The model $N_{I<I_{\rm c}}$ for a particular $i$th subfield is calculated from

$\displaystyle N_{I<I_{\rm c},i}^{\rm mod}$
$\displaystyle=\Omega_{i}\int_{3{\rm kpc}}^{16{\rm kpc}}n_{I<I_{\rm
c}}(l_{i},b_{i},s)s^{2}ds\hskip 19.91684pt(I_{\rm c}=18,21),$ (33)

where $\Omega_{i}$ is the sky area, and $n_{I<I_{\rm c}}(l_{i},b_{i},s)$ is
the number density of stars with $I<I_{\rm c}$ at a distance $s$ from the Sun
toward $(l_{i},b_{i})$, which is expressed as

$\displaystyle n_{I<I_{\rm c}}(l_{i},b_{i},s)=[n_{\rm d}(l_{i},b_{i},s)+n_{\rm
B}(l_{i},b_{i},s)]\int_{I<I_{\rm c}}L_{{\cal M}_{I}}\left(I-{\rm
DM}(s)-A_{I}(l_{i},b_{i},s)\right)dI$ (34)

with the stellar number density for disk and bulge stars, $n_{\rm
d}(l_{i},b_{i},s)$ and $n_{\rm B}(l_{i},b_{i},s)$, respectively, and the
luminosity function for $I$-band absolute magnitude $L_{{\cal M}_{I}}({\cal
M}_{I})$ calculated using the PARSEC isochrone models (Bressan et al., 2012;
Chen et al., 2014; Tang et al., 2014) for a given IMF. We use a 3D extinction
distribution of

$\displaystyle A_{I}(l_{i},b_{i},s)=A_{I,{\rm
RC}}(l_{i},b_{i})\frac{1-\exp[-s/(h_{\rm dust}/\sin|b_{i}|)]}{1-\exp[-s_{\rm
RC}(l_{i},b_{i})/(h_{\rm dust}/\sin|b_{i}|)]},$ (35)

where $A_{I,{\rm RC}}(l_{i},b_{i})$ and $s_{\rm RC}(l_{i},b_{i})$ are the mean
extinction for the RC and mean distance to RC, respectively, both taken from
the Nataf et al. (2013) extinction catalog, and we use $h_{\rm
dust}=164~{}{\rm pc}$ based on the measurement by Nataf et al. (2013). Because
the extinction distribution is not uniform even inside a subfield of $\sim
0.044~{}{\rm deg}^{2}$, we further divide each subfield into 8 sub-subfields
and then calculate the $N_{I<I_{\rm c},i}^{\rm mod}$ value for each subfield
by summing the values over the 8 sub-subfields with different $A_{I,{\rm RC}}$
values.

$\chi^{2}$ for the star count data is given by

$\chi^{2}_{N_{I<I_{\rm c}}}=\sum_{i}\left(\frac{N_{I<I_{\rm c},i}^{\rm
mod}-N_{I<I_{\rm c},i}^{\rm obs}}{N_{I<I_{\rm c},i}^{\rm
err}}\right)^{2}\hskip 19.91684pt(I_{\rm c}=18,21),$

where we apply $N_{I<21,i}^{\rm err}=0.12\,N_{I<21,i}^{\rm obs}$ because Mróz
et al. (2019) estimated the errors as 10 to 15 % by comparing to star count in
HST images. We use a larger uncertainty of $N_{I<18,i}^{\rm
err}=0.14\,N_{I<18,i}^{\rm obs}$ because evolved stars like red giants mainly
contribute to $N_{I<18}$ and the isochrone model uncertainty for such evolved
stars is likely larger than that for main-sequence stars contributing to
$N_{I<21}$.

The model value $N_{I<I_{\rm c}}^{\rm mod}$ depends on the density model, as
well as the IMF, because the luminosity function $L_{{\cal M}_{I}}({\cal
M}_{I})$ depends on the IMF. As a quantity that is less dependent on the
density model but more on the IMF, we defined $f_{18/21}\equiv
N_{I<18}/N_{I<21}$. $\chi^{2}$ for this quantity is calculated by

$\chi^{2}_{f_{18/21}}=\sum_{i}\left(\frac{f_{18/21,i}^{\rm
mod}-f_{18/21,i}^{\rm obs}}{f_{18/21,i}^{\rm err}}\right)^{2},$

where we used $f_{18/21,i}^{\rm err}=0.16\,f_{18/21,i}^{\rm obs}$. This
uncertainty is taken slightly smaller than just a square root of the sum of
the $N_{I<21,i}^{\rm err}$ and $N_{I<18,i}^{\rm err}$ considering a positive
correlation between $N_{I<21,i}^{\rm obs}$ and $N_{I<18,i}^{\rm obs}$.

####  4.1.5 Number of microlensing events and $t_{\rm E}$ distribution by
OGLE-IV

Figure 4: Far left: Data distribution of six quantities used to constrain the
bulge density ($\rho_{\rm B}$) model (Nataf et al., 2013; Mróz et al., 2019).
The right six rows: The same distributions from each indicated model.
$\chi^{2}$ values defined for each quantity are shown. In the bottom panels
indicated by $N_{\rm eve}$, total number of (expected) event detections in
$i$th subfield, $\sum_{j=1}^{20}N_{{\rm eve},i}(t_{{\rm E},j})$, are plotted.
Figure 5: Residuals corresponding to Fig. 4. In $i$th subfield, $(p_{i}^{\rm
mod}-p_{i}^{\rm obs})/p_{i}^{\rm err}$, is plotted where $p$ is each indicated
parameter on the far right. Figure 6: Far left: Data distribution of five
quantities used to constrain the bulge velocity ($v_{\rm B}$) model (Clarke et
al., 2019; Rich et al., 2007; Kunder et al., 2012). The right six rows: The
same distributions from each indicated model. $\chi^{2}$ values defined for
each quantity are shown. Figure 7: Residuals corresponding to Fig. 6. In $i$th
subfield, $(p_{i}^{\rm mod}-p_{i}^{\rm obs})/p_{i}^{\rm err}$, is plotted
where $p$ is each indicated parameter. Figure 8: Left: The original OGLE-IV
$t_{\rm E}$ distributions of 2212 events in the high-cadence fields by Mróz et
al. (2017) (top) and 5788 events in the low-cadence fields by Mróz et al.
(2019) (bottom). Right: Modified $t_{\rm E}$ distributions by applying the
Golovich et al. (2020) correction with the Gaussian Process model (G20 option,
considered in Appendix A). The black open circles and the red curves are
indicated data and models, respectively.

Mróz et al. (2019) also provided the microlensing optical depth and the event
rate for each OGLE-IV field. However, it is statistically easier to deal with
another equivalent quantity, $N_{{\rm eve},i}(t_{\rm E})$, which is the number
of detected events as a function of $t_{\rm E}$ for each $i$th subfield,
because $N_{{\rm eve},i}(t_{\rm E})$ simply follows a Poisson distribution.
For a particular $i$th subfield and a $j$th $t_{\rm E}$ bin, the expected
number of event detections is given by

$\displaystyle N_{{\rm eve},i}^{\rm mod}(t_{{\rm E},j})$
$\displaystyle=\frac{2N_{I<21,i}^{\rm mod}T_{\rm o}}{\pi\left\langle t_{\rm
E}\right\rangle_{i}}\tau_{I<21,i}\,f_{i}(t_{{\rm E},j})\epsilon_{i}(t_{{\rm
E},j})\Delta t_{\rm E},$ (36)

where $T_{\rm o}$ is the survey duration that takes 2741 or 2011 days
depending on the field, $\tau_{I<21,i}$ is the average optical depth over the
monitored $N_{I<21,i}^{\rm mod}$ stars, which is given by

$\displaystyle\tau_{I<21,i}=\frac{1}{N_{I<21,i}^{\rm mod}}\int_{3{\rm
kpc}}^{16{\rm kpc}}\tau(l_{i},b_{i},s)\,n_{I<21}(l_{i},b_{i},s)s^{2}ds$ (37)

with

$\displaystyle\tau(l_{i},b_{i},s)=\frac{4\pi
G}{c^{2}}\int_{0}^{s}\rho(l_{i},b_{i},s_{\rm L})\frac{s_{\rm L}(s-s_{\rm
L})}{s}ds_{\rm L},$ (38)

$f_{i}(t_{\rm E})$ is the probability density function of $t_{\rm E}$,
$\left\langle t_{\rm E}\right\rangle_{i}$ is the mean $t_{\rm E}$ value given
by $\left\langle t_{\rm E}\right\rangle_{i}=\int f_{i}(t_{\rm E})\,t_{\rm
E}\,dt_{\rm E}$, $\epsilon_{i}(t_{\rm E})$ is the detection efficiency for an
event with $t_{\rm E}$ in the subfield, and $\Delta t_{\rm E}$ is the $t_{\rm
E}$ bin size. We have been provided with the detection efficiency values in
addition to the number of detected events as a function of $t_{\rm E}$ for
each subfield through a private communication with P. Mróz.

We determine $\chi^{2}$ for the number of microlensing event detections with

$\displaystyle\chi^{2}_{N_{\rm eve}}=\sum_{i,j}-2\ln P[N_{{\rm eve},i}^{\rm
obs}(t_{{\rm E},j});N_{{\rm eve},i}^{\rm mod}(t_{{\rm E},j})],$ (39)

where $P(x;\lambda)=\lambda^{x}e^{-\lambda}/x!$ is the Poisson probability of
$x$ with the mean value $\lambda$. A $j$th $t_{\rm E}$ bin includes events
with $0.125(j-1)\leq\log[t_{\rm E}/{\rm days}]<0.125j$ ($j=$ 1, …, 20) for
$i\in{\rm O_{\rm low}}$, while it includes events with
$-0.30+0.14(j-1)\leq\log[t_{\rm E}/{\rm days}]<-0.30+0.14j$ ($j=$ 1, …, 20)
for $i\in{\rm O_{\rm hi}}$, where $i\in{\rm O_{\rm low}}$ and $i\in{\rm O_{\rm
hi}}$ indicate when $i$th subfield is located in the low- and high-cadence
fields, respectively. Note that 4,057 bins have $N_{{\rm eve},i}^{\rm
obs}(t_{{\rm E},j})>0$ out of the total 29,120 bins (1456 subfields $\times$
20 bins) while the remaining have $N_{{\rm eve},i}^{\rm obs}(t_{{\rm
E},j})=0$.

The probability density function of $t_{\rm E}$, $f_{i}(t_{\rm E})$, is
calculated through a Monte Carlo simulation of $\sim 10^{8}$ microlensing
events for each $i$th subfield with a given model. This is computationally
expensive and time-consuming when $f_{i}(t_{\rm E})$ is updated for every
proposed model in the MCMC fitting procedure described in Section 4.3. Thus,
in a fitting run, we instead use a representative $f_{i}(t_{\rm E})$
calculated with a tentative best-fit model because we found that the
$\chi^{2}_{N_{\rm eve}}$ value had negligible variance when the $f_{i}(t_{\rm
E})$ used was calculated with a model showing similar $\chi^{2}_{t_{\rm E}}$
introduced below.

Mróz et al. (2017) and Mróz et al. (2019) separately provided $t_{\rm E}$
distributions consisting of 2212 events in the high-cadence fields and 5788
events in the low-cadence fields, respectively, which are denoted by $N_{\rm
O_{\rm hi}}^{\rm obs}(t_{{\rm E},j})$ and $N_{\rm O_{\rm low}}^{\rm
obs}(t_{{\rm E},j})$, respectively. These are given by

$N_{o}^{\rm obs}(t_{{\rm E},j})=\sum_{i\in o}N_{{\rm eve},i}^{\rm obs}(t_{{\rm
E},j})\hskip 19.91684pt(o={\rm O_{\rm hi}},{\rm O_{\rm low}}),$

where the summation is taken over all subfields located in $i\in o$ ($o={\rm
O_{\rm hi}}$ or ${\rm O_{\rm low}}$). The range of summation is different from
the case for $\chi^{2}_{N_{\rm eve}}$ where only the subfields covered by
Nataf et al. (2013) are considered because a shape of $t_{\rm E}$ distribution
is not sensitive to variation in extinction values. A compared model value for
a particular $j$th $t_{\rm E}$ bin is

$\displaystyle N_{o}^{\rm mod}(t_{{\rm
E},j})=\eta_{o}\left\langle\epsilon_{o}\right\rangle_{j}\sum_{i\in
o}\omega_{i}f_{i}(t_{{\rm E},j})\Delta t_{\rm E}\hskip 19.91684pt(o={\rm
O_{\rm hi}},{\rm O_{\rm low}}),$ (40)

where the mean detection efficiency
$\left\langle\epsilon_{o}\right\rangle_{j}$ and weight for each subfield
$\omega_{i}$ are

$\left\langle\epsilon_{o}\right\rangle_{j}\equiv\frac{\sum_{i\in o}N_{{\rm
eve},i}^{\rm obs}(t_{{\rm E},j})}{\sum_{i\in o}N_{{\rm eve},i}^{\rm
obs}(t_{{\rm E},j})/\epsilon_{i}(t_{{\rm E},j})},$

and

$\omega_{i}=\sum_{j}\frac{N_{{\rm eve},i}^{\rm obs}(t_{{\rm
E},j})}{\epsilon_{i}(t_{{\rm E},j})},$

respectively, with $\eta_{o}$ being an arbitrary constant expected to be
approximately 1. Although $\eta_{o}=1$ provides $\sum_{j}N_{o}^{\rm
mod}(t_{{\rm E},j})=\sum_{j}N_{o}^{\rm obs}(t_{{\rm E},j})$, $\eta_{o}$ can
take a random value because $\sum_{j}N_{o}^{\rm mod}(t_{{\rm E},j})$ is simply
an expected number of the total event detections, and $\sum_{j}N_{o}^{\rm
obs}(t_{{\rm E},j})$ can differ from it. Thus, $\eta_{o}$ is chosen such that
it minimizes $\chi^{2}_{t_{\rm E},o}$ introduced below in every step of the
MCMC fitting.

We consider $\chi^{2}$ for the two $t_{\rm E}$ distributions by

$\displaystyle\chi^{2}_{t_{\rm E}}$ $\displaystyle=\chi^{2}_{t_{\rm E},{\rm
O_{\rm hi}}}+\chi^{2}_{t_{\rm E},{\rm O_{\rm low}}}$
$\displaystyle=\sum_{j}\left(\chi^{2}_{t_{{\rm E},j},{\rm O_{\rm
hi}}}+\chi^{2}_{t_{{\rm E},j},{\rm O_{\rm low}}}\right),$ (41)

where

$\chi^{2}_{t_{{\rm E},j},o}=\begin{cases}-2\ln P[N_{o}^{\rm obs}(t_{{\rm
E},j});N_{o}^{\rm mod}(t_{{\rm E},j})]&\text{ when $N_{o}^{\rm obs}(t_{{\rm
E},j})\leq 15$}\\\ \left[\frac{N_{o}^{\rm mod}(t_{{\rm E},j})-N_{o}^{\rm
obs}(t_{{\rm E},j})}{N_{o}^{\rm err}(t_{{\rm E},j})}\right]^{2}&\text{ when
$N_{o}^{\rm obs}(t_{{\rm E},j})>15$}\end{cases}$ (42)

with

$N_{o}^{\rm err}(t_{{\rm
E},j})=\left\langle\epsilon_{o}\right\rangle_{j}\sqrt{\sum_{i\in
o}\frac{N_{{\rm eve},i}^{\rm obs}(t_{{\rm E},j})}{\epsilon_{i}^{2}(t_{{\rm
E},j})}}.$

Recall that $P(x;\lambda)=\lambda^{x}e^{-\lambda}/x!$ is the Poisson
probability of $x$ with the mean value $\lambda$. The uncertainty $N_{o}^{\rm
err}(t_{{\rm E},j})$ is not simply a square root of the sum of the Poisson
errors of each $N_{{\rm eve},i}^{\rm obs}(t_{{\rm E},j})$, but instead, it is
weighted by $1/\epsilon_{i}^{2}(t_{{\rm E},j})$. Therefore, this makes the
$N_{o}(t_{{\rm E},j})$ distribution different from a simple Poisson
probability distribution; however, we consider $\chi^{2}_{t_{{\rm E},j},o}$
with the Poisson distribution when $N_{o}^{\rm obs}(t_{{\rm E},j})$ is small
($\leq 15$) because the dependency of $\epsilon_{i}(t_{{\rm E},j})$ on $i$
(i.e., subfield) is much smaller than that on $j$ (i.e., $t_{\rm E}$), and the
distribution is expected to remain similar to the Poisson distribution,
especially when number of subfields contributing to $N_{o}^{\rm obs}(t_{{\rm
E},j})$ is small. Because the Poisson distribution can be approximated to the
Gaussian distribution when $N_{o}^{\rm obs}(t_{{\rm E},j})$ becomes larger, we
consider $\chi^{2}_{t_{{\rm E},j},o}$ with the Gaussian distribution with
$N_{o}^{\rm err}(t_{{\rm E},j})$ to consider the weight of
$1/\epsilon_{i}^{2}(t_{{\rm E},j})$.

###  4.2 Prior Constraints and $\chi^{2}$ Penalty

Table 3: Priors used in fits for the bulge model and posterior values.

| Parameter | Reference | Value/Prioraa $x\pm x_{e}$ indicates Gaussian prior and $\chi^{2}_{\rm pena}=\sum_{x}a_{x}\times(\frac{p_{x}-x}{x_{e}})^{2}$ is applied in the fit, where $a_{x}=5$ for $M_{\rm VVV}$ and $a_{x}=1$ for the others, and $p_{x}$ is the value of parameter $x$ in a given model. | PosteriorbbPosterior values are from the best-fit G+GX model in Table 5, and the uncertainties are combinations of statistic errors and systematic errors due to model choice.
---|---|---|---|---
Sun |  $(R_{\odot},z_{\odot})$ [pc] | (1), (2) | (8160, 25) | (8160, 25)
|  $(v_{\odot,x},v_{\odot,y},v_{\odot,z})$ [km/s] | (1) | (-10, 243, 7) | (-10, 243, 7)
Bulge |  $\alpha_{\rm bar}$ [deg.] | (1) | 27 | 27
|  $M_{\rm VVV}$ [$10^{10}M_{\odot}$] | (3) | $1.32\pm 0.08$ | $1.14^{+0.10}_{-0.11}$
|  $\left\langle\sigma_{v_{\rm B},x^{\prime}}\right\rangle$ [km/s] | (1), (4) | $135\pm 5$ | $141.0^{+1.7}_{-4.4}$
|  $\left\langle\sigma_{v_{\rm B},y^{\prime}}\right\rangle$ [km/s] | (1), (4) | $105\pm 5$ | $113.6^{+3.4}_{-0.9}$
|  $\left\langle\sigma_{v_{\rm B},z^{\prime}}\right\rangle$ [km/s] | (1), (4) | $96\pm 5$ | $108.3^{+0.2}_{-1.7}$
|  $\Omega_{\rm p}$ [km/s/kpc] | (3) | $39.0\pm 3.5$ | $45.9^{+4.0}_{-5.4}$
IMF |  $M_{\rm br}$ [$M_{\odot}$] | (5), (6) | $0.60\pm 0.10$ | $0.90^{+0.05}_{-0.14}$
| $\alpha_{\rm hm}$ | (5), (6) | $2.30\pm 0.10$ | $2.32^{+0.14}_{-0.10}$
| $\alpha_{\rm ms}$ | (5), (6) | $1.30\pm 0.15$ | $1.16^{+0.08}_{-0.15}$
| $\alpha_{\rm bd}$ | (5) | $0.30\pm 0.70$ | $0.22^{+0.20}_{-0.55}$

References. — (1) Bland-Hawthorn & Gerhard (2016); (2) Gravity Collaboration
et al. (2019); (3) Portail et al. (2017); (4) Portail et al. (2015); (5)
Kroupa (2001); (6) Calamida et al. (2015).

We have prior information about certain fundamental parameters in our Galaxy
from several previous studies, and reasonable prior constraints on such
parameters help us efficiently examine the huge parameter space or disentangle
degeneracies among fit parameters. Particularly, in this case where our
parametric model lacks a dynamical consistency, prior constraints on the bulge
mass or kinematic parameters help us avoid converging into a completely
unphysical model.

Table 4.2 lists the fundamental parameters on which we apply prior constraints
during the fit. Other than the solar position, velocity and bar angle that are
fixed, there are 9 parameters: the model integrated mass within the VVV bulge
box, $M_{\rm VVV}$, mass-weighted velocity dispersions inside the bulge half
mass radius, $\left\langle\sigma_{v_{\rm B},i}\right\rangle$
($i=x^{\prime},y^{\prime},z^{\prime}$), bar pattern speed, $\Omega_{\rm p}$,
three IMF slopes, ($\alpha_{\rm hm}$, $\alpha_{\rm ms}$, $\alpha_{\rm bd}$),
and break mass, $M_{\rm br}$, where the VVV bulge box is defined the central
region inside $(x^{\prime},y^{\prime},z^{\prime})=(\pm 2.2,\pm 1.4,\pm 1.2)$
kpc (Wegg & Gerhard, 2013). Practically, we add a $\chi^{2}$ penalty depending
on a deviation from the applied prior value for each parameter, as described
below.

Portail et al. (2017) found a total dynamical mass inside the VVV bulge box of
$(1.85\pm 0.05)\times 10^{10}~{}M_{\odot}$, where the mass budget was
$(0.32\pm 0.05)\times 10^{10}~{}M_{\odot}$ for dark matter, $0.2\times
10^{10}~{}M_{\odot}$ for a nuclear stellar disk in the very center, and the
remaining $(1.32\pm 0.08)\times 10^{10}~{}M_{\odot}$ was for stellar objects
traceable by the RC stars. Because the data introduced in Section 4.1 is not
sensitive to the central nuclear stellar disk with a $\sim 50~{}{\rm pc}$
scale height, we use $(1.32\pm 0.08)\times 10^{10}~{}M_{\odot}$ as a prior for
$M_{\rm VVV}$ and apply a $\chi^{2}$ penalty of $5\times(\frac{M_{\rm
VVV}\,[10^{10}\,M_{\odot}]-1.32}{0.08})^{2}$ in the fit. Our best-fit models
presented in Section 4.3 all have $M_{\rm VVV}<1.32\times
10^{10}~{}M_{\odot}$, and the large factor of 5 was multiplied by the
$\chi^{2}$ penalty to avoid needlessly reducing the bulge stellar mass,
because our parametric Galactic model does not ensure dynamical consistency by
itself. Note that the light stellar mass in the bulge region does not violate
dynamics in itself; rather it indicates an additional stellar mass in a non-
sensitive region and/or a larger dark matter mass fraction compared to the
Portail et al. (2017) model. We discuss the implied dark matter mass in our
model later, in Section 5.3.

For the mass-weighted velocity dispersions inside the bulge half mass radius,
$\left\langle\sigma_{v_{\rm B},i}\right\rangle$
($i=x^{\prime},y^{\prime},z^{\prime}$), we apply a $\chi^{2}$ penalty of
$\sum_{i=x^{\prime},y^{\prime},z^{\prime}}\left(\frac{\left\langle\sigma_{v_{\rm
B},i}\right\rangle-\left\langle\sigma_{v_{\rm B},i}\right\rangle_{\rm
P15}}{5{\rm km/s}}\right)^{2}$, where $(\left\langle\sigma_{v_{\rm
B},x^{\prime}}\right\rangle,\left\langle\sigma_{v_{\rm
B},y^{\prime}}\right\rangle,\left\langle\sigma_{v_{\rm
B},z^{\prime}}\right\rangle)_{\rm P15}=$ (135, 105, 96) km/s that were derived
by Bland-Hawthorn & Gerhard (2016) using the Portail et al. (2015) dynamical
bulge model. Further, a $\chi^{2}$ penalty of $(\frac{\Omega_{\rm
p}-39.0}{3.5})^{2}$ is applied for the pattern speed of $\Omega_{\rm p}$ in
unit of km/s/kpc based on the value derived by Portail et al. (2017).

For the IMF parameters, the following $\chi^{2}$ penalties were applied for
each slope and break mass; $(\frac{M_{\rm br}/M_{\odot}-0.60}{0.10})^{2}$,
$(\frac{\alpha_{\rm hm}-2.30}{0.10})^{2}$, $(\frac{\alpha_{\rm
ms}-1.30}{0.15})^{2}$, and $(\frac{\alpha_{\rm bd}-0.30}{0.70})^{2}$, where
the prior values were determined based on the local IMF by Kroupa (2001)
combined with the break mass at $\sim 0.56\,M_{\odot}$ found by Calamida et
al. (2015).

Hereafter, we denote the sum of the $\chi^{2}$ penalty values for the above
nine parameters by $\chi^{2}_{\rm pena}$.

###  4.3 Fitting

####  4.3.1 Fitting procedure

We use the Markov Chain Monte Carlo (MCMC) methods (Metropolis et al., 1953)
for our fitting. The bulge model consists of 39 fit parameters, i.e., 7–16,
19, and 4 parameters for the density ($\rho_{\rm B}$), velocity ($v_{\rm B}$),
and IMF models, respectively. In the previous sections, we defined the
following eight $\chi^{2}$ values, $\chi^{2}_{\rm RC}$ for the OGLE-III RC
star count data in 9019 lines of sight (Section 4.1.1), $\chi^{2}_{\mu}$ for
the VIRAC proper motion measurements in 586 VVV sub-tiles (Section 4.1.2),
$\chi^{2}_{\rm BRA}$ for the BRAVA radial velocity data in 82 fields (Section
4.1.3), $\chi^{2}_{N_{I<I_{\rm c}}}(I_{\rm c}=18,21)$ and
$\chi^{2}_{f_{18/21}}$ for the OGLE-IV star count data in 1456 subfields
(Section 4.1.4), $\chi^{2}_{N_{\rm eve}}$ for the number of microlensing event
detections as a function of $t_{\rm E}$ by the OGLE-IV survey in the 1456
subfields (Section 4.1.5), and $\chi^{2}_{t_{\rm E}}$ for the two $t_{\rm E}$
distributions in the OGLE-IV high- and low-cadence fields divided into 20
$\log[t_{\rm E}/{\rm days}]$ bins (Section 4.1.5), in addition to the
$\chi^{2}$ penalty from priors on the nine parameters summarized in Table 4.2
(Section 4.2). Of the eight $\chi^{2}$ values, $\chi^{2}_{\mu}$ and
$\chi^{2}_{\rm BRA}$ depend on the $\rho_{\rm B}$ and $v_{\rm B}$ models,
$\chi^{2}_{N_{I<I_{\rm c}}}(I_{\rm c}=18,21)$ depends on the $\rho_{\rm B}$
and IMF models, and $\chi^{2}_{N_{\rm eve}}$ and $\chi^{2}_{t_{\rm E}}$ depend
on the $\rho_{\rm B}$, $v_{\rm B}$, and IMF models. Thus, ideally, all the 39
fit parameters are simultaneously fitted because of the correlations among the
$\rho_{\rm B}$, $v_{\rm B}$, and IMF models in the parameter space.

However, such a simultaneous fit is difficult for the following reasons:
First, the volume of the parameter space is very large and difficult to
examine in a reasonable amount of time. Second, the probability density
function of $t_{\rm E}$, $f_{i}(t_{\rm E})$, is needed to calculate
$\chi^{2}_{N_{\rm eve}}$, which involves time-consuming calculations, for a
given combination of $\rho_{\rm B}$, $v_{\rm B}$, and IMF models. As described
in Section 4.1.5, a representative $f_{i}(t_{\rm E})$ can be used for models
showing similar $\chi^{2}_{t_{\rm E}}$, but it is difficult to ensure the
similarity in such a simultaneous fit. Third, even if time permitted, the
“best-fit” model depends on what $\chi^{2}$ value was minimized in the fit.
Because the number of data points in the used datasets ranges from 40 bins for
$\chi^{2}_{t_{\rm E}}$ to $9019\times 3$ for $\chi^{2}_{\rm RC}$, agreement to
a small dataset could be neglected depending on arbitrarily selected weights
among the eight $\chi^{2}$ values. The weights are ideally 1 if every data
does not contain any systematic error, and the model selection is perfect;
however, since neither of these statements is true, the “best-fit” model would
be dominated by an arbitrary choice of weights.

Table 4: Each step in the iterative process to find best-fit models.  Step | Fit parameters | Minimized $\chi^{2}$ | InputaaInput models or functions are fixed during each step other than $\rho_{\rm 0,B}$ in step 1. | Output
---|---|---|---|---
1 | $M_{\rm br},\alpha_{\rm hm},\alpha_{\rm ms},\alpha_{\rm bd},\rho_{\rm 0,B}$ | ${\tilde{\chi}^{2}_{\rm IMF}}\tablenotemark{b}+\chi^{2}_{\rm pena}$ | Tentative best $\rho_{\rm B}$ and $v_{\rm B}$ models |  $L_{{\cal M}_{I}}({\cal M}_{I})$, $f_{i}(t_{{\rm E},j})$
2 |  $\rho_{\rm 0,B},\bm{p}_{r_{s}},R_{c}$ | ${\tilde{\chi}^{2}_{\rho_{\rm B}}}\tablenotemark{b}+5\chi^{2}_{\rm pena}$ |  $L_{{\cal M}_{I}}({\cal M}_{I})$ & $f_{i}(t_{{\rm E},j})$ from 1 |  $\rho_{\rm B}$ model
| \+ $x_{\rm KM,th}$ for KM option |  |  |
| \+ $f_{\rm 0,X},\bm{p}_{r_{s},{\rm X}},R_{\rm c,X},b_{\rm X}$ for 2-comp. model |  |  |
3 | $\Omega_{\rm p},v_{0}^{\rm str},y_{0}^{\rm str},\sigma_{v_{\rm B},i,0},\sigma_{v_{\rm B},i,1},\bm{p}_{r_{s},\sigma_{R}},\bm{p}_{r_{s},\sigma_{z^{\prime}}}$ | $2\chi^{2}_{v_{\rm B}}\tablenotemark{b}+5\chi^{2}_{\rm pena}$ |  $\rho_{\rm B}$ model from 2 |  $v_{\rm B}$ model
4 | Same as step 2 | ${\tilde{\chi}^{2}_{\rho_{\rm B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$ |  $v_{\rm B}$ model from 3, $L_{{\cal M}_{I}}({\cal M}_{I})$ & $f_{i}(t_{{\rm E},j})$ from 1 |  $\rho_{\rm B}$ model
5 | Same as step 3 | $2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$ |  $\rho_{\rm B}$ model from 4 |  $v_{\rm B}$ model

Because of these difficulties associated with the simultaneous fit, we perform
modeling by iterating the following step-by-step procedure, where steps 1–5
are also summarized in Table 4:

  1. 0.

At the beginning, tentative best-fit $\rho_{\rm B}$ and $v_{\rm B}$ models are
determined to use in step 1 by performing steps 2–5 without using
$\chi^{2}_{t_{\rm E}}$, $\chi^{2}_{N_{I<21}}$, $\chi^{2}_{N_{I<18}}$,
$\chi^{2}_{f_{18/21}}$, or $\chi^{2}_{N_{\rm eve}}$. The Kroupa (2001) IMF is
used when an IMF is needed.

  2. 1.

An IMF model is determined by minimizing ${\tilde{\chi}^{2}_{\rm
IMF}}+\chi^{2}_{\rm pena}$, where $\chi^{2}$ to evaluate an IMF model is
defined as ${\tilde{\chi}^{2}_{\rm IMF}}\equiv\chi^{2}_{t_{\rm
E}}+0.2(\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+\chi^{2}_{f_{18/21}})/3$, and
the best-fit $\rho_{\rm B}$ and $v_{\rm B}$ models from the previous run or
from step 0 are fixed and used. Subsequently, the probability density function
of $t_{\rm E}$ in $i$th subfield, $f_{i}(t_{\rm E})$, is calculated using the
$\rho_{\rm B}$, $v_{\rm B}$, and IMF models. The luminosity function $L_{{\cal
M}_{I}}({\cal M}_{I})$ is calculated using the IMF model. A $\chi^{2}_{t_{\rm
E}}$ value, $\chi^{2}_{t_{\rm E},{\rm ini}}$, is calculated using the
$\rho_{\rm B}$, $v_{\rm B}$, and IMF models.

  3. 2.

A $\rho_{\rm B}$ model is determined by minimizing
${\tilde{\chi}^{2}_{\rho_{\rm B}}}+5\chi^{2}_{\rm pena}$, where $\chi^{2}$ to
evaluate a $\rho_{\rm B}$ model is defined as ${\tilde{\chi}^{2}_{\rho_{\rm
B}}}\equiv\chi^{2}_{\rm
RC}+\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+2\chi^{2}_{N_{\rm eve}}$, and the
$f_{i}(t_{\rm E})$ and $L_{{\cal M}_{I}}({\cal M}_{I})$ from step 1 are fixed
and used.

  4. 3.

A $v_{\rm B}$ model is determined by minimizing $2\chi^{2}_{v_{\rm
B}}+5\chi^{2}_{\rm pena}$, where $\chi^{2}$ to evaluate a $v_{\rm B}$ model is
defined as $\chi^{2}_{v_{\rm B}}\equiv\chi^{2}_{\mu}+\chi^{2}_{\rm BRA}$, and
the best-fit $\rho_{\rm B}$ model from step 2 is fixed and used.

  5. 4.

The $\rho_{\rm B}$ model from step 2 is updated by minimizing
${\tilde{\chi}^{2}_{\rho_{\rm B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm
pena}$, where the best-fit $v_{\rm B}$ model from step 3 is fixed and used.

  6. 5.

The $v_{\rm B}$ model from step 3 is updated by minimizing $2\chi^{2}_{v_{\rm
B}}+5\chi^{2}_{\rm pena}$, where the best-fit $\rho_{\rm B}$ model from step 4
is fixed and used.

  7. 6.

A $\chi^{2}_{t_{\rm E}}$ value, $\chi^{2}_{t_{\rm E},{\rm fin}}$, is
calculated using the $\rho_{\rm B}$, $v_{\rm B}$, and IMF models from steps 4,
5, and 1, respectively. If $|\chi^{2}_{t_{\rm E},{\rm ini}}-\chi^{2}_{t_{\rm
E},{\rm fin}}|<4$, the iteration is stopped; otherwise, steps 1–6 are
repeated.

It is difficult to systematically determine weights among various datasets
suffering from different systematic errors; therefore, a factor multiplied by
each $\chi^{2}$ value is subjectively selected in our attempt to find a proper
balance between large and small datasets. We multiply 2 by $\chi^{2}_{N_{\rm
eve}}$ in ${\tilde{\chi}^{2}_{\rho_{\rm B}}}$ because our interest is in
microlensing study, and we want to have a better agreement on it. Further, we
multiply 2 by $\chi^{2}_{v_{\rm B}}$ because the number of data points
contributed is $\sim 30$ times fewer than ${\tilde{\chi}^{2}_{\rho_{\rm B}}}$.
$\chi^{2}_{\rm pena}$ is multiplied by 5 in steps 2–5 because the
$\chi^{2}_{\rm pena}$ is the sum of the $\chi^{2}$ penalty set on the nine
quantities, as the effect is otherwise easily diminished by a small
improvement of $\chi^{2}_{v_{\rm B}}$ or ${\tilde{\chi}^{2}_{\rho_{\rm B}}}$,
which could be falsely caused due to systematic errors in data.

In step 1, we fit the four parameters for the IMF and the mass normalization
factor, $\rho_{\rm 0,B}$, by minimizing ${\tilde{\chi}^{2}_{\rm IMF}}$, which
concerns not only agreement with the $t_{\rm E}$ distribution by
$\chi^{2}_{t_{\rm E}}$, but also agreement with star count data via the
$\chi^{2}_{N_{I<21}}$, $\chi^{2}_{N_{I<18}}$, and $\chi^{2}_{f_{18/21}}$
values. This is because the model values for the star count data given by Eq.
(33) are calculated using the luminosity function $L_{{\cal M}_{I}}({\cal
M}_{I})$, and $L_{{\cal M}_{I}}({\cal M}_{I})$ is calculated using a given
IMF. The sum of the three $\chi^{2}$ values for the OGLE-IV star count data,
$\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+\chi^{2}_{f_{18/21}}$, are multiplied
by 0.1 $\times$ 2/3 in ${\tilde{\chi}^{2}_{\rm IMF}}$ because the number of
data values are 1456 for each of the three datasets compared with 40 for the
$t_{\rm E}$ distribution, and we desire to place a larger weight on $t_{\rm
E}$ because $\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}$ is also included in
${\tilde{\chi}^{2}_{\rho_{\rm B}}}$. The factor 2/3 is multiplied because only
two of the three datasets, $N_{I<21}$, $N_{I<18}$, and $f_{18/21}$
($=N_{I<18}/N_{I<21}$), are independent data.

As mentioned in Section 4.1.5, the $\chi^{2}$ for the number of microlensing
event detections, $\chi^{2}_{N_{\rm eve}}$, is insensitive to a variation of
the probability density function for $t_{\rm E}$ in $i$th subfield,
$f_{i}(t_{\rm E})$, till the models used to calculate $f_{i}(t_{\rm E})$ show
similar values of $\chi^{2}_{t_{\rm E}}$, the $\chi^{2}$ for the shape of the
$t_{\rm E}$ distribution. In the fitting procedure, we use the same
$f_{i}(t_{\rm E})$ from step 1 throughout steps 2–6, where $\chi^{2}_{N_{\rm
eve}}$ contributes to $\chi^{2}_{\rm\rho_{\rm B}}$ in steps 2 and 4. Using the
same $f_{i}(t_{\rm E})$ can be justified only when $\chi^{2}_{t_{\rm E},{\rm
ini}}$ is similar to $\chi^{2}_{t_{\rm E},{\rm fin}}$, and thus a similarity
between the two values is used as the condition to end the iteration. The
differences among different iteration runs are $f_{i}(t_{\rm E})$ and the IMF
model, and the condition ensures convergences of not only $f_{i}(t_{\rm E})$,
but also the IMF, because the IMF parameters are primarily determined by shape
of $t_{\rm E}$ distribution when the $\rho_{\rm B}$ and $v_{\rm B}$ models are
sufficiently constrained by other datasets. Note that $\chi^{2}_{N_{\rm eve}}$
is only used to find the best-fit $\rho_{\rm B}$ models because a model value
for the number of microlensing events with $t_{{\rm E},j}$ in $i$th subfield,
$N_{{\rm eve},i}^{\rm mod}(t_{{\rm E},j})$, given by Eq. (36) is independent
of the $v_{\rm B}$ and IMF models once $f_{i}(t_{\rm E})$ is fixed.

Figure 9: Surface density distributions along the bar-axes of the two one-
component models (E and G) and two two-components models (E+EX and G+GX).
Parameters of each model are listed in Table 5. Note that the plotted surface
density does not include the contribution from very centered or very thin
components such as a nuclear stellar disk (Launhardt et al., 2002; Nishiyama
et al., 2013; Portail et al., 2017) or thin bar components (Wegg et al.,
2015). We neither confirm nor confute existence of such components in this
paper because we lack data in $|b|\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2^{\circ}$ or in $|l|>10^{\circ}$.

We use a Monte Carlo simulation of $3\times 10^{6}$ microlensing events to
calculate the model values for the $t_{\rm E}$ distribution, $N_{o}^{\rm
mod}(t_{{\rm E},j})$ ($o={\rm O_{\rm hi}},{\rm O_{\rm low}}$), which takes
$\sim$20 s to calculate a $\chi^{2}_{t_{\rm E}}$ value under our computational
environment. This causes a dispersion of $\chi^{2}_{t_{\rm E}}$ values $\sim
1.3$ depending on the random seed value. Thus, we conduct MCMC fits in step 1
with a temperature of $T_{\rm MCMC}\geq 4$ to avoid being stuck at a local
minimum created due to the random noise, where a proposed parameter set with
$\Delta\chi^{2}>0$ is accepted with a probability
$\exp[{-0.5\,\Delta\chi^{2}/T_{\rm MCMC}}]$ in an MCMC run. As a result of
this, we use $|\chi^{2}_{t_{\rm E},{\rm ini}}-\chi^{2}_{t_{\rm E},{\rm
fin}}|<4$ as the condition in step 6. Note that the $\chi^{2}_{t_{\rm E},{\rm
ini}}$ and $\chi^{2}_{t_{\rm E},{\rm fin}}$ values are calculated with a Monte
Carlo simulation of $3\times 10^{7}$ events, which is 10 times the one during
the MCMC calculations.

It is important to ensure that a solution from the MCMC fit with the
$\chi^{2}_{t_{\rm E}}$ dispersion of $\sim 1.3$ is not plagued by the random
noise. In Section 4.4, we compare the best-fit model from our fitting
procedure with the best-fit model from a grid search where $\chi^{2}_{t_{\rm
E}}$ at each grid is calculated with $3\times 10^{7}$ events. We find a
consistency between them, and confirm that the settings for the fitting
procedure have successfully determined a converged solution.

In the Monte Carlo simulation of microlensing events, we also consider binary
systems for lens objects, which are generated by following the binary
distribution developed by Koshimoto et al. (2020) based on Duchêne & Kraus
(2013). Because Mróz et al. (2017) and Mróz et al. (2019) did not include
binary lens events in the OGLE-IV $t_{\rm E}$ distributions, we reject a
detectable binary system whose central caustic size $w_{\rm LC}$ is larger
than the impact parameter $u_{0}$, where $u_{0}$ is randomly generated
uniformly from 0 to 1 in each trial of the Monte Carlo simulation. We use
$w_{\rm LC}=4q_{\rm LC}/(s_{\rm LC}-s_{\rm LC}^{-1})^{2}$ (Chung et al., 2005)
with the mass ratio, $q_{\rm LC}$, and the separation in unit of the angular
Einstein radius $\theta_{\rm E}$, $s_{\rm LC}$. Although the formula for
$w_{\rm LC}$ is an approximation for the planetary mass-ratio of $q_{\rm
LC}\ll 1$, Koshimoto et al. (2020) found that the criterion $w_{\rm LC}<u_{0}$
itself works reasonably well even for a stellar mass-ratio because of the
extension of $\theta_{\rm E}$ by a factor $\sqrt{1+q_{\rm LC}}$. Section
10.3.1 of Koshimoto et al. (2020) presents a more detailed discussion. Because
a tight close binary system with $s_{\rm LC}\ll 1$ can be undetectable under
this criterion, the resulting $t_{\rm E}$ distribution includes their
contribution, as well as single star lens systems. Such a binary system has a
lens mass equal to the total system mass, and their fraction for each $t_{\rm
E}$ bin ranges from $\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}3\%$
for $t_{\rm E}<10~{}{\rm days}$ to 10–14% for $t_{\rm E}>50~{}{\rm days}$.

####  4.3.2 Best-Fit Models

Table 5: Best-fit parameters for each model.  Model |  | IMF model |  |  $\rho_{\rm B}$ model |
---|---|---|---|---|---
|  ${\tilde{\chi}^{2}_{\rm sum}}$aa ${\tilde{\chi}^{2}_{\rm sum}}\equiv{\tilde{\chi}^{2}_{\rho_{\rm B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$, ${\tilde{\chi}^{2}_{\rm IMF}}\equiv\chi^{2}_{t_{\rm E}}+0.2(\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+\chi^{2}_{f_{18/21}})/3$, ${\tilde{\chi}^{2}_{\rho_{\rm B}}}\equiv\chi^{2}_{\rm RC}+\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+2\chi^{2}_{N_{\rm eve}}$, and $\chi^{2}_{v_{\rm B}}\equiv\chi^{2}_{\mu}+\chi^{2}_{\rm BRA}$. |  ${\tilde{\chi}^{2}_{\rm IMF}}$aa ${\tilde{\chi}^{2}_{\rm sum}}\equiv{\tilde{\chi}^{2}_{\rho_{\rm B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$, ${\tilde{\chi}^{2}_{\rm IMF}}\equiv\chi^{2}_{t_{\rm E}}+0.2(\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+\chi^{2}_{f_{18/21}})/3$, ${\tilde{\chi}^{2}_{\rho_{\rm B}}}\equiv\chi^{2}_{\rm RC}+\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+2\chi^{2}_{N_{\rm eve}}$, and $\chi^{2}_{v_{\rm B}}\equiv\chi^{2}_{\mu}+\chi^{2}_{\rm BRA}$. | $M_{\rm br}$ | $\alpha_{\rm hm}$ | $\alpha_{\rm ms}$ | $\alpha_{\rm bd}$ |  |  ${\tilde{\chi}^{2}_{\rho_{\rm B}}}$aa ${\tilde{\chi}^{2}_{\rm sum}}\equiv{\tilde{\chi}^{2}_{\rho_{\rm B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$, ${\tilde{\chi}^{2}_{\rm IMF}}\equiv\chi^{2}_{t_{\rm E}}+0.2(\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+\chi^{2}_{f_{18/21}})/3$, ${\tilde{\chi}^{2}_{\rho_{\rm B}}}\equiv\chi^{2}_{\rm RC}+\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+2\chi^{2}_{N_{\rm eve}}$, and $\chi^{2}_{v_{\rm B}}\equiv\chi^{2}_{\mu}+\chi^{2}_{\rm BRA}$. | $\rho_{\rm 0,B}$ |  | $x_{0}$ | $y_{0}$ | $z_{0}$ | $C_{\perp}$ | $C_{\parallel}$ | $R_{\rm c}$
|  |  | [$M_{\odot}$] |  |  |  |  |  | [$M_{\odot}/$pc3] |  | [kpc] | [kpc] | [kpc] |  |  | [kpc]
E | 80173 | 355 | 0.84 | 2.31 | 1.10 | 0.18 |  | 75481 | 9.72 |  | 0.67 | 0.28 | 0.24 | 1.4 | 3.3 | 2.8
G | 91319 | 391 | 0.90 | 2.40 | 1.18 | 0.17 |  | 86105 | 2.43 |  | 1.03 | 0.46 | 0.40 | 2.0 | 4.0 | 4.8
E+EX | 76261 | 342 | 0.86 | 2.32 | 1.13 | 0.18 |  | 72147 | 4.12 |  | 0.93 | 0.37 | 0.24 | 1.2 | 4.1 | 2.6
G+GX | 76500 | 337 | 0.90 | 2.32 | 1.16 | 0.22 |  | 72585 | 0.88 |  | 1.56 | 0.72 | 0.49 | 1.2 | 3.1 | 2.8
Model |  $\rho_{\rm B}$ model |  |  $v_{\rm B}$ model |
| $f_{\rm 0,X}$ | $b_{\rm X}$ | $x_{\rm 0,X}$ | $y_{\rm 0,X}$ | $z_{\rm 0,X}$ | $C_{\rm\perp,X}$ |  | $C_{\rm\parallel,X}$ | $R_{\rm c,X}$ |  |  $\chi^{2}_{v_{\rm B}}$aa ${\tilde{\chi}^{2}_{\rm sum}}\equiv{\tilde{\chi}^{2}_{\rho_{\rm B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$, ${\tilde{\chi}^{2}_{\rm IMF}}\equiv\chi^{2}_{t_{\rm E}}+0.2(\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+\chi^{2}_{f_{18/21}})/3$, ${\tilde{\chi}^{2}_{\rho_{\rm B}}}\equiv\chi^{2}_{\rm RC}+\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}+2\chi^{2}_{N_{\rm eve}}$, and $\chi^{2}_{v_{\rm B}}\equiv\chi^{2}_{\mu}+\chi^{2}_{\rm BRA}$. | $\Omega_{\rm p}$ | $v_{0}^{\rm str}$ | $y_{0}^{\rm str}$ | $\sigma_{v_{\rm B},x^{\prime},0}$ | $\sigma_{v_{\rm B},y^{\prime},0}$
|  |  | [kpc] | [kpc] | [kpc] |  |  |  | [kpc] |  |  | [km/s/kpc] | [km/s] | [pc] | [km/s] | [km/s]
E | – | – | – | – | – | – |  | – | – |  | 2203 | 49.5 | 49 | 393 | 64 | 75
G | – | – | – | – | – | – |  | – | – |  | 2382 | 40.5 | 12 | 20 | 76 | 68
E+EX | 1.44 | 1.38 | 0.28 | 0.18 | 0.29 | 1.3 |  | 2.2 | 1.3 |  | 1908 | 47.4 | 43 | 407 | 64 | 76
G+GX | 3.00 | 1.38 | 0.76 | 0.31 | 0.40 | 1.2 |  | 1.3 | 5.2 |  | 1831 | 45.9 | 28 | 11 | 64 | 75
Model |  $v_{\rm B}$ model |
| $\sigma_{v_{\rm B},z^{\prime},0}$ | $\sigma_{v_{\rm B},x^{\prime},1}$ | $\sigma_{v_{\rm B},y^{\prime},1}$ | $\sigma_{v_{\rm B},z^{\prime},1}$ | $x_{\rm 0,\sigma_{R}}$ | $y_{\rm 0,\sigma_{R}}$ |  | $z_{\rm 0,\sigma_{R}}$ | $C_{\rm\perp,\sigma_{R}}$ |  | $C_{\rm\parallel,\sigma_{R}}$ | $x_{\rm 0,\sigma_{z^{\prime}}}$ | $y_{\rm 0,\sigma_{z^{\prime}}}$ | $z_{\rm 0,\sigma_{z^{\prime}}}$ | $C_{\rm\perp,\sigma_{z^{\prime}}}$ | $C_{\rm\parallel,\sigma_{z^{\prime}}}$
| [km/s] | [km/s] | [km/s] | [km/s] | [kpc] | [kpc] |  | [kpc] |  |  |  | [kpc] | [kpc] | [kpc] |  |
E | 72 | 156 | 84 | 86 | 0.82 | 9.29 |  | 0.86 | 3.8 |  | 1.0 | 0.51 | 2.90 | 2.19 | 3.0 | 1.0
G | 75 | 136 | 109 | 101 | 1.03 | 2.15 |  | 0.73 | 4.9 |  | 1.0 | 0.52 | 1.44 | 1.10 | 2.3 | 1.0
E+EX | 71 | 152 | 78 | 82 | 0.86 | 3.22 |  | 0.95 | 4.3 |  | 1.0 | 0.56 | 2.00 | 3.82 | 3.7 | 1.1
G+GX | 70 | 155 | 78 | 83 | 0.94 | 4.23 |  | 0.88 | 4.6 |  | 1.0 | 0.70 | 1.73 | 2.03 | 4.8 | 1.0

As described in Section 2.3.1, we consider four different shapes for the
$\rho_{\rm B}$ model: two one-component models (E and G models) and two two-
components models (E+EX and G+GX models). For each of the E, G, E+EX, and G+GX
models, we consider four options, corresponding to the four disk models in
Table 2. We derived the best-fit parameters for each option with each model by
the fitting procedure described in Section 4.3.1. We define another $\chi^{2}$
value, ${\tilde{\chi}^{2}_{\rm sum}}={\tilde{\chi}^{2}_{\rho_{\rm
B}}}+2\chi^{2}_{v_{\rm B}}+5\chi^{2}_{\rm pena}$, and use its difference,
$\Delta{\tilde{\chi}^{2}_{\rm sum}}$, to indicate the better model among them.
The $\chi^{2}$ for the $t_{\rm E}$ distributions, $\chi^{2}_{t_{\rm E}}$, is
not considered in ${\tilde{\chi}^{2}_{\rm sum}}$ because the information
contained in $\chi^{2}_{t_{\rm E}}$ is included in the $\chi^{2}$ for the
number of microlensing event detections as a function of $t_{\rm E}$,
$\chi^{2}_{N_{\rm eve}}$, in a more statistically proper style. Although this
inclusion is only for the fields covered by Nataf et al. (2013), we
conservatively avoid the partial overlap. In this subsection, we determine the
best-fit model for each of the E, G, E+EX, and G+GX models by selecting the
best option from the four based on comparisons using
$\Delta{\tilde{\chi}^{2}_{\rm sum}}$.

Among the four disk models in Table 2, we decide to use the all-$z$ $+$ flat
$z_{\rm d}^{\rm thin}$ model based on the following two comparisons. The first
comparison is between the all-$z$ and low-$z$ models, where we find
$\Delta{\tilde{\chi}^{2}_{\rm sum}}<100$ between the two, which is not
considered significant. Thus, we select the all-$z$ models because it is more
broadly applicable compared to the low-$z$ models optimized for bulge sky. The
second comparison is between the flat- and linear-scale height models, where
we find that the best-fit model with the flat-scale height model is preferable
to the linear-scale height model by $\Delta{\tilde{\chi}^{2}_{\rm sum}}>1000$.
Half of this preference comes from $\chi^{2}_{N_{I<21}}+\chi^{2}_{N_{I<18}}$
while the other half comes from $2\chi^{2}_{v_{\rm B}}$. These results barely
depend on the choice of E, G, E+EX, and G+GX models, and we select the all-$z$
$+$ flat $z_{\rm d}^{\rm thin}$ model as the fiducial disk model for all of
them.

In Appendix A, we consider another option of using a different dataset of the
OGLE-IV $t_{\rm E}$ measurements. In this dataset, the original light curve
data are the same as those by Mróz et al. (2017, 2019), but the Gaussian
Process model developed by Golovich et al. (2020) is applied for the $t_{\rm
E}$ measurements. Appendix A presents a comparison of $\chi^{2}$ values
between the best-fit models with the original $t_{\rm E}$ distributions by
Mróz et al. (2017, 2019) and the best-fit models with the $t_{\rm E}$
distributions with the Golovich et al. (2020) model, where we find that the
original Mróz et al.’s $t_{\rm E}$ distributions are favored with respect to
the ${\tilde{\chi}^{2}_{\rm sum}}$, $\chi^{2}_{N_{\rm eve}}$, and
$\chi^{2}_{t_{\rm E}}$ values.

Table 5 lists all the best-fit parameters for the fiducial four models. Figs.
4–8 present comparisons between the data and model values using the four
models from Table 5, where we show $\chi^{2}_{p}$ values for each parameter
$p$ in each panel of Figs. 4, 6, and 8, while Figs. 5 and 7 show the residuals
corresponding to Figs. 4 and 6, respectively. Fig. 9 shows the surface density
distributions of each best-fit model.

Among the four models, we find that the two-components models are more
favorable than the one-component models by $\Delta{\tilde{\chi}^{2}_{\rm
sum}}>3600$. The E model is preferred over the G model by
$\Delta{\tilde{\chi}^{2}_{\rm sum}}\sim 11000$ in the one-component models. By
contrast, the two two-component models (E+EX vs. G+GX) show very similar
${\tilde{\chi}^{2}_{\rm sum}}$ values, as described in Section 2.3.1.

###  4.4 Uncertainty Assessments for Fundamental Parameters

The far-right column in Table 4.2 presents the posterior values for the
fundamental parameters on which we applied the prior in the fits. The fiducial
values are from the best-fit G+GX model, while the uncertainties are
combinations of statistical and systematic errors for each parameter.

The systematic errors are taken from the variation of each value depending on
the model choice listed in Table 5. These are the stellar mass within the VVV
bulge box ($\pm 2.2\times\pm 1.4\times\pm 1.2$ kpc), $M_{\rm
VVV}=1.14^{+0.00}_{-0.11}|_{\rm sys.}\,\times 10^{10}\,M_{\odot}$, mass-
weighted velocity dispersions inside the bulge half mass radius along the bar
axes, $(\left\langle\sigma_{v_{\rm
B},x^{\prime}}\right\rangle,\left\langle\sigma_{v_{\rm
B},y^{\prime}}\right\rangle,\left\langle\sigma_{v_{\rm
B},z^{\prime}}\right\rangle)=(141.0^{+0.0}_{-3.8}|_{\rm
sys.},113.6^{+3.2}_{-0.0}|_{\rm sys.},108.3^{+0.0}_{-1.5}|_{\rm sys.})\,{\rm
km/s}$, bar pattern speed, $\Omega_{\rm p}=45.9^{+3.6}_{-5.4}|_{\rm
sys.}\,{\rm km/s/kpc}$, break mass in the IMF, $M_{\rm
br}=0.90^{+0.00}_{-0.06}|_{\rm sys.}\,M_{\odot}$, and IMF slopes for three
different mass regions, $(\alpha_{\rm hm},\alpha_{\rm ms},\alpha_{\rm
bd})=(2.32^{+0.08}_{-0.01}|_{\rm sys.},1.16^{+0.02}_{-0.06}|_{\rm
sys.},0.22^{+0.00}_{-0.05}|_{\rm sys.})$. The representative values are for
the best-fit G+GX model, and the 0 error values indicate that the G+GX model
has the largest or lowest values among the four models in Table 5.

Statistical errors for the four velocity parameters are determined using the
posterior distribution for the G+GX model taken from the MCMC calculation in
step 5 in our fitting procedure that is described in Section 4.3.1. These are
$(\left\langle\sigma_{v_{\rm
B},x^{\prime}}\right\rangle,\left\langle\sigma_{v_{\rm
B},y^{\prime}}\right\rangle,\left\langle\sigma_{v_{\rm
B},z^{\prime}}\right\rangle)=(141.0^{+1.7}_{-2.2}|_{\rm
stat.},113.6^{+1.0}_{-0.9}|_{\rm stat.},108.3^{+0.2}_{-0.9}|_{\rm
stat.})\,{\rm km/s}$ and $\Omega_{\rm p}=45.9^{+0.6}_{-0.2}|_{\rm stat.}\,{\rm
km/s/kpc}$.

Although we have a posterior distribution from the MCMC calculation in step 1
for the other five parameters, $M_{\rm VVV}$, $M_{\rm br}$, $\alpha_{\rm hm}$,
$\alpha_{\rm ms}$, and $\alpha_{\rm bd}$, the distribution is contaminated by
the dispersion of $\chi^{2}_{t_{\rm E}}$ values of $\sim 1.3$ due to the
limited sample size of simulated events ($3\times 10^{6}$ events) in each step
of the MCMC calculation. Thus, we use the $\Delta\chi^{2}$ distributions from
a grid search with a ten times larger sample size of $3\times 10^{7}$ events
to determine statistical errors for the five parameters. Fig. 10 shows the
resultant $\Delta({\tilde{\chi}^{2}_{\rm IMF}}+\chi^{2}_{\rm pena})$ maps
calculated with the best-fit G+GX $\rho_{\rm B}$ and $v_{\rm B}$ models. In
the grid search, a grid is specified by a combination of $(M_{\rm
br},\alpha_{\rm hm},\alpha_{\rm ms},\alpha_{\rm bd})$ distributed uniformly,
and the mass normalization factor $\rho_{\rm 0,B}$, which is represented by
$M_{\rm VVV}$ here, is adjusted such that the ${\tilde{\chi}^{2}_{\rm
IMF}}+\chi^{2}_{\rm pena}$ for each grid is minimized.

The best-fit parameters from the grid search are $(M_{\rm VVV},M_{\rm
br},\alpha_{\rm hm},\alpha_{\rm ms},\alpha_{\rm bd})$ = ($1.16\times
10^{10}\,M_{\odot}$, $0.90\,M_{\odot}$, 2.30, 1.15, 0.20), which are almost
identical to the best-fit G+GX model parameters of ($1.14\times
10^{10}\,M_{\odot}$, $0.90\,M_{\odot}$, 2.32, 1.16, 0.22) considering the
limited resolution of the grid search. These two sets of parameters can be, in
principle, different because the grid search uses the $\rho_{\rm B}$ and
$v_{\rm B}$ models obtained from the last iteration of our fitting procedure,
whereas the IMF parameters of the best-fit G+GX model are determined using the
tentative best-fit $\rho_{\rm B}$ and $v_{\rm B}$ models from the last but one
iteration. Therefore, the identity of these two sets of parameters confirms a
convergence of our iterative fitting procedure. This also ensures that the
number of simulated events used in the MCMC fit, $3\times 10^{6}$, was
sufficient to find the best-fit parameters.

We use $\Delta({\tilde{\chi}^{2}_{\rm IMF}}+\chi^{2}_{\rm pena})=4$ as the
threshold value to determine the statistical errors for the five parameters
because the grid search is based on the best-fit $\rho_{\rm B}$ and $v_{\rm
B}$ models determined using the IMF parameters from step 1 in which the MCMC
temperature of $T_{\rm MCMC}=4$ is used. This yields $(M_{\rm VVV},M_{\rm
br},\alpha_{\rm hm},\alpha_{\rm ms},\alpha_{\rm bd})$ =
($1.14^{+0.10}_{-0.02}|_{\rm stat.}\times 10^{10}\,M_{\odot}$,
$0.90^{+0.05}_{-0.13}|_{\rm stat.}\,M_{\odot}$, $2.32^{+0.07}_{-0.10}|_{\rm
stat.}$, $1.16^{+0.08}_{-0.14}|_{\rm stat.}$, $0.22^{+0.20}_{-0.55}|_{\rm
stat.}$). By combining the statistical and systematic errors above, we derive
the posterior values listed in Table 4.2. Note that these error estimates
could be influenced by our non-standard fitting method, an iterative step-by-
step procedure, described in Section 4.3.1.

Figure 10: $\Delta({\tilde{\chi}^{2}_{\rm IMF}}+\chi^{2}_{\rm pena})$ map
calculated through a grid-search with the best-fit G+GX $\rho_{\rm B}$ and
$v_{\rm B}$ models. $M_{\rm VVV,pri}$ indicates prior value for $M_{\rm VVV}$
of 1.32 $\times 10^{10}\,M_{\odot}$.

##  5 Discussion

In this section, we discuss several fundamental parameters in our Galaxy
measured through our modeling process throughout Sections 3 and 4 by
comparison to previous studies. In Section 5.1, we discuss the 10 best-fit
parameters for the $v_{\rm d}$ models in Table 2. Section 5.2 examines how the
IMF parameters are constrained in our model framework, and compares these IMF
parameters with a local IMF by Kroupa (2001) and a bulge IMF by Zoccali et al.
(2000). In Section 5.3, we derive the stellar mass-to-light ratio
corresponding to our IMF and estimate the dark matter fraction inferred from
our results.

###  5.1 Disk Velocity Parameters

A review paper by Bland-Hawthorn & Gerhard (2016) summarizes estimates on the
local velocity dispersion values and the scale lengths of the distributions
from several previous studies (see their Table 5). The local velocity
dispersion values for the thin disk range from 34 (Piffl et al., 2014) to 48
km/s (Sanders & Binney, 2015) for $\sigma_{R,\odot}^{\rm thin}$ and from 20
(Binney, 2012) to 31 km/s (Sanders & Binney, 2015) for $\sigma_{z,\odot}^{\rm
thin}$, and our best-fit values in Table 2 are consistent with the estimates
for these two parameters in any of the four models considered.

Binney (2012) is the only one who did a fit to the scale length of the
velocity dispersion along $R$ for the thin disk, $R_{\sigma_{R}}^{\rm thin}$,
in the summary by Bland-Hawthorn & Gerhard (2016). Their estimate is
$R_{\sigma_{R}}^{\rm thin}=$ 3.3–20 kpc, which includes our estimated values
for all models other than the all-$z$ $+$ linear $z_{\rm d}^{\rm thin}$ model.
The value for the all-$z$ $+$ linear $z_{\rm d}^{\rm thin}$ model is
$R_{\sigma_{R}}^{\rm thin}=21.4$ kpc, and it is only slightly outside the
range presented by Binney (2012).

No study has presented a fit to the scale length of the velocity dispersion
along $z$ for the thin disk, $R_{\sigma_{z}}^{\rm thin}$, in the summary by
Bland-Hawthorn & Gerhard (2016), but some studies fixed the value at 9.0
(Piffl et al., 2014), 7.8 (Sanders & Binney, 2015), and 7.4 kpc (Binney &
Piffl, 2015) in their analysis. Our estimates on $R_{\sigma_{z}}^{\rm thin}$
ranging from 5.9–10.4 kpc are consistent with these values. The all-$z$ $+$
flat $z_{\rm d}^{\rm thin}$ model, which is our fiducial disk model selected
in Section 4.3.2, has $R_{\sigma_{z}}^{\rm thin}=5.9$ kpc which is shorter
than the fixed values in the previous papers, though.

Further, our estimates on the slope of age-velocity dispersion relation are
$\beta_{R}=0.22$–0.34 and $\beta_{z}=0.77$–0.82, which are fully and fairly
consistent with the estimates by Yu & Liu (2018) of $\beta_{R}=0.28\pm 0.08$
and $\beta_{z}=0.54\pm 0.13$, respectively.

By contrast, our estimates for some parameters for the thick disk show
inconsistency with previous studies. In Table 5 of Bland-Hawthorn & Gerhard
(2016), the local velocity dispersion values for the thick disk range from 25
(Binney, 2012) to 53 km/s (Binney & Piffl, 2015) for $\sigma_{R,\odot}^{\rm
thick}$ and from 33 to 65 km/s (Binney, 2012) for $\sigma_{z,\odot}^{\rm
thick}$, while our estimate for these parameters are $\sigma_{R,\odot}^{\rm
thick}=$ 68–75 km/s and $\sigma_{z,\odot}^{\rm thick}=47.8$–61.4 km/s.