Dynamical Stability of Bulk Viscous Isotropic and Homogeneous Universe

: In this paper, we study the phase space portrait of homogeneous and isotropic universe by taking different coupling functions between dark energy models and bulk viscous dark matter. The dimensionless quantities are introduced to establish an autonomous set of equations. To analyze the stability of the cosmos, we evaluate critical points and respective eigenvalues for different dynamical quantities. For bulk viscous matter and radiation in tachyon coupled ﬁeld, these points show stable evolution when γ (cid:29) δ but accelerated expansion of the universe for δ > 19 . The stability of the universe increases for some stationary points which may correspond to the late-time expansion for the coupled phantom ﬁeld.


Introduction
Recent astronomical observations such as type Ia supernova and cosmic microwave background radiations appear to be fairly compatible with the current cosmic accelerated expansion. These observational pieces of evidence have categorized the universe into four basic epochs. In the vacuum energy era, Planck magnitude of the universe transformed to a macroscopic immensity through a rapid expansion named as inflation. Upon reducing the temperature up to 10 3 K, the universe underwent matter-dominated epoch followed by the cosmic radiation state. Ultimately, in the current de Sitter era, the universe is in the phase of accelerating expansion. According to the latest statistics, our cosmos is composed of 95% dark sectors (dark energy (DE) and dark matter) and 5% baryonic matter [1,2]. We can observe the baryonic matter with the naked eye while the gravitational effects of other cosmic components confirm their presence.
Dark energy as an anti-gravitational force tends to accelerate the expansion of the cosmos. It has huge negative pressure with very low energy density in contrast to other cosmic components. The hidden characteristics of these dark sectors are still unknown even after the appropriate results of the ΛCDM model. There are two alternative DE model-one is to modify the geometric part of the Einstein-Hilbert action while the other one is to allocate the appropriate form of the energy-momentum tensor like scalar field models [3][4][5][6][7]. Another technique for cosmic evolution can be achieved through the transformation of the barotropic equation of state in Chaplygin gas [8].
Although a ΛCDM model is compatible with current cosmic observations but still no satisfactory argument has been achieved for the coincidence problem. Interaction of DE with dark matter is the novel technique that might address this issue. The interacting DE models have been proposed by several authors [9][10][11][12][13]. It has been found that the interaction between dark matter and DE models can assist to elaborate on the current cosmic evolution. If the DE models have Ω DE Ω DM of the order 1 and an accelerated stable solution, then the coincidence problem can be alleviated. The non-interacting DE models [14][15][16] showed late-time accelerated attractors with Ω DE = 1, which did not provide an appropriate solution for the coincidence problem. For the coupling, there does not exist any elementary theory that decides a particle form of interaction. It has been shown [17,18] that the coincidence problem can be resolved by selecting the appropriate coupling factors. There are two forms of coupling in the literature named as local and non-local coupling terms. In the local form of interactions, energy density has a direct relation with it [19], whereas non-local forms are directly proportional to Hubble parameter H and energy density. Yang et al. [20] investigated the observational consequences of stable interacting DE models in the light of recent cosmological data.
In the field of cosmology, the autonomous system is another notable technique to study the evolution of the universe. However, this approach does not provide any proper solution but measures the stability of the system efficiently [21]. Phase space analysis also plays a vital role in the inspection of stability for different models by describing momentum as well as a position at each point of the system. Xiao and Zhu [22] investigated dynamical properties of DE models through barotropic fluid in the background of quantum gravity. Shahalam et al. [23] analyzed the stability as well as accelerated behavior of the universe through the interaction of scalar field models with dark matter in the gravitational background of the FRW universe. We have extended their work to study the impact of nonlinear electrodynamics on the stability of FRW universe model through phase space portrait [24]. Recently, Sharif and Mumtaz [25] studied the phase space analysis of locally rotationally symmetric (LRS) Bianchi type I universe by taking different coupling functions between DE models and dark matter. In the present work, we study the effect of bulk viscosity on the stability as well as the expansion of the cosmos by using the dimensionless quantities.
The insertion of the bulk viscosity coefficient in the matter has become a focus of attention due to its remarkable consequences in astrophysics. This is a continuation of the idea of fluid velocity in the surrounding of dissipation and is accepted as a favorable phenomenon in inflation [26][27][28]. It has also been contemplated that bulk viscous fluid as a source of accelerated expansion can account for dark matter as well as late-time expansion of the universe [29]. The cosmological aspects of viscosity can be observed as a measure of the pressure required to reinstate equilibrium position when the expansion of fluids occurs in an expanding universe. Acquaviva and Beesham [30] analyzed the stability of FRW universe dominated by nonlinear bulk viscous matter. Sasidharan and Mathew [31] discussed the asymptotic properties of the FRW universe for three different bulk viscous coefficients and predicted the current cosmic expansion through their solutions. Sharif and Mumtaz [32] studied the effects of viscosity on the expansion of LRS Bianchi type I universe through phase space analysis.
In this paper, we explore stability as well as expanding behavior of FRW universe, by considering the interaction of DE models with bulk viscous dark matter, through phase space analysis. The paper is organized as follows. In Section 2, we provide a basic formalism for dynamical equations and bulk viscous parameter. We establish an autonomous set of equations through the dimensionless variables to analyze the pattern of cosmic evolution. Section 3 contains the phase space portrait of three specific forms of interaction between phantom energy and bulk viscous fluid. The characteristics of tachyon field with Q = γσ φ is discussed in Section 4. In the final section, we summarize the results.

