  In this paper we present COMET, a Gaussian process emulator of the galaxy
power spectrum multipoles in redshift-space. The model predictions are based on
one-loop perturbation theory and we consider two alternative descriptions of
redshift-space distortions: one that performs a full expansion of the real- to
redshift-space mapping, as in recent effective field theory models, and another
that preserves the non-perturbative impact of small-scale velocities by means
of an effective damping function. The outputs of COMET can be obtained at
arbitrary redshifts (up to $z \sim 3$), for arbitrary fiducial background
cosmologies, and for a large parameter space that covers the shape parameters
$\omega_c$, $\omega_b$, and $n_s$, as well as the evolution parameters $h$,
$A_s$, $\Omega_K$, $w_0$, and $w_a$. This flexibility does not impair COMET's
accuracy, since we exploit an exact degeneracy between the evolution parameters
that allows us to train the emulator on a significantly reduced parameter
space. While the predictions are sped up by at least two orders of magnitude,
validation tests reveal an accuracy of $0.1\,\%$ for the monopole and
quadrupole ($0.3\,\%$ for the hexadecapole), or alternatively, better than
$0.25\,\sigma$ for all three multipoles in comparison to statistical
uncertainties expected for the Euclid survey with a tenfold increase in volume.
We show that these differences translate into shifts in mean posterior values
that are at most of the same size, meaning that COMET can be used with the same
confidence as the exact underlying models. COMET is a publicly available Python
package that also provides the tree-level bispectrum multipoles in
redshift-space and Gaussian covariance matrices.
