Dynamical friction in dark matter superfluids: The evolution of black hole binaries

The theory of superfluid dark matter is characterized by self-interacting sub-eV particles that thermalize and condense to form a superfluid core in galaxies. Massive black holes at the center of galaxies, however, modify the dark matter distribution and result in a density enhancement in their vicinity known as dark matter spikes. The presence of these spikes affects the evolution of binary systems by modifying their gravitational wave emission and inducing dynamical friction effects on the orbiting bodies. In this work, we assess the role of dynamical friction for bodies moving through a superfluid core enhanced by a central massive black hole. As a first step, we compute the dynamical friction force experienced by bodies moving in a circular orbit. Then, we estimate the gravitational wave dephasing of the binary, showing that the effect of the superfluid drag force is beyond the reach of space-based experiments like LISA, contrarily to collisionless dark matter, therefore providing an opportunity to distinguish these dark matter models.


Introduction
The standard Λ Cold Dark Matter (ΛCDM) model, wherein dark matter (DM) consists of nonrelativistic, collisionless particles, offers an exquisite fit to various observations across different range of scales, from galaxies to large-scale structures of the Universe.Whether this empirical success persists down to galactic scales, where baryons play an important role in shaping the DM distribution (e.g., [1]), remains hotly debated [2].Alternatively, these potential small-scale shortcomings may indicate that an extension of the model is warranted, for example endowing DM with self-interactions that are strong enough to modify the galactic core [3,4].Another wellmotivated extension is to consider ultralight DM particles, called fuzzy dark matter [5,6], which have a de Broglie wavelength of order kiloparsec and behave as a Bose-Einstein condensate on astrophysical scales [7].
Yet another possibility is that DM consists of self-interacting, sub-eV bosonic particles that can achieve superfluidity.Over a certain mass range and scattering length, these particles can generate a condensate in sufficiently dense environments and low enough temperature [8].In some realizations of the model, long-range interactions between baryons can be mediated by phonons, thereby affecting galactic dynamics of ordinary matter [9][10][11][12][13].The simplest superfluid theory, based on quartic self-interactions, was initially ruled out in Ref. [14] from Bullet Cluster constraints and observations of galactic rotation curves [15,16], while recent weak-lensing observations constrain the model when a Modified Newtonian Dynamics (MOND)-like behaviour on galactic scales is approached thanks to a direct coupling to baryons [17].However, it has been recently revitalized in Refs.[18][19][20], by relaxing the assumptions of global thermal equilibrium and spherical symmetry.
It has been long established that the DM distribution affects the galaxy formation history [21,22].Density profiles obtained from numerical fits or N -body cosmological simulations, such as the Navarro-Frenk-White (NFW) profile [23], suggest that the mass distribution is peaked near the galactic center and dies off at larger distances.The predicted distribution is, however, affected in the presence of a central black hole (BH), depending both on the properties of the DM particles and the formation history of the central object.In an early study, Gondolo and Silk developed a Newtonian scheme to study how collisionless DM redistributes around the BH starting from an initial power-law distribution [24].They showed the formation of a spike, with a density peak near the radius of marginally bound circular orbits, below which the DM particles are accreted into the BH.The analysis has been later improved in Ref. [25], where a fully general relativistic calculation showed that the spikes can get significantly higher and closer to the central BH.On the other hand, the gravitational scattering of stars in the inner region can heat up the DM fluid and give rise to a softened profile [26][27][28][29] or even to disruption [30], while DM annihilations may induce a smoother cusp [31,32].Similarly, the presence of DM self-interactions [33][34][35][36] may wash out any initial and/or adiabatically altered particle distribution and give rise to smoother spikes [34].
Recently, two of us calculated the density profile of superfluid DM around supermassive BHs assuming different superfluid equations of state, for example when it is predominantly characterized by two-body and three-body interactions [37].It was found that, depending on the distance from the central object, the DM density is characterized by different functional behaviors: at large distances the DM appears as a superfluid core, which disappears into a cusp as we approach smaller distances, with an overall enhancement by a factor of order 10 2 − 10 3 relative to the asymptotic core density.The characteristics of this profile distinguish it from the standard predictions for collisionless DM [38].
The formation of these superfluid spikes may have important consequences for the formation of binaries near the galactic cores, where a compact object can become gravitationally bound to the central BH.These binaries have a mass ratio that spans values of order 10 −4 to 10 −3 for intermediate mass-ratio inspirals [39], and 10 −6 to 10 −5 for extreme mass-ratio inspirals [40].They can be observed with upcoming third-generation gravitational wave (GW) experiments like the Einstein Telescope [41][42][43] and Cosmic Explorer [44], or space-based interferometers like LISA [45], respectively.The impact of a DM spike on these binaries is twofold: it modifies the gravitational well felt by the secondary compact object, and it induces dynamical friction drag forces that slow down the binary's trajectory [46].Dynamical friction is a phenomenon wherein an object (perturber) moving through a discrete or continuous medium generates a density fluctuation (wake) in that medium.The gravitational attraction from this wake has the effect of slowing down the motion of the object.The pioneering study of Chandrasekhar [47] considered a perturber linearly moving in a collisionless medium, while Refs.[48][49][50] focused on linear motion in a gaseous medium.The case of circular binaries was also investigated in Refs.[51][52][53][54][55][56].More diverse environments, such as fuzzy dark matter, selfinteracting or superrandiance-driven scalar clouds, and Bose-Einstein condensates, have also been explored recently [6,[57][58][59][60][61][62][63][64].
For the case of superfluid DM, Ref. [65] computed the dynamical friction force experienced by a perturber moving through a uniform core either at constant velocity or accelerated due to an external field.For the steady motion, it highlighted the relevance of higher-gradient corrections in the phonon dispersion relation, which are responsible for the so-called quantum pressure in the hydrodynamical description.Our goal in this work is to understand the impact of a density spike on the evolution of a binary system orbiting within a superfluid environment.To this end, we compute the dynamical friction force experienced by a small perturber moving in circular orbits around a massive black hole surrounded by a DM spike, and then discuss the relevance of the drag force relative to the emission of GWs.
The paper is organized as follows.In Sec. 2 we review the basic notion of the superfluid DM model and the corresponding spikes around BHs.In Sec. 3 we calculate the dynamical friction force for circular binaries.In Sec. 4 we put these results together to discuss the evolution of black hole binaries.We conclude and summarize our results in Sec. 5. Two Appendices are devoted to technical details.In the following, we use natural units ℏ = c = 1.

Excursus on superfluid dark matter
In this Section, we summarize the basic description of the theory of superfluid DM [9][10][11][12][18][19][20].In this model, DM particles are sufficiently light and strongly self-interacting to thermalize and condense at the center of galaxies, creating a superfluid core.

Theory of dark matter superfluidity
The simplest realization of the superfluid DM model consists of a non-relativistic massive complex scalar field ψ, with mass m DM and quartic self-interactions, described by the Hamiltonian: We assume that this bosonic sector is minimally-coupled to gravity.The strength of the contact interactions is parameterized by the coupling g 2 , which is responsible for Bose-Einstein condensation and for generating a stable superfluid state, if it is positive, g 2 > 0. In other words, repulsive selfinteractions are essential for the system to exhibit superfluidity.The theory has a global U(1) symmetry, with the corresponding Noether charge given by the particle number density n.
From an effective field theory perspective, the superfluid condensate is described by a classical field configuration with finite number density, which spontaneously breaks the U(1) global symmetry.At zero temperature, the superfluid equation of state is given by in terms of the mass density ρ = m DM n.The corresponding spectrum of gapless perturbations (i.e., phonons) around the homogeneous solution is given by [65,66] where is the adiabatic sound speed of the superfluid, corresponding to the first derivative of the pressure with respect to density.The non-relativistic approximation inherent in Eq. ( 1) is valid as long as c s ≪ 1.The first term in Eq. ( 3) describes the tachyonic contribution associated with the system's gravitational instability; the second term is the energy cost to excite a sound wave of momentum k; and the third term represents the kinetic energy k 2 2m DM of a constituent of the system.In the following, we ignore the thermal corrections to the equation of state and sound speed, proportional to the temperature of the system [12].
From the above dispersion relation, one can show that modes softer than the Jeans scale k J are unstable [20,66] The corresponding Jeans length scale, λ J = 2π/k J , sets an upper bound on the size of a gravitationally stable superfluid core.One can identify two different regimes depending on the magnitude of the dimensionless parameter ξ, which encodes the nature of the pressure sustaining the core.
• For negligible self-interactions, ξ ≪ 1, the "quantum pressure" of the system prevents the gravitational collapse, giving rise to a Jeans scale of This is the Jeans scale which is characteristic of fuzzy DM models, a regime we are not considering in this paper.
• On the other hand, in the regime ξ ≫ 1 self-interactions are strong enough to be the main contribution in stabilizing against the collapse of the superfluid.In this case, the superfluid core is sustained by the interaction pressure, with the corresponding Jeans scale This is the regime of interest to describe a superfluid state, and we henceforth refer to it as interaction pressure case.
Let us mention that different theories can be considered as well.For example, one can also envisage a theory beyond quartic self-interactions, where the superfluids are mostly described by three-body interactions [10].The corresponding interaction Hamiltonian would involve a hexic potential of the form in terms of the coupling constant g 3 .The corresponding zero-temperature equation of state and sound speed read in terms of the cutoff scale Λ = 1/ 8g 3 m 3 DM .In the non-relativistic limit, the phonon dispersion relation for the three-body case is equivalent to the one for quartic interactions, shown in Eq. ( 3).The interested reader can find additional details in Appendix A.

Formation of the superfluid phase
Self-interacting bosons can generate a superfluid state through Bose-Einstein condensation, which occurs if the conditions of degeneracy and thermalization are satisfied.However, it is important to clarify the role of these two conditions in forming the dark matter superfluid phase.
Degeneracy occurs when the de Broglie volume of the system is highly occupied, and it can be satisfied also by systems which are out-of-equilibrium.Examples are models of fuzzy dark matter or cosmological axions, in which self-interactions are never sufficiently efficient to establish equilibrium.In general, one can define the coherence scale λ c of the system, describing the scale below which the gas of DM particles can be considered as an effective condensate.This is set by the minimum between the scale of gravitational stability and the de Broglie wavelength [67] Here, v is the velocity dispersion of the system.Let us investigate what properties DM should have to determine a coherence length of order kpc.Without thermalization (so without a mechanism that would push all particles in the zero momentum state), the typical velocity dispersion of DM particles in galaxies is v ≃ 10 −3 and extremely light masses (m DM ≃ 10 −21 eV) are required to have coherence over kpc scales.In this case the de Broglie wavelength is comparable to the Jeans length, and a coherent, kpc-size, degenerate soliton is formed.
On the other hand, if one focuses on heavier DM particles, it is the de Broglie wavelength that sets the scale of coherence, with sub-kpc extent.For these masses, in order to obtain a superfluid state which is coherent over the Jeans length and extends to kpc scales, both degeneracy and equilibrium are necessary (the latter condition responsible for decreasing the velocity dispersion of the system).These requirements can be translated into bounds on the parameters of the model.
The condition of degeneracy requires that the de Broglie wavelengths λ dB of the particles overlap in galaxies.By requiring that λ dB is larger than the average interparticle separation, one gets an upper bound on the DM mass [9,10] eV , (11) in terms of the DM halo mass M DM at virialization.This condition therefore implies that only light enough particles form a Bose-Einstein condensate, while heavy particles do not.This inequality can be milder if one considers halos that are only partially degenerate, but we will not be concerned with this possibility.
The second requirement, DM thermalization within galaxies, is achieved when the bosonic particles reach their maximal entropy state via self-interactions.Central regions of the halo, characterized by higher densities with respect to its outskirts, are expected to catalyze thermalization, by enhancing interaction rates and facilitating equilibration.This defines a second important length scale, the thermal radius R T .This is the radius within which dark matter particles experienced at least one interaction over the galaxy lifetime t g .Therefore, the region of the halo which is in thermal equilibrium is defined by the condition Here, σ = m 2 DM g 2 2 /4π is the two-body scattering cross-section, while ρ is the density of the halo, which is assumed to be NFW-like.We also defined the Bose-enhancement factor, N = nλ 3 dB .If the gas is degenerate, N ≫ 1, then interactions are enhanced by the highly occupied phase space. 1 The thermal radius enters in the non-linear equation (12) through the density profile and the velocity.
To summarize, if sub-eV dark matter particles are considered, then the halo is degenerate and a macroscopically large number of particles occupy the ground state and give rise to a Bose-Einstein condensate in the region within R T .Of course, if the thermal radius is larger than the virial radius itself, the whole halo will undergo the superfluid phase transition.We discuss these two scenarios in more detail in the next section.
The self-interaction cross section is also subject to an upper bound from merging galaxy clusters, such as the Bullet Cluster [15,16].In such merging events, the gas component is observed to be displaced with respect to the DM component.This observation requires that the DM particles must experience a negligible amount of scatterings while passing through the target cluster, giving rise to the bound under the assumption of 2-body scattering.Here N is evaluated for typical cluster density distributions.Notice that this is not the "standard" Bullet Cluster bound σ m DM ≲ cm 2 g [15].As discussed in [18,20], the interaction rate is enhanced by the Bose-enhancement factor for degenerate particles, making the bound way more stringent.If we consider particles heavier than the eV scale, 1 If the gas is in equilibrium, the equipartition theorem for ideal gases allows to rewrite the condition N ≫ 1 as where Tc is the critical temperature of the superfluid, and ζ the Riemann zeta function.
the bound relaxes to the well-known Bullet Cluster bound.It is important to emphasize that, as demonstrated in [20], this enhanced bound reverts to the standard non-degenerate counterpart also if the majority of dark matter in the halo is in the superfluid phase (i.e. when thermalization is also achieved).This is because the derivation of ( 14) relies on the assumption that particles are generally scattered into highly occupied energy levels.However, if most of the dark matter is in the condensed phase, this assumption no longer holds, and the use of N becomes invalid.

