Enabling matter power spectrum emulation in beyond-$\Lambda$CDM cosmologies with COLA

We compare and validate COLA (COmoving Lagrangian Acceleration) simulations against existing emulators in the literature, namely Bacco and Euclid Emulator 2. Our analysis focuses on the non-linear response function, i.e., the ratio between the non-linear dark matter power spectrum in a given cosmology with respect to a pre-defined reference cosmology, which is chosen to be the Euclid Emulator 2 reference cosmology in this paper. We vary three cosmological parameters, the total matter density, the amplitude of the primordial scalar perturbations and the spectral index. By comparing the COLA non-linear response function with those computed from each emulator in the redshift range $0 \leq z \leq 3$, we find that the COLA method is in excellent agreement with the two emulators for scales up to $k \sim 1 \ h$/Mpc as long as the deviations of the matter power spectrum from the reference cosmology are not too large. We validate the implementation of massive neutrinos in our COLA simulations by varying the sum of neutrino masses to three different values, $0.0$ eV, $0.058$ eV and $0.15$ eV. We show that all three non-linear prescriptions used in this work agree at the $1\%$ level at $k \leq 1 \ h$/Mpc. We then introduce the Effective Field Theory of Dark Energy in our COLA simulations using the $N$-body gauge method. We consider two different modified gravity models in which the growth of structure is enhanced or suppressed at small scales, and show that the response function with respect to the change of modified gravity parameters depends weakly on cosmological parameters in these models.