Condensation of Galactic Cold Dark Matter

We consider the steady-state regime describing the density profile of a dark matter halo, if dark matter is treated as a Bose-Einstein condensate. We first solve the fluid equation for"canonical"cold dark matter, obtaining a class of density profiles which includes the Navarro-Frenk-White profile, and which diverge at the halo core. We then solve numerically the equation obtained when an additional"quantum pressure"term is included in the computation of the density profile. The solution to this latter case is finite at the halo core, possibly avoiding the"cuspy halo problem"present in some cold dark matter theories. Within the model proposed, we predict the mass of the cold dark matter particle to be of the order of M_chi c2 = 10^-24 eV, which is of the same order of magnitude as that predicted in ultra-light scalar cold dark matter models. Finally, we derive the differential equation describing perturbations in the density and the pressure of the dark matter fluid.


Introduction
Observations on the rotational curves of spiral galaxies show that the velocities of the virialized material lying farther than the extent of the luminous matter from the galactic center reach a constant value [1]. Various theories aim at explaining this discrepancy between observations and Newton's virial theorem, including a modification of the gravitational potential [2,3] or of the Poisson equation [4,5], conformal gravity [6,7], and the metric skew tensor gravity [8,9,10].
Nowadays, the most promising way to explain the observations of the galactic rotation curves [11,12] consists in postulating the existence of non-luminous (dark) matter, distributed in a halo which extends much farther than the luminous component of a galaxy. Further, this dark component is supposed to be non-relativistic (Cold Dark Matter, CDM [13,14,15], see also [16,17]), since it is usually assumed to consist of massive particles with very low thermal velocities. Work on colliding galaxy clusters seem to confirm the existence of dark matter dominating the mass content of spiral galaxies and galaxy clusters [18,19]. An indirect confirmation also comes from the success of the concordance cosmological model, or Λ Cold Dark Matter (ΛCDM) model, in reproducing the anisotropies observed in the cosmic microwave background [20,21]. Among the most promising candidates for the CDM component are the Weakly Interacting Massive Particle (WIMP, [22]) or a population of zero-momentum axions [23,24,25,26,27]. However, these "canonical" CDM models usually feature problems in reproducing some observable properties of galaxies, most remarkably the under-abundance of small scale structure with respect to what predicted from simulations (the "missing satellite" problem [28,29,30]), the absence of a central density cusp rather predicted in numerical simulations (the "cusp" problem [31,32,33,34]), and too little dense subhalos compared with theoretical expectations (the "too big to fail" problem [35]). In more detail, observations of both nearby dwarf galaxies and low surface brightness galaxies show that the density profile of the CDM halo at the core reaches a constant value [36,37,38,39,40,41,42]. In contrast, various N-body simulations predict that the CDM density distribution steepens at the center of the halo [43,44,45,46,47].
Among the solutions proposed to overcome these issues, it has been suggested that dark matter could consist of a coherent scalar field with long range correlation, whose quanta are very light particles [48,49,50,51,52]. In fact, on short length scales, light scalar fields do not behave as perfect CDM and would inhibit cosmological structure growth [53,54,55]. This solution would consist of a viable mean to suppress low mass galaxies and provide cored profiles in CDM-dominated galaxies [51,52,56]. This peculiar form of dark matter might form a Bose-Einstein Condensate (BEC), described by the Gross-Pitaevskii or non-linear Schrödinger equation [57,58,59,60,61,62,63,64,65,66,67,68,69,70,71]. Alternatively, axions can also be modeled as a coherent BEC with small spatial gradient [72,73,74,75,76,77]. Recent 3D simulations of the gravitational collapse of wavelike cold dark matter (ΨDM) with a mass of the order of 10 −22 eV show the effects on structure formation due to their large Jeans length [78,79,80], which improved over previous work on the subject [57,58,59,60]. Other alternative models embed dark matter condensation into space-time with torsion [81]. A general review of the models proposed is given in Ref. [82]. The cosmological evolution of a BEC dark matter component has also been extensively explored [83,84,85,86].
At galactic scales, the evolution of a self-gravitating CDM system can be described as a fluid following the equation of continuity and the Navier-Stokes Equation (NSE). When these equations are derived for a Bose-Einstein fluid, an additional "Quantum Pressure" (QP) term appears [59].
In this paper, we treat cold dark matter as a pressure fluid with a dynamics described by the NSE, and we derive the equations for the zeroth-and first-order perturbations in the density and pressure of "canonical" and BEC cold dark matter. For this, we assume a rotating halo in which the proper velocity of dark matter is treated as a first order perturbation in the motion.
The paper is organized as follows. After the short review of fluid dynamics in Sec. 2, we discuss some popular halo models fitting numerical simulations in Sec. 3. In Sec. 4 we show that, at the lowest order, the dark matter density in the halo follows the Lane-Emden equation in a rotating frame. When quantum pressure is included, the Lane-Emden equation modifies as discussed in Sec. 5, and results for the halo profile predict a finite core. In Appendix 6, a generic expression for density perturbations and the proper velocity of dark matter is derived and will be part of future work.

