Boson Stars from Self-Interacting Dark Matter

We study the possibility that self-interacting bosonic dark matter forms star-like objects. We study both the case of attractive and repulsive self-interactions, and we focus particularly in the parameter phase space where self-interactions can solve well standing problems of the collisionless dark matter paradigm. We find the mass radius relations for these dark matter bosonic stars, their density profile as well as the maximum mass they can support.


Introduction
Bosonic degrees of freedom arise generically and naturally in theories of fundamental physics, both in the Standard Model and beyond. The Higgs boson is of paramount importance, being the only fundamental scalar in the Standard Model [1,2], but many other scalar degrees of freedom have been proposed to extend particle physics to high energy scales. These include (among many others) the axion of QCD [3] or the scalar which drives the expansion of the universe in quintessence models [4].
These bosonic particles often make good Dark Matter (DM) candidates as well. One reason for this is that unlike the Higgs, many of these new scalars would be stable or long-lived enough that they could coalesce into DM halos which constitute the seeds of galaxy formation. Unlike the usual collisionless cold DM picture, however, we are interested in the scenario where large collections of these bosons form bound states of macroscopic size due to their self-gravitation (and self-interaction generically). For this picture to be consistent, the scalars are taken to be sufficiently cold so that they may coalesce into a Bose-Einstein Condensate (BEC) state, and can thus be described by a single condensate wavefunction. These wavefunctions can indeed encompass an astrophysically large volume of space and have thus been termed "boson stars" [5].
It was shown many years ago that objects of this type are allowed by the equations of motion, first by Kaup [6] and subsequently by Ruffini and Bonazzola [7] in non-interacting systems. They found a maximum mass for boson stars of the form M max ≈ 0.633M 2 P /m, where M P = 1.22 × 10 19 GeV is the Planck mass and m is the mass of the individual bosons. (This is very different from the analogous limit for fermionic stars, termed the Chandrasekhar limit, which scales as M 3 P /m 2 ). Later, it was shown by Colpi et al. [5] that self interactions in these systems can cause significant phenomenological changes. In particular, they examined systems with repulsive self-interactions, and show that the upper limit on the mass is M max ≈ 0.02 √ λM 3 P /m 2 , where λ is a dimensionless φ 4 coupling. 1 This extra factor of M P /m as compared to the noninteracting case makes it more plausible that boson stars can have masses even larger than a solar mass. A different method of constraining the boson star parameter space, which fits the coupling strength using data from galaxy and galaxy cluster sizes, has been considered in [8,9].
The situation for attractive self-interactions is slightly more complex. The simplest case involves a self-interaction of the form λφ 4 , where λ < 0 for attractive interactions. If this were the highest-order term in the potential, then it would not be bounded below, and so one typically stabilizes it by the addition of a positive φ 6 term. We will assume that the contribution of such higher-order terms is negligible phenomenologically (we address the validity of that assumption in Section 3.3). Furthermore, in this scenario the typical sizes of gravitationally bound BEC states is significantly smaller than the repulsive or non-interacting cases. This is because the only force supporting the condensate against collapse comes from the uncertainty principle. Gravity and attractive self-interactions tend to shrink the condensate. We will see in Section 3 that the maximum mass for an attractive condensate scales as M max ∼ M P / |λ|. This result was originally derived using an approximate analytical method [10], and was later confirmed by a precise numerical calculation [11].
DM self-interactions have already been proposed and studied in different contexts [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]. One of the main reasons why DM self-interactions can play an important role is due to the increasing tension between numerical simulations of collisionless cold DM and astrophysical observations, the resolution of which (for the moment) is unknown. The first discrepancy, known as the "cusp-core problem", is related to the fact that dwarf galaxies are observed to have flat density profiles in their central regions [32,33], while N-body simulations predict cuspy profiles for collisionless DM [34]. Second, the number of satellite galaxies in the Milkly Way is far fewer than the number predicted in simulations [35,36,37,38,39,40]. Last is the so-called "too big to fail" problem: simulations predict dwarf galaxies in a mass range that we have not observed, but which are too large to have not yet produced stars [41].
The solution of these problems is currently unknown, but a particularly well-motivated idea involves self-interacting DM (SIDM). Simulations including such interactions suggest that they have the effect of smoothing out cuspy density profiles, and could solve the other problems of collisionless DM as well [42,43,19]. These simulations prefer a self-interaction cross section of 0.1 cm 2 /g σ/m 10 cm 2 /g. There are, however, upper bounds on σ/m from a number of sources, including the preservation of ellipticity of spiral galaxies [44,45]. The allowed parameter space from these constraints nonetheless intersects the range of cross sections which can resolve the small-scale issues of collisionless DM, in the range 0.1 cm 2 /g σ/m 1 cm 2 /g.
Self-gravitation and additionally extra self-interactions among DM particles can lead in some cases to the collapse of part of the DM population into formation of dark stars. The idea of DM forming star-like compact objects is not new. Dark stars that consist of annihilating DM might have existed in the early universe [46,47,48]. Dark stars have been also studied in the context of hybrid compact stars made of baryonic and DM [49,50,51,52] as well as in the context of mirror DM [53,54,55,56]. Additionally some of the authors of the current paper studied the possibility of dark star formation from asymmetric fermionic DM that exhibits Yukawa type self-interactions that can alleviate the problems of the collisionless cold DM paradigm [57]. Unlike the dark stars of annihilating DM, asymmetric dark stars can be stable and observable today. [57] displays the parameter space where it is possible to observe such dark stars, providing mass radius relations, corresponding Chandrasekhar mass limits and density profiles. Self-interactions in dark stars have also been considered in [58] for fermionic particles, as well as in [59] for bosonic ones.
In this paper we examine the dark stars composed of asymmetric self-interacting bosonic DM. The study is fundamentally different from that of [57] because unlike the case of fermionic DM where the stability of the star is achieved by equilibrium between the Fermi pressure and gravitation, bosonic DM does not have a Fermi surface. They form a BEC in the ground state and it is the uncertainty principle that keeps the star from collapsing. We are going to demonstrate how DM self-interactions affect the mass radius relation, the density profile and the maximum mass of these DM bosonic stars in the context of the self-interactions that reconcile cold DM with the observational findings.
Note that we set = c = 1 in what follows.

