On the stability of fermionic non-isothermal dark matter halos

The stability of isothermal dark matter halos has been widely studied before. In this paper, we investigate the stability of non-isothermal fermionic dark matter halos. We show that in the presence of temperature gradient, the force due to the pressure has both inward and outward components. In some regions of halos, the inward force that provides stability is due to the pressure rather than gravity. Moreover, it is shown that higher temperature gradients lead to halos with lower mass and size. We prove that if the temperature is left as a free positive profile, one can place no phase-space lower bound on the mass of dark matter. For halos that are in the low degeneracy classic domain, we derive an analytic expression of their temperature in terms of their mass density and place an upper bound on the mass of dark matter by requiring that temperature is not negative. We then use the Burkert mass profile for the Milky Way to show that if the central temperature of the halo is a few Kelvins, the mass of dark matter cannot exceed a few keV.


Introduction
A variety of independent observations leave no doubt on the existence of dark matter (DM). These include the early measurements of galaxies velocity dispersion in the Coma cluster [1], the rotation curves in galaxies [2], the recent measurements of the gravitational lensing [3,4], the Bullet cluster [5], the anisotropies in the CMB [6], and the large scale structures [7].
Since no candidate particle has been observed in any of the experimental searches [8][9][10][11], direct measurement of DM properties such as its mass is not possible at present. However, galaxies are giant DM laboratories where we can investigate the properties of the constituting particles by relating them to the observable features of DM halos such as their stability. The involved forces are gravity which can be explored a e-mail: ahmad_borzou@baylor.edu (corresponding author) by observation and the gradient of the pressure of DM. The latter depends on the equation of state (EOS) of DM and as a result, provides an invaluable possibility to investigate the foundations of these long-sought particles.
Since the visible mass of galaxies is made of fermions, DM may also have fermionic nature. In this case, the pressure of the halo is described by the Fermi-Dirac statistics. This scenario has been studied before in [12][13][14] for isothermal halos. Recently, degenerate models of fermionic DM halo have attracted attention mainly due to their potentials for addressing the core-cusp problem [15][16][17][18].
In the Fermi-Dirac statistic, the pressure is a function of the temperature and the fugacity-or equivalently the chemical potential. Therefore, the force due to the pressure has in general two components. The first force of the pressure is due to the derivative of the fugacity, which is investigated in references above within the isothermal halos. The second force of the pressure is due to the temperature gradient. As we will see in this paper, one of the two forces of the pressure can be inward, adding to the pulling gravitational force and deepening the potential well.
In general, fermions can be compressed in arbitrarily small volumes despite the restriction on their phase-space. Black holes are the living examples of compressed fermionic systems. Even though the phase-space of fermions is limited, there is an unlimited momentum space available to them. By reducing their configuration volume, they will occupy higher momentum states to meet the Pauli exclusion principle. It is only the balance of the forces due to their pressure and gravity that determines the size of fermionic DM halos. As we will see in this paper, the inward force of the pressure can be the dominant pulling force in a vast region of the halo that confronts the pushing force of the pressure and maintains stability.
In this paper, we discuss that a constant temperature across DM halos is not validated through observations. Moreover, even for collision-less DM models, there are a variety of heat generation mechanisms that can lead to at least slight devi-ations from isothermal models. We review the gravitational and non-gravitational frictions and gravitational contraction as the main sources of heat generation, and radiation and convection as the means of heat transfer in DM halos. The frictional effects may play crucial roles in satellite galaxies that move with relatively high speeds through their host galaxies. Also, the gravitational contraction mechanism of heat generation is more important in these compact halos than in larger dilute halos. Therefore, as we will show in this paper, it is possible that the compactness of satellite galaxies is due to their different temperature profiles.
We show that if DM distribution in the halos is Maxwellian, i.e. fermions are entirely at low degeneracy level, the temperature profile can be expressed analytically in terms of the mass profile. By requiring non-negative temperatures, we put an upper bound on the mass of DM. Also, we analytically prove that the temperature profile is irrelevant if DM halo is entirely in the infinite degeneracy level.
To study non-isothremal fermionic DM halos, we derive the most general hydrostatic equilibrium equation for a spherical Fermi-Dirac system. Computer software is developed to solve the field equations numerically. We validate the software by reproducing the known solutions of different degeneracy levels. Since the software is meant to be general, rather than assuming a specific DM and galactic model to derive the temperature as a function of the radius, we reserve the assumption for directly inserting a temperature profile into the software.
We study non-isothermal models with a generic temperature profile of the form T = T 0 1 + r , in a range of low central temperature of T 0 ≤ 1 K. We show that the isothermal solution r 0 = ∞ leads to the largest halo size and mass. By increasing the temperature gradient, the mass and the size of the halo decrease. We present solutions that are more compressed than their corresponding infinitely degenerate halos.
The current lower-bound on the mass of DM is derived using highly compressed infinitely degenerate isothermal solutions [18]. However, in the presence of temperature gradient, more compressed halos are possible; allowing for the possibility of lighter DM masses. We discuss that since (i) the pressure and the mass density of Fermi-Dirac systems are functions of two independent profiles of temperature and fugacity, and (ii) there is only one stability equation in terms of the mass density and pressure, any arbitrarily light dark matter mass can explain any observationally supported mass profile if the temperature is left as a free positive profile. Therefore, phase-space lower bound on the mass of DM requires modeling the non-universal temperature profile of DM halos. This paper is organized as follows. In Sect. 2, the statistics of Fermi-Dirac systems and their stability criteria are reviewed. In Sect. 3, we derive the most general stability equation for halos made of fermions and present computer software to solve it. In the same section, we study a class of solutions whose temperatures decrease with distance from the center. Also, we place an upper bound on the mass of DM using the stability of low degenerate halos. In Sect. 4, we discuss the phase-space lower bounds on the mass of DM and show the relevance of temperature profiles for them. A conclusion will be drawn in Sect. 5.

