Viscous effect in the late time evolution of phantom universe

We investigate the cosmological implications of a phantom dark energy model with bulk viscosity. We explore this model as a possible way to resolve the big rip singularity problem that plagues the phantom models. We use the latest type Ia supernova and Hubble parameter data to constrain the model parameters and find that the data favor a significant bulk viscosity over a non-constant potential term for the phantom field. We perform a dynamical analysis of the model and show that the only stable and physical attractor corresponds to a phantom-dominated era with a total equation of state that can be greater than $-1$ due to the viscosity. We also study the general effect of viscosity on the phantom field and the late time evolution of the universe. We apply the statefinder diagnostic to the model and find that it approaches a nearby fixed point asymptotically, indicating that the universe can escape the big rip singularity with the presence of bulk viscosity. We conclude that bulk viscosity can play an important role in affecting the late-time behavior as well as alleviating the singularity problem of the phantom universe.


I. INTRODUCTION
Cosmological observations in recent decades have revealed more details about the evolutional history of the universe.One of the major discoveries is the late time acceleration of cosmic expansion, which can be generally understood from at least two different perspectives.Either Einstein's theory of gravity is incomplete and calls for modifications (for recent reviews, see, e.g., Refs.[1,2]), or there is an unknown form of energy dubbed dark energy (DE) that exhibits repulsive behavior.The simplest model of DE is the cosmological constant Λ model with cold dark matter (ΛCDM), which assumes that DE consists of a constant energy density that fills the space uniformly and has a constant equation of state (EoS) w DE = −1, where w DE is the ratio of pressure to energy density of the DE content.
However, the cosmological constant has no clear physical origin and faces several theoretical challenges (see, e.g., Ref. [3]).Therefore, many alternative DE models have been proposed and explored (see, e.g., Refs.[4][5][6]).These models can be roughly divided into two categories: the quintessence model with w DE > −1 [7,8] and the phantom model with w DE < −1 [9,10], which have different implications for the ultimate fate of the universe.
Quintessence models generally predict that the universe will undergo an eternal expansion, while phantom models usually indicate that the universe will enter a super-accelerated expansion phase and end in a finite time with a cosmic singularity known as big rip, where all structures will be torn apart.Recent observations seem to favor the phantom model over the quintessence model [9,[11][12][13][14].Nonetheless, such discussions about w DE rest upon the assumption that it is simply given by p DE /ρ DE .This representation is purely phenomenological and lacks underlying physics.To provide a more comprehensive understanding of the nature of DE including its temporal evolution or perturbation, a more fundamental framework should be considered.One approach involves introducing a dynamic scalar field as part of either the quintessence or phantom models, depending on the specific properties of its potential and kinetic terms [9].In particular, the phantom scalar field model has attracted significant interests within the literature (see, e.g., Refs.[4,13,15] and the references therein).
The studies of phantom models typically lead to discussions about the big rip singularity [16][17][18][19][20][21][22][23][24][25].To address this issue, researchers have introduced modifications to DE based on quantum effects [26], geometric effects derived from modified gravities [27,28], or interactions between cosmic content [29,30].One way to introduce interactions involves considering the viscosity of each component.This approach accounts for the dissipative properties of real fluids.In particular, bulk viscosity is the most relevant for cosmology since we assume a homogenous and isotropic universe.This can be incorporated into the standard cosmological scheme by redefining the effective pressure p eff = p − Π with a viscous pressure term Π to restore thermal equilibrium [31][32][33].
Invoking viscosity in cosmology is proven to be useful to resolve or soften the cosmic singularity problem in different models [34][35][36][37][38][39][40][41].Following this approach, it is shown that the singularity problem can be alleviated in the anisotropic phantom universe with viscosity [25].
Moreover, an interacting phantom DE with dark matter induced by the viscous approach is shown to be able to avoid the big rip singularity [17,18].
In this work, we will study viscous phantom scalar field DE model and explore its late time behavior.The paper is structured as follows.In Sec.II, we review the phantom scalar field model of DE and introduce the viscous cosmology framework.Section III contains the late time observation fit to constrain the model parameters.Dynamical system analysis and statefinder diagnostic of the model are given in Secs.IV and V, respectively.And we conclude our study in the last section.

II. VISCOUS PHANTOM SCALAR FIELD DE MODEL
The action for a phantom field minimally coupled to gravity is given by where V (ϕ) is the potential of the phantom field ϕ.The energy density and pressure of the phantom field are [9] We assume a spatially flat Friedmann-Lemaître-Robertson-Walker metric for the homogeneous and isotropic universe, given by where a(t) is the scale factor.Then, the Friedmann equations read where ρ m and p m are the energy density and pressure of dust matter, respectively, H = ȧ a is the Hubble parameter, and the dot represents the derivative with respect to the cosmic time t.
As mentioned in the Introduction, the bulk viscosity dissipation in the cosmic phantom field fluid can be represented by a pressure term −Π added to p ϕ .For bulk viscosity related to the cosmic expansion, we assume that Π ∝ H, which implies that the effective pressure of the phantom field is where ξ ϕ is the viscosity coefficient.Then, the evolutionary equations for dust matter and the phantom field are given by ρm + 3Hρ m = 0 , Using Eq. ( 2), we obtain the equation for the phantom field where the prime denotes the derivative of the potential V with respect to the field ϕ.
Generally, the viscosity of a fluid may depend on energy density, pressure, spacetime geometry and so on.For simplicity, in the current work, we assume the viscosity is proportional to the changing rate of the phantom field, i.e., where the viscosity coefficient ξ 0 is a constant parameter.As for the potential of the phantom field, it is shown that an exponential form of potential can match or account for the acceleration of expansion [42][43][44][45][46][47], which is then adopted in the current work and is given by where α 0 and V 0 are constants.

III. OBSERVATIONAL CONSTRAINTS
We constrain the viscous phantom model using the late-time observational data sets based on Markov Chain Monte Carlo method.We use the Pantheon compilation of Type Ia supernova (SNIa) [5] and Hubble parameter (H(z)) data points [48]  One can see that for Model P, the phantom potential cannot be constant at 1σ confidence level as the phantom field is the sole mechanism in operation for the late-time cosmic acceleration.However, if both the viscosity and the phantom field with its potential are taken into consideration, as in Model vP, the fitting result favors a significant viscosity and does not exclude the constant potential case with α 0 = 0.Where γ = V 0 /H 2 0 .

IV. DYNAMICAL ANALYSIS
By using the dimensionless variables [42,43,51] the partial differential equations of the phantom scalar field can be recast into an autonomous system as where N = ln a = − ln(1 + z).The Friedmann equation can also be rewritten as a dimensionless constraint where Ω i = ρ i 3H 2 is the density parameter of the corresponding component.The critical points of the above system and their existence conditions are listed in Table II.The stability of each critical point can be inferred from the eigenvalues of the linearized ) all ξ 0 and α 0 ̸ = 0 1 TABLE II.The critical points and their existence conditions of autonomous system (11).
perturbation matrix M of the autonomous given by •For critical point A 1 , the eigenvalues of the linearized perturbation matrix are η (1) where In Fig. 2, the colored areas mark the possible regions of parameters for the critical point A 1 to exist.The stable (pink) and saddle (blue) regions are divided by the line Ω ϕ = 1.
•For critical point A 2 , the eigenvalues of the linearized perturbation matrix are η (1) A 2 exists as long as α 0 ̸ = 0.The sign of the first eigenvalue η A 2 is bifurcated by the hyperbola α 0 ξ 0 = −1/2, while the second eigenvalue η (2) A 2 is always negative for any real value of parameter ξ 0 .Therefore, the stable (pink) and saddle (blue) regions of A 2 are shown as Fig. 3.However, we note that A 2 only exists when Ω ϕ < 0. This scenario is also known as the perfect fluid supra-dominated era [52], which is beyond the scope of the current work and is considered unphysical.We constrain the parameters α 0 and ξ 0 as such that 0 ≤ Ω m , Ω ϕ ≤ 1 and this point A 2 will not exist.•For critical point A 3 , the eigenvalues of the linearized perturbation matrix are η (1) A 3 also exists for α 0 ̸ = 0 and any real ξ 0 .The second eigenvalue η A 3 is always negative, while the first one can have either sign.The stable (pink) and saddle (blue) regions of A 3 are shown in Fig. 4. We note that at this point the phantom field energy density parameter Ω ϕ = 1, which is invariant under the variation of the parameters α 0 and ξ 0 .So, this point represents a phantom field dominated universe and is always physical.
The saddle and physically relevant regions of parameters for A 1 are subsets of the stable regions of A 3 , but only under the circumstance that A 1 exists (i.e., the pair of parameters lies in the colored region depicted in Fig. 2).In this scenario, the system will evolve from A 1 to A 3 in the cosmic history since A 2 does not exist for these parameters.However, it is also possible that the parameters do not allow for the existence of A 1 (i.e., they fall outside the colored region depicted in Fig. 2), yet satisfy the physical conditions 0 ≤ Ω m , Ω ϕ ≤ 1.Then, A 3 becomes the sole existing critical point and the system will simply converge towards it.
The best-fit result of Model vP corresponds to the second scenario.Moreover, the total EoS of the cosmic fluid can be expressed in terms of x, y and ζ as The universe is in accelerating phase if −1 < w tot < −1/3, and decelerating if −1/3 < w tot < intercepts with the blue hyperbola y 2 − x 2 = 1 only at one point (x, y) = (0, 1), which means that either the model needs severe fine-tuning or it will encounter big rip singularity.This is independent of the fitting.For Model vP, the regions w tot ∈ [−1, −1/3] and w tot ∈ [−1/3, 0] calculated from the best-fit result are painted cyan and yellow in Fig. 5, respectively.One can see that the cyan region covers a portion of the blue hyperbola that represents all possibilities of A 3 .The best-fit A 3 is in this region.Therefore, the presence of bulk viscosity can indeed allow the model evolve into an attractor that is still in the w tot > −1 region, and help avoid the cosmic singularity.
A heteroclinic orbit representing the evolutionary history of our universe most likely starts somewhere near the green hyperbola (matter dominated era) in the yellow region (decelerating phase with −1/3 ≤ w tot ≤ 0), crosses the purple hyperbola (present time) in the cyan region (accelerating phase), and eventually reaches some point on the blue hyperbola (attractor A 3 ).However, the preceding section of such a heteroclinic orbit seems to cross from the Ω ϕ < 0 region to Ω ϕ > 0 region in early time and encounter a singularity of the phantom field EoS w ϕ = p eff /ρ ϕ .This issue arises because our analysis relies on late time observation data, which do not capture the early stages of cosmic phase transitions.
Further research may be needed to address this typical concern rooting from the negative kinetic energy of phantom field.

