Cosmological singularities and analytical solutions in varying vacuum cosmologies

We investigate the dynamical features of a large family of running vacuum cosmologies for which Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document} evolves as a polynomial in the Hubble parameter. Specifically, using the critical point analysis we study the existence and the stability of singular solutions which describe de-Sitter, radiation and matter dominated eras. We find several classes of Λ(H)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (H)$$\end{document} cosmologies for which new analytical solutions are given in terms of Laurent expansions. Finally, we show that the Milne universe and the Rh=ct\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{h}=ct$$\end{document} model can be seen as perturbations around a specific Λ(H)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (H)$$\end{document} model, but this model is unstable.


Introduction
Over the last two decades, most studies in cosmology strongly indicate that the universe is spatially flat and it consists of ∼ 4% baryonic matter, ∼ 26% dark matter and ∼ 70% of dark energy (hereafter DE), thought to be driving the phenomenon of cosmic acceleration [1][2][3][4][5][6]. Although there is a general consensus regarding the main properties of DE, namely it has a negative pressure, the origin of this unexpected component of the universe has yet to be identified. This has given rise to a plethora of alternative cosmological scenarios which mainly generalize the nominal Einstein-Hilbert action of general relativity using either a new field in nature [7][8][9][10], or a non-standard gravity theory that increases the number of degrees of freedom [11][12][13][14][15][16]. The introduction of a cosmological constant, , is probably the simplest modification of the Einstein-Hilbert action which can be considered. In the framework of the so called CDM model, the cosmological constant coexists with cold dark matter (CDM) and ordinary baryonic matter (see [17] for review). Although the CDM model fits accurately, the current cosmological data suffers from two problems [18][19][20][21]. The first is the 'old' cosmological problem, namely the expected (Planck natural unit) vacuum energy density is ∼ 120 orders of magnitude larger than the presently observed value of . The second problem is the coincidence problem: that the density of dark energy is so similar to the matter density today (the two were equal when the universe had expanded to about 75% of it present expansion scale).
There was a great effort to explore the (H ) models both analytically and observationally but a dynamical analysis based on critical points is still missing. The purpose of our work is to bridge this gap, by determining the de Sitter phases of a general family of (H ) models to search for singular solutions of the form a (t) ∝ t p which secure the existence of radiation ( p = 1/3) and matter ( p = 2/3) dominated eras, respectively. We will investigate the stability of the critical points in order to understand the dynamical and cosmological properties of the (H ) models. Here, the main mathematical tool that we use is that of the singularity analysis of differential equations, and more specifically we apply the ARS (Ablowitz, Ramani and Segur) algorithm [57][58][59]. This algorithm provides a method to construct the analytical solution of a differential equation which is expressed as a Laurent expansion around the singular leading-order term (for some applications on gravitational theories see [60][61][62][63][64][65] and references therein). Information regarding the stability of the trajectories close to the singularity can be extracted directly from the ARS algorithm.
The structure of the manuscript is as follows. In Sect. 2 we briefly introduce the concept of the running (H ) cosmologies. In Sects. 3 and 4 we present the main results of our work, namely we study the critical points and their stability as well as we provide the corresponding analytical solutions. Finally, in Sect. 5 we draw our conclusions.

