One line to run them all: SuperEasy massive neutrino linear response in $N$-body simulations

We present in this work a novel and yet extremely simple method for incorporating the effects of massive neutrinos in cosmological $N$-body simulations. This so-called "SuperEasy linear response" approach is based upon analytical solutions to the collisionless Boltzmann equation in the clustering and free-streaming limits, which are then connected by a rational function interpolation function with cosmology-dependent coefficients given by simple algebraic expressions of the cosmological model parameters. The outcome is a {\it one-line modification} to gravitational potential that requires only the cold matter density contrast as a real-time input, and that can be incorporated into any $N$-body code with a Particle--Mesh component with no additional implementation cost. To demonstrate its power, we implement the SuperEasy method in the publicly available \gadgetcode{} code, and show that for neutrino mass sums not exceeding $\sum m_\nu \simeq 1$ eV, the total matter and cold matter power spectra are in sub-1% and sub-0.1% agreement with those from state-of-the-art linear response simulations in literature. Aside from its minimal implementation cost, compared with existing massive neutrino simulation methods, the SuperEasy approach has better memory efficiency, incurs no runtime overhead relative to a standard $\Lambda$CDM simulation, and requires no post-processing. The minimal nature of the method allows limited computational resources to be diverted to modelling other physical effects of interest, e.g., baryonic physics via hydrodynamics.