Radial oscillations of boson stars made of ultralight repulsive dark matter

We compute the lowest frequency radial oscillation modes of boson stars. It is assumed that the object is made of pseudo-Goldstone bosons subjected to a scalar potential that leads to a repulsive self-interaction force, and which is characterized by two unknown mass scales $m$ (mass of the particle) and $F$ (decay constant). First we integrate the Tolman-Oppenheimer-Volkoff equations for the hydrostatic equilibrium of the star, and then we solve the Sturm-Liouville boundary value problem for the perturbations using the shooting method. The effective potential that enters into the Schr{\"o}dinger-like equation as well as several associated eigenfunctions are shown as well. Moreover, we found that the large frequency separation, i.e. the difference between consecutive modes, is proportional to the square root of the mass of the star and the cube of the mass scale defined by $\Lambda \equiv \sqrt{m F}$.


Introduction
Since the pioneer work of F. Zwicky about the dynamics of the Coma galaxy cluster in the 30's [1], and the observations made by V. Rubin to determine the rotation curves of galaxies a few decades later [2], we are convinced that most of the non-relativistic matter in the Universe is dark, usually referred to as cold dark matter. In modern times current wellestablished data coming from many different sources confirm the existence of dark matter [3], although its nature and origin still remains a mystery. For a review on dark matter see [4,5], and for recent reviews on dark matter detection searches see [6,7,8].
Usually dark matter in the standard parametrization of the Big-Bang cosmological model is assumed to be made of weakly interacting massive particles, a conjecture which works very well at large (cosmological) scales (≥ M pc), but encounters several problems at smaller (galactic) scales, like the corecusp problem, the diversity problem, the missing satellites problem and the too-big-to-fail problem [9]. These problems may be tackled in the context of selfinteracting dark matter [10,11], as any cuspy feature will be smoothed out by the dark matter collisions. In addition, if dark matter consists of ultralight scalar particles with a mass m ≤ eV , and with a small repulsive quartic self-interaction a Bose-Einstein condensate (BEC) may be formed with a long range correlation. This scenario has been proposed as a possible solution to the aforementioned problems at galactic scales [12,13,14].
An excellent dark matter candidate with selfinteractions is the Quantum Chromodynamics (QCD) axion [15,16,17,18,19], a hypothetical spin 0 particle proposed by the Peccei-Quinn mechanism [20,21] to solve dynamically the so-called strong CP problem [22], which can be stated as "why is the θ parameter in QCD so small?". The QCD axion, well studied in the literature, is one of the good dark matter candidates [23], although it has not been detected yet. What is more, within String Theory [24,25], other ultralight axions with a wide range of masses, have been shown to be good dark matter candidates [26,27,28].
The QCD axion is a member of a more general class of particles, called pseudo-Goldstone bosons, which are very common in Particle Physics, and they arise in the following way: An additional global symmetry is introduced, which is spontaneously broken at some high scale F with soft explicit symmetry breaking at a lower scale Λ ≪ F . Since the additional symmetry is a spontaneously broken global symmetry, there must be a Goldstone boson associated with this symmetry. However, due to the soft explicit symmetry breaking, this boson acquires a small non-vanishing mass, given by m = Λ 2 /F , instead of being massless. QCD axions have been used in Cosmology to accelerate the Universe, either in a concrete inflationary scenario called Natural Inflation [29,30,31,32], or to describe the current cosmic acceleration [33,34]. Its self-interaction, however, is found to be attractive, while tackling the cusp/core problem, as already mentioned before, requires a repulsive force.
Axion stars are gravitationally bounded collections of axions. The possibility of the formation of axion stars was considered for the first time in [35]. Since then axion stars have been studied in [36,37,38,39,40,41]. Boson stars have been studied in [42,43,44,45,46,47,48,49,50,51]. See also [52,53] for Newtonian self-gravitating Bose-Einstein condensates. The maximum mass for bosons stars in non-interacting systems was found in [42,43], while in [44,45] it was pointed out that self-interactions can cause significant changes. In [48,49] the authors constrained the boson star parameter space using data from galaxy and galaxy cluster sizes.
It is well-known that the properties of stars, such mass and radius, depend crucially on the equationof-state, which unfortunately is poorly known. Given the recent advances in Helioseismology and Asteroseismology in general, studying the oscillations of stars and computing the frequency modes offer us the opportunity to probe the interior of the stars and learn more about the equation-of-state, since the precise values of the frequency modes are very sensitive to the thermodynamics of the internal structure of the star [54]. In this work we develop an identical strategy to study the internal properties of these hypothetical axions stars which are on the list of a future astronomical observational programs [55]. For previous works on radial oscillations of stars see [56,57,58,59,60,61,62,63,64] and references therein.
It is the goal of the present work to study radial oscillations of Dark BEC stars made of ultralight repulsive scalar particles identified as pseudo-Goldstone bosons. Our work is organized as follows: after this introduction, we present the equation-of-state in section two, and the equations for the hydrostatic equilibrium in the third section. The equations for the perturbations are discussed in section four, and our numerical results in the fifth section. Finally we conclude in the last section. We work in natural units in which c = = 1. In these units all dimensionfull quantities are measured in GeV, and we make use of the conversion rules 1m = 5.068 × 10 15 GeV −1 , 1kg = 5.610 × 10 26 GeV and 1sec = 1.519 × 10 24 GeV −1 [65].

Hydrostatic equilibrium 2.1 Equation-of-state
The perturbative Lagrangian of a relativistic real scalar field φ is given by where the scalar potential is of the form and where we consider renormalizable theories only, ignoring all higher order terms. In this work the scalar field is identified with any pseudo-Goldstone boson. The sign of the quartic self-interaction is taken to be positive since we assume a repulsive self-interaction for the Dark pseudo-Goldstone boson. Therefore, the model assumed here is characterized by two unknown mass scales, namely the mass of the scalar particle, m, as well as the decay constant, F ≫ m, arising from the spontaneous breaking of some global symmetry. Unfortunately, it turns out that it is not easy to obtain scalar field models with a tiny mass and a repulsive force within known Particle Physics, although some attempts have been made [66]. In the following, without relying on concrete Particle Physics models, we shall assume that this is possible, and we shall study radial oscillations of objects made of Dark pseudo-Goldstone bosons.
The above scalar potential combined with the Gross-Pitaevskii equation [67,68,69], also known as non-linear Schrödinger equation, leads to the following equation-of-state for the ultralight pseudo-Goldstone boson [40,66]: where the constant K is computed to be [40,66] where a new mass scale Λ ≡ √ mF has been introduced.
The interested reader can find more details on the properties of the Bose-Einstein condensate in the appendix in the end of the article.

Structure equations
We briefly review relativistic stars in General Relativity. The starting point is Einstein's field equations without a cosmological constant with G N being Newton's constant, which is the following will be set to unity. In the exterior of the star the matter energy momentum tensor T µν vanishes. For matter we assume a perfect fluid with pressure P , energy density ǫ and an equation of state P (ǫ). For the metric in the case of static spherically symmetric spacetimes, we consider the following ansatz ds 2 = −e A(r) dt 2 +e λ(r) dr 2 +r 2 (dθ 2 +sin 2 θdϕ 2 ), (6) with two unknown metric functions of the radial distance A(r) and λ(r). For the exterior problem r > R, with R being the radius of the star, one obtains the well-known Schwarzschild solution where M is the mass of the star. For the interior solution r < R, we introduce the mass function, λ(r) = − ln (1 − 2m(r)/r). The continuity of the model is guaranteed by matching the two solutions at the surface of the star, for which we obtain m(R) = M . The Tolman-Oppenheimer-Volkoff (TOV) equations for the interior solution [70,71], with a vanishing cosmological constant read and where the prime denotes differentiation with respect to r (i.e., m ′ ≡ dm/dr). The first two equations are to be integrated with the initial conditions m(r = 0) = 0 and P (r = 0) = P c , where P c is the central pressure. The radius of the star is determined requiring that the energy density vanishes at the surface, P (R) = 0, and the mass of the star is then given by M = m(R). Finally, the other metric function can be computed using the third equation together with the boundary condition A(R) = ln (1 − 2M/R).
To avoid confusion a couple of remarks are in order here. First, since the wavefunction of the Dark BEC star extends to infinity, the radius of the star is formally infinite. This is the reason why in the literature it is used the symbol R 99 (it contains the 99 % of the mass of the star) rather that R. In the following for simplicity we drop the subindex 99. Furthermore, in [72] in solving the Klein-Gordon equation coupled to Einstein's field equations for any scalar field potential the authors noticed that the boson star was not isotropic, since the corresponding radial pressure component p r and the tangential pressure component p t were different, p r = p t , and so there was a nonvanishing anisotropy p t − p r . In more recent works, however, bosonic dark matter inside the star was modelled as a Bose-Einstein condensate characterized by a polytropic EoS of the form p = Kρ 2 , and the standard TOV equations for isotropic fluids were considered [73,74]. This is possible within the Thomas-Fermi approximation where the fractional pressure anisotropy is small, see e.g. [46] for its application to self-gravitating bosons stars.

Radial oscillations 3.1 Equations for the perturbations
If δr is the radial displacement and δP is the perturbation of the pressure, the equations governing the dimensionless quantities ξ = δr/r and η = δP/P are the following [60,64,75] and η ′ (r) = ξ ω 2 r(1 + ǫ/P )e λ−A − 4P ′ (r) P −ξ 8π(P + ǫ)re λ − r(P ′ (r)) 2 P (P + ǫ) where e λ = (1 − 2m/r) −1 and e A are the two metric functions, ω is the circular frequency oscillation mode (with ν = ω/(2π) being the ordinary frequency), and γ is the relativistic adiabatic index defined to be Here, however, we prefer to work equivalently with a second order differential equation used in [63] − supplemented with the boundary conditions at the origin r = 0 and at the surface of the star r = R: ζ(r = 0) = 0 and δp(r = R) = 0. In the previous equation ζ = r 2 e −A/2 ξ and the background functions f , G an H are given by and finally the perturbation of the pressure can be computed as Therefore, contrary to the previous hydrostatic equilibrium problem, which is an initial value problem, this is a Sturm-Liouville boundary value problem, and as such the frequency ν is allowed to take only particular values, the so-called eigenfrequencies ν n . Each ν n corresponds to a specific radial oscillation mode of the star. Accordingly, each radial mode is identified by its ν n and by an associated eigenfunction ζ n . The Sturm-Liouville boundary value problem at hand can be treated equivalently as a quantum mechanical problem by recasting the second order differential equation for ζ into a Schrödinger-like equation [76] of the form where the new variables τ and ψ are defined as acoustic radius τ = r 0 f −1 (z)dz and ψ(τ ) = ζ/u. The effective potential is found to be where the function Π is given by Π = −f ′ −G/f while u is determined by the condition u ′ /u = −Π/(2f ).

Numerical results
We consider for the ultralight scalar boson a mass m = 10 −8 eV , and a decay constant of the order of the GUT scale, F ∼ 10 16 GeV , which correspond to a mass scale Λ ≡ √ mF = 330 M eV , and central pressures that correspond to a boson star with a mass M 1 = 1.86 × 10 −14 M ⊙ (light object) and M 2 = 9.3 × 10 −10 M ⊙ (heavy object) and radius R = 6.9 km. Therefore, the order of magnitude of the fundamental mode is expected to be [77] ω 0 ∝ M/R 3 = 2.74 mHz for the light star and ω 0 ∝ 0.62 Hz for the heavy one, while the spacing in frequencies ν n = ω n /2π for the highly excited modes is computed by [78,79] where c s is the sound speed given by c 2 s = dP/dǫ. Conveniently, ν 0 reads where a o is a numerical constant, such that a o = π 0 z/sinz dz ≈ 5.8993. In particular, we notice that ν 0 is proportional to Λ 3 and it scales as ∼ M 1/2 . For a given mass one can easily compute the constant spacing at higher excited modes, and for the masses considered here we find ν 0 = 0.73 mHz for the light star and ν 0 = 0.16 Hz for the heavy star.
We computed the lowest oscillation modes summarized in tables 1 and 2. We see that the fundamental mode is of the expected order of magnitude. Introducing τ 0 = τ (R) = 11.4 min for the light object and τ 0 = 3.07 sec for the heavy one, the effective potential U that enters into the Shrödinger-like equation versus dimensionless τ /τ 0 is shown in Figures 1 and 2 for the light and the heavy star, respectively, while the corresponding eigenfunctions of several modes are shown in Figures 3 and 4. We recall that in a Sturm-Liouville boundary value problem the number of zeros of the eigenfunctions corresponds to the overtone number n, namely the first excited mode corresponds to n = 1 and has only one zero, the second excited mode corresponds to n = 2 and has two zeros, while the fundamental mode does not have zeros and corresponds to n = 0. In Fig. 3 we show the eigenfunctions ψ corresponding to the fundamental mode, n = 0 (blue), the first two excited modes, n = 1 (orange), and n = 2 (red), and the last two (9th excited:brown and 10th excited:magenta). In Fig. 4 we show the eigenfunctions corresponding to the fundamental and to the first four excited modes. Finally, in Figures 5 and 6 we show the large frequency separation ∆ν n = ν n+1 − ν n versus frequen- cies themselves for Λ = 330 M eV for the light and the heavy star, respectively. At higher order n the spectrum exhibits a series of equally spaced modes, precisely as expected.
Before we finish a couple of remarks are in order about radial and non-radial oscillations of the boson stars considered here. First, regarding detection and oscillations of Dark BEC stars, when these objects collide with neutron stars, or with other objects with a strong magnetic field, depending on how ordinary matter interacts with the ultralight pseudo-Goldstone bosons, fast radio bursts may be emitted in the surface of the neutron star. Such a possibility has been pointed out in the literature for axion stars [80]. In such collisions, the tidal forces of the neutron star may produce large variations on the structure of these light Dark BEC stars (that stretch in the direction of the neutron star), which will excite its radial and non-radial oscillations modes [80]. These perturbations in the boson stars will then lead to significant variations in their electromagnetic radiation. From the event rate in a galaxy we can determine the mass of the boson star [81]. Furthermore, current and future radio telescopes on Earth such as Arecibo or SKA, and the next generation of stellar missions such as "PLAnetary Transits and Oscillation of stars" (PLATO) [82], have the potential to put important constraints in their properties, or find such stars if they actually exist.
On the other hand, it would be useful if there were other complementary ways to test the existence of the boson stars and probe their interior. To this end, the Love number formalism can be of great help to characterize the deformation of the star due to the gravitational forces of the companion object. In binary systems the tidal forces may have observable effects on the emission of electromagnetic and gravitational radiation. In particular, the quadrupolar deformation of the star, Q ij , is proportional to the external tidal tensor E ij [83,84] where λ is the tidal de-formability, and k 2 is the Love number. These two parameters depend entirely on the EoS of the star, and the Love formalism has been proven to be a powerful tool to discriminate between regular black holes, for which k 2 = 0, and other alternatives, such as Exotic Compact Objects [85]. In the Newtonian limit, P ≪ ǫ, β → 0, with β = M/R being the compactness of the star, and for a polytropic power γ = 2 (polytropic index unity) the Love number is computed to be k 2 = (π 2 −5)/(2π 2 ) ≃ 0.26 [83], which agrees with the tabulated values of [86]. This value for a BEC Dark star (more general for polytropes with n = 1) is larger than the highest one for neutron stars with realistic equations-of-state, ∼ 0.14 [83]. Moreover, it has been proposed that very light dark matter clumps may be detected by Pulsar Timing Arrays [87,88], while space-based detectors, such as eLISA, will probe gravitational waves (GWs) in the mHz range [89,90]. The amplitude of the GW from such light objects is expected to be extremely low. If, however, a huge number of boson stars pulsates simultaneously, the amplitude of the stochastic gravitational wave background may be significantly enhanced.

Conclusions
We have studied radial oscillations of Dark BEC stars made of ultralight repulsive scalar particles. Although it is highly non-trivial to obtain a scalar potential with a tiny mass and a repulsive selfinteraction, we have assumed in this work that this is possible without relying on concrete Particle Physics models. The model is characterized by two free, unknown mass scales, namely the mass of the scalar boson and the decay constant arising from the spontaneous breaking of a global symmetry. First we solved the Tolman-Oppenheimer-Volkoff equations for the hydrostatic equilibrium of the star, and then using the known background solutions we solved the Sturm-Liouville boundary value problem for the perturbation with the shooting method. We have computed the fundamental as well as several excited modes for two different star masses, and we have shown graphically i) several eigenfunctions corresponding to the first three and two highly excited oscillation modes, and ii) how the large frequency difference varies with the frequencies themselves. In addition, we have reformulated the boundary value problem equivalently by writing down a Schrödinder-like equation, and we have shown the effective potential together with the first five values of ω 2 .

