The Formulation of Scaling Expansion in an Euler-Poisson Dark-fluid Model

We present a dark fluid model described as a non-viscous, non-relativistic, rotating, and self-gravitating fluid. We assumed that the system has spherical symmetry and the matter can be described with the polytropic equation of state. The induced coupled non-linear partial differential equation system was solved by using a self-similar time-dependent ansatz introduced by L. Sedov and G. I. Taylor. These kinds of solutions were successfully used to describe blast waves induced by an explosion since the Guderley-Landau-Stanyukovich problem. We showed that these kinds of solutions can provide new solutions that are consistent with the Newtonian cosmological framework. We have found that such solutions can be applied to describe normal-to-dark energy on the cosmological scale.


Introduction
In the second half of the 20th century, various self-similar solutions have been found following Gottfried Gurderley's famous discovery of spherically symmetric self-similar solutions that describe an imploding gas that collapses to the center [1].In this paper, we use the kinds of self-similar solutions that were independently found by Leonid Ivanovich Sedov and Sir Geoffrey Ingram Taylor during the 1940s [2,3].Despite the fact that such models have been well-known for decades, they have recently received attention again.This ansatz has been already applied successfully in several hydrodynamical systems, such as in the context of the three-dimensional Navier-Stokes and Euler equations [4], heat equations [5,6], and star formation [7].In addition, the concept of self-similarity has a wide range of applications in general relativity.Homothetic solutions were first introduced by Cahill and Taub [8], and have been studied extensively for such different topics as asymptotic solutions in cosmology [9] and the gravitational collapse of black holes [10].
The existence of the dark matter was first proposed by the Dutch astronomer Jacobus Cornelius Kapteyn [11], and became widely known through Zwicky's famous work from 1933 [12].During the second half of the century, solid experimental evidence was provided by Vera Rubin, Ken Ford, and others [13,14].However, the general existence and specified properties of dark matter remain one of the most disputed topics in theoretical astrophysics.Dark fluid is a theoretical attempt to describe the properties of dark matter and its unification with dark energy into one hypothesized substance [15].
Our goal is to use the Sedov-Taylor ansatz to describe the time evolution of a dark fluid-like material characterized by a coupled nonlinear partial differential equation system.
In our model, we study one of the simplest dark fluid materials, described by a polytropic (linear) equation of state.The dynamical evolution of the dark fluid is governed by the Euler equation and the gravitational field is described by the corresponding Poisson equation.We find time-dependent scaling solutions of the velocity flow, density flow, and gravitational fields, which can be good candidates to describe the evolution of the gravitationally coupled dust-like dark matter in the universe.We show that these kinds of solutions are consistent with the Newtonian Friedmann equations.They satisfy the mass conservation and acceleration equations and provide a similar expansion rate of the universe as the traditional solutions.The aim of this study is to broaden the knowledge of time-dependent self-similar solutions in these dark fluid models by improving upon and extending our previous model [16].We tested our model on cosmological scales; moreover, we studied how the obtained evolution of the Hubble parameter and the ratio of the density parameters resemble other physical models.

