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

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


INTRODUCTION
The base Λ-cold dark matter (Λ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 Λ 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 Λ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 Λ and also can account for the small, but significant, deviations from the ΛCDM model as persistently suggested by the high * akarsuo@itu.edu.tr † nihan.katirci@itu.edu.tr ‡ nozdemir@itu.edu.tr § javazquez@icf.unam.mx 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 Λ 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 ≈ −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-α 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 theo-ries, 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 ΛCDM model implies that the late time acceleration is not driven by Λ 1 and of DE yielding w < −1 and/or ρ < 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 2 , which can be realized if we relax the isotropic space assumption of the Λ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 Λ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 Λ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., (∆T /T ) =2 ∼ 10 −5 . And, in the simplest anisotropic generalization of the standard Λ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, ρ σ 2 ∝ 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 Ω σ 2 ,0 ∼ 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, Ω σ 2 ,0 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 Ω σ 2 ,0 10 −4 for z ∼ 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 ω. 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 σ 2 ≤ 5 3 , −1.14 w < −1 and 0.50 ≤ z PDL < 0.65 (see [100] for a similar result) for positive values of the BD parameter ω. 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 ω = −1 and similarly d-brane models yield BD gravity with ω ∼ −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 ω, determining the deviations from GR, which is recovered in the limit of |ω| → ∞.
Theoretically ω ≥ − 3 2 is necessary to avoid the Jordan/scalar field (ϕ) 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, ω 40000 at the 2σ 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 (ϕ), as in this case one should consider the corresponding mass (M ) of the Jordan field along with the BD parameter, viz., the {ω, M } parameter region [115]. It is shown that ω = O(1) is allowed for the Solar System constraints provided that M 2 × 10 −17 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 H0 10 −33 eV, which leaves the constraint from the Cassini spacecraft still valid.
The cosmological constraints on ω 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 (ϕ) = const. and considers CMB data from the first Planck release along with the constraints from Big Bang Nucleosynthesis (BBN) light element abundances, gives ω > 890 at the 2σ C.L. [116]. The constraints from Planck 2015 (2013) and a compilation of baryon acoustic oscillations (BAO) data are ω > 333 (ω > 208) at the 2σ C.L., weakly dependent on the exponent n > 0 for a power-law potential U (ϕ) ∝ ϕ n mimicking Λ in the late universe [117,118] (also see [119] for a recent work considering a more complicated potential). The BBN leads to ω 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 ω at a level comparable to the local tests [121][122][123].
The typical mass scale of M ∼ 10 −33 eV imposed by cosmological observations is so small that the local constraint ω 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 Λ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 ω 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, ω = ω(ϕ), in addition to the presence of potential U (ϕ) and this leads to the possibility of having an ω 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 secondorder 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 |ω| = 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 ΛCDM model, which is characterized by two additional free parameters, namely, the BD parameter ω and the corresponding present day density parameter of the expansion anisotropy Ω σ 2 ,0 . The role of the cosmological constant is taken over by the Jordan field potential of the form U (ϕ) ∝ ϕ 2 and the physical ingredient of the Universe is considered as in the standard Λ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 ΛCDM.

ANISOTROPIC MASSIVE BRANS-DICKE GRAVITY EXTENSION OF THE ΛCDM MODEL
We consider the BD action [106,110] written in the Jordan frame in the following form: where ϕ = ϕ(t) is the Jordan scalar field (function of cosmic time t only) and ω = const. is the Brans-Dicke parameter, R is the Ricci scalar, g is the determinant of the metric g µν , and S Matter is the matter action, which is independent of ϕ 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. 3 We assume M 2 = const. so that, as can also be seen from the action, it stands like a cosmological constant as 2ωM 2 ≡ Λ 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 ΛCDM-type cosmologies, such that the mass of the Jordan field alone can play the role of positive cosmological constant like in the standard ΛCDM cosmology provided that 2ωM 2 ≡ Λ > 0, and the Jordan field ϕ varying slowly enough on the top of this can account for small deviations from the standard ΛCDM model in a particular way, which in turn may lead to an improved fit to the observational data w.r.t. the standard ΛCDM model. In line with that, we intend to study the BD extension of the standard ΛCDM model as a correction, and therefore we demand the term 2ωM 2 to be positive definite, which requires ω > 0 as long as we keep M 2 > 0 to avoid the Jordan field from having an imaginary mass. 4 Hence, in this study, unless otherwise 3 We consider the bare mass as it is done in [100], where it formally appears as the mass of a minimally coupled canonical scalar field when curvature scalar is dropped. On the other hand, when the effective mass of the Jordan field is considered there are various definitions in the literature (see, e.g., [130] for a recent discussion and further references), though one cannot say that these would not fail from a strict particle physics point of view. 4 If this is not the case then 2ωM 2 = Λ will contribute to the field equations like a negative cosmological constant, which may be compensated by a rapidly changing Jordan field. Within the standard BD gravity (i.e., massless BD gravity) in the presence of pressureless source, there exist cosmological solutions with accelerating expansion if −2 < ω < −1, and Jordan field is real (i.e., the effective cosmological gravitational coupling is positive definite) as well if − 3 2 < ω < − 4 3 leading to deceleration parameter in the range 0 > q > −1, and to ϕ 2 ∝ (1 + z) 2 and ϕ 2 ∝ (1 + z) 3 in the two boundaries of this range, correspond-is mentioned, we carry out our investigations by assuming ω > 0, which is already stronger than the assumption ω ≥ −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; where S = S(t) is the mean scale factor. Here, the exponent β = β(t) satisfies the relationβ 2 = 1 6 σ 2 , where σ 2 = σ ij σ ij and σ 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 S S +Ṡ 2 S 2 − σ 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) where ρ 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: ingly [131]. Hence, obviously, if we consider ω < 0 and M 2 > 0 leading to 2ωM 2 < 0 (i.e., effectively negative cosmological constant), then it will be necessary to confine ourselves to the range − 3 2 < ω < − 4 3 , and moreover the decelerating effect of the term 2ωM 2 < 0 would be compensated by bringing the value of the BD parameter closer to −4/3, which in turn, implies a faster Jordan field [viz., for ω ∼ −4/3 Jordan field changes as fast as the pressureless source, i.e., ϕ 2 ∼ (1+z) 3 ] whereas we are looking for slowly changing Jordan field, say, ϕ 2 ∼ const. since we demand it to do only small modifications on the ΛCDM dynamics.
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 We consider the standard cosmological fluids, namely, pressureless fluid (p m = 0) and radiation/relativistic fluids (p r = ρ r /3), so that, here, ρ = ρ m + ρ r and p = p r . If the Jordan field is approximately constant then the first Friedmann equation (4) becomes 3Ṡ 2 S 2 ≈ 4 ϕ 2 ρ+ σ 2 2 +2ωM 2 (as like in GR), where the effective cosmological gravitational coupling strength is then given by G = with zero subscript denoting today's value (throughout the paper). 5 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 ) 8πG 0 = 4/ϕ 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 µν,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: where ρ DE is the energy density and p DE,x and p DE,y are the principal pressures of the effective DE along the xaxis and the y-and z-axes, respectively, and which read 5 Note however that for bound systems in the quasi-static regime, e.g., our solar system, the effective Newton's constant is G N = 2ω+3 2ω+4 G 0 [110], thus observers in a bound system which formed today would measure the same cosmological and local gravita- 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 This, in turn, leads to a modification in the evolution of the expansion anisotropy provided that ϕ 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 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]).