An overview of the Fermi-Dirac statistics and non-isothermal halo stability
The visible matter in galaxies is made of fermions. It is quite likely that DM also has the same nature. In general, both temperature and mass density are functions of the radial distance from the center of halos. Therefore, the degeneracy of fermionic DM can be at the opposite extreme levels at two locations of one same halo. Therefore, the most general EOS of the fermionic matter and its corresponding stability equation is of interest. In this section, we review the statistics of non-isothermal fermionic systems and the corresponding stability equation. We also briefly discuss the sources that can establish temperature gradients in DM halos.

EOS of Fermi-Dirac statistics
Since the mass density, and the pressure are local quantities measured in the free-falling frame at a given distance r , it is easiest to derive the EOS of fermionic DM from its energymomentum tensor in the same frame. Calculations in this frame are particularly advantageous if the corrections to the EOS due to possible interactions between DM are of interest. The energy density and the pressure operators of the system are the temporal and the spatial components of the energy-momentum tensor associated with the Lagrangian of a non-interactive Dirac spinor field where the integral is over microscopic scales and V is a local volume.
In a free-falling frame, the gravitational effects are absent, and the metric is Minkowskian. Working in this frame, we can insert the free field expansion of the fermionic field into equation above to derive the energy density and pressure operators. The pressure, the energy density, and the mass density are the ensemble averages of the corresponding quantum operators, which using the non-relativistic approximation, read P = 2(kT ) 5 2 where z is the fugacity, β = 1 kT , and α ≡ h √ 2π m in terms of the Planck constant. Also, the Fermi-Dirac integrals are defined as with (ν) being the gamma function. Moreover, it is useful to know that the derivative of the Fermi-Dirac integrals read The most general EOS for non-interacting fermionic DM can be easily read from Eq. (2) where When the degeneracy is high, i.e. z 1, we can use the Sommerfeld approximation In the case of very low degeneracy level, we can use the following approximation when h(z) 1, and we recover the EOS of the classic ideal gas.

Conservation of energy-momentum, stability equation, and temperature profile
In a stable solution, the net force on an arbitrarily small volume of mass is zero at any point of the halo. The hydrostatic equilibrium equation can be systematically derived using the conservation of the energy-momentum tensor together with Newton's field equation and reads where G is Newton's constant, and M(r ) is the mass enclosed within the distance r.
If DM annihilation and creation are not significant, the conservation of the energy-momentum also implies that This and Eq. (9) can be combined into a second order differential equation where the pressure and density are local, i.e., measured in the free falling frame at distance r . Also, the equation is valid regardless of the statistics of DM. Finally, we would like to mention that the equilibrium equation above can be equivalently derived in a frame attached to the center of halo, rather than to a free-falling observer. In this frame, the energy of DM particles receives an extra gravitational potential energy which should be absorbed by the chemical potential, i.e. the latter is also different in the two frames. The disadvantages of this frame are that the equations contain both astronomical and microscopic lengths, the equations are not manifestly covariant anymore, and in the case of strong gravitational effects, cumbersome calculations are needed to arrive at Tolman-Oppenheimer-Volkoff stability equation.