-Varying cosmology
In this section we briefly present the main points of the running vacuum cosmology. If we model the expanding universe as a perfect fluid with density ρ, and corresponding pressure 1 p = wρ, then the energy-momentum tensor is given by T μν = −p g μν + (ρ + p)U μ U ν . In this context, the term g μν can be absorbed by the total energy momentum tensor T μν ≡ T μν + g μν ρ , where ρ = /(8π G) is the vacuum energy density which is related to , while the corresponding equation of state is p = −ρ . Combining the above expressions we arrive at and thus the Einstein's field equations become For the spatially flat FLRW metric the field equations boil down to those of Friedmann equations, namely where for simplicity we have set 8π G ≡ 1. Although, a constant term is the simplest possibility, it is interesting to mention that the Cosmological Principle embodied in the FLRW spacetime allows to be a function of the cosmic time, or of any collection of homogeneous and isotropic dynamical variables, i.e.
In order to find more general solutions for this scenario, ones in which there is evolution of the form of the expansion rate over time, we need to explore a more general functional form for (H ). The case of viable running vacuum can be placed in the general framework of quantum field theory in curved spacetime [50,66,67]. Specifically, in ref. [26] the following expression (for recent review see [68]) was proposed for the functional form of (H ), where k > 0. Notice, that this formula is well motivated in the general context of the re-normalization quantum field theory (QFT) in curved spacetime and thus only even powers in H are allowed by the general covariance of the theory (see [50], [51][52][53]). In this family of running vacuum models the spacetime emerges from a pure non-singular de Sitter vacuum stage, allowing a graceful exit from inflation to a radiation-dominated phase, followed by dark matter and vacuum regimes, followed by evolution to a late-time de Sitter phase [28,29,68]. From the observational viewpoint, the current vacuum model agrees with the latest expansion data and it predicts a growth rate of clustering which is in agreement with the observations (for more details see [55,68,69] and references therein). Inspired by Ref. [26] and based only on phenomenological arguments we introduce in our analysis the following form of the running vacuum where the leading term c 0 describes in good approximation the current universe and the other terms introduce a dynamical evolution which can be important under special conditions. Unlike Eq. (6), here the exponent n is allowed to take positive and negative values as well as it can be either even or odd. From the phenomenological point of view, the parametrization (7) describes a large family of dynamical models of the vacuum energy (see also [68]). As expected, for the special case where n = −2k the current form of (H ) reduces to that of Eq. (6).
Furthermore, for a general quadratic vacuum model, with where F(H ) is an arbitrary function. Introducing the following parametrizations, Equation (8) then becomes This equation is just the original master equation (8) with a different value for w, which means that the quadratic term H 2 of Eq. (7) can be easily absorbed into the Friedmann equation (8). In this context, if we decide to use Eq. (10) in order to study the current class of decaying vacuum models, then (H ) takes the form where the new constants λ 0 and λ 2 are related with those of (7) by c 0 = 3−c 1 3 λ 0 and c 2 = 3−c 1 3 λ 2 . At this point, it is important to mention that the equation of state parameter for the matter source has been changed as given by (9).
In order to understand the possible evolutionary tracks of the cosmic fluid, we consider the following two cases: Case A: Assume that the matter source is that of an ideal fluid for which the equation of state parameter w 1 is constant. Hence, from Eq.(10), we deduce that (12) or, with the aid of H =ȧ/a, we find At this point we would like to discuss the following issue: for the vacuum model with n > 0 is it possible to have a super-accelerated expansion of the universe in the far future?
In order to answer this question we need to compute the deceleration parameter q = −(1 +Ḣ H 2 ) by numerically integrating Eq. (13) for non-relativistic matter w 1 = 0. Concerning the initial conditions we use the present values, mainely, a (t 0 ) = 1 andȧ (t 0 ) = 70, λ 0 = 3 m0 H 2 0 − λ 2 H −n 0 , with m0 0.3. We find that in the far future the deceleration parameter tends to that of de-Sitter cosmology (q ≈ − 1) and thus we exclude the possibility of a super-accelerated expansion of the universe. In Fig. 1 we plot the evolution of the deceleration parameter for various λ 2 and n values. Practically, H −n with n > 0 does not really affect the cosmic expansion at late times. Moreover, one can easily check that in this case the term H −n plays no role in the early inflation. Under the above conditions in order to provide a viable running vacuum model one may introduce other additional terms, namelyḢ n in the form of (H ) (see [68] and references therein). Case B: Here we assume that the cosmic fluid contains two ideal fluids (nominally case matter and radiation), hence and where p 1 = w 1 ρ 1 and p 2 = w 2 ρ 2 . Moreover, assume the second fluid does not interact either with the first fluid or with (H ). Hence, the conservation law for this component where we have set w 2 = const. Therefore, it is easy to show that we need to introduce the evolution of the second fluid into Eq. (12), namely Notice, that the pair (w 1 , w 2 ) lies in the region (−1, 1).
In the rest of the paper we proceed to determine analytical solutions of the Friedmann equations (13) and (18). Since these equations have a nonlinear nature we apply the algorithmic method of singularity analysis due to Ablowitz, Ramani and Segur (ARS) [57][58][59]. This method provides information about the existence and stability of singular cosmological solutions in different eras. In the limit where the constant λ 0 vanishes with n < 0, the closed-form solution of (12) was found by Perico et al. [26] (see also [70][71][72]). However, the existence of singular solutions close to the singularity has not yet been studied. As we shall see below, the constant term, λ 0 , modifies the solution close to the singularity and so affects the entire cosmic history.

