Dynamical systems methods and statefinder diagnostic of running vacuum (interacting dark energy) models

We study three interacting dark energy models within the framework of four-dimensional General Relativity and a spatially flat Universe. In particular we consider running vacuum models, where dark energy interacts with dark matter, while relativistic matter as well as baryons are treated as non-interacting fluid components. We compute the statefinders parameters versus red-shift as well as the critical points and their nature applying dynamical systems methods. Our main findings indicate that i) significant differences between the models are observed as we increase the strength of the interaction term, and ii) all the models present an unique attractor corresponding to acceleration.

is one of the approaches to quantum gravity, inspired by the well-known Brans-Dicke theory [19,20], where Newton's constant is replaced by an dynamical scalar field following the identification φ → G −1 . In SD cosmological models the cosmological constant becomes time dependent similarly to the running vacuum dynamics, although in the SD scenario Newton's constant, too, acquires a time dependence.
Remarkably, measurements of the expansion rate based on Hubble-diagram of high-redshift objects [50,51] suggest that a rigid Λ term is ruled out by a statistical significance of ∼ 4σ, accounting for deviations from ΛCDM model. The aforementioned deviations allow the possibility of both dynamical and interacting DE, which is realizable within the framework of running vacuum scenario [46]. More generically, interacting DE models are interesting for several reasons. First of all, it is a possibility that should not be ignored, under the assumption that DE and DM do not evolve separately but interact with each other non-gravitationally. For an extensive review on DE and DM interactions, see [52] an references therein. In addition, the "why now problem" may be addressed if our current Universe sits at a stable fixed point (attractor) of the corresponding dynamical system, and this attractor corresponds to acceleration and to 0 < Ω m,0 < 1, with Ω m,0 being today's normalized density of matter. Thus, the system will always reach its attractor at late times irrespectively of the initial conditions. It can be easily shown that this scenario cannot be realized if there is no interaction between DE and matter [53].
As several DE models predict very similar expansion histories, all of them are still in agreement with the available observational data. It thus becomes clear that it is advantageous to introduce and study new appropriate quantities capable of discriminating between different dark energy cosmological models at least at background level. Hence, in order to compare different dark energy models we can introduce parameters in which derivatives of the scale factor beyond the second-order appear. To this end, one option would be to study the so-called statefinder parameters, r, s, defined as follows [54,55] r ≡ ... a aH 3 , where the dot denotes differentiation with respect to the cosmic time t, H =ȧ/a is the Hubble parameter, and q = −ä/(aH 2 ) is the decelerating parameter. We see that the statefinder parameters are expressed in terms of the third derivative of the scale factor with respect to the cosmic time, contrary to the Hubble parameter and the decelerating parameter, which are expressed in terms of the first and the second time derivative of the scale factor, respectively. It is straightforward to verify that for the ΛCDM model without radiation the statefinder parameters take constant values, r = 1, s = 0. These parameters may be computed within a certain model, their values can be extracted from future observations [56,57], and the statefinder diagnostic has been applied to several dark energy models [58][59][60][61][62]. As we will see later on, r, s can be very different from one model to another even if they predict very similar expansion histories. Considering that running vacuum (RVM) models offers an interesting framework to study phenomenology beyond to ΛCDM model, the main goal of the present work is to analyse three scenarios of running vacuum DE models [46] in two respects: On the one hand, by applying the dynamical systems methods, we compute the critical points for each scenario and study their stability. On the other hand, in other to discriminate between the several running vacuum DE models and ΛCDM, we perform the statefinder diagnostic by means computing the statefinder parameters as a function of the redshift, studying their high and low-redshift limits. Our work is organized as follows: after this introduction, we present the formalism in Section 2. In the third section, upon the dynamical system analysis we compute the corresponding critical points, and we also discuss the statefinder parameters. Finally we summarize our findings and present our conclusions in Section 4. We adopt the mostly positive metric signature, (−, +, +, +), and we work in natural units where c = = 1.

