On the road to percent accuracy V: the non-linear power spectrum beyond $\Lambda$CDM with massive neutrinos and baryonic feedback

In the context of forthcoming galaxy surveys, to ensure unbiased constraints on cosmology and gravity when using non-linear structure information, percent-level accuracy is required when modelling the power spectrum. This calls for frameworks that can accurately capture the relevant physical effects, while allowing for deviations from $\Lambda$CDM. Massive neutrino and baryonic physics are two of the most relevant such effects. We present an integration of the halo model reaction frameworks for massive neutrinos and beyond-$\Lambda$CDM cosmologies. The integrated halo model reaction, combined with a pseudo power spectrum modelled by HMCode2020 is then compared against $N$-body simulations that include both massive neutrinos and an $f(R)$ modification to gravity. We find that the framework is 5% accurate down to at least $k\approx 3 \, h/{\rm Mpc}$ for a modification to gravity of $|f_{\rm R0}|\leq 10^{-5}$ and for the total neutrino mass $M_\nu \equiv \sum m_\nu \leq 0.15$ eV. We also find that the framework is 4(1)% consistent with the Bacco (EuclidEmulator2) emulator for $\nu w$CDM cosmologies down to at least $k \approx 3 \, h$/Mpc. Finally, we compare against hydrodynamical simulations employing HMCode2020's baryonic feedback modelling on top of the halo model reaction. For $\nu \Lambda$CDM cosmologies we find 2% accuracy for $M_\nu \leq 0.48$eV down to at least $k\approx 5h$/Mpc. Similar accuracy is found when comparing to $\nu w$CDM hydrodynamical simulations with $M_\nu = 0.06$eV. This offers the first non-linear and theoretically general means of accurately including massive neutrinos for beyond-$\Lambda$CDM cosmologies, and further suggests that baryonic effects can be reliably modelled independently of massive neutrino and dark energy physics. These extensions have been integrated into the publicly available ReACT code.