De Sitter phases
Consider the scenario with the two ideal fluids, ρ 1 (t) and ρ 2 (t). We choose the dimensionless variables 1 3H 2 , and after some algebra the field equations reduce to the following system of first-order equations: while the constraint Eq. (3) becomes In order to have a de Sitter phase the corresponding critical point P = H P, P of the dynamical system (19), (20) needs to obeyḢ (P) = 0. This condition implies Therefore, if we consider that H P = 0, which is always true for λ 0 = 0, then 1 = 0, and so we have Moreover, combining the equation with Eq. (21) we obtain 2 = 0, which implies that the vacuum component (H ) ensures the existence of de Sitter points. For example, in the simplest scenario where n = −1, Eq. (22) provides the following two critical points In order to study the stability of the solution around the critical points P 1 , and P 2 , we linearize the dynamical system by substituting in the system (19), (20), with ε 2 → 0. Next, we determine the eigenvalues of the linearized system and conclude that the point is stable if all the eigenvalues have negative real parts.
For n = −1, the linearized system around the critical points iṡ The eigenvalues of the linearized system are functions of the parameters w 1 , w 2 and λ 0, λ 2 . In Figs. 2 and 3 we present the regions of the parameter space, {λ 0 , λ 2 }, for which the points P 1 , P 2 are real and the eigenvalues have negative real values. The plots are for specific values of the parameters w 1 , and w 2 . The common area of the regions provides the appropriate pairs, {λ 0 , λ 2 } , for which the de-Sitter points are stable. Now for n = 1, Eq. (22) admits the three critical points Fig. 2 Plot of the region in the space {λ 0 , λ 2 } where the eigenvalues of the linearized system (19 ), (20) are negative. The first line is for w 1 = 0 and w 2 = 1 3 , while the second line is for w 1 = 1 3 and w 2 = 0. The plots in the left column corresponds to the point P 1 , and the plots in the right column to the point P 2 . The points are stable only when the colored areas overlap Otherwise, when < 0, only P 1 is a real point and describes always a de Sitter Universe.
Furthermore, for n = 2, the algebraic Eq. (22) admits the four solutions, which are de Sitter points for values of λ 0 , λ 2 for which the solutions are real. Easily, it follows that when λ 0 , λ 2 are positive, then points P 2± are always real.

Singularities and analytical solutions
Technically , in order to study the existence of matterdominated eras, that is, singular solutions of the form a (τ ) = a 0 τ p , τ = t − t 0 , we apply the ARS algorithm, hence the analytical solution can be written as a Painlevé Series around the (movable) singularity. Moreover, the ARS algorithm provides important information regarding the stability of the singular solution a (τ ) = a 0 τ p . The approach of the ARS algorithm is inspired by the work of Kowalevskaya [73], where the solution of a differential equation is described by a power-law function when the differential equation admits a movable singularity. Hence, a necessary condition for the approach to succeed is the existence of at least one movable singularity. The position of the singularity is t 0 , and the singulariy is characterized movable when t 0 depends on the initial conditions.
The ARS algorithm is developed via a three-step process, a brief description of which is Step A: The determination of the leading-order behavior. The leading-order term needs to be a negative integer, or at least a non-integral rational number. It is important to point out that the coefficient of the leading-order term may or may not be determined explicitly.
Step B: The proof of the existence of sufficient number of (arbitrary) integration constants by determining the resonances (Fuchs indices). The value minus one should always appear always of one of the resonances. This is important for the singularity to be movable and hence the ARS algorithm to be valid.
Step C: The consistency test is performed by substituting an expansion of the power-law series which describes the solution of the original equation, to test that it is indeed a true solution.
From the leading-order term and the resonances, we can determine the step of the Laurent expansion (Painlevé series) which describes the actual solution of the original equation. Whereas, the coefficients of the Laurent expansion are determined from the consistency test. For a right Painlevé series the resonances must be positive, for a left Painlevé series the resonances must be negative, while for a full Laurent expansion the resonances have to be mixed. Clearly, in the case of an autonomous second-order ordinary equation, since there is only one free resonance, the possible Laurent expansions are either left or right Painlevé series. For more details on the method we refer the reader to the review article of Ramani et al. [74].

Case A: one ideal fluid
In this case we consider one ideal fluid. In order to determine the movable singularities of Eq. (13) we perform the change of variables, a (t) → A −1 (t), and multiply the equation with the term A 2 (t)Ȧ n (t) and we finḋ Then we insert A (τ ) = A 0 τ p in the above equation, where τ = t − t 0 , and we search for the dominant terms in order to determine the power p. Notice, that t 0 is an integration constant and it denotes the position of the singularity.