Sources of heat generation and heat transfer
Since the stability equation depends also on the gradient of temperature, it is important to model the temperature profile of DM halos. In this section, we show that unlike the mass density with the equation for hydrostatic equilibrium, the temperature profile does not have a unique differential equation. It not only depends on the model of DM but also is a function of galactic properties such as relative speeds, the position of other galaxies, and the mass densities. We specifically discuss the means of heat transfer like radiation and convection and the origins of the heat, such as gravitational contraction, tidal forces, and friction. Heat generation Any phenomenon through which the macroscopic kinetic energy of galaxies is transformed, due to the conservation law, to the microscopic kinetic energy of the constituting particles, falls under the category of friction. The drag forces felt by compact dwarf galaxies that are moving through their hosts may set up considerable temperature gradients. Such forces usually depend on the relative speed, the mass densities, and other environmental properties. In collision-less models of DM, the Chandrasekhar friction is the most significant drag force. In interacting models of DM, a variety of frictional forces are possible. For QED-like interactions, the drag force can be approximated by the familiar form of the friction in the fluid dynamics.
Conversion of the gravitational potential energy to the kinetic energy of DM is another means for heating galaxies. The well-known Kelvin-Helmholtz mechanism belongs to this category.
Tidal stirring is another means of heating a halo. Most of the satellite galaxies orbiting their host experience a tidal force due to the difference between the host's gravitational force at the front and rear edges of the satellite. Heat transfer In collision-less DM models, only gravitational sources are available for radiation. The Kelvin-Helmholtz mechanism is an example. If DM interacts through a longrange force, in SIDM category, for example, DM can cool down through dark radiations, for instance see [19][20][21]. The Eddington equation for temperature gradient due to radiation reads [22] dT dr where κ is a constant and L is the luminosity of the radiation at distance r from the center. Convection is the most efficient means of heat transfer and is available to both collision-less and interacting DM models. It can be shown that the condition for the convection in a Fermi system reads where the right hand side is the adiabatic temperature gradient for a general fermionic gas. If the heat transfer is fast enough and the galaxy under study is convective, one can take the equality in this equation to determine the temperature gradient.

Stability solutions of fermionic DM halo
In this section, we would like to study the stability of DM halo under any degeneracy level. We first discuss the extreme degeneracy scenarios analytically and then present computer software for the sake of investigating the degeneracy regions that cannot be explored analytically.

Temperature profile of a DM halo with Maxwellian distribution
In a wide range of models, DM particles are assumed to have Maxwellian distribution. Therefore, their halos are described by the EOS of a classic ideal gas which reads where naught refers to the values at the center, m is the mass of DM, k is the Boltzmann constant, and T ≡ T 0 y is the temperature at arbitrary distance from the center. This EOS combined with Eq. (9) implies that If the temperature is constant everywhere in the halo, yρ = ρ, and a direct substitution confirms the stability of the wellknown isothermal solution This solution is singular at the center, i.e. does not satisfy the initial conditions, and also is in contradiction with observations of mass density at around the center. The temperature profile can be derived through the integration of differential equation above where y 0 = 1 is used. Since temperature cannot be negative at any distance from the center, we can derive an upper bound on the mass of dark matter Since possible dark matter interactions are not strong, T 0 is not expected to be high, and we should be able to find a fair maximum value for it. Therefore, the upper mass for dark matter can be set by observing the mass profile of different halo types.
As an example, in [23], the preferred mass profile for the halo of the Milky Way is found to be the Burkert profile with ρ 0 2.7 × 10 −21 (kg m −3 ) and r 0 10 (kpc). Therefore, r 2 dr rise with the distance until reaches a flat plateau of ∼ 0.58 kg 2 m −4 -hence the maximum value, at around 30 (kpc). Inserting these into Eq. (18), we find that If T 0 is around a few Kelvin, the mass of DM cannot exceed a few keV. This can be converted to a lower bound for the dispersion velocity at the center of the halo The predicted temperature profile of the Milky Way with the largest m T 0 is plotted in Fig. 1. Finally, we would like to emphasize the role of visible matter that has been neglected in the preceding analysis. In the central regions of most of the galaxies, the visible matter can have a higher mass density than DM, see for example the simulations in [24,25]. Therefore, the gravitational force due to visible matter in the interior of galaxies is greater than the force due to DM. Consequently, M(r ) in Eq. (9) should be replaced by M(r ) + M * (r ) where asterisk refers to the visible matter. This means that the force of pressure of DM should now confront the force of gravity due to both DM and visible matter. Therefore, the upper bound on the mass of DM in Eq. (18) should be modified by M(r ) → M(r ) + M * (r ) in the denominator. Due to the extra gravitational force, DM should be even lighter than when the visible matter is neglected.

Temperature profile of rotating halos
In general, the stability equation for a halo rotating with constant angular velocity ω reads where x 1 is the radial coordinate in the plane normal to the rotation axis in a cylindrical coordinate system and hat refers to the unit vector. The last term on the right hand side is the centrifugal force while the rest of the terms are the same as before.
In an isothermal halo, this stability equation can be met only if the halo takes a non-spherical shape of an ellipsoid. However, the shape of a non-isothermal halo does not necessarily depend on its rotation. The reason is that (i) the tem- perature of a halo is not necessarily a function of its shape and can be changed by halo's environment, and (ii) temperature gradient adds an extra component to the force due to the pressure.
As an example, in this section, we present the temperature profiles of four perfectly spherical but rotating halos with rotation periods of 0.8, 1, 1.2, and 1.7 ×10 18 year . For the sake of comparison, we choose their mass density, DM mass, and the rest of the parameters to be the same as the halo presented in Fig. 1. The temperature profiles that enforce the stability are shown in Fig. 2. The very bottom overlying solid lines are the profiles along the axis of rotation x 3 and as expected are the same as the temperature profile of the nonrotating halo. The non-solid lines are temperature profiles along the radial coordinate x 1 on the x 3 = 0 plane. Although the mass density is still spherical, the temperature profile has deviated from spherical symmetry. The equatorial edges of these rotating halos are heated up to confront the centrifugal force while the central regions are left intact.
In summary, halo rotation breaks the spherical symmetry of the stability equation, which can be accommodated by breaking the spherical symmetry of mass density, or temperature profile, or both. The complication is that the temperature is also a function of the environment and therefore is not fixed by mass density. This leads to a dramatically high number of possible stable solutions that can be ruled out only after the temperature profile is modeled. For mainly this reason, we only consider non-rotating spherical halos in the rest of this paper and leave the more general stability solutions for later.

