  We present the first simulations evolving resolved spectra of cosmic rays
(CRs) from MeV-TeV energies (including electrons, positrons, (anti)protons, and
heavier nuclei), in live kinetic-MHD galaxy simulations with star formation and
feedback. We utilize new numerical methods including terms often neglected in
historical models, comparing Milky Way analogues with phenomenological
scattering coefficients $\nu$ to Solar-neighborhood (LISM) observations
(spectra, B/C, $e^{+}/e^{-}$, $\bar{p}/p$, $^{10}$Be/$^{9}$Be, ionization,
$\gamma$-rays). We show it is possible to reproduce observations with simple
single-power-law injection and scattering coefficients (scaling with rigidity
$R$), similar to previous (non-dynamical) calculations. We also find: (1) The
circum-galactic medium in realistic galaxies necessarily imposes a
$\sim10\,$kpc CR scattering halo, influencing the required $\nu(R)$. (2)
Increasing the normalization of $\nu(R)$ re-normalizes CR secondary spectra but
also changes primary spectral slopes, owing to source distribution and loss
effects. (3) Diffusive/turbulent reacceleration is unimportant and generally
sub-dominant to gyroresonant/streaming losses, which are sub-dominant to
adiabatic/convective terms dominated by $\sim0.1-1\,$kpc turbulent/fountain
motions. (4) CR spectra vary considerably across galaxies; certain features can
arise from local structure rather than transport physics. (5) Systematic
variation in CR ionization rates between LISM and molecular clouds (or Galactic
position) arises naturally without invoking alternative sources. (6) Abundances
of CNO nuclei require most CR acceleration occurs around when reverse shocks
form in SNe, not in OB wind bubbles or later Sedov-Taylor stages of SNe
remnants.