Dynamical Equations
In this section, we compute dynamical equations for bulk viscous matter and DE in the background of isotropic and homogeneous model of universe given by such that a(t) depicts the scale factor. The viscous energy-momentum tensor [33] has the form where σ, p, ζ and C γ correspond to the energy density, pressure, bulk viscosity and dynamical velocity (canonical momentum) of the fluid, respectively. The dynamical velocity is defined by C γ = ξu γ and ξ is named as specific enthalpy of a fluid element, with ψ its rest mass density. The dynamical velocity has been contemplated as a generalization of the fluid's velocity in the presence of dissipation [34]. The bulk viscosity can be measured by ∇ γ C γ which is compatible with the idea that it vanishes for an incompressible fluid, i.e., ∇ γ C γ = 0. Moreover, the rest mass density remains conserved along the lines of flow, i.e., ∇ γ (ψu γ ) = 0 such that ψ = ψ 0 a −3 , where ψ 0 is an integration constant. Now, we consider the interaction between the phantom field and the dark matter through the coupling Q. This will be helpful to analyze how interaction terms can influence the qualitative behavior of models containing scalar fields and bulk viscous matter. These terms have been studied widely in literature in the context of reheating and inflation [35,36]. The density of the universe remains preserved throughout the evolution, but, for each constituent, it does not hold the conservation law: where H =˙a a , while dot represents the temporal derivative. In addition, p φ , p m correspond to pressure and σ φ and σ m exhibit the energy density of phantom field and matter (dust), respectively with σ tot = σ m + σ φ and p tot = p m + p φ . In the considered framework, the energy between the components can transform in both ways depending upon the signature of the coupling function. For Q > 0, the energy outflows from the matter part towards the phantom field and vice versa. In literature, several interaction terms have been studied to observe their effect on the cosmological issues [37][38][39][40]. The field equations, in a spatially flat FRW universe, take the form where [41]. To transform a set of equations into autonomous system, we introduce the dimensionless quantities as which satisfy the following relations where U = ln a. Through exponential potential, the parameter l becomes a constant value. In this scenario, Equations (4) and (7) can be expressed aṡ The effective equation of state (EoS) and field density parameter [41] are given by

Coupled Phantom Field
In this section, we examine the impact of different interaction terms (between viscous dark matter and phantom field) on the stability of FRW universe.