Superfluids in galaxies
Once the conditions of degeneracy and thermalization are fulfilled, a superfluid phase of size R T is generated at the center of the dark matter halo.As long as the interaction pressure dominates over the degeneracy pressure, we have implying that the superfluid phase is, in general, susceptible to Jeans instability.The net result is that the core of size R T fragments into self-gravitating superfluid configurations.The profile of these lumps can be computed by solving the hydrostatic equilibrium equation under the assumption of spherical symmetry.Using the equation of state in Eq. ( 2), the solution to the coupled Jeans-Poisson equations is in terms of the core central density ρ 0 , and diameter λ J .Notice that in the case of a two-body interacting superfluid, the latter coincides with the Jeans scale in the interaction pressure case (7).
To have an idea of its dependence on the parameters of the model, we can rewrite it in terms of the scattering cross section as Let us stress that, in order to focus on the regime λ J ≳ 2 kpc, we have to consider dark matter halos which are mostly in the superfluid phase.In this way, the enhanced bullet cluster bound (14) relaxes to its non-degenerate analog.If this were not the case, these values of the Jeans scale would be in conflict with the enhanced bound and therefore be excluded. 2 However, the collection of solitons that would reside in the halo is not expected to be stable.For solitons residing outside the central region of the halo, the assumption of self-gravitation breaks 2 For the three-body interacting case, with the equation of state given by Eq. ( 9), the core profile takes the form [10] Notice that, in contrast to the two-body case, the soliton size R3 depends on the central density ρ0.Following Ref. [10], throughout the paper we assume the fiducial values mDM = eV and Λ = meV for the model's parameters, which are possibly in tension with the Bullet cluster bound [20], but show indicatively how the picture changes for different equations of state.This region includes a central superfluid soliton, surrounded by a collection of weakly interacting streams of superfluid debris orbiting around.This debris is the result of the tidal disruption of the population of solitons that formed in the outskirts of the thermal core.Depending on whether the halo is in thermal equilibrium or not, an outskirt of out-of-equilibrium degenerate particles is present outside the thermal core, up to the virial radius of the halo.
down.In [20], it was shown tidal effects are effective in destroying them.The net result of this process is the formation of weakly-interacting superfluid streams, orbiting around a leftover central soliton. 3See Fig. 1 for a pictorial representation.Neglecting baryons, this complex distribution is then expected to match a coarse-grained NFW profile [23,69] due to the weakly-interacting nature of the debris.We refer the reader to the paper [20] for more details on the formation of tidal streams.
Finally, we can quantify an upper bound on the size of the central soliton.In order to see how comfortable we are in the choice of the model parameters to get a superfluid soliton of tens of kpcs in size, it is useful to express the equation of state at matter-radiation equality in terms of the Jeans scale [20] P ρ eq ≃ 10 One could therefore admit sizes of about 30 kpc to get a ratio of pressure and density at matterradiation equality of (P/ρ)| eq ≃ 0.01, which is consistent with CMB observations [70].

