Ongoing and upcoming cosmological surveys will significantly improve our
ability to probe the equation of state of dark energy, $w_{\rm DE}$, and the
phenomenology of Large Scale Structure. They will allow us to constrain
deviations from the $\Lambda$CDM predictions for the relations between the
matter density contrast and the weak lensing and the Newtonian potential,
described by the functions $\Sigma$ and $\mu$, respectively. In this work, we
derive the theoretical prior for the joint covariance of $w_{\rm DE}$, $\Sigma$
and $\mu$, expected in general scalar-tensor theories with second order
equations of motion (Horndeski gravity), focusing on their time-dependence at
certain representative scales. We employ Monte-Carlo methods to generate large
ensembles of statistically independent Horndeski models, focusing on those that
are physically viable and in broad agreement with local tests of gravity, the
observed cosmic expansion history and the measurement of the speed of
gravitational waves from a binary neutron star merger. We identify several
interesting features and trends in the distribution functions of $w_{\rm DE}$,
$\Sigma$ and $\mu$, as well as in their covariances; we confirm the high degree
of correlation between $\Sigma$ and $\mu$ in scalar-tensor theories. The
derived prior covariance matrices will allow us to reconstruct jointly $w_{\rm
DE}$, $\Sigma$ and $\mu$ in a non-parametric way.