The Model
We consider a set of coupled nonlinear partial differential equations which describe the non-relativistic dynamics of a compressible fluid with zero thermal conductivity and zero viscosity [17], These equations are the continuity, Euler equation, and equation of state (EoS), respectively.We assume that the system has spherical symmetry, and we are interested in solving it in one dimension.If we assume the fluid is ideal and the system has spherical symmetry, we can reduce the multi-dimensional partial differential equation (PDE) system into the one-dimensional radial-dependent one: Here, the dynamical variables are ρ = ρ(r, t), u = u(r, t), and P = P(r, t), respectively referring to the density, radial velocity flow, and pressure field distributions, with g being the radial component of an exterior force density.As we noted briefly in the introduction, we use the following general linear equation of state: Several forms of the EOS are available in astrophysics, and polytropic ones have been successfully used in the past; see Emden's famous book [18].A great variety of applications can be found [19].In Equation ( 3), the w parameter can vary depending on the type of matter that governs the system's evolution.Traditionally, w = 0 is used; this value corresponds to the EoS for ordinary non-relativistic matter or cold dust.For our case, we choose a negative value for w which leads to different kinds of dark-fluid scenarios, as was presented in detail by Perkovic [20].In this paper, we choose w = −1 which represents the simplest case of an expanding universe governed by dark matter.Smaller values could cause a "big rip".The adiabatic speed of sound can be evaluated from Equation (3), and it is easy to show that it is constant: which is a necessary physical condition.Furthermore, let us assume that we have an additional self-gravitating term in Equation (2b).In this case, the exterior force density g can be expressed in the following way: where Φ = Φ(r, t) is the Newtonian gravitational potential and satisfies the Poisson equation, which is coupled to the previously proposed PDE system [21] ∇ where G is the universal gravitational constant, which is set to unity in our further calculations.Note that we can add an additional constant term Λ to Equation (1b) which plays a similar role to the cosmological constant in Einstein's equations: Below, we show that this constant cannot be used, as it does not lead to a consistent self-similar solution, which is what we are looking for here.Note that this observation in our model can be an indirect proof of the nonexistence of the static Universe picture.We can extend the exterior force density further with a rotating term.In this case, we add a phenomenological rotation term to Equation (2b), meaning that the equation takes the following modified form: where ω is a dimensionless parameter that describes the strength of the rotational effect and θ is the polar angle.We assume that the rotation is slow; therefore, we can expect that the spherical symmetry is not broken.This statement is satisfied if the ω parameter is sufficiently small, implying that the rotational energy is negligible compared to the gravitational energy.The self-similar analysis of various rotating and stratified incompressible ideal fluids was investigated in two Cartesian coordinates in [22].Note that the geometrized unit system (c = 1, G = 1) was applied for the calculations below, which can be converted to other units as well.See Appendix A for more details and unit conversations.

Scaling Solution and Sedov-Taylor Ansatz
We would like to find and study analytic solutions of equations by applying the long-established self-similar ansatz by Sedov and Taylor [2,3], which can be expressed in the following form: where r means radial and t means time dependence.Note that the so-called shape functions ( f , g, h) only depend on rt −β ; thus, we introduce a new variable where ζ is a dimensionless quantity in geometrized units.The (not yet determined) exponents are called similarity exponents (α, β, γ, and δ), and have physical relevance; β describes the rate of spread of the spatial distribution during the time evolution (if the exponent is positive) or contraction (if β < 0).In addition, the other exponents describe the rate of decay of the intensity of the corresponding field.Solutions with integer exponents are called self-similar solutions of the first kind, while the second kind denotes the non-integer ones.Self-similarity is based on the concept that physical quantities preserve their shape during time evolution.A general description of the properties of these types of scaling solutions can be found in our previous publication [16].
We assume that the shape functions are sufficiently smooth and continuously differentiable at least twice in ζ over the entire domain.Thus, we calculate the relevant time and space derivatives of the shape functions and substitute them into Equation (2).As a consequence, we usually obtain an overdetermined algebraic equation system for the similarity exponents.Other possible scenarios may play out as well; these have been presented in detail in [16].We obtain the following numerical value for the exponents for both the non-rotating and rotating cases: α = 0, β = 1, γ = 2, and δ = 0.If we add the Λ constant to the Euler equation, such a solution for the similarity equation cannot be found.From these results, it is evident that the dynamical variables such as the velocity, gravitational potential, and density flow have spreading properties.Our physical intuition says that spreading is somehow similar to expansion, which is a basic property of the universe at astronomical or cosmological scales.By substituting the obtained numerical values of the similarity exponents, we can reduce the induced PDE system into an ordinary differential equation (ODE) system that depends only on the ζ independent variable.We find that the obtained equation system has the following form: It can easily be noticed that the ordinary differential equation system presented in Equation ( 11) cannot be solved analytically.For linearized nonautonomous ordinary differential equation systems, the stationary point of the phase space can be found, and we can say something about the general asymptotic behavior of the solutions as well [23].Nonetheless, there is no generally known method for nonlinearized nonautonomous differential equation systems.Moreover, the existence and uniqueness of smooth solutions has not yet been proven in multiple dimensions.Therefore, it is a reasonable approach to solve the obtained ordinary differential equation system in Equation ( 11) numerically for a large number of parameter sets (based on physical considerations) in order to explore the behavior of the solution of the system with different boundary and initial conditions.One example of such a numerical solution can be seen in Figure 1.As an example, at a specific parameter and initial condition set, the shape function of the velocity f (ζ) is almost linear and increasing after a short decrease, providing a hint of Hubble expansion-like behavior.The shape function g(ζ) is asymptotically flat after a quick ramp-up, corresponding to the conservation of matter.The last shape function h(ζ) has an increasing polynomial trend with a slight positive exponent, which is connected to the gravitational potential.To obtain a sufficiently smooth numerical solution, we solved the ODE system using an adaptive numerical integration provided by Wolfram Mathematica 13.1 [24].For all of our calculations, the integration limits were ζ 0 = 0.001 and ζ max = 40, the same as in [16].As has been said before, we established initial conditions to obtain the numerical solution; for this reason, we use ranges R of f (ζ 0 ) = 0.005 − 0.5 and g(ζ 0 ) = 0.001 − 0.1, while for the second-order differential equation we have This choice of initial conditions reflects that, first it is physically reasonable that the density flow range is R(g) ⊂ R + and finite.Recent results suggest that dark fluid could possibly have negative mass [25]; however, in the present model this leads to singular solutions.Second, our choice of initial velocity flow R( f ) ⊂ R + makes for an initially radially expanding fluid.We have seen that if the initial value for f (ζ) and g(ζ) is set outside of the previously given range, the solution of the differential equation becomes singular.Moreover, we have seen that the variation in the initial condition corresponding to the shape function of the gravitational potential does not affect the trend of the time evolution of the system, as it only causes vertical shifts.Therefore, we set the initial numerical value equal to zero.
We are interested in finding the solution of the ODE system as a function of the spatial and time coordinates.We transform our single-variable numerical solutions into two-variable functions, for which we use the inverted form of Equation (10).It can easily be noticed that if we look at the shape of the ansatz, the solution has a singularity at t = 0. Thus, we use the 0.001 ≤ t ≤ 25 and 0.001 ≤ r ≤ 25 domains to obtain the space-and time-dependent initial dynamical functions u(r, t), ρ(r, t), and Φ(r, t).