Superfluid density profile around black holes
Galaxies are expected to host a supermassive BH at their center.The presence of a supermassive BH modifies the density profile of the central superfluid soliton, due to its gravitational potential.One expects that the density profile steepens in its vicinity, from the otherwise nearly homogeneous core, and that DM self-interactions are correspondingly enhanced.In this Section, we briefly summarize the main results derived in Ref. [37], where the interested reader can find additional details.We consider a non-spinning, central black hole with mass M BH and radius given by Note that we have in mind a supermassive BH which can potentially be observed with LISA.The BH gravitational well is expected to modify the DM distribution within the characteristic radius By definition, r h is the radius at which the DM potential equals the potential induced by the BH or, equivalently, where the total enclosed DM mass and the BH mass are comparable.To simplify the problem, we assume a spherically-symmetric DM distribution, which is also isotropic in velocity and has relaxed to a near-equilibrium state.We further assume that the BH mass is much smaller than the total integrated DM halo mass M DM , while it dominates the total mass of all particles bound to it in the cusp.We describe the density profile in three different regions.In the outer region, beyond the BH sphere of influence r h , the BH has negligible effects on the superfluid DM.The latter is therefore described by a superfluid core and its match to the standard NFW envelope, as discussed in the previous Section.In the intermediate region, at radii smaller than r h , the BH starts to modify the DM density profile because of its gravitational well.It continues until relativistic effects to the DM motion become so important that they modify the superfluid properties and equation of state.At this point, an inner region starts, where the DM profile is further modified.The interested reader can find a more detailed discussion in Ref. [37].
In the following, we will briefly review the DM profile under the influence of the BH.At distances r ≲ r h , the BH starts to modify the density profile due to its potential well.Within this region, the DM particles are gravitationally bounded to the BH, with characteristic velocity Assuming energy equipartition, the DM particles can be described as having a temperature Within this region, we will assume that the BH influence is, however, not strong enough to modify the equation of state of the DM fluid.The DM density profile can then be computed by solving the hydrostatic equilibrium equation The first term of the middle expression comes from the superfluid self-gravity, while the second term from the BH potential Φ BH = −GM BH /r.In the last step, we have made the simplifying assumption that at radii smaller than r h the BH dominates the gravitational potential.At sufficiently small radii, when DM velocities approach the speed of light, the non-relativistic approximation breaks down.Neglecting the small contribution to the DM stress-energy tensor, one can adopt the Schwarzschild metric to describe the spherically-symmetric space-time as such that the relativistic generalization of Eq. ( 24) reads [34] dP (r) The solution to this equation provides the density profile of the superfluid DM in the intermediate region.
From Eq. ( 23), the temperature of the superfluid is connected to the distance from the black hole.Because of this, there exists a specific scale r deg [37] at which the temperature of the system may exceed the critical temperature of the superfluid.In this case, the condition of degeneracy for Bose-Einstein condensation breaks down, and DM is not sufficiently cold to remain in a superfluid state.Thus the DM is no longer degenerate, and its equation of state is instead approximated by the ideal gas law because of the large number of interactions around the BH [34] In the two-body case, it turns out that the r deg is smaller than the BH horizon, such that the DM remains in the superfluid state all the way to the BH horizon.In the three-body case, however, the degeneracy radius is larger than the BH horizon, and the evolution of the DM profile is correspondingly modified.(Note, however, that thermalization never breaks down as we approach the BH.)Following Refs.[34,37], one can compute the DM density profile in this regime by solving the hydrostatic equilibrium equation together with a heat equation within the gravothermal fluid approximation [34,[71][72][73].These equations must be solved down to the radius r mb = 4GM BH [25,34], which corresponds to marginally-bound circular orbits in the Schwarzschild geometry.At r mb = 4GM BH , DM particles are completely accreted into the BH, and the DM density profile plummets.
Figure 2 shows the density profile for the superfluid DM models discussed above, assuming a supermassive BH of mass M BH = 10 6 M ⊙ .The green and blue lines show the profile for the twoand three-body interacting case, respectively.In both cases, the profile becomes increasingly steep within the BH sphere of influence, growing by orders of magnitude compared to the case without the BH, until it drops at the accretion radius for the three-body case.Comparing the two interacting models, one can appreciate that the profile for P ∝ ρ 3 increases slightly less than in the P ∝ ρ 2 case.This is due to the fact that, as we go to smaller radii, the pressure in the three-body case is significantly larger, resulting in an overall milder density growth.
The role of pressure (or self-interactions) crucially impacts the growth of the density profile.This is manifest when one considers the standard case of collisionless dark matter (CDM), as  studied in the pioneering work by Gondolo and Silk [24].From an initially isotropic profile of the form ρ = ρ 0 (r/r h ) −γ , they studied the evolution of a point mass at the center of the distribution, which grows adiabatically by accreting particles and gives rise to a massive BH.The dark matter feels the BH potential and evolves until it reaches the final configuration [24] In Fig. 2 we show the behavior of the cusp for collisionless DM assuming no tilt γ = 0 for the initial profile.As one can appreciate, the absence of self-interactions results into a much larger growth (about fifteen orders of magnitude) with respect to the case of superfluid DM.The relevance of interactions is further confirmed if one considers a model of self-interacting dark matter (SIDM), studied in Ref. [34].Assuming a velocity-independent cross section, it was found that the density profile has an approximate slope of ρ SIDM ∝ r −3/4 .This model is shown in brown in the Fig. 2, highlighting the overall trend that weaker interactions further enhance the profile compared to models with larger pressure, as is the case of superfluids.The huge hierarchy between these different DM models will manifest in the strength of the dynamical friction force experienced by a perturber moving through the fluids.
3 Dynamical friction: circular orbits Dynamical friction is a dissipative phenomenon experienced by a massive object in motion within a medium.It corresponds to the gradual depletion of the body's momentum due to its gravitational interactions with the wake it generates in the surrounding medium.Since the pioneering work of Chandrasekhar [47], the standard setup for studying dynamical friction corresponds to considering a body in linear motion with constant velocity.This setup has been widely studied for different media, such as collisional gases [50], superfluids [65] and fuzzy dark matter [57].
Although the case of linear motion may provide an accurate approximation of dynamical friction for general trajectories, this approximation breaks down when the timescale over which friction significantly alters the dynamics of the massive object becomes comparable to the timescale over which the general trajectory deviates from linear motion.Binaries enter this last category, as the effect of dynamical friction may be important over several orbital cycles.In this Section, we derive the dynamical friction of bodies moving in a circular motion in a homogeneous superfluid medium, following the hydrodynamic set-up employed in Refs.[55,56].We sum up the preliminary steps and refer the reader to these papers for more details.
Consider an external probe of mass µ (which in the next Section will stand for the binary reduced mass) moving in a homogeneous superfluid medium of density ρ 0 .The linear fluid response of the medium to the orbiting probe is determined by solving the Euler and continuity equations for the fluid density and velocity [74], which represent the hydrodynamical form of the nonrelativistic Gross-Pitaevskii equation, where the last term in the second equation identifies the so-called "quantum pressure" term.One perturbs these hydrodynamical equations around a static background of homogeneous density ρ 0 as where the density inhomogeneity α is created by the probe moving in the fluid.This density inhomogeneity, together with the external perturber ρ ext = µ δ (3) ⃗ r−⃗ r p (t) , both source a perturbation ϕ in the gravitational potential: with ⃗ r p (t) denoting the probe position at time t.It is easy to show that the perturbed equations can be recast as a second-order equation for the density inhomogeneity [55,56], On the left-hand side, we have neglected the gravitational mass term coming from self-gravity, −m 2 g α = −4πGρ 0 α [65].This term is relevant for modes with sufficiently low momenta in the dispersion relation.Its effects on dynamical friction have been described in Ref. [65] for a linear trajectory.We briefly discuss its role in Sec.3.3, but leave its complete inclusion to future work.Furthermore, we stress that adopting this expression for the dispersion relation,

DM
, allows us to study also the region of momenta smaller than the chemical potential, discussed in Appendix A, with a redefinition of the effective sound speed and particle mass.Equation ( 32) can be solved in terms of the retarded Green function, through a convolution with the source term [55,56].Dynamical friction is then derived as the gravitational attraction between the wake α(⃗ r, t) and the probe as where ⃗ r ′ p (t) = ⃗ r p (t − τ ) is the position of the probe evaluated at the retarded time t − τ and ϵ > 0 enforces causality.Notice that the formula ( 33) is general and assumes only the homogeneity of the unperturbed medium.
Specializing now to a probe moving on a circular orbit of radius r 0 and constant angular velocity Ω, we parameterize ⃗ r p (t) in spherical coordinates, choosing the motion to lie in the equatorial plane: Using the standard spherical harmonics decomposition, the exponentials can be expanded as Meanwhile, in the polarization basis (ê + , ê− , ẑ), with ê± = 1 √ 2 (iŷ ∓ x), the momentum vector ⃗ k is decomposed as Using these, the friction force can be simplified to Since the circular motion lies in the plane x − y, spherical symmetry implies that the force along the ẑ direction vanishes.Furthermore, since the drag force is a real vector, one can further simplify the integrals in the ŷ direction and integrate over the solid sphere using Wigner 3j symbols, such that the expression reduces to a 1 + 1 dimensional problem.The resulting force is given by [55,56] where we have introduced a dimensionless force ⃗ F DF in the form of a sum over angular multipoles (ℓ, m): The cutoff ℓ max will be discussed below.The unit vectors r and φ point, respectively, in the directions perpendicular and tangential to the circular orbit.Equation ( 39) thus provides the perpendicular and tangential components of the drag force experienced by a probe moving on circular orbits.Notice that the assumption of circular orbit is justified as long as the acceleration induced by the drag force is smaller than the orbital acceleration.We will discuss the consistency of this assumption in the next Section.The above expression involves two functions, γ ℓm and S m ℓ,ℓ−1 , discussed below.
The first function, γ ℓm , is the result of the angular integration of the original integral (33).It is given by in terms of Gamma functions Γ.Notice that γ ℓm is not affected by changing the microscopic properties of the medium, and it is completely determined by the trajectory and the assumed homogeneity of the background.
The second function, S m ℓ,ℓ−1 , is the 1 + 1 dimensional analog of the starting integral in Eq. ( 33).For steady-state motion, it reads Thus the problem of evaluating dynamical friction boils down to evaluating S m ℓ,ℓ−1 .The next subsections are dedicated to solve the integral in different regimes of dark matter interactions and masses.

Full result and discussion
In this Section, we evaluate the function S m ℓ,ℓ−1 by applying the Cauchy integral formula.The procedure is detailed in Appendix B, and we report here the result and main steps of the computation.First of all, we notice that in the procedure of identifying the poles of Eq. ( 41), the following dimensionless quantities emerge.The first quantity is the Mach number which characterizes whether the motion of the perturber through the medium is subsonic (M < 1) or supersonic (M > 1).The second quantity is This describes the value of the azimuthal quantum number m above which the k 4 term in the dispersion relation becomes important.Similar to the case of linear motion, the more the probe moves supersonically, the more important is the quantum pressure in characterizing the friction.
In the process of identifying the poles, ℓ q and M enter through the following two functions With these dimensionless quantities and functions at hand, we may evaluate the integral (41) as described in Appendix B. The result is where i ℓ and k ℓ are the modified spherical Bessel functions of the first and second kind, and θ(m) denotes the Heaviside step-function.In particular, for vanishing azimuthal number, this reduces to The last necessary ingredient to evaluate the drag force in Eq. ( 39) is the maximum multipole ℓ max up to which the sum over ℓ is performed.It is determined by the size R probe of the probe, and radius r 0 of the orbit: Both radial and tangential components converge as ℓ max → ∞.For fixed Mach number, the tangential force starts converging to its asymptotic value only after the multipole moment ℓ q , which carries the information of quantum pressure, is reached.The growth in the tangential component for low multipoles as ℓ max increases is reminiscent of the Coulomb logarithm encountered in the case of a gaseous medium, to be discussed below, which is regulated only when multipoles probing the k 4 /4m 2 DM part of the dispersion relation are included.In other words, when ℓ max > ℓ q , quantum pressure provides the regulator of the theory and the tangential component converges, whereas if ℓ max < ℓ q , the size of the probe provides the effective cutoff.
Figure 3 shows the behavior of the dynamical friction force in terms of the Mach number M, for different choices of the size multipole ℓ max , fixing the orbital radius to r 0 = 3/m DM c s in order to get ℓ q ≃ O(1).In the left panel, one can see that the radial component of the friction force converges after summing up a few multipoles, and that is suppressed for subsonic velocities M < 1.On the other hand, in the right panel, the tangential component of the friction is shown to grow as more multipoles are included.As discussed above, this is reminiscent of the Coulomb logarithm encountered in the sound regime [55] and which is regulated by introducing quantum pressure in the dispersion relation [56].Differently from the sound regime depicted in [55], where the radial friction vanishes in the subsonic regime, neither the radial nor the tangential friction component vanishes.The reason is that quantum pressure provides a non-vanishing subsonic radial component, which is non-negligible in the case of low ℓ q .

Sound regime
In the limit of negligible quantum pressure, the superfluid dispersion relation simplifies to ω k ≃ c s k.This limit recovers the standard sound result and can be obtained if keeping the Mach number M fixed.In this case the momentum integral simplifies to [55] S m ℓ,ℓ−1 = iπ 2c 2 s j ℓ (mM) h (1) with the case m = 0 still determined by Eq. ( 46), albeit without the second term.Let us comment on the limit (48).While this is the correct limit to recover the sound regime, where all multipoles are dominated by sound modes and the dispersion relation is led by the c 2 s k 2 term, in the case of subsonic motion the condition ℓ q ≳ O(1) serves as a reasonable approximation for computing the tangential drag force. 4This approximation works fairly well due to the rapid convergence of dynamical friction in the subsonic limit, once a few multipoles are summed up.Notice that this fast convergence is not there in the supersonic regime, due to the presence of the Coulomb logarithm.Finally, let us mention that in the intermediate case 1 ≪ ℓ q < ℓ max , the setup can be modelled as a sound regime with an effective cutoff given by a size ∼ (m DM c s ) −1 .In this regime, multipoles between ℓ q and ℓ max are quantum-pressure dominated and do not contribute to the Coulomb logarithm.
Figure 4 shows the drag force in the sound regime.In the subsonic limit, the tangential force is always non-vanishing for M ≲ 1 and larger than the radial component.The tangential component is characterized by a large bump at M ≃ 1.4, analogous to the linear-trajectory case discussed in Ref. [50].The radial component presents instead a series of bumps that are due to the overlapping of the Mach cone and the sonic sphere, since the perturber is able to enter its own wake on a circular orbit.These results match the ones previously obtained in Refs.[51,55].The right panels exhibit the convergence of the force when higher multipoles are taken into account, similarly to the full case shown in Fig. 3.
In the subsonic and sound regime, M < 1 and ℓ q → ∞, the tangential component of the force can be summed up in the limit of infinite ℓ max , since the tangential component of the friction is convergent in the limit of subsonic velocities.The resummed version takes the exact form Notice that this resummation is exact up to O (mΩ/m DM c s ) corrections and breaks down only for M ≥ 1.Moreover, if we expand this result for small Mach numbers, the tangential component of the friction reduces to the drag created by a linear trajectory perturber turned on at t = 0 [50], as expected.
In contrast, the radial component does not admit an analytical summed-up version.However, the numerical summation of all multipoles confirms that this term vanishes for small Mach numbers M ≪ 1, as illustrated in the central panel of Fig. 4.

Quantum pressure regime
The second limit we may consider of the full result (45) is the one in which quantum pressure is the main contribution to the dynamics of the wake.This scenario is obtained if ℓ q is small, Once this condition is saturated, the Green function determining S m ℓ,ℓ−1 is mostly dominated by quantum pressure for every m. 5 In this limit, dynamical friction reduces to the case of fuzzy dark matter, with S m ℓ,ℓ−1 determined by Notice that this is the same result as obtained for fuzzy DM in Ref. [56], albeit with a caveat.While the exact fuzzy DM limit is infrared divergent, with the divergence localized in the (ℓ = 1, m = 0) term, the full result (45) is not.This implies that the tail of corrections is important in characterizing the divergence of this specific contribution.
We close this Section with a discussion of IR divergences.In general, there is a natural infrared cutoff for the theory, given by the Jeans scale λ J , which contributes to the dispersion relation of the fluid fluctuations as shown in Eq. ( 3).However, in the case of superfluids, there is a second emergent scale in the problem, determined by the healing length (m DM c s ) −1 .In regularizing the expanded result ( 52), a natural cutoff would be determined by the minimum length between the Jeans scale and the healing length It is worth noting that this cutoff is anticipated to naturally arise in the computation if self-gravity is reintroduced.By a comparison of these scales, it is straightforward to see that the choice of the cutoff boils down to the relative hierarchy of the Jeans mass m g = √ 4πGρ and the chemical potential m DM c 2 s .If the first dominates, then the dynamics of all mode functions within the fluid region are mostly determined by quantum pressure.This is the case in which the full result of Eq. ( 45) reduces to the fuzzy dark matter limit depicted in Ref. [56].
On the other hand, if m DM c 2 s ≫ m g , the infrared cutoff is determined by the healing length.In this case, the dynamics of the fluid fluctuations on the scale k ∼ r −1 0 is dominated by the quantum pressure, but, as the typical momentum of the fluctuation gets softer, there is an intermediate scale λ −1 J ≪ k ≪ r −1 0 at which fluctuations become sound modes.In general, these modes are not important in characterizing the dynamics of the wake, as the modes relevant for dynamical friction are typically on scales of size r 0 .Only the (ℓ = 1, m = 0) term of Eq. ( 45) is sensitive to them due to its infrared divergent behavior.

Summary of the Section
We summarise the result of this Section.The value of the parameter ℓ q allows us to understand which regime of the theory we are probing.
Large ℓ q ≫ 1 indicate the dominance of the sound speed term c s k in the phonon dispersion relation.As ℓ q increases, more modes with wavelengths shorter than the size of the orbit fall within the sound regime.While these modes may not be important in the subsonic regime, where modes comparable to the size of the orbit primarily determine the friction, for supersonic motion also short wavelengths contribute significantly.This is a reminiscence of the Coulomb logarithm that is encountered in the analysis of dynamical friction in the gaseous medium.Therefore, we distinguish the following two sub-cases: • If ℓ q > ℓ max , all modes which are relevant for the problem are in the sound regime.This corresponds to the case where quantum pressure is negligible at all scales and in which the friction is approximated by [55].
• If ℓ q < ℓ max , quantum pressure regulates the Coulomb logarithm at scales bigger than the peculiar size of the perturber.Due to the convergence of dynamical friction, multipoles which are bigger than ℓ q are not important in characterizing the friction.This scenario can be modelled as a sound regime with a perturber of effective size ∼ (m DM c s ) −1 .
Conversely, small values ℓ q ≲ 1 suggest a quantum-pressure dominated regime.In this case, modes with a wavelength larger than the size of the orbit provide an important contribution to the dynamical friction.This is because, in the limit of zero self-interactions, there is an infrared divergence localized in the (ℓ = 1, m = 0) term, which is regulated by the healing length in the case of the superfluid medium.However, if the Jeans scale of the problem is smaller than the healing length, the first provides the regulator and the friction reduces to the result of Ref. [56].Finally, intermediate values indicate regimes where both contributions effectively contribute in the dispersion relation.

Evolution of black hole binaries
In this Section, we can finally discuss the evolution of a black hole binary system surrounded by a superfluid dark matter environment.We focus on extreme mass-ratio inspirals, i.e., with tiny mass ratio m BH /M BH ≪ 1.Such systems can be modelled within BH perturbation theory by studying the quasi-adiabatic orbital motion of a point-particle with mass m BH (the secondary) around a much heavier BH with mass M BH (the primary, taken to be non-rotating).To simplify the setup, we assume that the dark matter generates an enhanced spike only around the heavier BH, as described in the previous Sections, with the orbiting lighter BH moving through the superfluid and subject to dynamical friction.As we will show, the drag force will compete with gravitational radiation to determine the binary's evolution.

Setting the stage
In order to understand the role of dynamical friction on the evolution of binaries of compact objects, we consider an extreme mass-ratio inspiral system of two BHs, with masses M BH = 10 6 M ⊙ and m BH = 10M ⊙ as reference values, moving in the superfluid core within a binary of radius r.The choice for the mass M BH of the central supermassive BH is motivated by the possibility of observing such systems with space-based interferometers like LISA.Choosing smaller values would give rise to intermediate mass-ratio inspirals, for example in the mass range 10 2 M ⊙ ≲ M BH ≲ 10 3 M ⊙ , which can eventually be observed with future third-generation GW experiments, like the Einstein Telescope and Cosmic Explorer.As discussed in the previous Section, in order to evaluate the dynamical friction force experienced by the lighter component of the binary, let us estimate the main parameters determining the drag force.In Fig. 5, we show in the left panel the Mach number for the two and three body interacting models of superfluid DM.As one can appreciate, the evolution of the sound speed at small radii implies that we are always in the subsonic regime, reaching M ≈ 1 only when the BH strongly modifies the DM density profile.On the other hand, before the BH starts to modify the profile around 10 −4 r h , we are well within the subsonic regime M ≪ 1.In the central panel, we plot the radial evolution of the parameter ℓ q , which quantifies the relevance of quantum pressure.It is manifest that it always satisfies ℓ q ≫ 1, such that within the BH sphere of influence the relevant regime for dynamical friction is the sound regime, ℓ q ≫ 1, as discussed in Sec.3.2.
Finally, in order to provide an analytical estimate of the drag force, we can compute the ratio of ℓ q and the finite-size cutoff ℓ max , defined in Eq. (47).The latter determines the maximum ℓ in the sum over multipoles.Their ratio is given by where in the second equality we have substituted the size of the perturber R probe = 2Gm BH and the tangential velocity v = Ωr = GM BH /r.The estimate, supported by observing the radial behavior shown in the right panel of Fig. 5, proves that the ratio is always larger than unity, implying that the size provides the cutoff of the theory and that quantum pressure effects can be safely ignored.
There, the reference values m DM ≃ µeV (two-body case) and m DM ≃ eV (three-body case) have been assumed for the masses of the dark matter candidate.Let us stress that this conclusion holds only for heavy enough DM particles, in particular for the range of masses m DM ≫ 10 −9 eV (twobody interactions) and m DM ≫ 10 −3 eV (three body interactions).For lighter particles, it is not possible to ignore the influence of quantum pressure because, as the orbit of the binary shrinks, it would eventually reach a scale at which quantum pressure becomes the dominant contribution to the dynamics.The analysis of this general case requires a numerical investigation, since one must evolve the orbit using the complete formula for dynamical friction.We postpone this discussion to future work.
To first approximation, as long as we are in the sound regime (ℓ q ≫ 1) and subsonic limit (M ≪ 1), we can focus on the tangential component of the friction force, with resummed expression given by Eq. ( 50).An important caveat is that our derivation of dynamical friction relied on a homogeneous background, while the presence of the central BH with mass M BH induces a nontrivial density profile, which can in principle invalidate the computation.On the other hand, as long as the growth in density is sufficiently mild, in a sense made precise below, then an "adiabatic" approximation is justified for dynamical friction.
To see that this is indeed the case here, first note that the radial drag component for a subsonic perturber in the sound regime is suppressed.Furthermore, along the radial direction and over the minimum length scale relevant to estimate the drag force in the sound regime, set by the characteristic size of the overdensity Gµ/c 2 s [50], the radial growth of the density profile is indeed negligible, Thus we are justified in studying the evolution of an orbit of radius r by replacing the homogeneous density profile ρ 0 in the dynamical friction formula with the non-homogeneous density profile ρ(r), evaluated at a given r.In this way, Eq. ( 50) for the tangential component of the friction becomes where we have expanded for small Mach numbers and kept the leading order term in M ≪ 1.This approximate expression offers a good to the exact force, as shown in the left panel of Fig. 6.
For completeness, we can compare the drag force obtained for superfluids with the one for a collisionless medium, given by [47] using the Gondolo-Silk density profile shown in Fig. 2. As shown in the right panel of Fig. 6, the dynamical friction force for superfluid dark matter is orders of magnitude smaller than for collisionless dark matter.The reason for this large hierarchy is two-fold: i) collisionless DM exhibits a stronger enhancement in density in the vicinity of the supermassive BH, as shown in Fig. 2; ii) the further suppression of the superfluid force by M 3 in the subsonic regime.This result highlights the stronger role of dynamical friction for collisionless DM than for superfluids, which will be confirmed in the next Section.Finally, comparing the superfluid drag forces in the right panel of Fig. 6, one notices that at large distances the force is larger in the three-body case because of its larger Mach number, while at smaller distances the two-body superfluid force dominates because of its further enhanced density profile.Before discussing the role of dynamical friction in the evolution of a binary system, let us comment on the consistency of assuming a circular orbit for the motion of a perturber under the influence of a drag force.This assumption is justified as long as the change of the orbital radius in a period is smaller than the radius itself.By expressing it in terms of the ratio between the acceleration induced by the drag force and the orbital one, one gets because of the large hierarchy in the BH masses and the fact that, within r h , the mass of the DM halo is much smaller than the one of the central BH.

