Dynamical analysis and statefinder of Barrow holographic dark energy

Based on the holographic principle and the Barrow entropy, Barrow holographic dark energy had been proposed. In order to analyze the stability and the evolution of Barrow holographic dark energy, we, in this paper, apply the dynamical analysis and statefinder methods to Barrow holographic dark energy with different IR cutoff and interacting terms. In the case of using Hubble horizon as IR cutoff with the interacting term $Q=\frac{\lambda}{H}\rho_{m}\rho_{D}$, we find this model is stable and can be used to describe the whole evolution of the universe when the energy transfers from the pressureless matter to the Barrow holographic dark energy. When the dynamical analysis method is applied to this stable model, an attractor corresponding to an accelerated expansion epoch exists and this attractor can behave as the cosmological constant. Furthermore, the coincidence problem can be solved in this case. Then, after using the statefinder analysis method to this model, we find this model can be discriminated from the standard $\Lambda$CDM model. Finally, we have discussed the turning point of Hubble diagram in Barrow holographic dark energy and find the turning point does not exist in this model.


I. INTRODUCTION
Observations [1][2][3][4][5][6] indicate that the universe is currently undergoing a phase of accelerated expansion. Nevertheless, Einstein's General Relativity with only radiation and matter cannot result in an accelerated expansion of the universe. In order to explain this observational result, a mysterious matter called dark energy was introduced to explain this accelerated phase. Since dark energy has not been detected directly, there is no reason to assume dark energy resembles known forms of matter or energy, and it becomes one of the biggest mysteries in modern cosmology. The simplest dark energy model is the cosmological constant (ΛCDM) in which the constant vacuum energy density can drive the accelerated expansion. Although it is well in agreement with the current observations [7], it suffers from the fine-tuning problem and the coincidence problem.
Based on the generic holographic principle and using the Barrow entropy [76] instead of the Bekenstein-Hawking entropy, Barrow holographic dark energy (BHDE) was proposed by considering the future event horizon as IR cutoff [77,78]. When the IR cutoff was chosen as the Hubble horizon, a new BHDE was proposed [79]. Then, the thermodynamics of this new BHDE with an interaction term was studied [80]. However, this new BHDE is unstable since its squared sound speed is negative [79]. It is noted that an interaction term between dark energy and dark matter can be used to solve the classical instability in HDE [51] and THDE [63], and this interaction term can also alleviate the cosmic coincidence problem [16,50,51]. So, in BHDE models, it is unclear whether the classical instability and the cosmic coincidence problem can be solved when an interaction term is taken into account.
When a dark energy model is proposed, it should fit the conventional standard cosmology as well as explain the current accelerated expansion. For a viable cosmological model, it should describe the whole evolution history of the universe. Namely, the universe stems from a radiation dominated epoch, and then enters into a matter dominated epoch to enable the formation of largescale structures, and eventually evolves into a dark energy dominated epoch. To analyze the evolution of the universe, the dynamical analysis method is introduced, which is an excellent method to investigate the qualitative behavior of a given cosmological model and can avoid the difficulty in solving non-linear cosmological equations. In this method, the dynamics of the universe can be described by analyzing the behavior of critical points of the dynamical system of a model, and the critical points always indicate the main evolution epoch of the universe. The stable point is used to describe a late time epoch dominated by the dark energy, the saddle point can denote the matter dominated epoch, and the unstable point corresponds to the early radiation dominated epoch. This method was used with great success in analyzing the evolution of the universe within HDE models [63,68,[81][82][83][84][85]. These positive results in HDE lead to some interesting questions in BHDE and motivate us to study BHDE: Whether a stable BHDE model can describe the whole evolution of the universe. How to discriminate the BHDE model from the standard ΛCDM model if the BHDE model describes the cosmic evolution successfully? The main task of this paper is to answer these questions.
Furthermore, a turning point of Hubble diagram in HDE was studied in Ref. [86] which shows that HDE model may be at odds with the cosmological principle. Here, the free parameter region 0 < c < 1, which was constrained by the current data, was considered. For 0.5 < c < 1, this turning point exists in the future, while it occurs in the past and is observable for c < 0.5. In BHDE with the future event horizon as IR cutoff, the combination H(z) + SN Ia place constraints on the free parameter C = 3.421 +1.753 −1.611 and the parameter of Barrow entropy ∆ = 0.094 +0.094 −0.101 [78]. It is interesting to discuss whether this turning point exists in BHDE with the future event horizon as IR cutoff when the constraints C = 3.421 +1.753 −1.611 and ∆ = 0.094 +0.094 −0.101 are considered. Regarding that BHDE may offer an interesting framework to study phenomenology beyond to HDE, the main goal of this paper is to discuss six BHDE models: we first consider BHDE with different IR cutoffs, and then, two kinds of interaction terms are used. The analysis is performed in four respects: Firstly, by analyzing the squared sound speed, we study the stability of BHDE models. Secondly, we use the dynamical analysis method to analyze the phase space behavior of the stable BHDE models and discuss the cosmic coincidence problem. Thirdly, in order to discriminate the stable BHDE model from the standard ΛCDM model, we per-form the statefinder diagnostic by depicting the evolution trajectories of statefinder parameters. Finally, we plot the evolution curves of H(z) and study the turning point of Hubble diagram in BHDE models.
The paper is organized as follows. In Section II, we investigate the evolution and stability of the universe in BHDE with different IR cutoff and interacting terms. In Section III, we use the dynamical analysis method to analyze the phase space behavior of the stable BHDE model. In Section IV, the statefinder diagnostic pairs are used to discriminate the stable BHDE model from the standard ΛCDM model. In Section V, we discuss the turning point of Hubble diagram in BHDE. Finally, our main conclusions are presented in Section VI.

II. THE UNIVERSE EVOLUTION
Using the Barrow entropy [76], the energy density of BHDE is given as [77] ρ D = CL △−2 . (1) Here, C is a parameter with the dimension [L] −2−△ and ∆ satisfies the relation 0 ≤ ∆ ≤ 1. For the case ∆ = 0, ρ D provides a smooth spacetime structure, i.e. the standard holographic dark energy, and ∆ = 1 corresponding to the most intricate structure. We consider a homogeneous and isotropic flat Friedmann-Robertson-Walker universe where a(t) is the scale factor with t being cosmic time, and setting κ 2 = 8πG. The Friedmann equation is given as where ρ r , ρ m and ρ D denote the energy density of radiation, pressureless matter and BHDE, respectively. The radiation, pressureless matter and BHDE conservation equations take the forṁ in which w D = pD ρD denotes the equation of state parameter of BHDE and Q represents the energy exchange between pressureless matter and BHDE. For Q > 0, energy transfers from BHDE to pressureless matter, and energy transfers from pressureless matter to BHDE for Q < 0. In this paper, we consider two interacting cases Q = H(αρ m +βρ D ) [63,83,87] and Q = λ H ρ m ρ D [87,88]. Here, α, β and λ are coupling constants.
To discuss the stability of BHDE model, we need to analyze its squared sound speed For v 2 s > 0, this model can be stable against perturbations. Otherwise, it is unstable.
To find a stable and suitable model for BHDE, we will analyze six models in the following: (i)Model I: Non-interacting BHDE with future event horizon as IR cutoff (BHDEF); (ii)Model II: Interacting BHDE with future event horizon as IR cutoff and Q = H(αρ m + βρ D ) (IBHDEFA); (iii)Model III: Interacting BHDE with future event horizon as IR cutoff and Q = λ H ρ m ρ D (IBHDEFL); (iv)Model IV: Non-interacting BHDE with Hubble horizon as IR cutoff (BHDEH); (v)Model V: Interacting BHDE with Hubble horizon as IR cutoff and Q = H(αρ m + βρ D ) (IBHDEHA); (vi)Model VI: Interacting BHDE with Hubble horizon as IR cutoff and Q = λ H ρ m ρ D (IBHDEHL).

A. BHDE with future event horizon
Considering the future event horizon as IR cutoff, the energy density of BHDE takes the form [77] with which satisfies the relationṘ h = HR h − 1. For the case ∆ = 0 with C = 3c 2 , it reduces to the standard holographic dark energy.
In the following, we assume the present scale factor a 0 = 1 so that the redshift z satisfies the relation z = 1 a − 1. Since the combination H(z) + SN Ia constrain C = 3.421 +1.753 −1.611 and ∆ = 0.094 +0.094 −0.101 [78], we consider C = 3.5 and ∆ = 0.1 within BHDE with future event horizon as IR cutoff.

Non-interacting
When there is no energy exchange between the pressureless matter and BHDE, we obtain Q = 0 and σ = 0. By differentiating Eq. (14) and using Eq. (6), one can obtain the equation of state parameter of BHDE with F = 1 HR h and F ′ satisfies Using Eqs. (14) and (16), the squared sound speed can be written as Substituting Eq. (17) into the above relation, v 2 s > 0 leads to These results mean the condition v 2 s > 0 cannot be satisfied during the whole evolution of the universe. So, this model is unstable.
Throughout this paper, we choose Ω   conditions, by solving Eqs. (11), (12) and (17) numerically, we get the evolution curves of Ω D , Ω m , w D and q which are plotted in Figs. (1) and (2). And the black dotted line in Figs. (1) and (2) represents the situation of ∆ = 0 which corresponds to the standard entropy. Form Fig. (1), we can see Ω D → 1 and Ω m → 0 at the late time, and the universe can be dominated by BHDE at the late time evolution. The left panel of Fig. (2) shows that the BHDE behave as the quintessence and w D can cross the phantom line. The right panel of Fig. (2) indicates a suitable range for the transition redshift is obtainable. But this model is unstable since v 2 s > 0 cannot be satisfied, and some evolution curves of v 2 s are depicted in Fig. (3).

Interacting with Q = H(αρm + βρD)
When the energy exchange Q = H(αρ m + βρ D ) between the pressureless matter and BHDE exists, the equation of state parameter w D has the form And the squared sound speed v 2 s is given as Here, F ′ is given in Eq. (17). Since the expression of v 2 s is too complicated, we cannot get the analytical solution for v 2 s > 0 and then resort to numerical solution. Using the initial conditions with ∆ = 0.1, C = 3.5 and solving Eqs. (11), (12) and (17), we obtain the evolution curves for Ω D , Ω m , w D , q and v 2 s which are plotted in Figs. (4), (5) and (6). The specific case of ∆ = 0 is plotted by the black dotted line, and the values of α, β and C are 0.1, 0.1 and 0.8, respectively.
It can be seen from Fig. (4) where Ω D is larger than 1 and Ω m is less than 0 for a negative β. And for the positive β, the universe can be dominated by BHDE at the late time. The evolution curves of w D and q is shown in Fig. (5). The left panel of Fig. (5) shows that BHDE can behave as quintessence and cross the phantom line for a negative α, and the right one shows that the late time acceleration can be achieved and the transition redshift is obtainable. These results mean that α < 0 and β > 0 can be used to describe the evolution of the universe. However, the result from Fig. (6) means this model is unstable since the evolution curve of v 2 s is not always positive during the whole evolution epoch.

Interacting with Q = λ H ρmρD
When the energy exchange Q = λ H ρ m ρ D between the pressureless matter and BHDE is adopted, w D and v 2 s are written as and Here, F ′ is given by Eq. (17). Then, solving Eqs. (11), (12) and (17) with ∆ = 0.1 and C = 3.5, we get the evolution curves of Ω D , Ω m , w D , q and v 2 s which are plotted in Figs. (7), (8) and (9). The black dotted line in Figs. (7), (8) and (9) denotes the specific case of ∆ = 0 with λ = 0.1 and C = 0.8. From these figures, one can see that this model can produce suitable results if proper values are chosen, but it is unstable.

B. BHDE with Hubble horizon
Using the Hubble horizon as IR cutoff, the energy density of BHDE leads to For ∆ = 2(δ − 1), the energy density of BHDE has the same form as that of Tsallis holographic dark energy (THDE) in mathematics [56]. Different from THDE, the parameter ∆ in BHDE is limited to 0 ≤ ∆ ≤ 1.

Non-interacting
For the case Q = 0, differentiating Eq. (23) and using Eq. (6), the equation of state parameter w D and the squared sound speed v 2 s lead to and Using the initial conditions and solving Eqs. (11) and (12), we obtain the evolution curves of Ω D , Ω m , w D , q and v 2 s which have been plotted in Figs. (10), (11) and (12). In Figs. (10), (11) and (12), the case ∆ = 0 is depicted by the black dotted line. It can be seen from above figures that the current acceleration can be realized and the whole evolution of universe can described by this model. The universe can evolve into the era depicted by the standard ΛCDM model since w D approach to −1 at the late time, and a suitable range for the transition      redshift is obtainable for the large value of ∆. For BHDE with Hubble horizon as IR cutoff, the observable Hubble data favor a large value of ∆ which can be seen in the right panel of Fig. (27). Thus, for BHDE with Hubble horizon as IR cutoff, we set ∆ = 0.8. The right panel of Fig. (12) which is plotted for ∆ = 0.8 shows that this model can describe the evolution of universe. However, this model is also unstable since the evolution curves of v 2 s in the left panel of Fig. (12) means v 2 s < 0.
2. Interacting with Q = H(αρm + βρD) For the situation Q = H(αρ m + βρ D ), we can write w D and v 2 s as follows and In Figs. (13), (14) and (15), the evolution curves of Ω D , Ω m , w D , q and v 2 s have been plotted with the initial conditions and ∆ = 0.8. Here, we have plotted the case ∆ = 0 with α = 0.7 and β = 0.2 by the black dotted line. From Fig. (13), we can see that Ω D is larger than 1 and the value of Ω m becomes negative at the late time era for the negative β. The left panel of Fig. (15) shows that this model can be stable for the negative α. So, a negative α with a positive β can realize the late time acceleration and ensure this model is stable. Unfortunately, the result from the right panel of Fig. (15) shows that Ω r dominates the evolution of the universe at the high redshift region z > 2 × 10 11 in which a approaches to 0.

Interacting with Q = λ H ρmρD
For the case Q = λ H ρ m ρ D , after some tediously calculations, w D and v 2 s lead to and The evolution curves of Ω D , Ω m , w D , q and v 2 s are plotted in Figs. (16), (17) and (18) with the initial conditions and ∆ = 0.8. The black dotted line in Figs. (16), (17) and (18) represents the case ∆ = 0 with λ = 0.5. From these figures, we can see that the current acceleration and the whole evolution of universe can be realized in this model, i.e. the universe evolves from the radiation dominated era to the pressureless matter dominated era, and then enters into the BHDE dominated era. Since both w D and q approach to −1, BHDE behaves as the cosmological constant at the late time evolution, and the universe will eventually evolve into an epoch described by the standard ΛCDM model. In addition, the left panel of (18) shows that this model is stable for some negative coupling constant λ, while it is unstable for a positive coupling constant. These results means when the energy transfers from the pressureless matter to BHDE, this model is stable. Thus, this model can be stable under some specific conditions.

III. DYNAMICAL ANALYSIS OF THESE STABLE MODELS WITH HUBBLE HORIZON AS IR CUTOFF
In previous section, we have discussed the cosmic evolution and stability of BHDE models, and it is found that only IBHDEHA and IBHDEHL model can be stable against perturbations under some specific conditions. In the following, we will analyze the dynamical behaviour of these two stable models.
In order to investigate the complete asymptotic behaviour of IBHDEHA and IBHDEHL models, we will use the dynamical system techniques to study the dynamical behaviour of these models. To achieve this goal, following Ref. [63,87,[89][90][91], we can obtain the critical points by solving the autonomous system equations By linearizing the autonomous system equations, we obtain the corresponding first order differential equations. Then, the stability of critical points can be determined by the eigenvalues of the coefficient matrix of the first order differential equations. In the case for taking the Hubble horizon as IR cutoff, the autonomous system consists of Eqs. (11) and (12). Then, solving Ω ′ m = Ω ′ D = 0, we obtain three critical points given in Table I. The value of Ω r can be obtained by Eq. ( 8). For q < 0, the corresponding critical point    represents an acceleration phase, otherwise it is a deceleration one.
Thus, point A 1 denotes the radiation dominated deceleration epoch, point A 2 represents the pressureless matter dominated deceleration epoch, and point A 3 is an acceleration epoch determined by the value of α and β. For the case β = 0, point A 3 becomes (0, 1) and w D = −1 is obtained. Thus, A 3 can behave as the cosmological constant Λ dominated acceleration era under the condition β = 0.
To analyze the stability of these critical points, we linearize this autonomous system and then obtain the eigenvalues of the corresponding critical points. These results are also shown in Table I. Under the constraints conditions 0 < ∆ < 1 and −1 < α ≤ 1, point A 1 is unstable and P 2 is a saddle point, the stability of A 3 is determined by α and β. The stability of these points indicates that the universe stems from the radiation dominated era (A 1 ) and then evolves into the pressureless matter dominated era (A 2 ), and eventually enters into the BHDE dominated late time acceleration epoch (A 3 ).
Since A 3 is a stable point with Ω m = β 3−α+β and Ω D = 3−α 3−α+β , the ratio of energy density r = ρm ρD = ΩM ΩD can approach a constant, which means the coincidence problem can be solved. An example of this case is shown in Fig. (19). The behavior of attractor A 3 is shown in the left panel of Fig. (19), and the trajectories are shown in the right one. For the case β = 0, the universe evolves from the radiation dominated era (A 1 ) into the pressureless matter dominated era (A 2 ), and eventually enters the BHDE dominated late time acceleration epoch (A 3 ) which can behave as the cosmological constant Λ. The phase space and evolution trajectories of this case have been plotted in Fig. (20).
Although the IBHDEHA model is stable under some specific conditions and can describe the complete evolution epoch of the universe, Ω r cannot dominate the evolution of the universe in the suitable redshift regions. Some examples are plotted in Fig. (21) and the right panel of Fig. (15). So, this model cannot describe the whole evolution history of the universe.

B. Interacting with Q = λ H ρmρD
Solving Ω ′ m = Ω ′ D = 0, we obtain four critical points shown in Table II. From this table, it can be seen that point P 1 denotes a deceleration epoch dominated by radiation, point P 2 represents a deceleration epoch dominated by the pressureless matter, and point P 3 denotes the BHDE dominated acceleration era. Point P 4 is fully determined by the value of λ. It is easily seen that P 3 can behave as the cosmological constant Λ dominated acceleration era since w D = −1 and q = −1.
By linearizing this autonomous system, we obtain the eigenvalues of the corresponding critical points. These results are also shown in Table II. Under the constraint condition 0 < ∆ < 1, point P 1 is unstable and P 2 is a saddle point, the stability of P 3 and P 4 are determined by λ. For λ < 1, P 3 is a stable point, while a stable point P 4 requires λ > 1. Although P 3 and P 4 both behave as an attractor, they cannot be stable simultaneously.
Thus, for λ < 1, these results show that the universe evolves from the radiation dominated era (P 1 ) into the pressureless matter dominated era (P 2 ), and eventually enters the BHDE dominated late time acceleration epoch (P 3 ) which can behave as the cosmological constant Λ. Some examples of this situations are plotted in Fig. (22). The left panel describes the behavior of attractor P 3 . And the right panel depicts the phase space trajectories on the Ω m − Ω D plane in which all trajectories stem from P 1 , approach to P 2 , and eventually converge into P 3 . In this situation, the whole history of the universe can be described by this model. For 1 < λ < 2, the stability of these points indicates that the universe stems from the radiation dominated era (P 1 ) and then evolves close to the pressureless matter dominated era (P 2 ), and eventually enters into the BHDE dominated late time acceleration epoch (P 4 ) in which the coincidence problem can be solved. An example of this case is shown in Fig. (23). The behavior of attractor P 4 is shown in the left panel of Fig. (23), and the trajectories is shown in the right one. In this case, the coincidence problem is solved.
For λ ≥ 2, the universe will finally evolve into an accelerated epoch dominated by the pressureless matter rather than the BHDE. Although the coincidence problem can be solved in this case, the acceleration in this case is caused by the pressureless matter.

IV. STATEFINDER ANALYSIS
In previous section, we have discussed the dynamical evolution of the universe in IBHDEHA and IBHDEHL models. We find only IBHDEHL model can explain the current accelerated expansion and describe the whole evolution of the universe. And an attractor behaving as the cosmological constant Λ exists in this model. Then, how to discriminate this model from the standard ΛCDM model becomes another problem. To solve this problem, we use the statefinder diagnostic pairs {r, s}, which was introduced by Sahni et al. [92], to analyze this model.
The statefinder parameters r and s are defined as [92] The statefinder parameter is a geometrical diagnostic since it only depends on a. By differentiating Eq. (9), we can express the statefinder parameters r and s by Ω ′ m and Ω ′ D . Then, solving Eqs. (11) and (12) numerically, we can obtain the evolution curves of the universe in the parameter space r − s.

Label
Critical P oints (Ωm, ΩD) wD q Eigenvalues Conditions P oints Current value  Table I  In Fig. (24), we have plotted an example for the evolution curves of the statefinder diagnostic pairs {r, s}. In the left panel, we take some positive λ, and use the negative one in the right panel. The black dotted line in Fig. (24) represents HDE model. From Fig. (24), one can see that the statefinder pairs {r, s} start from the left side of the ΛCDM fixed point(0, 1), i.e. s < 0 and r > 1 which is the characteristic of Chaplygin gas, and then enter into the region s > 0, r < 1 where the trajectories of quintessence and phantom stay in these regions. After passing through these regions, they finally tend to the ΛCDM fixed point(0, 1). It is noticeable that the behaviour of the statefinder pairs {r, s} in these model is similar to the quintom dark energy [93] and interacting vacuum energy model [94]. In addition, for λ = −0.5, the statefinder pairs start from the Chaplygin gas region and then approach to the ΛCDM fixed point. For λ = −0.3, the present value of {r, s} is close to the ΛCDM fixed point, and it is not for the other cases. These trajectories in r − s plane demonstrate the contrasting behavior of this model and the standard ΛCDM model strikingly.
For the sake of complementarity, we have plotted the evolution trajectories for another statefinder diagram {r, q} in Fig. (25). From Fig. (25), one can see that the standard Λ model starts from the standard cold dark matter fixed point (0.5, 1), while the evolution curves for IBHDEHL have a smaller deviation from this fixed point. And both models can evolve into the de Sitter expansion fixed point (−1, 1) in the future.

V. TURNING POINT IN HUBBLE DIAGRAM
In HDE, it was found that there exists a turning point in the Hubble diagram H(z), which leads to HDE model conflicting with the cosmological paradigm [86]. This turning point is determined by the constant parameter c in HDE model, and this point occurs in the future for 0.5 ≤ c < 1, while it becomes observable for c ≤ 0.5.     This result was obtained in [86] and it was also shown in the left panel of Fig. (26). When the cosmological observational data are used to constrain the free parameter c in HDE model, the combination of different data provides different best fitting values of c, and some examples are shown in Table III where most results indicate that the turning point can be postponed to the future. Then, in order to discuss the turning point in BHDE models, we have plotted the evolution curves of Hubble parameter H in Fig. (26) and (27). The error bar in Fig. (26) and (27) represent the observational Hubble parameter data [46,96]. In the right panel of Fig. (  and C = 3.421 +1.753 −1.611 , which provide us an approach to solve the turning point in Hubble diagram. In BHDEA and BHDEF models, the constraints indicate that a small ∆ is favored by the current observations. However, in BHDEH model, the right panel of Fig. (27) indicates the observable Hubble data favor ∆ > 0.7, which is different from the values of ∆ in BHDEA and BHDEF models. To obtain a best fitting value of the parameters in BHDEA, BHDEF and BHDEH models, one require to combine more observational data to constrain the parameters. It  is interesting to compare BHDEA, BHDEF and BHDEH models and combine more current observations to constrain the parameters in these models, and this is what we shall do in the future.

VI. CONCLUSION
Based on the holographic principle and the Barrow entropy, which results from the modification of the black hole surface due to quantum-gravitational effects, a new HDE model named BHDE has been proposed. In this paper, by considering the future event horizon and the Hubble horizon as IR cutoff, we analyze the evolution and stability of BHDE models with different interaction terms. We find all of these BHDE models can explain the current accelerated expansion, but only IBHDEHA and IBHDEHL can be stable under some specific conditions in which the coupling constant requires to be negative. A negative coupling constant may result in Q < 0, which indicates that the pressureless matter is giving energy to BHDE. When we study the evolution of the universe in IBHDEHA and IBHDEHL models, we find only IBHDEHL model can describe the whole evolution of the universe, i.e. the universe stems from the radiation dominated epoch, passes through the pressureless matter dominated epoch, and eventually approaches the dark matter dominated acceleration epoch. Since both w D and q can approach to −1 at the late time evolution of the universe, the BHDE in IBHDEHL model can behave as the cosmological constant.
Then, we apply the dynamical analysis techniques to IBHDEHA and IBHDEHL models which are stable against perturbations. Although IBHDEHA model can describe the entire evolution epoch of the universe, Ω r cannot dominate the evolution of the universe in the corresponding redshift regions. Thus, only IBHDEHL model can describe the whole evolution of the universe. In this model, the results show that there exist three different critical points for λ < 1, namely the unstable point represents the radiation dominated epoch, the saddle point denotes the pressureless matter dominated epoch, and the stable point corresponds to the BHDE dominated epoch. For the stable point, since both w D and q equal to −1, this point can behave as the cosmological constant. The attractor behavior and evolution curves of the stable point show that IBHDEHL model can describe the thermal history of the universe. Thus, IBHDEHL model is a viable cosmological model. In addition, for 1 < λ < 2, the stable point can alleviate the coincidence problem since the ratio of the energy densities can approach to a constant.
In order to discriminate IBHDEHL model from the standard ΛCDM model, we apply the statefinder analysis method to IBHDEHL model. The statefinder diagrams show that the evolution curves of IBHDEHL model will eventually approach to the ΛCDM fixed point in r − s phase plane and the de Sitter expansion fixed point in r − q phase plane. Furthermore, the statefinder diagrams indicate that the interaction between the pressureless matter and BHDE has a significant effect on the evolution of the universe, and the IBHEDHL model can be distinguished from the standard ΛCDM model by this method. Finally, we discuss the turning point of Hubble diagram in BHDE models. For BHDEF model, the turning point is nonexistent under the current observational results. For BHDEH model, this point vanishes for ∆ > 0.7, and this point can also be avoided for ∆ ∼ 0.1 in which the turning point occurs in the region z ∼ −1 which means a → ∞. When ∆ ≤ 0.02, this point disappears. These results indicate that the turning point in Hubble diagram does not exist in BHDE models under the acceptable conditions.