Results
Here, we present the solutions of the self-gravitating non-relativistic dark fluid.First, we provide a detailed introduction to the global properties of the solutions in a non-rotating system.Second, we show the effect of slow rotation on the solutions.In addition, we compare the results from the two cases with each other and with the previous results in [16].Note that the spherical symmetry of the system was kept conserved for all the cases.

Non-Rotating System
In the first case, we set the ω parameter to zero and used the obtained numerical values for the similarity exponents (α, β, γ, and δ) to obtain the exact ordinary differential equation.For the numerical integration, we used the same initial conditions f (ζ 0 ) = 0.5 and g(ζ 0 ) = 0.01 applied in our previous paper for the velocity and density flows, respectively.First, we used the time and radial projection of the unknown functions to obtain a better understanding.Figure 2 illustrates the spatial and time projections of the obtained velocity, density, and gravitational potential.These results for the velocity and density flows are consistent with our initial statement that these kinds of solutions of dark fluids can be used as a model to describe the exploding system (e.g., the universe).Similar behaviors can be seen for the radial velocity and the density, which both have fast decays in time at all distances.In addition, they both have a real singularity at t = 0 due to the shape of the ansatz.However, the radial distribution has a different nature; the density increases excessively at distances near the center of the explosion, and becomes linear at larger distances.On the contrary, the velocity grows polynomially with the radial distance.It can be seen that the gravitational potential decreases hyperbolically over time and becomes asymptotically flat.The gravitational potential has a natural singularity at the origin.The local existence of global weak solutions on a domain outside of the origin in spite of the existence of a singularity has been shown by Tsunge and others.Here, we find a different radial velocity profile than in the previous non-rotating model with two equations presented in [16].Furthermore, the similarity exponents are different (α = 0, β = 1, and γ = −1 in the two-equation model).This is most likely due to the new smoother solution forming as a consequence of a result of the second derivative appearing in the Poisson equation.It can be seen that the solution depicted above in Figure 2 is numerically stable in the specified initial and boundary condition range.It is more relevant to investigate the dynamics of the complete fluid in time and space in order to understand certain general trends or physical phenomena as the function of the initial conditions.For this reason, we now evaluate the related energy densities, which are as follows: Figure 3 further illustrates the fact that the kinetic energy density has a singularity at t = 0, as we have already seen in the case of the radial velocity.It has linearly enhancing maxima at larger distances, and has a fast decay in time for all radial distances.As mentioned above, the gravitational potential is a negative polynomial in time; therefore, we can obtain the total energy density of the system.From the total energy density distribution, it is apparent that the short-time behavior of the system is dominated by the initial explosion, while the long-range structure is regulated by the gravitational potential.