Gravitational wave dephasing
After setting the stage, we can estimate the change in the trajectory of binaries, orbiting for example in the LISA band, due to torques generated by the drag force induced by the presence of the DM environment.This change could eventually be used to detect and measure the effect of DM with GWs, as discussed for example in Refs.[75][76][77].Following the prescription laid out in Refs.[75,78,79], one can estimate the dephasing cycles, which account for changes in the gravitational wave phase induced by dynamical friction within the DM spike.For simplicity, we approximate the orbit of the binary as circular, which is justified by the fact that dynamical friction circularizes inspiralling orbits within a DM spike [80].(See, however, Ref. [81] for a discussion of its effects on eccentricity, whose increase could lead to cusps and O(few) enhancement in the GW waveform.Nevertheless, this increase is not expected to change our main conclusions.) The radial equation of motion for circular binaries in the Newtonian limit is given by [79] where Ω is the source's orbital frequency, and M DM (r) is the mass of the DM spike enclosed within a sphere of radius r.Following the assumptions made in Sec. 2, one can check that the DM halo mass M DM (r) ≪ M BH within r < r h for the regions of parameter space we consider, so that the quantity ϵq(r) ≪ 1.
Because ϵq(r) is small, the orbital evolution can be considered adiabatic, and the r term in Eq. ( 59) can be neglected.Rearranging and expanding to first order in ϵ, one obtains the following expression for r(Ω) [79], where we have defined r 0 = (GM BH /Ω 2 ) 1/3 .The adiabatic evolution of the circular orbit during the inspiral allows us to study the change in the orbital radius by computing how the orbital energy changes as the radius decreases, following the energy balance equation: The orbital energy loss Ėorbit can be expressed using the standard Newtonian expression for a circular orbit as where µ = m BH M BH /(m BH + M BH ) ≃ m BH indicates the reduced mass of the binary and the prime denotes differentiation with respect to r.At the lowest post-Newtonian order, the GW energy loss ĖGW is given by the quadrupole formula, which for circular orbits reads [82] Lastly, the energy loss ĖDF due to dynamical friction is sourced as the inspiraling perturber interacts with the surrounding DM superfluid, which creates a wake behind the massive object that decelerates it.The drag force on the orbiting object is given by Eq. ( 56), such that the energy loss becomes where only the tangential direction contributes to the evolution in the subsonic regime 6 .The relative strength between the gravitational and drag contributions is given, at zeroth order in ϵ, by the ratio in terms of the GW frequency f = Ω/π.This ratio is plotted in the left panel of Fig. 7 for the superfluid two-and three-body interacting cases, spanning a frequency range from the frequency at which the binary semi-major axis is equal to the BH sphere of influence, f r h , to the frequency f ISCO of the innermost stable circular orbit.The latter is computed taking into account self-force effects due to the lighter component of the binary [86], and is approximately given by 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 1 10 -22 7: Left panel: Ratio between the dynamical friction and gravitational radiation contributions as a function of the GW frequency f , for the superfluid models with two-(green) and three-body (blue) interactions, and for a collisionless/CDM medium (red line, obtained using the profile in Eq. ( 28) with γ = 0).The vertical dashed lines indicate, respectively, the frequency fr h at which r = r h , the frequency where the BH profile starts deviating from the core, ρ(r) ≳ ρ0, and the frequency f5yr corresponding to an observation time of five years for LISA.Right panel: GW dephasing for the models discussed above, assuming an extreme mass-ratio inspiral with BH masses MBH = 10 6 M⊙ and mBH = 10M⊙.
Figure 7 shows that, at small GW frequency, the binary evolution is dominated by dynamical friction.In particular, within the superfluid core, GW emission is negligible and dynamical friction provides the main channel through which the binary shrinks.However, once the binary semi-major axis becomes comparable to the radius at which the dark matter profile develops the spike due to the black hole influence, which takes place when ρ(r) ≳ ρ 0 at around r ≃ 10 −4 r h as visualized in Fig. 2, then GW emission starts to dominate the binary's inspiral.For larger frequencies, GW emission becomes the predominant driving channel leading to the binary merger. 7o ascertain the role of dynamical friction in the evolution of binaries observable with LISA, which is characterized by an observation time of at most five years, we have indicated by a dashed line in Fig. 7 the minimal frequency f 5yr for which binary inspirals lie within this observation time.As one can appreciate, the relative strength between gravitational and drag contributions satisfies ĖDF / ĖGW ≪ 1 in the frequency range (f 5yr , f ISCO ), showing that dynamical friction for superfluid dark matter is not sizeable enough to leave an imprint in the binary's evolution.
Using Eq. ( 60), we can invert Eq. ( 61) to get the frequency evolution equation in terms of the functions where one can easily recognize the contributions from GW emission and dynamical friction by checking the sound-speed dependence of a specific contribution.From Eq. ( 67), one can find the orbital phase in the stationary phase approximation, from which the number of GW cycles N cycles (f ) can be obtained as One can then split the contributions from dynamical friction and GW evolution by computing the GW dephasing whose behavior is shown in the right panel of Fig. 7. Values smaller than unity, ∆N cycles ≲ 1, indicate an efficient role of dynamical friction in driving the binary's inspiral.Similarly to what was discussed for the energy ratio ĖDF / ĖGW , one can appreciate that the GW dephasing is manifest only for wide enough binaries, where dynamical friction represents the main force dictating the binary's evolution.This conclusion is robust for both models of superfluid dark matter assumed in the paper, irrespective of their interactions.Moreover, in agreement with Fig. 7, in the range of frequencies that can be probed by LISA, the role of dynamical friction on the binaries is not large enough to be detected.The results obtained for the GW dephasing are based on a series of assumptions in the derivation of Eq. ( 56) for the drag force.Let us briefly comment on the consistency of our approach.First, our derivation assumes subsonic motion for the lighter component in the binary.This is justified at distances r ≳ 10 −4 r h (see Fig. 5), where dynamical friction dominates the evolution, but it breaks down at smaller distances, where the Mach number reaches values O(1).Even though the force may receive corrections, we do not expect our conclusions to dramatically change, both for the good agreement of the analytical expression as shown in the left panel of Fig. 6, and because at smaller distances GW emission already provides the dominant channel in driving the binary inspiral.
Second, our derivation is based on the phonon dispersion relation discussed in Appendix A and shown in Eq. ( 3), which assumes nonrelativistic DM particles and small sound speed c s ≪ 1.These assumptions, however, break down at small radii, where the DM becomes relativistic (see Ref. [88] for a detailed derivation of the phonon dispersion relation for a relativistic Bose gas and Appendix A for few details).This occurs approximately at the orbital frequency While a dedicated computation is needed to have a clear picture of the relativistic case, which may eventually change the details of the drag force expression, we do not expect our main conclusions to be affected, since in the regime f ≳ f rel ≃ 10 −2 f ISCO dynamical friction is already negligible compared to GW emission.Finally, following the comparison and discussion of Fig. 6, one can contrast these results with the ones for the case of collisionless dark matter (also discussed in Ref. [79], where relativistic effects are taken into account).By looking at the red lines in Fig. 7, one can appreciate that dynamical friction for a collisionless rather than superfluid medium has a stronger role in driving the binary inspiral, also in the frequency range probed by LISA.The reason behind this difference lies in the enhancement of the density profile for the two models, see Fig. 2. In particular, as already discussed above, the different nature of the interactions in the superfluid model compared to a collisionless medium is responsible for the suppressed growth of the profile as the BH horizon is approached.The large difference in the enhancement translates into a much weaker dynamical friction force for superfluids rather than collisionless dark matter.

