  Dwarf spheroidal (dSph) galaxies are believed to be strongly dark matter
dominated and thus are considered perfect objects to study dark matter
distribution and test theories of structure formation. They possess resolved,
multiple stellar populations that offer new possibilities for modeling. A
promising tool for the dynamical modeling of these objects is the Schwarzschild
orbit superposition method. In this work we extend our previous implementation
of the scheme to include more than one population of stars and a more general
form of the mass-to-light ratio function. We tested the improved approach on a
nearly spherical, gas-free galaxy formed in the cosmological context from the
Illustris simulation. We modeled the binned velocity moments for stars split
into two populations by metallicity and demonstrate that in spite of larger
sampling errors the increased number of constraints leads to significantly
tighter confidence regions on the recovered density and velocity anisotropy
profiles. We then applied the method to the Fornax dSph galaxy with stars
similarly divided into two populations. In comparison with our earlier work, we
find the anisotropy parameter to be slightly increasing, rather than
decreasing, with radius and more strongly constrained. We are also able to
infer anisotropy for each stellar population separately and find them to be
significantly different.
