  In the next decade, cosmological surveys will have the statistical power to
detect the absolute neutrino mass scale. N-body simulations of large-scale
structure formation play a central role in interpreting data from such surveys.
Yet these simulations are Newtonian in nature. We provide a quantitative study
of the limitations to treating neutrinos, implemented as N-body particles, in
N-body codes, focusing on the error introduced by neglecting special
relativistic effects. Special relativistic effects are potentially important
due to the large thermal velocities of neutrino particles in the simulation
box. We derive a self-consistent theory of linear perturbations in Newtonian
and non-relativistic neutrinos and use this to demonstrate that N-body
simulations overestimate the neutrino free-streaming scale, and cause errors in
the matter power spectrum that depend on the initial redshift of the
simulations. For $z_{i} \lesssim 100$, and neutrino masses within the currently
allowed range, this error is $\lesssim 0.5\%$, though represents an up to $\sim
10\%$ correction to the shape of the neutrino-induced suppression to the cold
dark matter power spectrum. We argue that the simulations accurately model
non-linear clustering of neutrinos so that the error is confined to linear
scales.