Coupling Q = γσ m
We take the coupling function Q = γσ m , where γ is an arbitrary constant. For this coupling, Equation (12) leads tö where δ = Hζa 3 is the viscosity factor. The autonomous Equations (9) and (10) reduce to To analyze the stability around a critical point, we evaluate its Eigenvalues (EVs) by converting the autonomous set of equations into a Jacobian matrix. Eigenvalues with negative or positive real part correspond to the sink/stable or source/unstable critical nodes, respectively, whereas the EVs with different signatures correspond to the saddle points.
For the stationary node C 1 = (x, 0), with the respective EVs are not given due to lengthy expressions. To analyze the impact of γ, δ and l on the stability of stationary points, we portray their dynamical characteristics in Figures 1 and 2 by taking different values of the model parameters. When l > 0, C 1 shows stable behavior for δ < 0 along with all values of γ as both of its EVs become negative for this range of parameters. This critical point depicts unstable behavior for other possible values of γ and δ as shown in Figure 1. It is mentioned here that the universe experiences accelerated expansion when the relation (3δ + 1)x 2 > Γ=5, Β=10 For C 2 = (x, 0) with the EVs correspond to the stable node for negative values of δ and l as all the EVs have negative signature.
For positive values of δ and l, this point fluctuates between the unstable and saddle points.
It is noted that, for negative values of δ, C 3 depicts stable results for positive as well as negative values of γ. For l < 0, we get a stable scenario for the positive values of δ.
When C 4 = (x, 0) with the point gives the opposite results to the point C 3 . We achieve more stable results as well as an accelerating era for particular values of parameters as compared to [23]. Most of the stationary points behave as the stable nodes for diverse ranges of parameters as displayed in Figures 1 and 2. It is mentioned that the coincidence problem cannot be resolved through the considered coupling function as density parameter does not attain the proper value, i.e., Ω DE < 0. The evolutionary results are summarized in Table 1 and acceleration phenomenon is discussed for the parameters γ = ±20 and δ = ±1. Table 1. Stability analysis of phantom universe for Q = γσ m .

Coupling Q = γσ φ
In this scenario, Equation (12) reduces tö A system of autonomous Equations (9) and (10) takes the form For the stationary node C 1 = (x, 0) with we get stable solutions for positive as well as negative values of γ and δ under some restrictions. For the positive values of dynamical quantities, stable nodes can be achieved with the condition δ γ while, for negative values, the condition is just reversed. It is noted that, for positive values of parameters, we cannot determine the accelerating phase of the universe as the value of effective EoS becomes complex for γ = 20 and δ = 1. It is observed that this critical point can alleviate the coincidence problem as the density parameter fulfills the criterion, i.e., Ω DE Ω DM = O(1) for γ < 0, δ > 0 and γ > 0, δ < 0.
When C 2 = (x, 0) with the system fluctuates between the unstable and saddle points for negative values of δ, whereas the stable results can be evaluated for δ > 0 corresponding to all values of l and γ. By increasing the parameter δ, we find more negative EVs depicting stable future attractor. For C 3 = (x, 0), where This fixed point shows unstable or saddle trends except for γ, δ > 0 and l > 0 with γ δ. In this case, the universe undergoes accelerated expansion for negative values of γ. When C 4 = (x, 0) with the system shows varying behavior. When l > 0, saddle trends can be observed for several combination of the parameters. The results are shown in Figures 3 and 4 and collective behavior in Table 2. The acceleration phenomenon is discussed for the parameters γ = ±20 and δ = ±1. For some choices of model parameters, the critical points become imaginary which can not be displayed on the phase portrait; therefore, their results are summarized in Table 2.

Ranges of Model Parameters Stability Acceleration
Now, we take the coupling as a linear function ofσ m andσ φ which reduces Equation (12) tö The autonomous system becomes For the critical point C 1 = (x 1 , 0), the value of x 1 is not given due to lengthy-expression. The corresponding EVs show stable solutions for the positive values of δ and l while the accelerating behavior cannot be determined as the effective EoS gives complex values. However, we obtain the stable as well as saddle zones for l < 0 and δ < 0 corresponding to complete range of γ. For γ < 0 and δ < 0, we find the accelerated expansion of the universe with oscillating behavior between unstable and saddle nodes. For C 2 = (x 2 , 0), the critical point varies between the unstable and saddle points for negative values of δ with γ > 0 as well as γ < 0 while the cosmic expansion has been observed in the latter case only. For δ > 0, the system moves towards more stability corresponding to l > 0 and l < 0.
When C 3 = (x 3 , 0), we observe mostly saddle and an unstable trend for the positive vales of l, whereas dynamical stability of C 3 corresponds to l < 0. We get a stable node for δ, l < 0 and γ > 0 which also lies in the accelerating phase. For C 4 = (x 4 , 0), the conditions are totally reversed as that of C 3 . Here, we have stable solutions for l > 0 while the negative values of l lead the universe towards instability. The analysis is summarized in Table 3 and the numerical plots are given in Figures 5 and 6.  Figure 6. Phase plane portrait of phantom universe for Q = γ(σ m +σ φ ) with l < 0. Table 3. Stability analysis of the phantom universe for Q = γ(σ m +σ φ ).

