Stirred and shaken: dynamical behavior of boson stars and dark matter cores

Bosonic fields can give rise to self-gravitating structures. These are interesting hypothetical new"dark matter stars"and good descriptions of dark matter haloes if the fields are very light. We study the dynamical response of Newtonian boson stars (NBS) when excited by external matter (stars, planets or black holes) in their vicinities. Our setup can describe the interaction between a massive black hole and the surrounding environment, shortly after the massive body has undergone a"kick", due to the collapse of baryonic matter at the galactic center, or dark matter depletion as a reaction to an inspiralling binary. We perform the first self-consistent calculation of dynamical friction acting on moving bodies in these backgrounds. Binaries close to coalescence"stir"the NBS core, and backreaction affects gravitational waveforms at leading $-6PN$ order with respect to the dominant quadrupolar term; the coefficient is too small to allow detection by next-generation interferometers. We also show that the gravitational collapse to a supermassive black hole at the center of a NBS is accompanied by only a small change in the surrounding core. The NBS eventually gets accreted, but for astrophysical parameters this occurs only after several Hubble times.

Bosonic fields can give rise to self-gravitating structures. These are interesting hypothetical new "dark matter stars" and good descriptions of dark matter haloes if the fields are very light. We study the dynamical response of Newtonian boson stars (NBS) when excited by external matter (stars, planets or black holes) in their vicinities. Our setup can describe the interaction between a massive black hole and the surrounding environment, shortly after the massive body has undergone a "kick", due to the collapse of baryonic matter at the galactic center, or dark matter depletion as a reaction to an inspiralling binary. We perform the first self-consistent calculation of dynamical friction acting on moving bodies in these backgrounds. Binaries close to coalescence "stir" the NBS core, and backreaction affects gravitational waveforms at leading −6P N order with respect to the dominant quadrupolar term; the coefficient is too small to allow detection by next-generation interferometers. We also show that the gravitational collapse to a supermassive black hole at the center of a NBS is accompanied by only a small change in the surrounding core. The NBS eventually gets accreted, but for astrophysical parameters this occurs only after several Hubble times.
Introduction. The nature and properties of dark matter (DM) are arguably among the most important open issues in science. Interesting candidates for DM include light bosonic fields, of which the axion is a prototypical example [1], or generalizations thereof such as axionlike particles [2,3], ubiquitous in string-inspired scenarios [4,5]. Due to moduli compactifications, the spectrum of these particles could be populated uniformly down to the Hubble scale, m H ∼ 10 −33 eV. The phenomenology of bosons with masses eV can be dramatically different from that of weakly-interacting massive particles at the GeV-TeV scale. Ultralight bosons can form macroscopic Bose-Einstein condensates which provide a natural alternative to the standard structure formation through DM seeds and to the cold DM paradigm [6][7][8][9]. It has been recently recognized that ultralight boson fields with masses of the order of 10 −22 eV are a compelling candidate for cold DM [9][10][11][12][13]. A similar, albeit much wider, phenomenology arises in models of ultralight vector fields, such as dark photons, also a generic prediction of string theory [14]. All experiments probing the interaction between DM and standard model particles have so far been inconclusive. Given that DM interacts gravitationally, and that practically all galaxies are permeated in this elusive matter, it is natural to study further its gravitational interaction. This prospect is made all the more attractive with the advent of gravitational-wave (GW) astronomy [15][16][17].
Light bosons can clump to form self-gravitating, stellar-mass or supermassive "boson stars" [18][19][20]. The study of the dynamics of such objects is important for a number of reasons, ranging from stability to the way they interact with surrounding bodies (stars, BHs, etc) [21]. It is crucial to understand scalar and GW emission, and possible constraints imposed via observations. Along these lines, we have in mind the understanding of local changes in the density triggered by the presence of a massive BH or star, the drag exerted by the bosonic clump on stars moving within it, the flux of energy and momentum induced by coalescing binaries, etc. These issues are relevant for large scale DM structures and GW physics [21][22][23][24][25], but also from the perspective of the interaction between DM stars and neutron stars or black holes (BHs) [26]. These questions can also be of direct interest for theories with screening mechanisms, where new degrees of freedom -usually scalars -can be nonlinearly screened on some scales [27]. Such mechanisms give rise to nontrivial profiles for the new degrees of freedom, for which many of the questions above apply.
We study the response of localized scalar configurations to bodies moving in their vicinities. The objects themselves -Newtonian boson stars (NBSs) -have been studied for decades [18][19][20]. Despite this and the recent activity at the numerical relativity level [28][29][30][31][32][33][34], their interaction with smaller objects (describing, for example, stars piercing through or orbiting such NBSs) has hardly been studied. The variety and disparity of scales in the problem makes it ill-suited for full-blown numerical techniques, but ideal for perturbation theory.
We use units where the speed of light, Newton's constant and reduced Planck's constant are all set to unity, c = G = = 1. Setup: background. We consider a general U (1)invariant, self-interacting, complex scalar field Φ(x µ ) minimally coupled to gravity, described by the action where R is the Ricci scalar of the spacetime metric g µν , g ≡ det(g µν ). The scalar field stress-energy tensor For a spherically symmetric, stationary NBS where Ψ(r) is a real function, regular everywhere. One can perform the Newtonian limit of the theory, with a spacetime (|U (r)| 1) Assuming that the scalar Φ is non-relativistic, one finds the leading-order Schrödinger-Poisson equations of motion (∇ 2 is the Laplacian operator) The non-relativistic regime implies that Ω = µ − γ, in the limit 0 < γ µ. In this case, the energy Ω of the individual scalar particles forming the NBS is approximately given by their rest-mass energy µ. Remarkably, this system is left invariant under the transformation These relations are crucial: once a NBS solution is found, all others can be obtained through a rescaling of that solution. The numerical solution of system (3), with appropriate boundary conditions, is straightforward [35][36][37][38]. Accurate and simple fits for the potential and scalar are shown in an upcoming work [38]. At large distances, the scalar decays exponentially as Ψ ∼ e − √ 2µγr /r, whereas the Newtonian potential falls off as −M NBS /r. For fundamental NBSs one finds γ 0.162712M 2 NBS µ 3 , where M NBS is the NBS gravitational mass. All the fundamental NBSs satisfy the scaling-invariant mass-radius relation M NBS µ 9.1/(R NBS µ) where the NBS radius R NBS is defined as the radius of the sphere enclosing 98% of its mass [12,20,36]. We find the coupling GM NBS µ/(c ) = 7.5 × 10 −3 M NBS /10 10 M µ/10 −22 eV . Setup: perturbations. We are interested in the dynamical response to external perturbers, which disturb the spherically symmetric, stationary background above, with the assumption |δΨ| 1, where Ψ 0 is the radial profile of the undisturbed NBS. The fluctuation δΨ induces a change δT S µν in the stress-tensor which can be used to calculate physical quantities such as energy, linear and angular momenta radiated in a given process, obtained by computing the flux of certain currents, Here, ξ t = −∂ t , ξ ϕ = ∂ ϕ are timelike and spacelike Killing vector fields of the background spacetime, respectively. The limit is being taken at fixed retarded time.
The index i = {x, y, z} and e x , e y , e z are unit spacelike vectors in the x, y, z directions. Low-energy fluctuation modes, with energy ω 2 < (µ − Ω) 2 , are confined within the NBS. The interaction with an external perturber may deposit a considerable amount of energy and momentum in these modes. Together with fluxes at infinity, these determine, for example, the energy loss of the perturber or the dynamical friction it is subjected to. The deposition in the background object happens through excitation of its normal modes, which can be quantified once the dynamical equations for the perturbations have been established.
Non-relativistic perturbations of the form (5) and corresponding gravitational potential fluctuation δU satisfy the linearized system [38] i∂ t δΨ = − 1 2µ with the perturber treated as one or more pointlike sources, each of rest-mass m p , . We perform an expansion of δΨ, δU in spherical harmonics, and convert the equations to Fourier space with δΨ, δU ∼ e ±iωt . The end result is a set of coupled ordinary differential equations which can be robustly solved with standard techniques [38]. Note that this system is invariant under the re-scaling (U 0 , Ψ 0 , γ, ω) → λ 2 (U 0 , Ψ 0 , γ, ω), r → r/λ. Thus, one can map perturbations of a single NBS to the entire family of solutions. The perturbative scheme requires that |δΨ| 1, which can always be enforced by making m p as small as necessary. The background construction neglects higher-order PN contributions. A self-consistent perturbative expansion requires that such neglected terms ∼ U 2 0 do not affect the dynamics ∼ δU of small fluctuations, imposing Treating perturbers as pointlike is a successful and standard tool in BH perturbation theory [40][41][42], in seismology [43] or in calculations of gravitational drag by fluids [44,45]. In this approximation one loses smallscale information. For ultralight fields the smallest scale external to the pointlike particle is the Compton wavelength of the field, which is much larger than the size of stars, planets or BHs. In other words, we do not expect to lose important details of the physics at play. Free oscillations. The characteristic, non-relativistic oscillations of NBSs are regular solutions of the system (7)- (8) satisfying Sommerfeld conditions at large distances. The first few characteristic frequencies are listed in Table I. They are all normal mode solutions, with real frequencies, confined within the NBS. The characteristic frequencies cluster around γ. The fundamental monopolar mode in the Table agrees with one published result after proper normalization [39]. Modes of relativistic stars have been considered in the literature [46][47][48][49][50] and should smoothly go over to the numbers in Table I. An impurity at the center. Static perturbations of NBSs are interesting in their own right. For perturbers far away from the NBS center, the induced tidal effects can dissipate energy and lead to distinct signatures, both in GW signals and in the dynamics of objects close to such configurations [51][52][53]. The fluctuation in the density induced by a massive object sitting at the center is obtained in Ref. [38]. The particle attracts scalar field towards the center, where the gravitational potential corresponds solely to that of the point-like mass. We find an insignificant change in the local DM mass density, δρ M (0)/ρ M (0) ∼ 10 m p /M NBS . A BH eating its host NBS. On the other hand, studies of particle-like DM find that its density close to BHs increases significantly [54,55]. Such fact is in some tension with observations [10], but can be explained via scattering of DM by stars or BHs, or accretion by the central BH, induced by heating in its vicinities [56][57][58]. These results do not generalize to light scalars, at least when the configuration is spherically symmetric, since there are no stationary BH configurations with scalar "hair" [59,60]. If a (non-spinning) BH forms at the center of a NBS, the scalar is accreted by the BH, thus our previous results cannot be extrapolated to this situation, and describe the system only at intermediate times. What is the lifetime of such a system? When M BH M NBS , a perturbative calculation suffices. Consider a sphere of radius r + centered at the origin of the NBS. The NBS is stationary, and there is a flux of energyĖ in ≈ 10 −3 µ 7 r 2 + M 5 NBS crossing such a sphere inwards and outwards [38]. When a BH horizon exists at r + = 2M BH [61], a fraction of the incoming flux is absorbed by the BH. Low-frequency waves (µM BH 1) are poorly absorbed [62], and when ω ∼ μ 10 .
We have tested the above physics with a series of toy models, including a massive, non self-gravitating, scalar confined in a spherical cavity with a small BH at the center [38]. The models conform to the physics just outlined. WithĖ abs =Ṁ BH and fixed NBS mass, one finds the timescale where χ ≡ M NBS /M BH . In other words, the timescale for the BH to increase substantially its mass -a conservative indicative of the lifetime of the entire NBS -is larger than a Hubble timescale for realistic parameters. This timescale is the result of forcing the BH with the nearly monochromatic field of the NBS. When the NBS material is nearly exhausted, a new timescale is relevant, that of the quasinormal modes of a BH surrounded by a massive scalar ∼ M BH (M BH µ) −6 < τ , still typically larger than a Hubble time. When rotation is included, the entire setup may become even more stable; superradiance from the BH might even support nearly stationary, but non spherically-symmetric, configurations [63,64]. Plunges into NBSs. Baryonic matter tends to slowly accumulate near the center of a DM structure, where it may eventually collapse to a massive BH. Gravitational collapse can impart a recoil velocity v recoil to the BH, of the order of 300 Km/sec [65]. This will cause a BH to oscillate within the DM halo, on timescales of ∼ 10 6 yr 10 3 M pc −3 /ρ, and with an amplitude ∼ 69 pc v recoil 300 Km/s 10 3 M pc −3 /ρ, well within our nonrelativistic approximation. To model these processes, consider first a particle plunging head-on into a NBS, with velocity v = −v R e z (with v > 0) at the NBS surface. The particle crosses the NBS exactly once, and for small v R the process is always well-described by a Newtonian, non-relativistic approximation.
The motion triggers a flux of scalar particles at infinity. The radiation and momentum spectrum E rad , P rad are shown in Fig. 1 for different v R . The energy lost by the crossing particle, E lost , is also shown and is only a fraction of E rad , since the scalars have nonzero rest-mass. We find the following good description of our results accurate to within 5% of error for 0 v R 2.5M NBS µ. This interval spans over non-relativistic astrophysical rel-evant velocities (e.g 0 v R [km/s] 6000) for the DM core of the Milky Way.
The momentum lost by the plunging object is given by P lost = E lost /v R , and also gives us the dynamical friction exerted on the traveling body, dP/dt ∼ P lost v/(2R NBS ). Previous estimates, using an infinite, non self-gravitating scalar as an approximation [9,66], are plagued by near and far-distance divergences. By contrast, our results are finite and regular in the small-velocity limit, and are the first self-consistent calculation of the energy and momentum released when a massive body travels through a self-gravitating medium. Including the self-gravity of the scalar is a crucial ingredient for a correct calculation of dynamical friction (details to appear in Ref. [38]). At low velocities, our results predict a friction force one order of magnitude smaller than estimated without all the ingredients.
We can also calculate the relaxation timescale of a BH undergoing bound motion inside a NBS. For smallamplitude motion, we find that the amplitude A of the motion decays exponentially in tie as a consequence of dynamical friction, A = A 0 e −t/τs with the timescale τ s ∼ 10 10 yr 10 −22 eV µ This result described the timescale for a kicked BH (or star) to settle down at the center of an halo purely due to the dynamical friction caused by dark matter. We find that, excluding all but the interaction with the scalar, BHs of mass m p 10 5 M settle down in a timescale smaller than the Hubble time. This result can be compared with the timescale of damping due to dynamical friction caused by stars in the galactic core, estimated to be ∼ 0.1 τ s if we take the mass in the galactic core M c = M NBS [67]. Given the associated uncertainties, these results imply that DM needs to be taken into account when studying such processes. Low-energy binaries within NBSs. Consider now compact binaries within a DM core. These could be two BH binaries formed via core collapse or two neutron star binaries somewhere close to the galactic center. These systems have been observed recently via electromagnetic counterparts to GWs [68]. Binaries, either at an early or late stage in their life, stir the scalar field and produce disturbances in its profile. At high (ω orb γ, µU 0 ), but still non-relativistic frequencies (ω orb µ) equations (7) and (8) imply that |Ψ 0 δΨ| |δU |. Thus, equation (8) reduces simply to ∇ 2 δU = 4πP . In this regime, we find the following analytic solution, Here, M = 2m p for equal-mass binaries. The substitution (1 + (−1) m ) → 1 describes a single particle of mass m p revolving around a BH of mass M BH = M m p . The analytic expression above agrees with a full numerical solution within ∼ 2% in the entire regime of validity. At fixed orbital frequency the flux converges exponentially in l. Emission starts at ω = γ µ, and the coupling between gravity and the scalar implies that large frequency sources radiate less. The energy emitted in scalar waves by equal-mass binaries, with respect to their own quadrupole GW flux,Ė GW = (32/5)m 2 p M −2 (M ω orb ) 10/3 , reads aṡ T orb 16 yrs 9 2 , showing promising numbers, for low frequency emitters. High-energy binaries within NBSs. Rapidly moving binaries, such as those suitable for LIGO or LISA sources, do not fit into a non-relativistic regime. The relevant description of these systems for frequencies ω orb µ is [38] We find the following simple result for the flux in the dominant l = m modes, Again, (1 + (−1) m ) → 1 describes a single particle of mass m p revolving around a massive object.
To understand the impact of such energy loss on GW signals, we add to the quadrupole formula the energy flux (15), and compute the correction to the GW phase in the stationary approximation [69]. In Fourier space, we decompose the phase of the GW signalh(f ) = Ae iΥ(f ) as Υ(f ) = Υ (0) [1+δ Υ ], where Υ (0) = 3/128(Mπf ) −5/3 represents the leading term of the phase's post-Newtonian expansion, and f = ω orb /π. We find for the scalar contribution in equal-mass binaries. Such a correction corresponds to a −6P N order contribution, and we have normalized it against possible values for LISA and galactic cores. Even though it is a very low and negative PN correction, its magnitude is so small as to make its detection via GWs seem hopeless [23,24].
Discussion. This work shows how self-gravitating NBSs respond to time-varying, localized matter fluctuations. These are structures that behave classically: they are composed of N ∼ 10 100 10 −22 eV/µ 2 particles; a binary of two supermassive BHs in the last stages of coalescence emits more than 10 60 particles. Our results show unique features of bosonic ultralight structures. For example, they are not easily depleted by binaries. Even a supermassive BH binary close to coalescence would need a Hubble time or more to completely deplete the scalar in a sphere of ten-wavelength radius around the binary. In other words, the perturbative framework is consistent and robust. This study should be extended to eccentric motion, or to self-gravitating vectorial configurations or even other nonlinearly interacting scalars [70]. Our results should also be a useful benchmark for numerical relativity simulations involving boson stars in the extreme mass ratio regime, when and if the field is able to accommodate such challenging setups.

ACKNOWLEDGEMENTS
We are indebted to the Theory Institute at CERN and to Waseda University for warm hospitality while this work was being completed. We thank Ana Sousa Carvalho for advice and support, and Emanuele Berti for useful comments and suggestions.
V. C. acknowledges financial support provided under the European Union's H2020 ERC Consolidator Grant "Matter and strong-field gravity: New frontiers in Einstein's theory" grant agreement no. MaGRaTh-646597. R.V. was supported by the FCT PhD scholarship SFRH/BD/128834/2017. L.A. acknowledges financial support provided by Fundaçao para a Ciência e a Tecnologia Grant number PD/BD/128232/2016 awarded in the framework of the Doctoral Programme IDPASC-Portugal. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904. We thank FCT for financial support through Project No. UIDB/00099/2020. We acknowledge financial support provided by FCT/Portugal through grant PTDC/MAT-APL/30043/2017. The authors would like to acknowledge networking support by the GWverse COST Action CA16104, "Black holes, gravitational waves and fundamental physics."