1 Introduction

Over the last few years, different cosmological observations such as type Ia supernovae [1, 2], cosmic microwave background radiation (CMBR) [3, 4], Baryon acoustic oscillations (BAO) [5, 6], galaxy redshift survey [7] and large scale structure [8, 9] strongly suggest that our universe is currently undergoing a phase of accelerated expansion. Also, it is believed that the mysterious force known as dark energy (DE) with huge negative pressure is responsible for the current expanding universe with an acceleration. The present Planck data tells that there is 68.3% DE of the total energy contents of the universe. We know very well that the standard cosmology has been extraordinarily successful but it unable to solve some serious issues including the search for the best DE candidate. Still there is some uncertainty in the origin and composition of DE except some particular ranges of the equation of state (EoS) parameter \(\omega _{de}\). The approaches to answer this DE problem fall into two representative categories: one is to introduce dynamical DE in the right-hand side of the Einstein field equations in the framework of general relativity and the second one is to modify the left hand side of the Einstein equations, leading to a modified theories of gravity. In the literature, there are very nice reviews on both dynamical DE models and modified theories  [10,11,12,13,14,15,16,17]. In the absence of any strong argument in favor of DE candidate, a variety of DE models have been discussed. Cosmological constant is the primary DE candidate for describing DE phenomenon but it has some serious problems like fine tuning and cosmic coincidence. Due to this reason, several dynamical DE models include a family of scalar fields such as quintessence [18,19,20,21], phantom [22,23,24,25], quintom [26, 27], tachyon [28,29,30], K-essence [31], various Chaplygin gas models like generalized Chaplygin gas, extended Chaplygin gas and modified Chaplygin gas [32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47] have been developed.

The holographic DE (HDE) model [48] has been suggested in the context of quantum gravity with the help of holographic principle [49]. This holographic principle says that the bound on the vacuum energy \(\varLambda \) of a system with size L should not cross the limit of the black hole (BH) mass having the same size due to the formation of BH in quantum field theory [50]. The energy density of HDE is defined as

$$\begin{aligned} \rho _{de}=3d^2m_p^2L^{-2}. \end{aligned}$$
(1)

where \(m_p\) is the reduced Planck mass. This HDE model gives the relationship between the quantum fields energy density in vacuum to the various cut-offs such as infrared and ultraviolet. Granda and Oliveros [51] proposed a IR cut-off containing the local quantities of Hubble and time derivative Hubble scales. The advantages of HDE with Granda and Oliveros cutoff (new HDE model) is that it depends on local quantities and avoids the causality problem appearing with event horizon IR cutoff. The new HDE model can also obtain the accelerated expansion of the universe and also showed that the transition redshift from deceleration phase (\(q>0\)) to acceleration phase (\(q<0\)) is consistent with current observational data [51, 52]. Nowadays, HDE attracted attention as it can alleviate the issue of cosmic coincidence, i.e., why the energy densities due the dark matter and the DE should have a constant ratio for the present universe [53]. Also, various works show that the HDE model is in fairly good agreement with the observational data [54,55,56,57]. Nojiri and Odintsov [58] have proposed unifying approach to early-time and late-time universe based on generalized HDE and phantom cosmology, and recently generalized this idea as Hinflation [59]. In recent years, various entropy formalism have been used to investigate the gravitational and cosmological setups. Also, some new HDE models are constructed such as Tsallis HDE (THDE) [60, 61], Renyi HDE model (RHDE) [62] and Sharma-Mittal HDE (SMHDE) [63]. Among these models, RHDE based on the absence of interactions between cosmos sectors, and this model shows more stability by itself [62]. SMHDE is classically stable in the case of non-interacting cosmos. The THDE model based on the Tsallis generalized entropy, which is never stable at the classical level [60, 61, 63]. Hence, with this motivation, in this work we consider the HDE with another entropy formalism i.e., Tsallis HDE.

As mentioned above, another approach to explore the present accelerated expansion of the universe is the modified theories of gravity. Standard Einstein’s theory of gravitation may not completely explain the gravity at very high energy. In this situation cosmic acceleration would arise not from DE as a substance but rather from the dynamics of modified gravity [12]. The simplest alternative to general relativity is the scalar-tensor theory obtained by Brans and Dicke [64] (BD). But, in case of w parameter value, there is a major difference between the theoretical and observational data. It is observed that the theoretical values are much less than that of observational value which motivates many researchers to explore various aspects of universe in BD framework [65,66,67,68,69,70,71,72,73,74,75,76]. Nojiri et al. [77] have studied the properties of singularities in the phantom DE universe. They have, also, mentioned that the phantom-like behavior of EoS parameter \(\omega _{de}\) may appear from BD theory, either from the non-minimal coupling of a scalar Lagrangian with gravity, or from negative (non-standard) potentials, or even the usual matter may appear in phantom-like nature. Recently, Saridakis et al. [78] have discussed HDE through Tsallis entropy and its cosmological evolution through observational constraints. Barboza et al. [79] and Nunes et al. [80] have studied DE models through non-extensive Tsallis entropy framework and cosmological viability of non-gaussian statistics. Agostino [81] has investigated the holographic principle through the nonadditive Tsallis entropy, used to describe the thermodynamic properties of nonstandard statistical systems such as the gravitational ones. Zadeh et al. [82] have explored the effects of different infrared cutoffs, such as the particle, Ricci horizons and the Granda–Oliveros cutoffs, on the properties of THDE model. Sharma and Pradhan [83] have investigated diagnosing THDE models with statefinder and \(\omega _{de}-\omega _{de}^\prime \) plane analysis. Sadri [84] has studied observational constraints on interacting THDE model. Sharif and Saba [85] have reconstructed the THDE model with Hubble horizon within the framework of f(G, T) gravity. Nojiri et al. [86] have studied modified cosmology from extended entropy with varying exponent. Ghaffari et al. [87] have discussed interacting and non-interacting THDE models by considering the Hubble horizon as the IR cutoff within BD scalar theory framework, while Jawad et al. [88] have studied cosmological implications of THDE in modified version of BD scalar theory. In both the models the authors have considered the BD scalar field \(\phi \) as a power function of average scale factor a(t). Here, we are interested to extend the study of THDE models in BD theory with scalar field \(\phi \) as logarithmic function of average scale factor.

