Anisotropic massive Brans–Dicke gravity extension of the standard $\Lambda$CDM model

Abstract

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.

1 Introduction

The base $\Lambda$-cold dark matter ($\Lambda$CDM) model, relying on the inflationary paradigm [1,2,3,4,5,6,7,8], is the simplest and most successful cosmological model to describe the dynamics and the large scale structure in agreement with the most of the currently available observational data [9,10,11,12]. It is today credited as the standard cosmological model, yet it is probably not where the story has concluded but the hardest part has just begun. It suffers from severe theoretical issues relating to the cosmological constant $\Lambda$ being responsible for the late time acceleration of the Universe [13,14,15,16,17,18,19] and, on the observational side, from tensions of various degrees of significance between some existing data sets [20,21,22,23,24,25,26,27,28]. Besides, based on the most minimal a priori assumptions, model independent reconstructions of the evolution of the dark energy (DE) equation of state (EoS) parameter w [29,30,31] and also model independent diagnoses [32, 33] exhibit a dynamical behaviour of w(z).

Nevertheless, even small deviations from/corrections to the $\Lambda$CDM mostly imply/require profound, and sometimes highly non-trivial, modifications to the fundamental theories of physics [34]. Indeed, we still do not have a promising and concrete fundamental theory leading to DE models that are more general than $\Lambda$ and also can account for the small, but significant, deviations from the $\Lambda$CDM model as persistently suggested by the high precision data. Though, depending on the characteristics of the DE models favored by observations, we can decide whether it is more natural to consider DE as an actual physical source [15, 35] or as an effective source originated from a modification [36,37,38,39,40] to the standard theory of gravity, i.e., general relativity (GR). For instance, eliminating $\Lambda$ by detecting $w >-1$ would not be by itself illuminating to the nature of DE, but any detection of $w < -1$ would be very illuminating to the nature of gravitation. For a perfect barotropic fluid, the adiabatic sound speed $c^2_a$ is the physical propagation speed of perturbations, and therefore $w<-1$ (viz., phantom [41] or quintom [42, 43] DE models) typically accompanying by $c^2_a<0$, implies the instability of perturbations (and/or ghost instabilities, see, e.g., [44]), whereas, as shown in [45, 46], scalar-tensor theories of gravity, such as Brans–Dicke (BD) theory, can lead to an effective source with $w<-1$ that does not correspond to the change of the sign of $c^2_a$, and hence perturbations can still be stable. Similarly, any detection of DE with negative energy density would also be hugely informative, such that this for an actual physical source is of course physically ill, whereas it can be an effective source so that negative energy density does not lead to any pathology since this is not the true energy density. On the observational side, the constraints on the EoS parameter of DE persistently indicate that $w\approx -1$, and so do not exclude $w < -1$ [12] and it has recently been shown in [21, 32, 47] that DE models with energy densities passing below zero at high redshifts fit the data better and can address the tensions relevant to Lyman-$\alpha$ forest measurements [20]. Such DE sources are indeed possible, in general, in modified theories of gravity with an effective gravitational coupling strength weaker in the past, e.g., in scalar-tensor theories, if we collect all modifications to the usual Einstein field equations to define an effective DE. Consequently, it can be argued that any detection of a deviation from $\Lambda$CDM model implies that the late time acceleration is not driven by $\Lambda$¹ and of DE yielding $w< -1$ and/or $\rho <0$ implies that gravity is not minimally coupled (which eliminate the actual perfect-fluid models of DE and stands as a strong sign in favor of the effective DE models from modified gravity).

Generically, the effective DE models from modified gravity induce non-zero anisotropic stresses,² which can be realized if we relax the isotropic space assumption of the $\Lambda$CDM model along with the modification to GR. Of course, we are not able to observe anisotropic stresses directly, yet they can reveal themselves through their effect on the evolution of the expansion anisotropy as well as on the average expansion rate of the Universe. Then, a natural question is whether the observations allow or suggest an anisotropic space, unless otherwise anisotropic stresses would be irrelevant. Indeed, there are some clues for questioning the spatially maximally symmetric Universe assumption, i.e., Robertson–Walker (RW) background, of the $\Lambda$CDM model. This has been mainly motivated by hints of anomalies in the CMB distribution first observed on the full sky by the WMAP experiment [49,50,51,52] and so remained in Planck experiment [53,54,55,56]. So far, the local deviations from the statistically highly isotropic Gaussianity of the CMB in some directions (the so called cold spots) could not have been excluded at high confidence levels [52, 53, 57, 58]. Furthermore, it has been shown that the CMB angular power spectrum has a quadrupole power lower than expected from the best-fit $\Lambda$CDM model [59, 60]. Several explanations for this anomaly have been proposed [61,62,63,64,65,66,67,68,69,70] including the anisotropic expansion of the Universe. On the other hand, inflation (canonical) isotropizes the Universe very efficiently [71, 72], leaving a residual anisotropy that is negligible for any practical application in the observable Universe. This could be irrelevant if an anisotropic expansion is developed only well after the matter-radiation decoupling, for example during the domination of DE, say, by means of its anisotropic pressure acting as a late-time source of not insignificant anisotropy [69, 73,74,75,76,77,78,79,80,81,82,83] (see also, e.g., [84,85,86,87] for constraint studies on the anisotropy of DE). Indeed, the CMB provides very tight constraints on the anisotropy at the time of recombination [88,89,90] of the order of the quadrupole temperature, i.e., $(\Delta T/T)_{\ell =2} \sim 10^{-5}$. And, in the simplest anisotropic generalization of the standard $\Lambda$CDM (viz., replacing the spatially flat RW metric by Bianchi type I metric), the energy density corresponding to the expansion anisotropy scales as the inverse of the square of the comoving volume, $\rho _{\sigma ^2}\propto S^{-6}$, which implies an isotropization of the expansion from the recombination up to the present, leading to the typically derived upper bounds on the corresponding density parameter today as $\Omega _{\sigma ^2,0}\sim 10^{-20}$. However, this is true if the anisotropic expansion is not generated by any anisotropic source, say, by an anisotropic DE (effective or actual), arising after decoupling [69, 73,74,75,76,77,78,79,80,81,82,83]. Nevertheless, almost all of these strong constraints on the expansion anisotropy in the late Universe, $\Omega _{\sigma ^2,0}\lesssim 10^{-22}$, assume the non-existence of anisotropic sources (actual or effective) in the late universe, and are derived in fact from the contribution of the expansion anisotropy on the average expansion rate of the Universe [91,92,93,94,95]. On the other hand, direct observational constraints, e.g., from SNIa data, on the expansion anisotropy of the later Universe are much weaker, namely, on the order of $\Omega _{\sigma ^2,0}\lesssim 10^{-4}$ for $z\sim 0$ [96, 97].

It is in fact always possible to imitate any modification to the Einstein field equations as GR with an arbitrary mixture of scalars, vectors and tensors [98, 99]. Therefore, one may think of abandoning any attempt to distinguish between the actual physical DE and the effective DE from modified gravities. However, a modified theory of gravity, e.g., the Brans–Dicke theory, can modify both the average expansion rate of the Universe and the evolution of the expansion anisotropy in a distinctive way depending on its free parameter characterizing the model, e.g., the Brans–Dicke parameter $\omega$. So, of course, owing to Occam’s razor, looking for these two concurrent modifications predicted by a specific modified theory of gravity and testing them against the observational data can provide us a strong reason for favoring or disfavoring the modified theory of gravity under consideration over GR. In this paper, for instance, the effective EoS parameter corresponding to the expansion anisotropy, the present day EoS of the effective DE and the redshift at which it evolves into the phantom region range in specific intervals as $1< w_{\sigma ^2}\le \frac{5}{3}$, $-1.14\lesssim w<-1$ and $0.50 \le z_\mathrm{PDL}< 0.65$ (see [100] for a similar result) for positive values of the BD parameter $\omega$. Indeed, revealing the origin of the late time acceleration of the Universe considering such distinctive features of the DE models from modified theories of gravity is among the scientific themes underlying some upcoming experiments (e.g., [101]).

Finally, we have also independent motivations for modifying GR by fundamental theoretical physics. Namely, almost all of the attempts to quantize GR introduce deviations from it in the form of extra degrees of freedom, higher powers of the curvature in the action, higher order derivatives in the field equations, or non-local terms. The low-energy limits of the string theories typically yield BD gravity with $\omega =-1$ and similarly d-brane models yield BD gravity with $\omega \sim -1$ [102,103,104,105]. Indeed, among all possible alternatives to GR, it is found natural by many to first consider scalar-tensor theories [106,107,108,109,110], which only add a (usually massive) scalar degree of freedom, the BD-like scalar field, to the two massless, spin-2 polarizations (gravitons) contained in the metric tensor. In relevance with our focus in this paper, the field equations in such models can be regarded as effective Einstein field equations and so the terms originating from the BD-like scalar field and its derivatives, moved to the right hand side, can always be interpreted as an effective fluid. Within this approach, the correspondence between general scalar-tensor theories and this effective fluid has been worked out explicitly first in [111] (see also [112] for the preceding work), and has shown that, in general, the corresponding effective fluid is imperfect. A recent work [113] extended and completed the correspondence showing how a symmetry of BD theory translates into a symmetry of this fluid such that this correspondence is valid for any spacetime geometries. Brans and Dicke’s 1961 theory [106, 110], motivated first by the implementation of Mach’s principle in gravity, today provides a prototype of scalar-tensor theories [107,108,109].

The original BD theory [106, 110] has only one additional constant parameter $\omega$, determining the deviations from GR, which is recovered in the limit of $|\omega | \rightarrow \infty$. Theoretically $\omega \ge -\frac{3}{2}$ is necessary to avoid the Jordan/scalar field ($\varphi$) from yielding negative energy density values in the Einstein frame that leads, e.g., the Minkowski vacuum to be unstable. Observations on a wide range of scales further constrain BD theory around GR. The tightest constraints, to date, are imposed by the observations of gravitational phenomena in the Solar System; specifically, $\omega \gtrsim 40000$ at the 2$\sigma$ confidence level (C.L.) from observations of radio signals from the Cassini spacecraft as it passes behind the Sun [114]. It can be vastly relaxed for the massive BD theory, i.e., when the Jordan field is accompanied by a potential $U(\varphi )$, as in this case one should consider the corresponding mass (M) of the Jordan field along with the BD parameter, viz., the $\{\omega ,M\}$ parameter region [115]. It is shown that $\omega =O(1)$ is allowed for the Solar System constraints provided that $M\gtrsim 2\times 10^{-17}\mathrm eV$ [115]. The constraints on M, however, are typically subject to the cosmological observations, since, in particular, it becomes cosmologically significant at substantially lower scales, viz., at the Hubble mass scale $M_{H_{0}}\simeq 10^{-33}\mathrm eV$, which leaves the constraint from the Cassini spacecraft still valid.

The cosmological constraints on $\omega$ are complementary to the local constraints as they probe different length and time scales, as well as different epochs of the Universe. The strongest cosmological constraint, which assumes $U(\varphi )=\mathrm const.$ and considers CMB data from the first Planck release along with the constraints from Big Bang Nucleosynthesis (BBN) light element abundances, gives $\omega >890$ at the 2$\sigma$ C.L. [116]. The constraints from Planck 2015 (2013) and a compilation of baryon acoustic oscillations (BAO) data are $\omega >333$ ($\omega >208$) at the 2$\sigma$ C.L., weakly dependent on the exponent $n>0$ for a power-law potential $U(\varphi )\propto \varphi ^n$ mimicking $\Lambda$ in the late universe [117, 118] (also see [119] for a recent work considering a more complicated potential). The BBN leads to $\omega \gtrsim 277$ (for a universe containing dust, radiation and conventional vacuum energy), which assumes power-law solutions and, by using general solutions, can be somewhat altered depending on the behaviour of the Jordan field in the early universe [120]. It is forecast that future cosmological experiments will be able to constrain $\omega$ at a level comparable to the local tests [121,122,123].

The typical mass scale of $M\sim 10^{-33}\mathrm eV$ imposed by cosmological observations is so small that the local constraint $\omega \gtrsim 40000$ from the Cassini spacecraft must be satisfied at all times in all parts of the Universe, and leads to the conclusion that the massive BD extension of the standard $\Lambda$CDM model must be phenomenologically very similar to it throughout most of the history of the Universe, implying that the cosmological features of BD theory would be observationally irrelevant. On the other hand, it is still worth studying in detail and understanding the cosmological constructions within massive BD theory by considering the cosmological constraints on $\omega$ only, as BD theory is representative of some more general gravity theories that can evade the local constraints (see [37] and references therein for further reading). For instance, straightforwardly, it is possible to consider a further extension of BD theory by allowing the BD parameter to be some functions of the Jordan field, $\omega =\omega (\varphi )$, in addition to the presence of potential $U(\varphi )$ and this leads to the possibility of having an $\omega$ small enough to be significant in the past, while being large enough today to be consistent with the tight constraints imposed by observations of gravitational phenomena in the Solar System. Alternatively, for instance, in the Horndeski theory (the most general scalar–tensor theory having second-order field equations in four dimensions) [124, 125], while on cosmological scales BD theory emerges as an approximation, the derivative self-interactions of the Horndeski scalar can be large enough on small scales (higher curvature than cosmological environments) to screen the scalar resulting at the same time to the recovery GR. For example, it was proposed in [126] that the so called Galileon theories (a close cousin of BD theory) could also explain the late time acceleration of the Universe while evading solar system constraints. A different kind of example is that, in the presence of extra dimensions, it is possible to make BD theory (massive or massless) consistent with gravitational tests (including solar system tests) for $|\omega |=O(1)$ [127]. See also [128] for a recent overview on cosmological probes of gravity for testing deviations from GR in relevance with such modified gravity theories.

Motivated by the discussions above, we present a detailed theoretical and observational investigation of the anisotropic (LRS Bianchi type I metric) massive Brans–Dicke gravity extension of the standard $\Lambda$CDM model, which is characterized by two additional free parameters, namely, the BD parameter $\omega$ and the corresponding present day density parameter of the expansion anisotropy $\Omega _{\sigma ^2,0}$. The role of the cosmological constant is taken over by the Jordan field potential of the form $U(\varphi )\propto \varphi ^2$ and the physical ingredient of the Universe is considered as in the standard $\Lambda$CDM. We confine the present work, basically, to the non-negative values of the BD parameter since, here, we mainly consider the extension as a correction to the standard $\Lambda$CDM.

2 Anisotropic massive Brans–Dicke gravity extension of the $\Lambda$CDM model

We consider the BD action [106, 110] written in the Jordan frame in the following form:

$$\begin{aligned} \begin{aligned} S_\mathrm{JBD}=&\int \mathrm{d}^4 x \sqrt{-g}\bigg [\frac{\varphi ^2}{8}R-\omega \bigg (\frac{1}{2}\nabla _{\mu }\varphi \nabla ^{\mu }\varphi +\frac{1}{2}M^2\varphi ^2\bigg )\bigg ]\\&+S_\mathrm{Matter}, \end{aligned} \end{aligned}$$

(1)

where $\varphi =\varphi (t)$ is the Jordan scalar field (function of cosmic time t only) and $\omega =\mathrm const.$ is the Brans–Dicke parameter, R is the Ricci scalar, g is the determinant of the metric $g_{\mu \nu }$, and $S_\mathrm{Matter}$ is the matter action, which is independent of $\varphi$ so that the weak equivalence principle is satisfied. It is clear from the way of writing the action that the term $M^2$ stands as the bare mass-squared of the Jordan field.³ We assume $M^2=\mathrm{const.}$ so that, as can also be seen from the action, it stands like a cosmological constant as $2\omega M^2\equiv \Lambda$ and thereby can drive accelerated expansion. Hence, switching to massive BD from GR with a positive cosmological constant provides us with an opportunity to construct $\Lambda$CDM-type cosmologies, such that the mass of the Jordan field alone can play the role of positive cosmological constant like in the standard $\Lambda$CDM cosmology provided that $2\omega M^2\equiv \Lambda >0$, and the Jordan field $\varphi$ varying slowly enough on the top of this can account for small deviations from the standard $\Lambda$CDM model in a particular way, which in turn may lead to an improved fit to the observational data w.r.t. the standard $\Lambda$CDM model. In line with that, we intend to study the BD extension of the standard $\Lambda$CDM model as a correction, and therefore we demand the term $2\omega M^2$ to be positive definite, which requires $\omega >0$ as long as we keep $M^2>0$ to avoid the Jordan field from having an imaginary mass.⁴ Hence, in this study, unless otherwise is mentioned, we carry out our investigations by assuming $\omega >0$, which is already stronger than the assumption $\omega \ge -3/2$ to avoid the Jordan field from yielding negative energy density values in the Einstein frame that leads, for instance the Minkowski vacuum to be unstable.

