  Aiming at discriminating different gravitational potential models of the
Milky Way, we perform tests based on the kinematic data powered by the Gaia DR2
astrometry, over a large range of $(R,z)$ locations. Invoking the complete form
of Jeans equations that admit three integrals of motion, we use the independent
$R$- and $z$-directional equations as two discriminators ($T_R$ and $T_z$). We
apply the formula for spatial distributions of radial and vertical velocity
dispersions proposed by Binney et al., and successfully extend it to azimuthal
components, $\sigma_\theta(R,z)$ and $V_\theta(R,z)$; the analytic form avoids
the numerical artifacts caused by numerical differentiation in Jeans-equations
calculation given the limited spatial resolutions of observations, and more
importantly reduces the impact of kinematic substructures in the Galactic disk.
It turns out that whereas the current kinematic data are able to reject
Moffat's Modified Gravity (let alone the Newtonian baryon-only model),
Milgrom's MOND is still not rejected. In fact, both the carefully calibrated
fiducial model invoking a spherical dark matter (DM) halo and MOND are equally
consistent with the data at almost all spatial locations (except that probably
both have respective problems at low-$|z|$ locations), no matter which a tracer
population or which a meaningful density profile is used. Because there is no
free parameter at all in the quasi-linear MOND model we use, and the baryonic
parameters are actually fine-tuned in the DM context, such an effective
equivalence is surprising, and might be calling forth a transcending synthesis
of the two paradigms.