Rotating System
In this section, we analyzed the effect of slow rotation compared to the non-rotating case.First, we studied the effect of the variation of the maximal angular velocity (the ω parameter).We chose the polar angle (the θ parameter at the equatorial), as it provides the largest effect.As previously stated, we assume in our further analysis that the spherical symmetry is not broken.According to this assumption, we fix the gravitational force density to be at least one magnitude larger than the centrifugal force density at every time and space (∥ f grav ∥ ≫ ∥ f centr ∥).The numerical results show that we can find a range of ω within which the constraint is fulfilled as long as the previously specified initial condition set is valid.We demonstrate that the asymptotic behavior of the numerical solution includes a significant dependence of ω on the acceptable domain (0 < ω < 0.3) of the parameters and initial conditions.
Comparing the results shown in Figure 4 with the non-rotating case in Figure 2, it is evident that slow and constant rotation does not affect the time or spatial distribution of the gravitational potential.Moreover, it can be seen that the radial density profile of the system is nearly uniform, and is identical to the previous case in that it decreases rapidly over time.Thus, we can conclude that the rotation accelerates the even distribution of the material in space and speeds up inflation.The singular behavior close to t = 0 is not affected by the rotation, as was expected.However, a significant difference can be seen in the first graph (the top left panel of Figure 4).
It can be seen that the radial profile of the velocity flow starts from zero in the origin and shows exponential growth for the short-range behavior.
Figure 5 illustrates the critical influence of the ω on the long-range asymptotic behavior of the time evolution of the velocity flow.Moreover, an increase in the ω value causes significant modifications in the radial profile for both the velocity and density flows while leaving the time evolution unaltered.In our analysis of the behavior of the obtained numerical solutions, we found that the inspected initial value range shows similar behavior; ω < 1 can be found for every initial and boundary condition in which the long-range asymptotic structure alternates.An example of this can be seen in Figure 5.Likewise, we studied the properties of the relevant dynamic variables for the previous case.The energy density associated with the rotation and the total energy is these quantities are depicted in Figure 6.

Connection to Newtonian Friedman Equation
Next, we apply the obtained scaling solution to cosmology in order to explain the evolution of the universe.Our strategy is to describe the concept of a cosmic fluid in terms of Newtonian cosmology, which we use because relativistic effects are not of great significance in this context.In this case, there is no need for a global reference point; thus, a scale factor a(t) can be introduced which contains the entire time evolution that affects our chosen reference frame.Henceforth, the relative distances in time can be denoted in this term: where R(t) is a continuous function and l is real number.In this chapter, we use classical conservation equations to show the connection between the obtained self-similar solution and the traditional Friedmann equations.