Equation for fluid dynamics
Newton's equations for a parcel of density ρ and proper velocity v, written in a reference frame with theẑ axis in the direction of increasing altitude, reads Here, p is the pressure acting on the parcel, φ is the gravitational potential, and τ describes all additional external forces in the system, like the mean gravitational field generated by all nearby galaxies. In addition, Π is a rank-two tensor describing the dissipative phenomena in the fluid, The two constants appearing in Eq. (2.2) are known in the literature respectively as the dynamic viscosity η and the second viscosity coefficient ζ [87,88]. The total time derivative of the velocity field can be explicitly written as the sum of a partial time derivative and the dyadic (advection) term which introduces a non-linear component in Newton's equation, where v = |v| and we have defined the vorticity of the velocity field as We assume that the galactic halo rotates at a constant rate Ω, and we switch to the rotating frame by setting v → v + Ω × r, obtaining the NSE in the rotating frame Here, Ω×Ω×r and 2 Ω×v are respectively the Coriolis and the centrifugal acceleration terms.
The NSE couples to two additional equations which express flux conservation (continuity equation), dρ dt + ρ (∇ · v) = 0, (2.6) and the value of the gravitational potential generated by the matter density ρ (Poisson equation), In the following, we look for a solution to the set of Eqs. (2.5) and (2.6) in the steady-state regime, ∂v ∂t = ∂ρ ∂t = 0. (2.8) Furthermore, since we are treating the velocity as a first order term in the perturbation series, we neglect all advection terms. Under these conditions, Eq. (2.5) reads while Eq. (2.6) in the steady-state regime is rewritten as the incompressibility relation ∇·v = 0 for the DM flow.

Density profiles from N-body simulations
In the following, we define the profiles so that to match the local halo density ρ ⊙ = 0.4 GeV/cm 3 [89,90,91,92,93] (see also Ref. [94]) when the radius r equates the distance of the Sun from the center of the halo r ⊙ ≈ 8.5 kpc. The ΛCDM model is very successful in predicting various observational features such as the galaxy rotation curves [11,12], which are closely related to the nature of the CDM and the formation history of the halo. Observations show that the shape of rotation curves of CDM-dominated galaxies is universal across a wide mass range [37,38,39].
Early work on the analytical shape of the density profile [95,96,97,98,99,100], showed that an isothermal profile with ρ ∼ r −2 is to be expected in a flat Friedmann-Robertson-Walker Universe. In addition, the density profile would steepen when sharper spectra of the initial density perturbation are assumed. These claims were first validated in simulations [101,102,103,104,105]. The isothermal profile [100] is derived from balancing the gravitational and pressure forces, assuming that pressure is linearly dependent on the matter density, where q = r/r ⊙ , q s = r s /r ⊙ , and ρ iso and r s respectively parametrize density and radial size of the Galaxy. The isothermal profile has the characteristic to remain finite when r → 0, approaching the value However, the isothermal profile is too shallow when compared with results from numerical simulations of non-colliding dark matter. Departures from the power-law behavior were first reported in Refs. [106,107,108]. A recent review has been given in Ref. [109]. Using numerical simulations, Navarro, Frenk, and White [NFW, 31,43,44] fitted results for CDM halos with mass spanning over four orders of magnitude with the function thus predicting a rather dense central cusp growing with ∼ 1/q. Additional numerical simulations [46,47] motivated Moore [30,32] to propose the following form for the CDM profile which diverges at the core as ∼ 1/q 3/2 . The isothermal, NFW, and Moore profiles, together with other parametrizations such as the BE [110] and PISO [111] profiles, can collectively be described by the generalized expression [112] ρ CDM (q, q s α, β, γ) = ρ ⊙ Θ CDM (q, q s , α, β, γ), The function Θ CDM (q, q s , α, β, γ) is normalized so that In Table 3, we have summarized the values of the parameters α, β, and γ for these popular halo models, including references. More recent N-body simulations [113,114,115,116,117,118,119] favor three-parameter profile models like the Einasto profile [120,121], a generalization over the Sérsic model [122], where r s and δ are constants. The halo density at the core predicted by the Einasto profile is finite, with In Figure 1, we compare the plots of the isothermal (dot-dashed line), NFW (dotted line), Moore (solid line), and Einasto profiles (dashed line), with the normalization satisfying Θ CDM (r ⊙ ) = 1. For the profiles considered, we have set r s equal to the value given in Table 3, while the parameters for the Einasto profile are r s = 20 kpc and δ = 0.17 1 .

Equations for the halo profile neglecting quantum pressure corrections
We assume that the DM stream velocity is v ≪ Ω L, where L is the typical galactic length scale. To give a numerical example, for a period of rotation Ω −1 = 200 My and for a length scale L = 50 kpc, we obtain v ≪ 250km/s. For this reason, we first look at a numerical resolution of the set of Eqs. (2.5)-(2.7) in which we neglect the stream velocity. Under these conditions, the set of equations describing the balance between pressure and density in a galactic CDM halo is Eq. (4.2) expresses mass conservation around an infinitesimal volume. The curl of Eq. (4.1) yields ∇ρ × ∇ p = 0, which is a structural condition between the bulk pressure and density which is fulfilled by the barotropic relation which is known in the literature as the Lane-Emden equation [124,125,126]. For the generic barotropic relation expressed in Eq. (4.4), the Lane-Emden equation requires a numerical resolution [127,128,129]. An analytic solution exists when the relation is polytropic p ∝ ρ 1+1/n , with n = 1, 2, 5. When Ω = 0, Eq. (4.5) has often found applications in the study of collision-less systems such as globular clusters and primordial galaxies [1]. The rotating Lane-Emden equation in cylindrical coordinates has been discussed by Stodolkiewicz [130] and Ostriker [131] for the case of a non-rotating isothermal cylinder, by Schneider and Schmitz [132] for a generic polytropic fluid, and by Christodoulou and Kazanas [133] in the context of planetary formation for a linear polytropic relation p ∝ ρ (see also Refs. [134,135]).
Here, we consider the case in which the rotation is not neglected, assuming that both density and pressure of DM in the galactic disk do not depend on the azimuthal coordinate φ. We assume a cylindrical symmetry of the density, by setting ρ = ρ ⊙ Θ(r), (4.6) where Θ = Θ(r) is a function depending on the distance from the galactic center r only, normalized so that Θ(r ⊙ ) = 1. In cylindrical coordinates, Eq. (4.5) is expressed as To enforce the barotropic relation in Eq. (4.4), we assume that the dark matter BEC behaves as a self-interacting polytropic fluid, with pressure depending on density as [92,100,136,137] p = U 2 ρ, (4.8) where U is a constant with dimensions of a velocity, known as the isothermal sound speed [133].
Here, we parametrize the relation between the bulk density and pressure as Defining the galactic rotational period T = 2π/Ω, we introduce the parameter ω = √ 2 Ω τ = 374 My T , (4.10) whose order of magnitude of the parameter is approximately one. Defining and switching to the new variable q = r/r ⊙ , the Lane-Emden Eq. (4.7) is rewritten as In the following, we solve Eq. (4.12) with the boundary conditions that Θ(q) and its derivative match the corresponding quantities from the Einasto profile at the solar neighborhood: 4.1 Non-rotating halo ω = 0 When ω = 0, the solution to the differential Eq. (4.12) reads where q s is related to q U by (4.15) The expression in Eq. (4.14) ensures that Θ(q) > 0 for any value of q, and satisfies the normalization Θ(1) = 1. The halo profile in Eq. (4.14) can be rephrased in terms of the function defined in Eq. (3.6) as Θ(q, κ) = Θ CDM (q, q s , κ, 2 + κ, 2 − κ). Notice that, for κ = 1, Eq. (4.14) reduces to the NFW profile given in Eq. (3.3). We have thus obtained the important result that the NFW profile can be obtained as a solution to the non-rotating Lane-Emden equation. We search for the value of the parameter κ which is consistent with the boundary condition for the derivative of the halo profile, as given in Eq. (4.13). This yields to a transcendental expression for κ, which has to be solved numerically. Setting for example q s = (20kpc)/r ⊙ , we obtain κ = 0.81, which is the value we adopt later in Sec. 4.2. In Fig. 2, we plot the value of κ obtained from Eq. (4.17) as a function of q s . Fig. 3 shows the halo density profile, in units of ρ ⊙ , as a function of the distance form the galactic center q = r/r ⊙ . We set κ = 0.2 (blue solid line), κ = 0.5 (red line), κ = 1 (green line), and κ = 2 (dark green line). For any value of κ, all functions drop to zero for large values of q. At the halo core, the functions tends to infinity for 0 < κ < 2, with the steepness of the halo density near the galactic center decreases for increasing values of κ. In the limit κ = 2, the halo density profile at the halo core reaches a finite value at the galactic center, similarly to the prediction of the isothermal halo example in Eq. (3.1) where ρ iso (0) = ρ ⊙ . Despite avoiding the cuspy halo problem, the result for κ = 2 cannot be reconciled neither with observations nor with galactic simulation. The galactic rotational curves follow from the density profile in Eq. (4.14). In fact, inserting the density profile into the Poisson Eq. (2.7) written in the form and integrating twice from q ′ = 0 to q ′ = q, we obtain the expression for the gravitational field where Notice that, for large values of q ≫ q s , the potential in Eq. (4.19) has the analytical expression given in Ref. [138] φ ≈ v 2 rot ln q s q , for q ≫ q s . where ∇φ is the gradient of the potential expressed in cylindrical coordinates. We clearly obtain v ≈ v rot for q ≫ q s . In Fig. 4 we plot the virial velocity from Eq. (4.22) as a function of q for κ = 0.2 (blue solid line), κ = 0.5 (red line), κ = 1 (green line), and κ = 2 (dark green line). We obtain that the rotational curves rise from the value v = 0 at the origin, to reach a flat shape for r > r ⊙ .

Rotating halo ω = 0
For ω = 0, we solve Eq. (4.12) numerically with the boundary conditions given in Eq. (4.13).  Fig. 5 shows an unrealistic behavior for large radii r ≫ r ⊙ , where the halo density profile exponentially drops to the constant value ω 2 , instead of dropping to zero as for the case of the Einasto profile, the NFW profile, and the non-rotating solution. This behavior is due to the fact that, for r ≫ r ⊙ , the derivative term in Eq. (4.12) drops to zero as 1/r 2 , and the equation reduces to Θ = ω 2 . We have further checked that, in this limit, expanding the exact solution Θ into a perturbation series, the first order perturbation also explains the exponentially damped oscillations observed in the full solution in Fig. 5. To sum up, since Θ = ω 2 for large values of r, we have obtained an unrealistic result because the angular velocity Ω in Eq. (4.12) is assumed constant. To  circumvent this problem, we modify the Lane-Emden equation by assuming that the angular velocity depends on q as where Ω 0 is the value of Ω when q = 0, and f (q) is a generic function of q satisfying 2 f (0) = 1 and lim q→+∞ f (q) = 0. (4.24) The Lane-Emden Eq. (4.12) is thus modified as The results shown on the right panel of Fig. 5 are then obtained for the constant function f (q) = 1 . Here, we use the sample function which satisfies the conditions in Eq. (4.24), to obtain a solution to Eq. (4.25) that fits the numerical results at large radii, as shown in the right panel of Fig. 5. At the halo core, the numerical resolution of Eq. (4.12) grows indefinitely for q → 0, for ω < 1 and f (q) = 1, and for any value of ω and f (q) as in Eq. (4.26). This behavior is similarly to what obtained for the NFW, Moore, and the halo profile with ω = 0 shown in Fig. 3. Within the framework of Eq. (4.7), the divergence of ρ at the halo core has been addressed by Christodoulou and Kazanas [133], who suggest to use a composite model in which the hydrostatic solution only applies at large radii, while the halo distribution is truncated to a constant value at small radii. Here, we suggest a different solution which involves the addition of a quantum pressure term, as discussed in Section 5.