SIDM parameter space
As we mentioned above, galactic scale N -body simulations of cold, non-interacting DM indicate that the central regions of galaxies should have a "cuspy" density profile, contrary to the cored profiles one observes. This, along with the "missing satellites" and "too big to fail" problems, has led some to question the non-interacting DM paradigm. While some believe that the inclusion of baryonic physics could alleviate these issues [60,61,62,63], it remains an open question. On the other hand, the inclusion of self-interactions in the DM sector could resolve these issues without creating tension with other astrophysical constraints. These two conditions can be simultaneously satisfied if the cross section per unit mass for DM satisfies Assuming a velocity independent cross section, [19] found that σ/m = 1 cm 2 /g tends to over-flatten dwarf galaxy cores and that it is marginally consistent with ellipticity constraints of the Milky Way. On the other hand a value of 0.1 cm 2 /g satisfies all constraints and flattens dwarf galaxy cores sufficiently. Let us consider a potential of the form Note that λ > 0 (λ < 0) signifies a repulsive (attractive) interaction. The resulting DM-DM scattering cross-section is at tree level. Plugging this into Eq. (1), we get the constraint This matches the results of [64]. For perturbativity, we should restrict λ 4π, which would imply that our results are valid only for m 100 MeV. In this mass range, it is plausible that these DM particles coalese into boson stars at some point in early cosmology. If a large fraction of DM is contained inside boson stars, the derived parameter space may be significantly altered [65], since boson star-DM interactions and boson star self-interactions may become significant. We will however assume that boson stars are rather scarce and the DM self-interactions are dominated by DM-DM scattering.

