LAUSR

We provide an end-to-end exploration of a distinct modified gravitational theory in Jordan-Brans-Dicke (JBD) gravity, from an analytical and numerical description of the background expansion and linear perturbations, to the nonlinear regime captured with a hybrid suite of $N$-body simulations, to the parameter constraints from existing cosmological probes. The nonlinear corrections to the matter power spectrum due to baryons, massive neutrinos, and modified gravity are simultaneously modeled and propagated in the cosmological analysis for the first time. In the combined analysis of the Planck CMB temperature, polarization, and lensing reconstruction, Pantheon supernova distances, BOSS measurements of BAO distances, the Alcock-Paczynski effect, and the growth rate, along with the joint ($3\times2$pt) dataset of cosmic shear, galaxy-galaxy lensing, and overlapping redshift-space galaxy clustering from KiDS and 2dFLenS, we constrain the JBD coupling constant, $\omega_{\rm BD}>1540$ (95% CL), the effective gravitational constant, $G_{\rm matter}/G=0.997\pm0.029$, the sum of neutrino masses, $\sum m_{\nu}<0.12$ eV (95% CL), and the baryonic feedback amplitude, $B<2.8$ (95% CL), all in agreement with the standard model expectation. We show that the uncertainty in the gravitational theory alleviates the tension between KiDS$\times$2dFLenS and Planck to below $1\sigma$ and the tension in the Hubble constant between Planck and the direct measurement of Riess et al. (2019) down to ~$3\sigma$; however, we find no substantial model selection preference for JBD gravity relative to $\Lambda$CDM. We further show that the neutrino mass bound degrades by up to a factor of $3$ as the $\omega_{\rm BD}$ parameterization becomes more restrictive, and that a positive shift in $G_{\rm matter}/G$ suppresses the CMB damping tail in a way that might complicate future inferences of small-scale physics. (Abridged)