Cosmological singularities and analytical solutions in varying vacuum cosmologies

We investigate the dynamical features of a large family of running vacuum cosmologies for which $\Lambda$ 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 $\Lambda(H)$ cosmologies for which new analytical solutions are given in terms of Laurent expansions. Finally, we show that the Milne universe and the $R_{h}=ct$ model can be seen as perturbations around a specific $\Lambda(H)$ 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).
An alternative approach to resolving these two problems is to consider the so called Λ (t)CDM models, wherein Λ is allowed to vary with cosmic time (see [31,32,70] and references therein). This class of models [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] is based on a dynamical vacuum energy density that evolves as a power series of the Hubble rate (for review see [28], [74]). Powered by a decaying vacuum energy density, the spacetime can emerges from a pure non-singular initial de Sitter vacuum stage, "gracefully" exit from inflation to a radiation era followed by dark-matter and vacuum-dominated phases, before finally, evolving to a late-time de Sitter phase [32,61,70]. Recently, Sola et al. [72] tested the performance of the running vacuum models against the latest cosmological data and they found that the Λ(H) models are favored than the usual ΛCDM model at ∼ 4σ statistical level (see also [73]). These developments have led to growing interest in Λ(H) cosmological models.
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 [52][53][54]. This algorithm provides a method to construct the analytical solution of a differential equation which is expressed as a Laurent expansion around the singular leadingorder term (for some applications on gravitational theories see [55][56][57][58][59][60] 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 Section 2 we briefly introduce the concept of the running Λ(H) cosmologies. In sections 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 section 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 tensorT µν ≡ 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. Λ = Λ(χ(t)). Moreover, considering G = const. (for other scenarios see [74]) the Bianchi identity that insures the covariance of the formulation, implies an energy exchange between ρ Λ (or Λ) and ρρ +Λ + 3H (ρ + p) = 0.
The above equations (3)-(5) are not independent, since from Eq.(4) and Eq.(3) we can extract Eq.(5). We can see there exists a simple exact scaling solution in which all the terms in (3) and p = (γ − 1)ρ, with γ constant, are proportional to t −2 . This gives an exact solution of (3)- (5): Hence, the ratio of the dark energy density to the other matter densities are always constant, and given by Observationally, this ratio is about Λ ρ ≃ 7/3 which, in the dust (γ = 1) era, would require h = 20/9.This solution contains only one parameter (h) after the equation of state (γ) is fixed and it has the interesting 'everlasting lambda' property that the cosmological constant is always important because it is always of similar order to that of ρ. Another scenario of this sort has been proposed but turned out to make predictions about spatial fluctuations in conflict with current observations [23]. A different approach that does not have this problem is introduced in refs. [24,25].
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 [26,28,29]. Specifically, in ref.
[32] 2 the following expression (for recent review see [69]) was proposed for the functional form of Λ(H), 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 [69,70]. 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 [69,71,72] and references therein). Furthermore, for a general quadratic vacuum model, with Λ (H) = c 1 H 2 + F (H) , we find where F (H) is an arbitrary function. Introducing the following parametrizations, Eq.(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−c1 3 λ 0 and c 2 = 3−c1 3 λ 2 . 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 or, with the aid of H =ȧ/a, we find 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 impliesρ 2 + 3(1 + w 2 )Hρ 2 = 0, a solution of which is 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) [52][53][54]. 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. [32] (see also [75], [76]). 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 (t) = 3H 2 ρ 1 (t) , Ω 2 (t) = 3H 2 ρ 2 (t), 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. (1) and (2) 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.

SINGULARITIES AND ANALYTICAL SOLUTIONS
Technically, in order to study the existence of matter-dominated 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 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 [62], 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 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. 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. [63].

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.
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 or p ≃ − 2 3 (1 + w 1 ) .
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 substitute A (τ ) = A 0 τ − 2 3(1+w 1 ) + mτ − 2 3(1+w 1 ) +s , and we linearize around m = 0. From the remaining terms, we find that the coefficient of the leading-order term vanishes if and only if s (s + 1) = 0 (30) which means that the resonances are s 1 = −1 and s 2 = 0. The second resonance tells us that the second integration constant is the coefficient A 0 , of leading-order term. We show that equation (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 3 . 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 4.1.2. Negative power, n < 0 Now we study the case where the exponent n in Λ(H) = λ 0 + λ 2 H −n is negative, that is n = − |n|. 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 section 4.1.1 we find that the only possible leading term is 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:

38880
.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 [69]. 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 nonrelativistic 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 second-order correction to the solution, which means that there are 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. 3 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.   Fig.) and the relative difference of the energy density ΩΛ (t) for the solutions which presented in Fig. (3). For the numerical integration we considered λ2 = 1, n = −1, while the solid line provides the relative error of solution with λ0 = 10 −1 λ2, dash-dash line for λ0 = 10 −2 λ2, dot-dot line for λ0 = 10 −3 λ2 and the dash-dot line for λ0 = 10 −4 λ2. For the numerical integration we started close to the singular solution with w1 = 1 3 , that is, at the radiation epoch.
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. 4 we plot the evolution of the scale factors which are considered in Fig. 3, while in the right panel of Fig. 4 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 section 2). In this case using the transformation A (t) = a −1 (t), the second-order differential equation (18) becomesȦ Using τ = t − t 0 and substituting A (t) = A 0 τ p in (39), we find the expression 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 matter-dominated 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 (39) and from the coefficient of the dominant term in the linearized equation around the value m = 0 we obtain s (1 + s) = 0, hence the two resonances are −1 and 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 a (τ ) = a 1 τ 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 3 ). Again, if in the case of the dust we consider another fluid, for instance w 2 = − 1 3 , which corresponds to the curvature term, then the analytical solution is given by Eq.(34) with the following nonzero coefficients (for n = 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 leading-order terms of Eq. (40) are −3 + (3 + n) p and −1 + p [3 (2 + w 2 ) + n]. Keeping only these terms, Eq. (40) reduces to and so .
Therefore, that new singular solution describes the matter-dominated 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. 4 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 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 [77] 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. 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+w1) , 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. In this context, we found families of Λ(H) cosmologies for which new analytical solutions are given in terms of 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.