In this work, we are interested in studying the both non-interacting and interacting Tsallis holographic dark energy in Brans–Dicke scalar–tensor theory by considering homogeneous, isotropic FRW flat universe. Here, we focus our attention on the THDE models in BD theory with logarithmic expansion of average scale factor a(t) for BD scalar field \(\phi \). The present work is organized as follows: in Sect. 2, we derive BD field equations in the presence of THDE and pressure less dark matter (DM). Also, we have constructed non-interacting and interacting THDE models along with their complete physical discussion. Finally, in Sect. 3, we present the conclusions of this paper and also made a comparative analysis.

2 Tsallis holographic dark energy in BD theory

We consider the homogeneous, isotropic and flat FRW metric in the form

$$\begin{aligned} ds^2=dt^2-a^2\{{dr^2}+r^2(d\theta ^2+sin^2~\theta d\psi ^2)\}, \end{aligned}$$
(2)

where a(t) is the scale factor of the model. The spatial volume (V), Hubble parameter (H) and deceleration parameter (q) of this model are given by

$$\begin{aligned} V= & {} a^3 \end{aligned}$$
(3)
$$\begin{aligned} H= & {} \frac{{\dot{a}}}{a} \end{aligned}$$
(4)
$$\begin{aligned} q= & {} -\frac{{\dot{H}}}{H^2}-1. \end{aligned}$$
(5)

Different theories of gravitation have been proposed as alternatives to Einstein’s general theory of gravity. But, the scalar–tensor theory formulated by Brans and Dicke [64] is supposed to be the best alternative to Einstein’s theory. We consider the universe filled with pressure less DM with energy density \(\rho _m\) and DE with density \(\rho _{de}\). Hence, in this case the BD field equations for the combined scalar and tenor fields are given by

$$\begin{aligned}&R_{ij}-\frac{1}{2}Rg_{ij}=-\frac{8\pi }{\phi }(T_{ij}+{\overline{T}}_{ij})-\phi ^{-1}\left( \phi _{i;j}-g_{ij}\phi ^{,k}_{;k}\right) \nonumber \\&\quad -w\phi ^{-2}\left( \phi _{,i}\phi _{,j}-\frac{1}{2}g_{ij}\phi _{,k}\phi ^{,k}\right) , \end{aligned}$$
(6)
$$\begin{aligned}&\phi ^{,k}_{;k}=\frac{8\pi }{(3+2w)}(T+{\overline{T}}) \end{aligned}$$
(7)

and the energy conservation equation is

$$\begin{aligned} (T_{ij}+{\overline{T}}_{ij})_{;j}=0, \end{aligned}$$
(8)

which is a consequence of field equations (6) and (7). Here, R and \(R_{ij}\) are the Ricci scalar and Ricci tensor respectively, w is a dimensionless coupling constant. \(T_{ij}\) and \({\overline{T}}_{ij}\) are energy-momentum tensors for pressure less DM and THDE, respectively, which are defined as

$$\begin{aligned} T_{ij}= & {} \rho _{m} u_iu_j \end{aligned}$$
(9)
$$\begin{aligned} {\overline{T}}_{ij}= & {} (\rho _{de}+p_{de})u_iu_j-p_{de}g_{ij} \end{aligned}$$
(10)

where \(p_{de}\) and \(\rho _{de}\) are the pressure and energy density of DE respectively and \(\rho _m\) is energy density of DM.

The field equations (6) and (7) for the metric (2) are obtained as

$$\begin{aligned}&2\frac{\ddot{a}}{a}+\frac{{\dot{a}}^{2}}{a^{2}}+\frac{w}{2}\frac{{\dot{\phi }}^{2}}{\phi ^{2}}+2\frac{{\dot{a}}{\dot{\phi }}}{a\phi }+\frac{\ddot{\phi }}{\phi }=-\frac{\omega _{de}\rho _{de}}{\phi } \end{aligned}$$
(11)
$$\begin{aligned}&3\frac{{\dot{a}}^{2}}{a^{2}}-\frac{w}{2}\frac{{\dot{\phi }}^{2}}{\phi ^{2}}+3\frac{{\dot{a}}{\dot{\phi }}}{a\phi }=\frac{\rho _{de}+\rho _{m}}{\phi } \end{aligned}$$
(12)
$$\begin{aligned}&{\ddot{\phi }}+3{\dot{\phi }}\frac{{\dot{a}}}{a}=\frac{\rho _{de}(1-3\omega _{de})+\rho _{m}}{3+2w}\nonumber \\ \end{aligned}$$
(13)

and the energy conservation equation (8), leads to

$$\begin{aligned} {\dot{\rho }}_{de}+{\dot{\rho }}_{m}+3H(\rho _{de}(1+\omega _{de})+\rho _m)=0, \end{aligned}$$
(14)