V. STATEFINDER DIAGNOSTIC
The statefinder diagnostic uses the parameters {q, r, s} that are derived from higher derivatives of the scale factor, to differentiate various DE models.For the current model, the parameters are Using the best-fit parameters, we plot the r-q and r-s trajectories for both Model P and Model vP in Figs. 6 and 7, respectively.The present time is indicated by the rounded marker on the curves.As shown in Fig. 6, Model P exhibits a monotonically increasing deviation from the de Sitter point in the future and asymptotes to a point in the region with q < 0, r > 1.In contrast, Model vP will asymptotically approach a stable fixed point (q, r) = (−0.9794,0.9391) near the de Sitter point.Fig. 7 reveals that both models have r < 1 and s > 0 in the early universe, and they both cross the ΛCDM point (0, 1), but have different behaviors afterwards.Model P moves away from the ΛCDM point, while Model vP returns to it and passes it again before converging to a nearby fixed point (s, r) = (0.0137, 0.9391) in the far future.
The effective EoS parameters of the total cosmic fluid, w tot , for both Model P and vP are plotted in Fig. 8.As shown in the figure, w tot of Model vP remains above −1 and asymptotes to −0.9863 in the infinite future, which implies that the universe will avoid the big rip singularity.This result is consistent with the dynamical analysis in the previous section.Moreover, the effective EoS parameter, w ϕ , of the viscous phantom field is also We also investigate the effect of viscosity on the evolutionary history by computing the ratio K of the effective pressure to the intrinsic pressure of the phantom field, which is defined as Fig. 9 shows the evolution of K as a function of ln a.It is evident that viscosity has a significant impact on reducing the pressure of the phantom field, especially around z = 0.8475 (ln a = −0.6138),where K reaches a minimum value of about 0.4921 and the pressure of the field is only a half of the intrinsic pressure.At present (z = 0), the effective pressure is about one third (32.73%) lower than the intrinsic pressure due to viscosity.In the future, viscosity will continue to reduce the intrinsic pressure by about one quarter, as K asymptotes to 0.7672.