DM scattering with boson stars
To quantify how scarce boson stars have to be within this approximation, we assume that boson stars have a characteristic radius R, mass M and number density n BS . The mass, number density and selfinteraction cross section of free DM is taken to be m, n and σ. The mean free path a DM particle travels before hitting another DM particle or a boson star will be λ DM = (nσ) −1 and λ BS ∼ (n BS πR 2 ) −1 , respectively. Scattering with boson stars has to be much rarer than with other free DM in our approximation. Therefore we require λ DM λ BS . For the DM density we use the typical value of the solar system, i.e. ρ DM = M n BS + mn ≈ 0.3GeV/cm 3 . These requirements lead to the following condition Taking self-interactions to be that of Eq. (3), and the boson star radius to be comparable to the minimum radius (which scales the same for both signs of interaction) R ∼ |λ|M P /m 2 (see Eq. (26)), Eq. (5) becomes The maximum mass of a boson star with non-negligible attractive interactions is ∼ M P / |λ|. Since this scaling is only proportional to a single power of M P , the first term in the denominator of Eq. (6) tends to dominate. We obtain in the attractive scenario where AU is an astronomical unit. The minimum mean distance between attractive boson stars can therefore within this approximation be (n att BS ) −1/3 ≈ 40(|λ|m/MeV) −1/3 AU. In the scenario with repulsive interactions the maximum mass scales as √ λM 3 P /m 2 . Therefore the second term in the denominator of Eq. (6) dominates. The number density must satisfy The minimum mean distance between repulsive boson stars which leaves our approximation valid can at most be (n rep BS ) −1/3 ≈ 5 × 10 2 λ 2/3 (m/MeV) −2/3 pc.

Bosonic Dark Matter
An important property of light scalar particles that has been examined extensively in the literature [66,67] is that large collections (particle number N 1) can transition to a BEC phase at relatively high temperature, as compared to terrestrial experiments with cold atoms. The critical temperature for condensate occurs when the de Broglie wavelength is equal to the average interparticle distance, λ dB = [ζ(3/2)/n] 1/3 , where n is the average number density of the particles and ζ(x) is the Riemann Zeta function. This implies a critical temperature for transition to the BEC phase of the form In this paper, we will assume that all relevant scalar field particles are condensed, i.e. that the system is in its ground state, a perfect BEC. The effect of thermal excitations is examined in [68] and they are expected to be negligible as long as T < T c is satisfied.

Non-Interacting Case
It is instructive to begin with the case of boson stars bound only by gravity, first analyzed in [6]. In this seminal work, Kaup considers the free field theory of a complex scalar in a spacetime background curved by self-gravity. The equations of motion 2 were solved numerically. The maximum mass of these solutions was found to be M max ≈ 0.633M 2 P /m, the oft-quoted Kaup limit for non-interacting boson stars. This value was later confirmed by Ruffini and Bonazzola [7], who used a slightly different method by taking expectation values of the equations of motion in an N -particle quantum state.
Interacting field theories are more complex. In particular, for cross sections satisfying Eq. (1), the phenomenology of repulsive and attractive interactions are very different, and accordingly, the methods required to analyze them are different as well. We outline the relevant methods in the sections below.

Repulsive Interactions
If the self-interaction is repulsive, we can make use of the result of Colpi et al. [5]. Like Kaup, their method begins with the relativistic equations of motion for a boson star, the coupled Einstein and Klein-Gordon equations, but including a self-interaction term represented by Λ: where the rescaled variables are x = mr, σ = √ 4πGΦ (Φ the scalar field), Ω = ω/m (ω the particle energy), and Λ = λM 2 P /(4πm 2 ). In addition to the scalar field itself, A(r) and B(r) must be solved for; these represent the deviations from the flat metric due to the self-gravity of the condensate, In practice, one can trade the metric function A(r) for the mass M(x) by the relation In the limit that the interactions are strong (precisely, Λ 1), the system can be simplified significantly, as one can perform a further rescaling of the equations: σ * = σΛ 1/2 , x * = xΛ −1/2 , and M * = MΛ −1/2 . The relevant parameters of Section 2 suggest a value of Λ = O(10 40 ) or higher, so it is completely safe to neglect terms proportional to Λ −1 . In this limit the equations simplify to 2 The non-interacting equations of motion are equivalent to Eq.s (10) and (11) in the limit Λ → 0.
where the pressure p * and density ρ * are given by In this limit, the equations do not depend on Λ, and one finds numerically that there is a maximum (dimensionless) mass M * max ≈ 0.22. Restoring the appropriate dimensions, one finds This bound on the mass of repulsive boson stars was confirmed very precisely using a hydrodynamic approach as well [69]. Figures 1 and 2 show the mass-radius relation and selected density profiles, respectively. The branch to the left of the peak in Figure 1 represents unstable equilibria, where the ground state energy is higher than the equilibrium on the right branch with the same number of particles (and thus the same quantum numbers).
where M = 1.99 × 10 30 kg is the solar mass. The range of masses allowed by these inequalities are represented in Figure 3. Because of the strength of the repulsive interactions, these solutions can have masses several orders of magnitude above M . If there is a significant number of such objects in the Milky Way, it could have important observational signatures. However, a detailed analysis of the formation of these objects is required, in order to give some indication of whether DM boson stars in galaxies have masses close to the maximum value or lower.