Conclusions
The theory of superfluid dark matter, based on the condensation and thermalization of selfinteracting, sub-eV bosonic particles at the center of galaxies, provides a novel way to conciliate the triumph of the ΛCDM model on cosmological scales with the dynamics on galactic scales.Within this theory, the dark matter bosons are able to generate a superfluid core in the galactic center, surrounded by particles in the normal phase in the outskirts of the halo.The presence of massive black holes in these environments is however expected to modify the superfluid core, creating overdensity regions known as dark matter spikes.
In this work, we investigated how superfluid dark matter spikes affect the evolution of binary systems orbiting in their environment by quantifying the role that dynamical friction plays on the system's inspiral.In particular, we computed the drag force experienced by a perturber assembled in an extreme mass ratio inspiral with the central supermassive black hole and compared its contribution to the binary evolution with the emission of gravitational waves.We found that the amount of gravitational wave dephasing for a superfluid environment is not large enough to be detected with space-based experiments like LISA.This result can be used to disentangle the superfluid theory from the standard predictions of collisionless dark matter.
Our work provides a preliminary step to a full modelling of the evolution of a binary system within a dark matter superfluid and can be extended by including several additional effects.First of all, a dark matter spike induces a change in the mass enclosed within the orbit of the secondary body as the system inspirals, modifying the Keplerian frequency at a given orbital radius from that of a black hole in vacuum [78].The second (relativistic) effect is related to the amount of dark matter particles accreted by the lighter black hole as the binary evolves, which changes its mass and therefore impacts on the system's inspiral [61,62,89].The efficiencies of accretion and dynamical friction are finally impacted by a third effect, that is dark matter feedback, based on the fact that the dark matter distribution should be jointly evolved to determine consistently its effect on the emitted gravitational waves [76].For dynamical friction this amounts to determining the local density of dark matter particles around the secondary [76,77], while for accretion it consists of taking into account the loss of particles from the distribution function which are accreted by the secondary [61,62,90].Including these effects would represent a huge step towards a full characterisation of the evolution of a binary system in superfluids.
Furthermore, we stress that one could generalize our results to other gravitational wave experiments like Einstein Telescope and Cosmic Explorer, which would be sensitive to lighter central black holes, including as well the information of black hole spins, and study different regimes of interactions for superfluid dark matter, for example described by alternative fluid equations of state.
Finally, the suppression of dynamical friction for superfluids compared to cold dark matter may offer a natural explanation to the puzzle concerning the absence of mergers for the five globular clusters orbiting Fornax [91], which should have long ago spiraled to the center if the dwarf galaxy is composed of a cuspy cold dark matter profile.Superfluid dark matter, characterized by longer timescales for dynamical friction, may offer a solution to understand the survival of these clusters [65] (see Refs. [6,57,59,92] for alternative models).Moreover, superfluid dark matter could be accompanied by the formation of vortices [10], which could be stable when a central perturber (like a SMBH) is present [93,94].Their stability could lead to the transfer of angular momentum from the binary system to the fluid, therefore affecting its orbital evolution, and modify the dark matter dynamics around the compact objects.We leave these further refinements and prospects to future work.The first step consists of extending the momentum integration range down to minus infinity by exploiting the parity of the integrand.Second, we utilize the relation j ℓ (x) = 1 2 (h ℓ (x) + h (2) ℓ (x)), where h (1,2) are the spherical Hankel functions.In this way, we can exploit the exponential behavior where the indices i and j run over 1 and 2, identifying the first and second-type spherical Hankel functions.The bar on the indices means to exchange a given Hankel function with the orthogonal one.Explicitly, h ℓ and h ( 2) ℓ , following the same notation of Ref. [56].Contours in the complex plane are chosen depending on the value of i and j: if {i, j} = {1, 1}, we close the contour in the upper half of the complex plane; whereas if {i, j} = {1, 2}, {2, 1}, and {2, 2}, we close the contour in the lower half plane.To simplify the computation, we push the k 0 divergence on the positive part of the complex plane, so that we have to include it only for {i, j} = {1, 1}.In the following, we provide the solution of the integral in the regimes m > 0 and m < 0.

