New, Spherical Solutions of Non-Relativistic, Dissipative Hydrodynamics

We present a new family of exact solutions of dissipative fireball hydrodynamics for arbitrary bulk and shear viscosities. The main property of these solutions is a spherically symmetric, Hubble flow field. The motivation of this paper is mostly academic: we apply non-relativistic kinematics for simplicity and clarity. In this limiting case, the theory is particularly clear: the non-relativistic Navier–Stokes equations describe the dissipation in a well-understood manner. From the asymptotic analysis of our new exact solutions of dissipative fireball hydrodynamics, we can draw a surprising conclusion: this new class of exact solutions of non-relativistic dissipative hydrodynamics is asymptotically perfect.


Introduction
The nearly perfect fluid behaviour in high-energy heavy ion collisions has been observed by the Brookhaven National Laboratory's Relativistic Heavy Ion Collider (RHIC) by the four RHIC collaborations, BRAHMS [1], PHENIX [2], PHOBOS [3], and STAR [4]. Currently, the mainstream of hydrodynamical modelling of high-energy heavy ion reactions seems to focus mostly on the numerical solutions of dissipative relativistic hydrodynamics. However, several analytic exact solutions are famous or very well known: the status of the field including applications of exact solutions of fireball hydrodynamics has been reviewed recently in ref. [5], including a discussion of several of the open questions of this field.
On the topic of describing relativistic heavy ion collisions with exact solutions of hydrodynamics, the most well-known papers deal with perfect fluid hydrodynamics for a 1 + 1 dimensional, longitudinally expanding fireball: for example, the boost-invariant Hwa-Bjorken solution [6,7] or an accelerating solution with finite rapidity density by Landau and Belenkij [8]. One of the main focuses of seeking new exact solutions of relativistic hydrodynamics is to generalise the Hwa-Bjorken solution to a 1 + 3 dimensional, accelerating flow profile.
One of the first famous 1 + 3 dimensional, spherically symmetric exact solutions of fireball hydrodynamics is the Zimányi-Bondorf-Garpman (ZBG) solution [9]. In the ZBG solution, a finite density and temperature profile is described, which both vanish after a time-dependent cutoff distance given by the scale parameter R(t). The ZBG solution is accelerating, and the flow-field corresponds to a spherically symmetric Hubble flow: v(r, t) =Ṙ (t) R(t) r. The ZBG solution has been generalised to a Gaussian density and a homogeneous temperature profile in ref. [10] and subsequently to an arbitrary temperature and a matching density profile in ref. [11], in both cases keeping the spherically symmetric Hubble flow field. The present manuscript generalises these solutions for an arbitrary, temperature-dependent speed of sound, using a thermodynamically consistent equation of state as proposed in ref. [12], and for a rather general, temperature-dependent shear and bulk viscosities, using the non-relativistic Navier-Stokes equations. For the sake of simplicity and clarity, throughout the body of this work, we assume the validity of a spherically symmetric Hubble flow.
In recent simulations of relativistic viscous hydrodynamics, allowing the presence of small specific shear [13] and bulk [14] viscosities provides an acceptable description of the sQGP (strongly coupled quark-gluon plasma), so the importance of viscous corrections to the Hwa-Bjorken solution and other exact solutions of relativistic perfect fluid hydrodynamics also emerges.
The effect of shear viscosity with longitudinal acceleration has been investigated recently by a perturbative, relativistic solution of Navier-Stokes and Isreal-Stewart theory in ref. [15]. Dissipation in relativistic hydrodynamics has been addressed in an exact form as well by two recent manuscripts [16,17]. These works discuss analytic results on the effect of bulk viscosity on 1 + 3 dimensional exact solutions of the relativistic Navier-Stokes as well as the Israel-Stewart equations, generalising the Hwa-Bjorken solution for non-vanishing kinematic bulk viscosities in 1 + 3 dimensions [16,17]. The velocity field of these solutions is also a Hubble flow, which describes reasonably well the asymptotic velocity profile of small and exploding fireballs with possible applications even in high-energy heavy ion collisions [17].
It seems to us that the new solutions of refs. [16,17] provide a great tool to investigate the late time effects of bulk viscosity in the hydrodynamic evolution. However, the fundamental equations of dissipative, relativistic hydrodynamics are debatable, and so far, several schemes have been proposed, but without a fundamental and generally accepted method, that are able to provide a stable and causal theory, and at the same time have also the correct non-relativistic limit: the non-relativistic Navier-Stokes equations. This motivated us to investigate the non-relativistic and spherically symmetric limit of the relativistic exact solutions described in refs. [16,17]. It turned out that in this limit, we obtained new, simple, and exact solutions of the non-relativistic Navier-Stokes equations. This solution is described in this work. This solution also allows for a straightforward generalisation to a rotating and directional Hubble type flow-field. The discussion of these non-spherical solutions is clearly beyond the scope of the present manuscript: we plan to present these solutions elsewhere.
In this paper, we present a new family of exact and analytic solutions of non-relativistic, dissipative hydrodynamics. As a result of the non-relativistic regime, we do not seek for a direct application of these solutions to describe experimental data, although such applications are possible even in the perfect fluid limit, as demonstrated before in refs. [9,18]. In this work, we opted to remain in the entirely academic framework and do not compare our theoretical results to experimental data. Although we apply a spherically symmetric, three-dimensional Hubble flow, the scale parameter R = R(t) is a time-dependent parameter, and its second derivative, hence the acceleration of the expansion is non-vanishing in some of the parametric solutions that we present. We also discuss coasting, accelerationless solutions for a homogenous pressure profile.
The spherical symmetric Hubble profile excludes the possibility of the rotation of the fireball, and it also causes the cancellation of the shear viscosity term from the dynamical equations of non-relativistic hydrodynamics. To simplify further the dissipative dynamics, we neglect possible heat conductivity effects. These simplifications allow us to focus on the main aim of this manuscript: to investigate the asymptotic effects of the bulk viscosity in the simplest possible case of fireball hydrodynamics.