Attractive Interactions
If DM self-interactions are attractive, then the method of [5] does not apply. However, assuming relativistic corrections are negligible, one can instead solve the nonrelativistic equations of motion numerically where V is the trapping potential, which in our case is the gravitational potential of the BEC and satisfies the Poisson equation The s-wave scattering length a is related to a dimensionless φ 4 coupling λ by a = λ/(32πm).
Here, ρ(r) = m · n(r) = m · |ψ(r)| 2 is the mass density of the condensate, which is normalized such that d 3 rρ(r) = M , the total mass. The three terms on the right-hand side of Eq. (16) correspond to the kinetic, gravitational, and self-interaction potentials, respectively. As our notation signifies, we will assume that the density function is spherically symmetric, i.e. ρ( r) = ρ(r), which should be correct for a ground state solution.
Because the Gross-Pitäevskii + Poisson system (hereafter GP, defined by Eqs. (16) + (17)) cannot be solved analytically in general, we use a shooting method to integrate the system numerically over a large range of parameters. As boundary conditions, we choose the values of ψ(0) and V (0) so that both functions are regular as r → 0, and so that asymptotically ψ(r) → 0 and rV (r) → 0 exponentially as r → ∞. Some examples of integrated density functions are given in Figure 5. Our numerical procedure requires the following rescaling of the dimensionful quantities: where the dimensionless quantities on the RHS are denoted with a tilde. The equations take the form where∇ denotes a gradient with respect tor, and we have explicitly taken a < 0. These are the equations we solve. Similar rescaled equations were used in [71], but for repulsive interactions, and unlike [71], we also scale away the scattering length a. This makes our solutions valid for any generic a < 0. In Figure 4 we show the mass-radius relation for the bosonic stars, which agrees well with the results obtained in [11]. As in the repulsive case, there is a maximum mass for these condensates, but this mass is significantly smaller for attractive interactions. For parameters satisfying Eq. (4), our analysis shows that condensates of this type would be light and very dilute, having masses < 1 kg and radii R ∼ O(km). (Our assumption that the General Relativistic effects could be neglected in this case is therefore well supported a posteriori.) One can arrive at a good, order of magnitude analytic estimate on the size and mass of condensates by a variational method which minimizes the total energy. To this end, we follow the approach of [10] by using the GP energy functional, As input, we choose an ansatz for the wavefunction ψ(r), and subsequently compute the energy of the condensate by integrating Eq. (20) up to some maximum size R. Minimizing the energy with respect to R should give a good estimate for the size of stable structures. Note that the gravitational potential V (r) must be chosen self-consistently to satisfy Eq. (17) for a given choice of ψ(r).
In order to illustrate the salient features of the method, we will choose a simple ansatz for the wavefunction: which is normalized as above. Performing the energy integral gives the result where A ≡ 1/(2m) and B ≡ 6Gm 2 /5. Minimizing E(R) with respect to R gives two critical points In this calculation, a natural length scale X ≡ A/B emerges. For any a = 0 (repulsive or attractive), the minimum of the energy lies at the solution with the "+" sign, i.e.
In the case of attractive interactions, there is a critical number of particles N max ≡ X/(9|a|), above which the real energy minimum disappears and no stable condensate exists. Using M max = mN max , this analysis sets a value for the maximum mass for stable condensates with attractive interactions: The corresponding limit on the radius is a lower bound, attractive boson stars being stable only for Note that while the coefficient depends on the details of the wavefunction ansatz, the scaling relations M att max ∼ M P / |λ| and R att min ∼ √ λM P /m 2 are completely generic. Using Eq. (4), we find The range of masses allowed by these inequalities is given by the green band in Figure 6. We plot the maximum masses over many orders of magnitude, between 1 eV and 1 GeV, but the maximum mass of boson stars with such strong attractive self-interactions is still < 1 kg. Note that the numerical results agree well with the predictions of the variational method to within an order of magnitude, even for the naïve constant density ansatz in Eq. (21). These estimates can be improved further by a more robust ansatz for the wavefunction.
As an example of a physical model, field theories describing axions exhibit an attractive self-coupling through the expansion of the axion potential V (A) = m 2 f 2 1 − cos(A/f ) , where A is the axion field, m is the axion mass, and f is the axion decay constant. Gravitationally bound states, particularly in the context of QCD axions, have become the topic of much recent interest [72,73,74]. These states typically have maximum masses of roughly 10 −11 M , far below the bounds set in this section, because the couplings are typically many orders of magnitude smaller.
As we pointed out in the introduction, in the case of attractive interactions the potential is unbounded from below since λ < 0. Therefore there must exist higher dimensional operators suppressed by some cutoff. The first irrelevant operator with a Z 2 symmetry is φ 6 /µ 2 c where µ c is the cutoff scale. We will now set a lower limit for µ c by requiring that the φ 6 term is negligible with respect to the φ 4 term for typical boson star field values. Assuming that the kinetic energy of the field is negligible, the energy density is roughly equal to the potential. The maximum mass and minimum radius in Eqs. (25) and (26) can also be used to estimate the energy density as ρ ≈ M max /R 3 min ≈ m 6 /|λ| 2 M 2 P . Now we can estimate the field valueφ inside the boson star with attractive interactions to be Requiring |λ|φ 4 φ6 /µ 2 c we obtain the inequality µ c m/|λ|.  4), as a function of DM particle mass m. The green band is the region consistent with solving the small scale problems of collisionless cold DM. The blue region represents generic allowed interaction strengths (smaller than 0.1 cm 2 /g) extending up to the Kaup limit which is shown in black. The red shaded region corresponds to λ 4π. Note that the horizontal axis is measured in grams.