Connection to the Expansion Rate
Here, we introduce the total mass inside a radius r: where V (t) ⊂ R 3 is a sphere with radius R(t) and r ∈ (0, R(t)).Therefore, within the classical Newtonian framework, the energy and mass conservation principle are fulfilled separately: Replacing the derivation with the integral formula, we can assume that equality holds for every l: In addition, we impose the kinematic condition that We impose ρ(r, t) > 0 for ∀r ≥ R(t) and assume that the velocity u(r, t) and density ρ(r, t) are obtained from the self-similar ansatz in Equation ( 9).We expand the density and velocity flow into the following respective power series: After substituting ζ = R(t)t −β , Equation (18) becomes On the relevant space and time scale, we only consider the contributions of the finite term.The discussion of the rotating case in Section 4 requires that we use a high-order polynomial to approximate the f (η) function.Based on the analysis of the numerical results in Figure 4, we restrict ourselves to using polynomes up to the eighth order, which are well able to approximate the curves: Moreover, a further simplification can be applied within the domain of interest for ρ by describing it as a polynomial with a rational exponent ρ(ζ) ∼ Aζ κ , where κ = 6/7.Hence, Equation (20) can be rewritten as Consequently, the mass conservation equation in (15) becomes a nonautonomous firstorder differential equation that cannot be solved analytically due to its highly nonlinear terms.
We can find further simplifications in the non-rotating limit (ω → 0), as the higherorder terms are less and less significant; therefore, the velocity flow simply becomes u(r, t) ∼ u 1 ζ 1 + u 2 ζ 2 and the mass conservation equation is The above equation can be solved analytically; the solution is where µ = 1 − (α + β), ν = κ − βκ, α, β, γ are the similarity exponents and A, u 1 , u 2 the positive constants, the κ positive exponent is obtained from the solution of Equation ( 11), H 1 is an integration constant, and Γ() is the upper incomplete Gamma function [26].Equation ( 24) can be simplified by substituting the similarity exponents and constants from the non-rotating case: , where κ = 6 7 .
We can now use Hubble's law of expansion to determine the H 1 integration constant and the numerical integration for the rotating case: where H 0 = 66.6 +4.1 −3.3 km/s/Mpc is the experimental value of the Hubble constant [27].The expansion rates for the rotating (yellow line) and non-rotating (blue line) Universe are drawn in Figure 7 in the usual units of standard cosmology.The dashed line represents the present measured value of the Hubble parameter.

Connection to the Critical Density
Energy conservation has to be investigated as well, as it leads to the critical density parameter of the universe.We can apply the Newton equation to obtain the so-called acceleration equation [29]: Next, we multiply the equation by Ṙ(t) and integrate over time; then, Equation ( 27) becomes which is the energy conservation equation.In general relativity, U 0 is related to the curvature of spacetime.However, because we are working in the framework of Newtonian cosmology it is instead associated with a kind of mechanical energy.Equation ( 28) can be transformed into the following form using the definition H(t) := Ṙ(t)/R(t): where ), and we use the following definition for the critical density: The value of the dynamical constant U t determines the evolution of the system.If U t > 0, then the universe eventually re-collapses.On the other hand, if U t < 0 then the system expands to infinity.Thus, U t = 0 is a special case in which the universe shows similar expanding behavior but reaches zero velocity only in infinite time [29].Here, we are searching for an isentropic solution; therefore, we introduce the entropy conservation equation dE + pdV = 0 from [30].Using the non-steady flow internal energy formula E = Vρ(u + 1) in the c = 1 unit system, and using the equation of state from Equation (3) in [31], the entropy conservation [32] can be expressed in terms of with the solution Thus, in this model the relationship between the time-dependent Hubble parameter and the density parameter is which cannot be solved explicitly.Here, Ω 0,CDM or Ω DE,0 are density parameters related to the cold dark matter (CDM) and dark energy (DE), respectively.The relationship between the density parameters and the dynamical parameters of the self-similar solution is defined and explained in the next section.