Acknowlegements
We are grateful to the anonymous reviewer for a careful reading of the manuscript, for his/her constructive criticism, and useful comments and suggestions. It is a pleasure to thank K. Clough, D. Hilditch and C. Moore for enlightening discussions. The authors thank the Fundação para a Ciência e Tecnologia (FCT), Portugal, for the financial support to the Center for Astrophysics and Gravitation-CENTRA, Instituto Superior Técnico, Universidade de Lisboa, through the Grant No. UID/FIS/00099/2013.

A A cold dilute boson gas
We model the collection of ultralight scalar particles in the Dark BEC star as a cold dilute boson gas.
Under these conditions only binary collisions at low energy are relevant, and thus they are characterized by the s-wave scattering length l irrespectively of the details of the two-body potential. Therefore, we can replace the full potential by a short range repulsive delta-potential of the form [67] which implies a self-interaction cross section given by σ = 4πl 2 . The ground state properties of the boson gas are described by the mean-field Gross-Pitaevskii equation, also known as non-linear Schrödinger equation i ∂Ψ(t, r) ∂t = − ∇ 2 2m + m U + U 0 |Ψ(t, r)| 2 Ψ(t, r), and where we have introduced the quantum potential V Q = −(1/2m)∇ 2 √ ρ/ √ ρ, while ρ can be viewed as the density of the fluid. Almost all particles are in the condensate, the effective pressure of which is computed to be P = 2πl In the case, the scattering length, the mass of the particles and the scale Λ are related via [40,66] l m 3 = 1 32π(F m) 2 = 1 32πΛ 4 , and therefore we obtain for the constant K the final expression K = 1/(2Λ) 4 .