II. THEORETICAL FRAMEWORK
We consider a flat (k = 0) FLRW Universe and setting κ 2 = 8πG, with G being the Newton's constant, the scale factor a(t) satisfies the Friedmann equations where ρ A and p A denote the energy density and pressure of each individual fluid component, respectively. The equation-of-state parameter for each fluid component p A = w A ρ A takes the values: w = 0 for baryons and dark matter, w = 1/3 for radiation and w = −1 for DE. The system of cosmological equations also includes the conservation equations for the non-interacting fluids (baryons, as well as for the interacting components (DE and dark matter) Here, Q represents the source term, i.e. the energy exchange between DE and DM. Particularly, in the running vacuum cosmology scenario, the cosmological coupling Λ varies as Λ ≡ Λ(H 2 ) or Λ ≡ Λ(R) [63]. In such a models, the dynamics of vacuum is due to the energy exchange with some of the fluid components that participate to the evolution of the Universe. As running vacuum models (RVM) seem to perform better than the ΛCDM in some circumstances, in the present work we consider three scenarios, namely I, II and III, found e.g. in [46] and precisely labelled as "running vacuum model" (RVM): where the dimensionless parameters {ν i } measure the strength of the interaction term Q i . Following previous works [53,[64][65][66][67] we introduce normalized densities (dimensionless, positive quantities) where ρ cr = 3H 2 /κ 2 is the critical energy density. On the one hand, the first Friedmann equation is a constraint or where Ω m ≡ Ω dm + Ω b . Because of the constraint, there are either two or three independent normalized densities depending on the interacting model. In particular, in scenarios I and III there are two, namely Ω r , Ω Λ , while the third one Ω m = 1 − Ω r − Ω Λ , whereas in scenario II there are three, Ω r , On the other hand, the second Friedmann equation takes the form where we have defined the total equation-of-state parameter w T = p T /ρ T , which is given by Finally, instead of cosmological time t we introduce the number of e-folds N ≡ ln(a), and we define the time derivatives for any quantity A as followsȦ Using the definitions and the cosmological equations one can obtain first order differential equations for Ω A with respect to N . The equations for Ω r , Ω b are the same in all three scenarios since they are non-interacting components The equation for Ω Λ depends on the interaction term Q, and therefore there are 3 cases Finally, q and r are computed to be while s can be computed using its definitions once q and r are known. Thus, we can compute the statefinder parameters, {r, s}, as a function of the red-shift z ≡ a 0 /a − 1 (with a 0 being the present value of the scale factor a), after solving the system of differential equations given by Eqs. (21)- (27) in three different models for the dimensionless densities Ω A . Although a numerical integration of the cosmological equations to obtain {r, s} is possible, in the following we will obtain exact analytical expressions, see next section.

III. ANALYTICAL SOLUTIONS
The system of coupled equations may be directly integrated to obtain concrete expressions for the energy densities in terms of the scale factor, as was done e.g. in Ref. [46]. Although these solutions were previously reported, neither the statefinder diagnostic nor the phase space were analysed. This is precisely the goal of the present article, filling thus a gap in the literature. in this paper we want to complete analysis by including the statefinder diagnostic showing, in figures, how the set {r, s} evolves for different values of redshift, as well as the phase space of the above parameters. We start by considering the corresponding dark matter density ρ dm and dark energy density ρ Λ respect to the scale factor for each model, i.e.: Thus, for each particular model the above profile densities give the evolution of dark matter and dark energy respectively. It is important to point out that the models analysed here boil down to the ΛCDM model when ν i → 0. Finally, for convenience, we introduce the dimensionless Hubble rate E(z) ≡ H(z)/H 0 , where H 0 = 100h (km sec −1 )/Mpc) is the Hubble constant. Accordingly, the parameters {q, r, s} are computed as follows and s(z) is given by (2), where X z ≡ dX/dz for any quantity X. Using the expressions for the energy densities shown before, one can obtain exact analytical expressions for all quantities of interest versus red-shift, E(z), q(z), r(z), s(z), see section V.