Equations for the halo profile including quantum pressure corrections
In the literature, the set of Eqs. (2.5)-(2.7) describes classical fluid dynamics. An analogous set of equations might be derived in the context of BEC cold dark matter starting from the Gross-Pitaevskii equation [59,60,63,64,65,66,67]. The proof considers the Schroedinger equation for the wave function ψ(r, t) describing a CDM particle, where M χ is the mass of the dark matter particle, φ(r) the gravitational potential satisfying the Poisson Eq. (2.7), the potential giving the Coriolis and centrifugal forces is the potential yielding the viscosities is φ η = −η r · ∇v, and F (ρ) a generic function of the number density, ρ = |ψ(r, t)| 2 . From this latter expression, it follows that the wave function can be described as ψ(r, t) = √ ρ e i S(r,t) , where S(r, t) is the action of the particle. Once the expression for ψ(r, t) has been substituted into the steady-state regime of the Schroedinger Eq. (5.1), the real component of the expression obtained gives the Navier-Stokes Eq. (2.9) with an extra "quantum pressure" (QP) term appearing on its right-hand side [60,61,63]. Since ∇S is the momentum of the particle, setting v = ∇S/M χ gives where the pressure is given as a function of F (ρ) as The dependence of the QP term on the inverse of the CDM particle mass squared is then predicted by the Schroedinger equation [59]. The additional QP term does not modify the barotropic relation p = p(ρ), since it does not appear in the curl of Eq. (5.2). At the same time, the divergence of Eq. (5.2) (the Lane-Emden equation) contains an additional term with respect to Eq. (4.5), In cylindrical coordinates, assuming that density depends on the radial coordinate only, and switching to the adimensional variable q = r/r ⊙ , Eq. (5.4) rewrites as where β is an adimensional quantity given by (5.6) Solving Eq. (5.6) for the mass of the CDM particle and reinserting natural units, we find a result within the mass range given in Refs. [53,54,55,78,79,80]. In facts, a CDM particle with a mass in the range 10 −22 ÷10 24 eV is required to solve the cusp problem and to suppress the small-scale spectrum [51]. Ref. [62] independently obtain a mass range similar to the one given in this work, in the context of galactic CDM BEC. We have shown the value of the mass of the CDM particle as a function of β in Fig. 6. Eq. (5.5) is solved numerically with the boundary conditions in Eq. (4.13), plus the additional requirement that the second and third derivatives of the numerical resolution match the corresponding quantities in the Einasto profile at the solar neighborhood as well. For the numerical resolution, we have used the radius dependence of the angular momentum given in Eq. [65] that the BEC cold dark matter assumption alleviates the cuspy core problem appearing when simulating the evolution of dark matter cores. The stability of such halo, which in principle is not guaranteed [69], will be the subject of a future study. We fit the results for the value of the halo profile at the halo core from the numerical resolution to Eq. 5.5, obtaining the result in Eq. (5.8). The fit shows a kink atq = 0.0012, and reaches a constant value for β → +∞, corresponding to a decreasing mass of the DM particle. In fact, for β 1, the mass density at the halo core reaches the value We show the fit in Fig. (8).
In Appendix 6, we derive the expressions for density and pressure perturbations, assuming that they are of the same order of magnitude as the free-streaming velocity |v|. The solution and discussion to this set of equations will be the subject for further study.

Summary
Ultra-light scalar dark matter with long-range correlations has been considered as a valid alternative to the ordinary WIMP paradigm. In this paper, we have discussed the halo density profile obtained when treating CDM as a perfect fluid. Our work corroborates the  suggestion led by some authors [65] that the "cuspy" halo core, predicted in "canonical" CDM models, is removed when long-range correlations are included. For this, we have first solved the fluid equation for CDM, obtaining that the density profile diverges at the core, see Fig. 3. When an extra "quantum pressure" term arising from these long-range correlations is included into the fluid equation, the numerically-computed density profile remains finite at the halo core, see Fig. 7.
The curl of Eq. (6.5) results in the expression 2∇ × Ω × v = η ∇ 2 ξ, (6.12) where we introduced the vorticity ξ = ∇ × v. We parametrize the velocity in terms of three new adimensional functions u, v, w, depending on r only, as v = r s τ ur + vφ + wẑ . (6.13) Combining the three components of Eq. (6.12) with the continuity equation gives du dr + u r = 0, (6.14) A common solution to the set of Eqs. (6.14)-(6.16) that avoids a divergence at infinity is v = v 0 /r, for a constant vector v 0 . We now derive the expression for the divergence of Eq. (6.5). Using the barotropic relation in Eq. (4.8), the relation between pressure and density perturbations is 17) or, writing the series expansion of the function Θ = Θ 0 + Θ 1 , where Θ 0 is the solution to Eq. (5.5) and Θ 1 a small perturbation, we obtain ρ 1 = ρ ⊙ Θ 1 , and p 1 = U 2 ρ ⊙ Θ 1 . (6.18) Using the continuity Eq. (6.6), the divergence of the dissipation term is ∇ 2 (∇ · v) = 0. Once the Poisson Eq. (6.7) is taken into account, the divergence of Eq. (6.5) is a differential equation for Θ 1 , where the expressions for the coefficients O i are