In the context of forthcoming galaxy surveys, to ensure unbiased constraints on cosmology and gravity when using non-linear structure information, per cent-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 cold dark matter (Λ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 Λ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 4 per cent accurate down to at least k≈3hMpc−1 for a modification to gravity of |fR0| ≤ 10−5 and for the total neutrino mass Mν ≡ ∑mν ≤ 0.15 eV. We also find that the framework is 4 per cent consistent with EuclidEmulator2 as well as the Bacco emulator for most of the considered νwCDM cosmologies down to at least k≈3h Mpc−1. Finally, we compare against hydrodynamical simulations employing HMCode2020’s baryonic feedback modelling on top of the halo model reaction. For νΛCDM cosmologies, we find 2 per cent accuracy for Mν ≤ 0.48 eV down to at least k ≈ 5h Mpc−1. Similar accuracy is found when comparing to νwCDM hydrodynamical simulations with Mν = 0.06 eV. This offers the first non-linear, theoretically general means of accurately including massive neutrinos for beyond-ΛCDM cosmologies, and further suggests that baryonic, massive neutrino, and dark energy physics can be reliably modelled independently