Fixed point (Ωr, ΩΛ)
Eigenvalues  We briefly review the stability analysis based on the nature of the fixed points (FPs), see e.g. [53,[64][65][66][67]. Suppose that for a dynamical system with a two-dimensional phase space (x, y), its time evolution is determined by the following system of coupled first order differential equations First, the fixed point(s) is (are) computed setting dx/dt = 0 = dy/dt, and one has to solve the system of two algebraic equations F (x 0 , y 0 ) = 0 = G(x 0 , y 0 ). Then, to determine the nature of the fixed point(s) we linearise the equations around that point, x(t) = x 0 + δx, y(t) = y 0 + δy ignoring higher order terms. One obtains a system of two coupled linear equations of the formẊ where the column X contains the two functions δx(t), δy(t), while A is a two-dimensional matrix, the elements of which are given by Finally, we compute the eigenvalues λ 1 , λ 2 of A, the sign of which determines the nature of the fixed point(s). In particular, the critical point is stable (A) when both eigenvalues are negative, unstable (R) when both eigenvalues are positive, and a saddle point (S) if the eigenvalues are of opposite sign. Furthermore, if q(x 0 , y 0 ) < 0, w T (x 0 , y 0 ) < −1/3, the fixed point at hand corresponds to acceleration. The procedure may be easily generalized in a straightforward manner for a three-dimensional phase-space.
The fixed points and their nature (stability conditions) for all three models considered in this work are shown in the Tables I,II,III,IV,V, and VI.   Tables I and II. The critical point I.a is a dark-energy-dominated de Sitter solution with Ω Λ = 1 and w DE = w T = −1. By analysing the stability through the matrix of linear perturbations we get for this fixed point the eigenvalues So, it is a stable FP (attractor) for ν < 1 and a saddle point for ν > 1. Point I.b is a scaling solution with Ω Λ = ν, w DE = −1, and the total equation of state is given by w T = −ν. Since 0 < Ω Λ < 1, thus one has that 0 < ν < 1. Furthermore, for this solution the accelerated expansion occurs only when ν > 1/3, otherwise we have decelerated expansion. For point I.b we obtain the eigenvalues Thus, it is always a saddle point because µ 1 < 0 and µ 2 > 0 in the physical range of ν. It is worth noticing here that for model I this solution represents the matter solution during the dark-matter-dominated epoch once that we restrict the value of ν to be ν ≈ 0. For slightly higher values of ν we may have important consequences for structure formation as the matter fluctuations can be suppressed with respect to the ΛCDM model [68]. Finally, point I.c is also a scaling solution for the which Ω r = 1 + 3ν, Ω m = −4ν, and Ω Λ = ν, with total equation of state w T = 1/3. This point I.c is not physically acceptable since it does not satisfy the condition 0 < Ω A < 1.

Fixed point (Ωr, ΩΛ)
Eigenvalues In the case of model II we obtain four critical points which are shown in Tables III and (IV). Point II.a is a matter dominated solution representing ordinary baryonic matter, such that Ω b = 1, and with w T = 0. The eigenvalues for this critical point are and therefore it is always a saddle point. Point II.b corresponds to a radiation dominated solution, Ω r = 1, for which one has that w T = 1/3 and therefore there is not acceleration. For this fixed point we find the eigenvalues Fixed point Existence Acceleration Clearly, for ν dm < 1, point II.c is a stable node and therefore an attractor. The last solution for this model is the fixed point II.d which is a scaling solution with Ω Λ = ν dm , as the physical requirement implies 0 < ν dm < 1. Also, this solution is characterized by w T = −ν dm , with the decelerating and accelerating regimes satisfying 0 < ν dm < 1/3 and ν dm > 1/3, respectively. Similarly to the solution I.b of model I, point II.d behaves as a dark matter solution in the limit ν dm 1, with a small contribution of dark energy proportional to ν dm , during the matter dominated epoch, and thus suppressing the growth of matter perturbations. Stability analysis leads us to the eigenvalues Since we require 0 < ν dm < 1, the point II.d is always a saddle point. So, the thermal history of the Universe is successfully reproduced for model II provided it satisfies the restriction ν dm 1.

C. Model III
Model III, has three critical points which are shown in Tables V and VI. Point III.a is a dark matter dominated solution for the which Ω m = 1, and w T = 0. The eigenvalues associated with this critical point are and hence, one can see that it is always a saddle point for the value 0 < ν Λ < 1.
On the other hand, point III.b is a solution for which radiation component is the dominant one, being that Ω r = 1 and w T = 1/3. From stability analysis we obtain the eigenvalues which means that it is an unstable node for ν Λ < 4/3 and a saddle in the opposite case ν Λ > 4/3. Finally, the point III.c is a scaling solution with Ω Λ = 1 − ν Λ that, interestingly, can also be an attractor with w T = −1 + ν Λ , allowing to alleviate the so-called cosmological coincidence problem [28]. The physical requirement 0 < Ω Λ < 1 implies 0 < ν Λ < 1, and from the constraint w T < −1/3, the accelerated expansion occurs for the value ν Λ < 2/3. For this critical point one finds the eigenvalues So, point III.c is a stable FP for ν Λ < 4/3 and a saddle point for 4/3 < ν Λ < 1. Like model II, the present model III is also physically viable to successfully reproduce the thermal history of the Universe from the radiation dominated era, going through the standard matter dominated era, to late times when the dark energy component dominates the total energy density and pressure of the Universe.

V. STATEFINDER ANALYSIS
In the present analysis we have calculated for all three models I, II, and III, the analytical expression of the Hubble parameter E(z) as explicit function of the redshift z, and then we have obtained the corresponding functions q(z), r(z) and s(z). It is important to observe that in the s−r plane, the flat ΛCDM scenario correspond to the point (0, 1), while that in the q − r plane the point (−1, 1) is the asymptotic de Sitter solution [69]. Also, as it has been highlighted in the previous dynamical analysis, Model I has a serious drawback in relation with the other models because it fails in describing the standard radiation dominated era, since if we impose the Friedmann equation constraint on the energy densities, it leads to a negative fractional matter energy density for ν > 0. Since we are only focused on the evolution at late times, particularly in the transition from matter dominated era to the present time, we neglect the radiation component in the computation of E(z) and {q, r, s}. However, in order to produce the plots shown in Fig.(1) we take into account the contribution coming from radiation.
For Model I, we find the Hubble parameter as Therefore, by using Eqs. (28), (29) and (2) we obtain At this level, it is interesting to analyse the limiting cases, namely: i) when z → −1 and ii) when z 1. For at high-refshifts, z 1, Eq. (46) becomes Accordingly, at high-refshifts, the scale factor behaves as a(t) ∼ a 0 (t/t 0 ) α . So, from Eqs. (47), (48) and (49), the asymptotic values for q, r and s, respectively, when z 1, are given by This solution represents the power-law expansion of a non-relativistic matter-dominated universe which shows a small deviation from the standard matter-dominated era (α = 2/3) for ν 1.
On the other hand, in the limit z → −1, the Universe reaches the late-time de-Sitter attractor such that with q = −1, r = 1 and s = 0. So, at the present time z = 0, the values of the state-finder parameters are and thus, for the state-finder parameters we have Now, the Hubble rate at z 1 yields while the parameters q, r and s behaves as Therefore, in this case we recover the standard matter-dominated era with 0 < s < 1.
In the limit z → −1, the model predicts a de Sitter solution with such that q = −1, r = 1 and s = 0. Hence, at the present, z = 0, we obtain the values and we also we obtain the expressions Similarly as in model II, for z 1, it is recovering the standard matter-dominated era, where the Hubble rate satisfies the relation (62). Also, in this limit we find with 0 < s < 1.
On the other hand, in the limit z → −1, unlike models I and II, the model III behaves as a scaling solution with with a scale factor of the form a ∼ a 0 (t/t 0 ) β . It is straightforward to check that in this limit, the statefinder parameters become In Fig. 1 we have depicted the behaviour of the parameters q, r and s as functions of the red-shift z, along with the trajectories of evolution in the q − r and s − r planes, for Models I (blue dashed line), II (red-dotted), and III (orange dash-dotted). For a sake of comparison, we have also included the ΛCDM model (solid black line). Recall that, in order to produce the plots shown in Fig. 1, the contribution from radiation has been taken into account, and we also set the following values for the fractional energy densities: Ω 0 dm = 0.26, Ω 0 b = 0.04, and Ω 0 r = 9 × 10 −5 . Naturally, Ω 0 m ≡ Ω 0 dm + Ω 0 b and Ω 0 Λ = 1 − Ω 0 m − Ω 0 r . In these plots the behaviour of the parameters q(z), r(z) and s(z) is in agreement with the analytical results that we have obtained in the limit case of a negligible radiation component for z → −1. For all the three models I, II and III, it is seen that the pair (s, r) starts in the left-hand side of the ΛCDM fixed point, which is characteristic of the hybrid expansion law (HEL), Chaplygin gas and Galileon models, such that s < 0 and r > 1 [70]. Let us notice that this behaviour is very different from what occurs in the case of the quintessence model for which it is observed that the trajectory in the (s, r) plane starts in the region 0 < s < 1 and r < 1. On the other hand, the trajectory in the (q, r) plane starts in the region bounded by 0 < q < 1 and r > 1, being that in the case of models I, and II, the ΛCDM line is crossed at some red-shift in the past to then evolve towards the de Sitter fixed point (q = −1, r = 1) at the future. For model III, the behaviour becomes different because in this case the asymptotic fixed point is a scaling solution such that for ν Λ = 0.3, the benchmark values of the fractional energy densities are Ω Λ = 0.7 and Ω m = 0.3, as it has been depicted in Fig. 2.
Note added: One day before we received the referee report, a new work related to ours appeared [71], which indicates that the topic is interesting. Our work was carried out in complete independence from them, and vice versa. In that work, the authors have analysed the dynamics and evolution of several Λ-varying cosmological models. The critical points and their nature have been determined, and the corresponding phase space is shown, although they have not discussed at all the statefinder parameters. We have checked that where there is overlap their results are in agreement with ours, which is a confirmation that both calculations are error-free.

