  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.