We consider the simplest anisotropic generalization of the spatially flat and homogeneous spacetime, i.e., locally rotationally symmetric (LRS) Bianchi type I metric, which can be written as follows;

$$\begin{aligned} \mathrm{d}s^2 = -\mathrm{d} t^2+S^2\left[ e^{4\beta }\mathrm{d}x^2+ e^{-2\beta }(\mathrm{d}y^2+\mathrm{d}z^2)\right] , \end{aligned}$$

(2)

where $S=S(t)$ is the mean scale factor. Here, the exponent $\beta =\beta (t)$ satisfies the relation $\dot{\beta }^2=\frac{1}{6}\sigma ^2$, where $\sigma ^2=\sigma _{ij}\sigma ^{ij}$ and $\sigma ^{ij}$ are shear scalar and tensor, respectively. The scalar curvature for this metric (2) may be written in terms of the mean scale factor and shear scalar as $R=-6\left( \frac{\ddot{S}}{S}+\frac{\dot{S}^2}{S^2}\right) -\sigma ^2$. Throughout the paper a dot denotes derivative w.r.t. cosmic time t.

We consider all types of matter distribution (namely, the usual cosmological sources such as radiation, baryons, etc.) as isotropic perfect fluids, which are described by the following energy-momentum tensor (EMT)

$$\begin{aligned} {T_{\mu }}^{\nu }={\text {diag}}[-\rho ,p,p,p], \end{aligned}$$

(3)

where $\rho$ and p are the energy density and pressure, respectively.

The field equations for the action (1) within the framework of the metric (2) in the presence of (3) read:

$$\begin{aligned}&3\frac{\dot{S}^2}{S^2}-\frac{1}{2}\sigma ^2-2\omega \frac{\dot{\varphi }^2}{\varphi ^2}+6\frac{\dot{S}}{S} \frac{\dot{\varphi }}{\varphi }-2\omega M^2=\frac{4}{\varphi ^2}\rho , \end{aligned}$$

(4)

$$\begin{aligned}&2\frac{\ddot{S}}{S}+\frac{\dot{S}^2}{S^2}+\frac{1}{2}\sigma ^2-\sqrt{\frac{2}{3}}\left[ \dot{\sigma }+\left( 3\frac{\dot{S}}{S}+2\frac{\dot{\varphi }}{\varphi } \right) \sigma \right] \nonumber \\&\quad +(2 \omega +2) \frac{\dot{\varphi }^2}{\varphi ^2} +2\frac{\ddot{\varphi }}{\varphi }+4\frac{\dot{S}}{S}\frac{\dot{\varphi }}{\varphi }-2 \omega M^2=-\frac{4}{\varphi ^2}p, \end{aligned}$$

(5)

$$\begin{aligned}&2\frac{\ddot{S}}{S}+\frac{\dot{S}^2}{S^2}+\frac{1}{2}\sigma ^2+\sqrt{\frac{1}{6}}\left[ \dot{\sigma }+\left( 3\frac{\dot{S}}{S}+2\frac{\dot{\varphi }}{\varphi } \right) \sigma \right] \nonumber \\&\quad +(2 \omega +2) \frac{\dot{\varphi }^2}{\varphi ^2}+2\frac{\ddot{\varphi }}{\varphi }+4\frac{\dot{S}}{S} \frac{\dot{\varphi }}{\varphi }-2 \omega M^2=-\frac{4}{\varphi ^2}p, \end{aligned}$$

(6)

for which the latter two correspond to the pressure equations along the x-axis and the y- and z-axes, respectively. The Klein–Gordon equation for the Jordan field reads

$$\begin{aligned} \frac{\ddot{\varphi }}{\varphi }+3\frac{\dot{S}}{S}\frac{\dot{\varphi }}{\varphi } -\frac{3}{2\omega }\left( \frac{\ddot{S}}{S}+\frac{\dot{S}^2}{S^2}\right) -\frac{1}{2\omega }\frac{\sigma ^2}{2}+M^2=0. \end{aligned}$$

(7)

We consider the standard cosmological fluids, namely, pressureless fluid ($p_{\mathrm{m}}= 0$) and radiation/relativistic fluids ($p_{\mathrm{r}} = \rho _\mathrm{r}/3$), so that, here, $\rho =\rho _{\mathrm{m}}+\rho _{\mathrm{r}}$ and $p=p_{\mathrm{r}}$. If the Jordan field is approximately constant then the first Friedmann Eq. (4) becomes $3\frac{\dot{S}^2}{S^2}\approx \frac{4}{\varphi ^2}\rho +\frac{\sigma ^2}{2}+2\omega M^2$ (as like in GR), where the effective cosmological gravitational coupling strength is then given by $G=\frac{\varphi _0^2}{\varphi ^2}\, G_0$, where $G=\frac{1}{2\pi \varphi ^2}$ and $G_0=\frac{1}{2\pi \varphi ^2_0}$ with zero subscript denoting today’s value (throughout the paper).⁵

We consider the conventional representation of the true equation for scalar-tensor gravity interacting with matter in the Einsteinian form with the constant $G_0=G(t_0)$

$$\begin{aligned} R_{\mu \nu }-\frac{1}{2}R g_{\mu \nu }=8\pi G_0 (T_{\mu \nu ,\mathrm{m}}+T_{\mu \nu ,\mathrm{DE}}), \end{aligned}$$

(8)

$8\pi G_0=4/\varphi ^2_0$ following Refs. [45, 46], where it is discussed that modified gravity theories can be recast in the standard GR form such that (8) where all the new geometrical terms are grouped (on the r.h.s.) to form an effective DE contribution denoted as $T_{\mu \nu ,\mathrm{DE}}$ (for a detailed discussion, see introduction in Ref. [46]). We note that the shear scalar is kept on the left hand side as it is a part of the Einstein tensor, namely, the metric itself only. Accordingly, we define the effective DE in the following way:

$$\begin{aligned}&\begin{aligned} 3\frac{\dot{S}^2}{S^2}-\frac{1}{2}\sigma ^2=\frac{4}{\varphi ^2_0}\,(\rho +\rho _\mathrm{DE}), \end{aligned} \end{aligned}$$

(9)

$$\begin{aligned}&\begin{aligned} 2\frac{\ddot{S}}{S}+\frac{\dot{S}^2}{S^2}+\frac{1}{2}\sigma ^2-\sqrt{\frac{2}{3}}\left( \dot{\sigma }+3\frac{\dot{S}}{S}\sigma \right) = -\frac{4}{\varphi ^2_0}\,(p+p_{\mathrm{DE},x}), \end{aligned} \nonumber \\ \end{aligned}$$

(10)

$$\begin{aligned}&\begin{aligned} 2\frac{\ddot{S}}{S}+\frac{\dot{S}^2}{S^2}+\frac{1}{2}\sigma ^2+\sqrt{\frac{1}{6}}\left( \dot{\sigma }+3\frac{\dot{S}}{S}\sigma \right) = -\frac{4}{\varphi ^2_0}\,(p+p_{\mathrm{DE},y}), \end{aligned} \nonumber \\ \end{aligned}$$

(11)

where $\rho _{\mathrm{DE}}$ is the energy density and $p_{\mathrm{DE},x}$ and $p_{\mathrm{DE},y}$ are the principal pressures of the effective DE along the x-axis and the y- and z-axes, respectively, and which read

$$\begin{aligned} \rho _{\mathrm{DE}}= & {} \frac{\varphi _0^2}{4}\bigg [\left( \frac{4}{\varphi ^2}-\frac{4}{\varphi _0^2} \right) \rho + 2\omega \frac{\dot{\varphi }^2}{\varphi ^2}-6\frac{\dot{S}}{S}\frac{\dot{\varphi }}{\varphi }+2\omega M^2\bigg ], \nonumber \\ \end{aligned}$$

(12)

$$\begin{aligned} p_{\mathrm{DE},x}= & {} \frac{\varphi _0^2}{4}\bigg [\left( \frac{4}{\varphi ^2}-\frac{4}{\varphi _0^2} \right) p-2\sqrt{\frac{2}{3}}\frac{\dot{\varphi }}{\varphi }\sigma \nonumber \\&+(2 \omega +2) \frac{\dot{\varphi }^2}{\varphi ^2} +2\frac{\ddot{\varphi }}{\varphi }+4\frac{\dot{S}}{S}\frac{\dot{\varphi }}{\varphi }-2\omega M^2\bigg ],\end{aligned}$$

(13)

$$\begin{aligned} p_{\mathrm{DE},y}= & {} \frac{\varphi _0^2}{4}\bigg [\left( \frac{4}{\varphi ^2}-\frac{4}{\varphi _0^2} \right) p+\sqrt{\frac{2}{3}}\frac{\dot{\varphi }}{\varphi }\sigma \nonumber \\&+(2 \omega +2) \frac{\dot{\varphi }^2}{\varphi ^2} +2\frac{\ddot{\varphi }}{\varphi }+4\frac{\dot{S}}{S}\frac{\dot{\varphi }}{\varphi }-2\omega M^2\bigg ]. \end{aligned}$$

(14)

We note that the Jordan field and shear scalar are coupled and they together lead to an anisotropy in the pressure of the effective DE, which may be represented by

$$\begin{aligned} \Delta p_{\mathrm{DE}}=p_{\mathrm{DE},y}-p_{\mathrm{DE},x}=\sqrt{6}\frac{\varphi _0^2}{4}\frac{\dot{\varphi }}{\varphi }\sigma . \end{aligned}$$

(15)

This, in turn, leads to a modification in the evolution of the expansion anisotropy provided that $\varphi$ is dynamical and the space is not exactly isotropic as would be expected from a realistic cosmological model (see, e.g., [132]). We can restate this equation in terms of the cosmological gravitation coupling strength as

$$\begin{aligned} \Delta p_{\mathrm{DE}}=-\sqrt{\frac{3}{2}}\frac{1}{8\pi G_0}\,\frac{\dot{G}}{G}\,\sigma . \end{aligned}$$

(16)

This effective anisotropic pressure (in the presence of shear) induced by the Jordan field in the Jordan frame is absent in the Einstein frame (see, e.g., Ref. [132]).

3 Exact solution compatible with standard $\Lambda$CDM model

We follow the method, given in [100], for obtaining exact cosmological solutions of a scalar-tensor gravity theory compatible with the $\Lambda$CDM model, albeit we consider LRS Bianchi type I spacetime rather than spatially flat Robertson–Walker spacetime.⁶$^{,}$⁷ To do so, we first re-express the set of differential equations given by Eqs. (4)–(7) using $\frac{\mathrm{d}z}{\mathrm{d} t}=-H(1+z)$, where $H=\dot{S}/S$ is the average Hubble parameter and z is the cosmic redshift we define in terms of S as $z=-1+\frac{1}{S}$. Accordingly, we re-express the energy density equation (4) as

$$\begin{aligned}&3H^2\bigg [1-\frac{2\omega }{3}(1+z)^2\frac{{\varphi '}^2}{\varphi ^2}-2(1+z)\frac{{\varphi '}}{\varphi }\bigg ]\nonumber \\&\quad -\frac{\sigma ^2}{2}-2\omega M^2=\frac{4}{\varphi ^2}\rho , \end{aligned}$$

(17)

and obtain the pressure equation as

