*Article* **Big Rip Scenario in Brans-Dicke Theory**

**Sasmita Kumari Pradhan 1,2, Sunil Kumar Tripathy 3,\*, Zashmir Naik 1, Dipanjali Behera 3 and Mrutunjaya Bhuyan 4,\***


**Abstract:** In this work, we present a Big Rip scenario within the framework of the generalized Brans-Dicke (GBD) theory. In the GBD theory, we consider an evolving BD parameter along with a self-interacting potential. An anisotropic background is considered to have a more general view of the cosmic expansion. The GBD theory with a cosmological constant is presented as an effective cosmic fluid within general relativity which favours a phantom field dominated phase. The model parameters are constrained so that the model provides reasonable estimates of the Hubble parameter and other recent observational aspects at the present epoch. The dynamical aspects of the BD parameter and the BD scalar field have been analysed. It is found that the present model witnesses a finite time doomsday at a time of *tBR* 16.14 Gyr, and for this scenario, the model requires a large negative value of the Brans-Dicke parameter.

**Keywords:** cosmological constant; generalised Brans-Dicke theory; Big Rip

**1. Introduction**

Late-time cosmic acceleration is one of the most bizarre and unsolved problems in modern cosmology. In scalar field cosmological models, the late-time cosmic acceleration issue is predominantly attributed to an exotic dark energy (DE) form that corresponds to a cosmic fluid having low energy density, as well as negative pressure. This is usually understood through a quantity dubbed as the equation of state (EoS) parameter *ωD* = *pρ* , where *p* represents the DE pressure and *ρ* symbolises the dark energy density. The dark energy with a negative pressure corresponds to a negative EoS parameter. Despite several attempts made by astronomers and cosmologists, the experimental determination of *ωD* remains challenging. Its precise estimation at the present epoch along with the knowledge of its development over a long period may unravel the mystery of the dark energy whose nature and origin remains speculative so far. In the ΛCDM model, the cosmological constant Λ with *ωD* = 1 plays the role of dark energy. However, in canonical scalar field models, quintessence fields or phantom fields shoulder the burden for the late-time cosmic speed-up, while the EoS parameter for the quintessence field lies in the range − 2 3 ≤ *ωD* ≤ −13 [1–3], which for the phantom fields, becomes *ωD* < −1 [4]. However, the EoS parameter as constrained from recent observational data favours a phantom phase in the Universe with *ωD* < −1 [5], while constraints from the CMB data in the nine-year WMAP survey sugges<sup>t</sup> that *ωD* = −1.073+0.090 −0.089 [6], a combination of the CMB data with Supernova data, predicts *ωD* = −1.084 ± 0.063 [7]. Other constraints on the EoS parameter include *ωD* = −1.035+0.055 −0.059 from Supernova cosmology project [8], *ωD* = −1.03 ± 0.03 from recent Planck 2018 results [9] and from Pantheon data *ωD* = −1.006 ± 0.04 [10]. In phantom dark energy models, the energy conditions are usually violated, and the

**Citation:** Pradhan S.K.; Tripathy, S.K.; Naik, Z.; Behera, D.; Bhuyan, M. Big Rip Scenario in Brans-Dicke Theory. *Foundations* **2022**, *2*, 128–139. https:// doi.org/10.3390/foundations2010007

Academic Editor: Eugene Oks

Received: 16 December 2021 Accepted: 11 January 2022 Published: 17 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Universe may witness a blowing-up of the curvature of space-time at a finite time, leading to the dissolution of the whole material Universe into pieces. This picture of the finite future of the Universe, dubbed as the Big Rip singularity, concerns the recent cosmological research [11,12]. The finite-time future singularity leads to inconsistencies which led to different proposals in recent times, including the quantum effects to delay singularity, possible gravity modification, or the coupling of dark energy and dark matter [13,14].

