  Boltzmann solvers are an important tool for the computation of cosmological
observables in the linear regime. In the presence of massive neutrinos, they
involve solving the Boltzmann equation followed by an integration in momentum
space to arrive at the desired fluid properties, a procedure which is known to
be computationally slow. In this work we introduce the so-called generalized
Boltzmann hierarchy (GBH) for massive neutrinos in cosmology, an alternative to
the usual Boltzmann hierarchy, where the momentum dependence is integrated out
leaving us with a two-parameter infinite set of ordinary differential
equations. Along with the usual expansion in multipoles, there is now also an
expansion in higher velocity weight integrals of the distribution function.
Using a toy code, we show that the GBH produces the density contrast neutrino
transfer function to a $\lesssim 0.5\%$ accuracy at both large and intermediate
scales compared to the neutrino free-streaming scale, thus providing a
proof-of-principle for the GBH. We comment on the implementation of the GBH in
a state of the art Boltzmann solver.