Coupled Tachyon Dynamics
In this section, we consider tachyon field as a DE source and interact it with bulk viscous dark matter. In this scenario, the conservation equations turn out to bė where [42]. The evolution equations become Introducing the dimensionless quantities such that the autonomous set of equations turn out to be For the inverse square potential, the parameter l represents a constant value. We consider the interaction term Q = δσ φ which reduces Equations (26) and (27) tȯ φ The effective EoS is given by For C 1 = (0, 0), the corresponding EVs are It is found that EVs indicate a stable zone only for positive values of γ and δ with the restriction γ δ, while, for the reverse condition, the system moves gradually towards the saddle zone (Figures 7 and 8).
Moreover, the effective EoS becomes negative for δ > 1 9 , which corresponds to the accelerated expanding phase of the cosmos, but still the coincidence problem can not be resolved. It is also noticed that this critical point has no dependence on l and shows saddle and unstable trends for other possible values of parameters. When C 2 = (±1, 0) and C 3 = (±1, ± (−1+3δ) √ 3 l ), the metric is indeterminant. We cannot evaluate the corresponding EVs as the system becomes undefined at these fixed points for interacting tachyon field model. The results are given in Table 4.  Table 4. Stability analysis of tachyon universe for Q = δσ φ .

Conclusions
In this paper, the impact of bulk viscosity on the stability of FRW universe has been examined through different coupling functions between dark matter and DE models. We have constructed an autonomous set of equations by using dimensionless quantities to inspect the stability of the prescribed system. We have evaluated the stationary nodes by settingx =ý = 0 and their respective EVs to explore the effects of γ, δ and l on the stability as well as expansion of the universe.
Firstly, the interaction of dark matter with phantom field has been investigated through EVs for diverse range of parameters (Figures 1-6). We analyze the dynamical behavior of three different coupling functions through phase space analysis. For Q = γσ m , C 1 and C 2 depict a stable DE model for δ < 0 corresponding to l > 0 and l < 0, respectively, followed by late accelerated expansion of the cosmos. The critical points C 3 and C 4 correspond to stable points for δ < 0 and δ > 0, respectively, for all values of γ which completely lie in the DE dominated era as w e f f < − 1 3 . For this DE model, all points lie in the de Sitter cosmic phase if the condition (3δ + 1)x 2 > 1 3(1+9δ) is satisfied. For the second coupling Q = γσ φ , the critical point C 1 shows accelerated stable evolution for γ < 0, δ > 0 and γ > 0, δ < 0 which can also alleviate the coincidence problem.
When Q = γ(σ m +σ φ ), we have found stable solutions of C 1 and C 2 for positive values of δ. Finally, for the coupled tachyon field with Q = δσ φ , the cosmic portrait represents the stable behavior for positive values of γ and δ when γ δ. It is worth mentioning here that these nodes belong to the accelerating phase of the universe if δ > 1 9 . We conclude that the inclusion of a bulk viscosity coefficient in the coupling of dark sectors makes the FRW universe dynamically more stable as compared to [24]. For all the considered couplings, our analysis shows that the coincidence problem can be resolved through the interaction term Q = γσ φ in the coupled phantom field model in contrast to [23]. In [23], for the coupling Q = γσ φ , a stable fixed point has been found which behaves as an attractive node but did not depict the accelerating phase of the universe and hence did not solve the coincidence problem. In our work, for the same coupling, the critical point C 1 can alleviate the coincidence problem as it shows the accelerated stable evolution and density parameter fulfills the criterion, i.e., Ω DE Ω DM = O(1) for γ < 0, δ > 0 and γ > 0, δ < 0. The coupling Q = γσ φ alleviates the coincidence problem which shows the compatibility of our results with the current cosmological observations. For future directions, we can adopt this specific form of coupling for other cosmological issues like tension on the determination of H 0 via low-redshift analysis [43] and σ 8 [44], cosmic microwave background anomalies [45] in the light of recent cosmological data.