In recent times, because of the concern regarding the ultimate fate of the Universe for phantom accreted dark energy with *ωD* < −1, a number of cosmological models have been used to investigate the Rip cosmologies based on the general relativity (GR) and modified gravity. Darbowski et al. studied the Big Rip singularity and the ultimate cosmic fate [15,16], and Granda and Loaiza showed the occurrence of Big Rip for kinetic and Gauss-Bonnet coupling [17]. The classical and quantum fate of the Big Rip cosmology has been studied by Vasilev et al. [18]. Within the framework of *f*(*T*) gravity, Hanafy and Saridakis presented a cosmological model where the Universe may last forever in a Pseudo Rip scenario [19]. Recently, Ray et al. studied the Big Rip and some Pseudo Rip cosmological models in the context of an extended gravity theory [20], and Pati et al. investigated the possible occurrence of Rip scenarios in an extended symmetric teleparallel theory [21]. Big Rip singularity is a unique singularity that possibly occurs in a phantom scenario violating the Null Energy Condition. The occurrence of Big Rip singularity, for which the energy density and the scale factor of the Universe diverge, dissolutes the bound system and ultimately leads to the tearing up of the Universe in finite time. Such a scenario has become a major concern for cosmologists. As such, the present study is aimed at investigating the possibility of the occurrence of a Big Rip scenario within the framework of the generalised Brans-Dicke (GBD) theory. The GBD theory incorporates a self-interacting potential, as well as a dynamically varying Brans-Dicke(BD) parameter. Previously, Montenegro and Carneiro have investigated some cosmological models leading to Big Rip kinds of solutions within the Brans-Dicke theory in the presence of decaying vacuum density. In that work, they considered a time-independent negative BD parameter *ω* = −1 [22]. They obtained cosmological solutions with a time-varying deceleration parameter which may lead to negative energy density. In general, the BD theory is a most popular modified gravity theory proposed as an alternative to GR, where the gravity is mediated by a scalar field. The BD theory has already been used emphatically in addressing many issues in cosmology and astrophysics including the explanation for the inflationary scenario [23]. Over a period of time, the BD theory has been well-studied through different tests, including the gravitational radiation from gravitational wave bursts [24–27].

The manuscript is presented as follows: In Section 2, the basic field equations for the GBD theory are estimated for an anisotropic LRS Bianchi I (LRSBI) metric. Additionally, the dynamics of the GBD theory when incorporating a cosmological constant is appraised. Section 3 is devoted to a Big Rip scenario through a scale factor that diverges at a specific time, and discusses the time evolution of the BD scalar field, BD parameter, and the selfinteracting potential under the Big Rip scenario. Lastly, the conclusion is drawn, and a concise summary of the present investigation is given in Section 4. Throughout this work, we chose the natural unit system: 8*πG*0 = *c* = 1, where *G*0 symbolises the Newtonian gravitational constant at the present epoch and *c* represents the speed of light in vacuum.

### **2. Basic Formalism**

In a Jordan frame within the GBD theory with a self-interacting potential *V*(*φ*) and a time-dependent BD parameter *<sup>ω</sup>*(*φ*), we have the action as [28,29]

$$S = \int d^4x \sqrt{-g} \left[ \phi \mathcal{R} - \frac{\omega(\phi)}{\phi} \phi^{,\mu} \phi\_{,\mu} - V(\phi) + \mathcal{L}\_m \right],\tag{1}$$

where *R* defines the curvature and L*m* is the matter Lagrangian. The time-dependence aspect of the BD parameter emerges naturally in the Kaluza-Klein theory, string theory, or in the supergravity theory [30,31]. Different issues in cosmology have been investigated in recent times in the GBD framework with a time-dependent BD parameter [32–36]. Because of the GBD field equations may be obtained as [36,37]

$$G\_{\mu\nu} = \left. \frac{T\_{\mu\nu}}{\phi} + \frac{\omega(\phi)}{\phi^2} \left[ \Phi\_{,\mu}\Phi\_{,\nu} - \frac{1}{2} g\_{\mu\nu}\Phi\_{,\mu}\Phi^{,\mu} \right] + \frac{1}{\phi} \left[ \Phi\_{,\mu;\nu} - g\_{\mu\nu}\Box\phi \right] - \frac{V(\phi)}{2\phi} g\_{\mu\nu} \tag{2}$$