Temperature profile of a fermionic DM halo at full degeneracy level
At the full degeneracy level, z ∞, and using the Sommerfeld approximation, the Fermi-Dirac integrals read Therefore, the ratio of the Fermi-Dirac integrals can be written in terms of the mass density The EOS in Eq. (5) reads which unlike the EOS of the classical ideal gas, has no temperature dependence. Consequently, the temperature profile becomes irrelevant in the full degeneracy limit, and the stability Eq. (11) leads to the well-known Lane-Emden equation with numerically known solutions. So far, we have studied the solutions for when every location of DM halo is at extremely low or full degeneracy levels. Since DM halo is denser and hotter at the center and less dense and colder at the edge, the degeneracy level of DM can vary with the distance from the center.

Software
Study of the most general scenario calls for numerically solving the most general stability equation because (i) the dependence of pressure on the mass density and temperature, both of which functions of r , takes a complicated nature in the partial degeneracy level, and (ii) the known solutions of specific degeneracy levels still depend on the initial conditions that are not known when transiting from one solution to the other.
Inserting the pressure and mass density from Eq. (2) into Eq. (11), the most general hydrostatic equilibrium equation for non-interacting fermions reads where the dimensionless variables are defined as ξ ≡ r , and s ≡ z z 0 . The boundary conditions are s 0 = y 0 = 1 which are implied by the definitions of s and y and ds dξ | ξ =0 = dy dξ | ξ =0 = 0 which can be understood from Eq. (9) knowing that M(r ) approaches zero faster than r 2 when we move toward the center.
The numerical solutions of differential equation (26) are sought in terms of y and Ln(s), instead of s. The latter choice is because the fugacity can take computationally infinite values. We have written computer code in Python [26], to study the solutions to the most general stability Eq. (26) for fermionic non-interacting halos. The mass of DM, the central values for the number density and the temperature, and the temperature profile in terms of ξ are taken from the user as inputs. The latter is to preserve the generality of the software and allows linking it to data-driven optimization methods for the phenomenological investigation of halo temperatures. The reason for this approach is that the equation for temperature gradient not only depends on the DM model but also is a function of the environmental variables of a specific galaxy.
The software starts from the center of the galaxy and uses the input values to compute the rest of the parameters like α and z 0 . Next, the numerical step size is determined using the Richardson extrapolation at the center. Depending on the magnitude of the error, the step size ( r ) can be between one-hundredth of a parsec (pc) to one parsec.
The software moves toward the edge until the density reaches one-thousandth of its value at the center. To reduce the numerical errors to the order of ( r ) 4 while preserving the computation speed, at each iteration step we use the Verlet method to find Ln(s) and its first derivative using d 2 (Ln(s)) dξ 2 from Eq. (26).
At any iteration, the software determines if the system is in the high or partial or non-degenerate regime. It then uses the appropriate approximations to calculate the Fermi-Dirac integrals. In the case of partial degeneracy where no asymptotic behavior is known, the software uses the optimization methods to learn the fugacity from the Fermi-Dirac integrals at that point. The code will be in its slowest mode when encountering the partial degeneracy level because, unlike the extreme degeneracy levels, given the mass density and temperature, it is not trivial to eliminate f 5 2 (z) of the pressure in Eq. (2) in terms of f 3 2 (z) of the mass density in the same equation.
The returned result will be in the form of a set of five plots showing the mass density in units of the critical mass density ρ c 9 × 10 −27 kg m −3 , the temperature, natural logarithm of the fugacity, the mass of the galaxy, and the chemical potential in the free-falling frame μ ≡ kT Ln(z). The software also reports the dynamical time, the total gravitational potential energy, and the total kinetic energy of the halo. All of the non-isothermal solutions reported in this paper have total kinetic energy equal to approximately half of their gravitation potential energy and as a result, are in the Virial equilibrium. Since the full Fermi-Dirac EOS is exclusively used in the software, the transition between degeneracy levels is smooth and the Pauli exclusion principle is always sustained.  Every reported solution is re-derived using step sizes of onetenth of the original step size to assure their stability against the global accumulation of the numerical errors.

Software validation
Since our code can be used to study any possible solution for non-interacting fermionic DM halo, and for validation purposes, in this section, we re-derive a set of known stability solutions.
We start with an isothermal model with DM mass of 1 MeV at the temperature of 1300 K, and a mass density of 8 × 10 −20 kg m −3 at 1 parsec from the center. The initial parameters are chosen such that the system is in the classical regime where the hydrostatic equation has an exact solution given in Eq. (16). Figure 3 shows the numeric and exact solutions for the given parameters. The mass density is the first plot from the left, is in fair agreement with the exact solution even though the latter does not satisfy the initial conditions due to its singularity at r = 0. The solution also contradicts observations due to its cuspy nature at the center. The second plot from the left is the temperature profile showing that it is assumed constant over the halo. The third plot is the logarithm of the fugacity indicating that z 1 at all times and confirms the classical nature of the solution. The fourth plot is the mass of the halo. It does not reach a flat plateau indicating another problem of this solution. The last plot is the chemical potential in the free falling frame. Since the temperature is constant, it is proportional to the logarithm of the fugacity in the third plot. However, the x-axis is in the logarithm scale providing the lost information at around the center.
To validate the software in the opposite of the spectrum, we reproduce a fully degenerate DM halo presented in [16] where a lower DM mass limit of 200 eV is derived using the observed line of sight dispersion velocity of dwarf galaxies. We set the DM mass to 200 eV at the temperature of 10 −4 K, and mass density of 10 −20 kg m −3 at 1 parsec from the center. The profile of the system is shown in Fig. 4 where the mass density and total mass of the system are in agreement with those reported in [16]. Since the temperature is not exactly zero, the logarithm of the fugacity is decreasing instead of being infinite as in [16]. Nevertheless, as far as this logarithm is large enough, the full degeneracy regime is approximately valid and the results are stable.
In the two solutions above, the entire DM halo was either non-degenerate or highly degenerate. We now reproduce the double plateau isothermal solution in [12,14] where DM halo is highly degenerate at the center, partially degenerate in the middle, and non-degenerate close to the edge. To achieve such solution, we choose the mass of DM to be 200 eV, the density at the center to be ρ 0 = 10 −21 kg m −3 , and the temperature at the center to be T 0 = 0.0003 K. The solution is depicted in Fig. 5. It should be noted that, unlike in the previous two solutions, both of the axes of the mass density are transformed into the logarithmic scales to reproduce the looks of the corresponding solutions in the references. Also, to capture the second plateau, we did not terminate the code until the density became 10 −6 , instead of 10 −3 , times the density at the center.
So far, we have reproduced the known isothermal solutions. Since non-isothermal solutions are not well investigated before, such comparison is not possible in that domain. However, we still have the analytic evaluations of Sects. 3.1 and 3.3 that, as we will discuss later, can validate a subset of non-isothermal solutions. For the rest of the solutions, their numerical accuracy is validated by first running the code in the normal mode and re-running it with intervals of one-tenth of the interval in the normal mode. We only report the solutions if the two trials lead to the same numerical solutions. We would like to mention that for extremely steep temperature gradients, the latter validation may fail, these results are not presented in this paper. It is always possible to reduce the intervals of the numerical method until a valid solution is reached at the cost of slower computations. However, we postpone the study of such domains until more advanced numerical methods are implemented.
Finally, the reported solutions are stable against small changes in the input values. An example of such study is shown in Fig. 6 where a solution with central temperature of T 0 = 0.1 K, temperature profile of y = 1 + 50ξ 2 −1 , and central mass density of ρ 0 = 10 −22 (kg m −3 ) is plotted together with solutions of seven small perturbations 3.6 Non-isothermal DM halos with y = 1 + bξ 2 −1 as the temperature profile A physically acceptable temperature profile does not take negative values and also does not increase monotonically with the distance from the center. It also must satisfy the initial conditions y 0 = 1 and dy dξ 0 = 0. We confirm that with b controlling the level of non-isothermality, satisfies all these requirements. First, we study non-isothermal effects by choosing DM mass of 100 eV-lowest possible DM mass if halos are isothermal [18], with central temperature of T 0 = 0.1 K, and central mass density of ρ 0 = 10 −22 (kg m −3 ) with different values of b. The results are shown in Fig. 7.
Before analyzing the outputs of the software, we would like to discuss their validity. The curve of b = 0 is the classic isothermal solution. We validated the software in this domain in the previous section in Fig. 3. From the third subplot, we can see that the solutions of b = 0.0001 and b = 0.001 are in the domain of low degeneracy with z 1 at any distance from the center. Therefore, the analytic formula in Eq. (17) should explain the temperature profile. We validate the latter by inserting the numerical mass density from the software into Eq. (17) and comparing the predicted temperature profile The temperature profile y = 1 + bξ 2 −1 with b = 0.001 as the true curve and the predicted profile using Eq. (17) and assuming that the mass density is given in Fig. 7 with the assumed one. This comparison is presented in Fig. 8 where the analytic temperature profile of software's mass profile is predicted to be exactly the same as the assumed one. The other curves with higher b values enter the partial degeneracy level at some distance from the center. Their accuracy is validated by re-running the code with reduced distance intervals as discussed in the previous section.
The very left subplot of Fig. 7 shows the predicted mass density. For small temperature gradients with b ≤ 0.1 (top panel), the mass density starts to increase at some large distance from the center and then drops to zero. Even the classical analytically validated solution for b = 0.001 possesses a clear overdensity at around 10 (kpc). For larger temperature gradients with b ≥ 10 (bottom panel), the mass profile has the shape of a doughnut with a "hole" at the center with a negligible radius. This dilute region is located where the visible matter and the central black hole are and may not be detectable. The observable central density will be the peak of the curve and is higher than what we have inserted in the software.
The total masses of the halos in Fig. 7 reach a flat plateau in non-isothermal solutions while continues to infinity in the isothermal case. Therefore, non-isothermal halos are more consistent with expectations.
The very right subplot in Fig. 7 shows that even the lowest degree of non-isothermality leads to an increasing chemical potential in terms of r as opposed to the decreasing chemical potential of the isothermal solution. Figure 7 also indicates that by increasing the temperature gradient, the radius, and the mass of halo decrease allowing for more compressed solutions. As we will see below, the two solutions with b ≥ 10 (bottom panel) are more compressed than their corresponding fully degenerate halos-as the most compact scenarios in isothermal cases.
We would like to mention that, unlike the solutions with constant temperature, in non-isothermal scenarios, the dis- Fig. 9 The escape velocity √ −2 and kT m corresponding to b = 1 of Fig. 7. DM particles at the center have higher kinetic energy than the gravitational potential energy. Near the edge, however, their kinetic energy falls toward zero due to the decrease in their temperature. The dynamical time of the system is 1.8 × 10 8 years. If an imaginary mass bubble moves from the center toward the edge, this extended period assures us that it will lose its kinetic energy to the colder surrounding and cannot escape the system by the time it reached the edge persion velocity of DM particles decreases with the distance from the center as can be seen in Fig. 9. Although the kinetic energy of DM particles is higher than the gravitational potential energy at the center, the kinetic energy at the edge is negligible and the dispersion velocity is much less than the escape velocity.
To explore the temperature profile in the high degeneracy level scenarios, we lower the central temperature to T 0 = 10 −5 (K) but keep the DM mass and central mass density to be m = 100 eV and ρ 0 = 10 −22 (kg m −3 ) respectively. From Sect. 3.3, we expect no dependence on the temperature profile and consequently on the value of b. The results for different b settings are shown in Fig. 10, and confirm the validity of our software in this domain.
In isothermal scenarios, the high degeneracy level solution presented in Fig. 10 is the most compact possible halo. However, by comparing Figs. 7 and 10, we can observe that more compressed halos are possible in the presence of temperature gradient. Figure 7 represents some halos that are more compressed than their corresponding fully degenerate halos of Fig. 10 because there is a pulling force of pressure in the nonisothermal cases that is absent in the isothermal solutions. The forces that are involved in Eq. (9) are the gravitational force on the right-hand side and the force due to the pressure on the left-hand side. Taking a derivative of the most general pressure in Eq. (2), the force of pressure reads The first force of the pressure is proportional to dy dr and is absent in isothermal solutions. If the temperature gradient has a negative sign, which we expect it to have, this is a pushing force that decompresses the halo. If as in Fig. 1, the temperature rise with the distance, the force will be inward.
The second force of the pressure is proportional to d (Ln(z)) dr which has a varying sign in non-isothermal solutions. It is a pulling force that compresses the halo at the center and a pushing force at the edge. In isothermal halos, it is only an outward force. In Fig. 11, the forces are depicted for isothermal b = 0 and non-isothermal b = 0.001, b = 1, and b = 10 solutions of Fig. 7. We would like to emphasize that the first two solutions are analytically validated. In non-isothermal halos, it is the inward force of the pressure that compresses the halo at around the center, and the gravitational force is negligible. From the figure, one can observe that the strengths of the forces of the pressure are orders of magnitude higher than the strength of gravity in the presence of temperature gradient. Also note that, toward the center, d P 2 dr has opposite signs for b = 0 and b = 0. Therefore, with a temperature gradient of just 0.1 K over thousands of light-years, one can neglect the force of gravity for the first few kilo-parsecs of the halo. With such strong additional pulling forces, it is easy to overcome the Pauli blocking forces even if DM particles are extremely light. In [27], a cosmological model for such light fermions is presented that does not contradict the observations of large scale structures. Later in Sect. 4, we discuss that such additional attractive forces can significantly lower the so-called phase-space lower bounds on the mass of fermionic DM that is currently derived in the literature for isothermal halos in the absence of the attractive forces associated with nonisothermality.
The effects of increasing the central temperature are shown in Fig. 12 where all of the solutions have DM mass of 100 eV and central mass density of ρ 0 = 10 −22 (kg m −3 ). Again, we use the generic profile of the form y = 1 + bξ 2 −1 where different T 0 values for both nonisothermal case of b = 100 (top panel) and isothermal case of b = 0 (bottom panel) are presented. As can be seen from the figure, by increasing the central temperature, the size of halo decreases in non-isothermal solutions and increases in isothermal solutions. Interestingly, the smallest halo in Fig. 12 (top) has the same central temperature as the largest halo in Fig. 12 (bottom).
These two opposite behaviors root back to the different natures of their fugacity-or equivalently the chemical potential, profiles. When b = 0, higher T 0 leads to a steeper tem-  g·r and ∇ P2·r g·r respectively. The force of gravity is always equal to the net of the two forces of the pressure. In the central region of the non-isothermal halos, it is the pulling force of the pressure rather than gravity that maintains the stability