VI. CONCLUSION AND DISCUSSIONS
Recent cosmological observations suggest that the late-time acceleration of the universe is more likely driven by phantom DE rather than quintessence DE.However, phantom DE models typically encounter the big rip singularity problem.In this work, we explore the possibility of avoiding the big rip singularity by introducing bulk viscosity in the phantom field model.We constrain the model parameters using the latest SNIa and H(z) data.The results indicate that the data favor a significant bulk viscosity over a non-constant potential term for the phantom field.
We then perform a dynamical analysis of the viscous phantom field model using the best-fit values of the parameters.We find that the only stable and physical attractor of the autonomous system is A 3 , which corresponds to a phantom-dominated era.The other critical points, may that be either unstable, unphysical, or cannot exist for the best-fit parameters, do not represent different cosmological eras, since the two variables of the autonomous system only correspond to the kinetic and potential energy.We also plot the possible curves that represent the matter-dominated era and the present epoch on the phase portrait.The heteroclinic orbit that describes our universe is expected to cross these curves and converge to the A 3 attractor.However, we cannot trust the early part of the orbit or the critical points that precede the matter-dominated era, since we only fit the model for the latetime behavior.Due to bulk viscosity, some combinations of the parameters may allow the A 3 attractor have a total EoS w tot of cosmic fluid that falls within the accelerating region −1 < w tot < −1/3.This implies that the phantom universe may avoid the big rip singularity and end in a state with w tot > −1 with the presence of bulk viscosity.
We apply the statefinder diagnostic to the viscous phantom field model and compare it with the non-viscous model.We find that the statefinder parameters of the viscous model do not diverge monotonically from the ΛCDM point, but rather approach a nearby fixed point asymptotically.The bulk viscosity acts as a dissipative force that lowers the intrinsic pressure of the phantom field.The viscosity has its maximum effect around z = 0.8475, where it lowers the intrinsic pressure by about a half.At the present time (z = 0), the viscosity reduces the intrinsic pressure by a third, resulting in a lower effective pressure.In the asymptotic future, the viscosity will still lower the intrinsic pressure by a quarter.This reduction enables the universe to escape the big rip singularity.
For simplicity, we have limited our study to a specific form of the phantom field potential and the bulk viscosity, which may not be the most general or realistic choice.It would be worthwhile to explore other forms of potential and viscosity that can fit the observational data and avoid the big rip singularity.Due to the scope of our work, we have only focused on the late-time behavior of the viscous phantom field model, and ignored the early-time dynamics that may involve other mechanisms.A more comprehensive study should include the full history of the universe and investigate the transitions between different cosmological eras.We also use conventional methods of data analysis and dynamical analysis, which may have some limitations or biases.Future research could employ artificial intelligence techniques to improve the accuracy and efficiency of data fitting, parameter estimation, and model selection [53,54].