$$
\Box \phi\_- = \frac{T}{2\omega(\phi) + 3} - \frac{2V(\phi) - \phi \frac{\partial V(\phi)}{\partial \phi}}{2\omega(\phi) + 3} - \frac{\frac{\partial \omega(\phi)}{\partial \phi} \phi\_{,\mu} \phi^{,\mu}}{2\omega(\phi) + 3},\tag{3}
$$

where - represents the d'Alembert operator. A perfect fluid distribution with the energymomentum tensor *<sup>T</sup>μν* = (*ρ* + *<sup>p</sup>*)*<sup>u</sup>μuν* + *pgμν* is considered, so that the *T* = *<sup>g</sup>μνTμν* is the trace. In order to model the universe, we consider an LRSBI universe [38]

$$ds^2 = -dt^2 + A^2 dx^2 + B^2(dy^2 + dz^2). \tag{4}$$

Even though the present observable Universe is mostly isotropic and homogeneous and can be mostly described as an FRW metric, some of the observations obviously hint of a possible departure from isotropy [39–42]. The amount of cosmic anisotropy present may be very small, but we cannot simply rule out its possibility. Additionally, the LRSBI model resembles the flat FRW model, but allowed us to incorporate a small but finite anistropy in the model. The GBD field equations for the anisotropic model become [32,36]

$$(2k+1)\tilde{\zeta}^2 H^2 \quad = \underbrace{\rho}\_{\Phi} + \frac{\omega(\phi)}{2} \left(\frac{\dot{\phi}}{\Phi}\right)^2 - 3H \left(\frac{\dot{\phi}}{\Phi}\right) + \frac{V(\phi)}{2\Phi},\tag{5}$$

$$2\tilde{\xi}\dot{H} + 3\tilde{\xi}^2 H^2 = \begin{array}{c} -\frac{p}{\rho} - \frac{\omega(\phi)}{2} \left(\frac{\dot{\phi}}{\phi}\right)^2 - 2\tilde{\xi}H\left(\frac{\dot{\phi}}{\phi}\right) - \frac{\ddot{\phi}}{\phi} + \frac{V(\phi)}{2\phi}, \tag{6}$$

$$(k+1)\sharp \dot{H} + (k^2 + k + 1)\sharp^2 H^2 = \ -\frac{p}{\phi} - \frac{\omega(\phi)}{2} \left(\frac{\phi}{\phi}\right)^2 - (k+1)\sharp H \left(\frac{\phi}{\phi}\right) - \frac{\phi}{\phi} + \frac{V(\phi)}{2\phi}.\tag{7}$$

In the above equations, *ξ* = 3 *k*+2 is an anisotropic parameter where *k* is a positive constant that decides the relationship among the directional expansion rates: *A* ˙ *A* = *k B* ˙ *B* . One should note that the isotropic behaviour of the model may be obtained for *ξ* = 1. *H* = 13 *A*˙*A* + 2 *B*˙*B* = 1*ξ B*˙*B* is the Hubble parameter.

The BD scalar field satisfies the Klein-Gordon equation

$$\frac{\ddot{\Phi}}{\Phi} + 3H\frac{\dot{\Phi}}{\Phi} = \frac{\rho - 3p}{2\omega(\phi) + 3} - \frac{\frac{\partial\omega(\phi)}{\partial\phi}\dot{\Phi}^2}{2\omega(\phi) + 3} - \frac{2V(\phi) - \phi\frac{\partial V(\phi)}{\partial\phi}}{2\omega(\phi) + 3}.\tag{8}$$

The GBD theory may be recast as an effective GR picture by incorporating a cosmological constant Λ. In such a case, the total energy density and the total pressure respectively become *ρT* = *ρ* + Λ and *pT* = *p* − Λ, and the GBD field equations reduce to