VI. CONCLUDING REMARKS
In summary, in the present work we have applied phase-space dynamical techniques to three running vacuum dark energy models, and we have computed the statefinder parameters as functions of the red-shift.
From our dynamical analysis we have shown that the three models, I, II, and III, can explain the current accelerated expansion phase of the Universe, being that the corresponding fixed point, either a de Sitter solution (Model I and II) or scaling solution (Model III), is also an attractor in all the cases, provided that ν i = {ν, ν dm , ν Λ } < 1. Also, for Model II and III the thermal history of the universe can be successfully reproduced, from the radiation-dominated era, passing through the matter-dominated era, and finally reaching the dark energy-dominated phase. In this aspect, Model I has a serious drawback in relation with the other models because it does not satisfy the physical requirement of a positive fractional energy density such that 0 < Ω A < 1.
More interestingly, in the case of Model II, the fixed point representing the matter era is a scaling solution in which the dark energy density has a contribution to the total energy density during the dark matter era, with Ω Λ = ν dm . So, when the parameter ν dm is small this contribution is also small. Thus, the parameter ν dm indicates a slight deviation with respect to the standard mater era whose value may be more constrained from large-scale structure (LSS) data. It is due to the fact that when the universe enters in this fixed point the growing of matter density perturbations can be suppressed by the presence of dark energy, such that the dark matter density contrast grows less rapidly that the first power of the scale factor, depending on the amount of dark energy [68]. On the other hand, in the case of model III, the final attractor is a scaling solution with accelerated expansion which satisfies Ω Λ = 1 − ν Λ , w DE = −1 and w T = −1 + ν Λ . This class of solution is very interesting because it provides a natural mechanism in alleviating the fine-tuning problem, or cosmological coincidence problem of dark energy [28]. It can adjust the current values of the cosmological parameters such as Ω 0 Λ = 0.7 and Ω 0 m = 0.3, and at the same time explain the accelerated expansion. Regarding the statefinder analysis, first let us recall that the pair {r, s} is defined using third order time derivatives of the scale factor, and that the statefinder parameters have the potential to discriminate between different dark energy models. Now, we should differentiate between the three cases shown in Fig. 1. The first column was plotted for ν i = 0.3, the second column was depicted for ν i =0.07, whereas in the third column it is assumed that ν i = 0.007. Starting from the left panel, we observe evident discrepancies between the models. This difference is also natural because running vacuum models are quantum-inspired, which means that the deviations with respect to the classical counterpart should be, in general, small. Taking the above idea seriously, the parameter ν i , which encodes the quantum features, should be taken in such a way that the effect on the classical solution will be soft. We then claim that the last column in Fig. 1 should be taken as a more suitable situation. Furthermore, notice that when ν i is taken to be close to zero (third column), the deceleration parameter q(z) looks qualitative identical to those for ΛCDM. If we now analyze the parameter r, we observe once more a notorious similarity to the standard scenario, although the second parameter (i.e., s) exhibits a remarkable difference. Thus, although classically these models should be equivalent, at the level of the statefinder diagnostic, this is not the case.