Conclusions
In this paper we studied the possibility that self-interacting bosonic DM forms stars. We assumed that self-interactions are mediated by a λφ 4 interaction and we investigated what type of stars can be formed in the case of both attractive and repulsive self-interactions, giving particular emphasis to the parameter phase space of masses and couplings where the DM bosons alleviate the problems of collisionless DM. We have considered DM particles that populate the BEC ground state. We estimated the maximum mass where these dark stars are stable, the mass-radius relation and the density profile for generic values of DM mass and self-interacting coupling λ. We leave several things for future work. The first and most important is the mechanism of formation for these bosonic dark stars. Sufficiently strong self-interactions can lead to the gravothermal collapse of part or the whole amount of DM to dark stars [75]. In this case, DM self-interactions can facilitate the formation of bosonic stars because DM particles get confined to deeper self-gravitating wells simply by expelling high energetic DM particles out of the core. As the core loses energy, the virial theorem dictates that the core shrinks and heats up the same time. This leads to further energy loss and thus to the gravothermal collapse. Such a scenario could also explain why the black hole at the center of the Galaxy is so heavy, since DM bosonic stars could provide the initial seed required for the further growth of the supermassive black hole [76]. It is interesting to note that boson stars can coexist in equilibrium with black holes, as shown in [77,78]. One should also notice that if the whole density of DM collapses to dark stars, one does not have to be within the narrow band of parameter space depicted in Figures 3 and  6. Another possibility is the creation of high DM density regions due to adiabatic contraction, caused by baryons [79,80]. Moreover, bosonic DM particles can get trapped inside regular stars via DM-nucleon collisions. The DM population is inherited by subsequent white dwarfs that, in case of supernovae 1a explosions, can leave the bosonic matter intact, either alone or with some baryonic matter [81].
Asymmetric bosonic dark stars where no substantial number of annihilations take place will not be very visible in the sky, although present. Gravitational lensing could be one way to deduce the presence of such stars in the universe. Additionally, if the DM boson interacts with the Standard Model particles via some portal (e.g. kinetic mixing between a photon and a dark photon), thermal Bremmstrahlung could potentially produce an observable amount of luminosity. This is particularly interesting since such a photon spectrum would probe directly the density profile of the boson star. Bosonic stars could also disguise themselves as "odd" neutron stars. For example, it is hard to explain sub-millisecond pulsars with typical neutron stars. XTE J1739-285 could possibly be such a case, since it allegedly rotates with a frequency of 1122Hz [82]. Compact enough bosonic stars would have no problem to explain such high rotational frequencies. Another possibility is the observation of compact stars with masses higher than the maximum mass a neutron star can support. Such might be the case of the so-called "black widow" PSR B1957+20, with a mass of 2.4 solar masses [83]. Therefore, abnormal neutron stars can well be the smoking gun for the existence of asymmetric dark stars either with fermionic constituents like [57], or with the bosonic ones studied here.