5 FIG. 1 .
FIG. 1. Constraints on the viscous phantom model parameters from 1σ to 2σ confidence level.

FIG. 2 .
FIG. 2. The parameter regions for the existence of critical point A 1 .The shaded areas are considered unphysical due to either Ω ϕ > 1 (shaded pink area) or Ω ϕ < 0 (shaded blue area).

FIG. 3 .
FIG. 3. The parameter regions for the existence of critical point A 2 .

FIG. 4 .
FIG. 4. The parameter regions for the existence of critical point A 3 .

Figure 5
Figure5shows the phase portrait of the viscous phantom field in the late universe.The phase plane is divided into three regions: (i) Ω ϕ > 1 (the light gray region above the upper branch of the blue hyperbola y 2 − x 2 = 1); (ii) 0 < Ω ϕ < 1 (the region between the upper branch of the blue hyperbola y 2 − x 2 = 1 and red straight lines y = ±x); and (iii) Ω ϕ < 0 (the regions below the red straight lines).The blue hyperbola y 2 − x 2 = 1 represents all the possible attractors A 3 for different values of parameters {α 0 , ξ 0 }.The blue rounded marker indicates the specific attractor point A 3 for the best-fit parameters.The purple hyperbola corresponds to Ω ϕ0 = 1 − Ω m0 = 0.678, which is the best-fit density parameter of the phantom field at present time.The green hyperbola corresponds to Ω ϕ = 0.018, which is the value of Ω ϕ when Ω m reaches its maximum in the cosmic history, indicating the matter-dominated era calculated from the fitting result.

FIG. 5 .
FIG.5.The phase space evolution of the dynamical system for Model vP.

FIG. 6 .
FIG.6.Evolving trajectories of the statefinder pairs in the q-r plane.The black dot (−1, 1) represents the de Sitter phase, and the solid dots on the lines represent the current state of the models.

FIG. 7 .
FIG. 7. Evolving trajectories of the statefinder pairs in the r-s plane.The black dot (0, 1) represents the ΛCDM model, and the solid dots on the lines represent the current state of the models.

FIG. 8 .
FIG. 8.The evolution of the total EoSs, w tot , of Model P and Model vP, as well as the effective EoS, w ϕ of Model vP.

FIG. 9 .
FIG.9.Evolution of the ratio of effective pressure p ef f to scalar field pressure p ϕ .
to fit the model parameters.The SNIa data consist of 279 samples from Sloan Digital Sky Survey (SDSS) and Supernova Legacy Survey (SNLS) with redshifts 0.03 < z < 0.68, and 1048 samples with redshifts 0.01 < z < 2.3 including the Hubble Space Telescope (HST) samples and various low-z samples.The H(z) data include 26 data points from Baryon Acoustic Oscillations and 31 data points from the differential age method.For comparison, besides fitting the viscous phantom model under consideration (denoted as Model vP in the following), we also perform the same fitting procedure for the phantom model without viscosity (denoted as Model P in the following).The fitting results are summarized in TableI1 .Figure1shows the constraints on the parameters of Model vP at 2σ confidence level.χ 2 min /dof.1053.93/10991053.174/1098TABLEI.The best fitting values of model parameters and 1σ confidence level for Model P and Model vP under Pantheon+H(z) data sets.