EXACT SOLUTION COMPATIBLE WITH STANDARD ΛCDM MODEL
We follow the method, given in [100], for obtaining exact cosmological solutions of a scalar-tensor gravity theory compatible with the ΛCDM model, albeit we consider LRS Bianchi type I spacetime rather than spatially flat Robertson-Walker spacetime. 6 7 To do so, we first re-express the set of differential equations given by Eqs. 6 We note that in the GR limit (say, ϕ = const.) the anisotropic model under consideration here, (4)- (7), reduces mathematically to the standard ΛCDM+stiff matter model in the RW framework [133]. However, the role of the stiff matter with positive energy density in [133] is played by the shear scalar here, so that one can straightforwardly utilize the rich class of solutions given in [133] for the GR limit in our model. This presents a good example that one could find various solutions of the model under consideration here but yet we focus on the solution obtained by extending the method given in [100] to LRS Bianchi type I metric. 7 We followed the method given in [100] for obtaining an exact cosmological solution of BD gravity theory for the expansion rate H(z), with its ΛCDM model counterpart up to a large redshift, (4)-(7) using dz dt = −H(1 + z), where H =Ṡ/S is the average Hubble parameter and z is the cosmic redshift we define in terms of S as z = −1 + 1 S . Accordingly, we re-express the energy density equation (4) as and obtain the pressure equation as using Eqs. (5) and (6), the shear propagation equation as by subtracting Eq. (5) from Eq. (6), and finally reexpress the scalar field equation (7) as where denotes derivative with respect to redshift (d/dz). 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 provided that we set p = p m = 0, viz., the universe is filled with only pressureless matter. We next note that constants D 1 = D 2 = −2ωM 2 , and that the shear scalar σ 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) = σ 2 2 . In accordance with these points, we assume viz., for pressureless matter in a spatially flat, homogeneous and isotropic universe. A more general exact solution of the same setup was given in [134] (much earlier than [100]), was nevertheless given for the scale factor in cosmic time and very complicated for extracting an exact H(z) required for observational analyses. (18) and (20) and then solve for the rate of change of Jordan field in z as 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 where ϕ 0 = ϕ(z = 0) is the present time value of the Jordan field, and which in turn gives The integration of the shear propagation equation given in Eq. (19) gives where σ 2 0 = σ 2 (z = 0) is the present time value of the shear scalar, leading, together with (23), to whereas it is σ 2 GR ∝ (1 + z) 6 in GR (corresponding to the |ω| → ∞ limit of the BD gravity), and thereby the energy density corresponding to the shear scalar can be defined as 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 which is exactly the same as the one for the pressureless matter in the standard ΛCDM model. Finally using ϕ 2 , σ 2 and ρ 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 and the energy densities corresponding to the mass of the Jordan field and to the expansion anisotropy today read, respectively, One may check for consistency that we recover the H(z) of the standard ΛCDM model when we set ω → ∞ (viz., in the GR limit of BD gravity giving ϕ → constant) and ρ σ 2 ,0 = 0 (viz., isotropic expansion as in the RW spacetime metric). We also note that ρ M can be regularized to give a finite positive value consistent with the observations by setting M 2 ω = constant so that M 2 → 0 as ω → ∞ such that the term M 2 ω arising from their multiplication would always remain unaltered and finite.