Figure 1 :
Figure 1: Pictorial representation of the central superfluid region of size RT.This region includes a central superfluid soliton, surrounded by a collection of weakly interacting streams of superfluid debris orbiting around.This debris is the result of the tidal disruption of the population of solitons that formed in the outskirts of the thermal core.Depending on whether the halo is in thermal equilibrium or not, an outskirt of out-of-equilibrium degenerate particles is present outside the thermal core, up to the virial radius of the halo.

Figure 2 :
Figure 2: Comparison between the dark matter profile enhanced around a central BH with mass MBH = 10 6 M⊙ for the cases of collisionless (red), self-interacting (brown) and superfluid (green and blue) models.

Figure 3 :
Figure 3: Radial (left panel) and azimuthal (right panel) components of the dynamical friction force in terms of the Mach number M. The orbital radius is fixed to r0 = 3/mDMcs in order to get an intermediate case (ℓq ≃ O(1)) between the quantum pressure and sound regimes.The different lines correspond to different values of the maximum multipole ℓmax in the sum.

Figure 4 :
Figure 4: Left panel: drag force in the sound speed regime in terms of the Mach number M, assuming ℓmax = 10 [51] and r0 = 100/mDMcs to enforce ℓq ≫ 1. Central and right panels: radial and azimuthal components of the force assuming the orbital radius r0 = 30/mDMcs to describe the sound regime.The different lines correspond to different values of the maximum multipole ℓmax in the sum.