Positive power, n > 0
For n > 0 and in the context of A (τ ) = A 0 τ p Eq. (26) then takes the form from which it follows that the only possible dominant term is that of p n+1 [2 + 3 (1 + w 1 ) p]. Therefore if we assume that the leading-order behavior describes explicitly the solution at the singularity then the latter algebraic equation reduces to Hence, the singular solution of the ideal fluid, without the -varying term, describes the actual solution of the model close to the singularity.
For the second step of the ARS algorithm, we substi- The second resonance tells us that the second integration constant is the coefficient A 0 , of leading-order term. We show that Eq. (28) is not affected by A 0 , hence this constant is arbitrary. In this case the analytical solution of Eq. (26) is given by a right Laurent expansion, while its step depends on w 1 . Notice, that the presence of a right Painlevé series means that the matter-dominated era close to the singularity, is an unstable solution. 2 We note that the problem contains two integration constants, namely A 0 and w 1 , hence it is not necessary to perform the consistency test. Below, we complete the current analysis by presenting some analytical solutions which are physically interesting. (a) Analytic solution for a non-relativistic fluid: For a pressureless matter source, we have w 1 = 0. Therefore, from Eq. (29) we calculate p = − 2 3 , which means that the Laurent expansion has a step of 1 3 , that is, the solution of (26) is where we need to determine the coefficients A 1 , A 2 , .. etc.
In the case of the −varying model (11), with n = 1, the first six non-zero coefficients of (31) are: A 0 , which is an arbitrary integration constant, and while the rest of non-zero coefficients are those of A 21+3ν , ν ∈ N.On the other hand, for other values of the power n, the coefficients are different. For example, in the case of the (H ) model with n = 3, aside from A 0 , the first non-zero coefficients are while the other non-zero coefficients are A 18+3ν , ν ∈ N.
From the aforementioned solutions, we observe that the solution close to the matter-dominated era includes an additional term A 0 τ 4 3 , namely and so the scale factor is written as a Taylor expansion around τ = 0 as follows, (b) Analytic solution for a radiation fluid: Considering that the ideal fluid is radiation (w 1 = 1/3) from (29), we obtain p = − 1 2 , which means that the solution is expressed by a right Painlevé series with step 1 2 , that is, whereĀ 0 is an integration constant. For n = 1, the nonzero coefficients of the expansion (34) arē .N.
For n = 2, the coefficients of the expansion arē Therefore, the analytical solution prior to the radiation era is well approximated by the following expression: and so

Negative power, n < 0
Now we study the case where the exponent n in (H ) = λ 0 + λ 2 H −n is negative. We remind the reader that we are still using the variable A (τ ) = A 0 τ p . Under those conditions Eq. (26) gives Following similar notations to those of Sect. 4.1.1 we find that the only possible leading term is the one with coefficient 2 + 3 (1 + w 1 ) p as long as n > −2, which implies that the exponent p is given by Eq. (29). As expected, from the ARS algorithm, we find that the constant A 0 is arbitrary, hence the corresponding resonances (30) are s 1 = − 1 and s 2 = 0 respectively, whilst the analytical solution is written as a right Painlevé series. Specifically, we find the following interesting solutions.
(a) Analytic solution for a non-relativistic fluid: For w 1 = 0, the Laurent expansion is of the form of (31), where with n = −1, namely (H ) = λ 0 +λ 2 H , we obtain the following coefficients: .N This means that close to the singularity the scale factor is approximated by Notice that the second term of this solution is different with that of Eq. (33) for n > 0. The closed-form solution in the case of the (H ) = λ 0 + λ 2 H running-vacuum model has been found earlier in [68]. In this case the running vacuum model introduces a mild dynamical behavior of the vacuum energy at low energies, namely dark matter and dark energy eras. For more details regarding the dynamics as well as the corresponding observational constraints of the current vacuum model we refer the reader the work of [68]. On the other hand, terms of the form H 3 , H 4 (n < −2) they can be important for the early universe (see [26,76]). (b) Analytic solution for a radiation fluid: Here, the equation of state parameter (hereafter EoS) is w 1 = 1 3 , and the analytical solution is given by the series (34). For n = −1, the corresponding non-zero coefficients arē and so, prior to the radiation-dominated era (singular solution), the scale factor can be written as Again, this solution includes an extra term which differs from Eq. (35 ) for n > 0. At this point the following comments are appropriate. The current approximate solutions (37), (38), for non-relativistic matter and radiation contain extra terms which differ from those in (33) and (35) where n = 1. This is expected: for negative powers of n, the extra terms of the solutions are related to the λ 2 H −n term in Eq. (11), while for n > 0 the extra terms are affected by the constant vacuum term λ 0 . However, in this case the λ 0 term is introduced in the secondorder correction to the solution, which means that there are  Fig. 4. For the numerical integration we have considered λ 2 = 1, n = −1, while the solid line provides the relative error of the solution with λ 0 = 10 −1 λ 2 , the dash-dash line for λ 0 = 10 −2 λ 2 , the dot-dot line for λ 0 = 10 −3 λ 2 and the dash-dot line for λ 0 = 10 −4 λ 2 . We started the numerical integration close to the singular solution with w 1 = 1 3 , that is, at the radiation epoch differences between the two solutions λ 0 = 0 and λ 0 = 0, not very far from the singularity. In order to understand the differences between the two solutions, in Fig. 4 we present the evolution of the relative difference between the numerical simulations for n = − 1, and λ 0 = 0 and those with λ 0 = 0, for the scale factor and the Hubble expansion rate.
We observe that prior to the singularity the relative error is close to the error of the numerical integration; however, as the system evolves the relative error becomes larger. In the left panel of Fig. 5 we plot the evolution of the scale factors which are considered in Fig. 4, while in the right panel of Fig. 5 we show the relative difference in the energy density = (H ) /3H 2 between the different solutions. Recall that λ 0 does not correspond to the value of the cosmological constant at the present time. In particular, the latter is given by (H 0 ) = λ 0 + λ 2 H −n 0 , where H 0 is the value of the Hubble function at the present time. Therefore, since we know that (H 0 ) is small, we must have λ 0 of the order of λ 2 H −n 0 : that is, for n < 0, |λ 0 | > |λ 2 | . This makes clear that the correction terms which depend on λ 0 are important for the evolution of the system.
The singularity analysis fails for n < −2, and below we will carry out our analysis for the case of two ideal fluids.