Inclusion of radiation
The solution of the BD extension of standard Λ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, 4 ϕ 2 ρ m ∝ (1+z) 3+ 1 1+ω , is dominant over the effective cosmological constant 2ωM 2 and one may check that for this epoch, neglecting anisotropy, it reduces to the well known solution with S ∝ t 2+2ω 4+3ω and ϕ 2 ∝ t 2 4+3ω (see, e.g., [120] and references therein). We could not obtain analytically an exact cosmological solution -other than the trivial solution with ϕ = const.-using the same method when we include relativistic source, p r = ρ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 ρ r ∝ (1 + z) 4 and ρ m ∝ (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 eq = −1 + ρm,0 ρr,0 ∼ 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., ϕ = 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 eq , we have ϕ = ϕ 1 = const., implying G = G 1 = const. as long as anisotropy is negligible. Note that ϕ 1 (or G 1 ) will be different than ϕ 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 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 eq the Jordan field evolves approximately in accordance with (23) and thereby we can estimate that which implies where G 1 is then approximately the cosmological gravitational coupling strength for z > z eq (see [120] for further details). Thus, for ω > 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., . 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 σ 2 ∝ σ 2 0 (1 + z) 6 like in GR in the presence of isotropic perfect fluid, since we have ϕ ∼ 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 z BBN , expansion anisotropy will dominate over the terms due to the radiation and hence we obtain once again ρ σ 2 ∝ (1 + z) 6+ 2 1+ω (which can be obtained simply by setting ρ M = 0 and ρ σ 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, ρ σ 2 ,0 ρ r,0 ρ m,0 ∼ ρ M for ω 10 (see Sect. 3.3.1 for this choice on ω)] for z < z eq . Thus, we can write 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 ∼ 3380, and the recombination that leads to photon decoupling should have taken place when z ∼ 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 ϕ constant relying on its attractor behavior during the radiation-dominated phase, we can write for z eq < z < z eq,σ 2 ,r .
Note that, for consistency between (34) and (35), here in (35) we have ϕ 1 ∼ ϕ 0 (1 + z eq ) − 1 2+2ω in contrast to ϕ 0 in (34), in accordance with the above discussions leading to (32), and similarly ρ σ 2 ,1 ∼ ρ σ 2 ,0 (1 + z eq ) 1 1+ω . Using (35), we estimate (at least roughly) the radiation-expansion anisotropy equality redshift as z eq,σ 2 ,r ∼ −1 + ρr,0 ρ σ 2 ,1 , which can also be written as follows; Avoiding expansion anisotropy from spoiling the physical processes relevant to BBN will then roughly require z eq,σ 2 ,r > z BBN . Typically, in a viable cosmological model, we expect the matter-radiation equality redshift to be z eq = −1 + ρm,0 ρr,0 ∼ 3380, the BBN to take place at z BBN ∼ 3 × 10 8 and Ω r,0 ∼ 10 −4 . Accordingly, taking z eq,σ 2 ,r > 3 × 10 8 , we obtain 1+ω , namely, the following stringent constraint on the density parameter of the expansion anisotropy today: which reduces to Ω σ 2 ,0 ∼ 10 −21 in the GR limit (ω → ∞). 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 ω > 0 we consider in this study, BD theory leads to more stringent constraints on the expansion anisotropy (e.g., Ω σ 2 ,0 10 −25 for ω = 0), while the case ω < −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 1.6, as Ω σ 2 ,0 10 −4 (see also [97] for a recent work). We note that this value is achieved upon choosing ω = − 21 17 in (37), which can be done without worrying for theoretical issues associated with ω < − 3 2 , but this will cost large deviations from GR (viz. rapidly varying strength of the gravitational coupling) as well as the standard ΛCDM. We shall further discuss possible modifications on the expansion anisotropy for different values/ranges of ω and their comparisons with the standard Λ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 from our solution (see, e.g., Refs. [132,[139][140][141] for anisotropy dominated universes in the context of scalartensor theories). This epoch however is out of interest in the present work.

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+Λ 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 ϕ freezes, implying the BD gravity mimics GR, and thereby the effective DE will mimic nothing but Λ (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 (|ω| → ∞). (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 Λ, 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: where the energy density of the effective DE, ρ DE , reads which gives for z = 0. Dividing (39) by 3H 2 0 , we obtain where we define f (z) = ρDE ρDE,0 , reading which gives f (z = 0) = 1 consistently with Ω m,0 +Ω DE,0 + Ω σ 2 ,0 = 1. Note that we define the present time values of the density parameters as Ω i,0 = ρ i,0 /ρ c,0 , where ρ c,0 = 3H 2 0 ϕ 2 0 /4 is the present time value of the critical energy density, and that we eliminate Ω M,0 , that appears in f (z) originally, using the relation which can easily be obtained from (40). Using (44) along with Ω m,0 + Ω σ 2 ,0 = 1 − Ω 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 ω [viz., γ(ω)] as 8 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, which give for z = 0. The anisotropy of the EoS of the effective DE may be given as the difference between its principal EoS param- (50) Accordingly, we have which is approximately equal to 1 1+ω Ω σ 2 ,0 ΩDE,0 , provided that the effective DE is dominant today (z = 0). Additionally, we can write for dust domination which is positive definite as ω > 0, and for expansion anisotropy domination which depends only on ω, is negative definite as ω > 0, monotonic in ω > 0 and takes values as We finally, usingρ DE + 3Hρ DE (1 +w DE ) = 0, define a mean/volumetric effective EoS parameter for the effective DE as which explicitly reads and for z = 0. We note that, given Ω DE,0 > 0 and Ω m,0 and Ω σ 2 ,0 are already positive definite, for ω > 0, the contributions from both the pressureless matter and the expansion anisotropy to the present value of the mean EoS of the effective DE,w DE,0 , are negative, i.e., anisotropic BD extension of the standard ΛCDM model predicts phantom like effective DE in the present time Universe.

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 ω 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 + H H (1 + z). This reads, as we expect in a viable cosmology Ω σ 2 ,0 < Ω r,0 Ω m,0 , for z ∼ 0, and for z = 0 today. We confine the present study to small deviations from the standard ΛCDM model and hence it is reasonable to expect, in accordance with the most recent Planck results [12], that Ω m,0 ∼ 0.3 leading to q 0 ∼ −0.55 from q = −1+ 3 2 Ω m in the standard ΛCDM, which in turn implies from (58) that we should have |ω| 10, which leads to ω 10 since have already assumed ω ≥ 0. 9 Unless otherwise, for instance, if ω = 0, then q 0 ∼ −0.4, which could be decreased to a reasonable value, viz., q 0 ∼ −0.55, by decreasing Ω m,0 considerably, namely Ω m,0 ∼ 0.22, which implies very large deviations from the standard ΛCDM. Therefore, our assumption ω ≥ 0 already allows large deviations from the standard Λ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 ω, namely, when ω ∼ −1. Hence, in this section, we shall depict some key/interesting features of the anisotropic BD model and its deviation from the standard ΛCDM model assuming ω ≥ 0.
First of all, one may check easily that, similar to the standard ΛCDM, its anisotropic BD extension asymptotically approaches to de Sitter Universe in the far future, and accordingly, at this limit, the effective DE mimics exactly the cosmological constant as At low redshifts, viz., at which ρ m ρ m,0 and ρ σ 2 ρ m , the effective DE can approximately be described by Here, we first note that the non-negativity of the effective DE today (z = 0), ρ DE,0 ≥ 0, leads to a natural lower bound on the value of the energy density corresponding to the effective cosmological constant (viz., ρ M ) as ρ M ≥ 1 γ − 1 ρ m,0 , which allows zero lower bound only at the GR limit, i.e., at ω → ∞ (viz., γ → 1). We next note that γ takes values within the range 1 2 ≤ γ < 1 for ω ≥ 0 and hence the coefficient of ρ m is negative definite, γ − 1 < 0, which implies that, because ρ m decreases as the Universe expands, ρ DE will increase as the Universe expands at z ∼ 0, i.e., the effective DE will behave like a phantom fluid (w DE −1) at z ∼ 0 [see (56) for w DE (z = 0)]. Namely, we see from (40) or (55) that given ρ m = ρ m,0 (1 + z) 3 , the effective DE passes below the phantom-divide-line (PDL, w DE = −1) at and afterwards stays there forever, i.e.,w DE < −1 for z < z PDL . It is worth noting that the BD gravity predicts almost a specific redshift for the PDL crossing, such that and 0.63 ≤ z PDL < 0.65 for ω ≥ 10 (see [100] for a similar result). This is indeed an interesting result, such that z PDL ∼ 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 ∼ 0.6 with high confidence would imply a strong reason for favoring the BD gravity over GR, or vice versa. During this period, ∆w DE can approximately be given by (51) and we note that the effective DE will be almost isotropic, ∆w DE ≈ 0, for ω > 0 since Ω σ 2 ,0 Ω DE,0 in a realistic cosmology. At moderate redshifts (e.g., z PDL z z 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 This plateau is persistent all the way to the large redshift values at which either the radiation becomes dominant over pressureless matter (z ∼ z eq ), or the contribution from the expansion anisotropy becomes significant. The plateau ofw DE is located within the range 0 ≤w DE ≤ 1 3 depending on the value of ω ≥ 0. For ω = 0, the effective DE behaves like an extra relativistic degree of freedom asw DE = 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 ΛCDM model only slightly), we expect ω 10 and this places the plateau in a value in the range 0 <w DE 0.03. During this period, the evolution of the anisotropy of the effective DE can approximately be given by (52). We note that, although ∆w DE increases as (1 + z) 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; We note that, except the GR limit ω → ∞ (i.e., γ → 1), the mean/volumetric EoS parameter of the effective DE yields a second plateau placing in the range 1 w DE 5 3 for ω ≥ 0, and the effective DE behaves like the expansion anisotropy with a value of magnitude of ratio in the range 0 < |ρ DE /ρ σ 2 | 1 2 but yields negative energy density since the coefficient γ −1 is negative for ω ≥ 0. According to this the energy density of the effective DE changes sign and we calculate from (40) that this occurs at Ωm,0 (1 + z) 3+ 2 1+ω , which may be read off from (42), we obtain We note that, in a realistic cosmological setup, say, ω 10 (even for ω ≥ 0) and Ω m,0 Ω σ 2 ,0 , this ratio is obviously much larger than unity, implying that z 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 ϕ freezes out) so that the effective DE will be constant (i.e., mimics cosmological constant) while the expansion anisotropy keeps on growing as ρ σ 2 ∝ (1 + z) 6 as z increases. Eventually, the expansion anisotropy will dominate over radiation and the dynamical effective DE -viz., ρ DE ∝ ρ σ 2 (γ − 1) as may be seen from (40) by setting ρ M = 0 = ρ 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.

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 ρ σ 2 ∝ (1+z) 6 and thereby give Ω σ 2 ,0 10 −21 . On the other hand, the direct/model independent observational constraints from, e.g., type Ia SNe observations are rather weak, e.g., Ω σ 2 10 −4 for z ∼ 0 [96,97]. Such large values may be possible provided that the redshift dependence of ρ σ 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 ω ≥ 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 ω < −1 (since the expansion anisotropy in this case yields flatter redshift dependency) and even vastly for − 5 3 < ω < −1, which however would cost to large deviations from GR, hence from Λ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 Λ [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 ω < −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 We first note that it approaches Zeldovich fluid, i.e., w σ 2 → 1, at the GR limit |ω| → ∞ and, for ω ≥ 0 (finite), is in the range 1 < w σ 2 ≤ 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 ω > −1 as can be seen from (27). We depict w σ 2 versus ω in Fig. 1.
For ω ≤ − 4 3 , the EoS parameter spanning the range −1 ≤ w σ 2 < 1, which is familiar from the conventional cosmology: (i) ω = −2 (w σ 2 = 1 3 ), the expansion anisotropy mimics radiation, viz., plays the role of an extra relativistic species. mimics Λ > 0. It might be interesting to remind here that − 3 2 < ω < − 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 ω ∼ − 4 3 appears in d-branes string models [104,105] and hence it is conceivable that such models would predict σ 2 ∼ const. For − 4 3 < ω < −1, we have w σ 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 ω = −1 (the low energy effective string action limit of BD gravity [104,105]) is interesting that w σ 2 exhibits a pole, namely, lim ω → −1 + w σ 2 = +∞ and lim ω → −1 − w σ 2 = −∞, which imply that we must set σ 2 = 0 in this case (at least in our solution). In relevance with this, in the case of ω > −1 but ω ≈ −1, w σ 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 ω = 0 -the lowest value we consider in this study-we have ω σ 2 = 5 3 and for ω 10, which, as we discussed above, is required for small deviations from the ΛCDM model (by keeping isotropic RW metric), we have 1 < w σ 2 1.06. These imply that considering BD gravity with w 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 ΛCDM model in contrast to DE models that can be described by scalar field/s within GR.

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 Ω m,0 , the baryon density parameter today Ω b,0 h 2 , and the dimensionless Hubble constant h as well as the two free parameters ω and Ω σ 2 ,0 that account for the anisotropic BD extension of the ΛCDM model. We consider H(z) given in (34) rather than (39), since the former one contains radiation (viz., Ω 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 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 σ 2 , which contributes to the H(z) directly via the term ρ σ 2 ,0 (1 + z) 6+ 2 1+ω and also indirectly via its effect on the ρ DE as a result of BD theory).
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 I 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 eff = 3.046) with minimum allowed mass m ν = 0.06 eV theoretically well determined within the standard model of particle physics and the density parameter of radiation Ω r,0 = 2.469 ×  I: Constraints on the anisotropic Brans-Dicke extension of the standard ΛCDM model using the combined data sets PLK+BAO+SN+H. For two-tailed distributions, the results are given in 1σ and for one-tailed distributions, the results are given in 2σ. ΛCDM column corresponds to GR limit (Brans-Dicke gravity with ω → ∞) and the shear scalar σ 2 = 0 (isotropic). Parameters and ranges of the uniform priors assumed in our analysis and derived parameters are labeled with * .

Parameter
Anisotropic . The photon energy density today ρ γ,0 is not subject to our analysis since it is well constrained, such that it has a simple ρ γ = π 2 15 T 4 CMB relation with the CMB monopole temperature [148], which is very precisely measured to be T CMB,0 = 2.7255±0.0006 K [149]. Table I 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 ΛCDM model. Fig. 2 summarizes the constraints on the parameters that characterize the anisotropic BD extension to the standard ΛCDM model. Left panel of this figure displays 2D marginalized posterior distributions on the BD parameter ω (deviation from GR) and the density parameter corresponding to the expansion anisotropy today Ω σ 2 ,0 (deviation from isotropic expansion). We see from Table I that Table I. The extra two parameters from the anisotropic BD extension to the ΛCDM model may shift and relax some of the constraints, as it is the case of Ω m,0 = 0.3002 ± 0.0067, and therefore as a consequence the earlier redshift at matterradiation equality z eq = 3368.62 ± 30.89, compared to the ΛCDM model with z eq = 3359.45 ± 24.14 (see right panel of Fig. 3).
We depict, in Fig. 4 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 forw DE and from redshift z = 10 to z = 10 6 for ∆w DE and Ω DE . Though, note that the extensions to the redshift values z ∼ z 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σ and 2σ confidence intervals are plotted as black lines and red lines correspond to the best-fit values found over the analysis. At low redshifts, the upperleft panel shows that the crossing of the PDL occurs at z PDL = 0.64839 ± 0.00058, just as we concluded from the discussion that follows equation (62). For redshift values z < z PDL , the effective DE exhibits phantom behaviour w DE < −1 as may be seen from (56) or (61)   (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 ∆w DE increases only about an order of magnitude from its current value ∆w DE,0 10 −7 (see Table I) 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 Ω DE ∼ 0.1 for z ∼ 1. 75 and Ω DE ∼ 0 for z 10),w DE starts to noticeably climb up and settles in the first plateau ofw DE ∼ 0 (see (64) and paragraph covering it) lying between z ∼ 50 and z ∼ z eq (viz., throughout the pressureless matter dominated epoch). There is a period in this plateau during whichw DE > 0 [see (64)], viz., Ω 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, Ω DE during this period can never grow up to considerable values, viz., remains less than a percent at 1σ C.L. and few percents at 2σ 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., ∆w DE 0.07 at photon decoupling redshift z ∼ 1100) until the effective DE starts to leave this plateau as the expansion anisotropy starts to become dominant. We see thatw DE exhibits a pole at log 10 z DE,pole = 4.80 ± 0.58 then settles in a new plateau ofw DE ∼ 1 starting at z ∼ 10 5 during which ∆w DE exhibits a similar behavior and settles down at ∆w DE ∼ −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 ∼ z eq is basically an artifact of omitting radiation in Sect. 3.3.1 to have explicit expression of ρ 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 eq and be maintained all the way to the redshift values at which BBN took place z ∼ z 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 Λ (viz.,w DE = −1 and ∆w DE = 0 implying that it is isotropic and maintains its energy density value at z ∼ z eq for z > z eq ) and expansion anisotropy increases as ∝ (1 + z) 6 . However, although the data constrain z DE,pole to be about 1-2 orders of magnitude larger than z 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 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 eq , using the constraint on Ω σ 2 ,0 from Table I we calculate that the expansion anisotropy dominates over radiation at redshift values smaller than z 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 Ω σ 2 ,0 nor ω. These two parameters are not independent and hence we should further investigate the upper boundary on Ω σ 2 ,0 and lower boundary on ω, of course, particularly, by using the data providing information about the Universe at z z eq , such as the peak of the matter-power spectrum relevant to z = z eq and BBN relevant to z ∼ 10 8 .

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), readṡ which is negative definite considering ω ≥ 0 in this study, provided that the Universe is expanding (H > 0). The constraints we found on ω and H 0 (see Table I) predict Ġ G < 1.092 × 10 −12 yr −1 (68% C.L.) for z = 0, (70) similar to those given for z ∼ 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 |Ġ/G| 10 −8 yr −1 for z ∼ z eq . On the other hand -given that ϕ = const.
(and hence G = const.) is an attractor solution when radiation is dominant-the redshift dependence of G becomes flatter w.r.t. G ∝ (1+z) 1 1+ω [see (24)] as the radiation becomes more significant w.r.t. pressureless matter and thereby |Ġ/G| 10 −8 yr −1 for z ∼ z eq will never be achieved, but instead G will become almost constant for z ∼ z eq and remain so for z > z 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 eq and leads to G 1 ∼ G 0 (1 + z eq ) 1 1+ω [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 ∼ z eq w.r.t. its present time value as ∆G G0 | z∼zeq = 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: for z z eq . These constraints, of course, are valid also for the epoch of primordial nucleosynthesis that takes place when z BBN ∼ 3 × 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 ∆G/G 0 ).

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, Ω σ 2 ,0 = 10 −8.48 (see Table I), leads to a lower bound on the expansion anisotropy-radiation equality (Ω σ 2 = Ω r ) redshift as z eq,σ 2 ,r ∼ 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 (Ω m = Ω r ) redshift z eq ∼ 3369 (see Table I). Indeed, using the results given in Table I, we see that the Universe could be dominated by the expansion anisotropy as Ωr,eq+Ωm,eq = Ω σ 2 ,eq 2Ωm,eq 10 2 at the matter-radiation equality, and thereby conclude that the upper bound on Ω σ 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 eq , of a mode that enters the horizon, H eq a eq , at the matter-radiation transition [138,151].
In our model, k eq = H eq a eq = Heq (1+zeq) 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 Ω σ 2 ,0 . To do so, using (34), we write where we use a 0 = a(z = 0) = 1, and Ω M,0 = 1 − Ω m,0 − Ω r,0 −Ω σ 2 ,0 . We note that, for ω 50 and z eq ∼ 3369 (see Table I), the contribution to k eq from the term with Ω M,0 is negligible, k eq is by far more sensitive to Ω σ 2 ,0 and, for a given set of values of the parameters, k eq increases with increasing Ω σ 2 ,0 and decreases with increasing ω (viz., as the BD gravity approaches to GR). We obtain k eq = 0.01034 ± 0.00012 for Ω σ 2 ,0 = 0 (isotropic expansion), which is slightly larger than, but yet consistent with the isotropic ΛCDM model value k eq = 0.01024 ± 0.00007. Accordingly, switching to the massive BD (1) leads only to a slight increase in k eq , as expected from the constraint ω 50 (see Table I). It may be useful to note that these two values are consistent with the recent Planck release [12] value k eq = 0.010339 ± 0.000063 (TT,TE,EE+lowE+lensing+BAO) obtained for the base ΛCDM model. On the other hand, we see that the expansion anisotropy, although negligible today, shifts k eq to unrealistically large values, viz., we obtain k eq = 0.15159 ± 0.00327 for the upper bound Ω σ 2 ,0 = 10 −8.48 (see Table I). We then work out the values of Ω σ 2 ,0 that can shift k eq to reasonable values. See Fig. 5 showing explicitly k eq with respect to Ω σ 2 ,0 with errors to make a simple comparison between the models. We obtain k eq = 0.02824 ± 0.00057 by setting Ω σ 2 ,0 = 10 −10 , and k eq = 0.01067 ± 0.00013 by setting Ω σ 2 ,0 = 10 −12 , which are still inconsistent with the k eq values given for the isotropic ΛCDM model. But then, we obtain k eq = 0.01037 ± 0.00012 by setting Ω σ 2 ,0 = 10 −13 , and k eq = 0.01034 ± 0.00012 by setting Ω σ 2 ,0 = 10 −14 , where we notice that only the last decimals are different. We observe further that k eq does not change for Ω σ 2 ,0 10 −14 anymore and remains consistent with the isotropic ΛCDM model values obtained in this study as well as recent Planck release, which in turn implies that we cannot distinguish Ω σ 2 ,0 10 −14 from Ω σ 2 ,0 = 0 by means of k eq . We finally calculate using Ω σ 2 ,0 ∼ 10 −14 that the Universe is indeed dominated by matter+radiation at matter-radiation equality, viz., we now have Ω σ 2 ,eq 2Ωm,eq ∼ 10 −3 . Thus, by means of matterradiation transition, we conclude Ω σ 2 ,0 10 −14 along with ω 50, 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 I.

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 Λ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 ∼ 1 MeV to T ∼ 0.1 MeV and the age of the Universe from t ∼ 1 s to t ∼ 3 min corresponding to redshift z ∼ 3×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 ∼ 3 × 10 8 (implying z eq,σ 2 ,r 3 × 10 8 ) by using the constraints on the relevant parameters from Table I and found out that Ω σ 2 ,0 10 −21 for ω 50 in line with our preliminary investigation in Sect. 3.1 [cf., see Eq. (37)]. On the other hand, z eq,σ 2 ,r corresponds to the redshift at which ρ σ 2 /ρ 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 ρ σ 2 /ρ r ratio up to a few percents, viz., Hence, considering this condition, we can obtain a new constraint on Ω σ 2 ,0 stronger than the one given in (37) obtained from the condition z eq,σ 2 ,r 3 × 10 8 . To work this out, we first write, from (35), where we also used Ωm,0 Ωr,0 = ρm,0 ρr,0 = 1 + z eq . Then, we use in this equation the condition (74) along with the constraints on Ω m,0 and z eq given in Table I We note that, here, the contribution from the effective anisotropic pressure of the Jordan field during z < z eq , which is encoded in the term (1 + z eq ) 1 1+ω in (75), is not significant, such that it leads to only fifteen percent smaller upper bound value for Ω σ 2 ,0 , viz., (1+z eq ) − 1 1+ω 10 −0.07 = 0.85 for ω 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 ω 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 ∝ t 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: 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 ρ r = π 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 It is clear from (77) and (78) that we can set the following relationg * = G1 G0 g * implying that the altered G, i.e., G 1 , in the BD gravity can equivalently be treated as modified g * , namely,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 ≡ H HSBBN = g * g * , which can be adopted as S = G1 G0 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 η 10 (the number density of baryons n b to the number density of CMB photons n γ defined as η 10 = 10 10 n b /n γ ). 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 η 10 6.5 and 0.85 S 1.15 and, for 4 He (helium) and D (deuterium) these are 10 [155] Y p = 0.2381 ± 0.0006 + 0.0016[η 10 + 100(S − 1)], (79) y DP = 45.7(1 ± 0.06)[η 10 − 6(S − 1)] −1.6 , (80) where, Y p ≡ nHe n b and y DP ≡ 10 5 nD nH are 4 He and D mass fractions, respectively, and here in our study, η 10 = 273.9 Ω b,0 h 2 and S = G1 G0 (1 + z eq ) 1 2(1+ω) , see (33). We see that η 10 ≈ 6.1 for both models from the constraints on Ω b,0 h 2 given in Table I and find that the cosmological data we considered constrain the altered expansion rate of the Universe for z > z eq , hence during BBN, as 1 < S < 1.0665 at 68% C.L. and 1 < S < 1.0839 at 95% 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 p and y DP values, which in turn can be used for a further investigation of the constraint on the BD parameter ω.
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 p = 0.2449 ± 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-α systems studied in [158] giving y DP = 2.527 ± 0.030 (68 % C.L.). We note that the BBN-predicted Y p and y DP for both the Λ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 η 10 is fixed). We give a summary of our findings by depicting the 2D marginalized posterior distribution of the BBN-predicted Y P via (79)/y DP via (80) and ω in the top panel/bottom panel of Fig. 6. In which, for a comparison, we depict also the bands of the SBBN-predicted Y P via (79)/y DP via (80) for the Λ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 P and ω with decreasing ω as it approaches its lower bound, whereas for large values of ω the Y P − ω contour at 68% C.L. (dark green) approaches the band of the SBNN-predicted Y P at 68% C.L. (grey) as it should be. Besides these, more importantly, we see that the Y P − ω contour at 68% C.L. (dark green) stays above the band of the independent observational estimate band (red) for ω 250, namely, its consistency with the independent observational estimate of Y p = 0.2449 ± 0.0040 from [157] requires ω 250 at 68% C.L. providing us an improved constraint on the BD parameter ω. We see that the BBN-predicted y DP via (80) is insensitive to ω, y DP − ω contour at 68% C.L. (dark green) already covers the independent observational estimate band (red) for ω 250 (the improved constraint from our Y P investigation) but is much wider than it (due to the relatively large internal error in (80), viz., y DP ∝ 1 ± 0.06), and y DP − ω contour for BD gravity (dark green) approaches the SBBN-predicted y DP band at 68% C.L. (grey) for large ω values as it should be. These show that we are not able to deduce a new constraint on the BD parameter ω by comparing the BBN-predicted value of y DP via (80) with its independent observational estimate.
Thus, by means of the BBN, we reach the most stringent constraints on Ω σ 2 ,0 and ω, the two free parameters that determine the anisotropic BD extension to the standard ΛCDM model, as follows; Ω σ 2 ,0 10 −23 and ω 250.
Finally, considering these improved constraints, we find that the contribution to this constraint of Ω σ 2 ,0 from the anisotropy of the effective pressure of the Jordan field during z < z eq , which is encoded in the term (1 + z eq ) 1 1+ω in (75), is quite insignificant, such that, it leads to only a few percent smaller upper bound value for Ω σ 2 ,0 , viz., (1 + z eq ) − 1 1+ω 10 −0.01 = 0.977 for ω 250.

