Dynamical systems methods and statender diagnostic of interacting vacuum 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 first consider two vacuum models where dark energy interacts with dark matter, while relativistic matter as well as baryons are treated as non-interacting fluid components. Secondly, we investigate a third model where the gravitational coupling is assumed to be a slowly-varying function of the Hubble rate and dark energy and dark matter interact as well. We compute the statefinders parameters versus red-shift as well as the critical points and their nature applying dynamical systems methods. In the case of only an interaction term, 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. On the other hand, when we allow for a variable gravitational coupling, we find that (i) the deviation from the concordance model depends of both the strength of gravitational coupling parameter and the interaction term, and (ii) there is an unique attractor corresponding to acceleration.


Introduction
The origin and nature of dark energy (DE), the fluid component that currently accelerates the Universe [1][2][3], is one of the biggest mysteries and challenges in modern theoretical Cosmology. Clearly, Einstein's General Relativity [4] with radiation and matter only cannot lead to accelerating solutions. A positive cosmological constant [5] is the simplest, a e-mail: grigorios.panotopoulos@tecnico.ulisboa.pt (corresponding author) b e-mail: angel.rincon@pucv.cl c e-mail: giovanni.otalora@pucv.cl d e-mail: nelson.videla@pucv.cl most economical model in a very good agreement with a great deal of current observational data. Since, however, it suffers from the cosmological constant (CC) problem [6], other possibilities have been considered in the literature over the years. The CC problem, introduced by Zeldovich for the first time more than fifty years ago [7], may be summarized in a few words as follows: It is an impressive mismatch-by many orders of magnitude-between the observational value of vacuum energy, and the expected value from particle physics due to vacuum fluctuations of massive fields. Although some progress has been made up to now, see e.g. [8][9][10][11], the origin of the CC problem still remains a mystery.
Regarding the CC problem and possible alternatives to the CDM model, either a modified theory of gravity is assumed, providing correction terms to GR at cosmological scales, or a new dynamical degree of freedom with an equation-of-state (EOS) parameter w < −1/3 must be introduced. In the first class of models (geometrical DE) one finds for instance f (R) theories of gravity [12][13][14][15], brane-world models [16][17][18] and Scalar-Tensor theories of gravity [19][20][21][22], while in the second class (dynamical DE) one finds models such as quintessence [23], phantom [24], quintom [25], tachyonic [26] or k-essence [27]. For an excellent review on the dynamics of dark energy see e.g. [28].
Furthermore, regarding the value of the Hubble constant H 0 , there is nowadays a tension between high red-shift CMB data and low red-shift data, see e.g. [29][30][31][32]. The value of the Hubble constant extracted by the PLANCK Collaboration [33,34], H 0 = (67 − 68) km/(Mpc s), is found to be lower than the value obtained by local measurements, H 0 = (73 − 74) km/(Mpc s) [35,36]. This tension might call for new physics [37]. What is more, regarding large scale structure formation data, the growth rate from red-shift space distortion measurements has been found to be lower than expected from PLANCK [38,39].
Both tensions may be alleviated within the framework of running vacuum dynamics [40][41][42][43][44][45][46]. In this class of models, contrary to a rigid cosmological constant = const., vacuum energy density can be expressed as a function of the Hubble rate, i.e. ρ = ρ (H ), being of dynamical nature and at the same time it may interacts with dark matter, and also it accounts for a running of the gravitational coupling G [47]. We remark in passing that other alternatives approaches to running vacuum dynamics do exist, and one may mention for instance the scale-dependent (SD) scenario [48][49][50], in which it is assumed that the couplings of the original classical action acquire a scale-dependence. The SD scenario 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 [51,52] 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. Secondly, and perhaps the main motivation for an interaction in the dark sector, is currently motivated that this scenario can solve the current cosmological tensions in some data, see e.g. [53][54][55]. Additionally, other recent relevant results regarding the interaction between DE and DM were found in [56,57]. For an extensive review on DE and DM interactions, see [58] 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 [59].
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 [60,61] where the dot denotes differentiation with respect to the cosmic time t, H =ȧ/a is the Hubble parameter, and q = −ä/(a H 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 [62,63], and the statefinder diagnostic has been applied to several dark energy models [64][65][66][67][68]. 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 models within the running vacuum dynamics: we first consider two vacuum models where dark energy interacts with dark matter [46], and secondly, we investigate a third model where the gravitational coupling is assumed to be a slowly-varying function of the Hubble rate [47] and dark energy and dark matter interact as well. The analysis is performed 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 basic equations and analytical solutions for Models I and II in Sects. 2 and 3, respectively. In the fourth section, upon the dynamical system analysis we compute the corresponding critical points for Models I and II, while in the fifth section we discuss the statefinder parameters for the same models. In Sect. 6, we present the main results for Model III regarding the dynamical system and statefinder analysis. Finally we summarize our findings and present our conclusions in Sect. 7. We adopt the mostly positive metric signature, (−, +, +, +), and we work in natural units where c =h = 1.

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 equationof-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, radiation) as well as for the interacting components (DE and dark mat- 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) [69]. 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 first consider two scenarios, namely I and II, 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 . In the present paper we are interested in studying the late-times cosmology within the interacting vacuum energy scenarios. So, we have neglected the coupling to radiation and baryons because at lower redshifts, z 2, the contribution to the total energy density coming from these components is smaller than the dark energy and dark matter components. It is worth to mention that an eventual coupling between radiation and dark energy could have a significant effect on the dynamics of early Universe, see e.g. [28]. Particularly, within the interacting vacuum energy, the nucleosynthesis sets strong constraints on the strength of the coupling, which becomes much smaller than the unity [70].
Following previous works [59,[71][72][73][74] 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 scenario I there are three, r , , b , in contrast, in the scenario II there are two, namely r , , with the third one being m = 1 − r − , while the fourth dm = 1 − r − − b . 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 followṡ 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 the scenarios since they are non-interacting components The equation for depends on the interaction term Q, and therefore there are three 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. (20)- (26) in two 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.

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.: for Model II (26) 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 where H 0 = 100 hkm/(Mpc s) is the Hubble constant. Accordingly, the parameters {q, r, s} are computed as follows and s(z) is given by (2), where X z ≡ d X/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 Sect. 5.

Dynamical systems methods
We briefly review the stability analysis based on the nature of the fixed points (FPs), see e.g. [59,[71][72][73][74]. 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 forṁ 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 pro- cedure may be easily generalized in a straightforward manner for a three-dimensional phase-space.
The fixed points and their nature (stability conditions) for all two models considered in this work are shown in the Tables 1, 2, 3, and 4.

Model I
In this case we obtain four critical points, which are shown in Tables 1 and 2. Point I.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. This fixed point is not physical because it represents an era dominated by baryons. It is well known that cold dark matter constitutes the dominant component during the matter-dominated era at redshift 1 z 10 3 , and thus during this epoch it provides the main contribution for structure formation in the universe. Point I.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 which means that it is always an unstable FP for ν dm > 0.
On the other hand, point I.c is a de Sitter-dominated solution for which = 1 and w DE = w T = −1. So, this solution presents accelerated expansion for all values of ν dm . In this case, we find the eigenvalues Clearly, for ν dm < 1, point I.c is a stable node and therefore an attractor. The last solution for this model is the fixed point I.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. Point I.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 I.d is always a saddle point. The critical point I.d is a dark matter dominated solution for the model I, with a small contribution from dark energy density given by = ν dm 1 and total equation of state w T = −ν dm ≈ 0. Although this fixed point can provide accelerated expansion for ν dm > 1/3, the observational constraints on ν dm 1 [70] do not allow that this happens. So, the thermal history of the Universe is successfully reproduced for model I provided it satisfies the restriction ν dm 1. It is also important to note that in the present framework of dynamical systems the negative values for ν dm are excluded since this would imply a negative energy density = ν dm < 0. So, as several authors usually do, we have given preference to maintain the physical condition ρ ≥ 0 in agreement with the weak energy condition (WEC) [28].

Model II
Model II, has three critical points which are shown in Tables 3 and 4. Point II.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 II.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 II.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 II.c is a stable FP for ν < 4/3 and a saddle point for 4/3 < ν < 1. Like model I, the present model II 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. Considering the observational constraints as, for instance, those found in [75,76] (and references therein) the value of the strength of the coupling is quite smaller that unity, such that ν 0.01, and this is in agreement with our theoretical bounds. So, the results obtained from the dynamical analysis in this section should be supplemented by the observational bounds. In fact, although the scaling solution for model II allows us to adjust 0 0.73 and 0 m = 0.27 for ν = 0.3, actually this solution cannot reached at z = 0, but only asymptotically to reproduce the whole thermal history of the universe. In other words, in order to obtain the physical trajectory in the phase space consistently with observational data, this scaling solution only can be reached at the future, in such a way that we need to choice smaller values of ν ( 10 −2 ), allowing to have 0 0.73 at z = 0, and asymptotically

Statefinder analysis
In the present analysis we have calculated for both models I, and II, the analytical expression of the Hubble parameter E(z) as explicit function of the red-shift 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 [77]. 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 and thus, for the state-finder parameters we have r (z) = 1 + 9(1 − ν dm )ν dm 2 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 matterdominated 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 I, for z 1, it is recovering the standard matter-dominated era, where the Hubble rate satisfies the relation (47). Also, in this limit we find with 0 < s < 1.
On the other hand, in the limit z → −1, unlike models I, the model II 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 Thus, for ν = {0.01, 0.001} one gets q 0 ≈ {−0.985, −0.999}, r 0 ≈ {0.9555, 0.9955}, and s 0 = {0.01, 0.001}.
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 (red-dotted), and II (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, and II, it is seen ) for Model II and for ν = 0.01. Different trajectories correspond to different initial conditions. All of them meet at the attractor ∼ (0.0, 1.0) at late times labelled with a point. Notice that our solution does not reach the de-Sitter solution, but it is undistinguishable from it 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 [78]. 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, 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 II, the behaviour becomes different because in this case the asymptotic fixed point is a scaling solution such that for ν = 0.01, the benchmark values of the fractional energy densities are 1 and m 0, but at z = 0, we have 0.7 and m 0.3, as it has been depicted in Fig. 2.

Variable gravitational coupling
In this section we generalize our previous results by allowing a variable gravitational constant within the framework of interacting RVM's.
The generalized Friedmann equations are given by where G(t) is a function of the cosmic time t, and ρ m denotes the non-relativistic matter energy density, including both baryons and dark matter. The conservation law for this model can be written aṡ Following our analysis for interacting RVM's, this equation can also be splitted into a set of two separate evolution equations for the energy densities ρ m and ρ , according to Ref. [47], as followṡ Let us note that the two above equations are reduced to the Eqs. (8), (9), in the case when the gravitational coupling becomes constant. On the other hand, in order to compare this model with current observational data regarding dark energy, we rewrite the above Friedmann equations in the standard form with G 0 being is the constant gravitational coupling. In these equations we have introduced the effective energy density and pressure density of dark energy, defined as respectively. This effective dark energy density includes the effect of both vacuum energy density ρ (t) and the dynamically changing gravitational coupling G(t). It is straightforward to show that it satisfies the evolution equatioṅ Therefore, it is easy to see that the two evolution equations (70) and (75) together are consistent with the energy conservation law for ρ m and ρ DE . As it is usually done, we introduce the fractional and critical energy densities the EOS parameter of dark energy and the total EOS parameter With this, the condition for having accelerated expansion becomes w e f f < −1/3, or equivalently q < 0.

An example for variable G
In order to obtain concrete results, we assume the following phenomenological ansatz for the gravitational coupling depending on the Hubble rate, according to [47] G(X ) where X ≡ E 2 = (H/H 0 ) 2 . Also, in order to extend our previous analysis for interacting RVM's, we consider the coupling between DE and DM to be Q = 3ν m Hρ m (Model II). For this coupling function Q and ansatz (79), the equations (70) and (75) take the form whose solutions are given by and respectively. Upon replacement of the above solutions into Eq. (71), one obtains ρ 0 cr = ρ 0 m + ρ 0 de . On the other hand, doing the same but now in Eq. (72) we find After solving this equation, we obtain the dimensionless Hubble rate squared E 2 (N ), and then H (N ), or equivalently, H (z), by introducing N = − ln(1 + z). Also, from (76), one has that and m (X ) = 1 − de (X ). From Eqs. (77), (78), we find The Eq. (84) is an autonomous equation which can be treated following the same framework of dynamical systems. This equation has a single critical point X c which can be determined from For example, numerically we have found that for 0 de 0.73, ν dm 0.01 and ν G 1 × 10 −3 , the value X c is approximately X c 0.73. This result gives us the asymptotic (future) value of the Hubble rate by using the relation H c = H 0 √ X c . So, from the observational data of Planck 2018 Ref. [34], one has that H 0 = 67.4 km/(Mpc s) and therefore H 0 57.48 km/(Mpc s) at z → −1.
The stability of this fixed point can be studied by considering the solution where the perturbation δ X (N ) satisfies δ X (N ) 1. Thus, replacing it in Eq. (84) we obtain whose solution is with μ = −3(1 − ν m ). So, since ν m 1 then μ < 0 and accordingly, the fixed point is always an attractor. We note that the stability does not depend on ν G at leading order in perturbation.
Similarly, by using Eqs. (28) and (2), the statefinder parameters r and s are computed to be which reduce to those already obtained for Model II in the limit ν G → 0. In Fig. 3 (left panel), it is shown the evolution of the Hubble rate H (z) for the present model by solving the differential equation (84) for some values of ν dm and ν G . It is also added the Hubble rate H C DM of the CDM model, along with the current available data for H (z) from [79] and [80]. Also, we depict the behaviour of the exact relative difference with respect to the concordance model, for fixed ν dm and a pair of different values of ν G . We take ν dm = 0.01 and two different values of ν G , the first value ν G = 5 × 10 −4 and the second one ν G = 1 × 10 −3 , which are included into the physical range ν G ∈ 5 × 10 −4 , 1 × 10 −3 , obtained from observations, see for example Refs. [40,42,81]. It is observed an increasing of r H for higher red-shifts, z 2, and after the present time z = 0, in the future. Particularly, for z 10 we obtain r H 3%, whereas that for z = −1 we have r H 0.18%. In Fig. 4 we plot the variation with respect to z of some cosmological parameters such the fractional energy densities of Also, when z → −1, the model tends asymptotically to an attractor which is a de Sitter solution with de = 1, dm = 0, w de = −1 and w T = −1. For higher red-shifts we observe that w de becomes more sensitive to the values of ν G , taking smaller values than −1 and going deeper in the phantom regime for larger values of ν G . Nevertheless, let us note that the energy density of dark energy decays very quickly and the effective cosmic fluid behaves as nonrelativistic matter with w T 0, and therefore allowing the existence of the standard matter-dominated era [28,82].
Finally, in Fig. 5 we show the evolution of the statefinder parameters r (z) (left panel) and s(z) (right panel) as functions of red-shift, for the same set of values of parameters used in above plots. It can be seen that at z = 0, and for ν dm = 0.01, ν G = 5 × 10 −4 (short dashed line), these parameters take the values r 0.988 and s 3.75 × 10 −3 . For the larger value ν G = 1×10 −3 (large dashed line), at z = 0, we get r 0.988 and s 3.80 × 10 −3 . So, for lower red-shift, there is a small difference (of the order of 4% for z 2) between the results when varying ν G , and this difference is much smaller for the values of r than for s. When z → −1, the trajectory of the system in the plane of statefinder parameters tends toward the de Sitter expansion at the future, with r = 1 and s = 0.

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 two models, I, and II can explain the current accelerated expansion phase of the Universe, being that the corresponding fixed point, either a de Sitter solution (Model I) or scaling solution (Model II), is also an attractor in all the cases, provided that ν i = {ν dm , ν } < 1. Also, for Model I and II 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 energydominated phase.
More interestingly, in the case of Model I, 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 [83]. On the other hand, in the case of model II, 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 discriminate between the two cases shown in Fig. 1. The first column was plotted for ν i = 0.01, the second column was depicted for ν i = 0.001. 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 (second column), the deceleration parameter q(z) looks qualitative identical to those for CDM. If we now analyse 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.
Finally, it has been performed a further analysis regarding to a class of models with variable gravitational coupling within the interacting vacuum scenario, by using the tools dynamical system and statefinder analysis. In doing so, we have studied a particular model for G(H ) from Ref. [47], and Q = 3ν dm Hρ dm . It has been introduced an effective cosmic fluid for describing DE by defining, both effective energy density and pressure, which contain the contribution of the running of gravitational coupling. We have shown that the model has only fixed point which is an attractor and de Sitter solution, allowing to adjust the current data of H (z), along with the other cosmological parameters such that the fractional energy density of DE and the equation of state of dark energy at the present time. In particular, we found that the effect of ν dm becomes more significant for higher redshift, and at the future when z → −1 in comparison with the CDM model. Regarding the statefinder analysis we observe that the s parameter becomes more sensitive to higher redshifts than the r parameter, when we vary the couplingν G . As a final remark, an exhaustive study including several anszats for both the gravitational coupling and the interaction between DE and DM deserves a separated project, reason why these ideas will de addressed in a future work. Note added: One day before we received the referee report, a new work related to ours appeared [84], 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. tionally, N. V. would like to express his gratitude to the Instituto Superior T écnico of Universidade de Lisboa for its kind hospitality during the final stages of this work. The author G. O acknowldeges DI-VRIEA for financial support through Proyecto Postdoctorado 2019 VRIEA-PUCV.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The present work is a theoretical study, and therefore no experimental data has been listed.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .