Density Profiles of Gases and Fluids in Gravitational Potentials from a Generalization of Hydrostatic Equilibrium

It is well-known that equilibrium density profiles of gases and fluids in gravitational potentials have an r −1 dependence, where r is the radius. This is found in both astronomical observations and detailed simulations in spherical-ly-symmetric geometries. It is also well-known that the standard equation for hydrostatic equilibrium does not produce such solutions. This paper utilizes a Lagrangian formulation that produces a closed-form r −1 solution and identifies the needed terms to be added to the standard equation for hydrostatic equilibrium to obtain such a solution. Variants of the r −1 solution avoid a density singularity at the origin and a total-enclosed mass singularity at infinity. The resulting solutions are shown to be in good agreement with well-established density profiles of ordinary matter in galaxies, dark matter haloes, and the atmosphere of earth. Comparisons are made to solutions based on the standard hydrostatic equation, including solutions based on the Lane-Emden equation. The origins of differences are explained.


Introduction
Early numerical solutions for dark matter (DM) density profiles were found to have dependence inversely proportional to radius [1]. Astronomical observations of both ordinary matter and dark matter have also inferred density profiles closely resembling the r −1 profile, where r is the radius from the nominal center of the distribution, in a region surrounding the origin [2] [3] [4] [5] [6]. However, the profile has several serious drawbacks. First, it gives infinite density at the origin. Second, it gives infinite total enclosed mass. Third, it is not consistent with many observations of the shape of the inner profiles of low-spectral brightness galaxies [7] [8] [9]. Nonetheless, the r −1 profile is "ubiquitous" [10]. This paper will provide theoretical support that such a profile is indeed ubiquitous when there is thermal and mechanical equilibrium.
To analyze the spatial distribution of a gas of particles, N-body simulations [1] [ 11] or solutions to the Vlasov equation [12] are typically used. Two simplified governing equations are considered for DM density profiles in this paper. The first is the standard equation for hydrostatic equilibrium in a spherically-symmetric geometry: here ρ is the number density of DM, P is the pressure, r is the radius, m υ is the mass of a DM particle, ( ) enc M r is the enclosed mass, and G is the gravitational constant. This equation can be solved using the well-known Lane-Emden formulation [13], which uses the assumption that pressure is a function of density only: where c γ is a constant for a given polytropic exponent γ . An inhomogeneous form of the Lane-Emden equation can be used when ordinary matter (OM) is present as shown later in the paper. As shown in Section 4, calculations using Equations (1) and (2) [16]. This is found to be true for any polytropic exponent between about 1.1 and 2 in Section 4 below.
To address this, a generalization of Equation (1) is derived using a Lagrangian formulation. This derivation is given in Section 2. A partial validation of this derivation is given in Section 3 using the properties of a well-known density profile of a gas in a gravitational potential, namely, the earth's atmosphere. Section 4 presents behavior of solutions near the origin. Section 5 presents the behavior of solutions far from the origin and presents a sample calculation for a dark-matter halo. Section 6 compares the solutions to Einasto and de-projected Sersic profiles. Section 7 compares the solutions to results from the Lane Emden equation.
Section 8 provides a brief summary.

Derivation of a Generalized Equation for Hydrostatic Equilibrium
The assumed kinetic and potential energy densities for a medium in a spherically symmetric geometry under the influence of gravity is given by . . , . .
where the variables are defined as in Equation (1). In the topic of this paper, m rad r ρ . Note that the pressure is treated as a kinetic energy for the obvious reason. This will allow a variational approach with a Lagrangian density. Explicit relativistic effects and the possible impact of angular momentum are neglected in this initial analysis The pressure is assumed to be a function of density only. As an example, the number density for fermions in Equation (3) assuming thermodynamic equilibrium is where p is the fermion momentum, v m is its mass, s n is the number of spin states, v µ is the chemical potential, and ћ is Planck's constant. The other variables have been defined earlier. In Equation (4), the chemical potential is set to zero for non-interacting particles. With Equation (4) the pressure for fermions has the usual form in the non-relativistic limit, and is given by This equation should hold for arbitrary small variations To address the last term, assume Equation (2) The integral in the last term in Equation (8) can then be eliminated using Equation (7) to give The second term in Equation (9) is the change in the potential energy of the system at radius r due to mass at radius r, and can be dropped if the enclosed mass is understood to be the enclosed mass strictly less than r.
To address the term involving  the second line in Equation (9) can be integrated by parts, accounting for the volume integral of the Lagrangian density, as well as the r 2 factor that is associated with it. After considerable algebra, the result is This equation now has the required form so that the expression in curly brackets is zero: Equation (11) is the result for strict equilibrium for a medium with a pressure that depends on density only in a spherically-symmetric geometry, and with a single species of matter. The limiting case of the standard hydrostatic equation is obtained by neglect of the second terms on both sides of Equation (11).
Based on the usual derivation of Equation (1), the standard hydrostatic equation should be valid when mechanical equilibrium applies but the stricter action equilibrium of Equation (11) does not apply. The term involving ( )( ) It is also worth noting that Equation (11) can be written as This simplification indicates a solution for density of the form This result is both remarkable and perhaps expected. It arises because of the spherical geometry and the r −1 potential as discussed above. It matches the Navarro-Frenk-White (NFW) result [1] as well as other related results near the origin [2] [3] [5]. However, it has two serious drawbacks. First, it gives infinite density at the origin. Second, it gives infinite total enclosed mass. One may attempt to address these issues within the Lagrangian formalism by introducing a density constraint using a Lagrange multiplier ρ λ , and a total particle number (mass) constraint with multiplier M λ . With these two constraints, Equation (10) be- The last term arises from the variation of the total-mass constraint with respect to density. For r small (close to the origin) or large (far from the origin), one might expect this integral to be approximately zero, since the total mass constraint requires that the radial integral from zero to infinity of the density variations be identically zero. This version of the equation will be explored in subsequent sections of this paper.

Generalized Equation for Earth's Atmospheric Density Profile with Altitude
As a basic check of the results of Section 2, one may consider the case of earth's atmosphere. In the earth's atmosphere, it is well-known that (a) the temperature is not constant with altitude in the troposphere [17], and (b), that despite (a), a density profile computed with an isothermal atmosphere is a decent approximation. This apparent inconsistency is investigated in this section with various solutions and in the process provides some confirmation that the generalized hydrostatic equation has validity. First the version of the generalized equation for hydrostatic equilibrium given by Equation (11) where T(z) is the temperature, which is treated as an input from [17] This equation simplifies to (neglecting terms of order z/r e ), This can be solved analytically: here 0 ρ is the density at sea level, which is an input. The solution of Equation (19) gives poor agreement with measured data, with significantly larger density than measured at higher altitudes. This will be seen to arise from the underlying assumption of strict action (mechanical and thermal) equilibrium. Equation (19) may be compared to the isothermal approximation, which amounts to retention of only the first terms on both sides of Equation (17) and assuming constant temperature. One obtains the usual exponential atmosphere in this case, where T 0 is the temperature of the atmosphere at sea level.
There is also the adiabatic equation for pressure and density. The above Equations (17)- (19) assumes strict action equilibrium. The adiabatic equation assumes that hot parcels of air near the surface of the earth rise without significant heat exchange with the surrounding air. Hence it corresponds to a process in quasi-static equilibrium that is not in strict mechanical and thermal equilibrium at all layers. With this insight, the adiabatic result can be derived from Equation (17) by neglect of the last terms on both sides of the equation, and assuming the gas satisfies a polytropic relation (2) with a polytropic exponent ~1.4 γ . With these assumptions, one obtains the result from the standard hydrostatic equation, which is the well-known formula for density [18]: here c p is the specific heat at constant pressure of the atmosphere at sea level.
Note that K. The latter is consistent with the ideal gas law and the sea-level values of pressure and density given in [17]. For Figure 1, a polytropic exponent of 1.38 rather than 1.40 is used and is justifiable because of the presence of significant water vapor in the middle atmosphere between 1 and 15 km. The integration step size is 250 m. It should be noted that Equation (21)  To remedy the errors above 15 km, the additional terms in Equation (15) can be restored gradually above 15 km, representing the return to true equilibrium at higher altitudes. The altitude dependent factor used to restore the full Equation (15) is given by Equation (22): where z is in km in this expression and H(z) is the Heaviside function. This factor gives a gaussian onset to equilibrium and multiplies the terms in Equation (15). A numerical integration of (15) with this factor is shown in Figure 1, along with the other solutions given above. All are compared to the standard (measured) atmospheric pressure and density from [17].
Comparing the cyan and blue curves for both pressure and density in Fig

Solution near the Origin
As discussed at the end of Section 1, the unconstrained equilibrium solution gives infinite density at the origin. This defies the Fermi-Dirac equation, which implies a maximum density due to fermion degeneracy pressure. It also seems counter-intuitive, at least in the absence of a gravitational singularity. To remedy this, a density constraint (as well as a total particle constraint) is added to the Lagrangian density, with a result given by Equation (14). It is not difficult to see that Equation (14) can be written The last term of Equation (23), comprising the last two lines, can be estimated in several ways. The simplest approach is to assume that the density variations in the integral average to zero, yielding zero for this term. However, a more careful approach is appropriate, and this can be done assuming the zeroth-order solution, Equation (13). It is further assumed is much less than 1, to investigate behavior near the origin. Then the various factors in the last term of Equation (23) are The "big-oh" and "little-oh" notations are used in Equations (24b) and (24c). With Equations (24), the last term of Equation (23) can then be estimated as Note also that with constrained total particle number, the variations should which allows Equation (25) to be re-written as With bounded variations ( ) r δρ ′ in 0 to r, this expression tends to zero as r tends to zero, when the denominator is not zero (as should be the case).
Hence, however the last term is estimated, the summary result for Equation (23) near the origin is Note that this has two possible solutions for any given r. The first solution is obtained by the usual approach, setting the quantity inside curly brackets to zero. This gives the r −1 solution presented in Equation (13) This latter solution is consistent with an incompressible medium, as occurs when the density is limited by Fermi-Dirac statistics. It is also consistent with the observed features of most celestial bodies, in which the density is approximately constant at their core, and also often in spherical shells. It is also approximately consistent with observations for dark matter in central regions that are less than a kpc in diameter [7] [9]. It implicitly assumes that the properties of the constituent materials are constant within spherical regions.
It also is self-consistent with expression (27) equal to zero in such a region.
With these facts in mind, the solution (29) of constant density in a region nearest the origin is a candidate solution. It should be further noted that this same result of constant density can be applied far from the origin where particles achieve their nominal cosmic background temperature and density.
Next consider the standard hydrostatic equation. One can assume a solution of the form 1 0 0 r r ρ − for the density and substitute into Equation (1) with Equation (2) and check for consistency. One finds that ( )( ) One notes that the left-hand side varies as r γ − whereas the right-hand side is a constant under the assumption of   The solution is clearly approximate because for r greater than a r the density can be become complex and also because the density variation of DM in Equation (32) is not included in the enclosed mass in Equation (31). However, it does show that a constant density to order of ( ) 2 a r r is consistent with the standard hydrostatic equation.
To summarize Section 4, it is found that setting the variations of the Lagrangian density to zero yields an alternative, constant-density solution. This solution is consistent with the observed density profile of many celestial bodies near their origin, including observationally-inferred dark matter profiles.
The standard hydrostatic equation also shows solutions consistent with approximately constant density at the origin but is seen to be inconsistent with r −1 solutions. With these considerations in mind, the constant-density solution is used in a region of radius r 0 near the origin in the following sections.