Halo encounters
Rather than being isolated, DM halos are regularly encountering other galaxies. Such collisions are categorized depending on the sizes and velocities of the colliding galaxies. In general collisions lead to non-zero time derivative of the distribution function. The time evolution of the distribution function is governed by the Boltzmann equation where C[ f ] is a functional of the distribution function accounting for the collision, i.e., the gravitational forces due to the other system, is the self-gravitational potential, and v is the velocity of DM particles. When the encounter is long term, as in the case of major mergers, such that the time derivative of the distribution function is non-perturbatively different from zero, the system is not stable and stability solution as the subject of this paper becomes irrelevant.
When galactic collisions are fast, the distribution function deviates slightly from the initial form. Therefore, an approach similar to the one employed for studying the anisotropies in the cosmic microwave background should be viable. In this case, the distribution function of fermionic DM can be imagined to have a form similar to the Fermi-Dirac distribution where we assume that δz z and δβ β. After expanding around the initial Fermi-Dirac distribution and neglecting the second order terms, the distribution function reads 1 ze βε + 1 − e βε εzδβ + δz During the halo encounter, the collision term C[ f ] can be found from for example [28] using the cross-section of the gravitational interactions. After the encounter, the partial derivative of f with respect to time is still non-zero. The stability Eq. (9) is only valid at the very end of the process when ∂ f ∂t = 0. However, since both sides of the Boltzmann equation can be expanded in terms of δβ and δz, and since the zeroth and first order equations are independent, the stability solutions presented in this paper are always valid to the zeroth order.
In summary, the stability solutions of this paper are only applicable to halos that have not been involved in a major merger in the past few Giga-years. Recent studies show that the probability of such encounters is higher for more massive halos [29][30][31].
One interesting halo encounter is when the transferred energy is entirely spent to overcome the gravitation bound of the outer regions and leads to mass loss, such that C[ f ] and ∂ f ∂t are only non-zero in the outer regions of the halo. Therefore, the stability solutions presented in this paper are still valid in the interior regions because the gravitational force at any distance is only a function of the enclosed mass and the balance of the forces in the inner regions will not be changed significantly.
Tidally truncated subhalos submerged in a host are in general subject to at least dynamical friction which can significantly change their temperature profile and hence the density profile. However, due to Gauss's law, the gravitational force of the interior region is not affected by the mass loss and the stability solutions presented in this paper are well justified up until the limiting radius provided that the temperature profile reflects the encounter.