where overhead dot denotes ordinary differentiation with respect to time t. Here, \(\omega _{de}\) is the equation of state (EoS) parameter of dark energy and is given by

$$\begin{aligned} \omega _{de}=\frac{p_{de}}{\rho _{de}}. \end{aligned}$$
(15)

In literature it is also common to use a power-law relation between BD scalar field \(\phi \) and average scale factor ’a’ of the form \(\phi =\phi _0a^l\) [89, 90], where \(\phi _0\) is a constant and l is a power index. Various authors have discussed different aspects of this form of scalar field \(\phi \) and have shown that it leads to constant deceleration parameter [72, 91] and also time varying deceleration parameter [69, 92]. Recently, Kumar and Singh [93] have introduced a BD scalar field evolves as a logarithmic function of the average scale factor to study the evolution of holographic and new agegraphic DE models. The relation is given by

$$\begin{aligned} \phi =\phi _1~ln(\beta _1+\beta _2~a(t)) \end{aligned}$$
(16)

where \(\phi _1\), \(\beta _1>1\) and \(\beta _2>0\) are constants. Recently, Singh and Kumar [94], Sadri and Vakili [95], and Aditya and Reddy [96] have investigated holographic DE models in BD theory using this logarithmic law for scalar field.

The energy density of Tsallis holographic DE model is given by [61]

$$\begin{aligned} \rho _{de}=\eta L^{2\delta -4} \end{aligned}$$
(17)

where \(\eta \) is a parameter with dimensions \([L]^{-2\delta }\). For \(\delta =1\) the above equation gives the usual holographic DE \(\rho _{de}=3d^2m_p^2\), with \(\eta =3d^2m_p^2\) and \(d^2\) the model parameter. Also, it is interesting to mentioning that in the special case \(\delta =2\) the above equation gives the standard cosmological constant model \(\rho _{de}=constant=\varLambda \) (Saridakis et al. [78]).

By considering the Hubble horizon as the IR cutoff, \(L =H^{-1}\), in BD theory the energy density (17) takes the form

$$\begin{aligned} \rho _{de}=3d^2\phi H^{-2\delta +4}. \end{aligned}$$
(18)

The dimensionless density parameters are defined as

$$\begin{aligned} \varOmega _{m}=\frac{\rho _{m}}{\rho _{cr}};~~\varOmega _{de}=\frac{\rho _{de}}{\rho _{cr}};~~\varOmega _{\phi }=\frac{\rho _{\phi }}{\rho _{cr}} \end{aligned}$$
(19)

where \(\rho _{cr}=3m_p^2H^2\) is called the critical energy density and in BD theory it can be written as \(\rho _{cr}=3\phi H^2\).

In the following sections, we consider the two cases: non-interacting model and interacting model. We determine, in both the cases energy density of THDE \(\rho _{de}\), EoS parameter \(\omega _{de}\), deceleration parameter q, squared sound speed \(v_s^2\) and \(\omega _{de}-\omega _{de}^\prime \) plane by solving the BD field equations. We also study their physical behavior.

2.1 Non-interacting model

First we consider that there is no energy exchange between the two fluids (cosmic sectors), and hence, the energy conservation equation (8) leads us to the following separate conservation equations:

$$\begin{aligned} {\dot{\rho }}_{de}+3H\rho _{de}(1+\omega _{de})= & {} 0, \end{aligned}$$
(20)
$$\begin{aligned} {\dot{\rho }}_{m}+3H\rho _m= & {} 0. \end{aligned}$$
(21)

Taking differentiation with respect to time t for Eq. (20), and using BD scalar field Eq. (16), we have

$$\begin{aligned} {\dot{\rho }}_{de}=\rho _{de}H\left( \frac{\beta _2a}{(\beta _1+\beta _2 a)~ln(\beta _1+\beta _2 a)}+(4-2\delta )\frac{{\dot{H}}}{H^2}\right) .\nonumber \\ \end{aligned}$$
(22)

By taking the time derivative of Eq. (12), using Eqs. (16) and (18)–(21), we obtain