CONCLUSIONS
We have carried out an explicit detailed theoretical and observational investigation of an anisotropic mas-sive Brans-Dicke (BD) gravity extension of the standard ΛCDM model, wherein the extension is characterized by the two additional degrees of freedom: the so called BD parameter ω and the present day value of the density parameter corresponding to the shear scalar (a measure of the expansion anisotropy), Ω σ 2 ,0 . The role of the cosmological constant Λ is taken over by the Jordan field potential of the form U (ϕ) = 1 2 M 2 ϕ 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 11 . 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 ω, being the measure of the deviation from GR (|ω| → ∞), 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 ω, 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 σ 2 ∝ (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 ΛCDM model, so we have mainly confined our investigations to small deviations from this model via sufficiently small Ω σ 2 ,0 and large ω values. We have shown through some preliminary cosmological discussions that |ω| 10 (roughly) cannot be called as a small deviation, yet, for completeness we have extended our investigations to non-negative values of ω. Indeed, considering the combined data sets PLK+BAO+SN+H, we have obtained ω 50, M < 1.51 × 10 −34 eV and Ω σ 2 ,0 10 −8. 48 (Ω σ 2 ,0 10 −14 , when matter-radiation equality is considered). 12 Then by means of BBN, we have improved these constraints to ω 250 (in particular, from the comparison of the helium abundance prediction of the model with the direct measurements) and Ω σ 2 ,0 10 −23 . We have also found that the contribution of the anisotropy of the effective DE (viz., 11 In more complicated anisotropic spacetimes, the anisotropic spatial curvature contributes to the shear propagation equation and hence to the redshift dependency of the shear scalar. For instance, the most general spatially flat anisotropic spacetimes are of Bianchi type VII 0 , which, in addition to the simple expansionrate anisotropies, yield anisotropic spatial curvature that mimics traceless anisotropic fluid [48]. 12 We make use of a modified version of the simple and fast MCMC code that computes expansion rates and distances from the Friedmann equation, named SimpleMC [21,145]. The method we use constrains the expansion anisotropy through its contribution to the average expansion rate of the Universe. ∆w DE,0 < 4.23 × 10 −7 ) to this constraint on Ω σ 2 ,0 is insignificant, namely, led to only a few percent smaller (stronger) upper bound value. 13 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 ΛCDM model presents no significant deviations from it all the way to the BBN. The strongest cosmological constraints on ω present in the literature, e.g., ω 890 [116], would obviously just strengthen this conclusion. Moreover, we should further consider the strongest local constraint ω > 40000 [114] as well since the constraint M < 1.51×10 −34 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 ω, the simplest anisotropic massive BD gravity extension of the standard ΛCDM model cannot be distinguished from its simplest GR based anisotropic extension [95] when we consider the data from cosmological observations. 14 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 ω to be some functions of ϕ, the observations from the Cassini spacecraft place upon strong constraint on ω as in the original BD theory. This constraint, however, now only applies to the local value of ω(ϕ), i.e., of the present day value of ϕ in the Solar System. Similarly, it is possible that ω has spatial variation, so that, for instance, it could be constrained to be ω > 40000 locally, but could be in line with ω ∼ 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 13 It can be seen from H(z) given in (34) that the sufficiently large negative ω values also lead to small deviations from ΛCDM, but in this case we have slightly flatter redshift dependence of the expansion anisotropy, which in turn can relax the constraint on Ω σ 2 ,0 at few percents level compared to the one we have obtained along with large positive ω values. Nevertheless, for such large negative values of ω, being less than the scale invariant limit − 3 2 , one also should cope with stability issues. 14 For example, the deviation of the effective EoS parameter corresponding to the expansion anisotropy from w σ 2 = 1 (corresponding to the GR limit) becomes just ∼ 10 −5 for ω = 40000, whilst it is already very small, viz., ∼ 10 −3 , for the cosmological lower bound ω ∼ 250 (and ω ∼ 890).
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 Λ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 ω 10 2 (ω 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 Λ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 |ω| ∼ O(1), which, interestingly, is in principle the natural order of magnitude for the BD parameter (as, e.g., ω = −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 Λ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 ω, say, ω ∼ −1 but larger than the scale invariant limit − 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, ω = −1 corresponds to the low energy effective string action and ω ∼ −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 ω values, we carried out a discussion in Sect.3.3.2 showing how dramatically the redshift dependence of the expansion anisotropy can alter when ω ∼ −1. For instance, for ω = − 4 3 , expansion anisotropy becomes non-dynamical (mimicking Λ) and for ω = −1, expansion anisotropy should be set to zero. The region − 3 2 < ω < − 4 3 is also interesting that, in the presence of only dust (without Λ 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 ω consistent with local experiments and cosmo-logical 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 |ω| = 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 Λ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 ω-along with, for instance, the full Planck data -as in [93,94] providing constraints on the anisotropic expansion on the top of the Λ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 Λ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.