Phase-space mass bounds on DM mass
In 1979, Tremaine and Gunn derived the first lower limit on the mass of DM [32]. The derivation depends on a set of assumptions whose validities are not known yet. More specifically they assumed (I) a specific primordial phasespace density, (II) DM is collision-less, i.e., the maximum of its phase-space density is conserved, (III) galactic DM has a Maxwell-Boltzmann distribution, and (IV) DM halo is isothermal. The first two assumptions are not valid for interacting DM. The last two assumptions are also not valid in the degenerate non-isothermal models in which we are interested. The Tremaine-Gunn bound is more related to a knowledge of the primordial phase-space and its evolution over time (which are model dependent) than the Fermi-Dirac statistics of particles. The same bounds apply to some nonfermionic models of dark matter [33].
True models of galactic fermionic DM were later studied in for instance [12][13][14]. A lower bound on the mass of a genuinely fermionic DM can be derived using the lower limit of its dispersion velocity at the full degeneracy level, see for instance [18,34] If the Maxwell-Boltzmann distribution governs a coarsegrained fermionic DM, its dispersion velocity still needs to be larger than the minimum in Eq. (32), and the inequality leads to a lower limit on the mass of DM. In [16], it is discussed that the inequality is trivial if DM halo is made of a degenerate Fermi gas. As is mentioned by [18], the inequality is trivial as far as the Fermi-Dirac distribution is used to describe DM halo regardless of its degeneracy level, i.e. even if z 1. To see this, we start with the definition of the dispersion velocity 10 f 3 2 (z) 5 3 .
Even in the classical limit where the Fermi-Dirac distribution is effectively Maxwell-Boltzmann, the inequality is trivial. This can be seen by replacing the Fermi-Dirac integrals with their low-fugacity approximation. For this reason, our software exclusively (even in effectively Maxwellian DM halos) works with a Fermi-Dirac distribution such that the limitation of the phase-space is always respected. Also, as is mentioned in [18], if DM distribution is not entirely Maxwellian, it is not possible to use the dispersion velocity of entirely Maxwellian visible matter to learn the escape velocity of DM. Such learning becomes even less possible in non-isothermal halos because the dispersion velocity is a function of temperature profile and the mechanisms of heating the visible matter are due to the release of the stored potential energy of electric and strong forces, as well as gravity-based mechanisms. On the contrary, DM halos are most likely heated up by gravity-based mechanisms and friction. As can be seen from Fig. 9, in the outer regions of non-isothermal halos, the kinetic energy of DM particles are low enough that they cannot escape the gravitational well.
However, in [18], the lower bound of ∼ 100 eV for the mass of DM is derived if the observed dwarf galaxies are infinitely degenerate. The basis for the bound is that the decay period of satellite galaxies due to the Chandrasekhar friction has to be larger than 10 10 years. Chandrasekhar's estimation for the decay time is a function of (i) the velocity of the satellite, (ii) its distance from the center of the host galaxy, and (iii) mass and radius of the halo of the satellite.
Among the three enumerated factors above, only the mass and radius of a satellite's halo depend on the stability Eq. (26). If two solutions to the stability equation have the same halo mass and radius, their corresponding decay time will be the same. The lower bound of ∼ 100 eV in the reference above is derived for entirely degenerate halos-as the most compact possible solutions of isothermal scenarios. However, there is no observation confirming that such halos are infinitely degenerate. Due to the stronger frictional forces that such galaxies experience, it is likely that the halos are non-isothermal in which case the same halo size and mass are possible with lower DM mass. For example, in Fig. 13, we present a non-isothermal halo made of DM mass of 75 eV  Fig. 13-labeled true, are reproduced-labeled numeric, by numerically solving Eq. (40) and using Eq. (35) after assuming that the mass density in Fig. 13 with T 0 = 0.1 (K) is the true function tions, we can take this as another validation on top of the interval reduction validation method discussed in Sect. 3.5.
It is important to note that the escape velocity and the phase-space limitation arguments mentioned above are trivially satisfied by the solution to Eq. (40). Because, this is a stability equation whose solution guarantees that DM particles are trapped in the halo. Also, since the full Fermi-Dirac EOS is used for the derivation, the limitation of phase-space is naturally met.
Finally, we would like to emphasize that visible matter is often denser than DM at distances close to the center of halos. The extra gravitational force helps to contract DM even further. Since the gravitational force due to visible matter can be significant, fermionic DM halos are compacted much more than if the visible matter is neglected. Therefore, even in isothermal halos, the lower bound on the mass of DM is not as low as is currently derived.