$$\begin{aligned} \frac{{\dot{H}}}{H^2}= & {} -\biggl \{\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~ln(\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[ln(\beta _1+\beta _2 a)]^2}\nonumber \\&+9\varOmega _{de}(1+\omega _{de}+u)\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~ln(\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[ln(\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}\biggr \}\nonumber \\&\times \biggl \{6+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[ln(\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[ln(\beta _1+\beta _2 a)]^2}\biggr \}^{-1} \end{aligned}$$
(23)

where \(u=\frac{\rho _m}{\rho _{de}}=\frac{\varOmega _m}{\varOmega _{de}}\). On the other side, from Eqs. (18) and (20), we obtain

$$\begin{aligned} \frac{{\dot{H}}}{H^2}=\frac{3}{2(\delta -2)}\left[ 1+\omega _{de}+\frac{\phi _1\beta _2 a}{3(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\right] .\nonumber \\ \end{aligned}$$
(24)

From the above two Eqs. (23) and (24), we find that

$$\begin{aligned} \frac{{\dot{H}}}{H^2}= & {} -\biggl \{\frac{3\beta _2 a(2+\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}+9\varOmega _{de}u\biggr \}\nonumber \\&\times \biggl \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\biggr \}^{-1} \end{aligned}$$
(25)

Taking time derivative of Eq. (19), we obtain

$$\begin{aligned} {\dot{\varOmega }}_{de}=2\varOmega _{de} ~(1-\delta )~\frac{{\dot{H}}}{H}. \end{aligned}$$
(26)

In order to observe the behavior density parameter of THDE, we define \(\varOmega _{de}^\prime =\frac{{\dot{\varOmega }}_{de}}{H}\), where the prime denotes derivative with respect to ‘\(\ln ~a(t)\)’. Then from Eqs. (25) and (26) we have

$$\begin{aligned} \varOmega _{de}^{\prime }= & {} 2\biggl \{\frac{3\beta _2 a(2+\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}+9\varOmega _{de}u\biggr \}\nonumber \\&\times \varOmega _{de}(\delta -1)\nonumber \\&\times \biggl \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\biggr \}^{-1}. \end{aligned}$$
(27)

Also, from Eqs. (23) and (24) we obtain EoS parameter of THDE as

$$\begin{aligned} \omega _{de}= & {} -1-\Biggl \{\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}+9\varOmega _{de}u\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{\beta _2 a}{2(\delta -2)(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&\times \biggl \{6+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\biggr \}\Biggr \}\nonumber \\&\times \Biggl \{\frac{3}{2(\delta -2)}\biggl \{6+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}+9\varOmega _{de}\biggr \}\Biggr \}^{-1}. \end{aligned}$$
(28)
Fig. 1
figure 1

Plot of \(\varOmega _{de}\) of non-interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image

To study the behavior of dimensionless density parameter of THDE \(\varOmega _{de}\), we solve the Eq. (27) numerically corresponding to redshift z whose output is given in Fig. 1. In the evolution of dynamical parameters the THDE parameter \(\delta \) plays a crucial role and hence we have taken various values to \(\delta =0,~1,~1.1,~1.3,\) 1.5,  2. Here, we have used the initial value of \(\varOmega _{de}\) as \(\varOmega _{de}^0=0.73\) and additionally we have taken \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\) [87, 97] and \(u=0.3\) [98]. It can be seen that the \(\varOmega _{de}\) is positive and decreasing function throughout the evolution, and finally it approaches to different positive constant values for \(\delta =1,~1.1,~1.3,~1.5,~2\) whereas the density parameter is increasing function for \(\delta =0\). Also, we observe that the steepness of density parameter \(\varOmega _{de}\) increases with increase in \(\delta \). Observations from Planck’s data (2014) [99] have given the following constraints on the DE density parameter \(\varOmega _{de}=0.717^{+0.023}_{-0.020}\) (Planck \(+\) WP) and \(\varOmega _{de}=0.717^{+0.028}_{-0.024}\) (WMAP-9) and Planck’s data (2018) [100] have given the constrains on DE density parameter as \(\varOmega _{de}= 0.679 \pm 0.013\) (TT \(+\) lowE), \(0.699 \pm 0.012\) (TE \(+\) lowE), \(0.711^{+0.033}_{-0.026}\) (EE \(+\) lowE), \(0.6834 \pm 0.0084\) (TT, TE, EE \(+\) lowE), \(0.6847 \pm 0.0073\) (TT, TE, EE \(+\) lowE \(+\) lensing), \(0.6889\pm 0.0056\) (TT, TE, EE \(+\) lowE \(+\) lensing \(+\) BAO) by implying various combination of observational schemes at 68% confidence level. We observed that the THDE density parameter meets the above mentioned limits at present epoch, which shows that our results are consistent.

The EoS parameter is the relationship between pressure \(p_{de}\) and energy density \(\rho _{de}\) of DE whose expression is given by \(\omega _{de}=\frac{p_{de}}{\rho _{de}}\). The EoS parameter is used to classify the decelerated and accelerated expansion of the universe and it categorizes various epochs as follows: when \(\omega =1\), it represents stiff fluid, if \(\omega =1/3\), the model shows the radiation dominated phase while \(\omega =0\) represents matter dominated phase. In DE dominated accelerated phase, \(-1<\omega <-1/3\) shows the quintessence phase and \(\omega =-1\) shows the cosmological constant, i.e., \(\varLambda \)CDM model and \(\omega <-1\) yields the phantom era. In Fig. 2, we have plotted EoS parameter versus redshift z for different values of \(\delta \). We observe from Fig. 2 that the EoS parameter starts from matter dominated era, then goes towards quintessence DE era and finally approaches to vacuum DE era for different values of \(\delta =0,~1,~1.1,\) 1.3,  1.5. For \(\delta =2\), it is clear from the Eq. (28) that the EoS parameter becomes \(-1\), i.e., cosmological constant. This is in agreement with the results obtained by Saridakis et al. [78]. Also, it is interesting to mention here that as \(\delta \) decreases the transition from matter dominated phase to dark energy phase is delayed considerably. It may be noted that the EoS parameter, in this non-interacting case, never crosses the phantom divided line (PDL) \(\omega _{de}=-1\).

Fig. 2
figure 2

Plot of EoS parameter of non-interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image

Deceleration parameter is another important cosmological parameter by means of which one can distinguish current accelerated or early decelerated expansions. It is defined as

$$\begin{aligned} q=-1-\frac{{\dot{H}}}{H^2}. \end{aligned}$$
(29)

In this case, with the help of Eq. (25) this deceleration parameter can be obtained as

$$\begin{aligned} q= & {} -1+\biggl \{\frac{3\beta _2 a(2+\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}+9\varOmega _{de}u\biggr \}\nonumber \\&\times \biggl \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\biggr \}^{-1} \end{aligned}$$
(30)
Fig. 3
figure 3

Plot of deceleration parameter of non-interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image

The plot of above constructed deceleration parameter versus redshift z is shown in Fig. 3 for different values of \(\delta \) and assumed values of other parameters. It may be noted that for \(\delta =1,~1.1,~1.3,~1.5,~2\) the model exhibits a smooth transition from early deceleration era to the current acceleration era of the universe, whereas for \(\delta =0\) the model remains in the accelerated phase. Moreover, the transition redshift \(z_t\) is decreasing as \(\delta \) increases. We observed that the transition redshift \(z_t\) from a deceleration phase to an accelerated universe lies within the interval \(0.57<z_t<0.77\). This is in accordance with the recent cosmological observations [101, 102].

Fig. 4
figure 4

Plot of \(\omega _{de}-\omega _{de}^\prime \) plane of non-interacting THDE model for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image

The \(\omega _{de}-\omega _{de}^\prime \) (where \(\prime \) denotes differentiation with respect to \(\ln ~a\)) plane explains the accelerated expansion regions of the universe. Caldwell and Linder [103] have firstly proposed this plane for analyzing the quintessence scalar field. Two different classes have been characterized from this plane known as thawing region (\(\omega _{de}^\prime >0\) when \(\omega _{de}<0\)) and freezing region (\(\omega _{de}^\prime <0\) when \(\omega _{de}<0\)) on the \(\omega _{de}^\prime -\omega _{de}\) plane. Here, for non-interacting model, we have developed \(\omega _{de}-\omega _{de}^\prime \) plane for different values of \(\delta \) as given in Fig. 4. From Fig. 4, it can be observed that \(\omega _{de}-\omega _{de}^\prime \) plane corresponds to freezing region for all the values of \(\delta =0\), 1, 1.1, 1.3, 1.5. Observational data says that the expansion of the universe is comparatively more accelerating in freezing region. Also, for \(\delta =2\) the model coincides with the \(\varLambda CDM\) limit \((\omega _{de},\omega _{de}^\prime )=(-1,0)\). Hence, the behavior of \(\omega _{de}-\omega _{de}^\prime \) plane is consistent with the present day observations.

In order to check the stability against small perturbations of our non interacting THDE model in this scenario, we obtain the squared speed of sound. Sign of square of speed of sound plays a crucial role, i.e., its negativity (\(v_s^2<0\)) represents instability and vice versa. It can be defined as follows

$$\begin{aligned} v_s^2=\frac{{\dot{p}}_{de}}{{\dot{\rho }}_{de}}. \end{aligned}$$
(31)

By differentiating the relation \(\omega _{de}=\frac{p_{de}}{\rho _{de}}\) with respect to time t and dividing by \({\dot{\rho }}_{de}\), we get

$$\begin{aligned} v_s^2=\omega _{de}+\frac{\rho _{de}}{{\dot{\rho }}_{de}}~{\dot{\omega }}_{de} \end{aligned}$$
(32)
Fig. 5
figure 5

Plot of squared sound speed of non-interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image

Now, using Eq. (22) in above Eq. (32), we get

$$\begin{aligned} v_s^2=\omega _{de}+\frac{{\dot{\omega }}_{de}}{H\left( \frac{\beta _2a}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}+(4-2\delta )\frac{{\dot{H}}}{H^2}\right) } \end{aligned}$$
(33)

where \(\omega _{de}\) and \(\frac{{\dot{H}}}{H^2}\) are, respectively, given in Eqs. (25) and (28).

In the present scenario, we develop the squared speed of sound trajectories in terms redshift z as shown in Fig. 5 using same values of the constants in the model. For all the values of \(\delta \), we can observe from Fig. 5 that initially \(v_s^2\) exhibits a decreasing behavior but with positive sign. This shows that the model stable at early epochs. Later, all the curves related to the parameter \(\delta \) exhibit negative behavior. Thus for the non-interacting THDE model is found the stability at initial epoch and then becomes unstable in later epochs.

2.2 Interacting model

Here, we consider the interaction between the two fluids. Since the nature of both DE and dark matter is still unknown, there is no physical argument to exclude the possible interaction between them. Some observational evidences of the interaction in dark sectors have been found recently [104, 105]. Abdalla et al. [106, 107] have investigated the signature of interaction between DE and dark matter by using optical, X-ray and weak lensing data from the relaxed galaxy clusters. So, it is reasonable to assume the interaction between DE and dark matter in cosmology.

For this purpose we can write the energy conservation equations as

$$\begin{aligned} {\dot{\rho }}_{de}+3H\rho _{de}(1+\omega _{de})= & {} -Q, \end{aligned}$$
(34)
$$\begin{aligned} {\dot{\rho }}_{m}+3H\rho _m= & {} Q. \end{aligned}$$
(35)

where the quantity Q represents interaction between DE components. From the Eqs. (34) and (35) we can say that the total energy is conserved i.e., \({\dot{\rho }}_{tot}+3H(\omega _{tot}+1)\rho _{tot}=0\), where \(\rho _{tot}=\rho _m + \rho _{de}\). Since there is no natural information from fundamental physics on the interaction term Q, one can only study it to a phenomenological level. Various forms of interaction term extensively considered in literature include \(Q=3cH\rho _m\), \(Q=3cH\rho _{de}\) and \(Q=3cH(\rho _m+\rho _{de})\). Here, c is a coupling constant and positive c means that DE decays into dark matter, while negative c means dark matter decays into DE. Also, there are many other forms for interaction term considered in literature which are defined as \(Q=\gamma {\dot{\rho }}_i\) and \(Q=3cH\gamma \rho _i+\gamma {\dot{\rho }}_i\) where \(i={m, de, tot}\). A detailed analysis of linear interacting terms can be found in Izaquierdo and Pavon [108], Ferreira et al. [109], Sadeghi et al. [47, 110] and Wei [111]. Most of these interactions are either positive or negative and can not give the possibility to change their signs. Cai and Su [112] have proposed a new sign-changeable interaction by allowing \({\dot{\rho }}\) and deceleration parameterq into interacting term as \(Q=q(\gamma {\dot{\rho }}+3cH\rho _i)\) where \(i={m, de, tot}\). Here \(\gamma \) and c dimensionless constants and deceleration parameter is defined in Eq. (5). Recently, Sadri et al. [113] have compared different phenomenological linear as well as non-linear interaction cases in the framework of the holographic Ricci DE model and they found that the linear interaction are the best cases among the others. Inspired by the works of Pavon and Zimdahl [114], Sadeghi et al. [115], Honarvaryan et al. [116], Zadeh et al. [82] and Sadri [84], in this work, we consider the linear interacting term as \(Q=3c^2H(\rho _m+\rho _{de})\). It is worthwhile to mention here that this work can be extended using sign-changeable interaction term, which will be done in other forthcoming articles.

Fig. 6
figure 6

Plot of density parameter of interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\), \(u=0.3\) and \(\delta = 1.3\)

Full size image

Now, combining the time derivative of Eq. (12) with Eqs. (16)–(19), (34) and (35), we easily get

$$\begin{aligned} \frac{{\dot{H}}}{H^2}= & {} -\biggl \{\frac{3\beta _2 a(2-\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}\nonumber \\&-9\varOmega _{de}\left( u(c^2-1)+c^2\right) \biggr \}\nonumber \\&\times \biggl \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\biggr \}^{-1} \end{aligned}$$
(36)

In view of Eq. (26), from Eq. (36) we obtain

$$\begin{aligned} \varOmega _{de}^\prime= & {} 2\varOmega _{de}(\delta -1)\bigg \{\frac{3\beta _2 a(2-\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}\nonumber \\&-9\varOmega _{de}\left( u(c^2-1)+c^2\right) \bigg \}\nonumber \\&\times \bigg \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\bigg \}^{-1} \end{aligned}$$
(37)
Fig. 7
figure 7

Plot of EoS parameter of interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\), \(u=0.3\) and \(\delta = 1.3\)

Full size image

We solve the differential equation (37) for dimensionless density parameter \(\varOmega _{de}\) in interacting THDE model and plot it against redshift z as shown in Fig. 6. Here, we have analyzed all the dynamical parameters through both \(\delta \) and coupling coefficient \(c^2\). Remaining constant parameters are same as used in previous section. It can be seen that the density parameter \(\varOmega _{de}\) shows decreasing behavior for all values of \(\delta \) except for \(\delta =0\). Also, we observe that increase of the coupling coefficient \(c^2\) and parameter \(\delta \) cause the \(\varOmega _{de}\) increases but, the effect of coupling coefficient vanishes near past and \(\varOmega _{de}\) becomes more steeper with the increase of \(\delta \).

From Eqs. (18), (34) and (36), we find that

$$\begin{aligned} \omega _{de}= & {} -1-c^2(u+1)-\frac{\beta _2a}{3(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{2(\delta -2)}{3}\Bigg \{\bigg \{\frac{3\beta _2 a(2-\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2 a)]^3}\nonumber \\&-9\varOmega _{de}\left( u(c^2-1)+c^2\right) \bigg \}\nonumber \\&\times \bigg \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\bigg \}^{-1}\Bigg \} \end{aligned}$$
(38)

We can get EoS parameter behavior from it’s plots which are shown in Fig. 7, for various values of \(c^2\) and \(\delta \). It can be observed from Fig. 7 that the EoS parameter starts from quintessence phase and turns towards phantom region by crossing vacuum dominated era (phantom divide line \(\omega _{de}=-1\)) of the universe for all the values of \(\delta \) and \(c^2\). At late times, it can also be observed that the EoS parameter attains high phantom region with the increase of coupling coefficient \(c^2\) but, attains a constant value in phantom region for various values of \(\delta \). For \(\delta =0\) the model varies from matter dominated phase to dark energy phase, whereas for \(\delta =2\) the model remains in the phantom region (\(\omega _{de}<-1\)). It is interesting to note here that as \(\delta \) and \(c^2\) increases, there is a delay in transition of the model from quintessence to phantom phase.

Fig. 8
figure 8

Plot of deceleration parameter of of interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image
Fig. 9
figure 9

Plot of \(\omega _{de}-\omega _{de}^\prime \) plane of interacting THDE model for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\), \(u=0.3\) and \(c^2=0.15\)