$$\begin{aligned}&-\frac{\mathrm{d}H^2}{\mathrm{d}z}\left[ 1+z-(1+z)^2\frac{\varphi '}{\varphi }\right] \nonumber \\&\quad +3H^2\bigg \{1+(1+z)^2\bigg [\frac{2}{3}\frac{\varphi ''}{\varphi }+\frac{2}{3}(1+\omega )\frac{{\varphi '}^2}{\varphi ^2}\bigg ]\nonumber \\&\quad -\frac{2}{3}(1+z)\frac{\varphi '}{\varphi }\bigg \}+\frac{\sigma ^2}{2}-2\omega M^2=-\frac{4}{\varphi ^2}p, \end{aligned}$$

(18)

using Eqs. (5) and (6), the shear propagation equation as

$$\begin{aligned} \frac{\sigma '}{\sigma }-\frac{3}{1+z}+2\frac{\varphi '}{\varphi }=0, \end{aligned}$$

(19)

by subtracting Eq. (5) from Eq. (6), and finally re-express the scalar field equation (7) as

$$\begin{aligned}&-\frac{\mathrm{d}H^2}{\mathrm{d}z}\left[ \omega (1+z)^2\frac{\varphi '}{\varphi }+\frac{3}{2}(1+z)\right] \nonumber \\&\quad +H^2\left[ -2\omega \frac{\varphi ''}{\varphi }(1+z)^2+4\omega (1+z)\frac{\varphi '}{\varphi }+6\right] \nonumber \\&\quad +\frac{\sigma ^2}{2}-2\omega M^2=0, \end{aligned}$$

(20)

where $'$ denotes derivative with respect to redshift ($\mathrm{d}/\mathrm{d}z$).

We note that Eqs. (18) and (20) have exactly the same mathematical form, that is, a first order linear differential equation in $H^2$, such as

$$\begin{aligned} A_i(z)\frac{\mathrm{d}H^2}{\mathrm{d}z}+B_i(z)H^2+C_i(z)+D_i=0,\,\, i=1,2, \end{aligned}$$

(21)

provided that we set $p=p_{\mathrm{m}}=0$, viz., the universe is filled with only pressureless matter. We next note that constants $D_1=D_2=-2\omega M^2$, and that the shear scalar $\sigma ^2$ is the term that differs in our system of equations from the one given in [100]. Fortunately, it contributes to (18) and (20) in the same way, namely, as $C_1(z)=C_2(z)=\frac{\sigma ^2}{2}$. In accordance with these points, we assume $A_1(z)=A_2(z)$ in (18) and (20) and then solve for the rate of change of Jordan field in z as

$$\begin{aligned} \frac{\varphi '}{\varphi }=-\frac{1}{2(1+\omega )}\frac{1}{1+z}, \end{aligned}$$

(22)

which in turn renders the coefficients of $H^2$ identical, i.e, $B_1(z)=B_2(z)$, as well. Thus, integrating (22), it turns out that the solution of the Jordan field is

$$\begin{aligned} \varphi ^2=\varphi _0^2(1+z)^{-\frac{1}{1+\omega }}, \end{aligned}$$

(23)

where $\varphi _0=\varphi (z=0)$ is the present time value of the Jordan field, and which in turn gives

$$\begin{aligned} G=G_{0} (1+z)^{\frac{1}{1+\omega }}. \end{aligned}$$

(24)

The integration of the shear propagation equation given in Eq. (19) gives

$$\begin{aligned} \sigma ^2=\sigma _0^2(1+z)^6\left( \frac{\varphi }{\varphi _0}\right) ^{-4}, \end{aligned}$$

(25)

where $\sigma ^2_0=\sigma ^2(z=0)$ is the present time value of the shear scalar, leading, together with (23), to

$$\begin{aligned} \sigma ^2=\sigma _0^2(1+z)^{6+\frac{2}{1+\omega }}, \end{aligned}$$

(26)

whereas it is $\sigma _\mathrm{GR}^2 \propto (1+z)^6$ in GR (corresponding to the $| \omega | \rightarrow \infty$ limit of the BD gravity), and thereby the energy density corresponding to the shear scalar can be defined as

$$\begin{aligned} \rho _{\sigma ^2}\equiv \frac{\sigma ^2}{2} \frac{\varphi _0^2}{4}=\frac{\sigma ^2_0}{2} \frac{\varphi _0^2}{4}(1+z)^{6+\frac{2}{1+\omega }}, \end{aligned}$$

(27)

by considering the present time value of the Jordan field, i.e., of the cosmological gravitational coupling strength. Because BD theory [110] satisfies the local conservation of EMT, the energy density of the pressureless matter can immediately be written as

$$\begin{aligned} \rho _{\mathrm{m}}=\rho _{\mathrm{m,0}}(1+z)^3, \end{aligned}$$

(28)

which is exactly the same as the one for the pressureless matter in the standard $\Lambda$CDM model. Finally using $\varphi ^2$, $\sigma ^2$ and $\rho _{\mathrm{m}}$ from Eqs. (23), (26) and (28), respectively, in Eq. (17) we reach the following modified anisotropic Friedmann equation for BD theory in this solution

$$\begin{aligned} \begin{aligned} H^2=&\frac{4}{3}\frac{\gamma }{\varphi _0^2}\left[ \rho _{\mathrm{M}}+\rho _{\mathrm{m},0}(1+z)^{3+\frac{1}{1+\omega }}+\rho _{\sigma ^2,0}(1+z)^{6+\frac{2}{1+\omega }}\right] , \end{aligned} \end{aligned}$$

(29)

where

$$\begin{aligned} \gamma =\frac{6(1+\omega )^2}{(3\omega +4)(2\omega +3)}, \end{aligned}$$

(30)

and the energy densities corresponding to the mass of the Jordan field and to the expansion anisotropy today read, respectively,

$$\begin{aligned} \begin{aligned} \rho _{\mathrm{M}}=2M^2\omega \frac{\varphi _0^2}{4}\quad \text {and} \quad \rho _{\sigma ^2,0}=\frac{\sigma _0^2}{2}\frac{\varphi _0^2}{4}. \end{aligned} \end{aligned}$$

(31)

One may check for consistency that we recover the H(z) of the standard $\Lambda$CDM model when we set $\omega \rightarrow \infty$ (viz., in the GR limit of BD gravity giving $\varphi \rightarrow \mathrm{constant}$) and $\rho _{\sigma ^2,0}=0$ (viz., isotropic expansion as in the RW spacetime metric). We also note that $\rho _{\mathrm{M}}$ can be regularized to give a finite positive value consistent with the observations by setting $M^2\omega =\mathrm{constant}$ so that $M^2\rightarrow 0$ as $\omega \rightarrow \infty$ such that the term $M^2\omega$ arising from their multiplication would always remain unaltered and finite.

3.1 Inclusion of radiation

The solution of the BD extension of standard $\Lambda$CDM model, for either isotropic or anisotropic spatially flat spacetimes, can be achieved by using the method given in Ref. [100] provided that $p=0$, viz., the Universe is filled with only pressureless matter. It covers also the epoch of the Universe when the contributions from the pressureless matter in (29), namely, $\frac{4}{\varphi ^2}\rho _{\mathrm{m}}\propto (1+z)^{3+\frac{1}{1+\omega }}$, is dominant over the effective cosmological constant $2\omega M^2$ and one may check that for this epoch, neglecting anisotropy, it reduces to the well known solution with $S\propto t^{\frac{2+2\omega }{4+3\omega }}$ and $\varphi ^2\propto t^{\frac{2}{4+3\omega }}$ (see, e.g., [120] and references therein). We could not obtain analytically an exact cosmological solution – other than the trivial solution with $\varphi =\mathrm const.$ – using the same method when we include relativistic source, $p_{\mathrm{r}}=\frac{\rho _\mathrm{r}}{3}$. Fortunately the effects of the relativistic source on the background dynamics can be neglected in the late Universe. Namely, because local energy conservation holds in BD gravity in the Jordan frame, we have $\rho _{\mathrm{r}}\propto (1+z)^4$ and $\rho _{\mathrm{m}}\propto (1+z)^3$, implying that the effects of the relativistic source on the dynamics of the Universe will be significant for large redshift values, in particular, for the redshift values larger than the matter-radiation equality redshift, $z_{\mathrm{eq}}=-1+\frac{\rho _{\mathrm{m,0}}}{\rho _{\mathrm{r,0}}}\sim 3380$ [12], as in the standard cosmology [138]. The attractor solution of the BD gravity for radiation domination is well known that it is exactly the GR solution, i.e., $\varphi = \mathrm{const.}$ [120, 135,136,137,138]. Indeed, in general, BD cosmologies have exact solutions which show that they are dominated by the Jordan field at early times and by the perfect fluid matter sources at late times, which here can be considered as the period all the way down from matter-radiation equality to the earlier times covering the times relevant to Big Bang Nucleosynthesis (BBN). Hence, for $z>z_{\mathrm{eq}}$, we have $\varphi = \varphi _1= \mathrm{const.}$, implying $G=G_1= \mathrm{const.}$ as long as anisotropy is negligible. Note that $\varphi _1$ (or $G_1$) will be different than $\varphi _0$ (or $G_0$) as a consequence of the fact that the Jordan field has been evolving approximately in accordance with (23) for $z<z_{\mathrm{eq}}$, i.e., between the matter-radiation equality and today. Accordingly, our model will more approach a general relativistic cosmology as the radiation dominates over pressureless matter as we go to higher redshifts, yet the most part of period of the Universe during which $z<z_{\mathrm{eq}}$ the Jordan field evolves approximately in accordance with (23) and thereby we can estimate that

$$\begin{aligned} \begin{aligned} \varphi _1\sim \varphi _0(1+z_{\mathrm{eq}})^{-\frac{1}{2+2\omega }}\quad&\text {at}\quad z\sim z_{\mathrm{eq}} \\ \quad \quad \quad \quad \text {as well as} \quad \varphi \sim \varphi _1\quad&\text {for}\quad z\gtrsim z_{\mathrm{eq}}, \end{aligned} \end{aligned}$$

(32)

which implies

$$\begin{aligned} G_1\sim G_0(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }} \quad \text {for}\quad z\gtrsim z_{\mathrm{eq}}, \end{aligned}$$

(33)

where $G_1$ is then approximately the cosmological gravitational coupling strength for $z>z_{\mathrm{eq}}$ (see [120] for further details). Thus, for $\omega >0$, the cosmological gravitational coupling strength will remain constant (provided that anisotropy is negligible) during the radiation dominated epoch as just like in the standard cosmology, but will be larger than its present time value, e.g., $G_1/G_0\sim 3380^{\frac{1}{1+\omega }}$. This implies an enhanced expansion rate of the Universe during this epoch, which in turn will affect, for instance, the primordial element abundances. We shall further discuss this point in Sect. 5. We next see from (19) that we approximately have $\sigma ^2\propto \sigma _0^2(1+z)^6$ like in GR in the presence of isotropic perfect fluid, since we have $\varphi \sim \mathrm{constant}$ when the radiation is dominant over pressureless matter and expansion anisotropy is still insignificant. As we go to even larger redshift values, say, $z\gg z_{\mathrm{BBN}}$, expansion anisotropy will dominate over the terms due to the radiation and hence we obtain once again $\rho _{\sigma ^2}\propto (1+z)^{6+\frac{2}{1+\omega }}$ (which can be obtained simply by setting $\rho _{\mathrm{M}}=0$ and $\rho _{\sigma ^2,0}=0$ in our solution).

In the light of the above discussion, although we do not have exact analytical solution when we include the radiation as well, we can find out the contribution from the radiation. First of all, it is obvious that during radiation domination the Jordan field is constant, then the inclusion of radiation to the model will in general slow-down the Jordan field, namely, the model in general will deviate less from GR. Of course, that effect will be insignificant for a viable and realistic cosmological model [namely, cosmological model satisfying, roughly, $\rho _{\sigma ^2,0}\ll \rho _{\mathrm{r,0}}\ll \rho _{\mathrm{m,0}}\sim \rho _{\mathrm{M}}$ for $\omega \gtrsim 10$ (see Sect. 3.3.1 for this choice on $\omega$)] for $z< z_{\mathrm{eq}}$. Thus, we can write

$$\begin{aligned} H^2\approx & {} \frac{4}{3}\frac{\gamma }{\varphi _0^2}\bigg [\rho _{\mathrm{M}}+\rho _{\mathrm{m},0}(1+z)^{3+\frac{1}{1+\omega }}+\rho _{\mathrm{r},0}(1+z)^{4+\frac{1}{1+\omega }}\nonumber \\&+\rho _{\sigma ^2,0}(1+z)^{6+\frac{2}{1+\omega }}\bigg ]\quad \text {for}\quad z<z_{\mathrm{eq}}, \end{aligned}$$

(34)

which can be safely used all the way to the recombination redshift (CMB release), such that it is well known that the Universe should have transited from radiation to matter dominated era when $z\sim 3380$, and the recombination that leads to photon decoupling should have taken place when $z\sim 1100$, at which the pressureless matter is still dominant over radiation. Therefore, we note that H(z) given in (34) can safely be used for constraining the model by using CMB radiation data as well.

For the times before the matter-radiation equality, setting $\varphi$ constant relying on its attractor behavior during the radiation-dominated phase, we can write

$$\begin{aligned} H^2\approx & {} \frac{4}{3}\frac{\gamma }{\varphi _1^2}\left[ \rho _{\mathrm{r},0}(1+z)^{4}+\rho _{\sigma ^2,1}(1+z)^{6}\right] \nonumber \\&\quad \text {for}\quad z_{\mathrm{eq}}<z<z_\mathrm{eq,\sigma ^2,r}. \end{aligned}$$

(35)

Note that, for consistency between (34) and (35), here in (35) we have $\varphi _1\sim \varphi _0(1+z_{\mathrm{eq}})^{-\frac{1}{2+2\omega }}$ in contrast to $\varphi _{0}$ in (34), in accordance with the above discussions leading to (32), and similarly $\rho _{\sigma ^2,1}\sim \rho _{\sigma ^2,0}(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }}$. Using (35), we estimate (at least roughly) the radiation-expansion anisotropy equality redshift as $z_\mathrm{eq,\sigma ^2,r}\sim -1+ \sqrt{\frac{\rho _{\mathrm{r,0}}}{\rho _{\sigma ^2,1}}}$, which can also be written as follows;

$$\begin{aligned} z_\mathrm{eq,\sigma ^2,r}\sim -1+ \sqrt{\frac{\rho _{\mathrm{r,0}}}{\rho _{\sigma ^2,0}}\left( \frac{\rho _{\mathrm{m,0}}}{\rho _{\mathrm{r,0}}}\right) ^{-\frac{1}{1+\omega }}}. \end{aligned}$$

(36)

Avoiding expansion anisotropy from spoiling the physical processes relevant to BBN will then roughly require $z_\mathrm{eq,\sigma ^2,r}>z_{\mathrm{BBN}}$. Typically, in a viable cosmological model, we expect the matter-radiation equality redshift to be $z_{\mathrm{eq}}=-1+\frac{\rho _{\mathrm{m,0}}}{\rho _{\mathrm{r,0}}}\sim 3380$, the BBN to take place at $z_{\mathrm{BBN}}\sim 3\times 10^{8}$ and $\Omega _{\mathrm{r,0}}\sim 10^{-4}$. Accordingly, taking $z_\mathrm{eq,\sigma ^2,r}> 3\times 10^{8}$, we obtain $\frac{\rho _{\sigma ^2, 0}}{\rho _\mathrm{r, 0}}=\frac{\Omega _{\sigma ^2, 0}}{\Omega _\mathrm{r, 0}}\lesssim 10^{-17-\frac{4}{1+\omega }}$, namely, the following stringent constraint on the density parameter of the expansion anisotropy today:

$$\begin{aligned} \Omega _{\sigma ^2,0}\lesssim 10^{-21-\frac{4}{1+\omega }}, \end{aligned}$$

(37)

which reduces to $\Omega _{\sigma ^2,0}\sim 10^{-21}$ in the GR limit ($\omega \rightarrow \infty$). It is noteworthy that, in the GR limit, at which anisotropy in the pressure also disappears, our estimation on the expansion anisotropy for today is already of the same order of magnitude with the ones given in [94] from a general test of isotropy using CMB temperature and polarization data from Planck and in [91] from primordial nucleosynthesis, both of them consider GR and only isotropic fluids. This also shows the reliability of our estimation. We note that, for the case $\omega >0$ we consider in this study, BD theory leads to more stringent constraints on the expansion anisotropy (e.g., $\Omega _{\sigma ^2,0}\lesssim 10^{-25}$ for $\omega =0$), while the case $\omega <-1$, which we do not consider in the present study, can relax it. For instance, Campanelli et al. [96] have put model-independent upper bounds on the expansion anisotropy, from type Ia SNe observations, in the late Universe, namely, for $z\lesssim 1.6$, as $\Omega _{\sigma ^2,0}\lesssim 10^{-4}$ (see also [97] for a recent work). We note that this value is achieved upon choosing $\omega =-\frac{21}{17}$ in (37), which can be done without worrying for theoretical issues associated with $\omega <-\frac{3}{2}$, but this will cost large deviations from GR (viz. rapidly varying strength of the gravitational coupling) as well as the standard $\Lambda$CDM. We shall further discuss possible modifications on the expansion anisotropy for different values/ranges of $\omega$ and their comparisons with the standard $\Lambda$CDM model in Sect. 3.3.2.

Finally, for the times when the expansion anisotropy dominates over radiation (and hence obviously on pressureless matter), we can write

$$\begin{aligned} H^2\propto (1+z)^{6+\frac{2}{1+\omega }}\quad \text {for}\quad z>z_\mathrm{eq,\sigma ^2,r}, \end{aligned}$$

(38)

from our solution (see, e.g., Refs. [132, 139,140,141] for anisotropy dominated universes in the context of scalar-tensor theories). This epoch however is out of interest in the present work.

3.2 Effective anisotropic dark energy

The discussions in the previous subsections reveal that the massive BD gravity modification on the gravity sector, manifests itself when radiation becomes subdominant, say, in the relatively late Universe, that is, switching from GR+$\Lambda$ to the massive BD introduces corrections on the redshift dependencies of both the average expansion rate and expansion anisotropy. Although we consider the presence of radiation when we constrain the model using the cosmological data, while we are discussing the effective DE in what follows, we shall neglect its presence due to the following reasons: (i) Its effect on the effective DE will obviously be negligible in the late universe. (ii) When it is dominant the Jordan field $\varphi$ freezes, implying the BD gravity mimics GR, and thereby the effective DE will mimic nothing but $\Lambda$ (not to mention that DE would be subdominant at this epoch). (iii) A practical reason, such that we don’t have explicit analytical solution when we include radiation into the model. On the other hand, although we expect anisotropy to be negligible (even with respect to radiation) in the late universe, we keep expansion anisotropy due to the following reasons: (i) To see its effect on the effective DE explicitly. (ii) It leads to an anisotropy in the pressure of the effective DE, which in turn modifies the redshift dependency of the expansion anisotropy w.r.t. GR ($|\omega |\rightarrow \infty$). (iii) In relevance with (ii), its modified redshift dependency w.r.t. the one in GR (provided that it can be detected) would be a strong signal in favor of modified gravities (such as BD theory), rather than the presence of a DE (such as $\Lambda$, scalar fields) ingredient of the Universe assuming that GR is the true theory of gravity.

We implement the solution given in Sect. 3 and accordingly we re-write (9) as follows:

$$\begin{aligned} \begin{aligned} 3H^2=\frac{4}{\varphi ^2_0}\,\left( \rho _{\mathrm{m}}+\rho _{\mathrm{DE}}+\rho _{\sigma ^2}\right) , \end{aligned} \end{aligned}$$

(39)

where the energy density of the effective DE, $\rho _{\mathrm{DE}}$, reads

$$\begin{aligned} \rho _{\mathrm{DE}}=\gamma \rho _{\mathrm{M}}+\rho _{\mathrm{m}}\left[ \gamma \left( \frac{\rho _{\mathrm{m}}}{\rho _{\mathrm{m,0}}}\right) ^{\frac{1}{3(1+\omega )}}-1\right] +\rho _{\sigma ^2}(\gamma -1),\qquad \end{aligned}$$

(40)

which gives

$$\begin{aligned} \begin{aligned} \rho _{\mathrm{DE,0}}=&\,\gamma \rho _{\mathrm{M}}+\rho _{\mathrm{m,0}}\left( \gamma -1\right) +\rho _{\sigma ^2,0}(\gamma -1) \end{aligned} \end{aligned}$$