Discussion
According to current scientific understanding, dark matter and dark energy make up about 95% of the total energy density of the observable universe today.The dark fluid theory suggests that a single substance may explain both dark matter and dark energy.The behaviour of the hypothetical dark fluid is believed to resemble that of cold dark matter on galactic scales while exhibiting similar characteristics to dark energy at larger scales [25].Predictions can be obtained from our Sedov-von Neumann-Taylor blast wave-inspired nonrelativistic dark fluid model on galactic and cosmological scales.A useful feature of this model is that the initial value problem of the reduced ordinary differential equation system is easier to handle than the boundary and initial condition problem of the original partial differential equations.To provide a reliable practical basis for our dark fluid model, we have tested the theoretical results on astrophysical scales in order to demonstrate the similar nature of the solutions on the cosmological scale.
Our solution was developed based on cosmological observations, using the Hubble law to scale the expansion of the universe.Our model includes various scaling mechanisms through the use of a Sedov-type self-similar ansatz, which allows for describing different time decay scenarios [33].In the case of a non-rotating system, we can conclude that the radial velocity profile of the solution provides Hubble-like expansion, while the nonrotating model provides inflation-like behavior u(r, t) > 1 on a long-range timescale, which cannot be physical (causal).We have seen that the high initial velocity of the dark fluid relaxes to a small, constant, and non-relativistic value on long timescales (t > 6 billion years).Another interesting feature of the model is that the gravitation potential is a negative polynomial in time, which is consistent with the distribution of the density of the dark fluid.
One notable aspect in the case of the rotating model is that it does not show superluminous behavior in the expected time range, contrary to the non-rotating model.Simultaneously, we have found that the radial profile of the density becomes saturated and close to flat at far distances from the initial point; see the top right panel of Figure 6.In addition, it is possible to set the initial conditions such that the universe today is observed as flat Euclidean with the density parameter being the critical density corresponding to the Hubble parameter H 0 .The flatness of the universe is indicated by the recent measurements of WMAP [34].The sum index uses the baryonic energy (B), dark energy (DE), and cold dark matter (CDM).We can define the matter part of Ω M = Ω CDM + Ω B along with the full Ω = ∑ i Ω i and i ∈ {B, CDM, DE} [35].If the Ω CDM ≫ Ω B relation is correct, then we can assume the following identity: Accordingly, we can determine the relevant time and radial coordinates from the obtained results, which correspond to this specific energy relation evaluated at r = R(t 0 ) and t = t 0 .In the absence of dark fluid, the universe continues to expand indefinitely, though at a gradually slowing rate that eventually approaches zero.This causes an opentopology universe.In this case, the ultimate fate of the universe is that the temperature asymptotically approaches absolute zero in a so-called "big freeze".At the same time, it is essential to mention that a weakness of our non-relativistic model is that its results are not as precise as those of relativistic Friedmann equation-based models [25].

Conclusions and Outlook
In this paper, we have studied the behaviour of self-similar time-dependent solutions in a coupled nonlinear partial differential equation system describing a non-viscous, nonrelativistic, and self-gravitating fluid (Euler-Poisson system).The reason behind the applied self-similar solutions is that they provide a very efficient method to analyze various kinds of physical systems, particularly as concerns the hydrodynamical description of systems that involve collapse and explosion.
The analysis presented in this paper is a Euler-Poisson extension of our previous model [16].We found that a Sedov-Taylor type solutions exists and that the algebraic equation obtained for the similarity exponents has only one unique solution.We have used the obtained solution to describe the behaviour of a non-relativistic dark fluid on cosmological scales, and have presented the relevant kinematical and dynamical quantities.
In addition, we have shown that our model based on self-similarity is in agreement with the Newtonian Friedmann equations.
It can easily be noticed that this model has certain limitations due to its classical nature, although it does provide relatively adequate results on the cosmological scale.We have shown that the obtained quasi-analytical solution for the evolution of the Hubble parameter is in agreement with the standard cosmological model, e.g., [28].Moreover, in the previous section we have shown that the energy ratio from our model closely resembles the ratio of Ω M (t)/Ω(t), i.e., the relevant density parameters as of today.
This model has the additional practical benefit that the calculation does not require high computing performance and resources.Therefore, it can be used to estimate the physical value of the initial and boundary values when more sophisticated theoretical or numerical simulations are used.Moreover, it can provide a reliable basis for comparing for two-and three-dimensional hydrodynamical simulations.In future research, through reasonable effort it would be possible to improve this model in order to describe relativistic matter [36], Chaplygin gas [37], and two-fluid models.

Figure 3 .
Figure 3. Numerical solutions of the velocity flow u(r, t), density flow ρ(r, t), and gravitational potential Φ(r, t) as a function of the spatial and time coordinates in the case of a non-rotating system, additionally showing the distribution of the total and kinetic energy densities.We used ζ 0 = 0.001 for numerical integration and initial conditions f (ζ 0 ) = 0.5, g(ζ 0 ) = 0.01, h(ζ 0 ) = 0, and h ′ (ζ 0 ) = 1.

Figure 5 .
Figure5.Dependence of the space and time evolution on the maximal angular velocity ω; the different lines correspond to different values of the angular velocity ω.The curves were evaluated at a particular time (left), with radial coordinates provided on the vertical axis (right).A detailed explanation is provided in the main text.

Table A2 .
Relevant astronomical quantities and their corresponding values in SI units.