An accurate treatment of scattering and diffusion in piecewise power-law
  models for cosmic ray and radiation/neutrino transport
  A popular numerical method to model the dynamics of a 'full spectrum' of
cosmic rays (CRs), also applicable to radiation/neutrino hydrodynamics (RHD),
is to discretize the spectrum at each location/cell as a piecewise power law in
'bins' of momentum (or frequency) space. This gives rise to a pair of conserved
quantities (e.g. CR number and energy) which are exchanged between cells or
bins, that in turn give the update to the normalization and slope of the
spectrum in each bin. While these methods can be evolved exactly in
momentum-space (considering injection, absorption, continuous losses/gains),
numerical challenges arise dealing with spatial fluxes, if the scattering rates
depend on momentum. This has often been treated by either by neglecting
variation of those rates 'within the bin,' or sacrificing conservation --
introducing significant errors. Here, we derive a rigorous treatment of these
terms, and show that the variation within the bin can be accounted for
accurately with a simple set of scalar correction coefficients that can be
written entirely in terms of other, explicitly-evolved 'bin-integrated'
quantities. This eliminates the relevant errors without added computational
cost, has no effect on the numerical stability of the method, and retains
manifest conservation. We derive correction terms both for methods which
explicitly integrate flux variables (two-moment or M1-like) methods, as well as
single-moment (advection-diffusion, FLD-like) methods, and approximate
corrections valid in various limits.