Behavior of Solutions away from the Origin
The mass-constrained Equation (14) is most relevant here. The last term in the equation is assumed negligible in this case, since the integral of any physical density variations from r to infinity should tend to zero for large r in the case of constrained total particle number. The result is Inspection shows that this equation can be simplified considerably, giving, as in Equation (28), This result shows that the solution is independent of the density constraint ρ λ with the neglect of the last term of Equation (14). This equation may further be put into dimensionless form by defining , where is c ρ the cutoff density and c r is the cutoff radius where the density drops due to the constraint on the total particle count. The resulting equation is This equation may be solved numerically or perturbatively. For a perturbative treatment, assume ( )( ) c c r r ρ ρ is much less than 1, corresponding to radii near the origin but not at the origin. Then the zeroth order solution is given by Equation (12), where the "big-oh" notation is used to indicate the order of the approximation.
In the last expression, It should be noted that the results of Figure 2 are not explicitly dependent on the polytropic exponent, particle mass, or temperature. The result does depend on the latter 2 parameters implicitly through the initial density at or near the origin via Equation (4), for example. One may observe from Equation (35)  It is well known that logarithmic exponents less than −3 correspond to solutions that have convergent enclosed mass [1] [2].
In both this section and Section 4, the last term of Equation (14)

Comparison to Einasto and De-Projected Sersic Profiles
Taking all the above sections into consideration, one finds that is always less than 1 where there is appreciable density, so the resulting solution is This solution has a central core region that is qualitatively consistent with the size of inner profiles of low-spectral brightness galaxies [2] [8] [9]. Further, the constant density in a core region avoids a density singularity and so is consistent with a Fermi-Dirac density distribution at finite temperature. Moreover, the finite density at the origin is consistent with all observed astronomical bodies except black holes. Lastly, the "core" solution offers a potential path to resolution to the "core-cusp" problem in dwarf galaxies.
For the region just outside the core, the logarithmic slope given by Equation α β γ = model is given in the previous section. Equation (40) can also be compared to the well-known de-projected Sersic and Ei-nasto models [4] [5] [6]. The de-projected Sersic (dpS) density distribution for radiant matter is approximated by [2]: The parameter 0 rad ρ is obtained by setting the volume integral of Equation (41) equal to the measured or inferred ordinary mass of the galaxy. The variable e R denotes the radius which encloses 1/2 the total light of the galaxy. The dpS profile is also used with some success for characterization of DM density profiles versus radius. The other two parameters in Equation (41) are given conveniently and approximately from Equations (19) and (27)  (43) The parameter n ranges from roughly 2 to 4 for DM haloes of galaxies and clusters, from the same reference. With these values of n, the leading exponent n p in Equation (42) is −0.7 to −0.85, which is not too different from −1, and [2] notes that the logarithmic slope is almost always approximated by −1 in the vi- However, as with the dpS model, the onset of steeper logarithmic slopes, slopes less than −1 but greater than −3, occurs for smaller r than for Equation (40).
This different behavior of the model solutions can be addressed with the possible solutions or extensions of Equation (14). One approach is to add angular momentum. The net effect of angular momentum would be to limit the extent of the galactic halo and would result in a more rapid decline in density than given in Equation (40). However, faster declines are found in simulations without appeal to angular momentum. As a second approach, one may consider the results of Section 3, which indicate that if the medium is not in a strict equilibrium, then solutions with faster declines occur. The outer perimeters of a large halo are not expected to be in strict equilibrium because of the process of aggregation of subhaloes [10]. Proceeding as in Section 3, one may neglect the term  in the first line of Equation (14) Consider the region r which is significantly greater than 0 r but also signifi- For related states of matter, such as neutron stars, a polytropic exponent γ ranging from 3/2 to 2 is often used [20] [21] [22]. With these choices for γ , one obtains familiar forms of solutions, One may also compare the logarithmic derivative of Equation (48) The parameter 5 3 L is shown in Table 1 as a function of particle mass for temperatures of T = 100, 8, and 2 K respectively at the origin. These temperatures set the density 0 ρ at the origin via Equation (4) [26]. Masses of 20 eV/c 2 or more for these temperatures give 5 3 L of 3 kpc or less, corresponding to the diameter of dwarf galaxies. One may note that the temperatures and masses shown correspond to non-relativistic conditions.  Table 1 shows results for 5 3 r , assuming r 0 = 1 kpc. The table highlights the rapid variability of this scale parameter with mass. Masses of 1.2 to 1.6 eV/c 2 give scale sizes corresponding to the peak of the large-scale mass power spectrum, at a nominal background temperature of 2 K, and masses of 10 to 50 eV/c 2 yield scale sizes that correspond to radii of large to dwarf galaxies. Note that the masses shown in this table are not consistent with recent estimates of dark-matter particle mass [27] [28], which identify masses of the order of 5 keV. This implies that some underlying assumption of this paper is not consistent with the assumptions of those results. The only key assumptions of Table 1 are (a) mechanical or energy equilibrium applies, (b) the particles are fermions, and (c) ordinary gravity applies. The value of 5 3 r also depends on r 0 , the region in which the density is approximately constant and equal to 0 ρ . This region is a corollary of the assumption of fermionic matter in the absence of a gravitational singularity. Figure 3 shows some density profiles corresponding to Table 1. Plot (a) of Figure 3 shows results for the Lane-Emden equation, Equation (56). Plot (b) of Figure 3 shows results for the generalized hydrostatic equation, Equation (40). The assumed conditions are shown in Table 2. For the numerical integration, 4000 radial steps are used, with each step equal to 77 kpc. The seed mass at r = 0 is set to the estimated mass of OM of a galaxy similar to the Milky Way, which is chosen to be 9 × 10 10 solar masses [19]. The chemical potential is set to zero for the results of Figure 3 because the particles are assumed to be free and non-interacting rather than bound (however, the results shown in plot (b) of Figure 3 are also consistent with bound states of matter in fluid or gaseous form). The temperature T at the origin for plot (b) of Figure 3 is set to 25 K in order to remain above a number density floor of 10 8 m −3 out to a radius r c of 300 Mpc. This choice of number density floor is based on the estimated neutrino density in standard cosmology [29]. The average temperature of the profiles on the right is about 0.15 K in all cases, using Equation (4) to relate density to temperature. Note also that plot (b) of Figure 3 derives the r −1 solution. For this solution, the scale parameter r γ of Equation (50) and Table 1 is not directly relevant.  In summary, this Section compares the results of the Lane-Emden equation with that of the generalized hydrostatic equation. It is seen that the former produces a different dependence of scale size on temperature and particle mass that does the generalized hydrostatic equation, and the shapes of the solutions are significantly different near the origin; the former is flat, the latter has a cusp at or near the origin. However, in both cases, a narrow range of masses and temperatures yield scale sizes that are consistent with observed galactic and large-scale structure as seen from Table 1.

Summary
In summary, Section 2 shows a derivation of a generalized hydrostatic equation in spherically symmetric geometries with a single species of matter satisfying a polytropic relation. Section 3 provides a partial validation of the generalized hydrostatic Equation (10) using the earth's atmosphere. Section 4 shows that there is a constant-density solution for the Lagrangian formulation in a spherical region about the origin (or in spherical regions about the origin), and that the r −1 solution near the origin is not self-consistent with the standard hydrostatic equation, as expected, when the added terms of the generalized equation are not included. Section 5 explains how Equation (14) gives a convergent enclosed mass far from the origin and shows an example solution of the differential equation.
Section 6 demonstrates that the solutions can match the radial behavior of accepted models of OM and DM density versus radius away from the origin. Section 7 presents solutions of the generalized equation for the density and compares it to the standard hydrostatic equation with a polytropic exponent in a gravitational potential. The overall conclusion is that a generalized equation for hydrostatic equilibrium, Equation (40), gives a cuspy solution for DM that is more consistent with observations and simulations (e.g., [1] [2]) than the standard hydrostatic equation near the origin and more consistent with finite density and lower slope in the immediate vicinity of the origin [7] [9] than the de-projected Sersic or Einasto models. Because this paper shows that the r −1 den-