Anisotropic massive Brans-Dicke gravity extension of the standard
  $\Lambda$CDM model
  We present an explicit detailed theoretical and observational investigation
of an anisotropic massive Brans-Dicke (BD) gravity extension of the standard
$\Lambda$CDM model, wherein the extension is characterized by two additional
degrees of freedom; the BD parameter, $\omega$, and the present day density
parameter corresponding to the shear scalar, $\Omega_{\sigma^2,0}$. The BD
parameter, determining the deviation from general relativity (GR), by alone
characterizes both the dynamics of the effective dark energy (DE) and the
redshift dependence of the shear scalar. These two affect each other depending
on $\omega$, namely, the shear scalar contributes to the dynamics of the
effective DE, and its anisotropic stress --which does not exist in scalar field
models of DE within GR-- controls the dynamics of the shear scalar deviating
from the usual $\propto(1+z)^6$ form in GR. We mainly confine the current work
to non-negative $\omega$ values as it is the right sign --theoretically and
observationally-- for investigating the model as a correction to the
$\Lambda$CDM. By considering the current cosmological observations, we find
that $\omega\gtrsim 250$, $\Omega_{\sigma^2,0}\lesssim 10^{-23}$ and the
contribution of the anisotropy of the effective DE to this value is
insignificant. We conclude that the simplest anisotropic massive BD gravity
extension of the standard $\Lambda$CDM model exhibits no significant deviations
from it all the way to the Big Bang Nucleosynthesis. We also point out the
interesting features of the model in the case of negative $\omega$ values; for
instance, the constraints on $\Omega_{\sigma^2,0}$ could be relaxed
considerably, the values of $\omega\sim-1$ (relevant to string theories)
predict dramatically different dynamics for the expansion anisotropy.