Full size image

The deceleration parameter is given by

$$\begin{aligned} q= & {} -1+\bigg \{\frac{3\beta _2 a(2-\varOmega _{de})}{(\beta _1+\beta _2 a)~\ln (\beta _1+\beta _2 a)}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&-\frac{3\beta _2^2 a^2}{(\beta _1+\beta _2 a)^2~\ln (\beta _1+\beta _2 a)}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^2~[\ln (\beta _1+\beta _2 a)]^2}\nonumber \\&+\frac{w\beta _2^3a^3}{(\beta _1+\beta _2 a)^3~[\ln (\beta _1+\beta _2a)]^3}\nonumber \\&-9\varOmega _{de}\left( u(c^2-1)+c^2\right) \bigg \}\nonumber \\&\times \bigg \{6+6(\delta -2)\varOmega _{de}+\frac{6\beta _2 a}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]}\nonumber \\&-\frac{w\beta _2^2a^2}{(\beta _1+\beta _2 a)~[\ln (\beta _1+\beta _2 a)]^2}\bigg \}^{-1} \end{aligned}$$
(39)

Figure 8 represents the behavior of deceleration parameter for interacting THDE model for different values of \(\delta \) and \(c^2\). In both cases, it shows that the deceleration parameter starts from decelerating phase, then goes towards accelerating phase and finally approaches to \(q=-1\). Also, the transition from deceleration to acceleration occurred between \(z=0.6\) to 0.8. It may be noted here that the transition redshift of the model decreases as \(\delta \) increases and it is quite opposite in case of coupling coefficient \(c^2\), i.e., the transition is delayed as coupling coefficient increases.