Navier-Stokes Equations of Non-Relativistic, Viscous Hydrodynamics
The dynamical equations of non-relativistic, viscous hydrodynamics are discussed below. The continuity equation of the conserved particle density n reads as where n(r, t) ≡ n is a function of the spatial coordinates r = r x , r y , r z and time t, while v(r, t) ≡ v stands for the velocity field. We allow for a compressible expansion, ∇v = 0, keeping in mind that the fireball formed in the Little Bangs of heavy ion collisions corresponds to a quickly expanding fluid. We start from the Navier-Stokes equation for the conservation of energy, which is expressed as where ε(r, t) ≡ ε is the energy density, the pressure is denoted by p(r, t) ≡ p, T(r, t) ≡ T is the temperature, the heat conductivity coefficient is denoted by λ, while the shear and bulk viscosity coefficients are denoted by η and ζ, respectively, and D is Cauchy's strain tensor, which is defined as We consider the case of vanishing heat conduction, corresponding to the assumption of λ ≡ 0. The Navier-Stokes form of the Euler equation stands for the momentum conservation: where the single particle mass is denoted by m. Here, we assume that the bulk viscosity ζ may depend on position, so care must be taken when evaluating the first term on the right-hand side of the above equation.
The microscopic properties of the flowing matter are its thermodynamical features. These are characterised by an equation of state, which closes the above set of partial differential equations of non-relativistic hydrodynamics. In this paper, we follow refs. [12,[19][20][21][22] and assume that the equation of state can be written as These equations of state are thermodynamically consistent [12]. In the next section, we detail the temperature dependence of the speed of sound that follows from these equations of state. Although similar equations of state have been utilised in exact analytic solutions of both non-relativistic and relativistic hydrodynamics before, the corresponding temperaturedependent speed of sound has not been described or detailed in earlier papers, as far as we know.

Temperature Dependence of the Speed of Sound
Now, let us clarify the meaning of the κ(T) function. The simplest case is if this function is a constant: if κ(T) becomes independent of the temperature T, it is denoted as κ, and its relationship to the adiabatic index γ can be easily found as where C p is the heat capacity at constant pressure and C V stands for the heat capacity if the volume of the system is kept constant. This is also the reason why the adiabatic index γ is also called the heat capacity ratio. The well-known formula of the speed of sound reads as where σ stands for the entropy density. Thus, in general, the speed of sound may be dependent on both the temperature and on the chemical potential of conserved charges. In case of ideal gases as well as in case of our selected equations of states, Equations (5) and (6), the connection between κ and the speed of sound can be obtained from Equation (5), for a constant, temperature independent κ as follows: In the case of a temperature-dependent κ(T) function, the formula for the temperature dependence of speed of sound yields a more general expression: Based on Equations (7) and (9), we can define a temperature-dependent adiabatic index, and we have found that γ(T) still can be expressed as the ratio of generalised, temperature-dependent heat capacities: A further generalisation is possible for the case when the particle mass depends on only the temperature. The derivation is similar, and the formulae remains nearly unchanged. The speed of sound with a temperature-dependent κ(T) and particle mass m(T) is obtained as: We have also investigated the case of a multi-component matter, which is composed of a mix of various hadrons with temperature-dependent masses. In such a medium, the speed of sound depends on the average mass of the hadrons instead of the single particle mass: where Equation (13) is expressed in a general form for multi-component hadronic matter with temperature-dependent κ parameter and particle masses. It is clear that this form of the speed of sound is independent of the temperature (in)dependence of both the adiabatic index (hence κ as well) and the (average) particle mass: in the case of constant, temperatureindependent mass m, and κ or adiabatic index γ, Equation (13) reduces to the form for ideal gases, as shown in Equation (9).

Scale and the Continuity Equations for Spherically Symmetric Hubble Flow
We search for spherically symmetric, exact solutions of the Navier-Stokes equations, where the velocity field is a Hubble flow with a time-dependent Hubble parameter: where R(t) ≡ R is the scale of the expanding fireball. In a Hubble flow, R is a function of the time only. We seek self-similar solutions; thus, we introduce a scaling variable s that satisfies the scale equation: For the above defined spherically symmetric Hubble flow, a spherically symmetric solution of the scale equation is: and the D matrix, defined in Equation (3), becomes diagonal: Using Equation (15) for the velocity field, the solution of the continuity equation is: where V is an arbitrary function of the s scale variable, n 0 stands for n(0, t 0 ), the initial value of the R scale is R 0 ≡ R(t 0 ), and d = 3 is the number of spatial dimensions. The spherically symmetric Hubble flow field given by Equation (18) and the spherically symmetric, generic solution for the density of the conserved charge, given by Equation (19) will be common in various classes of exact solutions of the Navier-Stokes equations detailed in the subsequent parts of this manuscript. Due to the spherical symmetry of the selected velocity profile, the effects of shear viscosity cancel.

New, Exact Solutions for a Generic, Temperature-Dependent Speed of Sound
In this section, we describe our general results, gradually introducing more and more simplifying assumptions. First of all, let us note that even in the κ(T) = κ 0 constant case, the considered equation of state leads to a temperature-dependent speed of sound, as given in Equation (9). For a generic, temperature-dependent pressure to energy density ratio κ(T), the temperature dependence of the speed of sound may be even more complicated; see for example, Equation (12) and further details in Section 3. In this subsection, we collect results that are obtained for a generic, temperature-dependent pressure to energy density ratio κ(T).
When searching for solutions with a temperature-dependent κ(T) function, we utilise earlier results of solving the equations of non-relativistic perfect fluid hydrodynamics with a temperature-dependent κ(T) and corresponding temperature-dependent speed of sound, in particular the properties of the solutions in refs. [12,20,21,[23][24][25]. Based on these perfect fluid solutions, our ansatz for the self-similar temperature profile is the following, factorised form: where T (s) ≥ 0 is an arbitrary non-negative function of s, normalised as T (0) = 1. The non-trivial time dependence is described by the factor f T (t), and we denote the initial temperature at s = 0 at the initial time t 0 by T 0 . For a non-zero conserved charge or bariochemical potential, the pressure is thus determined by the equation of state, p = nT as: where p 0 = p(0, t 0 ) = n 0 T 0 . Using these considerations, the energy equation, Equation (2), can be rewritten as follows: This equation can be further simplified and solved if the temperature is spatially homogeneous, T (s) = 1, keeping a generic κ(T) function. For such a spatially homogeneous temperature profile, the left-hand side of the above equation depends on time only. For such a temperature-dependent pressure to energy density ratio κ(T), the temperature-dependent adiabatic index γ(T) is determined in Section 3, as given by Equation (11), and the temperature-dependent speed of sound by Equation (12). Equation (22) is a non-linear but ordinary differential equation. Its left-hand side depends on time, while its right-hand side depends not only on time but also on the ratio of the bulk viscosity to pressure. Thus, this equation can be solved if this ratio is assumed to be a constant: The resulting equation, valid for spatially homogeneous temperature profiles, T (s) = 1 reads as follows: In the Euler equation, the second derivatives cancel because of the special form of the velocity field, so we get back the same equation that corresponds to perfect fluids: Using a spherical Hubble profile and our ansatz for the temperature field, introduced by Equation (20), the Euler equation is reduced to the following second order, ordinary differential equation: Here, a constant of integration is denoted by C E , which is a constant parameter that also enters the following differential equation for the scaling functions: where the prime in V (s) and T (s) denotes a derivation with respect to s. Equation (27) indicates that the T (s) and V (s) profile functions for the density and for the temperature are not independent from each other: they must be chosen corresponding to a matching initial condition. The differential equation for the scaling functions can be integrated, and the scaling function for the density can be expressed with the help of the constant of integration and the scaling function of the temperature: Thus, it is clear that spatially homogeneous temperature profiles with T (s) = 1 correspond to Gaussian density profiles, V (s) = exp − C E 2 s 2 . Thus, we have reduced the complex set of partial differential equations of nonrelativistic, dissipative hydrodynamics with a temperature-dependent pressure to energy ratio κ(T) and the corresponding temperature-dependent speed of sound given by Equation (12), for spatially homogeneous temperature profiles, to a family of coupled and non-linear, but ordinary differential equations, Equations (22) and (26). These equations can be easily solved with broadly accessible, numerical software packages such as MATLAB, Maple, or Mathematica. Note that these solutions describe an accelerating fluid, where the second time derivative of the scale parameter,R(t), is a non-vanishing function of time, as follows from Equation (26) for a non-vanishing constant of integration, C E = 0. Accelerationless solutions are also obtained, according to Equation (26), they correspond to the C E = 0 special case.
For the sake of completeness, let us also mention the equation for the σ ≡ σ(r, t) entropy density, which describes the entropy production due to dissipation. This equation reads as:

Solutions for a Temperature-Independent Pressure to Energy Density Ratio
From now on, let us consider a temperature-independent energy density to pressure ratio κ(T) = κ 0 . In this manuscript, we are considering three-dimensional, finite fireballs; however, such equations of states are frequently used, for example in the case of the conformal Gubser flows, that correspond to the κ 0 = 3 case as detailed in [26,27]. For such a temperature-independent energy density to pressure ratio κ 0 , Equation (22) simplifies as

Analytic Solutions for a Spatially Homogeneous Pressure Distribution
We can provide exact analytic solutions of the above system of differential equations if the pressure distribution is spatially homogeneous. Although such a pressure distribution is academic, this solution is interesting, as its certain properties may be inherited by solutions with spatially inhomogenous pressure profiles as well. Let us consider two important remarks. First, consider that in this case, when the pressure depends only on time, the coefficient of the bulk viscosity can be chosen to be any pressure-dependent function because of the energy equation, Equation (30): Second, the condition of homogeneity requires that V (s)T (s) = 1. So, according to Equation (28), the constant of integration has to vanish, C E = 0. Then, the Euler equation is simplified toR = 0.
This equation reflects the fact that without a pressure gradient, the fireball expands with a constant rate,Ṙ =Ṙ 0 . Thus, the time-dependent scale of the fireball is expressed by the function R(t) =Ṙ(t − t 0 ). Using this result, the energy equation is reduced to the following: If we assume that the ratio of the bulk viscosity to pressure is a constant (ζ ∝ p), then we obtain an explicit solution for the f T (t) function: where the ratio ζ 0 /p 0 is the proportionality factor between ζ and p. Utilising this form, exact results can be obtained for the temperature, the conserved charge, and the pressure: This result is not only an analytic solution of the non-relativistic Navier-Stokes equations, but it turns out that this exact solution is the non-relativistic limit of a new, relativistic solution of dissipative hydrodynamics with accelerationless Hubble profile, as detailed in ref. [17]. Note that in Equation (35), an additional integration constant is allowed that controls the value of R(t) at t 0 , but in this case, we fixed this constant to 0. This result allows for an asymptotic analysis, which is given in more details in Section 7.1. The key observation is coming from the analysis of the second factor of Equation (34), which tends to a constant multiplicative factor for late times, t t 0 . According to Equation (20), the time dependence of the temperature is governed by f T (t), which is illustrated in Figure 1 for s = 0. In the left panel, we have shown five different curves with the same initial temperature, and the initial kinematic bulk viscosity is varied. The effect of entropy production is clearly visible compared to the black curve, which corresponds to perfect fluid. In the right panel, we varied not only the initial value of kinematic bulk viscosity but also the initial temperature. In this case, each coloured curve tends to the same asymptote shown by a solid black line. Thus, these curves approach the same final state. For a detailed explanation of this phenomenon, see Section 7. In the left panel, we set the same initial temperatures, but in the right panel, the curves start from different initial conditions, and each of them approach the solid black line, the perfect fluid asymptote.

Analytic Solutions for a Spatially Inhomogeneous Pressure Profile
Now, let us discuss the case where the pressure also depends on the spatial coordinates through the scale variable s. In Equation (30), one can realise that the left side of the equation depends only on time, but the right side has coordinate dependence, too, because both the pressure and the shear viscosity may depend on the scaling variable s. However, this s dependence of the right hand side is eliminated if the bulk viscosity coefficient is proportional to the pressure. Let us assume that the bulk viscosity has a special form: In this case, the energy equation, Equation (30), becomes an ordinary differential equation that depends on time only. In addition, let us factor out a trivial time dependence that is due to the decrease of the temperature due to the radial expansion. Let us rewrite the time-dependent factor f T (t) in the following form: where g T (t) stands for a dissipative correction in the time dependence of the temperature field. For a vanishing bulk viscosity coefficient, g T (t) ≡ 1. This ansatz is based on known non-relativistic, perfect fluid solutions with Hubble flow, as detailed in refs. [10,12,21,[28][29][30]. Using this ansatz, the energy equation of Equation (30) can be rewritten in terms of g T (t) and R(t) as follows:ġ The Euler equation can be also rewritten in terms of the dissipative correction g T as follows: Equations (42) and (43) are time-dependent, coupled, and non-linear but ordinary differential equations for R(t) and g T (t). These equations can be readily solved with the help of generally available mathematical packages such as Maple, MATLAB, or Mathematica. Furthermore, their structure is particularly clear and allows for an asymptotic analysis that corresponds to their late time, t t 0 , behaviour. This is the subject of the analysis of the next section.

Discussion: Attractor Behaviour in Other Hydrodynamical Solutions
Before we discuss the asymptotic properties of our newly described exact viscous solutions, let us mention that these asymptotic properties have a kind of attractor behaviour. Finding attractors in hydrodynamic systems as well as finding hydrodynamical attractors in non-equilibrium systems is a broad topic of great current research interest that is impossible to fully review in the current, rather academic manuscript with more limited scope.
In ref. [31], the Boltzmann equation was solved in a relaxation time approximation for an azimuthally symmetric radially expanding boost-invariant conformal system that is undergoing a longitudinally boost-invariant, radial Gubser flow [27], and a hydrodynamic attractor solution was found in various approximations of relativistic viscous hydrodynamics. Analytic solutions of dissipative relativistic spin hydrodynamics were studied recently based on a Gubser expansion. The evolution of spin potential turned out to be important in the future studies of spin polarisation in these solutions [32]. There were several other important studies of dissipative relativistic hydrodynamics based on the perturbations or other modifications of the relativistic and boost-invariant, hence infinite Gubser flow. Instead of detailing these results, due to the limited focus of our work, let us quote some recent review papers on these topics. For a very recent and rather pedagogical introduction to the effective descriptions relevant for attractors, in particular hydrodynamical attractors in high-energy heavy ion physics, holography, and kinetic theory, together with highlights of some recent advances in dissipative hydrodynamics, we recommend ref. [33]. The status of dissipative relativistic hydrodynamics has been reviewed less recently but more extensively in refs. [5,34]. This latter review [5] discussed and presented also the topics of exact solutions of perfect fluid hydrodynamics, so it is relevant background for our current manuscript.
For the context of our exact dissipative solutions and their relation to perfect fluid solutions, it is useful to mention that some exact non-relativistic as well as relativistic perfect fluid solutions seem to form a basis or a background solution that are recovered by our results in the limit of vanishing dissipation. For example, our results for the special case of κ = κ 0 generalise the Zimányi-Bondorf-Garpman [9] spherically symmetric perfect fluid solutions for non-vanishing bulk viscosity and for a cancelling shear viscosity coefficient, as evident from the scale equation: substituting ζ 0 = 0 to Equations (42) and (43), one recovers the scale equation for R(t) obtained already in ref. [9].
The scaling function of the density as well as for the temperature profile are related with a matching initial condition, as shown in Equation (28). In ref. [11], this equation was also obtained before but for a perfect fluid and a corresponding perfect fluid interpretation of the constant of integration (denoted by C E in the present manuscript and by C φ ref. [11]). It is also inspiring to note that if this constant of integration is vanishing, C E = 0, we obtain an inverse relationship between the scaling functions of the density and the temperature profiles, V (s)T (s) = 1. This matching condition is not limited to the non-relativistic kinematic domain, as it is found to be valid also for relativistic perfect fluid solutions, as detailed in ref. [35].
It is also inspiring to note that the above-mentioned class of perfect fluid solutions has a straightforward generalisation to spheroidal and ellipsoidal flows. The first triaxial, finite fireball solution, as far as we know, has been found by De, Garpman, Sperber, Bondorf, and Zimányi [36]. This oblate, ellipsoidal perfect fluid solution has been generalised to arbitrarily positive definite temperature profiles in refs. [12,37]. This suggests that a generalisation of the spherically symmetric, viscous solutions to oblate shaped, ellipsoidal fireballs may become feasible.
Furthermore, these expanding spheroidal or triaxial ellipsoids may not only expand and cool but rotate, too. Such exact and rotating fireball solutions have been obtained for non-relativistic perfect fluids in refs. [30,[38][39][40]. Due to the limited scope of this manuscript, we do not provide a complete list of references on exact rotating solutions of fireball hydrodynamics, as the already quoted papers, in particular ref. [40], seem to be sufficient to conjecture that viscous corrections of the non-relativistic Navier-Stokes equations can be obtained even for rotating, triaxially expanding ellipsoidal fireballs, too.
Thus, there seems to be a very interesting and deep connection between perfect fluid solutions and dissipative solutions of the Navier-Stokes equations. In the next section, we detail their relationship for the spherically symmetric, non-relativistic flows, which is the main topic of the current manuscript.

Asymptotically Perfect Fluid Behaviour
During the understanding of these new solutions, we have found that at late times, an expanding fireball with Hubble-type velocity profile is proceeding toward "perfection", or in other words, these new exact solutions of dissipative hydrodynamics are asymptotically equal to an exact solution of a perfect fluid hydrodynamics in the following well-defined sense: two functions of time t, called x(t) and y(t) are asymptotically equal, if lim t→∞ x(t) y(t) = 1.
Such an asymptotic equality is denoted by x ∼ y throughout this manuscript. The asymptotic equality of a dissipative solution of hydrodynamics with a perfect fluid solution implies that the effect of bulk viscosity can be scaled out at late times, corresponding to low temperatures. This result supports the conclusion of [17], which was obtained for relativistic kinematics, for a similar set of equations of state. In other words, this result implies that we cannot decide from final state measurements that the medium evolved as a perfect fluid with higher initial temperature and entropy content or as a viscous fluid with lower initial temperature. This behaviour is also perceived in ref. [16], which as far as we know presents the first exact solutions of the relativistic Navier-Stokes equations for non-zero bulk viscosities.
A similar but different phenomenon is reported in ref. [41], which presents an analytic solution of a set of differential equations that describe the transition from kinetic theory to hydrodynamics in the Bjorken expansion. Another interesting feature is provided by ref. [42]: for large times, a unique, non-singular, asymptotic solution of the non-linear Einstein equations in the bulk is found to be a perfect fluid solution of relativistic hydrodynamics. The conclusion of refs. [41,42] also highlight the importance of the identification of attractor solutions in dynamical equations and the theory of asymptotic series.
In the following subsection, we show how the approach to "perfection" can be described with analytic tools in the simplest case. Subsequently, we discuss this behaviour in more general cases and illustrate such a behaviour in two figures as well.

Asymptotic Analysis for the Case of a Homogeneous Pressure and Linear Bulk Viscosity
We have presented a new family of solutions with homogeneous pressure and arbitrary, pressure-dependent bulk viscosity in Section 6.1. At the end of that section, we found an analytic form for the time-dependent part of the temperature profile if we considered the case of ζ ∝ p: namely, a bulk viscosity coefficient that is a linear function of the pressure: where ζ 0 = ζ(p 0 ). In this case, the analytic form of the temperature, the conserved charge density, and the pressure have been obtained in Equations (37)- (39). For late times, t t 0 , the exponential terms in the T(t, s) temperature and the p(t) pressure are approaching a constant value. Thus, these fields are asymptotically are given as follows: where we stress that the ∼ sign stands for the asymptotic equality. The initial temperature and initial pressure are rescaled by the asymptotic dissipative correction and absorbed by the asymptotic perfect fluid initial temperature (T A 0 ) and pressure (p A 0 ): The effect of bulk viscosity is also absorbed into these parameters, so the bulk viscosity effect can be absorbed into the normalisation of the initial pressure and temperature, and it has no other asymptotic influence on the time evolution of the fireball. Accordingly, this asymptotic limit corresponds to a perfect fluid fireball hydrodynamics with T(t 0 ) = T A 0 and p(t 0 ) = p A 0 initial conditions. This family of perfect fluid solutions is already known and published in ref. [30]. However, that reference discusses a wider family of perfect fluid solutions: the spherically symmetric fireballs are generalised to spheroidally symmetric, rotating fireballs.
An interesting property of the kinematic bulk viscosity (ζ 0 /n 0 ) emerges from Equation (48). Using the p 0 = n 0 T 0 ideal gas approach, one can realise that the initial value of the kinematic bulk viscosity is a non-monotonic function of T 0 : and this function is illustrated in Figure 2, where ζ 0 /n 0 is divided by the unobservable initial time t 0 . Equation (48) sets an upper limit on ζ 0 /n 0 , and it is written as: where e denotes the Euler number. T A is an observable quantity by describing the experimental data with the perfect fluid solution of ref. [30]. For a fixed asymptotic solution with fixed T A 0 , the initial value of the kinematic bulk viscosity is a function of T 0 , and in this figure, the initial time is scaled out. This non-monotonic behaviour is described by Equation (48), and the maximum of the curve is given by Equation (51).

Asymptotic Analysis for an Inhomogeneous Pressure and Linear Bulk Viscosity
Now, let us examine the asymptotic behaviour of the more general solution presented in Section 6.2. It is a parametric solution, described by two coupled differential equations, namely Equations (42) and (43). We have numerically solved this set of differential equations for a perfect fluid and also for viscous fluids with different values of ζ 0 /p 0 . We obtained the R(t) scale of the fireball and the g T (t) function of the dissipative correction of temperature. These calculations are performed twice by two different approaches. The first approach is the conventional one, when we start the hydro evolution from the same initial conditions and vary the amount of bulk viscosity. The results of these calculations are shown in Figure 3, where one can see that the heating effect of the entropy production and that the same initial conditions lead to different asymptotic final states. This is a conventional feature of our solution. The second approach is when the initial conditions are varied, but the asymptotic equality of the solutions with different bulk viscosity coefficients are required. This is a non-conventional approach. We illustrated these cases in Figure 4. Looking at this figure, one can realise that the black curves, related to perfect fluid solutions, behave as asymptotic attractors. Thus, this figure provides the same conclusion that we have shown in the previous subsection for the case of homogeneous pressure with analytical tools: the exact solutions of dissipative fireball hydrodynamics are asymptotically equal to a perfect fluid solution with increased initial temperature and entropy content. If the asymptotic perfect fluid solution is fixed, then varying the bulk viscosity coefficients can be co-varied with the initial conditions, so that asymptotically, the same perfect fluid solution governs or attracts the late time behaviour of these dissipative solutions. In Figure 4, we characterised the curves with the particle mass m, the energy density to pressure ratio κ that characterises the equations of state, and the initial properties of the asymptotic perfect fluid attractor. We find that only the following physical quantities are important in determining the physical characteristics of the perfect fluid asymptotic solution: the initial temperature of the asymptotic perfect fluid solution T A 0 , the initial value of the scale R A 0 , and the initial velocity of the scaleṘ A 0 . These quantities, together with the initial density n 0 and the particle mass m, determine a perfect fluid solution that becomes asymptotically equal with the dissipative exact solution of the non-relativistic Navier-Stokes equations. Scheme 4 of Ref. [1] is replaced with the following corrected Figure 4. Note, that the caption of the new Figure 3 is unchanged, it is the same as that of Ref. [1]. Fortunately, the captions of both Figures 3 and 4 as well as the rest of the manuscript, including the conclusion, remain unchanged.
The authors would like to apologize for any inconvenience caused to the readers by these changes.
Author Contributions: Writing-original draft, G.K.; Writing-review and editing, G.K., T.C. and L.P.C. All authors have read and agreed to the published version of the manuscript. Scheme 4 of Ref. [1] is replaced with the following corrected Figure 4. Note, that the caption of the new Figure 3 is unchanged, it is the same as that of Ref. [1]. Fortunately, the captions of both Figures 3 and 4 as well as the rest of the manuscript, including the conclusion, remain unchanged.
The authors would like to apologize for any inconvenience caused to the readers by these changes.  R A 0 = 0 initial parameters, utilising an m = 940 MeV for the particle mass and a constant, temperature-independent κ = 3. The solid black line stands for a perfect fluid solution, and this perfect fluid curve labelled by zero bulk viscosity is approached by each of the shown exact viscous solutions asymptotically, T(t)~T A (t). The dashed blue, the dotted-dashed green, the dashed yellow, and the dotted-dashed red lines correspond to our new viscous solution of non-relativistic Navier-Stokes equations for different values of ζ 0 /p 0 , but for the same asymptotic solutions.

Summary
We have found a new family of analytic, exact solutions of non-relativistic Navier-Stokes hydrodynamics in 1 + 3 dimensions. In these solutions, the velocity field is spherically symmetric Hubble flow, so the effect of shear viscosity cancels, and we assumed that the heat conduction is negligible. With these features, we provide exact and analytic solutions for the cases of homogeneous and inhomogeneous pressure as well. Let us comment that even ignoring heat flow and assuming that the temperature and pressure are homogeneous is still appropriate for analysis of non-relativistic, dissipative hydrodynamics, and even for the analysis of non-relativistic heavy ion collisions, as shown for example in refs. [9,18]. The solutions for homogeneous pressures correspond to the non-relativistic limit of a recently found exact solution of the relativistic Navier-Stokes equation [17]. This result strengthens the conclusion of refs. [16,17]. In this class of exact solutions of dissipative hydrodynamics, for asymptotically late times, the effects of bulk and shear viscosity coefficients can be scaled out, and asymptotically, the same hadronic final state or attractor solution can be reached using different bulk and shear viscosity coefficients, with suitably co-varied initial conditions. These solutions share a common perfect fluid asymptotic solution.
The non-relativistic kinematic domain may make the significance of this paper rather academic. On one hand, the assumption of spherical symmetry, the absence of heat flow, and assuming homogeneous pressure and temperature field may obliterate the physical properties of dissipative fluids. On the other hand, we have found the effect of approach to "perfection" first in the relativistic kinematic domain [16,17], before writing it up in the current, clear-cut non-relativistic form. Although the domain of applicability is limited to spherical symmetry in the present manuscript, we conjecture that a similar process of asymptotic "perfection" of dissipative hydrodynamical solutions is valid also in spheroidal and ellipsoidal, as well as rotating, parametric solutions of the Navier-Stokes equations. Thus, we conjecture here that the phenomena of approach to "perfection" is a rather general phenomena. In the present manuscript, we have proven that the phenomena of approach to "perfection", namely the asymptotic equality of a dissipative hydrodynamical solution with a perfect fluid solution exists, and its domain of validity includes the non-relativistic, spherical fireballs that are governed by the Navier-Stokes equations. This domain of validity of the asymptotically perfect fluidity of dissipative flows is to be studied and extended to other kinds of fireball hydrodynamics in future, subsequent studies.
Thus, we have shown that the approach to "perfection" is not specific to relativistic kinematics [16,17] but seems to be a more general phenomenon, valid for certain exact, 1 + 3 dimensional, spherically symmetric, parameteric solutions of the Navier-Stokes equations. Further generalisations of these results are in progress but will be discussed in separate manuscripts.