  We present population-orbit superposition models for external galaxies based
on Schwarzschild's orbit-superposition method, by tagging the orbits with age
and metallicity. The models fit the density distributions, as well as
kinematic, age and metallicity maps from Integral Field Unit (IFU) spectroscopy
observations. We validate the method and demonstrate its power by applying it
to mock data, similar to those obtained by the Multi-Unit Spectroscopic
Explorer (MUSE) IFU on the Very Large Telescope (VLT). These mock data are
created from Auriga galaxy simulations, viewed at three different inclination
angles ($\vartheta=40^o, 60^o, 80^o$). Constrained by MUSE-like mock data, our
model can recover the galaxy's stellar orbit distribution projected in orbital
circularity $\lambda_z$ vs. radius $r$, the intrinsic stellar population
distribution in age $t$ vs. metallicity $Z$, and the correlation between
orbits' circularity $\lambda_z$ and stellar age $t$. A physically motivated
age-metallicity relation improves recovering the intrinsic stellar population
distributions. We decompose galaxies into cold, warm and hot + counter-rotating
components based on their orbit circularity distribution, and find that the
surface density, mean velocity, velocity dispersion, age and metallicity maps
of each component from our models well reproduce those from simulation,
especially for projections close to edge-on. These galaxies exhibit strong
global age vs. $\sigma_z$ relation, which is well recovered by our model. The
method has the power to reveal the detailed build-up of stellar structures in
galaxies, and offers a complement to local resolved, and high-redshift studies
of galaxy evolution.