The \(\omega _{de}-\omega _{de}^\prime \) plane is used to present the dynamical property of dark energy models, where \(\omega _{de}^\prime \) is the evolutionary form of \(\omega _{de}\) with prime indicates derivative with respect to \(\ln \) a. The \(\omega _{de}-\omega _{de}^\prime \) plane for the constructed interacting THDE model is developed for various values of \(\delta \) and \(c^2\) as shown in Fig. 9. We observed from these plots that \(\omega _{de}-\omega _{de}^\prime \) plane corresponds to the freezing region only. It is concluded that \(\omega _{de}-\omega _{de}^\prime \) plane analysis for the present scenario gives consistent results with the accelerated expansion of the universe.

Fig. 10
figure 10

Plot of squared sound speed of interacting THDE model versus redshift z for \(\beta _1=7.1\), \(\beta _2=2.15\), \(w=1000\), \(\varOmega ^0_{de}=0.73\) and \(u=0.3\)

Full size image

Using Eqs. (36) and (38) in Eq. (33) we can obtain the expression for squared sound speed \(v_s^2\) of interacting THDE model in BD theory. In Fig. 10, we plot its obtained expression in terms of redshift z for same values of the other constants involved in the equations as for above figures. It can be seen that initially \(v_s^2\) shows a decreasing behavior but with positive sign which shows the stable model. As the model evolves, these trajectories bear a negative behavior. Thus for the interacting THDE model in BD theory is stable at initial epoch and then exhibits instability at late times.