Figure 5 :
Figure 5: Left panel: The Mach number M = v/cs as a function of the radial distance, for the two-and three-body interacting models of superfluid dark matter.Central panel: The parameter ℓq, which characterizes the relevance of quantum pressure, for the same models.Since ℓq ≫ 1, we are deep in the sound regime.Right panel: Comparison between ℓq and the finite-size cutoff ℓmax.The dashed lines indicate the accretion radius r mb and the degeneracy radius r deg for the three-body case.We have assumed the reference values mDM = µeV and mDM = eV for the DM mass, respectively, and ρ0 = 10 −25 g/cm 3 for the homogeneous core density.

Figure 6 :
Figure 6: Left panel: Comparison between the numerical value of the azimuthal component of the dynamical friction force (blue line) and its analytical behavior found by resumming the multipole series in Eq. (50) (green line for the resummed analytical expression, and red line for its Taylor expansion).Right panel: Ratio between the superfluid and collisionless friction forces as a function of the distance.
Appendix B Evaluating the S m ℓ,ℓ−1 functionIn this Appendix we calculate the function S m ℓ,ℓ−1 , introduced in Eq. (41), by evaluating its defining momentum integral

2 )
(x) ∼ e −ix (B.2)and apply the Cauchy integral formula.One can identify the poles in the integrand of Eq. (41) ask 0 = 0 ; k 1,2 = ±im DM c s f + m ; k 3,4 = ±m DM c s f − m , (B.3)where the functions f ± m are defined in Eq. (44).The corresponding residues are given by Res