Dark Matter as a Non-Relativistic Bose-Einstein Condensate with Massive Gravitons

We confront a non-relativistic Bose--Einstein Condensate (BEC) model of light bosons interacting gravitationally either through a Newtonian or a Yukawa potential with the observed rotational curves of $12$ dwarf galaxies. The baryonic component is modelled as an axisymmetric exponential disk and its characteristics are derived from the surface luminosity profile of the galaxies. The purely baryonic fit is unsatisfactory, hence a dark matter component is clearly needed. The rotational curves of five galaxies could be explained with high confidence level by the BEC model. For these galaxies, we derive: (i) upper limits for the allowed graviton mass; and (ii) constraints on a velocity-type and a density-type quantity characterizing the BEC, both being expressed in terms of the BEC particle mass, scattering length and chemical potential. The upper limit for the graviton mass is of the order of $10^{-26}$ $\text{eV/c}^2$, three orders of magnitude stronger than the limit derived from recent gravitational wave detections.


Introduction
The universe is homogeneous and isotropic at scales greater than about 300 Mpc. It is also spatially flat and expanding at an accelerating rate, following the laws of general relativity. The spatial flatness and accelerated expansion are most easily explained by assuming that the universe is almost entirely filled with just three constituents, namely visible matter, Dark Matter (DM) and Dark Energy (DE), with densities ρ vis , ρ DM and ρ DE , respectively, such that ρ vis + ρ DM + ρ DE = ρ crit ≡ 3H 2 0 /8πG ≈ 10 −26 kg/m 3 (where H 0 is the current value of the Hubble parameter and G the Newton's constant), the so-called critical density, and ρ vis /ρ crit = 0.05, ρ DM /ρ crit = 0.25 and ρ DM /ρ crit = 0.70 [1,2]. It is the large amount of DE which causes the accelerated expansion. In other words, 95% of its constituents is invisible. Furthermore, the true nature of DM and DE remains to be understood. There has been a number of promising candidates for DM, including weakly interacting massive particles (WIMPs), sterile neutrinos, solitons, massive compact (halo) objects, primordial black holes, gravitons, etc., but none of them have been detected by dedicated experiments and some of them fail to accurately reproduce the rotation curves near galaxy centers [3,4]. Similarly, there has been a number of promising DE candidates as well, the most popular being a small cosmological constant, but any computation of the vacuum energy of quantum fields as a source of this constant gives incredibly large (and incorrect) estimates; another popular candidate is a dynamical scalar field [5,6]. Two scalar fields are also able to model both DM and DE [7]. Extra-dimensional modifications through a variable brane tension and five-dimensional Weyl curvature could also simulate the effects of DM and DE [8]. In other theories, dark energy is the thermodynamic energy of the internal motions of a polytropic DM fluid [9,10]. Therefore, what exactly are DM and DE remain as two of the most important open questions in theoretical physics and cosmology.
Given that DM pervades all universe, has mass and energy, gravitates and is cold (as otherwise it would not clump near galaxy centers), it was examined recently whether a Bose-Einstein condensate (BEC) of gravitons, axions or a Higgs type scalar can account for the DM content of our universe [11,12]. While this proposal is not new, and in fact BEC and superfluids as DM have been considered by various authors , the novelty of the new proposal was twofold: (i) for the first time, it computed the quantum potential associated with the BEC; and (ii) it showed that this potential can in principle account for the DE content of our universe as well. It was also argued in the above papers that, if the BEC is accounting for DE gravitons, then their mass would be tightly restricted to about 10 −32 eV/c 2 . Any higher, and the corresponding Yukawa potential would be such that gravity would be shorter ranged than the current Hubble radius, about 10 26 m, thereby contradicting cosmological observations. Any lower and unitarity in a quantum field theory with gravitons would be lost [35].
In this paper, we discuss the possibility of a BEC formed by scalar particles, interacting gravitationally through either the Newton or Yukawa potential. Such a BEC, interacting only through massless gravitons has been previously tested as a viable DM candidate by confronting with galactic rotation curves [30,36].
In this paper, we solve the time-dependent Scrödinger equation for the macroscopic wavefunction of a spherically symmetric BEC, where in place of the potential we plug-in a sum of the external gravitational potential and local density of the condensate, proportional to the absolute square of the wavefunction itself, times the self-interaction strength. The resultant non-linear Schrödinger equation is known as the Gross-Pitaevskii equation. For the self-interaction, we assume a two-body δ-function type interaction (the Thomas-Fermi approximation), while we assume that the external potential being massive-gravitational in nature, satisfying the Poisson equation with a mass term. The BEC-forming bosons could be ultra-light, raising the question of why we use the non-relativistic Schrödinger equation. This is because, once in the condensate, they are in their ground states with little or no velocity, and hence non-relativistic for all practical purposes. Solving these coupled set of equations, we obtain the density function, the potential outside the condensate and also the velocity profiles of the rotational curves. We then compare these analytical results with observational curves for 12 dwarf galaxies and show that they agree with a high degree of confidence for five of them. For the remaining galaxies, no definitive conclusion can be drawn with a high confidence level. Nevertheless, our work provides the necessary groundwork and motivation to study the problem further to provide strong evidence for or against our model. This paper is organized as follows. In the next section, we set the stage by summarizing the coupled differential equations that govern the BEC wavefunction and gravitational potential and find the BEC density profiles. In Section 3, we construct the corresponding analytical rotation curves. In Section 4, we compare these and the rotational curves due to baryonic matter with the observational curves for galaxies. In Section 5, we find most probable bounds on the graviton mass, as well as derive limits for a velocity-type and a density-type quantity characterizing the BEC.

Self-Gravitating, Spherically Symmetric Bec Distribution in the Thomas-Fermi Approximation
A non-relativistic Bose-Einstein condensate in the mean-field approximation is characterized by the wave function ψ(r, t) obeying known as the Gross-Pitaevskii equation [37][38][39]. Here,h is the reduced Planck constant, r is the position vector; t is the time; ∆ is the Laplacian; m is the boson mass; is the probability density; the parameter λ > 0 measures the atomic interactions and is also related to the scattering length [40], characterizing the two-body interatomic potential energy: and finally V ext (r) is an external potential. For a stationary state, where µ is a chemical potential energy [40,41]. When µ is constant, Equation (1) reduces to present works [22], [30] where V Q is the quantum correction potential energy: We mention that Equation (5) is valid in the domain where ρ (r) = 0. The quantum correction V Q has significant contribution only close to the BEC boundary [21], therefore it can be neglected in comparison to the self-interaction term λρ. This Thomas-Fermi approximation becomes increasingly accurate with an increasing number of particles [42].
We assume V ext (r) to be the gravitational potential created by the condensate. In the case of massive gravitons, it is described by the Yukawa-potential in the non-relativistic limit: with ρ BEC = mρ, gravitational constant G, and characteristic range of the force R g carried by the gravitons with mass m g . The relation between R g and m g is R g =h/ m g c , where c is the speed of light andh is the reduced Planck constant. The Yukawa potential obeys the following equation: Contrary to Equation (5), Equation (8) is also valid in the domain where ρ (r) = 0. In the massless graviton limit, we recover Newtonian gravity, in particular Equations (7) and (8) reduce to the Newtonian potential and Poisson equation.

Mass Density and the Gravitational Potential inside the Condensate
The Laplacian of Equation (5) using Equation (8) gives For a spherical symmetric matter distribution, Equations (8) and (9) become where we introduced the notation 1 This system gives the following fourth order, homogeneous, linear differential equation for rρ BEC : with In the case of massless gravitons, πR * gives the radius of the BEC halo [30]. To have a real Λ, R g > R * should hold, constraining the graviton mass from above. Typical dark matter halos have πR * of the order of 1 kpc which gives the following upper bound for the graviton mass: m g c 2 < 4 × 10 −26 eV. Then, the general solution of Equation (13) is with integration constants A 1 , B 1 , C 1 and D 1 . This is why we impose the reality of Λ. For the imaginary case the general solution would contain runaway hyperbolic functions. This is also the solution of the system in Equations (10)-(11). Requiring ρ BEC to be bounded, we have D 1 = −B 1 . Then, the core density of the condensate is and the solution can be written as Substituting ρ BEC (r) in Equation (11), the gravitational potential is Being related to the mass density by Equation (5) gives The BEC mass distribution ends at some radial distance R BEC (above which we set ρ BEC to zero), allowing to express C 1 in terms of ρ (c) , R BEC and Λ as Finally, we consider the massless graviton limiting case m g → 0. Then, R g → ∞ implies Λ = √ 4πGm 2 /λ = 1/R * and C 1 = 0 (by Equation (19)). Then, ρ BEC (r) coincides with Equation (40) [22].

Gravitational Potential Outside the Condensate
The potential U is determined up to an arbitrary constant A 2 , i.e.
Here, U out Y satisfies Equation (8) Since an exponentially growing gravitational potential is non-physical, C 2 = 0 and The constants A 2 and B 2 are determined from the junction conditions: the potential is both continuous and continuously differentiable at r = R BEC : In the next section, we see that the continuous differentiability of the gravitational potential coincides with the continuity of the rotation curves.

Rotation Curves in Case of Massive Gravitons
Newton's equation of motions give the velocity squared of stars in circular orbit in the plane of the galaxy as Here, R is the radial coordinate in the galaxy's plane and U is the gravitational potential. In the case of massive gravitons, U is given by U = U Y + A, where U Y satisfies the Yukawa-equation with the relevant mass density and A is a constant.
The contribution of the condensate to the circular velocity is for r ≤ R BEC and for r ≥ R BEC .
In the relevant situations, the stars orbit inside the halo and their rotation curves are determined by the parameters: ρ (c) R 2 * , R BEC and Λ. In the limit m g → 0, the v 2 of the BEC with massless gravitons is recovered, given as Böhmer proposed [22] v 2 for r ≤ R BEC and for r ≥ R BEC .

Contribution of the Baryonic Matter in Newtonian and in Yukawa Gravitation
The baryonic rotational curves are derived from the distribution of the luminous matter, given by the surface brightness S = F/∆Ω (radiative flux F per solid angle ∆Ω measured in radian squared of the image) of the galaxy. The observed S depends on the redshift as 1/(1 + z) 4 , on the orientation of the galaxy rotational axis with respect to the line of sight of the observer, but independent from the curvature index of Friedmann universe. Since we investigate dwarf galaxies at small redshift (z < 0.002), the z-dependence of S is negligible. Instead of S given in units of solar luminosity L per square kiloparsec (L /kpc 2 ), the quantity µ given in units of mag/arcsec 2 can be employed, defined through where R is the distance measured the center of the galaxy in the galaxy plane and M is the absolute brightness of the Sun in units of mag. The absolute magnitude gives the luminosity of an object, on a logarithmic scale. It is defined to be equal to the apparent magnitude appearing from a distance of 10 parsecs. The bolometric absolute magnitude of a celestial object M , which takes into account the electromagnetic radiation on all wavelengths, is defined as M − M = −2.5 log(L /L ), where L and L are the luminosity of the object and of the Sun, respectively. The brightness profile of the galaxies µ(R) was derived in some works [43][44][45] from isophotal fits, employing the orientation parameters of the galaxies (center, inclination angle and ellipticity). This analysis leads to µ(R) which would be seen if the galaxy rotational axis was parallel to the line-of-sight. We used this µ(R) to generate S(R).
The surface photometry of the dwarf galaxies are consistent with modeling their baryonic component as an axisymmetric exponential disk with surface brightness [46]: where b is the scale length of the exponential disk, and S 0 is the central surface brightness. To convert this to mass density profiles, we fitted the mass-to-light ratio (Υ = M/L) of the galaxies. In Newtonian gravity, the rotational velocity squared of an exponential disk emerges as Freeman proposed [46]: with I and K the modified Bessel functions, evaluated at R/2b. In Yukawa gravity, a more cumbersome expression has been given in the work of De Araujo and Miranda [47] as where λ = h/m g /c = 2πR g is the Compton wavelength. For b/λ 1, the Newtonian limit is recovered.

Testing Pure Baryonic and Baryonic + Dark Matter Models
We chose 12 late-type dwarf galaxies from the Westerbork HI survey of spiral and irregular galaxies [43][44][45] to test rotation curve models. The selection criterion was that these disk-like galaxies have the longest R-band surface photometry profiles and rotation curves. For the absolute R-magnitude of the Sun, M ,R = 4.42 m [48] was adopted. Then, we fitted Equation (32) to the surface luminosity profile of the galaxies, calculated with Equation (31) from µ(R). The best-fit parameters describing the photometric profile of the dwarf galaxies are given in Table 1.
We derived the pure baryonic rotational curves by fitting the square root of Equation (33) to the observed rotational curves allowing for variable M/L. The pure baryonic model leads to best-fit model-rotation curves above 5σ significance level for all galaxies (the χ 2 -s are presented in the first group of columns in Table 1), hence a dark matter component is clearly required.
Then, we fitted theoretical rotation curves with contributions of baryonic matter and BEC-type dark matter with massless gravitons to the observed rotational curves in Newtonian gravity. The model-rotational velocity of the galaxies in this case is given by the square root of the sum of velocity squares given by Equations (29) and (33) with free parameters Υ, ρ (c) and R * . The best-fit parameters are given in the second group of columns of Table 1. Adding the contribution of a BEC-type dark matter component with zero-mass gravitons to rotational velocity significantly improves the χ 2 for all galaxies, as well as results in smaller values of M/L. The fits are within 1σ significance level in five cases (UGC3851, UGC6446, UGC7125, UGC7278, and UGC12060), between 1σ and 2σ in three cases (UGC3711, UGC4499, and UGC7603), between 2σ and 3σ in one case (UGC8490), between 3σ and 4σ in one case (UGC5986) and above 5σ in two cases (UGC1281 and UGC5721). We note that the bumpy characteristic of the BEC model results in the limitation of the model in some cases, the decreasing branch of the theoretical rotation curve of the BEC component being unable to follow the observed plateau of the galaxies (UGC5721, UGC5986, and UGC8490). The theoretical rotation curves composed of a baryonic component plus BEC-type dark matter component with massless gravitons are presented on Figure 1.  Table 1. Parameters describing the theoretical rotational curve models of the 12 dwarf galaxies. Best-fit parameters of the pure baryonic model in the first group of columns: central surface brightness S 0 , scale parameter b, M/L ratio Υ, along with the χ 2 of the fit. This model results in best-fit model-rotation curves above 5σ significance level for all galaxies. Best-fit parameters of the baryonic matter + BEC with massless gravitons appear in the second group of columns: M/L ratio Υ, characteristic density ρ (c) , distance parameter R * , along with the χ 2 of the fit and the respective significance levels. Constraints on the parameter m 2 /λ are also derived. In five cases, the fits χ 2 are within 1σ and marked as boldface. The fits are between 1σ and 2σ in three cases, between 2σ and 3σ in one case, between 3σ and 4σ in one case and above 5σ in two cases. Best-fit parameters of the baryonic matter + BEC with massive gravitons are given in the third group of columns only for the well-fitting galaxies: the range for R BEC and the upper limit on m g are those for which the fit remains within 1σ. Corresponding constraints on the parameter m/µ are also derived.

Pure Baryonic
Baryonic + BEC with m g = 0 Baryonic + BEC with m g > 0  We attempted to distinguish among galaxies to be included in well-fitting or less well-fitting classes based on their baryonic matter distribution. Several factors affect the goodness of the fits, as follows. The best-fit falls outside the 1σ significance level in the case of seven galaxies. Among these galaxies, UGC8490 and UGC5721 have (a 1 ) steeply rising rotational curve due to their centralized baryonic matter distribution (b < 0.5 kpc, v max > 50kms −1 ) with (a 2 ) long, approximately constant height observed plateau. Joint fulfilment of these criteria does not occur for the well-fitting galaxies, as b 0.5 kpc for them. The rest of the galaxies with best-fits falling outside the 1σ significance level have (b 1 ) slowly rising rotational curve due to their less centralized baryonic matter distribution (b > 0.5 kpc, v max < 50kms −1 ) with (b 2 ) short, variable height observed plateau, holding relatively small number of observational points (N ≤ 15, a small N lowers the 1σ significance level). The well-fitting galaxies do not belong to this group, as either they hold more observational points, or have a longer, approximately constant height observed plateau. We expect that for the galaxies not falling in the classes with baryonic and observational characteristics summarized by either properties (a 1 )-(a 2 ) or (b 1 )-(b 2 ) the BEC dark matter model represents a good fit. Finally, we note the galaxy UGC3711 represents a special case due to the lack of sufficient observational data. Although the shape of its rotational curve is very similar to that of the best-fitting galaxy, UGC12060, it is based on just six observational points, lowering the 1σ level. Its points also have smaller error bars, which increases the χ 2 . This results in the best-fit rotational curve of UGC3711 falling outside out the 1σ significance level.
Finally, we fitted the theoretical rotational curves given by both a baryonic component and a non-relativistic BEC component with massive gravitons, employing Yukawa gravity. The parameters Υ, ρ (c) and R * were kept from the best-fit galaxy models composed of baryonic matter + BEC with massless gravitons. The model-rotational velocity of the galaxies arises as the square root of the sum of velocity squares given by Equations (27) and (34) with free parameters R BEC and R g . Adding mass to the gravitons in the BEC model leads to similar performances of the fits.

Discussion and Concluding Remarks
We estimated the upper limit on the graviton mass, employing first the theoretical condition of the existence of the constant Λ, then analyzing the modelfit results of those five dwarf galaxies for which the fit of the BEC model with massive gravitons to data was within 1σ significance level.
Keeping the best-fit parameters ρ (c) , R * , we varied the value of R BEC and R g and calculated the χ 2 between model and data. The upper limit on the graviton mass m g has been estimated from the values of R g , for which χ 2 = 1σ has been reached. The results are given in Table 1. We plotted the theoretical rotation curves given by a baryonic plus a non-relativistic BEC component with massive gravitons with limiting mass in Figure 1. As shown in Table 2, the fit with the rotation curve data has improved the limit on the graviton mass in all cases. Table 2. Constraints for both the upper limit for the mass of the graviton (first from the existence of Λ, second from the rotation curves) and for the velocity-type and density-type BEC parameters (related to the mass of the BEC particle, scattering length and chemical potential) in the case of the five well-fitting galaxies. Comparing the theoretical rotation curves derived in our model with the observational ones, we found the upper limit to the graviton mass to be of the order of 10 −26 eV/c 2 . We also note that the constraint on the graviton mass imposed from the dispersion relations tested by the first three observations of gravitational waves, 7.7 × 10 −23 eV/c 2 [49], is still weaker than the present one.

ID
For the BEC, we could derive two accompanying limits: (i) first m 2 /λ has been constrained from the corresponding values of R * arising from the fit with the massless gravity model; and then (ii) m/µ has been constrained from the constraints derived for the graviton mass and our previous fits through Equations (19) and (20). These are related to the bosonic mass, chemical potential and scattering length, but only two combinations of them, a velocity-type quantitȳ and a density-type quantityρ were restricted, both characterizing the BEC. Their values are also given in Table 2 for the set of five well-fitting galaxies.
If the BEC consists of massive gravitons with the limiting masses m = m g determined in Table  2, the chemical potential µ and the constant characterizing the interparticle interaction λ can be determined as presented in Table 3. Table 3. Constraints on µ and λ assuming m = m g in case of the five well fitting galaxies. With this, we established observational constraints for both the upper limit for the mass of the graviton and for the BEC.