3 Conclusions

The current scenario of accelerated expansion of the universe has become more fascinating with the passage of time. In order to address this problem, different approaches have been adopted through lot of dynamical DE models and modified theories of gravity. In this work, we investigate this phenomenon by assuming the Tsallis Holographic DE within Brans–Dicke scalar–tensor theory of gravity by taking the BD scalar field as a logarithmic function of average scale factor. We have studied both scenarios i.e., interaction and non-interaction between DE and pressure less DM. We have summarized our results as follows:

\(\bullet \) In order to explore the viability of constructed THDE models, we have obtained different dynamical cosmological parameters. In their evolution (non-interacting and interacting) the THDE parameter \(\delta \) and coupling coefficient \(c^2\) play a significant role. Hence, we have explored the dynamics of physical parameters in terms of redshift z for the values of \(\delta =0,~1,~1.1,~1.3,\)  1.5,  2 and \(c^2=0.05,~0.15,~0.25\). We solve THDE density parameters (Eqs. (27) and (38)), in both non-interacting and interacting models, numerically corresponding to redshift z and their outputs are plotted in Figs. 1 and 6. It can be seen that the \(\varOmega _{de}\) is positive and decreasing function throughout the evolution, and at present epoch it approaches to 0.73 for all the three values of \(\delta \) and \(c^2\). Also, \(\varOmega _{de}\) increases with \(\delta \) and \(c^2\). But for \(\delta =0\), the energy density parameter increases. Planck’s observations (2014, 2018) [99, 100] have put the following constraints on the DE density parameter \(\varOmega _{de}=0.717^{+0.023}_{-0.020}\) (Planck \(+\) WP), \(0.717^{+0.028}_{-0.024}\) (WMAP-9), \(0.679 \pm 0.013\) (TT \(+\) lowE), \(0.699 \pm 0.012\) (TE \(+\) lowE), \(0.711^{+0.033}_{-0.026}\) (EE \(+\) lowE), \(0.6834 \pm 0.0084\) (TT, TE, EE \(+\) lowE), \(0.6847 \pm 0.0073\) (TT, TE, EE \(+\) lowE \(+\) lensing), \(0.6889\pm 0.0056\) (TT, TE, EE \(+\) lowE \(+\) lensing \(+\) BAO). We observed that, in both non-interacting and interacting, the THDE density parameter meets the above mentioned limits which shows that our results are consistent.

\(\bullet \) The behaviors of EoS parameter versus redshift for non-interacting and interacting THDE models are respectively given in Figs. 2 and 7. We observe that the non-interacting model starts from matter dominated era, then goes towards quintessence DE era and finally approaches to vacuum DE era for all three values of \(\delta \) (Fig. 2). For \(\delta =2\), the EoS parameter becomes \(-1\) i.e., cosmological constant, which is in agreement with the results obtained by Saridakis et al. [78]. Also, we observed that as \(\delta \) decreases the transition from matter dominated phase to dark energy phase is delayed is delayed considerably. It may be noted that the non-interacting model never crosses the PDL (\(\omega _{de}=-1\)). Also, the EoS parameter shows translation from matter dominated era and goes towards phantom region by crossing the phantom divide line (Fig. 7). In interacting THDE case, the model starts from accelerated phase and then goes towards phantom DE era by crossing the quintessence era as well as PDL (i.e. vacuum DE) (Fig. 7). This type of behavior corresponds to quintom. It can also be observed that as \(\delta \) and \(c^2\) increases the transition from quintessence to phantom phase of the interacting THDE model is delayed. The present Planck collaboration data (2018) [100] gives the limits for EoS parameter as