$$(2k+1)\mathfrak{J}^2 H^2 \quad = \begin{array}{c} \rho \text{T} \ + \rho\_{\Phi'} \end{array} \tag{9}$$

$$\left(\frac{k+3}{2}\right)\xi H + \left(\frac{k^2+k+4}{2}\right)\xi^2 H^2 \quad = \quad - (p\_T + p\_\phi)\_\prime \tag{10}$$

where,

$$
\rho\_{\phi} \quad = \quad (2k+1)\ddot{\xi}^2 H^2 \triangle \phi - \Im H \phi + \frac{\omega(\phi)}{2} \left(\frac{\phi^2}{\phi}\right),
\tag{11}
$$

$$p\_{\phi} = -\left(\frac{k^2 + k + 4}{2}\right) \mathfrak{J}^2 H^2 \triangle \phi - \left(\frac{k + 3}{2}\right) \mathfrak{J}H \triangle \phi + \phi + \left(\frac{k + 3}{2}\right) \mathfrak{J}H\phi + \frac{\omega(\phi)}{2} \left(\frac{\phi^2}{\phi}\right),\tag{12}$$

where *φ* = 1 − *φ*, which ensures the reduction to GR behaviour for *φ* = 0. It is interesting to note here that the GBD theory provides an extra cosmic fluid which may shoulder the burden of late-time acceleration. The corresponding effective EoS parameter becomes

$$\omega\_{\it eff} = \frac{p\_{\theta} - \Lambda}{\rho\_{\theta} + \Lambda} = -1 + \frac{-\left(\frac{k+3}{2}\right)\xi H \triangle \phi + \left[ (2k+1) - \left(\frac{k^2 + k + 4}{2}\right) \right] \xi^2 H^2 \triangle \phi + f(\phi, \dot{\phi}, \ddot{\phi})}{\Lambda + (2k+1)\xi^2 H^2 \triangle \phi + g(\phi, \dot{\phi})},\tag{13}$$

with *f*(*φ*, *φ*˙, *φ*¨) = *φ*¨ + *<sup>k</sup>*+<sup>3</sup> 2 *ξ* − 3*Hφ*˙ + *ω*(*φ*)*<sup>φ</sup>*˙ 2*φ* and *g*(*φ*, *φ*˙) = −3*Hφ*˙ + *ω*(*φ*) 2 *φ*˙ 2*φ* . Inthelowredshiftepoch,therecanbesmallvaluesof|*φ*|,andconsequently,

 Λ becomes the dominant term in the denominator of the second term in Equation (13). Additionally, in the numerator of Equation (13), the contribution of *f*(*φ*, *φ*˙, *φ*¨) may be neglected compared to other terms. In the limit of |*φ*| → 0, the effective EoS may be expressed as

$$
\omega\_{\varepsilon f\uparrow} \simeq -1 - \frac{\triangle \phi}{\Lambda} \left[ -\left(\frac{k+3}{2}\right) \zeta \dot{H} + \left[ (2k+1) - \left(\frac{k^2+k+4}{2}\right) \right] \zeta^2 H^2 \right], \tag{14}
$$

which may be reduced to

$$
\omega\_{eff} \simeq -1 + 2\frac{\triangle \phi}{\Lambda} \xi \dot{H}.\tag{15}
$$

for a small departure from cosmic anisotropy. The above equation tells us that we ge<sup>t</sup> a quintessence-like phase for *φ* > 0 and a phantom-dominated phase for *φ* < 0, at least in the low redshift epochs. In an earlier work, it has been shown from the reconstruction of the BD scalar field from observational *<sup>H</sup>*(*z*) data that *φ* < 0. This shows phantom-like behaviour in the GBD theory [38].

### **3. A Big Rip Scenario**

An explanation of the late-time cosmic speed-up issue with dark energy requires that the EoS parameter should be *ωD* < −1/3. The cosmological constant corresponds to *ωD* = −1. For dark energy, cosmological models favouring *ωD* < −1 are usually dominated with phantom energy for which the energy density goes up with time violating the dominant energy condition. The energy density in a phantom dominated dark energy model is proportional to the scale factor. As a consequence, the scale factor blows up at a finite time *tBR* − *t*0 2 <sup>3</sup>|<sup>1</sup>+*ωD*| 1 *<sup>H</sup>*0√<sup>1</sup>−Ω*m* , where *H*0 represents the present value of the Hubble parameter and Ω 0.3 is the matter density parameter [43]. Such a scenario is termed as the Big Rip, whose occurrence dissolutes the bounded system [43–45]. The Big Rip scenario in Phantom models leads to a unique singularity in the Universe and can be associated with the fundamental quantum gravity formalism [46].

We consider a Big Rip scenario with the scale factor evolving as

$$a(t) \cong \left(t\_{BR} - t\right)^{\alpha},\tag{16}$$

where *tBR* is the epoch where the scale factor blows up. *α* is a constant parameter related to the EOS parameter as

$$\alpha = \frac{2}{3(1+\omega\_D)}.\tag{17}$$

In a phantom field-dominated Universe, the EoS parameter is usually less than unity, that is, *ωD* < −1, which requires that the constant parameter *α* appearing in the scale factor should be negative, that is, *α* < 0. For the given Big Rip scenario, the Hubble parameter and the deceleration parameter (DP) are, respectively, *H*(*t*) = − *α tBR*−*t* and *q* = −1 − *H* ˙ *H*<sup>2</sup> = −1 + 1*α* . Additionally, we have *H*˙ = − *H*<sup>2</sup>*α* and *H*¨ = 2*α*<sup>2</sup> *H*3. One should note that, while the Hubble parameter contains two adjustable parameters, the deceleration parameter contains only one parameter, *α*. The deceleration parameter in the present Big Rip scenario comes out to be a constant quantity which should be negative to provide

an accelerating model. Because of this, we may constrain the parameter *α* from some recent observational constraint on the deceleration parameter. In a recent work, Camarena and Marra used the observational data from supernovae in a redshift range to constrain the DP as *q*0 = −1.08 ± 0.29 [47]. The central value of the deceleration parameter immediately fixes up the scale factor parameter as *α* = −12.49. Consequently, the EoS parameter may be constrained as *ωD* = −1.0533. This value is in close agreemen<sup>t</sup> with some recent measurements, as mentioned earlier. Once *α* is fixed, the other parameter *tBR* may be obtained from the present value of the Hubble parameter. In Figure 1, we show the evolution of the Hubble parameter for the constrained value of *α* that predicts a finite-time future singularity. The Hubble parameter increases with the cosmic expansion and blows up at a finite future. Assuming *H*0 = 74.3 km s<sup>−</sup><sup>1</sup> Mpc−1, the present model predicts a Big Rip occurring at a cosmic time *tBR* 16.14 Gyr.

**Figure 1.** Hubble parameter as a function redshift. The Hubble parameter is in km s<sup>−</sup><sup>1</sup> Mpc−<sup>1</sup> units.

### *3.1. Brans-Dicke Scalar Field*

Considering the Big Rip scenario, the evolutionary aspect of the Brans-Dicke scalar field may be studied within the framework of the GBD theory. Algebraic simplification of the field Equations (6) and (7) leads to the evolution equation for the BD scalar field as

$$-\frac{\dot{H}}{H} - 3H = \frac{\dot{\Phi}}{\Phi}.\tag{18}$$

In terms of DP, we may express the evolution equation as

$$(q-2)H = \frac{\dot{\phi}}{\phi}.\tag{19}$$

In our model, we obtained the deceleration parameter to be a constant quantity, and consequently, the BD scalar field may be obtained by integrating Equation (19) as

$$
\phi = \frac{\phi\_0}{a\_0^n} a^n,\tag{20}
$$

where *n* = *q* − 2 = 1 *α* − 3. *a*0 and *φ*0 are respectively the value of the scale factor and the BD scalar field at the present epoch. This relation clearly articulates a power-law behaviour of the BD scalar field with respect to the scale factor. It is worth mentioning here that the use of power-law functional behaviour of the scalar field is quite common in the literature. Since the deceleration parameter is a negative quantity in our model, the BD scalar field should decrease with the cosmic expansion.

The scale factor may be expressed as *a* = − *αH α*, and consequently, the BD scalar field becomes 

$$\frac{\Phi}{\Phi\_0} = \left(\frac{H\_0}{H}\right)^{\text{na}},\tag{21}$$

which ultimately leads to

$$\frac{\Phi}{\Phi} = -3nH^2. \tag{22}$$

In terms of the redshift defined as *z* = *a*0*a* − 1, the BD scalar field reduces to

$$\frac{\Phi}{\Phi\_0} = (1+z)^{-n}.\tag{23}$$

In Figure 2, the evolutionary aspect of the BD field is shown. The BD scalar field shows a decreasing trend from large positive values to vanishingly small values at late cosmic times. One should note from the figure that the BD scalar field behaves like ∼(<sup>1</sup> + *z*)3 in low redshift epochs.

**Figure 2.** Evolutionary behaviour of the BD scalar field.

### *3.2. Brans-Dicke Parameter*

The BD parameter is a dynamical quantity in the GBD theory, and its behaviour depends on the dynamical behaviour of the BD scalar field. The BD parameter may be obtained from the GBD field Equations (5)–(7) as

$$
\omega(\phi) = \frac{1}{\left(\frac{\dot{\phi}}{\dot{\phi}}\right)^2} \left[ -\frac{\rho + p}{\phi} - \frac{\ddot{\phi}}{\dot{\phi}} + k\xi^\jmath H \frac{\dot{\phi}}{\dot{\phi}} - 2\xi \dot{H} + 2(k-1)\xi^2 H^2 \right]. \tag{24}
$$

It should be mentioned here that, for a given Big Rip scenario as specified by a Hubble parameter *H* and a given cosmic anisotropy *ξ*, the dynamics of the BD scalar field is suitably obtained. Once the BD scalar field is obtained, the time-dependent aspect of the BD parameter now requires an equation of state *p* = *<sup>ω</sup>Dρ*, a relationship between the pressure *p* and the energy density *ρ*. Replacing *ρ* + *p* by (1 + *<sup>ω</sup>D*)*ρ* and using the fact that *φ*˙ *φ* = *nH* in Equation (24), we ge<sup>t</sup>

$$
\omega(\phi) = n^{-2} H^{-2} \left[ -\left( \frac{1 + \omega\_D}{\phi} \right) \rho + \left( (k\_\star^\pi + 3)n + 2(k - 1)\xi^2 \right) H^2 \right]. \tag{25}
$$

For the given Big Rip scenario, the conservation equation for the cosmic fluid

$$
\rho + 3H(\rho + p) = 0\tag{26}
$$

can be reduced to

$$\frac{\dot{\rho}}{\rho} = -\frac{2}{\alpha}H.\tag{27}$$

The energy density may be obtained from the integration of the conservation equation as

$$
\rho = \rho\_0 \left(\frac{H}{H\_0}\right)^2,\tag{28}
$$

where *ρ*0 is the present value of the energy density. Here, we used the fact that *a a*0 = *H*0 *H α* . It is now straight-forward to obtain the pressure as

$$p = \left(\frac{2}{3a} - 1\right) \rho\_0 \left(\frac{H}{H\_0}\right)^2,\tag{29}$$

so that *p* + *ρ* = 2 3*α ρ*0 *H H*0 2. The Brans-Dicke parameter may now be expressed in terms of the Hubble parameter as

$$
\omega(H,\xi,\mathfrak{a}) = \omega\_0(H\_0,\xi,\mathfrak{a}) - \chi \left[ \left(\frac{H}{H\_0}\right)^{(na-2)} - 1 \right],\tag{30}
$$

where *χ* = <sup>2</sup>*ρ*0*H*((<sup>3</sup>+*n*)*<sup>α</sup>*−<sup>1</sup>) 0 <sup>3</sup>*αφ*0*n*<sup>2</sup> and *<sup>ω</sup>*0(*<sup>H</sup>*0, *ξ*, *α*) = *kξ*+3 *n* + 2(*k* − 1) *ξ n* 2 − *χ* is the present value of the BD parameter.

One may note that the anisotropy affects the Brans-Dicke parameter, but it does not participate in its evolution. Only the non-evolving part of the Brans-Dicke parameter is modified in the presence of cosmic anisotropy. In fact, as is evident from Equation (30), the cosmic anisotropy brings about a change in the required present value of the BD parameter. Similar results have already been observed in an earlier work [32], where it was shown that for a power-law expansion and an exponential expansion law of the scale factor, the cosmic anisotropy affects only that part of the BD parameter which does not evolve with time. The similarity between this Big Rip model and that of the power-law behaviour and the exponential expansion model is that the deceleration parameter is non-evolving. Because of this, we may infer that, for the time-independent deceleration parameter, the cosmic anisotropy will not contribute to the evolution of the Brans-Dicke parameter. However, for models with a time-dependent deceleration parameter, the cosmic anisotropy affects the BD parameter as a whole [36,38]. In Figure 3, the BD parameter (normalized to its value at the present epoch) is shown for a representative value of the cosmic anisotropy *ξ*. It is observed that the BD parameter increases with the cosmic expansion. Since the Hubble parameter blows up at a time *t tBR*, the BD parameter also blows up at that epoch. In Figure 4, the evolutionary aspect of the BD parameter with respect to the BD scalar field is shown. With an increase in the BD scalar field, *<sup>ω</sup>*(*<sup>H</sup>*, *ξ*, *α*) decreases from a much higher value to low positive values. However, the decrement in *<sup>ω</sup>*(*<sup>H</sup>*, *ξ*, *α*) slows down for higher values of the scalar field. One may estimate the BD parameter at the present epoch within the purview of the present formalism to ge<sup>t</sup> a rip scenario at a finite future. In fact, for

a finite time future singularity such as the one discussed in the present work within GBD theory, we require the BD parameter to be approximately *ω*(*z* = 0) −3.6 × 10<sup>9</sup> at the present epoch.

**Figure 3.** Brans-Dicke parameter as a function of redshift.

**Figure 4.** Evolutionary behaviour of the Brans-Dicke parameter with respect to the BD field.

### *3.3. Self-Interacting Potential*

Within the GBD formalism, the time-dependence aspect of the self-interacting potential may be obtained as

$$V(\phi) = 2\phi \left[ \left\{ (2k+1)\ddot{\xi}^2 + 3n - \frac{\omega(\phi)}{2}n^2 \right\} H^2 - \frac{\rho}{\phi} \right]. \tag{31}$$

In order to assess the effect of the cosmic anisotropy, we may express the self-interacting potential as

$$\frac{V(\phi)}{2\phi} = \left\{ (2k+1)\zeta^2 - \frac{\omega\_0}{2}n^2 \right\} H^2 + \left[ 3n + \frac{\chi n^2}{2} \left( \left( \frac{H}{H\_0} \right)^{(na-2)} - 1 \right) \right] H^2 - \frac{\rho}{\phi}. \tag{32}$$

In the above Equation (32), the first term in the right-hand side bears the contribution of the cosmic anisotropy. In Figure 5, we show the evolutionary behaviour of the selfinteracting potential for a given cosmic anisotropy. Since in the Big Rip scenario, we discussed in the present work a large negative value of *ω*0 = *ω*(*z* = 0) is being required, the effect of the cosmic anisotropy becomes negligibly small. However, in other scenarios such as a bouncing one, we may ge<sup>t</sup> a substantial effect of the cosmic anisotropy on the self-interacting potential [36]. Because of this, we chose a representative value *k* = 1.0001 or a corresponding *ξ* = 0.999967 to plot the figure. The value of the self-interacting potential in the figure is normalized to its present value. One should note that, with the growth of cosmic expansion, *V*(*φ*) decreases from a higher value to small values at late times. This behaviour may be translated in terms of the scalar field to infer that the self-interacting potential shows an increasing trend with the BD scalar field.

**Figure 5.** Evolutionary behaviour of the self-interacting potential.

### **4. Summary and Conclusions**

In the present work, we studied a Big Rip scenario within the framework of a generalized Brans-Dicke theory. An LRSBI Universe is considered to incorporate directional anisotropy in the expansion rates. Such a model provides a more general approach compared to the FRW model. The generalized Brans-Dicke theory having a cosmological constant can be recast as a GR-like theory with cosmic fluid dominated by dark energy. The effective cosmological constant for such effective cosmic fluid dominated by dark energy may provide a quintessence-like or phantom-like behaviour depending on the nature of the scalar field. On the basis of the nature of the BD scalar field which has been reconstructed from the observational *<sup>H</sup>*(*z*) data, we showed that the present model favours a phantom dark energy model. In Phantom models, the energy density and the scale factor may grow sharply within a finite time, leading to a Big Rip situation. We considered a Big Rip scale factor with an arbitrary parameter, that was fixed from a recent observational value of the deceleration parameter, which ultimately fixes the effective EoS parameter as *ωD* = −1.05

in close conformity with recent observational estimates. Additionally, the present model witnesses a Big Rip scenario at a time *tBR* 16.14 Gyr. We studied the dynamical aspects of the scalar field, Brans-Dicke parameter, and the self-interacting potential. In conformity with other scalar field models, the BD scalar field shows a decreasing behaviour with cosmic time. The Brans-Dicke parameter, on the other hand, increases with the cosmic expansion. To the BD scalar field, the BD parameter (as normalized to its present value) decreases from a high value to almost constant values for higher values of the scalar field. In some recent tests concerning the gravitational radiation from gravitational wave bursts on the BD theory, the bounds on the BD parameter may be *ω* ≥ 104∼10<sup>5</sup> [27] or *ω* ≥ 10<sup>6</sup> [24]. Montenegro et al. used *ω* = −1 to obtain Big Rip kinds of solutions in the BD theory [22]. However, our present model requires a high negative value of the BD parameter at the present epoch of the order of ∼−10<sup>9</sup> to witness a Big Rip scenario in the finite future.

In our model, the anisotropy parameter does not affect the BD scalar field. It affects the BD parameter, but does not contribute to its time evolution aspect. In principle, only the non-evolving part of the BD parameter is affected by the anisotropy parameter. This is a feature usually observed for BD gravity models with a constant deceleration parameter [32]. However, the anisotropy affects the self-interacting potential. Since we require quite a large magnitude of the present epoch value of the BD parameter to ge<sup>t</sup> a viable Big Rip scenario, a small variation of the cosmic anisotropy does not substantially bring about a change in its numerical value.

**Author Contributions:** Conceptualization, S.K.T. and Z.N.; Methodology, S.K.T. and S.K.P.; Software, S.K.P.; Formal analysis, S.K.P., S.K.T. and D.B.; Investigation S.K.P., Z.N. and M.B.; Writing—original draft preparation, S.K.P. and S.K.T.; Writing—review and editing, M.B. and Z.N.; Supervision, S.K.T. and M.B.; Project administration, S.K.T. and M.B.; Funding acquisition, M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by FOSTECT Project No. FOSTECT.2019B.04, FAPESP Project No. 2017/05660-0, and the CNPq—Brasil.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The Data are contained within the article.

**Acknowledgments:** SKT thanks IUCAA, Pune (India) for providing support through the visiting Associateship programme and MB thanks T. M. Joshua for the discussion throughout the work.

**Conflicts of Interest:** The authors declare no conflict of interest.
