  We present a Markov chain Monte Carlo pipeline which can be used for unbiased
large-scale tests of gravity using galaxy cluster observations. The pipeline,
which currently uses cluster number counts to constrain the present-day
background scalar field $f_{R0}$ of Hu-Sawicki $f(R)$ gravity, fully accounts
for the effects of the fifth force on cluster properties including the
dynamical mass, the halo concentration and the observable-mass scaling
relations. This is achieved using general models which have been calibrated
over a wide and continuous mass range ($10^{11}M_{\odot}\lesssim
M\lesssim10^{15}M_{\odot}$) using a large suite of cosmological simulations,
including the first to simultaneously incorporate both screened modified
gravity and full baryonic physics. We show, using mock cluster catalogues, that
an incomplete treatment of the observable-mass scaling relations in $f(R)$
gravity, which do not necessarily follow the usual power-law behaviour, can
lead to unbiased and imprecise constraints. It is therefore essential to fully
account for these effects in future cosmological tests of gravity that will
make use of vast cluster catalogues from ongoing and upcoming galaxy surveys.
Our constraint framework can be easily extended to other gravity models; to
demonstrate this, we have carried out a similar modelling of cluster properties
in the normal-branch Dvali-Gabadadze-Porrati model (nDGP), which features a
very different screening mechanism. Using our full-physics simulations, we also
study the angular power spectra of the thermal and kinetic Sunyaev-Zel'dovich
effects in $f(R)$ gravity and nDGP, and demonstrate the potential for precise
constraints of gravity using data from upcoming CMB experiments. Finally, we
present a retuned baryonic physics model, based on the IllustrisTNG model,
which can be used for full-physics simulations within large cosmological
volumes.