$$\begin{aligned} \omega _{de}= & {} -1.56^{+0.60}_{-0.48} ~ \text {(Planck + TT + lowE)} \\ \omega _{de}= & {} -1.58^{+0.52}_{-0.41} ~ \text {(Planck + TT, TE, EE + lowE)} \\ \omega _{de}= & {} -1.57^{+0.50}_{-0.40} ~ \text {(Planck + TT,\;TE,\;EE + lowE + lensing)} \\ \omega _{de}= & {} -1.04^{+0.10}_{-0.10} ~ \\&\text {(Planck\, +\, TT,\;TE,\;EE\,+\,lowE\,+\,lensing\,+\,BAO)}. \end{aligned}$$

It can be seen that the trajectories of EoS parameter of both interacting and non-interacting THDE models (Figs. 2, 7) coincide with this Planck collaboration (2018) results.

\(\bullet \) The trajectory of deceleration parameter versus redshift for non-interacting and interacting models are shown in Figs. 3 and 8, respectively. It can be seen that, in both the cases, the model exhibits a smooth transition from early deceleration era to the current acceleration era of the universe. Also, the transition redshift \(z_t\) from a deceleration to acceleration lies within the interval \(0.58<z_t<0.8\). Also, the observational data from \(H(z)+SN~Ia\) for \(\varLambda \)CDM model given a range for transition redshift as \(z_t=0.682 \pm 0.082\) and \(q\rightarrow -1\) as \(z\rightarrow -1\). Hence, we can say that the above results are in accordance with the recent cosmological observations [101, 102].

\(\bullet \) In this present scenario, we develop the squared speed of sound \(v_s^2\) trajectories for both non-interacting and interacting THDE models (Figs. 5,  10). In both the models, \(v_s^2\) varies in positive region initially which shows the stability of the models. However, it has become negative within short interval of time and remains negative forever exhibits instability of the models later epochs. Many authors [60, 82] in the literature have shown that the non-interacting or interacting THDE models with different IR cutoff’s are unstable. However, it is interesting to mention that for increasing value of \(\delta \) and decreasing value of interaction parameter \(c^2\), we expect stability of the both models in near future. The \(\omega _{de}-\omega _{de}^\prime \) plane for the non-interacting and interacting THDE models is also developed as given in Figs. 4 and 9. It may be observed from these trajectories that \(\omega _{de}-\omega _{de}^\prime \) plane corresponds to freezing only. Many researchers have concluded that the expansion of the universe is comparatively more accelerating in freezing region. Also, in both the models \(\omega _{de}-\omega _{de}^\prime \) plane meets the observational data from Planck collaboration (Planck \(+\) WP \(+\) BAO, [99]) which is given by \(-0.89\le \omega _{de}\le -1.38\) and \(\omega _{de}^\prime <1.32\).

Now it will be interesting to compare our THDE models in BD theory with logarithmic scalar field with the other THDE models in literature with regard to the energy density (\(\varOmega _{de}\)), EoS (\(\omega _{de}\)), deceleration (q) parameters and the stability analysis \(v_s^2\). Zadeh et al. [82] has explored the effects of different IR cutoffs on the properties of THDE model. It seems that our findings are coincide with the results of Ghaffari et al. [98], Zadeh et al. [82] and Tavayef et al. [60]. Ghaffari et al. [87] have investigated THDE in BD theory, while Jawad et al. [88] have studied THDE in modified BD theory with scalar field as power function of average scale factor. Ghaffari et al. [98] have shown that the density parameter is increasing function and converges to \(\varOmega _{de}=1\) at late times, but in our models (both non-interacting and interacting) the density parameter has opposite behavior, i.e., \(\varOmega _{de}\) is decreasing positive function and finally tends to a positive constant value. This is due to because of logarithmic BD scalar field in our model. However, the present value of \(\varOmega _{de}\) is same as in [87] and also coincide with the observational results. Our THDE non-interacting model never crosses the PDL but, in the work of Jawad et al. [88] the EoS parameter behaves like quintom. We, also, observed from the findings of Ghaffari et al. [87] that the transition redshift of each dynamical parameter is same, but in our models it varies through \(\delta \) and \(c^2\). Our both THDE models are unstable for any value of \(\delta \) or coupling constant \(c^2\), a result the same as that of the models in [87, 88]. The above results leads to the conclusion that our THDE models in BD theory with logarithmic form of scalar field are in good agreement with the observational data and also we hope that the above investigations will help to have a deep insight into the behavior of THDE universe in BD cosmology.

It is interesting to mention here that, for \(\delta =0\), the dynamical behavior of properties such as \(\varOmega _{de}\) (Figs. 1, 6), \(\omega _{de}-\omega _{de}^\prime \) plane (Figs. 4, 9) and squared sound speed (Fig. 10) have some peculiar behavior when compared with the other values of \(\delta \). For \(\delta =0\), it is clear from the Figs. 1 and 6 that the density parameter \(\varOmega _{de}\) positive and increasing function which is quite opposite to the behavior of \(\varOmega _{de}\) when \(\delta \) is non-zero. Also, for \(\delta =0\), the interacting THDE model is unstable near present epoch whereas for other non-zero values of \(\delta \) the model becomes unstable near past. The interacting THDE model varies in DE region only for non-zero values of \(\delta \) whereas for \(\delta =0\) the interacting model varies in both matter dominated and DE regions.