(41)

for $z=0$. Dividing (39) by $3H_0^2$, we obtain

$$\begin{aligned} \begin{aligned} \frac{H^2}{H_0^2}=\Omega _{\mathrm{DE,0}}f(z)+\Omega _{\mathrm{m,0}}(1+z)^3+\Omega _{\sigma ^2,0}(1+z)^{6+\frac{2}{1+\omega }}, \end{aligned} \end{aligned}$$

(42)

where we define $f(z)=\frac{\rho _{\mathrm{DE}}}{\rho _{\mathrm{DE,0}}}$, reading

$$\begin{aligned} f(z)=\,1+\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}\left( \gamma -1\right) \left[ (1+z)^{6+\frac{2}{1+\omega }}-1\right] \nonumber \\ +\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}\left\{ (1+z)^3\left[ \gamma (1+z)^{\frac{1}{1+\omega }}-1\right] +1-\gamma \right\} , \end{aligned}$$

(43)

which gives $f(z=0)=1$ consistently with $\Omega _{\mathrm{m,0}}+\Omega _{\mathrm{DE,0}} + \Omega _{\sigma ^2,0}=1$. Note that we define the present time values of the density parameters as $\Omega _{i,0}=\rho _{i,0}/\rho _\mathrm{c,0}$, where $\rho _\mathrm{c,0}=3H_0^2\varphi _0^2/4$ is the present time value of the critical energy density, and that we eliminate $\Omega _\mathrm{M,0}$, that appears in f(z) originally, using the relation

$$\begin{aligned} \begin{aligned} \Omega _{\mathrm{M},0}=\gamma ^{-1}\left[ \Omega _{\mathrm{DE},0}-(\gamma -1)(\Omega _{\mathrm{m,0}}+\Omega _{\sigma ^2,0})\right] , \end{aligned} \end{aligned}$$

(44)

which can easily be obtained from (40). Using (44) along with $\Omega _{\mathrm{m,0}}+\Omega _{\sigma ^2,0}=1-\Omega _{\mathrm{DE,0}}$, we see that the present time density parameter of the effective DE is determined by the Jordan field’s mass and the BD parameter $\omega$ [viz., $\gamma (\omega )$] as⁸

$$\begin{aligned} \Omega _{\mathrm{DE},0}=\omega \frac{2}{3}\left( \frac{M}{H_0}\right) ^2+\frac{\gamma -1}{\gamma }. \end{aligned}$$

(45)

We obtain – upon straightforward arrangements in the equations obtained by using (29), (44) in the equations obtained by substituting (23), (30), (31) in (13), (14) – the principal EoS parameters of the effective DE along the x-axis and y- and z-axes in terms of the present day values of the density parameters as follows, respectively,

$$\begin{aligned} w_{\mathrm{DE},x}\equiv&\, \frac{p_{\mathrm{DE},x}}{\rho _{\mathrm{DE}}}=-1+\frac{1}{f(z)}\bigg [\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}\frac{2\omega +2}{2\omega +3}(1+z)^{3+\frac{1}{1+\omega }}\nonumber \\&\quad -\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}(1+z)^{3}-\frac{2}{2\omega +3}\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}(1+z)^{6+\frac{2}{1+\omega }}\nonumber \\&\quad -\frac{2}{3(1+\omega )}\sqrt{\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}} (1+z)^{3+\frac{1}{1+\omega }}\nonumber \\&\quad \times \sqrt{1+\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}(1+z)^{3}+\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}(1+z)^{6+\frac{2}{1+\omega }}} \bigg ], \end{aligned}$$

(46)

$$\begin{aligned} w_{\mathrm{DE},y}\equiv&\,\frac{p_{\mathrm{DE},y}}{\rho _{\mathrm{DE}}}=-1\nonumber \\&+\frac{1}{f(z)}\bigg [\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}\frac{2\omega +2}{2\omega +3}(1+z)^{3+\frac{1}{1+\omega }}-\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}(1+z)^{3}\nonumber \\&-\frac{2}{2\omega +3}\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}(1+z)^{6+\frac{2}{1+\omega }}\nonumber \\&+\frac{1}{3(1+\omega )}\sqrt{\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}}(1+z)^{3+\frac{1}{1+\omega }}\nonumber \\&\times \sqrt{1+\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}(1+z)^{3}+\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}(1+z)^{6+\frac{2}{1+\omega }}} \bigg ], \end{aligned}$$

(47)

which give

$$\begin{aligned} w_{\mathrm{DE},x,0}&=-1-\frac{\Omega _{\mathrm{m,0}}+2\,\Omega _{\sigma ^2,0}}{(2\omega +3)\Omega _{\mathrm{DE,0}}} -\frac{2\sqrt{\Omega _{\sigma ^2,0}}}{3(1+\omega )\Omega _{\mathrm{DE,0}}},\nonumber \\ \end{aligned}$$

(48)

$$\begin{aligned} w_{\mathrm{DE},y,0}&=-1-\frac{\Omega _{\mathrm{m,0}}+2\,\Omega _{\sigma ^2,0}}{(2\omega +3)\Omega _{\mathrm{DE,0}}} +\frac{\sqrt{\Omega _{\sigma ^2,0}}}{3(1+\omega )\Omega _{\mathrm{DE,0}}}\nonumber \\ \end{aligned}$$

(49)

for $z=0$.

The anisotropy of the EoS of the effective DE may be given as the difference between its principal EoS parameters ($\Delta w_{\mathrm{DE}}=w_{\mathrm{DE},y}-w_{\mathrm{DE},x}$) as

$$\begin{aligned} \begin{aligned} \Delta w_{\mathrm{DE}}=\,&\frac{1}{(1+\omega )}\frac{(1+z)^{3+\frac{1}{1+\omega }}}{f(z)} \sqrt{\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}} \\&\times \sqrt{1+\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{DE,0}}}(1+z)^{3}+\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}(1+z)^{6+\frac{2}{1+\omega }}}. \end{aligned} \end{aligned}$$

(50)

Accordingly, we have

$$\begin{aligned} \Delta w_{\mathrm{DE,0}}=\frac{1}{1+\omega }\,\frac{\sqrt{\Omega _{\sigma ^2,0}}}{\Omega _{\mathrm{DE,0}}}\quad \text {for}\quad z=0, \end{aligned}$$

(51)

which is approximately equal to $\frac{1}{1+\omega }\sqrt{\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{DE,0}}}}$, provided that the effective DE is dominant today ($z=0$). Additionally, we can write for dust domination

$$\begin{aligned} \Delta w_{\mathrm{DE}}\sim \frac{1}{(1+\omega )\gamma }\sqrt{\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{m,0}}}}\,(1+z)^{\frac{3}{2}}, \end{aligned}$$

(52)

which is positive definite as $\omega >0$, and for expansion anisotropy domination

$$\begin{aligned} \Delta w_{\mathrm{DE}}\sim -\frac{1}{1+\omega }\,\frac{1}{1-\gamma }, \end{aligned}$$

(53)

which depends only on $\omega$, is negative definite as $\omega >0$, monotonic in $\omega >0$ and takes values as $\Delta w_{\mathrm{DE}}=\{-2,-1.2\}$ for $\omega =\{0,\infty \}$.

We finally, using $\dot{\rho }_{\mathrm{DE}}+3H\rho _{\mathrm{DE}}(1+{\bar{w}}_{\mathrm{DE}})=0$, define a mean/volumetric effective EoS parameter for the effective DE as

$$\begin{aligned} \begin{aligned} {\bar{w}}_{\mathrm{DE}}=-1+\frac{1+z}{3}\frac{\rho _{\mathrm{DE}}'}{\rho _{\mathrm{DE}}}=-1+\frac{1+z}{3}\frac{f'(z)}{f(z)}, \end{aligned} \end{aligned}$$

(54)

which explicitly reads

$$\begin{aligned} \begin{aligned} {\bar{w}}_{\mathrm{DE}}=&-1+\frac{\Omega _{\mathrm{m},0}(1+z)^3\left[ \gamma \frac{3\omega +4}{3(1+\omega )}(1+z)^{\frac{1}{1+\omega }}-1\right] +\Omega _{{\sigma ^2},0}\frac{2(\gamma -1)(3\omega +4)}{3(1+\omega )}(1+z)^{6+\frac{2}{1+\omega }}}{\Omega _{\mathrm{DE,0}}+\Omega _{\mathrm{m,0}} \left\{ (1+z)^3\left[ \gamma (1+z)^{\frac{1}{1+\omega }}-1\right] +1-\gamma \right\} +\Omega _{\sigma ^2,0}(\gamma -1)\left[ (1+z)^{6+\frac{2}{1+\omega }}-1\right] }, \end{aligned} \end{aligned}$$

(55)

and

$$\begin{aligned} \begin{aligned} {\bar{w}}_{\mathrm{DE,0}}=&-1-\frac{1}{2\omega +3}\,\frac{\Omega _{\mathrm{m},0}}{\Omega _{\mathrm{DE,0}}}-\frac{2(5\omega +6)}{3(2\omega +3)(\omega +1)}\,\frac{\Omega _{\mathrm{\sigma ^2},0}}{\Omega _{\mathrm{DE,0}}} \end{aligned} \end{aligned}$$

(56)

for $z=0$. We note that, given $\Omega _{\mathrm{DE,0}}>0$ and $\Omega _{\mathrm{m,0}}$ and $\Omega _{{\sigma ^2},0}$ are already positive definite, for $\omega >0$, the contributions from both the pressureless matter and the expansion anisotropy to the present value of the mean EoS of the effective DE, ${\bar{w}}_{\mathrm{DE,0}}$, are negative, i.e., anisotropic BD extension of the standard $\Lambda$CDM model predicts phantom like effective DE in the present time Universe.

3.3 Preliminary investigations

3.3.1 Features of effective anisotropic dark energy

In this section, before the observational analysis, we discuss some of the features of the model, in particular those corresponding to the effective DE and the evolution of the expansion anisotropy. To do so, we start with justifying the positivity assumption on $\omega$ in our study by making use of one of the important parameters about the kinematics of the Universe, the deceleration parameter. We define the mean/volumetric deceleration parameter in terms of the average Hubble parameter as $q=-1+~\frac{H'}{H}(1+z)$. This reads, as we expect in a viable cosmology $\Omega _{{\sigma ^2},0}<\Omega _{\mathrm{r,0}}\ll \Omega _{\mathrm{m},0}$,

$$\begin{aligned} \begin{aligned} q\approx -1+\frac{4+3\omega }{2+2\omega }\,\frac{\Omega _{\mathrm{m},0}(1+z)^{3+\frac{1}{1+\omega }}}{1+\gamma \Omega _{\mathrm{m},0}\left[ (1+z)^{3+\frac{1}{1+\omega }}-1\right] } \end{aligned} \end{aligned}$$

(57)

for $z\sim 0$, and

$$\begin{aligned} \begin{aligned} q_0\approx -1+\frac{4+3\omega }{2+2\omega }\Omega _{\mathrm{m},0}, \end{aligned} \end{aligned}$$

(58)

for $z=0$ today. We confine the present study to small deviations from the standard $\Lambda$CDM model and hence it is reasonable to expect, in accordance with the most recent Planck results [12], that $\Omega _{\mathrm{m},0}\sim 0.3$ leading to $q_0\sim -0.55$ from $q=-1+\frac{3}{2}\Omega _{\mathrm{m}}$ in the standard $\Lambda$CDM, which in turn implies from (58) that we should have $|\omega | \gtrsim 10$, which leads to $\omega \gtrsim 10$ since have already assumed $\omega \ge 0$.⁹ Unless otherwise, for instance, if $\omega =0$, then $q_0\sim -0.4$, which could be decreased to a reasonable value, viz., $q_0 \sim -0.55$, by decreasing $\Omega _{\mathrm{m},0}$ considerably, namely $\Omega _{\mathrm{m},0}\sim 0.22$, which implies very large deviations from the standard $\Lambda$CDM. Therefore, our assumption $\omega \ge 0$ already allows large deviations from the standard $\Lambda$CDM, in other words, we do not oversimplify the model but yet basically avoid the cases most likely non-viable that would lead to extended discussions since the model exhibits various complicated dynamics particularly at small negative values of $\omega$, namely, when $\omega \sim -1$. Hence, in this section, we shall depict some key/interesting features of the anisotropic BD model and its deviation from the standard $\Lambda$CDM model assuming $\omega \ge 0$.

First of all, one may check easily that, similar to the standard $\Lambda$CDM, its anisotropic BD extension asymptotically approaches to de Sitter Universe in the far future,

$$\begin{aligned} H^2\rightarrow \frac{4}{3}\frac{\gamma }{\varphi _0^2}\rho _{\mathrm{M}}=\frac{2\gamma }{3}\omega M^2,\,\,\, \rho _{\sigma ^2}\rightarrow 0\,\,\,\text {and} \,\,\, \rho _{\mathrm{m}}\rightarrow&0 \nonumber \\ \,\,\,\text {for}\,\,\, z\rightarrow&-1, \end{aligned}$$

(59)

and accordingly, at this limit, the effective DE mimics exactly the cosmological constant as

$$\begin{aligned} \rho _{\mathrm{DE}}\rightarrow \gamma \rho _{\mathrm{M}},\,\,\, \bar{w}_{\mathrm{DE}}\rightarrow -1\quad \text {and}\quad \Delta w_{\mathrm{DE}}\rightarrow 0. \end{aligned}$$

(60)

At low redshifts, viz., at which $\rho _{\mathrm{m}}\simeq \rho _{\mathrm{m},0}$ and $\rho _{\sigma ^2}\ll \rho _{\mathrm{m}}$, the effective DE can approximately be described by

$$\begin{aligned} \begin{aligned}&\rho _{\mathrm{DE}}\simeq \gamma \rho _{\mathrm{M}}+(\gamma -1)\rho _\mathrm{m }\quad \text {and}\quad \\&\bar{w}_{\mathrm{DE}}\lesssim -1+\frac{(\gamma -1)\rho _\mathrm{m }}{\gamma \rho _{\mathrm{M}}+(\gamma -1)\rho _\mathrm{m }}\quad \text {for}\quad z\simeq 0. \end{aligned} \end{aligned}$$

(61)

Here, we first note that the non-negativity of the effective DE today ($z=0$), $\rho _{\mathrm{DE,0}}\ge 0$, leads to a natural lower bound on the value of the energy density corresponding to the effective cosmological constant (viz., $\rho _{\mathrm{M}}$) as $\rho _{\mathrm{M}}\ge \left( \frac{1}{\gamma }-1\right) \rho _{\mathrm{m},0}$, which allows zero lower bound only at the GR limit, i.e., at $\omega \rightarrow \infty$ (viz., $\gamma \rightarrow 1$). We next note that $\gamma$ takes values within the range $\frac{1}{2}\le \gamma < 1$ for $\omega \ge 0$ and hence the coefficient of $\rho _{\mathrm{m}}$ is negative definite, $\gamma -1< 0$, which implies that, because $\rho _{\mathrm{m}}$ decreases as the Universe expands, $\rho _{\mathrm{DE}}$ will increase as the Universe expands at $z\sim 0$, i.e., the effective DE will behave like a phantom fluid ($\bar{w}_{\mathrm{DE}}\lesssim -1$) at $z\sim 0$ [see (56) for $\bar{w}_{\mathrm{DE}}(z=0)$]. Namely, we see from (40) or (55) that given $\rho _{\mathrm{m}}=\rho _{\mathrm{m,0}}(1+z)^3$, the effective DE passes below the phantom-divide-line (PDL, $w_{\mathrm{DE}}=-1$) at

$$\begin{aligned} z=z_\mathrm{PDL}\equiv \left( \frac{2\omega +3}{2\omega +2}\right) ^{1+\omega }-1, \end{aligned}$$

(62)

and afterwards stays there forever, i.e., $\bar{w}_{\mathrm{DE}}<-1$ for $z<z_\mathrm{PDL}$. It is worth noting that the BD gravity predicts almost a specific redshift for the PDL crossing, such that

$$\begin{aligned} \frac{1}{2}\le z_\mathrm{PDL}\le e^{\frac{1}{2}}-1=0.65\quad \text {for}\quad \omega \ge 0, \end{aligned}$$

(63)

and $0.63\le z_\mathrm{PDL}<0.65$ for $\omega \ge 10$ (see [100] for a similar result). This is indeed an interesting result, such that $z_\mathrm{PDL}\sim 0.6$ can be taken as the signature of BD gravity. In other words, observations suggesting the presence of a DE source passing below PDL at $z\sim 0.6$ with high confidence would imply a strong reason for favoring the BD gravity over GR, or vice versa. During this period, $\Delta w_{\mathrm{DE}}$ can approximately be given by (51) and we note that the effective DE will be almost isotropic, $\Delta w_{\mathrm{DE}}\approx 0$, for $\omega >0$ since $\Omega _{\sigma ^2,0}\ll \Omega _{\mathrm{DE,0}}$ in a realistic cosmology.

At moderate redshifts (e.g., $z_\mathrm{PDL} \ll z \lesssim z_{\mathrm{eq}}$), namely, when the contribution from the pressureless matter to the EoS parameter of the effective DE given in (55) becomes significant but those from the expansion anisotropy and the mass of the Jordan field are not, the mean/volumetric EoS parameter of effective DE yields a plateau where

$$\begin{aligned} \begin{aligned} \rho _{\mathrm{DE}}\sim \rho _{\mathrm{m,0}}\gamma \left( \frac{\rho _{\mathrm{m}}}{\rho _{\mathrm{m,0}}}\right) ^{\frac{3\omega +4}{3\omega +3}} \quad \text {as}\quad {\bar{w}}_{\mathrm{DE}}\sim \frac{1}{3\omega +3}\\ \quad \text {for}\quad z_\mathrm{PDL} \ll z \lesssim z_{\mathrm{eq}}. \end{aligned} \end{aligned}$$

(64)

This plateau is persistent all the way to the large redshift values at which either the radiation becomes dominant over pressureless matter ($z\sim z_{\mathrm{eq}}$), or the contribution from the expansion anisotropy becomes significant. The plateau of ${\bar{w}}_{\mathrm{DE}}$ is located within the range $0\le {\bar{w}}_{\mathrm{DE}}\le \frac{1}{3}$ depending on the value of $\omega \ge 0$. For $\omega =0$, the effective DE behaves like an extra relativistic degree of freedom as ${\bar{w}}_{\mathrm{DE}}=\frac{1}{3}$ during this period and hence requires special attention. On the other hand, as it is discussed above, to obtain a viable model (viz., which deviates from the standard $\Lambda$CDM model only slightly), we expect $\omega \gtrsim 10$ and this places the plateau in a value in the range $0<{\bar{w}}_{\mathrm{DE}}\lesssim 0.03$. During this period, the evolution of the anisotropy of the effective DE can approximately be given by (52). We note that, although $\Delta w_{\mathrm{DE}}$ increases as $(1+z)^{\frac{3}{2}}$ during this period, the EoS parameter of the effective DE would not achieve a significant anisotropy during this period even under the weakest observational constraints we give in Sect. 4.

At very/sufficiently large redshifts, viz., when the only dominant contribution is from the expansion anisotropy in (55), the ratio of the energy densities of the effective DE and the expansion anisotropy freezes out;

$$\begin{aligned} \begin{aligned} \rho _{\mathrm{DE}}\sim (\gamma -1)\rho _{\sigma ^2}\quad \text {and}\quad \bar{w}_{\mathrm{DE}}\sim \frac{3\omega +5}{3\omega +3}. \end{aligned} \end{aligned}$$

(65)

We note that, except the GR limit $\omega \rightarrow \infty$ (i.e., $\gamma \rightarrow 1$), the mean/volumetric EoS parameter of the effective DE yields a second plateau placing in the range $1\lesssim \bar{w}_{\mathrm{DE}}\lesssim \frac{5}{3}$ for $\omega \ge 0$, and the effective DE behaves like the expansion anisotropy with a value of magnitude of ratio in the range $0<|\rho _{\mathrm{DE}}/\rho _{\sigma ^2}|\lesssim \frac{1}{2}$ but yields negative energy density since the coefficient $\gamma -1$ is negative for $\omega \ge 0$. According to this the energy density of the effective DE changes sign and we calculate from (40) that this occurs at

$$\begin{aligned} z=z_\mathrm{DE, pole}\sim \left( \frac{\gamma }{1-\gamma }\,\, \frac{\Omega _{\mathrm{m},0}}{\Omega _{{\sigma ^2},0}}\right) ^ {\frac{\omega +1}{3\omega +4}}-1, \end{aligned}$$

(66)

at which its EoS parameter exhibits a pole since at that redshift its energy density becomes zero. For $z>z_\mathrm{DE, pole}$, $\Delta w_{\mathrm{DE}}$ asymptotically approaches to (53) with increasing redshift and the EoS of the effective DE will be significantly anisotropic, viz., we have $1<\bar{w}_{\mathrm{DE}}\le \frac{5}{3}$ and $-2\le \Delta w_{\mathrm{DE}}<-1.2$ for $\omega \ge 0$. However, for $\omega \ge 0$, this interesting behavior is always irrelevant to the observable Universe since it occurs when the Universe is strongly dominated by the expansion anisotropy: Substituting $z=z_\mathrm{DE, pole}$ into $\frac{\rho _{\sigma ^2}}{\rho _{\mathrm{m}}}=\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{m,0}}} (1+z)^{3+\frac{2}{1+\omega }}$, which may be read off from (42), we obtain