Conclusion
We have studied non-isothermal non-interacting fermionic spherical dark matter halos. Using the full EOS of Fermi-Dirac statistical systems, we derive the most general stability equation and present computer software to numerically solve it. Since the full Fermi-Dirac EOS is used in the software, the transitions between degeneracy levels are smooth and the limitation of fermionic phase-space is never violated in the numerical solutions. From non-degenerate to highly degenerate Fermi halos with any temperature profile can be investigated with the software.
We have studied non-isothermal halos using a generic temperature profile of the form T = T 0 1 + r r 0 2 −1 , and shown that their chemical potential profile is substantially different from that of the corresponding isothermal halos. We show that the mass and radius of such non-isothermal halos decrease by increasing the temperature gradient. We have shown that the force due to the pressure has inward as well as outward components, and at the central regions of the studied halos, it is the inward force of the pressure, rather than gravity, that maintains the stability. We have discussed the phase-space lower bounds on the mass of DM as well as the importance of modeling the temperature profile of DM halos for deriving them. We have shown that if the temperature is left as a free parameter, any arbitrarily light DM mass can explain the observed mass profile of DM halos.
It has been discussed that the limitation of the phase-space of fermions does not restrict their configuration volume. It is the inward force acting on the fermions that determines the size of DM halos. In the presence of temperature gradient, the inward force due to the pressure adds to the inward force of gravity and maintains the stability of compressed fermionic halos. We have shown examples where the former force is orders of magnitude stronger than the latter for a rather vast region. The inward component of the force of pressure is absent in isothermal halos.
We have shown that if the quantum nature of DM is irrelevant in the halos, the temperature profile is analytically given in terms of the mass profile. By requiring that the temperature is not negative, we place an upper bound on the mass of DM. We find that if the central temperature of DM halo is only a few Kelvins, the mass of DM cannot be larger than a few keV.