  We present a Bayesian method to identify multiple (chemodynamic) stellar
populations in dwarf spheroidal galaxies (dSphs) using velocity, metallicity,
and positional stellar data without the assumption of spherical symmetry. We
apply this method to a new Keck/DEIMOS spectroscopic survey of the Ursa Minor
(UMi) dSph. We identify 892 likely members, making this the largest UMi sample
with line-of-sight velocity and metallicity measurements. Our Bayesian method
detects two distinct chemodynamic populations with high significance
($\ln{B}\sim33$). The metal-rich ($[{\rm Fe/H}]=-2.05\pm0.03$) population is
kinematically colder (radial velocity dispersion of $\sigma_v=4.9\pm0.8 \, {\rm
km \, s^{-1}}$) and more centrally concentrated than the metal-poor ($[{\rm
Fe/H}]=-2.29\pm0.05$) and kinematically hotter population ($\sigma_v
=11.5\pm0.9\, {\rm km \, s^{-1}}$). Furthermore, we apply the same analysis to
an independent MMT/Hectochelle data set and confirm the existence of two
chemodynamic populations in UMi. In both data sets, the metal-rich population
is significantly flattened ($\epsilon=0.75\pm0.03$) and the metal-poor
population is closer to spherical ($\epsilon=0.33_{-0.09}^{+0.12}$). Despite
the presence of two populations, we are unable to robustly estimate the slope
of the dynamical mass profile. We found hints for prolate rotation of order
$\sim 2 \, {\rm km \, s^{-1}}$ in the MMT data set, but further observations
are required to verify this. The flattened metal-rich population invalidates
assumptions built into simple dynamical mass estimators, so we computed new
astrophysical dark matter annihilation (J) and decay profiles based on the
rounder, hotter metal-poor population and inferred
$\log_{10}{(J(0.5^{\circ})/{\rm GeV^{2} \, cm^{-5}})}\approx19.1$ for the Keck
data set. Our results paint a more complex picture of the evolution of Ursa
Minor than previously discussed.