$$\begin{aligned} \frac{\rho _{\sigma ^2}(z=z_\mathrm{DE,pole})}{\rho _{\mathrm{m}}(z=z_\mathrm{DE,pole})}=\left( \frac{\gamma }{1-\gamma }\right) ^{\frac{3\omega +5}{3\omega +4}}\left( \frac{\Omega _{\mathrm{m,0}}}{\Omega _\mathrm{\sigma ^2,0}}\right) ^{\frac{1}{3\omega +4}}. \end{aligned}$$

(67)

We note that, in a realistic cosmological setup, say, $\omega \gtrsim 10$ (even for $\omega \ge 0$) and $\Omega _{\mathrm{m,0}}\gg \Omega _{\sigma ^2,0}$, this ratio is obviously much larger than unity, implying that $z_\mathrm{DE,pole}$ occurs always at a redshift at which the expansion anisotropy is already dominant. If we consider the presence of radiation as well, then we have the following: During the radiation domination BD gravity mimics GR (as Jordan field $\varphi$ freezes out) so that the effective DE will be constant (i.e., mimics cosmological constant) while the expansion anisotropy keeps on growing as $\rho _{\sigma ^2}\propto (1+z)^6$ as z increases. Eventually, the expansion anisotropy will dominate over radiation and the dynamical effective DE – viz., $\rho _{\mathrm{DE}}\propto \rho _{\sigma ^2}(\gamma -1)$ as may be seen from (40) by setting $\rho _{\mathrm{M}}=0=\rho _{\mathrm{m}}$ – will show up again and then its EoS parameter will respectively realize a pole and a plateau that follows that pole as z increases. Thus, in a realistic setup, the pole and the second plateau features of the effective DE occur in the expansion anisotropy dominated universe that takes place before (in terms of time) the radiation dominated universe and hence they are irrelevant to the observable universe.

3.3.2 Modified expansion anisotropy

Almost all of the model dependent observational constraints on the expansion anisotropy in the literature, e.g., [91,92,93,94,95], consider GR+isotropic sources leading to the steep redshift dependency of the energy density corresponding to the shear scalar as $\rho _{\sigma ^2} \propto (1+z)^{6}$ and thereby give $\Omega _{\sigma ^2,0}\lesssim 10^{-21}$. On the other hand, the direct/model independent observational constraints from, e.g., type Ia SNe observations are rather weak, e.g., $\Omega _{\sigma ^2}\lesssim 10^{-4}$ for $z\sim 0$ [96, 97]. Such large values may be possible provided that the redshift dependence of $\rho _{\sigma ^2}$ is modified properly (viz., is made modest), which may be done, in principle, by either introducing an anisotropic source, e.g., anisotropic DE, in GR or considering usual cosmological sources in modified theories of gravity such as BD gravity that leads to an effective anisotropic source.

As we have shown in Sect. 3.1 that switching to the massive BD gravity (1) with $\omega \ge 0$ results in having more stringent constraints on the expansion anisotropy [see (37) and discussion that follows]. These constraints in fact can be weakened for $\omega <-1$ (since the expansion anisotropy in this case yields flatter redshift dependency) and even vastly for $-\frac{5}{3}<\omega <-1$, which however would cost to large deviations from GR, hence from $\Lambda$CDM model, such that in this case we will have rapidly varying cosmological gravitational coupling strength [see (24)] or, in the approach realized in (39), to an effective DE significantly deviating from $\Lambda$ [see (40) along with (27)]. It may be interesting to discuss a bit more on how the expansion anisotropy can be eased down in case $\omega <-1$ by considering (39) based on the approach that can be described by (8). We see from (65) – or directly from (26) – that the effective EoS parameter corresponding to the expansion anisotropy, viz., shear scalar, is

$$\begin{aligned} w_{\sigma ^2}=1+\frac{2}{3\omega +3}. \end{aligned}$$

(68)

We first note that it approaches Zeldovich fluid, i.e., $w_{\sigma ^2}\rightarrow 1$, at the GR limit $|\omega |\rightarrow \infty$ and, for $\omega \ge 0$ (finite), is in the range $1< w_{\sigma ^2}\le \frac{5}{3}$, which is stiffer than the Zeldovich fluid. Zeldovich fluid yields the most rigid EoS parameter ($w=1$) compatible with the requirements of relativity theory [144], which in turn implies that we can not have a minimally interacting source whose energy density redshift dependence is steeper than $(1+z)^{6}$ within GR, whereas we have it in BD gravity from the expansion anisotropy as an effective source provided that $\omega >-1$ as can be seen from (27). We depict $w_{\sigma ^2}$ versus $\omega$ in Fig. 1.

Fig. 1

Effective equation of state parameter corresponding to the expansion anisotropy $w_{\sigma ^2}$ versus Brans–Dicke parameter $\omega$. $w_{\sigma ^2}$ exhibits a pole at $\omega =-1$, and $w_{\sigma ^2}\rightarrow 1$ as $|\omega |\rightarrow \infty$

For $\omega \le -\frac{4}{3}$, the EoS parameter spanning the range $-1\le w_{\sigma ^2}<1$, which is familiar from the conventional cosmology: (i) $\omega =-2$ ($w_{\sigma ^2}=\frac{1}{3}$), the expansion anisotropy mimics radiation, viz., plays the role of an extra relativistic species. (ii) $\omega =-\frac{5}{3}$ ($w_{\sigma ^2}=0$), it mimics dust, viz., CDM. (iii) $\omega =-\frac{3}{2}$ –the scale invariant limit – ($w_{\sigma ^2}=-\frac{1}{3}$), it mimics constant negative spatial curvature (hyperboloid space). (iv) $\omega =-\frac{4}{3}$ ($w_{\sigma ^2}=-1$), it mimics $\Lambda >0$. It might be interesting to remind here that $-\frac{3}{2}<\omega <-\frac{4}{3}$ corresponds to the range in which the universe exhibits accelerating expansion in the presence of only pressureless matter (see footnote 4). The BD parameter $\omega \sim -\frac{4}{3}$ appears in d-branes string models [104, 105] and hence it is conceivable that such models would predict $\sigma ^2\sim \mathrm{const}$. For $-\frac{4}{3}<\omega <-1$, we have $w_{\sigma ^2}<-1$, i.e., the expansion anisotropy mimics phantom sources implying that the expansion anisotropy decreases/increases with increasing/decreasing redshift in contrast to the all other cases. The case $\omega =-1$ (the low energy effective string action limit of BD gravity [104, 105]) is interesting that $w_{\sigma ^2}$ exhibits a pole, namely, $\lim _{\omega \, \rightarrow \, -1^+}w_{\sigma ^2}=+\infty$ and $\lim _{\omega \, \rightarrow \, -1^-}w_{\sigma ^2}=-\infty$, which imply that we must set $\sigma ^2=0$ in this case (at least in our solution). In relevance with this, in the case of $\omega >-1$ but $\omega \approx -1$, $w_{\sigma ^2}$ will be extremely large, implying that in this case the presence of the expansion anisotropy today or at any moment in the observable Universe will require extreme fine tuning. Finally, for $\omega =0$ – the lowest value we consider in this study – we have $\omega _{\sigma ^2}=\frac{5}{3}$ and for $\omega \gtrsim 10$, which, as we discussed above, is required for small deviations from the $\Lambda$CDM model (by keeping isotropic RW metric), we have $1<w_{\sigma ^2}\lesssim 1.06$. These imply that considering BD gravity with $w\gtrsim 10$ rather than GR as the law of gravity will additionally lead to small modification in the evolution of the expansion anisotropy (viz., a slightly steeper redshift dependency of the shear scalar compared to the one in GR) in the anisotropic generalization of the $\Lambda$CDM model in contrast to DE models that can be described by scalar field/s within GR.

4 Constraints from recent cosmological data

In this section we perform a parameter estimation and provide observational constraints on the free parameters of the model: the pressureless matter density parameter today $\Omega _{\mathrm{m,0}}$, the baryon density parameter today $\Omega _\mathrm{b,0}h^2$, and the dimensionless Hubble constant h as well as the two free parameters $\omega$ and $\Omega _{\sigma ^2,0}$ that account for the anisotropic BD extension of the $\Lambda$CDM model. We consider H(z) given in (34) rather than (39), since the former one contains radiation (viz., $\Omega _{\mathrm{r,0}}$) as well and can be safely used all the way to the recombination redshift (covering last scattering surface), so that we could include CMB data in our analysis. For the derived parameters we present the ones relevant to the effective DE. We consider (39), where the effective DE is defined in accordance with the exact explicit solution in redshift given in Sect. 3 and the approach described by (8). In this way, even though this solution ignores the presence of radiation, we could go further and investigate the dynamics of the model under consideration in terms of redshift for $z<z_{\mathrm{eq}}$ (during which radiation is negligible) analytically in the light of the observational constraints. We keep in mind that H(z) throughout this study is averaged over the volumetric expansion rate of the Universe and hence the method we use constrains the expansion anisotropy through its contribution to the average expansion rate of the Universe (viz., the shear scalar $\sigma ^2$, which contributes to the H(z) directly via the term $\rho _{\sigma ^2,0}(1+z)^{6+\frac{2}{1+\omega }}$ and also indirectly via its effect on the $\rho _{\mathrm{DE}}$ as a result of BD theory).

Table 1 Constraints on the anisotropic Brans–Dicke extension of the standard $\Lambda$CDM model using the combined data sets PLK+BAO+SN+H. For two-tailed distributions, the results are given in $1\sigma$ and for one-tailed distributions, the results are given in $2\sigma$. $\Lambda$CDM column corresponds to GR limit (Brans–Dicke gravity with $\omega \rightarrow \infty$) and the shear scalar $\sigma ^2=0$ (isotropic). Parameters and ranges of the uniform priors assumed in our analysis and derived parameters are labeled with $^*$

In order to perform the parameter space exploration, we make use of a modified version of the simple and fast Markov Chain Monte Carlo (MCMC) code that computes expansion rates and distances from the Friedmann equation, named SimpleMC [21, 145]. For an extended review of cosmological parameter inference see [146]. The code uses a compressed version of the Planck data (PLK), where the CMB is treated as a “BAO experiment” at redshift $z=1090$, measuring the angular scale of the sound horizon at that time, a recent analysis of Type Ia supernova (SN) data dubbed Joint Light-curve Analysis (JLA) compressed into a piece-wise linear function fit over 30 bins spaced evenly in $\log z$, and high-precision Baryon Acoustic Oscillation measurements (BAO), from the comoving angular diameter distance only (viz., growth rate measurement is not considered), the Hubble distance and the volume averaged distance, at different redshifts up to $z=2.36$. For a more detailed description about the datasets used see [21]. We also include, independent of/separate from other data sets, a collection of currently available data on H(z) obtained from cosmic chronometers (H) (see [143] and references therein). Table 1 displays the parameters used throughout this paper along with the corresponding flat priors; derived parameters labeled with $^*$.

