title: Cosmic expansion with matter creation and bulk viscosity
authors: Víctor H. Cárdenas Miguel Cruz Samuel Lepe
text: 
  
  
   
  
  victor.cardenas@uv.cl
  
  
  miguelcruz02@uv.mx
  
  
  samuel.lepe@pucv.cl
  
  ^1Instituto de Física y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Valparaíso, Chile 
  
  ^2Facultad de Física, Universidad Veracruzana 91000, Xalapa, Veracruz, México
  
  ^3Instituto de Física, Facultad de Ciencias, Pontificia Universidad Católica de Valparaíso, Avenida Brasil 2950, Valparaíso, Chile
  
  
  
  
  We explore the cosmological implications at effective level of matter creation effects in a dissipative fluid for a FLRW geometry; we also perform a statistical analysis for this kind of model. By considering an inhomogeneous Ansatz for the particle production rate we obtain that for created matter of dark matter type we can have a quintessence scenario or a future singularity known as little rip; in dependence of the value of a constant parameter, η, which characterizes the matter production effects. The dimensionless age of this kind of Universe is computed, showing that this number is greater than the standard cosmology value, this is typical of universes with presence of dark energy. The inclusion of baryonic matter is studied. By implementing the construction of the particle production rate for a dissipative fluid by considering two approaches for the expression of the bulk viscous pressure; we find that in Eckart model we have a big rip singularity leading to a catastrophic matter production and in the truncated version of the Israel-Stewart model such rate remains bounded leading to a quintessence scenario. For a non adiabatic dissipative fluid, we obtain a positive temperature and the cosmic expansion obeys the second law of thermodynamics.      
  
  
  
  
  Cosmic expansion with matter creation and bulk viscosity
      Samuel Lepe^3
      December 30, 2023
  ========================================================
  
  
  
  
  
  
  § INTRODUCTION
   
  
  
  The nature of the late times acceleration of the observable Universe is still far from our actual understanding <cit.>. However, this challenge has not been a reason for the community to remain staidly; on the contrary, this has motivated an exhaustive search for models or scenarios beyond General Relativity that attempt to roughly describe the current stage of the Universe in order to elucidate the nature of the catalyst of such accelerated expansion, usually termed as dark energy. As is well known, the cosmological constant approach is a promising scenario, but it has yet to face its own battles. For instance, the origin of the cosmological constant must be at Planck scales but its effects are observed only at cosmological scales (the current accelerated expansion). This difference between scales where the cosmological constant becomes relevant has turned the problem unmanageable, at Planck scales the value of the cosmological constant is greater by approximately 120 orders of magnitude than expected. A reconciliation for the capricious behavior of the cosmological constant at different scales was proposed in Ref. <cit.>, but this description depends on a quantum formulation for the fluctuations of spacetime; nowadays there is no quantum theory of gravity, therefore the description made was at a semiclassical level, this scheme it lacks of a full characterization for these quantum fluctuations and consequently of their cosmic evolution.
  
  
  Another more recent but no less controversial problem of our current Universe is the so-called H_0 tension, being H_0 the Hubble constant, this value represents the rate of expansion of the Universe at present time. The problem lies in the discrepancy between the value reported for H_0 by the Planck collaboration <cit.> and the one reported in <cit.>. Such difference between both values is not attributable to systematic errors. This has originated for example that schemes where dark energy has the peculiarity of interact with dark matter may be a response to solve this tension among other possible scenarios, that is, we keep looking models beyond standard cosmology, see for instance <cit.>. An interesting review on the interacting scheme for the dark sector can be found in <cit.>. Other scenarios promote that the H_0 tension is the result of a tension in the value T_0 of the cosmic microwave background (CMB) temperature. Despite in some schemes the H_0 tension can be alleviated, this implies a new paradigm regarding the well established value for the temperature of the CMB <cit.>. Other proposals suggest that an enhancement in the geometric description of the Universe is a viable alternative to resolve the aforementioned tension with no need of make use of exotic components or a reformulation of our conception of the gravitational interaction. In Ref. <cit.> was found that the H_0 tension can be mitigated by considering only a non-vanishing torsional tensor to describe the dark matter sector.
  
  
  In this work we consider the bulk viscous effects plus matter production in the cosmological fluid as an alternative to describe the late times behavior of the observable Universe; it is well known that viscosity effects have a negative contribution in the pressure of the cosmological fluid, therefore such effects lead to accelerated cosmic expansion. On the other hand, in some works it is discussed that matter production is due to the expansion of the Universe <cit.>. For this latter scenario the production of particles is consequence only of the dynamical gravitational background, i.e., the created particles interact only with gravity, therefore their abundance is determined by their masses uniquely. Given the features of the component known as dark matter, this theoretical scenario has been the subject of various tests to explain its nature, despite the fact that dark matter only interacts with gravity, there could be various mechanisms for its production <cit.>. As we will see below, the bulk viscous pressure can be written in terms of the matter production rate if we take into account the adiabatic condition for the entropy, i.e., S =. We will focus on the cosmological implications of the model at effective level. In the first part of the work we assume an inhomogeneous matter production rate, in general grounds this term is assumed to be constant or can be given as a function of the Hubble parameter Γ∝ H, but in this work we consider a generalization for this production rate given as Γ = Γ(H,Ḣ). A similar form for this term was considered in Ref. <cit.> by the authors. We obtain that for this kind of model in order to have a consistent description we must include baryonic matter, once it is included, we obtain that a quintessence dark energy behavior is allowed by the model at effective level if we restrict the parameters of the resulting model with the use of recent observational data. In the second part of the work we drop the Ansatz philosophy for the matter production term, Γ, and we consider two approaches for the dissipative effects in the cosmic fluid, in this second part such approaches lead us to a specific form for the matter production term. As we will see later, in one case the description allows a phantom cosmology leading to a catastrophic production of matter. However, for the second approach the model only allows a quintessence dark energy behavior again. An interesting characteristic of a dissipative fluid is the generation of entropy, in the perfect fluid description the entropy production and heat dissipation are lost <cit.>. However, a recent work showed that at cosmological scales the second law of thermodynamics is still valid <cit.>. On the other hand, in the literature can be found that bulk viscous effects are not ruled out at all by the observational data a provide a framework in which the H_0 tension can be weakened <cit.>. As additional examples, see for instance the Ref. <cit.>, where was stated that the presence of dissipative effects in the fluid contribute significantly to reproduce the experimental measurements of the longitudinal polarization of hyperons produced in relativistic heavy-ion collisions. In Ref. <cit.> was found that when the inflationary process is attenuated by dissipative causes, the inflaton interchanges its energy density with a emerging radiation component, which can be associated to the CMB and at late times such model leads to a cosmological constant type evolution.
  
  
  Over the years the approach in which matter production effects are employed has changed slightly. At first it was thought that these effects were important only in the early Universe providing a natural explanation for the reheating phase of the inflationary process <cit.> or depending on the rate of matter production, the origin of the Universe could be free of an initial singularity <cit.>. However, in the meantime it has also been shown that the matter production effects can play an important role in cosmic evolution, in Ref. <cit.>, supported by cosmological observations the authors showed that this scheme leads to an Universe in which the dark energy sector can be emulated by the particle production. Other relevant scenarios in which the production of matter has an important character can be found in <cit.>, this work explores the possibility that the current state of the Universe is transitory and will eventually present a decelerated phase, this transition is possible by considering matter creation in their model. Another example is given in <cit.> where the cosmological model encloses entropic forces with matter creation. Other interesting works can be found in <cit.>.
  
  
  The outline of this work is as follows: in Section <ref> we provide the cosmological equations for a model with matter creation effects and viscosity. Under the adiabatic condition for the entropy we construct the effective parameter state and we explore an inhomogeneous Ansatz for the particle production rate. We obtain the corresponding Hubble parameter and we calculate the age of this kind of Universe and we consider the inclusion of baryons. In Section <ref> we perform the statistical analysis of the model with baryonic matter and the use of current observational data. In Section <ref> we proceed in the inverse order, with the consideration of the dissipative and matter creation effects we construct the form of the particle production rate, Γ, we do not use the Ansatz philosophy. We explore the well known Eckart model and the truncated version of the causal Israel-Stewart description, where a cosmological constant evolution can be emulated. In Section <ref> we provide a description of the cosmological model when the entropy production is not adiabatic, we find that the temperature of the fluid is positive and obeys the second law of thermodynamics. Under this approach the cosmic expansion does not have the problem of negative entropy or temperature. Finally, in Section <ref> we give the final comments of our work. We will consider 8π G=c=k_B=1 units throughout this work.  
  
  
  
  
  § MATTER CREATION AND BULK VISCOSITY
  
  
  
  
  For a dissipative fluid the local equilibrium scalars such as the particle number density and its energy density are not altered by the dissipative effects. However, the pressure deviates from the local equilibrium pressure
  
      p_eff = p + Π,
  
  where Π is the bulk viscous pressure. From now on the quantity, p+Π, will denote effective pressure for the dissipative fluid. In this description the Friedmann equations for a flat FLRW spacetime can be written as follows
  
      3H^2 = ρ,    Ḣ+H^2 = -1/6[(1+3ω)ρ + 3Π],
  
  being ρ the energy density of the dissipative fluid, H denotes the Hubble parameter and the dot stands for derivatives with respect to time. In the previous equations we have considered a barotropic equation of state between the density and the pressure given as, p = ωρ, where ω is commonly known as parameter state and it is constrained to the interval [0,1). Note that the consideration of the Friedmann equations given in (<ref>) leads to a non conservation equation for the energy density
  
      ρ̇+3H(1+ω)ρ = -3HΠ.
  
  On the other hand, if matter creation exists, i.e., gravitational particle production, then the continuity equation for the particle number density takes the form 
  
      ṅ + 3H n = n Γ,
  
  where the possibilities Γ >0, Γ < 0 denote source or sink of particles, respectively. Γ it is also known as the particle production rate. From the Gibbs equation <cit.>
  
      T dS = d (ρ/n) + p d (1/n),
  
  we can write
  
      nTṠ   =   ρ̇ - ρ(1+ω)ṅ/n, 
         =    - 3HΠ - ρΓ (1+ω),
  
  where the Eqs. (<ref>) and (<ref>) were considered. Note that for this approach the entropy is not longer a constant. However, if we assume Ṡ=0, i.e, an adiabatic dissipative fluid with particle creation in order to be in agreement with the standard cosmological model, we obtain the following condition
  
      Π = - ρ(1+ω)Γ/3H,
  
  we will have a negative contribution from the viscous pressure to the non equilibrium pressure if the particle production rate is positive, i.e., no annihilation. If we insert the previous equation in (<ref>) we can write the continuity equation for the density in its standard form, ρ̇ + 3H(1+ω_eff)ρ = 0, where the effective parameter state has the form
  
      ω_eff = ω - (1+ω)Γ/3H.
  
  It is worthy to mention that if the dissipative fluid behaves as standard dark matter we have, ω = 0 (p=0), then the effective parameter only depends of the particle production effects, this case was considered in Ref. <cit.>, in such case the condition Γ > 3H must be fulfilled in order to have an effective phantom behavior.  
  
  
  
   §.§ Ansatz for the particle production rate
  
  
  The most simple assumption for Γ is given by a constant production rate. However, in order to study the implications of matter production on the cosmic expansion, some specific functions of the Hubble parameter for the particle production rate have been studied, for instance Γ∝ H^α, being α an appropriate constant <cit.>. A more general form for Γ can be found in Ref. <cit.>, where
  
      Γ = Γ(ρ, p(ρ), H, Ḣ, Ḧ,...),
  
  and the function p(ρ) is viable by means of the equation state. The form given in Eq. (<ref>) for Γ it is known as inhomogeneous particle production rate. In our case we will consider Γ given as
  
      Γ = Γ(H,Ḣ) = γ'(H)Ḣ = dγ(H)/dt,
  
  where the prime denotes derivatives with respect to the Hubble parameter by means of the chain rule. Note that the form of Γ given above together with Eq. (<ref>) lead us to a direct integration integration of Eq. (<ref>), yielding
  
      ρ = ρ_0 a^-3(1+ω)exp[(1+ω)γ(H)],
  
  then if γ(H) ∝ln (H^δ) and δ being a constant, we could write consistently an expression for the Hubble parameter by means of the first Friedmann equation. We focus on the following expression
  
      γ(H) = 2(1-1/η)ln(H/H_0),
  
  being η a constant value and H_0 the value of the Hubble parameter at present time. Using the above γ function and (<ref>), we can write
  
      Γ/3H = -2/3(1-1/η)(1+q),
  
  where q is the deceleration parameter defined as 1+q := -Ḣ/H^2. As can be seen, the particle production rate can be written as a function of the deceleration parameter. On the other hand, using the acceleration equation (<ref>) together with the expression for the viscous pressure obtained from the adiabatic condition (<ref>), one gets
  
      Γ/3H = 1 -2/3(1+q)/(1+ω),
  
  by equating both expressions for the quotient Γ/3H, we arrive to the following result
  
      1+q = 3/2{1/1+ω - (1-1/η)}^-1,
  
  then, the deceleration parameter takes a constant value if matter creation effects are introduced in a dissipative fluid. Note that if in the previous expression we consider the value, η = 2(1+ω)/(3+5ω), we have q=0, which represents a Dirac-Milne universe. This kind of Universe expands at constant rate since ä = 0. For standard dark matter we have, q = - (1-3η/2), in this case the accelerated cosmic expansion will take place only if the condition η < 2/3 is satisfied and no accelerated expansion for η = 2/3. Using the Eq. (<ref>) we can write for the effective parameter (<ref>)
  
      ω_eff = -1 + 2/3(1+q) = -1 + {1/1+ω - (1-1/η)}^-1.
  
  As expected, the effective parameter state has a contribution from the matter creation effects. If the following conditions are satisfied η > (1+ω)/ω or 0 < η < 2(1+ω)/(3+5ω), the effective parameter (<ref>) will behave as a phantom or quintessence fluid, therefore accelerated cosmic expansion can be obtained for a dissipative fluid with matter creation effects in presence of ordinary matter. On the other hand, for ω = 0 we have ω_eff = -1 + η, thus from the condition 0 < η < 2/3, the dissipative fluid will behave as quintessence and will exhibit phantom behavior with η < 0. For a Dirac-Milne Universe we will have ω_eff = -1/3.
   
  Using the Friedmann constraint (<ref>) we have, 3H^2 = ρ = ρ_0a^-3(1+ω)exp(γ(H)), rearranging properly the terms appearing in this expression after inserting the Eq. (<ref>); we obtain for normalized Hubble parameter
  
      E(z) = Ω_ρ^1/2(1-Δ)(1+z)^3(1+ω)/2(1-Δ),
  
  where we have defined Δ := (1+ω)(1 - 1/η), Ω_ρ corresponds to the standard definition of the density parameter, Ω_ρ := ρ_0/3H^2_0 and we also used the relation between the scale factor and the redshift, 1+z = a^-1, besides E(z) := H(z)/H_0. Note that for η = 1 we have a null contribution from matter creation effects and we recover the standard cosmology. In the standard dark matter case we have
  
      E(z) = Ω_ρ^1/2(1-Δ)(1+z)^3/2(1-Δ),
  
  for the quintessence scenario we have 0 < η < 2/3, therefore - ∞ < Δ < - 1/2, then as we approach to the far future, we have for the normalized Hubble parameter, E(z→ -1) → 0. For the phantom regime we have η < 0, thus Δ > 1, which leads to a singular scenario as we approach to the far future, E(z→ -1) →∞. A similar behavior is obtained for the energy density as we approach to the far future, for quintessence ρ→ 0 and ρ→∞ for phantom. It is worthy to mention that in the phantom scenario the singular nature takes place only at the far future and not for a specific value of the redshift, therefore this kind of singularity corresponds to a little rip. From the current observational data it is not possible to determine if the final fate of the Universe is a future singularity or not but if consistency with the supernova is demanded, then this kind of singular model could represent a viable alternative to the ΛCDM model <cit.>. If we consider the Eq. (<ref>) and H(a) = dln a/dt, we can compute the dimensionless age of the Universe
  
      H_0t(a) = ∫^a_0da'/a'H(a'),
  
  as stated in Ref. <cit.>, the value for the dimensionless age at present time, H_0t_0, is around 1 for the ΛCDM model and this is independent of the transition from the decelerated to accelerated stage, the authors call this fact as synchronicity problem given that it appears we are living in a special time. In Fig. (<ref>) we show the behavior of the dimensionless age of the Universe if we consider the Eq. (<ref>) with a quintessence behavior. We also consider the value ω = 1/3 and we give the value 0.315 ± 0.007 <cit.> to the energy density parameter Ω_ρ, according to the latest Planck collaboration results, this case corresponds to the dashed line in the figure. On the other hand, if we consider that created matter is the only component of the Universe we must have, Ω_ρ = 1, this corresponds to the solid line of the plot, this latter case is closer to the value 1 at present time (a=1). Note that at present time the value of the dimensionless age of the Universe is increased in both cases, however, the obtained values are close to 1; the augmentation in the age of the Universe at present time by dark energy was discussed in Ref. <cit.> for phantom dark energy.  
  
  
  
  
  
  
  
  
   §.§ Including baryonic matter
  
  
  
  It is worthy to mention that if we consider only the matter creation effects we are left with a constant deceleration parameter in our description, this can be seen in Eq. (<ref>). Despite that dark energy behavior can be obtained for certain values of the parameters η and ω, the model itself is not consistent with the ΛCDM model and the type Ia supernovae data; this issue can be alleviated with the inclusion of baryons as done in Ref. <cit.>. In this new scenario the form of the normalized Hubble parameter can be expressed as follows
  
      (H/H_0)^2 = Ω_B(1+z)^3 + (1-Ω_B)(1+z)^3(1+ω)(H/H_0)^2Δ,
  
  where we take into account the normalization condition, Ω_B + Ω_CM = 1, i.e., for the late times description of the Universe we only consider the contribution from the gravitationally created matter and baryons. From the previous expression it is not possible to find the normalized Hubble parameter analytically. However, if we consider the acceleration equation (<ref>) with the assumption, p_B = 0, and the equations (<ref>), (<ref>) together with the standard relationship between the redshift and the scale factor given above; we can write the following differential equation for the normalized Hubble parameter
  
  
      (1+z)EdE/dz-E^2-1/3(1+ω)(1-Ω_B)(1+z)(1-1/η)d ln E/dz = 1/2[(1+3ω)(1-Ω_B)+Ω_B(1+z)^3],
  
  
  from which E(z) can be obtained by numerical integration. We will use this equation to constraint the parameters of the model with the use of current cosmological observations later. In this we will have at effective level from the contribution of baryons and matter creation effects
  
      ω_eff = η - 1/1+Ω_B/3H^2_0(1-Ω_B)^η(1+z)^3(1-η),
  
  where we have considered the equation (<ref>) for the pressure together with the assumption that created matter behaves as standard dark matter, ω=0, and p_B=0.
  
  
  
  
  § OBSERVATIONAL CONSTRAINTS
  
  
  
  
  In this section we study to what extend the evolution implied by Eq. (<ref>) describes appropriately the observations. In particular, we use the latest type Ia supernovae sample called Pantheon <cit.> consisting in 1048 data points. The data gives us the apparent magnitude m at maximum from which we can compute the distance modulus μ = m-M with M the absolute magnitude for type Ia supernovae. Here we compute the residuals μ - μ_th and minimize the quantity
  
      χ^2 = (μ - μ_th)^T C^-1(μ - μ_th),
  
  where μ_th = 5 log_10( d_L(z)/10pc) gives the theoretical distance modulus, d_L(z) is the luminosity distance given by
  
      d_L(z)= (1+z)C/H_0∫dz/E(z),
  
  C is the covariance matrix released in <cit.>, and the observational distance modulus takes the form
  
      μ = m - M + α_1 X - α_2 Y,
  
  where m is the maximum apparent magnitude in band B, X is related to the widening of the light curves, and Y corrects the color. usually, the cosmology – specified here by μ_th – is constrained along with the parameters M, α_1 and α_2. These nuisance parameters are then marginalized to obtain the posterior probabilities for our parameters of interest: w, Ω_B and η. In the statistical study we use a prior for Ω_B and we assume w=0 assuming the contribution be as dark matter. Then, the free parameters to fit are Ω_B and η. Using the value Ω_B=0.0493 ± 0.0006 based on <cit.>. In this case it is possible to obtain as best fit values Ω_B = 0.049 ± 0.001, and η=0.174 ± 0.015. The confidence contours for 1σ and 2σ are shown in Fig. (<ref>). If we insert the constrained values for Ω_B and η in Eq. (<ref>) together with the Hubble constant reported in Ref. <cit.> one gets that at present time the effective parameter lies in the interval [-0.840997,-0.810997]; this interval is within the quintessence region determined by the DES collaboration (Dark Energy Survey) for the parameter state at present time, ω _de,0 = -0.95^+0.33_-0.39 <cit.>. On the other hand, as we approach the far future, from Eq. (<ref>) one gets ω_eff(z → -1) →η -1 = -0.826 ± 0.015, thus the value of the parameter state at the far future is very close to its present time value. Therefore this model represents a quintessence dark energy behavior along the cosmic evolution and is in agreement with some sets of observational data such as the latest Planck results (ω_de,0= -1.028 ± 0.031) and DES collaboration, where the quintessence behavior for dark energy is allowed.   
  
  
  
  
  
  
  § CONSTRUCTING Γ FROM VISCOUS MODELS
  
  
  
  
  In this section we construct the particle production rate from bulk viscous considerations. We will not consider a specific Ansatz for this term. According to the Israel-Stewart formalism, the bulk viscous pressure must obey the following transport equation <cit.>
  
      Π = -3ξ(ρ)H - τΠ̇ - ζ/2τΠ[3H + τ̇/τ - ξ̇/ξ - Ṫ/T],
  
  where ξ(ρ) is the bulk viscosity coefficient and given that it is a function of the energy density we must have ξ(ρ) ≥ 0 a typical form for the bulk viscosity coefficient can be found in the literature as ξ∝ρ^s being s a constant; T is the barotropic temperature which is also a function of the energy density by means of the integrability Gibbs condition. In general, an interesting parametrization for the bulk viscosity coefficient is given in terms of the scalar expansion of the fluid, i.e, ξ∝Θ where Θ = H, thus ξ∝ρ^1/2, i.e., s=1/2; physically this represents an increasing viscosity in the case of increasingly movements in the fluid <cit.>[The introduction of a 4-velocity field u^α provides the 1+3 decomposition of spacetime. The metric tensor splits as g_αβ = h_αβ - u_αu_β, where h_αβ is the projection tensor and is orthogonal to u^α. The kinematics of the 4-velocity field takes the form <cit.>
  
      ∇_αu_β = 1/3Θ h_αβ + σ_αβ + ω_αβ - u̇_αu_β.
  
  where Θ = H generalizes the Newtonian expansion and σ_αβ and ω_αβ generalize the Newtonian shear and Newtonian vorticity, respectively. u̇_α represents a 4-acceleration (the dot represents a time derivative). Spatial homogeneity and isotropy demand u̇_α = σ_αβ =  ω_αβ = 0. The scalar expansion determines the volume evolution.]. On the other hand, τ represents the relaxation time for bulk viscous effects. In Ref. <cit.> was shown with a perturbative method that the Israel-Stewart theory fulfills the causality and stability conditions. The energy functional can be written as ∑_A=1^8Ω_A(δ Z_A)^2 where δ Z_A represents certain combination of the perturbation functions, then the positivity of this functional is guaranteed for Ω_A≥ 0. A specific term is given as
  
      Ω_3(λ) := (ρ + p){ 1- λ^2[(∂ p/∂ρ)_S + ξ/τ(ρ + p)]}≥ 0,
  
  this requirement was shown to hold for all λ where, 0 ≤λ≤ 1, by taking the case λ = 1 one gets
  
      (∂ p/∂ρ)_S + ξ/τ(ρ + p)≤ 1,
  
  note that the first term corresponds to the adiabatic speed of sound, usually denoted as c^2_s, and the second term is identified as the speed of the bulk viscous perturbations, c^2_b, therefore from the above condition we have v^2 = c^2_s + c^2_b≤ 1 and this latter expression guarantees the causality of Israel-Stewart model. For a barotropic fluid, p=ωρ, and using the Eq. (<ref>) we obtain 
  
      c^2_b :=  ξ/τρ(1 + ω)≤ 1 - ω,
  
  usually the r.h.s. of the previous equation is written as ϵ(1-ω), where ϵ is a constant parameter. Therefore, from the previous equation we write for the relaxation time τ = ξ/c^2_bρ(1 + ω) = ξ/ϵ(1-ω^2). From the Friedmann equations (<ref>), we can write for the bulk viscous pressure
  
      Π = -2Ḣ - 3(1 + ω)H^2.
  
  In order to illustrate some results, in the following sections we will study two cases separately.
  
  
  
   §.§ Eckart model
  
  
  If we consider τ = 0 in Eq. (<ref>) we obtain the Eckart model. This approach have been studied exhaustively despite its superluminal propagation problem. However, represents a manageable framework for viscous effects, see for instance the references contained in <cit.>, where the Eckart framework is studied in the late and early Universe and tested with recent observational data. From the Eq. (<ref>) with τ = 0, the Friedmann equation (<ref>) and the standard definition for the bulk viscous coefficient we can write 
  
      Π = -3ξ(ρ)H = -3^s+1ξ_0H^2s+1,
  
  if we consider the previous expression in equation (<ref>) we arrive to an explicit form for the particle production rate given as
  
      Γ(H)/3H = 3^sξ_0/(1+ω)H^2(s-1/2),
  
  note that in order to maintain gravitational particle production, Γ > 0, we must have ξ_0 > 0 for null viscous effects we have not particle production. The case s=1/2 simplifies the above equation leading to a constant rate for the production of particles and it is given as Γ = √(3)ξ_0/(1+ω) or Γ = √(3)ξ_0 for the standard dark matter case. Now, using the Eq. (<ref>) and (<ref>) we can write
  
      Γ/3H = 1+2/3(1+ω)(Ḣ/H^2),
  
  and by means of (<ref>) the previous expression takes the following form
  
      Ḣ/H^2 = -3/2{ 1 + ω - 3^sξ_0H^2(s-1/2)}.
  
  If we consider the standard dark matter case, ω = 0, together with the case s=1/2, then the Eq. (<ref>) can be integrated straightforwardly, yielding
  
      H(t) = H_0[1 + 3√(3)/2(1/√(3) - ξ_0)H_0(t - t_0) ]^-1,
  
  where H_0 denotes evaluation of H(t) at t = t_0. The expansion of the Universe is guaranteed for 0 < ξ_0 < 1/√(3). However, for ξ_0 > 1/√(3) the Hubble parameter becomes negative, therefore we can write it in the following convenient form
  
      H(t) = 2/3√(3)(1/ξ_0-1/√(3))(t_s-t)^-1,
  
  where we have defined
  
      t_s = t_0 + 1/H_0[2/3√(3)(ξ_0-1/√(3))].
  
  Notice that t_s is a constant value and will represent a value for the time at the future at which the Hubble parameter becomes singular. According to the classification for future singularities given in <cit.>, we have a big rip singularity at t = t_s. Thus inserting the Hubble parameter (<ref>) in Eq. (<ref>) with s=1/2 and ω = 0, one gets
  
      Γ(t) = 3√(3)ξ_0H = (2ξ_0/ξ_0-1/√(3))(t_s-t)^-1.
  
  As can be seen, the particle production rate becomes singular as we approach to t_s. This result is contradictory since in the phantom scenario as the Universe reaches the future singularity, matter or any structure must be disaggregated. However, we must keep in mind that have been shown that phantom scenarios are unstable, i.e., at quantum level they have an unbounded negative energy that leads to the absence of a stable vacuum state. On the other hand, in order to avoid this kind of problem in Ref. <cit.> was proposed the interaction of phantom particles and standard matter, at least gravitationally. In consequence, the gravitational interaction allows processes as the spontaneous production from the vacuum of a pair of phantom particles and a pair of photons, to mention some. In this case the phase space integral is divergent, this indicates a catastrophic instability. This instability can be avoided by imposing a non Lorentz invariant momentum space cutoff, but despite this correction the density number of photons and phantom particles can be written as, n ∝ (Γ×), and according to the diffuse gamma ray background observations the production rate of photons, Γ, leads to higher values from the typical ones for the energy of the produced photons. Then, any cosmic ray experiment at earth should be detecting more events than normal and more energetic than those detected until now. These bounds for Γ and n together with other considerations of Ref. <cit.> suggest that the origin of phantom must come from new physics beyond the standard model of particles. See also Refs. <cit.> for similar discussions on the topic. 
  
  
  
   §.§ Israel-Stewart approach
  
  
  Setting ζ = 0 in the transport equation (<ref>) leads to the truncated version of the Israel - Stewart theory <cit.>. This effective model has been widely studied at cosmological level in several contexts, see the references in <cit.>. Therefore, if we consider the truncated transport equation (<ref>) together with the Friedmann equations (<ref>) and the continuity equation for the energy density (<ref>), we obtain a second order differential equation for the Hubble parameter
  
      Ḧ+ 3HḢ(1+ω)+9/2ϵ(1-ω^2){(1+ω)/3^sξ_0H^1-2s - 1 } H^3
          + ϵ (1-ω^2)/3^s-1ξ_0ḢH^2(1-s)=0.
  
  Taking the value s=1/2, as discussed in Ref. <cit.>, if dissipative effects are included in the ΛCDM model, the case s=1/2 is the only value that allows to write the standard form of the well known de Sitter solutions, H_0 = ±√(Λ/3), in the limit ξ→ 0; and if we also consider the change of variable x = ln (a/a_0), we can solve for the Hubble parameter in the standard dark matter case
  
      H(z) = H_0√(cos(β[c_2 + ln(1+z) ] )/cos (β c_2))(1+z)^α,
  
  where c_2 is an integration constant given as 
  
      c_2 = 1/βarctan(2{α - (1+q_0)}/β).
  
  The form written in Eq. (<ref>) for the Hubble parameter is obtained from the condition H(z=0) and the first derivative H'(z=0) given that we are solving a second order differential equation, for simplicity in the notation we have defined α = 3/4 + √(3)ϵ/4ξ_0 together with β = 1/2√(6ϵ/ξ_0(√(3)-6ξ_0-3ξ_0/2ϵ - ϵ/2ξ_0)), which are also constants. On the other hand, q_0 is the deceleration parameter evaluated at z=0. Note that at present time z=0, the Hubble parameter (<ref>) takes the value H_0.
  
  
  In Fig. (<ref>) we depict the normalized Hubble parameter, H(z)/H_0. In Ref. <cit.> the velocity of bulk viscous perturbations was constrained to the interval, 10^-11≪ c^2_b(ϵ) ≲ 10^-8, in order to obtain a similar behavior as the ΛCDM model at perturbative level, therefore we will consider these values for the parameter ϵ and the interval [10^-6,10^-2] for the constant ξ_0. As commented before, the solution (<ref>) depends on the value of the deceleration parameter at present time, we will consider its definition coming from the ΛCDM model q(z) = -1 + 3(2[1+(Ω_Λ,0/Ω_m,0)(1+z)^-3])^-1 together with the normalization condition Ω_Λ,0 + Ω_m,0 =1 and the recent value reported by the Planck collaboration for Ω_m,0 <cit.>. With these considerations we obtain -0.538 ≤ q_0≤ -0.517. In the upper panel of Figure (<ref>) we perform a comparison with the ΛCDM model (shaded region in blue) and the dashed lines corresponds to our viscous model with ξ_0 = 10^-2, ϵ = 10^-11, 10^-8 and the aforementioned values for q_0, as can be seen the model coincides with the ΛCDM model from the present time (z=0) to the far future (z=-1) and also in the past (z>0). In the lower panel we show the behavior of the normalized Hubble parameter when we vary the value of the constant ξ_0 to 10^-6, in this case we can see that the normalized Hubble parameter tends to zero (dotted lines). As the value of the parameter ϵ grows the model deviates significantly from the ΛCDM model, but, this could represent an unstable scenario for the viscous model.      
  
  
  
  
  Then, if we consider the Eq. (<ref>) with ω = 0, we can obtain for the particle production rate
  
      Γ(z)    =    3H(z)-2(1+z)dH(z)/dz, 
         =    H_0(1+z)^α√(cos(β[c_2 + ln(1+z) ] )/cos(β c_2))[3-2α. 
         +   . βtan(β[c_2 + ln(1+z) ] )].
  
  In this case we have that if at some stage of the cosmic evolution the argument β(c_2 + ln(1+z)) of the tangent function takes the value π(1+2n)/2 with n=0,± 1, ± 2,..., the particle production rate diverges. From the previous equation we can write the quotient
  
      Γ/3H = 1-2/3α + β/3tan(β[c_2 + ln(1+z) ]).
  
  We show the behavior of the above expression in Fig. (<ref>) with ξ_0 = 10^-2, ϵ = 10^-11, 10^-8 and q_0 = -0.538, -0.517, each case is represented by the dashed and solid lines. The behavior of this quotient is related to the effective parameter as given in Eq. (<ref>). For ω = 0 the quotient remains positive and tends to 1 at the far future (z=-1). This means that along the cosmic evolution we will always have dark matter production (Γ > 0) and as shown in the plot, the matter production does not become catastrophic. Relating the obtained behavior for the quotient with ω_eff given in (<ref>), we can see that the cosmological viscous fluid with matter production effects from the present time (z=0) to the future will behave as quintessence dark energy and at the far future will emulate a cosmological constant behavior, at this stage the matter production stops. If we consider the value ξ_0 = 10^-6 the resulting behavior for Γ/3H is similar as the one shown in Fig. (<ref>). Therefore, the truncated Israel-Stewart model with matter production does not allow a phantom scenario.  
  
  
  
  
  
  To end this section, we discuss the following Ansatz for the Hubble parameter
  
      H(t > t_s) = |A|/t_s(t/t_s-1 )^-1,
  
  this Ansatz was proposed in Ref. <cit.> for the full Israel-Stewart model in order to study some of its thermodynamics properties. One interesting feature of this Ansatz is that the cosmic evolution starts from an initial singularity given at, t = t_s, that posses the characteristics of a big bang and represents an expanding universe since |A| > 0. Inserting the Ansatz (<ref>) in the truncated differential equation of the Israel-Stewart model (<ref>), we obtain a quadratic equation for |A|, the solutions will be given as
  
      |A|_± = (3+3√(3)ϵ/ξ_0) ±√(9 + 3ϵ^2/ξ^2_0-6√(3)ϵ/ξ_0+36ϵ)/3√(3)ϵ/ξ_0-9ϵ,
  
  then |A| is a constant given in terms of the parameters ϵ and ξ_0. Thus, using the equation (<ref>) with the Ansatz (<ref>), we obtain the following expression for the particle production rate
  
      Γ(t) = 1/t_s(t/t_s-1 )^-1{ 3|A|_±-2},
  
  always that the condition 3|A|_± > 2 it is fulfilled, we will have that the cosmic evolution starts from an initial singularity and will be driven by a dissipative fluid with infinite matter production at the beginning (Γ > 0) and such production decays as the Universe expands. In this case we have the following simple expression for the quotient between the particle production rate and the Hubble parameter
  
      Γ/3H = 1-2/3|A|_±.
  
  Adopting the values for the parameters ξ_0 and ϵ as in the previous case, we obtain some different cases for (<ref>). With ξ_0 = [10^-6,10^-2] and ϵ = [10^-11, 10^-8] we have that Γ/3H ≈ 1 if we consider the solution |A|_+, then by means of Eq. (<ref>) with ω = 0, we find that the dissipative fluid with matter production behaves as a cosmological constant along the cosmic evolution. On the other hand, for |A|_- with the same set of values for ξ_0 and ϵ we have that depending on the election of these values, the quotient can be positive but ≪ 1, negative or zero. Using the Eq. (<ref>) with ω = 0 we can see that for a positive quotient the effective parameter is excluded from the quintessence region and for negative quotient ω_eff turns positive leading to a decelerated expansion. For the case of null particle production we have q=1/2 by means of Eq. (<ref>) and this result also represents a decelerated cosmic expansion. Therefore, the interesting case is given by the solution |A|_+, but as in the case discussed previously, no phantom regime is allowed. 
  
  
  
  
  § NON ADIABATIC EXPANSION
  
  
  
  
  The results obtained in the previous sections emerge from the adiabatic condition for the entropy, i.e. Ṡ = 0, leading to a constant entropy in time. However, if we consider non adiabaticity the thermodynamics description of the cosmological model becomes more consistent <cit.>. From Eqs. (<ref>), (<ref>) and (<ref>) we can compute the particle number density, yielding
  
      n(H) = n_0/V(H/H_0)^2(1-1/η),
  
  where V is the Hubble volume given as (a/a_0)^3. Note that for H=H_0, we simply have the constant density number, n_0/V, and for an expanding Universe this density is always positive. With the inclusion of dissipative effects we must consider the following evolution equation for the temperature <cit.>
  
      Ṫ/T = -3H(∂ p/∂ρ) + nṠ(∂ T/∂ρ),
  
  the integrability condition
  
      ∂^2S/∂ T ∂ n = ∂^2S/∂ n ∂ T,
  
  still holds, therefore we can write
  
      n∂ T/∂ n+(ρ +p)∂ T/∂ρ = T∂ p/∂ρ.
  
  If we consider these results together with Eqs. (<ref>), (<ref>) and a barotropic equation of state in the expression (<ref>), we obtain
  
  
      Ṫ/T = -3Hω1 + 3Π/ρ(1+ω)+Γ/3H/1 + Π/ρ(1+ω)( 1-Γ/3H)+Γ/3H( 1-Γ/3H),_Ξ(t)
  
  
  by considering the relationship between the redshift and the scale factor the equation for the temperature can be written as
  
      1/TdT/dz = 3ω(1+z)^-1Ξ(z),
  
  if we integrate
  
      T(z) = T_0exp[3ω∫^z_0Ξ(z')/1+z'dz'].
  
  It is worthy to mention that the obtained expression for the temperature is definite positive, in the case where we have null contribution from dissipative and matter creation effects, Π = Γ = 0, the temperature takes the standard definition obtained in a single fluid description, T(z) = T_0(1+z)^3ω and for the standard dark matter case the temperature takes a constant value, as in the ΛCDM model.
   
  
  On the other hand, if we consider the expression (<ref>), we can write
  
      Ṡ = 9H^3/nT[-2/3(1+q)+(1+ω)(1-Γ/3H)],
  
  and the expression (<ref>) was considered for the bulk viscous pressure. Now, inserting Eq. (<ref>) in the previous result one gets
  
      Ṡ = 9H^3/nT{2/3(1+q)[(1+ω)(1-1/η) - 1] + (1+ω)}.
  
  As can be seen the entropy production has contributions from matter creation and dissipative effects. For the standard dark matter case we have
  
      Ṡ ={[     9H^3/nT_0{ 1-2/3(1+q)/η},  ,; 9H^3/nT_0{ 1+2/3|1+q|/η},      , ].
  
  where we have considered that at effective level both cases can appear. In order to be in agreement with second law of thermodynamics, Ṡ > 0, the conditions T_0 > 0 and n > 0, must be satisfied together with the following cases: (i) in the quintessence scenario the conditions 1 > 2(1+q)/3η and η > 0 must be satisfied, the positivity of the entropy production is guaranteed for η < 0, (ii) for the phantom regime the fulfillment of the second law is achieved with η > 0 and for η < 0 we must have 1 > 2|1+q|/3η. Thus, under the non adiabatic condition for the entropy the cosmic expansion is free of the negative entropy or temperature problem <cit.>.   
  
  
  
  
  § FINAL REMARKS
  
  
  
  
  In this work we studied the cosmic evolution that emerges from the consideration of matter creation effects in a viscous fluid under the adiabatic condition for the entropy, i.e., constant entropy. As first approach we adopted an inhomogeneous Ansatz for the particle production rate, Γ. In this first scheme we obtained that the model at effective level could describe phantom or quintessence regimes depending on the values of the parameters η and ω, being the case of interest the dark matter type behavior for created matter denoted by ω = 0. For this special case phantom and quintessence behaviors are still present. On the other hand, we computed the dimensionless age of this kind of Universe for two cases: created matter represents the whole content of Universe and secondly created matter represents the actual matter of the Universe; for both cases we got that the age for these kind of universes deviates from 1 but remains close to this value, this is characteristic of universes with presence of dark energy.  However, when the deceleration parameter is computed in this first description we obtain a constant parameter as result, therefore a model of this kind it is not consistent with the cosmological observations or the ΛCDM model. In order to fix this issue we considered the inclusion of baryonic matter and performed the statistical analysis for the model in order to constraint the parameters Ω_B and η with ω = 0. Therefore, with the constrained parameters we proved that this model behaves as quintessence dark energy along the cosmic evolution.
  
  
  As second approach we dropped the Ansatz philosophy for the particle production rate and we constructed such term with the consideration of two descriptions for the bulk viscous pressure: the Eckart model and the truncated version of the Israel-Stewart formalism. The particle production rate has a relevant role in the cosmic evolution since
  
      ω_eff = - Γ/3H,
  
  for ω = 0, as discussed previously. If Γ > 3H, the phantom scenario is allowed and the case Γ < 3H could represent quintessence or decelerated expansion. In the Eckart description a big rip singularity appears, then in this case the Γ term also exhibits a singular behavior leading to a catastrophic matter production. This scenario contradicts the expected behavior of the phantom regime where matter or any structure must be diluted. We infer that this conduct is due to the fact that the Eckart formalism is non-causal, therefore it is not appropriate to describe late times in any type of Universe. On the other hand, if the truncated version of the Israel-Stewart model is considered we have several possibilities to study, however, in this work we focused only on two cases: an analytical solution for the Hubble parameter emerging from this approach and an Ansatz for the Hubble parameter that leads to a cosmic evolution with an initial singularity with the properties of a Big Bang. It is worthwhile to mention that in this scenario the catastrophic matter production is not present. For the analytical solution of the Hubble parameter we obtained that the corresponding Γ term leads to a quintessence dark energy evolution and at the far future the model could mimic the ΛCDM model. These scenarios are possible if the parameter ϵ (a constant parameter responsible of characterizing the velocity of bulk viscous perturbations) lies in a specific interval such that at perturbative level the viscous model is close to ΛCDM, if we consider other values for ϵ we could have other possibilities for the cosmic evolution of the model but our description could reveal an unstable behavior. When the Ansatz for the Hubble parameter is considered, we have two branches characterized by | A|_±; for the branch | A|_+ the cosmological constant expansion is emulated by our model and if we consider the branch | A|_- decelerated expansion is obtained. Then, this latter case is not favored and must be discarded.
  
  
  We must emphasize that it seems that the scenario of dark matter type production plus other cosmological effects does not allow the crossing to the phantom regime, i.e, we only have access to a quintessence dark energy scheme. The case of dark matter type production was discussed by the authors in Ref. <cit.>. We would like to comment that according to some recent tentative results the quintessence scenario seems to be the elected candidate by the observable Universe to lead its late times behavior, in Ref. <cit.> it is reported that the cosmic expansion could be slower than we think since the catalyst of such expansion is a quintessence component. Thus, the kind of model as studied in this work could provide a viable theoretical framework to describe our Universe.
  
  
  Finally, we discuss the cosmic expansion of our model if we implement the non adiabaticity condition. In this description we obtain that the temperature of the viscous fluid is positive definite and in the case of dark matter type behavior, ω=0, we got that such temperature takes a constant value as in the ΛCDM model even if Π and Γ are different from zero. For Π = Γ = 0, we recover the single fluid description at thermodynamics level. Then, the model obeys the second law if we take into account some conditions for the parameters of the model. These phantom or quintessence scenarios have positive entropy and temperature.
  
  
  
  § ACKNOWLEDGMENTS
  
  M.C. work was supported by S.N.I. (CONACyT-México).
  
  
  accel
  A. G. Riess et al. (Supernova Search Team), Astron. J. 116, 1009 (1998); S. Perlmutter et al. (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999); G. Hinshaw et al. (WMAP Collaboration), Astrophys. J. Suppl. 208, 19 (2013); N. Suzuki et al. (Supernova Cosmology Project), Astrophys. J. 746, 85 (2012).
  
  carlip
  S. Carlip, Phys. Rev. Lett. 123, 131302 (2019).
  
  planck
  N. Aghanim, et al. (Planck Collaboration), Astron. Astrophys. 641, A6 (2020).
  
  tension
  A. G. Riess, S. Casertano, W. Yuan, L. M. Macri and D. Scolnic, Astrophys. J. 876, 85 (2019); A. G. Riess, W. Yuan, S. Casertano, L. M. Macri and D. Scolnic, Astrophys. J. Lett. 896, L43 (2020).
  
  interact
  S. Pan, W. Yang and A. Paliathanasis, Mon. Not. R. Astron. Soc. 493, 3114 (2020); M. Lucca and D. C. Hooper, arXiv:2002.06127 [astro-ph.CO].
  
  wang
  B. Wang, E. Abdalla, F. Atrio-Barandela and D. Pavón, Rept. Prog. Phys. 79, 096901 (2016).
  
  temp
  M. M. Ivanov, Y. Ali-Haimoud and J. Lesgourgues, Phys. Rev. D 102, 063515 (2020).
  
  torsion
  F. Izaurieta, S. Lepe and O. Valdivia, Phys. Dark Univ. 30, 100662 (2020).
  
  prod1
  N. Herring and D. Boyanovsky, Phys. Rev. D 101, 123522 (2020).
  
  prod2
  Y. Ema, K. Nakayama and Y. Tang, J. High Energ. Phys. 2018, 135 (2018).
  
  us1 
  V. H. Cárdenas, M. Cruz, S. Lepe, S. Nojiri and S. D. Odintsov, Phys. Rev. D 101, 083530 (2020).
  
  maartens
  R. Maartens, arXiv:astro-ph/9609119.
  
  manuel
  M. Gonzalez-Espinoza and D. Pavón, Mon. Not. R. Astron. Soc. 484, 2924 (2019).
  
  andronikos
  W. Yang, S. Pan, E. Di Valentino, A. Paliathanasis and J. Lu, Phys. Rev. D 100, 103518 (2019).
  
  spin
  S. Bhadury, W. Florkowski, A. Jaiswal, A. Kumar and R. Ryblewski, arXiv:2008.10976 [nucl-th].
  
  cmb
  S. Kolekar, S. Shankaranarayanan and S. M. Chitre, arXiv:1912.06138 [gr-qc].
  
  reheat
  W. Zimdahl and D. Pavón, Mon. Not. R. Astron. Soc. 266, 872 (1994).
  
  lima
  L. R. W. Abramo and J. A. S. Lima, Class. Quantum Grav. 13, 2953 (1996).
  
  cosmography
  Y. L. Bolotin, V. A. Cherkaskiy, M. I. Konchatnyi, S. Pan and W. Yang, arXiv:2008.09602 [gr-qc].
  
  entropic
  H. Gohar and V. Salzano, arXiv:2008.09635 [gr-qc].
  
  lima1
  M. O. Calvao, J. A. S. Lima and I. Waga, Phys. Lett. A 162, 223 (1992).
  
  lima2
  J. A. S. Lima and A. S. M. Germano, Phys. Lett. A 170, 373 (1992).
  
  lima3
  J. A. S. Lima, R. Portugal and I. Waga, Phys. Rev. D 37, 2755 (1988).
  
  lima4
  R. Silva, J. A. S. Lima and M. O. Calvao, Gen. Rel. Grav. 34, 865 (2002).
  
  generalrate
  M. P. Freaza, R. S. de Souza and I. Waga, Phys. Rev. D 66, 103502 (2002).
  
  generalrate1
  R. O. Ramos, M. Vargas dos Santos and I. Waga, Phys. Rev. D 89, 083524 (2014).
  
  Nunes:2015rea 
  R. C. Nunes and D. Pavón, Phys. Rev. D 91, 063526 (2015).
    
  victor
  V. H. Cárdenas, Eur. Phys. J. C 72, 2149 (2012).
  
  victor1
  V. H. Cárdenas, arXiv:0812.3865[astro-ph].
  
  Nojiri:2005sr 
  S. Nojiri and S. D. Odintsov, Phys. Rev. D 72, 023003 (2005).
  
  little1
  P. H. Frampton, K. J. Ludwick and R. J. Scherrer, Phys. Rev. D 84, 063003 (2011).
  
  little2
  I. Brevik, E. Elizalde, S. Nojiri and S. D. Odintsov, Phys. Rev. D 84, 103508 (2011).
  
  little3
  F. Contreras, N. Cruz, E. Elizalde, E. González and S. D. Odintsov, Phys. Rev. D 98, 123520 (2018).
  
  avelino
  A. Avelino and R. P. Kirshner, Astrophys. J. 828, 35 (2016).
  
  caldwell
  R. R. Caldwell, Phys. Lett. B 545, 23 (2002).
  
  baryons
  G. Steigman, R. C. Santos and J. A. S. Lima, J. Cosmol. Astropart. Phys. 09, 033 (2009).
  
  panteon
  D. M. Scolnic, et al. ,
  
  Astrophys. J. 859, no.2, 101 (2018).
  
  
  
  PDG
  M. Tanabashi, et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018) and 2019 update.
  
  des
  M. A. Troxel, et al. (Dark Energy Survey Collaboration), Phys. Rev. D 98, 043528 (2018).
  
  is1
  W. Israel and J. M. Stewart, Ann. Phys. (N.Y.) 118, 341 (1979). 
  
  hiscock
  W. A. Hiscock and J. Salmonson, Phys. Rev. D 43, 3249 (1991).
  
  brevik
  I. Brevik and O. Gron, Astrophys. Space Sci. 347, 399 (2013).
  
  lindblom
  W. A. Hiscock and L. Lindblom, Ann. Phys. 151, 466 (1983).
  
  eckart
  C. Eckart, Phys. Rev. 58, 267 (1940)  ibid 58, 919 (1940); M. Cataldo, N. Cruz and S. Lepe, Phys. Lett. B 619, 5 (2005); A. Avelino, Y. Leyva and L. A. Ureña-López, Phys. Rev. D 88, 123004 (2013)   A. Avelino and U. Nucamendi, J. Cosmol. Astropart. Phys. 04, 006 (2009)  A. Avelino and U. Nucamendi, J. Cosmol. Astropart. Phys. 08, 009 (2010)   A. Montiel and N. Bretón, J. Cosmol. Astropart. Phys. 08, 023 (2011); T. Padmanabhan and S. M. Chitre, Phys. Lett. A 120, 433 (1987); N. D. Jerin Mohan, A. Sasidharan and T. K. Mathew, Eur. Phys. J. C 77, 849 (2017); A. Hernández-Almada, Eur. Phys. J. C 79, 751 (2019).  
  
  odintsov1
  S. Nojiri, S. D. Odintsov and S. Tsujikawa, Phys. Rev. D 71, 063004 (2005).
  
  phproduction
  J. M. Cline, S. Jeon and G. D. Moore, Phys. Rev. D 70, 043543 (2004). 
  
  phproduction1
  R. R. Caldwell, M. Kamionkowski and N. N. Weinberg, Phys. Rev. Lett. 91, 071301 (2003).
  
  phproduction2
  S. M. Carroll, M. Hoffman and M. Trodden, Phys. Rev. D 68, 023509 (2003).
  
  truncated
  D. Pavón, D. Jou and J. Casas-Vázquez, Ann. Inst. Henri Poincaré A 36, 79 (1982); L. P. Chimento and A. S. Jakubi, Class. Quantum Grav. 10, 2047 (1993); W. Zimdahl, Phys. Rev. D 53, 5483 (1996); A. A. Coley, R. J. van den Hoogen and R. Maartens, Phys. Rev. D 54, 1393 (1996); V. Méndez and J. Triginer, J. Math. Phys. 37, 2906 (1996); L. P. Chimento, A. S. Jakubi and D. Pavón, Class. Quantum Grav. 16, 1625 (1999); L. P. Chimento, A. S. Jakubi and D. Pavón, Int. J. Mod. Phys. D 9, 43 (2000).
  
  svalue
  N. Cruz, E. González, S. Lepe and D. Sáez-Chillón Gómez, J. Cosmol. Astropart. Phys. 12, 017 (2018).
  
  perturbations
  O. F. Piattella, J. C. Fabris and W. Zimdahl, J. Cosmol. Astropart. Phys. 05, 029 (2011).
  
  bang
  M. Cruz, N. Cruz and S. Lepe, Phys. Rev. D 96, 124020 (2017).
  
  grandon
  V. H. Cárdenas, D. Grandón and S. Lepe, Eur. Phys. J. C 79, 357 (2019); M. Cruz, S. Lepe and G. Morales-Navarrete, Class. Quantum Grav. 36, 225007 (2019).
  
  chemical
  M. Cruz, S. Lepe and S. D. Odintsov, Phys. Rev. D 98, 083515 (2018).
  
  prllast
  Y. Minami and E. Komatsu, Phys. Rev. Lett. 125, 221301 (2020).
  
  
  