Case B: two ideal gases
In this section we consider that the cosmic fluid consists two ideal fluids, where one of them is minimally coupled (see Case B in Sect. 2). In this case using the transformation A (t) = a −1 (t), the second-order differential equation (18) becomeṡ Using τ = t − t 0 and substituting A (t) = A 0 τ p in (39), we find the expression This provides the following three dominant behaviors:

Positive power, n > 0, solution I
For n > 0, if we assume that the dominant term close to the singularity is 2 + 3(1 + w 1 ) p, then p is given by Eq. (29). This implies that the singular solution describes the matterdominated era via the fluid ρ 1 (t). Of course we need to clarify that the above term is indeed the dominant one if and only if the corresponding EoS parameters of the two ideal fluids satisfy the inequality Within the framework of this latter restriction we present some special analytical solutions, below. It is important to mention that when w 2 = − 1 3 , Eq. (18) corresponds to that of one ideal effective fluid in a FLRW spacetime with nonzero spatial curvature k, such that k = −ρ 20 .
For the calculation of the resonances, we follow standard steps; namely, we substitute A (t) = A 0 τ Because the second resonance has the value zero, the second integration constant is the coefficient of the dominant term; hence, it is not necessary to perform the consistency test. However, in order to compare the current solutions with those of Case A we shall now investigate some applications. (a) Analytic solution for a non-relativistic fluid and curvature term: Here, we have (w 1 , w 2 ) = (0, − 1 3 ). Therefore, the leading-order term is p = − 2 3 and the step of the expansion is 1 3 . The analytical solution is given by expression (31) while, for n = 1, the nonzero coefficients of the Laurent expansion are given by Notice that the other non-zero coefficients of the series have the form A 10+2ν , ν ∈ N. Recall that the constraint (3) provides an algebraic relation between the free parameters of the model. As expected, for k = ρ 20 = 0 the current solution reduces to that of section 3.1, while close to the singularity, and for k = ρ 20 = 0, the approximated solution is or (b) Analytic solution for radiation and non-relativistic matter: Now, we assume that ρ 1 (t) and ρ 2 (t) are the radiation (w 1 = 1 3 ) and matter (w 2 = 0) densities respectively. In this case, the solution of the problem are the series (34) withĀ 0 arbitrary and for n = 1 the nonzero coefficients are found to bē while the other non-zero coefficients take the form A 10+2ν , ν ∈ N.
(c) Analytic solution for radiation and curvature term: Here we consider that the cosmic fluid contains radiation (w 1 = 1 3 ) and curvature (w 2 = − 1 where we clearly observe the differences with respect to that of the previous solution (radiation and non-relativistic matter).

Positive power, n > 0, solution II
The second class of solutions are those for which the leadingorder terms of Eq. (40) and so Therefore, that new singular solution describes the matterdominated era of the minimally coupled fluid ρ 2 (t).
In contrast to the aforementioned cases, the constant A 0 here is not arbitrary but it is given in terms of w 1 , w 2 and ρ 20 , by where w 2 = 0. Since A 0 is not arbitrary, we expect that the second resonance cannot be zero. Indeed, using our methodology we find the following resonances: Unlike the previous cases, where the corresponding solutions are written as right Painlevé series, we have a different situation here. In particular, if we select the EoS parameters (w 1 , w 2 ) such that s 2 < 0, then the solution is given by a left Painlevé series and thus the dominant term is an attractor. 3 (a) Analytic solution for radiation and curvature: Using w 1 = 1 3 and w 2 = − 1 3 Eqs. (45) and (47) provide p = −1 and s 1 = −1, s 2 = 1. Therefore, the analytical solution is 3 For more details, we refer the reader to [77][78][79] and references therein.
Notice, that the solution passes the consistency test. Moreover, the fact that the second resonance is positive implies that the solution (48) is unstable. Now, using the dominant term A A 0 τ −1 of the above solution, i.e. a(τ ) ∝ τ, we find that the Milne universe and the R h = ct model [80] can be seen as perturbations around the present (H ) model.
If the R h = ct or Milne models were somehow the effective behavior of the current (H ) model, then we could argue that the R h = ct and Milne models are both unstable.
To this end, for n = 1 the nonzero coefficients of (48) are, A 1 which is the second integration constant and where A 1 is a second integration constant.

Negative power, n < 0
The last case that we have to consider is when the power n is negative, i.e. n < 0. In this case it follows that the only possible dominant behavior is that with p = − 2 3(1+w 1 ) , while the corresponding resonances are s 1 = − 1 and s 2 = 0, which means that the second integration constant is the A 0 coefficient. Thus, since we know the two integration constants explicitly it is not necessary to perform step C of the ARS algorithm, i.e. the consistency test. However, for the completeness of our analysis we present the Laurent expansions for two cases of special interest. (a) Analytic solution for radiation and nonrelativistic matter If we assume that w 1 = 1 3 and w 2 = 0, then the dominant behavior is p = − 1 2 , and the solution is given as a right Painlevé series with step 1 2 , that is, series (34), where now the first nonzero coefficients for n = − 1 arē (b) Analytic solution for radiation and curvature: In a similar way, for w 1 = 1 3 and w 2 = − 1 3 , the solution is given by the right Painlevé series (34) with step 1 2 , where the first nonzero coefficients for n = − 1 arē In the following section we discuss our results and we draw our conclusions.

Conclusions
In the framework of running vacuum ( ) models we implemented a critical point analysis in order to study the existence and the stability of singular solutions which describe de Sitter, radiation and matter dominated eras. We used one of the most popular expressions for the vacuum parameter which describes a large family of running vacuum models, namely (H ) = c 0 + c 1 H 2 + c 2 H −n . We showed that this functional form of (H ) allows for a number of de Sitter eras, depending on the model parameters. The radiation and matter dominated eras are also investigated with the method of singularity analysis. In order to recover a sequence of radiation and matter eras, we found some interesting constraints on the cosmic fluid.
Furthermore, it is interesting to mention that a closed form solution exists only under special conditions. Indeed, analytical solutions (see [26,68]) are possible when the cosmic fluid contains only two components, namely one fluid (nonrelativistic or radiation) and running vacuum (H ). Also, in the case of [26] we have an additional restriction towards solving analytically the Friedamnn equations, namely we are dealing only with even powers in H , hence the term λ 2 H −n boils down to λ 2 H 2k with n = −2k and k ≥ 2 [see Eq. (7)]. Despite the fact that the above special solutions are available, we are unaware of any study concerning the critical points in running vacuum cosmologies. Moreover, if the cosmic fluid includes radiation, non-relativistic matter and running vacuum (H ) (this situation is close to reality) then analytical solutions are not yet available, implying that the approximated solutions of the present study are completely new and novel.
In this context, we found families of (H ) cosmologies for which new analytical solutions are given in terms of Laurent expansions. Finally, we showed that the Milne universe and the R h = ct model can be written as perturbations deviating from a special, but unstable, (H ) model.