In the analysis, we assume radiation content by including three neutrino species ($N_\mathrm{eff}=3.046$) with minimum allowed mass $\sum m_{\nu }=0.06\, \mathrm{eV}$ theoretically well determined within the standard model of particle physics and the density parameter of radiation $\Omega _{\mathrm{r},0}=2.469\times 10^{-5} h^{-2}(1+0.2271 N_\mathrm{eff})$, where dimensionless Hubble constant $h=H_0/100\, \mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}$ [147]. Throughout the analysis we assume flat priors over our sampling parameters: $\Omega _{\mathrm{m},0}=[0.05,1.5]$, $\Omega _{\mathrm{b},0} h^2=[0.02,0.025]$ and $h=[0.4,1.0]$. Whereas priors for the two parameters that characterize the anisotropic BD extension of the standard $\Lambda$CDM model are given by $\log _{10}\omega =[0,6]$ and $\log _{10}\Omega _{\sigma ^2,0}= [-12,0]$. The photon energy density today $\rho _{{\gamma },0}$ is not subject to our analysis since it is well constrained, such that it has a simple $\rho _{\gamma }=\frac{\pi ^2}{15} T_\mathrm{CMB}^4$ relation with the CMB monopole temperature [148], which is very precisely measured to be $T_{\mathrm{CMB},0}=2.7255\pm 0.0006\,\mathrm{K}$ [149]. Table 1 summarizes the observational constraints on the free parameters as well as the derived parameters of the model under consideration using the combined datasets PLK+BAO+SN+H; and for comparison those parameters used on the standard $\Lambda$CDM model.

Fig. 2

(Left panel): 2D marginalized posterior distribution of the parameters $\log _{10}\Omega _{\sigma ^2}$ and $\log _{10}\omega$ that determine the extension from the standard $\Lambda$CDM model. The 2D constraints are plotted with 1$\sigma$ and 2$\sigma$ confidence contours. (Right panel): the 3D posterior distribution in the $\{\bar{w}_{\mathrm{DE,0}}, ~\Delta w_{\mathrm{DE,0}}, ~\log _{10}\omega \}$ subspace, where the colour code indicates the value of $\log _{10}\omega$ using the colour bar

Fig. 3

(Left panel) 2D marginalised posterior distribution of the parameters corresponding to the mass of the Jordan field M and the Brans–Dicke parameter $\log _{10}\omega$. (Right panel) 1D constraints corresponding to the redshift of matter-radiation equality ($z_{\mathrm{eq}}$)

Fig. 4

The constraints on g(z) as a result of the data. These show the posterior probability Pr(g|z): the probability of g as normalized in each slice of constant z, with color scale in confidence interval values. The 1$\sigma$ and 2$\sigma$ confidence intervals are plotted as black lines and red lines correspond to the best-fit values found over the analysis, while dotted lines to the redshift of matter-radiation equality. (Top) $\bar{w}_{\mathrm{DE}}(z)$ at low z values (left) and at high redshift values in log-scale (right). (Middle) $\Delta w_{\mathrm{DE}}(z)$ at low z values (left) and at high redshift values in log-scale (right). (Bottom) $\Omega _{\mathrm{DE}}(z)$ at low z values (left) and at high redshift values in log-scale (right)

Figure 2 summarizes the constraints on the parameters that characterize the anisotropic BD extension to the standard $\Lambda$CDM model. Left panel of this figure displays 2D marginalized posterior distributions on the BD parameter $\omega$ (deviation from GR) and the density parameter corresponding to the expansion anisotropy today $\Omega _{\sigma ^2,0}$ (deviation from isotropic expansion). We see from Table 1 that current observations prefer small contributions from the expansion anisotropy $\log _{10}\Omega _{\sigma ^2,0} < -8.48$ (95% C.L.) with a large BD parameter (small deviation from GR) $\log _{10}\omega >1.69$ (95% C.L.), signaling that the anisotropic BD extension to the standard $\Lambda$CDM model should be studied indeed as a small correction in line with what we have been claiming and basically doing so far. Right panel of Fig. 2 displays 3D posterior distributions in the $\{\bar{w}_{\mathrm{DE,0}}, \Delta w_{\mathrm{DE,0}}$, $\log _{10}\omega \}$ subspace: the average/volumetric EoS parameter of the effective DE $\bar{w}_{\mathrm{DE,0}}$ and the measure of anisotropy of its EoS parameter $\Delta w_{\mathrm{DE,0}}$ coloured code with $\log _{10}\omega$. The former two are derived parameters that determine deviations of the effective DE from $\Lambda$ ($\bar{w}_{\mathrm{DE}}=-1$ and $\Delta w_{\mathrm{DE}}=0$) today and are controlled by $\omega$. We notice that for large values $\omega$ (color coded with red) we recover the $\Lambda$CDM case at the present time, where $\Delta w_{\mathrm{DE,0}}=w_{\mathrm{DE},y,0}-w_{\mathrm{DE},x,0} \rightarrow 0$ and $\bar{w}_{\mathrm{DE,0}} \rightarrow -1$, as expected by our previous analysis. We would like to remind here that in our observational analysis only positive values of $\omega$ were allowed. The left panel of Fig. 3 shows that current data impose upper bounds on the mass of the Jordan field given by $M<1.51 \times 10^{-34}$eV (95% C.L.). The correlation observed on the parameters $\omega$ and M is an interesting point to note. Even though both parameters have only one tail constraints, the effective cosmological constant $2\omega M^2$ –a component of the effective DE, see (31) and (40)– is well limited by $2\omega M{^2}=4.48\pm 1.22\times 10^{-66}$eV$^2$ and in agreement with the current value of $\Lambda =4.48 \pm 1.25\times 10^{-66}$eV$^2$, see Table 1. The extra two parameters from the anisotropic BD extension to the $\Lambda$CDM model may shift and relax some of the constraints, as it is the case of $\Omega _{\mathrm{m,0}}=0.3002 \pm 0.0067$, and therefore as a consequence the earlier redshift at matter-radiation equality $z_{\mathrm{eq}}=3368.62\pm 30.89$, compared to the $\Lambda$CDM model with $z_{\mathrm{eq}} = 3359.45\pm 24.14$ (see right panel of Fig. 3).

We depict, in Fig. 4, the behaviour of the effective DE (top and middle panels) as well as the evolution of its density parameter $\Omega _{\mathrm{DE}}=\rho _{\mathrm{DE}}/3H^2$ (bottom panel) in redshift, according to the constraints obtained from observational analysis (see Sect. 3.3.1 for the features of the effective DE). Left panel of Fig. 4 displays the behaviors of $\bar{w}_{\mathrm{DE}}$, $\Delta w_{\mathrm{DE}}$ and $\Omega _{\mathrm{DE}}$ for redshifts up to $z=2$. This panel was drawn by taking random samples from the posterior distribution of the parameter-space and colour coded by its likelihood (bluer colors represent regions of higher probability). Right panel of Fig. 4 displays the extended behaviors of the same functions; from redshift $z=1$ to $z=10^6$ for $\bar{w}_{\mathrm{DE}}$ and from redshift $z=10$ to $z=10^6$ for $\Delta w_{\mathrm{DE}}$ and $\Omega _{\mathrm{DE}}$. Though, note that the extensions to the redshift values $z\sim z_{\mathrm{eq}}$ and larger should be seen with a caution that these are obtained by making use of the effective DE derived in Sect. 3.3.1 without considering the presence of radiation. In this panel we depict the probability of a function normalized in each slice of constant z, with colour scale in confidence interval values. The 1$\sigma$ and 2$\sigma$ confidence intervals are plotted as black lines and red lines correspond to the best-fit values found over the analysis. At low redshifts, the upper-left panel shows that the crossing of the PDL occurs at $z_\mathrm{PDL}=0.64839\pm 0.00058$, just as we concluded from the discussion that follows equation (62). For redshift values $z<z_\mathrm{PDL}$, the effective DE exhibits phantom behaviour $\bar{w}_{\mathrm{DE}}<-1$ as may be seen from (56) or (61) and for higher redshift values, $z>z_\mathrm{PDL}$, $\bar{w}_{\mathrm{DE}}>-1$ as may be seen from (64) and (65). Evolving DE with an EoS parameter being below $-1$ at present, evolved from $w>-1$ in the past is named as quintom DE. We obtain the quintom DE as the upper-left panel shows, whereas, the explicit construction of quintom scenario is more difficult than other dynamical DE models, due to a no-go theorem which forbids the EoS parameter of a single perfect fluid or a single scalar field to cross the $w=-1$ boundary [43]. This property is distinctive that single-scalar-field models with canonical kinetic term are not allowed to satisfy, which also lead to that the Hamiltonian is unbounded from below. Interestingly this property is also achieved throughout some model independent analyses [23, 29, 30]. The middle left panel shows that, as can be deduced from (51) and (52), the effective DE becomes slightly more anisotropic with the increasing redshift, viz., the anisotropy of the EoS parameter $\Delta w_{\mathrm{DE}}$ increases only about an order of magnitude from its current value $\Delta w_{\mathrm{DE,0}}\lesssim 10^{-7}$ (see Table 1) to the one at $z=2$. Looking at the right panel we see that after few redshifts, as pressureless matter becomes dominant (see in the bottom panel that $\Omega _{\mathrm{DE}}\sim 0.1$ for $z\sim 1.75$ and $\Omega _{\mathrm{DE}}\sim 0$ for $z\gtrsim 10$), $\bar{w}_{\mathrm{DE}}$ starts to noticeably climb up and settles in the first plateau of $\bar{w}_{\mathrm{DE}}\sim 0$ (see (64) and paragraph covering it) lying between $z\sim 50$ and $z\sim z_{\mathrm{eq}}$ (viz., throughout the pressureless matter dominated epoch). There is a period in this plateau during which $\bar{w}_{\mathrm{DE}}>0$ [see (64)], viz., $\Omega _{\mathrm{DE}}$ increases (i.e., the energy density of the effective DE increases faster than of the pressureless matter) with increasing redshift. However, as can be seen in the bottom right panel, $\Omega _{\mathrm{DE}}$ during this period can never grow up to considerable values, viz., remains less than a percent at 1$\sigma$ C.L. and few percents at 2$\sigma$ C.L.. The anisotropy of the effective DE keeps on increasing with increasing redshift approximately in accordance with (52) during this plateau, but it remains positive definite and insignificant (e.g., $\Delta w_{\mathrm{DE}}\lesssim 0.07$ at photon decoupling redshift $z\sim 1100$) until the effective DE starts to leave this plateau as the expansion anisotropy starts to become dominant. We see that $\bar{w}_{\mathrm{DE}}$ exhibits a pole at $\log _{10}z_\mathrm{DE,pole}= 4.80\pm 0.58$ then settles in a new plateau of $\bar{w}_{\mathrm{DE}}\sim 1$ starting at $z\sim 10^5$ during which $\Delta w_{\mathrm{DE}}$ exhibits a similar behavior and settles down at $\Delta w_{\mathrm{DE}}\sim -1.2$, i.e., the effective DE eventually becomes highly anisotropic. However, note that this last stage of the effective DE starting just after $z\sim z_{\mathrm{eq}}$ is basically an artifact of omitting radiation in Sect. 3.3.1 to have explicit expression of $\rho _{\mathrm{DE}}(z)$, strickly speaking, is in fact unlikely to be realized at an observationally relevant past of the Universe (see the discussion starting with equation (65) in Sect. 3.3.1).

Instead, in a realistic setup, the radiation domination should start at $z=z_{\mathrm{eq}}$ and be maintained all the way to the redshift values at which BBN took place $z\sim z_{\mathrm{BBN}}$. During which the Jordan field is constant and hence, except an altered cosmological gravitational coupling strength, the Universe evolves exactly the same as in GR (see Sect. 3.1), such that the effective DE mimics $\Lambda$ (viz., $\bar{w}_{\mathrm{DE}}= -1$ and $\Delta w_{\mathrm{DE}}=0$ implying that it is isotropic and maintains its energy density value at $z\sim z_{\mathrm{eq}}$ for $z>z_{\mathrm{eq}}$) and expansion anisotropy increases as $\propto (1+z)^6$. However, although the data constrain $z_\mathrm{DE,pole}$ to be about 1-2 orders of magnitude larger than $z_{\mathrm{eq}}$, we see in Fig. 4 that, even within the 68% C.L. error region, the abrupt behaviour of the effective DE can start at redshift values less than $z_{\mathrm{eq}}$, which implies that the expansion anisotropy in this case is too large that it will never allow radiation to be dominant in the Universe, though over pressureless matter. Moreover, even there is a region where abrupt behaviour of the effective DE starts at redshift values larger than $z_{\mathrm{eq}}$, using the constraint on $\Omega _{\sigma ^2,0}$ from Table 1 we calculate that the expansion anisotropy dominates over radiation at redshift values smaller than $z_{\mathrm{BBN}}$, i.e., it will spoil the BBN processes. All these imply that the observational constraint method we performed here is able to put stringent enough constraints on neither $\Omega _{\sigma ^2,0}$ nor $\omega$. These two parameters are not independent and hence we should further investigate the upper boundary on $\Omega _{\sigma ^2,0}$ and lower boundary on $\omega$, of course, particularly, by using the data providing information about the Universe at $z\gtrsim z_{\mathrm{eq}}$, such as the peak of the matter-power spectrum relevant to $z=z_{\mathrm{eq}}$ and BBN relevant to $z\sim 10^8$.

5 Some further observational consequences and discussions

5.1 Variation of cosmological gravitational coupling strength

The normalized rate of change of the cosmological gravitational coupling strength in our exact solution (neglecting radiation), from (24), reads

$$\begin{aligned} \frac{\dot{G}}{G}=-\frac{H}{1+\omega }, \end{aligned}$$

(69)

which is negative definite considering $\omega \ge 0$ in this study, provided that the Universe is expanding ($H>0$). The constraints we found on $\omega$ and $H_0$ (see Table 1) predict

$$\begin{aligned} \bigg |\frac{\dot{G}}{G}\bigg | <1.092\times 10^{-12} \, \mathrm{yr}^{-1}\,\, (68\% \,\mathrm{C.L.}) \quad \mathrm{for}\quad z=0, \end{aligned}$$

(70)

similar to those given for $z\sim 0$ from various physical systems in which gravity is not negligible, such as the motion of the bodies of the Solar System, astrophysical and cosmological systems (see [150] for a comprehensive review). Assuming this relation (69) holds all the way to matter-radiation equality, we obtain $|\dot{G}/G|\lesssim 10^{-8} \, \mathrm{yr}^{-1}$ for $z\sim z_{\mathrm{eq}}$. On the other hand –given that $\varphi =\mathrm{const.}$ (and hence $G=\mathrm{const.}$) is an attractor solution when radiation is dominant– the redshift dependence of G becomes flatter w.r.t. $G\propto (1+z)^{\frac{1}{1+\omega }}$ [see (24)] as the radiation becomes more significant w.r.t. pressureless matter and thereby $|\dot{G}/G|\lesssim 10^{-8} \, \mathrm{yr}^{-1}$ for $z\sim z_{\mathrm{eq}}$ will never be achieved, but instead G will become almost constant for $z\sim z_{\mathrm{eq}}$ and remain so for $z>z_{\mathrm{eq}}$ as long as radiation continues to be dominant. Yet, in accordance with our detailed discussion in Sect. 3.1, (69) is mostly valid all the way to $z_{\mathrm{eq}}$ and leads to $G_1\sim G_0(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }}$ [see (33)]. Accordingly, we find that the constraints from the data predict the upper bound on the relative change in the strength of the cosmological gravitational coupling at $z\sim z_{\mathrm{eq}}$ w.r.t. its present time value as $\frac{\Delta G}{G_0}|_{z\sim z_{\mathrm{eq}}}=0.138\,(0.175)$ 68% C.L. (95% C.L.). This, in turn, implies that the strength of the cosmological gravitational coupling during the radiation dominance will be a constant $G_{1}$ satisfying the following constraints:

$$\begin{aligned} \begin{aligned} G_0\le G_{1}\lesssim 1.138\, G_0\,\,(1.175\, G_0)&\quad 68\% \,\mathrm{C.L.}\, (95\% \,\mathrm{C.L.}) \end{aligned} \end{aligned}$$

(71)

for $z\gtrsim z_{\mathrm{eq}}$. These constraints, of course, are valid also for the epoch of primordial nucleosynthesis that takes place when $z_{\mathrm{BBN}}\sim 3\times 10^8$, provided that the expansion anisotropy is still insignificant at that redshift. We finally note that these constraints are similar to those obtained from CMB and BBN (see Ref. [150] for a review on CMB and BBN constraints on $\Delta G/G_0$).

5.2 Matter-radiation transition

We find using (36) that, along with the constraints on the BD parameter and the radiation density parameter, the upper bound on the expansion anisotropy today, $\Omega _{\sigma ^2,0}=10^{-8.48}$ (see Table 1), leads to a lower bound on the expansion anisotropy-radiation equality ($\Omega _{\sigma ^2}=\Omega _{\mathrm{r}}$) redshift as $z_\mathrm{eq,\sigma ^2,r}\sim 10^3$, which is obviously not acceptable in a viable cosmological model since it implies that the expansion anisotropy dominates the Universe even at redshift values less than the matter-radiation equality ($\Omega _{\mathrm{m}}=\Omega _{\mathrm{r}}$) redshift $z_\mathrm{eq }\sim 3369$ (see Table 1). Indeed, using the results given in Table 1, we see that the Universe could be dominated by the expansion anisotropy as $\frac{\Omega _{\sigma ^2,\mathrm eq}}{\Omega _\mathrm{r,eq}+\Omega _\mathrm{m,eq}}=\frac{\Omega _{\sigma ^2,\mathrm eq}}{2\Omega _\mathrm{m,eq}}\lesssim 10^2$ at the matter-radiation equality, and thereby conclude that the upper bound on $\Omega _{\sigma ^2,0}$ given above should be further reduced to values guaranteeing that the Universe is radiation and matter dominated at the matter-radiation equality. Additionally, transition from radiation to matter domination is one of the most important epochs in the history of the Universe. This transition alters the growth rate of density perturbations: during the radiation era perturbations well inside the horizon are nearly frozen but once matter domination commences, perturbations on all length scales are able to grow by gravitational instability and therefore it sets the maximum of the matter power spectrum in GR as well as BD [138]. Namely, it determines the wavenumber, $k_{\mathrm{eq}}$, of a mode that enters the horizon, $H_{\mathrm{eq}}a_{\mathrm{eq}}$, at the matter-radiation transition [138, 151].

In our model, $k_{\mathrm{eq}}=H_{\mathrm{eq}}a_{\mathrm{eq}}=\frac{H_{\mathrm{eq}}}{(1+z_{\mathrm{eq}})}$ can be estimated analytically by assuming H(z) given in (34) holds all the way to matter-radiation equality. At this point, both the radiation and matter contribute equally to the total energy density. This gives us opportunity to reduce the constraints on $\Omega _{\sigma ^2,0}$. To do so, using (34), we write

$$\begin{aligned} k_{\mathrm{eq}}=\,&H_0\big [\Omega _\mathrm{M,0}(1+z_{\mathrm{eq}})^{-2}+2\Omega _{\mathrm{m,0}}(1+z_{\mathrm{eq}})^{1+\frac{1}{1+\omega }}\nonumber \\&+\Omega _\mathrm{\sigma ^2,0}(1+z_{\mathrm{eq}})^{4+\frac{2}{1+\omega }}\big ]^{1/2}, \end{aligned}$$

(72)

where we use $a_0=a(z=0)=1$, and $\Omega _\mathrm{M,0}=1-\Omega _{\mathrm{m,0}}-\Omega _{\mathrm{r,0}}-\Omega _\mathrm{\sigma ^2,0}$. We note that, for $\omega \gtrsim 50$ and $z_\mathrm{eq }\sim 3369$ (see Table 1), the contribution to $k_{\mathrm{eq}}$ from the term with $\Omega _\mathrm{M,0}$ is negligible, $k_{\mathrm{eq}}$ is by far more sensitive to $\Omega _\mathrm{\sigma ^2,0}$ and, for a given set of values of the parameters, $k_{\mathrm{eq}}$ increases with increasing $\Omega _\mathrm{\sigma ^2,0}$ and decreases with increasing $\omega$ (viz., as the BD gravity approaches to GR).

We obtain $k_{\mathrm{eq}}=0.01034 \pm 0.00012$ for $\Omega _{\sigma ^2,0}=0$ (isotropic expansion), which is slightly larger than, but yet consistent with the isotropic $\Lambda$CDM model value $k_{\mathrm{eq}}=0.01024 \pm 0.00007$. Accordingly, switching to the massive BD (1) leads only to a slight increase in $k_{\mathrm{eq}}$, as expected from the constraint $\omega \gtrsim 50$ (see Table 1). It may be useful to note that these two values are consistent with the recent Planck release [12] value $k_{\mathrm{eq}}=0.010339\pm 0.000063$ (TT,TE,EE+lowE+lensing+BAO) obtained for the base $\Lambda$CDM model. On the other hand, we see that the expansion anisotropy, although negligible today, shifts $k_{\mathrm{eq}}$ to unrealistically large values, viz., we obtain $k_{\mathrm{eq}}= 0.15159 \pm 0.00327$ for the upper bound $\Omega _{\sigma ^2,0}=10^{-8.48}$ (see Table 1). We then work out the values of $\Omega _{\sigma ^2,0}$ that can shift $k_{\mathrm{eq}}$ to reasonable values. See Fig. 5 showing explicitly $k_{\mathrm{eq}}$ with respect to $\Omega _{\sigma ^2,0}$ with errors to make a simple comparison between the models. We obtain $k_{\mathrm{eq}}=0.02824 \pm 0.00057$ by setting $\Omega _{\sigma ^2,0}=10^{-10}$, and $k_{\mathrm{eq}}=0.01067 \pm 0.00013$ by setting $\Omega _{\sigma ^2,0}=10^{-12}$, which are still inconsistent with the $k_{\mathrm{eq}}$ values given for the isotropic $\Lambda$CDM model. But then, we obtain $k_{\mathrm{eq}}=0.01037 \pm 0.00012$ by setting $\Omega _{\sigma ^2,0}=10^{-13}$, and $k_{\mathrm{eq}}=0.01034 \pm 0.00012$ by setting $\Omega _{\sigma ^2,0}=10^{-14}$, where we notice that only the last decimals are different. We observe further that $k_{\mathrm{eq}}$ does not change for $\Omega _{\sigma ^2,0}\lesssim 10^{-14}$ anymore and remains consistent with the isotropic $\Lambda$CDM model values obtained in this study as well as recent Planck release, which in turn implies that we cannot distinguish $\Omega _{\sigma ^2,0}\lesssim 10^{-14}$ from $\Omega _{\sigma ^2,0}=0$ by means of $k_{\mathrm{eq}}$. We finally calculate using $\Omega _{\sigma ^2,0}\sim 10^{-14}$ that the Universe is indeed dominated by matter+radiation at matter-radiation equality, viz., we now have $\frac{\Omega _{\sigma ^2,\mathrm eq}}{2\Omega _\mathrm{m,eq}}\sim 10^{-3}$. Thus, by means of matter-radiation transition, we conclude

$$\begin{aligned} \Omega _{\sigma ^2,0}\lesssim 10^{-14}\quad \text {along with}\quad \omega \gtrsim 50, \end{aligned}$$

(73)

where the upper bound on the expansion anisotropy is improved by reducing about six orders of magnitude w.r.t. the ones given in Table 1.

Fig. 5

The wavenumber ($k_{\mathrm{eq}}$) of a mode that enters the horizon ($H_{\mathrm{eq}}a_{\mathrm{eq}}$) at matter-radiation transition ($z_{\mathrm{eq}}$), viz., the maximum of matter power spectrum, for some fixed values of the expansion anisotropy today ($\Omega _{\sigma ^2,0}$). Anisotropic BD extension of the $\Lambda$CDM model, is labeled by green color, the standard base $\Lambda$CDM (in this study) is labeled by blue color, while the recent Planck release value for the standard base $\Lambda$CDM model is shown with red color

5.3 Big Bang Nucleosynthesis

Big Bang Nucleosynthesis provides a probe of the dynamics of the early Universe, which in turn gives us opportunity to further investigate the constraints on the anisotropic BD extension of the standard $\Lambda$CDM model. Such that, in the standard-BBN (SBBN) –assuming the standard model of particle physics is valid and the expansion of the Universe is governed by GR– the processes relevant to BBN take place when the temperature ranges from $T\sim 1\, \mathrm{MeV}$ to $T\sim 0.1\, \mathrm{MeV}$ and the age of the Universe from $t\sim 1\,\mathrm{s}$ to $t\sim 3\,\mathrm{min}$ corresponding to redshift $z\sim 3\times 10^8$ at which the Universe is radiation dominated. These, of course, should not be altered significantly in a viable cosmological model. Accordingly, we first looked for the condition of radiation dominance at $z\sim 3\times 10^8$ (implying $z_\mathrm{eq,\sigma ^2,r}\gtrsim 3\times 10^8$) by using the constraints on the relevant parameters from Table 1 and found out that $\Omega _{\sigma ^2,0}\lesssim 10^{-21}$ for $\omega \gtrsim 50$ in line with our preliminary investigation in Sect. 3.1 [cf., see Eq. (37)]. On the other hand, $z_\mathrm{eq,\sigma ^2,r}$ corresponds to the redshift at which $\rho _{\sigma ^2}/\rho _{\mathrm{r}}=1$, but as may be seen from the investigations in [91, 152] the expansion anisotropy does not lead to a considerable deviation from the SBBN for the $\rho _{\sigma ^2}/\rho _{\mathrm{r}}$ ratio up to a few percents, viz.,

$$\begin{aligned} \frac{\rho _{\sigma ^2}(z=z_{\mathrm{BBN}})}{\rho _{\mathrm{r}}(z=z_{\mathrm{BBN}})}\lesssim 10^{-2}. \end{aligned}$$

(74)

Hence, considering this condition, we can obtain a new constraint on $\Omega _{\sigma ^2,0}$ stronger than the one given in (37) obtained from the condition $z_\mathrm{eq,\sigma ^2,r}\gtrsim 3\times 10^8$. To work this out, we first write, from (35),

$$\begin{aligned} \frac{\rho _{\sigma ^2}}{\rho _{\mathrm{r}}}&\sim \frac{\rho _{\sigma ^2,0}(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }}(1+z)^6}{\rho _{\mathrm{r,0}}(1+z)^4}\nonumber \\&=\frac{\Omega _{\sigma ^2,0}}{\Omega _{\mathrm{m,0}}}(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }+1}(1+z)^2, \end{aligned}$$

(75)

where we also used $\frac{\Omega _{\mathrm{m,0}}}{\Omega _{\mathrm{r,0}}}=\frac{\rho _{\mathrm{m,0}}}{\rho _{\mathrm{r,0}}}=1+z_{\mathrm{eq}}$. Then, we use in this equation the condition (74) along with the constraints on $\Omega _{\mathrm{m,0}}$ and $z_{\mathrm{eq}}$ given in Table 1 by setting $z=3\times 10^8$ and reach the following constraint

$$\begin{aligned} \Omega _{\sigma ^2,0}\lesssim 10^{-23} \quad \text {along with}\quad \omega \gtrsim 50. \end{aligned}$$

(76)

We note that, here, the contribution from the effective anisotropic pressure of the Jordan field during $z<z_{\mathrm{eq}}$, which is encoded in the term $(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }}$ in (75), is not significant, such that it leads to only fifteen percent smaller upper bound value for $\Omega _{\sigma ^2,0}$, viz., $(1+z_{\mathrm{eq}})^{-\frac{1}{1+\omega }}\lesssim 10^{-0.07}=0.85$ for $\omega \gtrsim 50$. Thus, provided that the condition on the expansion anisotropy today given in (76) is satisfied, the expansion anisotropy will not have considerable effect on BBN, but the altered strength of the cosmological gravitation coupling depending on the BD parameter $\omega$ will, as we shall investigate in what follows.

Provided that the condition (76) is satisfied, the Universe is radiation dominated during the BBN epoch and hence, as it is discussed in Sect. 3.1, we have $a\propto t^{\frac{1}{2}}$ as in the SBBN based on GR except that the strength of the cosmological gravitational coupling during this epoch, $G_1$, can be slightly larger than its present time value, $G_0$, in our model based on BD gravity, see (71) and (33). Consequently, in accordance with (33), we can write the expansion rate of the Universe during the radiation dominance, hence during the BBN as well, as follows:

$$\begin{aligned} H^2=\frac{G_1 }{G_0}H^2_\mathrm{SBBN}\quad \text {with}\quad H^2_\mathrm{SBBN}=\frac{8\pi G_0 }{3}\frac{\pi ^2}{30}g_{*}T^4, \end{aligned}$$

(77)

where T is the temperature and $g_{*}(T)$ is the effective number of degrees of freedom counting the number of relativistic particle species determining the energy density in radiation as $\rho _\mathrm{r}=\frac{\pi ^2}{30}g_{*}T^4$. According to this, BD gravity can lead to a larger expansion rate for a given temperature since $G_1>G_0$ is allowed, see (71). We note that this is analogue of altering the expansion rate of the Universe during BBN by modifying $g_*$ (for instance, by introducing an extra massless degree of freedom such as sterile neutrino) within GR as

$$\begin{aligned} H^2=\frac{\tilde{g}_* }{g_*}H^2_\mathrm{SBBN} \quad \text {with}\quad H^2_\mathrm{SBBN}=\frac{8\pi G_0 }{3}\frac{\pi ^2}{30}g_{*}T^4. \end{aligned}$$

(78)

It is clear from (77) and (78) that we can set the following relation $\tilde{g}_*=\frac{G_1}{G_0} g_{*}$ implying that the altered G, i.e., $G_1$, in the BD gravity can equivalently be treated as modified $g_*$, namely, $\tilde{g}_*$, in GR. The effect of altered expansion rate of the Universe during BBN due to the modified $g_*$ within GR is well investigated in the literature [153,154,155] and is parametrized in terms of $S\equiv \frac{H}{H_\mathrm{SBBN}}=\sqrt{\frac{\tilde{g}_*}{g_*}}$, which can be adopted as $S=\sqrt{\frac{G_1}{G_0}}$ for our work within BD gravity by preserving $g_*$ as in the SBBN. Such that, deviations from $S=1$ (SBBN) will modify the neutron abundance and the time available for nuclear production/destruction, changing the BBN-predicted primordial element abundances. In general, it is necessary to access to a BBN code to study the BBN-predicted primordial abundances (viz., mass fractions) as functions of S and $\eta _{10}$ (the number density of baryons $n_\mathrm{b}$ to the number density of CMB photons $n_\mathrm{\gamma }$ defined as $\eta _{10}=10^{10} n_\mathrm{b}/n_\mathrm{\gamma }$ ). On the other hand, luckily, they have been identified (e.g., in [153,154,155]) by extremely simple but quite accurate analytic fits over a limited range in these variables as $5.5\lesssim \eta _{10}\lesssim 6.5$ and $0.85\lesssim S \lesssim 1.15$ and, for $^4$He (helium) and D (deuterium) these are¹⁰ [155]

$$\begin{aligned} Y_{\mathrm{p}}&=0.2381\pm 0.0006+0.0016[\eta _{10}+100 (S-1)], \end{aligned}$$

(79)

$$\begin{aligned} y_{\mathrm{DP}}&=45.7(1\pm 0.06)[\eta _{10}-6(S-1)]^{-1.6}, \end{aligned}$$

(80)

where, $Y_{\mathrm{p}}\equiv \frac{n_\mathrm{He}}{n_\mathrm{b}}$ and $y_{\mathrm{DP}}\equiv 10^5\frac{n_\mathrm{D}}{n_\mathrm{H}}$ are $^4$He and D mass fractions, respectively, and here in our study, $\eta _{10}=273.9\,\Omega _\mathrm{b,0} h^2$ and $S=\sqrt{\frac{G_1}{G_0}}\simeq (1+z_{\mathrm{eq}})^{\frac{1}{2(1+\omega )}}$, see (33). We see that $\eta _{10}\approx 6.1$ for both models from the constraints on $\Omega _\mathrm{b,0} h^2$ given in Table 1 and find that the cosmological data we considered constrain the altered expansion rate of the Universe for $z>z_{\mathrm{eq}}$, hence during BBN, as $1<S<1.0665$ at $68\% \mathrm{\, C.L.}$ and $1<S<1.0839$ at $95\% \mathrm{\, C.L.}$. These are well inside the validity intervals of (79) and (80) and hence they can be safely utilized here for obtaining the BBN-predicted $Y_{\mathrm{p}}$ and $y_{\mathrm{DP}}$ values, which in turn can be used for a further investigation of the constraint on the BD parameter $\omega$.

In the standard $\Lambda$CDM model [$S=1$ (GR), $\Omega _{\sigma ^2,0}=0$ and accommodating SBBN], considering the constraint on $\Omega _\mathrm{b,0}h^2$ given in Table 1, we obtain $\eta _{10}=6.1469 \pm 0.0404$ leading, from (79) and (80), to

$$\begin{aligned} Y_{\mathrm{P}}^\mathrm{SBBN}&=0.24793\pm 0.00065\quad (68\%\, \mathrm{C.L.}), \end{aligned}$$

(81)

$$\begin{aligned} y_{\mathrm{DP}}^\mathrm{SBBN}&=2.4999 \pm 0.1642\quad (68\%\, \mathrm{C.L.}). \end{aligned}$$

(82)

On the other hand, in the BD extension of the standard $\Lambda$CDM model, by assuming $\Omega _{\sigma ^2,0}\lesssim 10^{-23}$ [viz., expansion anisotropy is negligible during BBN, see (76)] and considering the constraints on $\Omega _\mathrm{b,0}h^2$, $\omega$ and $z_{\mathrm{eq}}$ given in Table 1, we obtain $\eta _{10}=6.1373\pm 0.0448$ and $1<S<1.0665$ at $68\% \mathrm{\, C.L.}$ leading, from (79) and (80), to

$$\begin{aligned} Y_{\mathrm{p}}&=0.24898\pm 0.00199\quad (68\%\, \mathrm{C.L.}), \end{aligned}$$

(83)

$$\begin{aligned} y_{\mathrm{DP}}&=2.5338\pm 0.1726\quad (68\%\, \mathrm{C.L.}). \end{aligned}$$

(84)

Fig. 6

The 2D constraints on the BBN predicted $^4$He (helium) mass fraction $Y_{\mathrm{p}}$ (top panel) and D (deuterium) $y_{\mathrm{DP}}$ mass fraction (bottom panel) in BD gravity are plotted with 1$\sigma$ (dark green) and 2$\sigma$ (light green) confidence contours. Grey bands are for the SBBN-predicted mass fractions with 1$\sigma$. Red bands are independent observational estimates with 1$\sigma$ from [157, 158]

Observations of helium and hydrogen recombination lines from metal-poor extragalactic H II regions provide an independent method (viz., a direct measurement) for determining the primordial helium abundance and a latest and widely accepted estimate comes from the data compilations of [157] giving $Y_{\mathrm{p}}=0.2449 \pm 0.0040$ (68 $\%$ C.L.). Similarly, the most recent estimate of the primordial deuterium abundance comes from the best seven measurements in metal-poor damped Lyman-$\alpha$ systems studied in [158] giving $y_{\mathrm{DP}}=2.527\pm 0.030$ (68 $\%$ C.L.). We note that the BBN-predicted $Y_{\mathrm{p}}$ and $y_{\mathrm{DP}}$ for both the $\Lambda$CDM model (81) and its anisotropic BD extension (83) obtained by using the cosmological data, led to consistent values with independent observational estimates. We notice, however, slightly larger mean values in the case of BD gravity, as suggested by (79) and (80) when $S>1$ (assuming $\eta _{10}$ is fixed). We give a summary of our findings by depicting the 2D marginalized posterior distribution of the BBN-predicted $Y_{\mathrm{P}}$ via (79)/$y_{\mathrm{DP}}$ via (80) and $\omega$ in the top panel/bottom panel of Fig. 6. In which, for a comparison, we depict also the bands of the SBBN-predicted $Y_{\mathrm{P}}$ via (79)/$y_{\mathrm{DP}}$ via (80) for the $\Lambda$CDM model accommodating SBBN ($S=1$) and of their independent observational estimates given in [157, 158]. We note that there is an increasing anti-correlation between $Y_{\mathrm{P}}$ and $\omega$ with decreasing $\omega$ as it approaches its lower bound, whereas for large values of $\omega$ the $Y_{\mathrm{P}}-\omega$ contour at 68% C.L. (dark green) approaches the band of the SBNN-predicted $Y_{\mathrm{P}}$ at 68% C.L. (grey) as it should be. Besides these, more importantly, we see that the $Y_{\mathrm{P}}-\omega$ contour at 68% C.L. (dark green) stays above the band of the independent observational estimate band (red) for $\omega \lesssim 250$, namely, its consistency with the independent observational estimate of $Y_{\mathrm{p}}=0.2449 \pm 0.0040$ from [157] requires $\omega \gtrsim 250$ at 68% C.L. providing us an improved constraint on the BD parameter $\omega$. We see that the BBN-predicted $y_{\mathrm{DP}}$ via (80) is insensitive to $\omega$, $y_{\mathrm{DP}}-\omega$ contour at 68% C.L. (dark green) already covers the independent observational estimate band (red) for $\omega \gtrsim 250$ (the improved constraint from our $Y_{\mathrm{P}}$ investigation) but is much wider than it (due to the relatively large internal error in (80), viz., $y_{\mathrm{DP}}\propto 1\pm 0.06$), and $y_{\mathrm{DP}}-\omega$ contour for BD gravity (dark green) approaches the SBBN-predicted $y_{\mathrm{DP}}$ band at 68% C.L. (grey) for large $\omega$ values as it should be. These show that we are not able to deduce a new constraint on the BD parameter $\omega$ by comparing the BBN-predicted value of $y_{\mathrm{DP}}$ via (80) with its independent observational estimate.

Thus, by means of the BBN, we reach the most stringent constraints on $\Omega _{\sigma ^2,0}$ and $\omega$, the two free parameters that determine the anisotropic BD extension to the standard $\Lambda$CDM model, as follows;

$$\begin{aligned} \Omega _{\sigma ^2,0}\lesssim 10^{-23}\quad \text {and}\quad \omega \gtrsim 250. \end{aligned}$$

(85)

Finally, considering these improved constraints, we find that the contribution to this constraint of $\Omega _{\sigma ^2,0}$ from the anisotropy of the effective pressure of the Jordan field during $z<z_{\mathrm{eq}}$, which is encoded in the term $(1+z_{\mathrm{eq}})^{\frac{1}{1+\omega }}$ in (75), is quite insignificant, such that, it leads to only a few percent smaller upper bound value for $\Omega _{\sigma ^2,0}$, viz., $(1+z_{\mathrm{eq}})^{-\frac{1}{1+\omega }}\lesssim 10^{-0.01}=0.977$ for $\omega \gtrsim 250$.

6 Conclusions

We have carried out 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 the two additional degrees of freedom: the so called BD parameter $\omega$ and the present day value of the density parameter corresponding to the shear scalar (a measure of the expansion anisotropy), $\Omega _{\sigma ^2,0}$. The role of the cosmological constant $\Lambda$ is taken over by the Jordan field potential of the form $U(\varphi )=\frac{1}{2}M^2\varphi ^2$ with M being the bare mass of the field. We have considered the LRS Bianchi type I metric, which generalizes the spatially flat RW metric simply by allowing a different scale factor along one of the three orthogonal axes, while preserving the spatial homogeneity and flatness, and the isotropic spatial curvature.¹¹ We have solved the field equations analytically and obtained the average Hubble parameter, H(z), explicitly by extending the method developed in [100]; described the anisotropic effective DE in accordance with the way of defining effective source given in [45, 46], and consistently included the radiation into the model.

The BD parameter $\omega$, being the measure of the deviation from GR ($|\omega |\rightarrow \infty$), by alone characterizes the dynamical behaviour of the effective DE as well as the redshift dependency of the expansion anisotropy. These two affect each other depending on $\omega$, such that the shear scalar contributes to the dynamics of the effective DE, and its anisotropic stress controls the dynamics of the shear scalar, in particular deviations from its usual form $\sigma ^2\propto (1+z)^6$ in GR, see H(z) in (34). We planned the current study as an extension in the sense of a correction to the standard $\Lambda$CDM model, so we have mainly confined our investigations to small deviations from this model via sufficiently small $\Omega _{\sigma ^2,0}$ and large $\omega$ values. We have shown through some preliminary cosmological discussions that $|\omega |\lesssim 10$ (roughly) cannot be called as a small deviation, yet, for completeness we have extended our investigations to non-negative values of $\omega$. Indeed, considering the combined data sets PLK+BAO+SN+H, we have obtained $\omega \gtrsim 50$, $M<1.51\times 10^{-34}\mathrm eV$ and $\Omega _{\sigma ^2,0}\lesssim 10^{-8.48}$ ($\Omega _{\sigma ^2,0}\lesssim 10^{-14}$, when matter-radiation equality is considered).¹² Then by means of BBN, we have improved these constraints to $\omega \gtrsim 250$ (in particular, from the comparison of the helium abundance prediction of the model with the direct measurements) and $\Omega _{\sigma ^2,0}\lesssim 10^{-23}$. We have also found that the contribution of the anisotropy of the effective DE (viz., $\Delta w_{\mathrm{DE,0}}<4.23\times 10^{-7}$) to this constraint on $\Omega _{\sigma ^2,0}$ is insignificant, namely, led to only a few percent smaller (stronger) upper bound value.¹³

All these have led us to conclude that, with the observations relevant to the dynamics of the Universe, the simplest anisotropic massive BD gravity extension of the standard $\Lambda$CDM model presents no significant deviations from it all the way to the BBN. The strongest cosmological constraints on $\omega$ present in the literature, e.g., $\omega \gtrsim 890$ [116], would obviously just strengthen this conclusion. Moreover, we should further consider the strongest local constraint $\omega >40000$ [114] as well since the constraint $M<1.51\times 10^{-34}\,\mathrm eV$ obtained here for the mass of the Jordan field leads to no relaxation on it [115]. This in turn implies that, against the local constraints, the cosmological features that arise from replacing GR by massive BD (such as the modified expansion anisotropy due to the anisotropic effective DE) are not significantly sensitive to the current cosmological observations. In other words, in view of the local constraints on $\omega$, the simplest anisotropic massive BD gravity extension of the standard $\Lambda$CDM model cannot be distinguished from its simplest GR based anisotropic extension [95] when we consider the data from cosmological observations.¹⁴

It is conceivable that our findings, when we consider only the observations relevant to the background dynamics of the Universe, are representative for the ones that would be obtained in the more general extensions, namely, through a gravity theory more general than the massive BD and therefore can evade the local constraints. For instance, promoting $\omega$ to be some functions of $\varphi$, the observations from the Cassini spacecraft place upon strong constraint on $\omega$ as in the original BD theory. This constraint, however, now only applies to the local value of $\omega (\varphi )$, i.e., of the present day value of $\varphi$ in the Solar System. Similarly, it is possible that $\omega$ has spatial variation, so that, for instance, it could be constrained to be $\omega >40000$ locally, but could be in line with $\omega \sim 10^2$ at cosmological scales [159, 160]. See, for instance, [37] and references therein for a further reading on such gravity theories closely related to the BD theory, of which the Horndeski/Galileon [124,125,126] theories are the most popular ones. In the anisotropic extensions through such gravity theories, however, one should deal with more complicated field equations, those may not be solved analytically and lead to the explicit expression of H(z), and hence the solutions and findings, presented in this work, are useful as they may be considered as the good approximations to such constructions. Consequently, the theoretical investigations we carried out here are, in general, instructive for the anisotropic extensions of the standard $\Lambda$CDM model replacing GR by a modified theory of gravity approximating massive BD at cosmological scales. When the observational constraints are considered as well, our findings against the cosmological constraints (the local constraints) with $\omega \gtrsim 10^{2}$ ($\omega \gtrsim 40000$) are instructive for the extensions through those gravity theories that can (cannot) evade local constraints. Yet, this also implies that the lesson learned from this study based on the massive BD would approximately be valid for such more general constructions, namely, many interesting features of such anisotropic extensions of the $\Lambda$CDM model that do not exist within its simple GR based anisotropic extension, would remain insignificant when the model is constrained by the observations.

We note that, in fact, these features become quite significant for $|\omega |\sim O(1)$, which, interestingly, is in principle the natural order of magnitude for the BD parameter (as, e.g., $\omega =-1$ appears in the low energy limit of string theories). This by alone, makes it quite interesting, at least theoretically, to further study this region in spite of that, as we have discussed above, it implies large deviations from the standard $\Lambda$CDM dynamics, and hence is not expected to be consistent with the real Universe. Such a study, however, is beyond the scope of the current paper, but yet, here, we would like to make a couple of relevant comments. In case of small negative values of $\omega$, say, $\omega \sim -1$ but larger than the scale invariant limit $-\frac{3}{2}$ (below which stability issues appear) the model exhibits several critical points that complicates the investigation in this region. So it is necessary to carry out several separate analyses within this small region. Meanwhile, this region is distinguished that it is relevant to string theories, namely, $\omega =-1$ corresponds to the low energy effective string action and $\omega \sim -1$ appears in d-brane constructions [102,103,104,105]. For these reasons in fact, although we have devoted the current work to non-negative $\omega$ values, we carried out a discussion in Sect. 3.3.2 showing how dramatically the redshift dependence of the expansion anisotropy can alter when $\omega \sim -1$. For instance, for $\omega =-\frac{4}{3}$, expansion anisotropy becomes non-dynamical (mimicking $\Lambda$) and for $\omega =-1$, expansion anisotropy should be set to zero. The region $-\frac{3}{2}<\omega <-\frac{4}{3}$ is also interesting that, in the presence of only dust (without $\Lambda$ or the mass of Jordan field), it is the region of BD theory leading to accelerated expansion of the Universe, and we notice that the expansion anisotropy also contributes to this acceleration since it also behaves like DE in this region. We do not know, so far, a concrete and successful mechanism making such small values (negative or positive) of $\omega$ consistent with local experiments and cosmological observations. Yet, involving extra spatial dimensions may be good place for looking for such possibilities as it was shown that it is possible to make BD theory (massive or massless) consistent with the gravitational tests (including solar system tests) for $|\omega |=O(1)$ in the presence of extra dimensions [127]. Noticing the presence of extra dimensions as part of string theories and the dynamical extra dimensions (internal space) contribute to the evolution of the four dimensional anisotropic spacetime like an anisotropic stress in a similar way that the Jordan field does in BD theory make this option more appealing and interesting [161,162,163,164,165]. This shows one of the many possibilities of rich behaviors which could be worthwhile to investigate in future studies on the extensions of the standard $\Lambda$CDM model wherein the expansion anisotropy exhibits non-trivial behaviors that may have interesting cosmological consequences.

A more rigorous observational investigation of the model may be another direction to follow. We have studied the observational constraints on the parameters of the model by considering the background dynamics of the Universe through the evolution of the comoving volume scale factor. These constraints can be improved when we consider the perturbation sector – as in [116,117,118] providing the most stringent cosmological constraints on $\omega$– along with, for instance, the full Planck data – as in [93, 94] providing constraints on the anisotropic expansion on the top of the $\Lambda$CDM at a level close to the ones from BBN by considering the CMB radiation temperature and polarization data-. Such a further observational investigation of the model requires considerable amount of work beyond the aim of the current paper, as it is necessary to study perturbations on the anisotropic background in the presence of an effective anisotropic source. This would strengthen the constraints on the parameters of the model, but would probably not change our conclusion that the model exhibits no significant deviations from the standard $\Lambda$CDM model all the way to the BBN. Yet, we cannot be sure about that unless we undertake this work in future and see the results.