Issue 
A&A
Volume 647, March 2021



Article Number  A70  
Number of page(s)  15  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202039865  
Published online  11 March 2021 
Dynamical friction in BoseEinstein condensed selfinteracting dark matter at finite temperatures, and the Fornax dwarf spheroidal
Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern 0315, Oslo, Norway
email: stian.hartman@gmail.com
Received:
6
November
2020
Accepted:
28
December
2020
Aims. The aim of the present work is to better understand the gravitational drag forces, also referred to as dynamical friction, acting on massive objects moving through a selfinteracting BoseEinstein condensate, also known as a superfluid, at finite temperatures. This is relevant for models of dark matter consisting of light scalar particles with weak selfinteractions that require nonzero temperatures, or that have been heated inside galaxies.
Methods. We derived expressions for dynamical friction using linear perturbation theory, and compared these to numerical simulations in which nonlinear effects are included. After testing the linear result, it was applied to the Fornax dwarf spheroidal galaxy, and two of its gravitationally bound globular clusters. Dwarf spheroidals are wellsuited for indirectly probing properties of dark matter, and so by estimating the rate at which these globular clusters are expected to sink into their host halo due to dynamical friction, we inferred limits on the superfluid dark matter parameter space.
Results. The dynamical friction in a finitetemperature superfluid is found to behave very similarly to the zerotemperature limit, even when the thermal contributions are large. However, when a critical velocity for the superfluid flow is included, the friction force can transition from the zerotemperature value to the value in a conventional thermal fluid. Increasing the mass of the perturbing object induces a similar transition to when lowering the critical velocity. When applied to two of Fornax’s globular clusters, we find that the parameter space preferred in the literature for a zerotemperature superfluid yields decay times that are in agreement with observations. However, the present work suggests that increasing the temperature, which is expected to change the preferred parameter space, may lead to very small decay times, and therefore pose a problem for finitetemperature superfluid models of dark matter.
Key words: dark matter / cosmology: theory
© ESO 2021
1. Introduction
When a massive object moves through a background medium, its gravitational field can cause the background to form an overdensity that trails it, and in turn exerts a gravitational force on the object that produced it. This is known as dynamical friction, and is a purely gravitational phenomenon. It can therefore also arise in systems in which the constituent components otherwise have very weak or no coupling to one another, or behave as collisionless particles, such as dark matter (DM) and stars. Many important processes in the formation of structure, the evolution of galaxies, and the dynamics of astrophysical systems, such as mergers (Jiang et al. 2008; BoylanKolchin et al. 2008), the sinking of satellites into their host halos (Colpi et al. 1999; Cowsik et al. 2009; Cole et al. 2012; Tamfal et al. 2020), the decay of orbiting black holes and binaries (Just et al. 2011; Pani 2015; Dosopoulou & Antonini 2017; Gómez & Rueda 2017), and bar–halo interactions in disk galaxies (Weinberg 1985; Debattista & Sellwood 2000; Sellwood 2014), therefore depend on the nature of this drag force.
The first detailed calculation of dynamical friction was carried out by Chandrasekhar (1943) in the context of stellar dynamics. Chandrasekhar considered the varying gravitational forces acting on a star as it moves through its stellar neighborhood, and found that it experiences a net average force opposite to its direction of motion, that is, a gravitational drag force. He treated the background of stars as an infinite homogeneous gas of collisionless particles following a MaxwellBoltzmann velocity distribution, an approach that can also be used for collisionless DM (Mulder 1983; Colpi et al. 1999; Binney & Tremaine 2008). However, for a collisional medium, pressure forces must be taken into account when computing the dynamical friction, and this has been done both analytically and numerically for various types of gases, such as ideal (Ostriker 1999; SánchezSalcedo & Brandenburg 1999; Lee & Stahler 2011, 2014; Thun et al. 2016), relativistic (Barausse 2007; Katz et al. 2019), and magnetized gases (SánchezSalcedo 2012; Shadmehri & Khajenabi 2012).
The nature of the dynamical friction due to DM is related to the nature of DM itself. The standard model of the universe, ΛCDM, includes cold and collisionless DM as the predominant matter component, making up about 80% of all matter. While extremely successful at explaining observables such as the microwave background radiation, largescale structure, the expansion history of the Universe, and important properties of galaxies (Davis et al. 1985; Percival et al. 2001; Tegmark et al. 2004; TrujilloGomez et al. 2011; Vogelsberger et al. 2014; Planck Collaboration XIII 2016; Riess et al. 2016), the identity of DM has remained elusive. Furthermore, there are discrepancies between simulations of structure formation at small scales, and observations (for reviews, see e.g., Weinberg et al. 2015; Del Popolo & Le Delliou 2017; Bullock & BoylanKolchin 2017). These discrepancies may have their solution within ΛCDM by including more realistic models of baryonic physics in simulations (SantosSantos et al. 2015; Sales et al. 2016; Zhu et al. 2016; Sawala et al. 2016), but the solution may also be in alternative models of DM (Hu et al. 2000; Spergel & Steinhardt 2000; Shao et al. 2013; Lovell et al. 2014; Schive et al. 2014; Elbert et al. 2015; Berezhiani & Khoury 2015; Khoury 2016; Schwabe et al. 2016; Mocz et al. 2017; Tulin & Yu 2018; Clesse & GarcíaBellido 2018; Boldrini et al. 2020). For these reasons, studies have also been carried out on dynamical friction in various DM models, such as fuzzy DM (Hui et al. 2017; BarOr et al. 2019; Lancaster et al. 2020), and selfinteracting BoseEinstein condensed (SIBEC) DM (Berezhiani et al. 2019), also known as superfluid DM. A number of studies have considered finitetemperature effects of interacting superfluid DM (Harko & Mocanu 2012; Slepian & Goodman 2012; Harko et al. 2015; Sharma et al. 2019). Of particular note is the one presented by Berezhiani & Khoury (2015), who suggested that superfluid DM, when provided with a special Lagrangian structure and coupling to the visible sector, can give rise to modified Newtonian dynamics (MOND; Milgrom 1983a,b,c; Famaey & McGaugh 2012) between baryons at galactic scales. This MONDian force is mediated by superfluid phonons, which cease to be coherent on scales larger than galaxies, resulting in the vanishing of the extra force and the preservation of the largescale success of CDM. For the fifth force to be MONDian, the DM particles need exotic threebody selfinteractions, and the DM fluid has to be above a certain temperature to be wellbehaved. Finitetemperature DM might arise through processes inside galaxies that transfer energy to the DM halo (Goerdt et al. 2010; Pontzen & Governato 2012; Read et al. 2019), possibly heating up the DM fluid.
Because the form of the dynamical friction experienced by visible matter embedded in DM halos depends on the properties of DM, observations of galaxies can be used to constrain DM. Dwarf spheroidal galaxies (dSph) are particularly wellsuited for this purpose. Being poor in visible matter, their dynamical behavior is dominated by their DM component and they therefore provide a testing ground for DM models (Battaglia et al. 2013; Walker 2013; Strigari 2018). One such system is the Fornax dSph and its five gravitationally bound globular clusters (GCs; Mackey & Gilmore 2003), with a sixth one recently found to likely be a genuine, albeit dim, GC (Wang et al. 2019a). The orbital decay times of these GCs, in particular the inner two (not counting the recently discovered sixth GC), due to dynamical friction from a CDM background, have been estimated to τ_{DF} ≲ 1 Gyr (Oh et al. 2000; Cole et al. 2012; Hui et al. 2017; ArcaSedda & CapuzzoDolcetta 2017), much shorter than the supposed age of the host system, τ_{age} ∼ 10 Gyr (del Pino et al. 2013; Wang et al. 2019b). Furthermore, there is no bright stellar nucleus at the center of Fornax dSph that would suggest the sinking of other GCs in the past. This apparent mismatch between theoretical prediction and observation suggests one of two scenarios; that we are witnessing these GCs just as they are about to fall into their host, implying a finetuning of their initial positions, which seems unlikely; or that there is some mechanism, or property of DM that stops the GCs from migrating towards the center of the Fornax dSph. This discrepancy between CDM estimates and observations is the socalled timingproblem, and a number of solutions have been proposed, such as massive black holes heating the system (Oh et al. 2000); assuming the CDM profile of Fornax to be cored instead of cuspy (Goerdt et al. 2006; Cole et al. 2012); inaccurate modeling of the Fornax system and the rate of the orbital decay (Cowsik et al. 2009; Kaur & Sridhar 2018; Boldrini et al. 2019; Leung et al. 2020; Meadows et al. 2020); or some exotic property of DM (Hui et al. 2017; Lancaster et al. 2020).
In this work we extend the analysis of dynamical friction in a zerotemperature superfluid to finite temperatures, where the fluid is in a mixed state of normal fluid–made up of thermal excitations–and superfluid. This type of system has pressure terms coming from both thermal excitations and selfinteractions, and can exhibit unique features due to the separate flow of the superfluid and normal fluid components. With an expression for the dynamical friction in SIBECDM, we estimate the time it takes for two of the GCs orbiting the Fornax dSph to sink into their host halo due to this gravitational drag, thereby inferring constraints on finitetemperature superfluid DM. The paper is organized as follows: in Sect. 2 the superfluid equations at both zero and finite temperatures are introduced, as well as some basic notions related to superfluidity. In Sects. 3 and 4 these equations are used to derive the dynamical friction at linear order, both in a steadystate and a finitetime scenario. The dynamical friction is also found using numerical simulations of the full superfluid hydrodynamic equations in Sect. 5, which is compared to the linear result in Sect. 6. The tools developed in the preceding sections are used in Sect. 7 to estimate the decay times of two of the GCs orbiting the Fornax dSph, and constraints on SIBECDM are inferred. In Sect. 8 a summary of this work and the main results are presented. Natural units are used throughout.
2. Hydrodynamics of finitetemperature superfluids
In the standard treatment, superfluids are often related to BoseEinstein condensates (BEC), which form when the temperature is sufficiently low and the particle density high enough that the de Broglie wavelengths of identical bosons overlap, creating a coherent state that can be described by a singleparticle wavefunction. This wavefunction is usually associated with the superfluid, and can therefore be regarded as a quantum mechanical effect at macroscopic scales. The wavefunction ψ at the meanfield level is governed by the GrossPitaevskii equation, a nonlinear Schrödinger equation with effective contact interactions parameterized by g;
The external potential V_{ext} can be a trapping potential, as is often used in cold atomic experiments, or a gravitational potential. The amplitude of ψ is related to the particle number density by n = ψ^{2}, and mass density ρ = mψ^{2}.
By inserting for the wavefunction
and defining the velocity field v = ∇S/m, the nonlinear Schrödinger equation can be reformulated in a hydrodynamic form. The real and imaginary parts of the Schrödinger equation give the set of equations
These are the socalled Madelung equations (Madelung 1926). The first is a continuity equation for mass, and the second is a quantum variant of the momentum equation, with the quantum potential
coming from the kinetic part of the Schrödinger equation that is present even in the absence of interactions. From the definition of the velocity field, we see that it is irrotational, because the curl of a gradient is zero. However, there can arise defects in the superfluid, around which the circulation is quantized as
because the complex wavefunction must be singlevalued. These special structures in superfluids are called quantum vorticies. Both the Schrödinger and Madelung formulations have been used in cosmology as models for DM in order to explain the absence of smallscale structure that is predicted in Nbody simulations of CDM (Schive et al. 2014; Mocz et al. 2017; Nori & Baldi 2018, 2021; Mina et al. 2020a,b).
At finite temperatures, the hydrodynamic formulation of a superfluid must take into account that the fluid is no longer completely superfluid. There is a thermal cloud of excitations in addition to the coherent superfluid state that carries entropy, gives a thermal contribution to the fluid pressure, and can be viscous and rotational. To complicate matters further, as the temperature of the fluid changes, the fraction of the fluid in this thermal cloud changes as well. This property of superfluids, to behave both as a superfluid (in the sense that we usually understand the term, as a fluid with zero viscosity, quantized circulation, and carrying no entropy) and a conventional fluid, has led to the development of a twofluid picture of superfluids. The hydrodynamic equations for a finitetemperature superfluid are (neglecting the quantum potential; Taylor & Griffin 2005; Chapman et al. 2014):
The thermal cloud, which we refer to as the “normal fluid”, has density ρ_{n}, velocity u_{n}, and transports both mass and thermal energy. The second component is the “superfluid”, with density ρ_{s}, a velocity field u_{s}, and carries no entropy. The total mass density is the sum of the two components, ρ = ρ_{n} + ρ_{s}, and likewise for momentum, j = ρ_{n}u_{n} + ρ_{s}u_{s}. The fluid pressure is P, the entropy density S, temperature T, and .
As previously mentioned, superfluids and BECs are related phenomena, but it is important to stress that they are not equivalent. The formation of a BEC does not automatically imply a superfluid. To see this we must consider the cocalled Landau criterion. Landau, in his seminal paper on superfluid liquid helium 4 (Landau 1941), made the following argument: Assume that dissipation and heating in a fluid takes place via the creation of elementary excitations. If these excitations become energetically unfavorable and cannot spontaneously appear, then dissipation and heating ceases, and the fluid becomes superfluid. The criterion for such a condition is for the relative velocity v between the superfluid and a scattering potential, such as an impurity or a container wall, to be smaller than a critical value,
where ϵ(p) is the energy of an elementary excitation with momentum p (Pitaevskii & Stringari 2016). This criterion shows that an ideal BEC, for which the excitation spectrum is ϵ(p) = p^{2}/2m, has v_{c} = 0 and is therefore not a superfluid. On the other hand, a Bose gas with weak interactions has–upon the formation of a BEC–an energy spectrum that is linear at small momentum, ϵ(p) = c_{s}p. Hence v_{c} = c_{s}, and weakly interacting BECs are superfluids.
The Landau criterion is usually derived with the velocity relative to an external scatterer in mind, but it also applies to the thermal excitations that make up the normal fluid. The critical value for the relative velocity w = u_{s} − u_{n} of the normal fluid and superfluid is smaller than the one determined by Eq. (11), but the difference is small at low temperatures and weak selfinteractions (Navez & Graham 2006).
The presence of the relative velocity w, because of the partially independent motion of the superfluid and normal fluid components in a finitetemperature superfluid, has important consequences for its behavior. The superfluid part does not carry heat, while the normal fluid does, allowing mass and entropy to flow separately. This becomes clear if we define the velocity field for the mass flux, v = j/ρ, and express the equations for mass and entropy conservation in terms of w and v;
For a finite superfluid fraction, the entropy has an additional flux term, and therefore entropy and mass can have different flow patterns. This property is called thermal counterflow. The equation for ∂w/∂t contains a driving term S∇T/ρ_{n}, and so the counterflow w tends to be directed towards regions of higher temperature, washing out thermal differences in the superfluid. As we see below, it is this property that makes the dynamical friction in a superfluid different from a corresponding fully normal fluid (i.e., a conventional fluid, ρ_{s} = 0, with the same pressure forces).
When the Landau criterion is broken, with w approaching and passing the critical velocity, the superfluid flow starts to decay as a tangle of quantum vortices form, and causes a mutual friction between the superfluid and normal fluid components (Skrbek 2011; Skrbek & Sreenivasan 2012; Barenghi et al. 2014). Such a dissipative effect is not present in the superfluid equations, but can be included with additional terms, as has been done in numerical studies of superfluid helium (Doi et al. 2008; Darve et al. 2012; Soulaine et al. 2017). However, to circumvent the need for extra parameters and the need to assume the functional form of the mutual friction, we instead follow the same approach used in a previous work (Hartman et al. 2020); the dissipative processes are assumed to take place instantaneously when the relative velocity w exceeds the critical velocity. The velocity field v_{s} is changed in such a way that the fluid momentum is conserved, and that only the magnitude of w is altered, not its direction, bringing it to w = v_{c}. In other words, we assume the mutual friction to be directed along w.
The critical temperature T_{c} is a central quantity in BEC superfluids. For T > T_{c}, a gas of identical bosons is a normal fluid, but for T < T_{c}, the particles begin accumulating in the ground state, forming a BEC, which in turn can form a superfluid. In the threedimensional, homogeneous, ideal Bose gas, this critical temperature is
where ζ(x) is the Riemann Zetafunction, and holds approximately for weakly interacting gases as well (Andersen 2004; Sharma et al. 2019).
For the thermodynamic quantities of a weakly interacting Bose gas, we again follow the approach used in a previous work (Hartman et al. 2020). The equation of state is approximated by an ideal gas with contributions from twobody interactions,
valid only for T < T_{c}. The fraction of particles in the condensate f_{0} and the superfluid f_{s} = ρ_{s}/ρ are both taken to be equal to the condensate fraction in the ideal case;
The critical velocity is approximated as
As long as the temperature is not too close to the transition point, and the interactions are sufficiently weak, these approximations work well.
3. Dynamical friction from steadystate linear perturbation theory
The starting point for computing the dynamical friction acting on an object, or a “perturber”, moving through the superfluid are Eqs. (7)–(10). The gravitational potential is sourced by both the background mass density ρ, and the mass distribution ρ_{pert} of the perturber:
The superfluid is assumed to be homogeneous, and so the fluid variables are expanded to linear order, ρ = ρ_{0} + δρ, S = S_{0} + δS, u_{s} = δu_{s}, and so on. The linear equations are
These can be combined into two coupled equations for δρ and δS;
The “0” subscript indicates that the quantity is evaluated at the background level. As expected, there are scaledependent pressure terms that inhibit the growth of mass density and entropy perturbations, but in the entropy equation there are additional effective pressure terms that further reduce entropy perturbations. These are due to thermal counterflow and depend on the superfluid fraction, vanishing in the fully normal fluid limit. It must be noted that the critical velocity v_{c} is not included in the present approach, but the effect of this on linear theory is considered further in Sect. 4, as well as in Sect. 5 using numerical simulations.
Writing δρ = αρ_{0}, and Fourier transforming into momentum (k) and frequency (k_{0}) space, the solutions of the kmodes α_{k} are found:
where the dispersion relation is
and
The dynamical friction is given by the change in the energy of the perturber,
where M and V are the mass and velocity of the perturber, and Φ_{α} is the gravitational potential of the background fluid,
This is readily found in kspace,
which can be Fourier transformed back into positionspace to give the dynamical friction,
Approximating the perturber as a point particle moving along the zaxis with constant velocity V,
or in kspace
yields the expression for the dynamical friction as
Equation (41) can be tackled by extending the k_{0}integral into the complex place and closing it in the upper half plane (assuming t > 0), meaning that contour integration can be used. The poles are pushed slightly off the real line by the prescription ω_{k±} → ω_{k±} + iϵ, and only the residual of the poles inside the contour contribute to the integral. Taking the limit ϵ → 0^{+} after integrating gives the dynamical friction as
Spherical polar coordinates are adopted for the integral over momentum, with the polar angle θ defined as the angle relative to the direction of propagation, the zaxis, and the force is evaluated at the position of the perturber, . The integrand is independent of the azimuthal angle, but depends on the polar angle through k_{z} = kcos θ. Integrating over the azimuthal angle therefore gives a factor 2π, while the polar angle in combination with the δfunction fixes the exponentials to one and places upper limits on the momentum, , where satisfies kV = ω_{k±}. Further constraints are placed on k: The remaining kintegral is bounded by the finite sizes of the perturber and the cloud it moves through, R_{max} = R_{cloud} and R_{min} = R_{pert}, otherwise both ultraviolet (UV) and infrared (IR) divergences may be encountered, because the perturber is modeled as a point particle, and the background fluid as infinite and uniform. We must also have , where is the minimum momentum for which ω_{k±} are real. At small k, or, equivalently, large scales, where ω_{k±} become complex or imaginary, the background cloud will be gravitationally unstable and deform. We denote as a general measure the upper limits in k for the two terms in Eq. (42) by , and the lower limits by . Inserting the expression for ω_{k±} and using that C_{4} − A = C_{1}, the dynamical friction becomes
There is an implicit criterion that , otherwise the integral is zero.
Equation (43) can be solved analytically, but its expression is not particularly enlightening. Instead, we focus on a few limiting cases for which the force reduces to a simplified form; zero temperature, the fully normal fluid, small velocities, and no selfgravitation.
3.1. Zerotemperature limit
Taking the limit T → 0 (under the assumption that terms such as S^{2}ρ_{s}/ρ_{n} go to zero as well) yields one band for the dispersion relation,
where the sound speed at zero temperature is
The dynamical friction becomes
with
3.2. Normal fluid limit
Taking the fully normal fluid limit ρ_{s} → 0 also gives one band for the dispersion relation,
with the sound speed in the fully normal fluid
The dynamical friction is again given by Eq. (46), but with
This is the same as the zerotemperature case, but with a different sound speed.
3.3. Smallvelocity limit
At sufficiently small velocities, , assuming that the finite sizes of the background cloud and perturber do not set the integral limits in Eq. (43), the dynamical friction becomes
This is equal to the friction force at T = 0 in the same limit, as opposed to when ρ_{s} = 0;
The dynamical friction of a superfluid therefore approaches the zerotemperature limit even when there is a significant thermal contribution. This happens because counterflow in the superfluid conspires against thermal perturbations, allowing the mass overdensity to behave similarly to a zerotemperature fluid. With only the interaction pressure that is present at zero temperature effectively damping density perturbations, the density contrast of the superfluid can grow larger (compared to a normal fluid at the same temperature) and hence produce a stronger net gravitational force acting on the perturber. However, we recall that this result does not include the effect of the critical velocity which would limit this thermal counterflow. In Sect. 4 we propose a scheme to include the critical velocity in linear perturbation theory, and then test the scheme using hydrodynamic simulations in Sect. 6.
3.4. Neglecting selfgravitation
The numerical results presented in Sect. 6, as well as the decay times of globular clusters in Sect. 7, are obtained when selfgravitation is neglected. It is therefore of interest to see what the steadystate linear theory predicts in this case as well.
Neglecting selfgravitation amounts to setting C_{2} = 0. The dispersion relation becomes
For the equation of state used throughout this work, and T/T_{c} ≲ 0.2, the above superfluid sound speeds can be accurately approximated by
We note that for c_{n} ≫ c_{T = 0}, we have . The dynamical friction takes the form
One feature that is clear in this limit is that F_{DF} jumps from zero as V becomes larger than c_{−}, and jumps again as it crosses c_{+}. It seems odd that the force should change value so dramatically when the velocity of the perturber crosses these thresholds, and indeed we find in the numerical simulations in Sect. 5 that it does not. The problem is that in the steadystate case, as considered in this section, the linear overdensity is symmetric upstream and downstream when the perturber moves at subsonic speeds, resulting in a zero net gravitational force at the position of the perturber. This is not an issue at supersonic speeds because the perturber moves faster than the background fluid can respond to the perturbation, which is at the speed of sound, resulting in a clear cone trailing the perturber (Ostriker 1999). At subsonic speeds, on the other hand, the fluid reacts faster than the perturber moves, and with an infinite amount of time to propagate this response, the firstorder perturbation of the background becomes symmetric. In order to overcome this shortcoming of steadystate linear perturbation theory, other studies have broken this symmetry by switching on the perturber for a finite time (Ostriker 1999; SánchezSalcedo 2012), or by going to secondorder perturbations (Lee & Stahler 2011; Shadmehri & Khajenabi 2012). In the following section, the finitetime approach is employed for a superfluid.
4. Dynamical friction from finitetime linear perturbation theory
For the finitetime calculation, Eqs. (26) and (27) are used without selfgravitation, and an approach very similar to the one used by Ostriker (1999) is followed.
The equations can be written in matrix form as
where
and
By diagonalizing matrix A, the coupled set of equations can be transformed into two decoupled wave equations of the form
which are solved using the retarded Green’s function for the wave equation in three dimensions:
For a point source switched on at the origin at t = 0 and moving at speed ,
where H(x) is the Heaviside function, the solution of χ becomes, upon defining s = z − Vt, ℳ_{i} = V/c_{i}, and R^{2} = x^{2} + y^{2},
The resulting overdensity δρ is a weighted sum of χ_{+} and χ_{−}, and the dynamical friction is obtained by integrating the gravitational force due to the overdensity over the whole volume, that is,
In spherical polar coordinates, s = r cos θ = rx and , we get
where we have again introduced an upper and lower cutoff of scales in the integral to avoid UV and IR divergences. The sound speeds c_{+} and c_{−} are the same as the ones given in Eq. (55), and
The dynamical friction from the finitetime calculation is compared to the steadystate result in Fig. 1. The discontinuities have been removed, with the force increasing with velocity V until it reaches a maximum near the sound speed, after which the perturber becomes supersonic and the friction force decreases with the same 1/V^{2} dependence as in the steadystate result. As time passes, the finitetime result approaches the steadystate result, as expected.
Fig. 1. Dynamical friction from linear perturbation theory using the finitetime approach (solid lines) and the steadystate approach (dotted lines) as a function of V. As time passes, the finitetime result, Eq. (69), approaches the steadystate result, Eq. (58). In the zerotemperature limit, we have T = 0, while in the normal fluid case we have ρ_{s} = 0. (a) t = 0.1R_{max}/V. (b) t = R_{max}/V. (c) t = 10R_{max}/V. 
Both approaches predict F_{DF} in the superfluid phase to be very close to the zerotemperature value, even when thermal pressure dominates over the contribution from selfinteractions. However, we must recall again that the Landau criterion is not included in linear perturbation theory, which will limit the thermal counterflow of the superfluid, making it behave more like a normal fluid, thus decreasing the dynamical friction as thermal pressure forces inhibit the growth of density perturbations. Let us therefore construct an ad hoc scheme to include the critical velocity in the linear theory.
The critical velocity is expected to have an effect when the relative velocity is of the order of the critical velocity and larger. Therefore, let us consider the linearized equation for the relative velocity,
The amplitude of δρ and δS, and hence δT, increases with M, driving w up to the critical value faster, causing the effect of the critical velocity on the system to be more prominent. Increasing M should therefore have a similar effect as lowering v_{c} in transitioning the superfluid dynamical friction from the T = 0 value to the fully normal fluid value.
We now assume that for an estimate of the characteristic counterflow of the system, there is an interpolating function with , , and a transitional region around , such that
where and is the dynamical friction from linear theory for the superfluid and fully normal fluid, respectively. Using Eq. (73) we can write
The length L and time Δt are characteristic scales over which the fluid attains the mass and entropy overdensity at the origin, δρ(0) and δS(0). The timescale can be estimated as Δt = L/v, where v is some characteristic velocity in the problem. The largest superfluid sound speed, c_{+}, which is essentially the fastest speed with which the superfluid can respond to disturbances, was found to work.
For the δfunction perturbation, the central values for the mass and entropy overdensities diverge in the linear theory, meaning that δρ(0) and δS(0) are not welldefined. Instead, these should be evaluated at some point near the origin, as was done for dynamical friction. With the equation of state used in this work, an estimate of the linear entropy contrast at R_{min}/2 is
The rough estimate of the counterflow is therefore
Only the form of the interpolating function remains to be specified. The simple but rather arbitrary choice,
was found to work well, as we see in Sect. 6.
5. Numerical simulation of dynamical friction
To test the calculations from linear perturbation theory, the full superfluid equations are integrated numerically. We use the frame comoving with the perturber, meaning that its gravitational field is static and centered at the origin, while the background fluid is moving. We take the perturber to be a sphere with uniform mass density . The system has rotational symmetry, and so cylindrical coordinates are employed; the axial distance is z, the distance along the axis of rotational symmetry, and the radial distance is r, the distance from the axis. The simulation volume is therefore a cylinder, and we take its domain to be −L < z < L and 0 < r < L.
The superfluid is initialized as a uniform fluid moving with velocity . More fluid is injected into the simulation volume with the same velocity at the z = +L boundary. The z = −L and r = L boundaries are taken to have zero gradients, while the inner boundary r = 0 has a reflecting boundary condition.
To numerically integrate the superfluid equations, a Godunov scheme similar to the one described in Hartman et al. (2020) is used. In the present work, the generation of entropy when the Landau criterion is broken is not included. Also, as we evolve the entropy instead of the energy, and we do not include any viscosity, the numerical scheme dissipates kinetic energy at shock fronts that is not converted into internal energy. In the absence of this shock heating, the total energy is not strictly conserved. Nevertheless, we have found that this inaccuracy is by and large negligible for the scenarios we consider here because the solutions are mostly in or near the linear regime.
The selfgravitation of the superfluid is neglected. The gravitational field it produces is computed only to find the resulting force on the perturber, that is, the dynamical friction. The initial fluid parameters are ρ = 2 × 10^{7} M_{⊙} kpc^{−3}, T = 0.2T_{c}, m = 500 eV, and g = 2 × 10^{−3} eV^{−2}. These parameters are chosen only to illustrate the basic features of dynamical friction in superfluids while keeping the simulation runtimes reasonably short. Unless stated otherwise, the size of the perturber is taken to be R_{pert} = 2 pc with mass M = 0.1 M_{⊙}, while the simulation size is L = 150 pc. The simulation is run until t = 10 pc/V, that is, until the background has moved 10 pc. This is small compared to the full simulated length, but is necessary for preventing boundary effects from interfering with the results. R_{min} is taken to be the size of the perturber, R_{min} = 2 pc, and R_{max} the radius of the cylindrical simulation volume, R_{max} = 150 pc, when compared to linear perturbation theory. The resolution of the simulated volume is 4096 × 2048 cells, in the z and r directions, respectively, for the superfluid case. In the zerotemperature and normal fluid limits, for which the numerical scheme was made secondorder in time and space using a MinMod slopelimited MUSCLHancock scheme (Toro 2006) without stability issues, a lower resolution of 2048 × 1024 is used.
An effective critical velocity , which is just v_{c} multiplied by some factor, is used in the simulations to show the effect of varying v_{c} without actually changing other parameters such as the particle mass and selfinteraction.
6. Comparison of perturbation theory and numerical simulation
In Fig. 2 the dynamical friction from the numerical simulations is compared to the linear result with the effect of the critical velocity included, Eqs. (69), (74), and (78). Even with the Landau criterion, given by Eq. (18), the dynamical friction in the superfluid can be very similar to the zerotemperature limit, as was shown in the linear theory. This similarity can also be seen in the mass density profile shown in Fig. 3. At T = 0, for which the pertuber is supersonic with V = 1.5c_{T = 0}, there is a welldefined supersonic cone that trails the pertuber, and the density contrast reaches about 4.5. The finitetemperature superfluid has a similar density contrast and supersonic cone, though not as welldefined, illustrating that the superfluid behaves like the T = 0 limit as thermal counterflow suppresses thermal perturbations. In the fully normal fluid case, the density contrast is much smaller, around 0.17, and the perturber is instead moving at subsonic speeds, because V = 1.5c_{T = 0} < c_{n}, hence there is no sonic cone. As is decreased, the relative velocity becomes increasingly limited and the thermal counterflow inefficient, causing the superfluid density profile to approach the fully normal fluid limit. The dynamical friction changes accordingly, as shown in Fig. 4. Furthermore, Fig. 5 shows the friction force as a function of the mass of the perturber M, confirming the expectation that increasing M causes a transition from superfluid to normal fluid behavior in a similar manner to decreasing v_{c}.
Fig. 2. Dynamical friction against velocity for the superfluid with varying , and for the zero temperature and the fully normal fluid limits. The results from finitetime linear perturbation theory are shown with dotted lines of the same colors. Even with the critical velocity included, the superfluid case gives a dynamical friction force of the same magnitude as the zero temperature limit. When is decreased, the superfluid approaches the fully normal fluid limit as thermal counterflow is increasingly limited. The sound speeds c_{T = 0} and c_{n} are indicated by the vertical dotted lines. 
Fig. 3. Density profiles and streamlines for V = 1.5c_{T = 0}. The mass density profiles are superimposed by the net mass density velocity, v = j/ρ, while the entropy density is superimposed by the relative velocity w = v_{s} − v_{n}. The perturber has mass M = 5 M_{⊙}, the simulation volume is L = 75 pc, and the time is t = 50 pc/V. (a) T = 0 limit. (b) ρ_{s} = 0 limit. (c) . (d) . (e) . 
Fig. 4. Dynamical friction against the effective critical velocity , for V = 1.5c_{T = 0}, with the results from finitetime linear perturbation theory included with dotted lines of the same colors. As is lowered, the dynamical friction goes from about the same value as the zero temperature limit to the value in the fully normal fluid limit, changing by about two orders of magnitude. 
The numerical results of Figs. 2, 4, and 5 show that the scheme to include v_{c} in the linear theory (Eqs. (74) and (78)) successfully captures the basic dependence on the perturber mass and critical velocity, though it fails at low velocities, V < c_{T = 0} ≈ c_{−}, suggesting that other factors might come into play at those speeds. However, as we see in the following section, this does not cause any problems when applied to the Fornax dSph, as in the relevant parameter space we have , which is far away from the transition between the superfluid and normal fluid phase, and therefore no interpolation is needed. No further attempt was therefore made to improve the scheme.
Fig. 5. Dynamical friction against the perturber mass M, for V = 1.5c_{T = 0}, with the results from finitetime linear perturbation theory included with dotted lines of the same colors. The departure from perturbation theory for the zerotemperature case at high M is due to nonlinear effects. Increasing the mass of the perturber causes the superfluid to behave increasingly like a normal fluid, similarly to the effect of decreasing . 
7. Application to the Fornax system
So far, only the physics of dynamical friction in superfluids has been discussed, with little reference to the real world. Now, armed with the expressions derived and tested in the previous sections, the parameter space of superfluid DM–the particle mass, selfinteraction, and temperature–can be explored by estimating the orbital decay time of GCs in the Fornax dSph, and checking whether the timing problem is alleviated for SIBECDM, or exacerbated.
The decay time can be defined as the time τ_{DF} it takes dynamical friction to reduce the angular momentum L of the GCs to zero;
where M, V, and r are the mass, circular orbital velocity, and the orbital radius of the GCs. The density profile of the Fornax dSph is modeled in the same way as in Cole et al. (2012), using a spherical doublepowerlaw profile^{1} of the form
The profile parameters, still following Cole et al. (2012), are listed in Table 1, and correspond to different models for the shape of the Fornax dSph. As SIBECDM, like many alternative theories for DM, is in part motivated by typically having a cored profile, we only focus on the Large core (LC) and Weak cusp (WC) models from Cole et al. (2012). It should also be noted that the density profile Eq. (80) models the total mass density, that is, both stellar and DM, but as DM is the dominant component, we use it as a pure DM profile. As illustrated in Fig. 6, subtracting a subdominant portion of the density ρ_{0} = ρ(r) in the computation of the SIBECDM dynamical friction in order to account for the presence of stellar mass does not significantly alter the value of the orbital decay time.
Fig. 6. Change in the orbital decay time as parameters related to the modeling of the Fornax dSph, the GCs, and the dynamical friction are varied. The reference values, which are for GC3 in the LC model, are labeled with the subscript “ref”. 
Halo mass profile parameters, using Eq. (80), with values from Cole et al. (2012).
Estimates of the masses, projected orbital radii r_{⊥}, and core radii r_{c} of the GCs, which we use as R_{min} in perturbation theory, are listed in Table 2. As in Lancaster et al. (2020) and Hui et al. (2017), is used as the “true” radial distance from the Fornax center. This larger radial distance leads to a longer estimate of the decay time τ_{DF}, as illustrated in Fig. 6. The orbital velocities of the GCs are assumed to be circular, determined by the total halo mass enclosed by their orbits, M_{encl},
Projected radial distances, masses, and core radii of the GCs (not including the sixth found by Wang et al. 2019b) in the Fornax dSph, taken from Mackey & Gilmore (2003).
Inside r_{s} the density profile of the Fornax dSph is approximately constant and cored for the LC and WC models. Hence, we use r_{s} as the “core size” of the Fornax, R_{c}, and the upper length scale when we use perturbation theory, R_{max}. The DM density is determined at the position of the GCs using Eq. (80).
There is a limited region of parameter space that is both physically relevant, and may provide a reasonable estimate of τ_{DF}. This region should satisfy the following:

The core radius of the halo obtained from hydrostatic equilibrium should not exceed the core radius of the dSph as modeled by Eq. (80).

The DM mass and selfinteraction should satisfy constraints from observations of the deceleration of DM in cluster collisions (Harvey et al. 2015).

The relaxation rate of DM should be higher than the rate of dynamical changes in the dSph, so that the system can thermalize and form a superfluid.

Perturbation theory is only properly valid for δρ/ρ ≪ 1.

The ad hoc scheme to include the critical velocity introduced at the end of Sect. 4 failed to reproduce the numerical results of Sect. 6 for velocities V < c_{−}. Therefore, we cannot trust the dynamical friction obtained from linear perturbation theory for these velocities. However, this should only be a problem near the transition , where the form of interpolation between the superfluid and normal fluid result is important.
While this list is not exhaustive, it provides a minimum set of criteria that should be fulfilled. Due to our ignorance of the general behavior of superfluid DM in a number of situations, we enforce relaxed variants of the above constraints.
As seen in the previous sections, and shown in an earlier work (Hartman et al. 2020), counterflow can effectively redistribute thermal energy in a superfluid. Therefore, the shape of the temperature profile of a realistic superfluid DM halo is unknown. The least constraining assumption is therefore made; that the counterflow has washed out any significant thermal differences, so that only the interaction part of the pressure (the only pressure present at T = 0) determines the hydrostatic profile. Demanding that the core radius of the halo be larger than the core radius obtained from hydrostatic equilibrium, which we define as ρ(R_{c})≈ρ(0)/2, gives
This relaxed constraint is only possible if it is physically feasible for the counterflow to transport a significant portion of the thermal energy away from the halo core. A supplementary criterion can be derived by demanding that the total entropy flux due to thermal counterflow at the core edge R_{c}, with w = v_{c}, be of the same order as the total entropy enclosed in R_{c}. This leads to
where Δt should be smaller than the age of the dSph, for example Δt ∼ 1 Gyr. As shown in Fig. 7, the difference between the T = 0 and the finite temperature treatment of the hydrostatic halo size can be very large, and we do not expect a realistic superfluid halo to be able to completely remove thermal differences, even if upper estimates of the thermal counterflow suggest it could. A realistic superfluid core radius is therefore expected to be larger than the zerotemperature estimate used to derive Eq. (82).
Fig. 7. Criteria listed in the text, and the orbital decay timescale for GC3 at T/T_{c} = 10^{−4} in the LC model for reference (solid black line). Permitted parameter space; the left side is from the constraint on the halo core radius in hydrostatic equilibrium, Eq. (82); the upper right side from the constraint from galaxy cluster collisions, Eq. (84); and the lower right side from the minimum relaxation rate needed to thermalize the fluid across the halo, Eq. (86) (solid blue line). V = c_{−}, with V < c_{−} on the left side (dotted blue line). Criterion for linear perturbation theory to be properly valid, with δρ/ρ_{0} < 1 satisfied on the left side (dashed blue line). Supplementary criterion for the T = 0 treatment of the hydrostatic halo size, with Eq. (83) satisfied on the left side (solid red line). , where the superfluid dynamical friction transitions from superfluid on the left side, to normal fluid on the right (dashed red line). Constraint on the halo core radius in hydrostatic equilibrium with thermal pressure included, with halo cores smaller than the core as modeled by Eq. (80) on the right side. 
By measuring the spatial offset of stars, gas, and DM in colliding galaxy clusters, a constraint on the selfinteraction cross section of DM, σ, can be established (Harvey et al. 2015). The lack of deceleration of DM and its proximity to the collisionless stars in these collisions places an upper limit of σ/m ≲ 0.5 cm^{−2} g^{−1}. In terms of the selfinteraction parameter g, this constraint reads (Pitaevskii & Stringari 2016)
While the above places upper limits on g, there is also a lower limit that must be considered given by the criterion that the DM superfluid should be thermalized across much of the halo core. For this we require the relaxation rate of DM, Γ_{DM}, to be higher than the rate of dynamical changes in the halo core, . For twobody interactions, the relaxation rate is Γ ∼ nσδv, where σ is the scattering crosssection and δv the velocity dispersion of the particles. In terms of g, as above, the crosssection is σ = m^{2}g^{2}/4π. However, for a condensed Bose gas, the relaxation rate is enhanced, that is, Γ ∼ 𝒩nσδv, where
because of the high occupation number of the ground state (Sikivie & Yang 2009). Using δv ∼ V, that is, that the DM velocity dispersion is of the same order as the GC orbital velocity, the criterion Γ_{DM} > Γ_{grav} becomes
It should be noted that the enhancement factor is included in this criterion, but not in the constraint from cluster collisions. This is another example of a relaxed constraint due to our ignorance of how the superfluid properties might change in the various situations. The characteristic speeds of cluster collisions are usually much higher than inside halos, which might result in a much larger disruption of the DM BEC. Furthermore, the DM fluid may not even be condensed throughout most of the cluster, only inside dense structures. We therefore choose the least restrictive constraint by including 𝒩 inside the dSph DM halo, but not outside.
The remaining constraints due to δρ/ρ ≪ 1 and V < c_{−} are readily obtained from perturbation theory and Eq. (57). The result from the finitetime approach, Eq. (69), with our proposed scheme for including the critical velocity, Eqs. (74) and (78), is used to compute the dynamical friction. The characteristic timescale t = r/V is used as the finite time the perturber has been active, though the results are not sensitive to this choice. A deficiency of the finitetime formalism is the lack of selfgravitation in the background fluid, and it might seem that a better choice is to instead use the steadystate expression, Eq. (43), which includes this property. However, that result assumes the perturber has acted on an otherwise static background for an infinite amount time, and does not take into account that the background can be rotationally supported, and therefore resist the largescale gravitational collapse induced by the perturber. Furthermore, numerical studies of dynamical friction in realistic halos show that linear theory can provide reasonable estimates of the gravitational drag force even with selfgravitation neglected if the mass of the perturber is significantly smaller than the mass of the host halo (Fujii et al. 2006; Aceves & Colosimo 2007; Binney & Tremaine 2008; Chapon et al. 2013; Antonini & Merritt 2011; Tamfal et al. 2020), as is the case here. However, because linear perturbation theory assumes a uniform background with an upper cutoff of scales to take into account the finite extent of the background, we focus on the GCs located inside r_{s}, where the density profile of the LC and WC models are approximately flat, which are GC3 and GC4. These are also the ones that the timingproblem usually applies too (Cole et al. 2012; Hui et al. 2017; ArcaSedda & CapuzzoDolcetta 2017).
The criteria on the parameter space listed above are illustrated in Fig. 7 for GC3 in the LC profile, with the estimated orbital decay time for reference. Some features that also hold for the other cases shown in Fig. 8 are worth pointing out. First, the transition point between superfluid and normal fluid behavior, , lies far away from the region δρ/ρ_{0} < 1 where perturbation theory is valid, meaning that we do not need to worry about the accuracy of the interpolation scheme described in Sect. 4. Second, the decay time becomes very large for V < c_{−}, because, as we have seen in the previous sections, the dynamical friction vanishes quickly for velocities below the lowest sound speed.
Fig. 8. Decay time of GC3 and GC4, as listed in Table 2, in the LC and WC models for the Fornax dSph density profile from Table 1 (solid line). Permitted parameter space; the left side is from the constraint on the halo core radius in hydrostatic equilibrium, Eq. (82); the upper right side from the constraint from galaxy cluster collisions, Eq. (84); and the lower right side from the minimum relaxation rate needed to thermalize the fluid across the halo, Eq. (86) (dotted line). Criterion for linear perturbation theory to be properly valid, with δρ/ρ_{0} < 1 satisfied on left side (dashed line). Limit due to a hydrostatic halo with thermal pressure included, with resulting core radii smaller than the core of the Fornax dSph as modeled by Eq. (80) to the right. Changing the temperature only changes the decay time of the normal fluid phase, as well as the crossover from superfluid to normal fluid. However, for the temperatures shown and lower, the normal phase is well outside the parameter space where perturbation theory is valid. (a) GC3 LC, T/T_{c} = 10^{−2}. (b) GC3 LC, T/T_{c} = 10^{−4}. (c) GC3 LC, T/T_{c} = 10^{−6}. (d) GC3 WC, T/T_{c} = 10^{−2}. (e) GC3 WC, T/T_{c} = 10^{−4}. (f) GC3 WC, T/T_{c} = 10^{−6}. (g) GC4 LC, T/T_{c} = 10^{−2}. (h) GC4 LC, T/T_{c} = 10^{−4}. (i) GC4 LC, T/T_{c} = 10^{−6}. (j) GC4 WC, T/T_{c} = 10^{−2}. (k) GC4 WC, T/T_{c} = 10^{−4}. (l) GC4 WC, T/T_{c} = 10^{−6}. 
The orbital decay time for a wide range of parameters is shown in Fig. 8 for the two GCs inside the core radius of the Fornax dSph, GC3 and GC4, in the two density profiles considered. τ_{DF} generally either attains a minimum value, τ_{DF, min}, or approaches infinity. The minimum values in the region V > c_{−} and δρ/ρ_{0} < 1 are summarized in Table 3, with τ_{DF} in the range 67 Myr−197 Myr. These timescales are considerably smaller than the CDM result assuming the same density profiles, 515 Myr−1327 Myr, with the dynamical friction given by Binney & Tremaine (2008)
where Λ ≈ rδv^{2}/GM, , δv is the velocity dispersion of CDM particles, taken to be δv ≈ V, and erf is the error function.
The decay time remains small even if parameters used to model the Fornax dSph, the GCs, and the dynamical friction are varied, as illustrated in Fig. 6. A notable exception is the position of the GC, for which τ_{DF} is considerably shorter when closer to the halo center, and likewise longer when further away. This implies that the value for τ_{DF} obtained from Eq. (79) overestimates the time it takes the GC to fully decay from its current position, but it also implies that the migration towards the halo center was slower in the past when the GCs were at larger radial distances. Indeed, estimates of τ_{DF, min} for GC1, GC2, and GC5, all of which are located at r ≳ 1 kpc, give decay times in excess of 4 Gyr. In the CDM case, the decay times for these GCs are even longer: 17 Gyr and more. These estimates do not suggest a timing problem for the outer GCs, even if their decay times are considerably shorter for SIBECDM compared to CDM. However, we note that these GCs are near or outside the radius r_{s}, where the density profile of the dSph falls sharply, and therefore we do not expect the result for the dynamical friction, nor Eq. (79), to necessarily provide a reasonable estimate of τ_{DF}. Nonetheless, the present results show that for a large region of the relevant parameter space of the SIBECDM model considered here, GC3 and GC4 are currently racing towards the center of their host halo in a SIBECDM universe.
Let us now consider τ_{DF} in light of constraints on the SIBECDM model from the literature. By fitting rotation curves of slowly rotating SIBECDM halos in hydrostatic equilibrium in 173 nearby galaxies from the Spitzer Photomery & Accurate Rotation Curves (SPARC) data (Lelli et al. 2016), Crǎciun & Harko (2020) estimated the properties of SIBECDM halos at T = 0, and found the best fit values for g/m^{2} to be between 2.7 × 10^{−4} eV^{−4} and 5.0 × 10^{−2} eV^{−4}. For reference, the estimated limit from hydrostatic equilibrium using Eq. (82) gives g/m^{2} of less than about 2 × 10^{−4} eV^{−4} or 10^{−3} eV^{−4}, depending on the profile used for the dSph. As the preferred values obtained by Crǎciun & Harko (2020) for zerotemperature SIBECDM satisfy V < c_{−} for the GCs and dSph profiles considered, leading to a vanishing dynamical friction, the T = 0 case does not have a timingproblem, a result that could also have been found using heuristic arguments; if the halo is largely supported by hydrostatic pressure, that is, its Jeans’ length is of the order of the DM halo core radius R_{c}, then density perturbations on smaller scales inside the halo will be highly suppressed, resulting in very weak dynamical friction, and therefore long decay times.
In a finitetemperature SIBECDM halo–for which we expect the preferred values for g/m^{2} obtained from fitting rotation curves to be lowered, as it provides additional pressure forces to support DM halos–the present results instead suggest that overly large orbital decay rates due to strong dynamical friction may arise. This is the opposite of what one would naively expect if the superfluid had been treated as a conventional thermal fluid, because an increased pressure generally leads to a smaller maximum dynamical friction. Instead, the superfluid essentially ignores the thermal contribution, and responds to a perturber as if it were at T = 0, which can yield a much stronger friction force.
8. Conclusion
We investigated the dynamical friction acting on an object due to a superfluid background, starting with steadystate linear perturbation theory. The wellknown issue of discontinuities in the friction force as the perturber’s velocity crosses the fluid sound speed was encountered. We therefore also employed a finitetime formalism, which removed these discontinuities, agreeing with previous studies that the dynamical friction increases with the velocity of the perturber until the sound speed is reached, after which the force decreases with the same V^{−2} dependence as the steadystate result. Both approaches predict the force in the superfluid phase to be very similar to the T = 0 limit, even when there are large thermal contributions, yielding a much stronger friction force than one might naively expect when compared to a conventional fluid at the same temperature. This happens because counterflow conspires against thermal perturbations, allowing the superfluid to respond to a perturbation as if it were at zero temperatures. However, the counterflow is only effective as long as it does not exceed the critical velocity v_{c}, which acts as an upper limit. For flows where the counterflow would normally exceed, but is limited by, the critical velocity, the superfluid instead behaves as a normal fluid. Therefore, decreasing v_{c} essentially causes a transition from a superfluid to a normal fluid, interpolating the dynamical friction from about the value at T = 0 to the value of the normal fluid, which can differ by several orders of magnitude. Numerical simulations were also used to investigate dynamical friction, confirming the general dependence of the force on the critical velocity and the mass of the perturber, which was found using linear perturbation theory. However, the linear theory failed to reproduce the shape of the superfluidnormal fluid transition for velocities smaller than the smallest sound speed, V < c_{−}.
Finally, the superfluid dynamical friction was applied to the Fornax dSph and two of its GCs. It was found that the relevant parameter space in which, among other things, perturbation theory is valid gives orbital decay times for these GCs that are much smaller than the age of the dSph, except in the region preferred in the literature (Crǎciun & Harko 2020). The present work therefore suggests that there is no timing problem for Fornax GCs in the SIBECDM model for the values of g/m^{2} obtained by Crǎciun & Harko (2020) by fitting rotation curves at T = 0. For a finitetemperature SIBECDM, on the other hand, for which the preferred parameter space of g/m^{2} is likely lowered, very large decay rates of Fornax GCs pose a problem.
The use of linear perturbation theory made it possible to probe a large region of parameter space that is difficult to explore with numerical simulations. The main limitations of the numerical scheme used in this work are the low order of the Godunov scheme used; the absence of entropy production, both when the critical velocity was enforced and in shock waves, which leads to the total energy not being strictly conserved; and the large difference between the superfluid sound speeds and dynamics, which results in very small timestepping and hence excessive diffusion of the numerical solution. All of these limit the parameters for which we can be confident that the numerical solution is correct, and therefore limits the range within which perturbation theory can be tested. Ideally, superfluid dynamical friction would have also been explored using simulations with realistic models for both the DM halo and perturber, as has been done for galaxies with standard CDM and gas (Chapon et al. 2013; Tamfal et al. 2020), but such a study requires an improved scheme for solving the superfluid hydrodynamics equations.
There appears to be a sign error in Eq. (1) in Cole et al. (2012) when comparing the resulting profiles to their own figures, as well as compared to other works that use the same type of profile (Zhao 1996; Walker et al. 2009; Hague & Wilkinson 2013).
Acknowledgments
We thank the Research Council of Norway for their support, and Benjamin Elder for the discussions that initiated this work. We also thank the anonymous referee for their helpful comments and suggestions that greatly improved this manuscript.
References
 Aceves, H., & Colosimo, M. 2007, Am. J. Phys., 75, 139 [Google Scholar]
 Andersen, J. O. 2004, Rev. Mod. Phys., 76, 599 [Google Scholar]
 Antonini, F., & Merritt, D. 2011, ApJ, 745, 83 [Google Scholar]
 ArcaSedda, M., & CapuzzoDolcetta, R. 2017, MNRAS, 464, 3060 [NASA ADS] [CrossRef] [Google Scholar]
 BarOr, B., Fouvry, J.B., & Tremaine, S. 2019, ApJ, 871, 28 [Google Scholar]
 Barausse, E. 2007, MNRAS, 382, 826 [CrossRef] [Google Scholar]
 Barenghi, C. F., Skrbek, L., & Sreenivasan, K. R. 2014, Proc. Natl. Acad. Sci., 111, 4647 [Google Scholar]
 Battaglia, G., Helmi, A., & Breddels, M. 2013, New Astron. Rev., 57, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhiani, L., & Khoury, J. 2015, Phys. Rev. D, 92, 103510 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhiani, L., Elder, B., & Khoury, J. 2019, JCAP, 2019, 074 [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press) [Google Scholar]
 Boldrini, P., Mohayaee, R., & Silk, J. 2019, MNRAS, 485, 2546 [Google Scholar]
 Boldrini, P., Miki, Y., Wagner, A. Y., et al. 2020, MNRAS, 492, 5218 [Google Scholar]
 BoylanKolchin, M., Ma, C.P., & Quataert, E. 2008, MNRAS, 383, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., & BoylanKolchin, M. 2017, ARA&A, 55, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, ApJ, 97, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S., Hoyos, C., & Oz, Y. 2014, J. High Energy Phys., 2014, 27 [CrossRef] [Google Scholar]
 Chapon, D., Mayer, L., & Teyssier, R. 2013, MNRAS, 429, 3114 [CrossRef] [Google Scholar]
 Clesse, S., & GarcíaBellido, J. 2018, Phys. Dark Univ., 22, 137 [Google Scholar]
 Cole, D. R., Dehnen, W., Read, J. I., & Wilkinson, M. I. 2012, MNRAS, 426, 601 [Google Scholar]
 Colpi, M., Mayer, L., & Governato, F. 1999, ApJ, 525, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Cowsik, R., Wagoner, K., Berti, E., & Sircar, A. 2009, ApJ, 699, 1389 [Google Scholar]
 Crǎciun, M., & Harko, T. 2020, Eur. Phys. J. C, 80, 735 [EDP Sciences] [Google Scholar]
 Darve, C., Bottura, L., Patankar, N. A., & Van Sciver, S. 2012, AIP Conf. Proc., 1434, 247 [Google Scholar]
 Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Debattista, V. P., & Sellwood, J. A. 2000, ApJ, 543, 704 [NASA ADS] [CrossRef] [Google Scholar]
 del Pino, A., Hidalgo, S. L., Aparicio, A., et al. 2013, MNRAS, 433, 1505 [NASA ADS] [CrossRef] [Google Scholar]
 Del Popolo, A., & Le Delliou, M. 2017, Galaxies, 5, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Doi, D., Shirai, Y., & Shiotsu, M. 2008, AIP Conf. Proc., 985, 648 [Google Scholar]
 Dosopoulou, F., & Antonini, F. 2017, ApJ, 840, 31 [Google Scholar]
 Elbert, O. D., Bullock, J. S., GarrisonKimmel, S., et al. 2015, MNRAS, 453, 29 [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Relativ., 15, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Fujii, M., Funato, Y., & Makino, J. 2006, PASJ, 58, 743 [Google Scholar]
 Goerdt, T., Moore, B., Read, J. I., Stadel, J., & Zemp, M. 2006, MNRAS, 368, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Goerdt, T., Moore, B., Read, J. I., & Stadel, J. 2010, ApJ, 725, 1707 [NASA ADS] [CrossRef] [Google Scholar]
 Gómez, L. G., & Rueda, J. 2017, Phys. Rev. D, 96, 063001 [Google Scholar]
 Hague, P. R., & Wilkinson, M. I. 2013, MNRAS, 433, 2314 [NASA ADS] [CrossRef] [Google Scholar]
 Harko, T., & Mocanu, G. 2012, Phys. Rev. D, 85, 084012 [CrossRef] [Google Scholar]
 Harko, T., Liang, P., Liang, S.D., & Mocanu, G. 2015, JCAP, 2015, 027 [Google Scholar]
 Hartman, S. T. H., Winther, H. A., & Mota, D. F. 2020, A&A, 639, A90 [EDP Sciences] [Google Scholar]
 Harvey, D., Massey, R., Kitching, T., Taylor, A., & Tittley, E. 2015, Science, 347, 1462 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [Google Scholar]
 Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541 [Google Scholar]
 Jiang, C. Y., Jing, Y. P., Faltenbacher, A., Lin, W. P., & Li, C. 2008, ApJ, 675, 1095 [NASA ADS] [CrossRef] [Google Scholar]
 Just, A., Khan, F. M., Berczik, P., Ernst, A., & Spurzem, R. 2011, MNRAS, 411, 653 [NASA ADS] [CrossRef] [Google Scholar]
 Katz, A., Kurkela, A., & Soloviev, A. 2019, JCAP, 2019, 017 [Google Scholar]
 Kaur, K., & Sridhar, S. 2018, ApJ, 868, 134 [Google Scholar]
 Khoury, J. 2016, Phys. Rev. D, 93, 103533 [NASA ADS] [CrossRef] [Google Scholar]
 Lancaster, L., Giovanetti, C., Mocz, P., et al. 2020, JCAP, 2020, 001 [CrossRef] [Google Scholar]
 Landau, L. 1941, Phys. Rev., 60, 356 [Google Scholar]
 Lee, A. T., & Stahler, S. W. 2011, MNRAS, 416, 3177 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, A. T., & Stahler, S. W. 2014, A&A, 561, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157 [Google Scholar]
 Leung, G. Y. C., Leaman, R., van de Ven, G., & Battaglia, G. 2020, MNRAS, 493, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Lovell, M. R., Frenk, C. S., Eke, V. R., et al. 2014, MNRAS, 439, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Mackey, A. D., & Gilmore, G. F. 2003, MNRAS, 340, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Madelung, E. 1926, Naturwissenschaften, 14, 1004 [Google Scholar]
 Meadows, N., Navarro, J. F., SantosSantos, I., BenítezLlambay, A., & Frenk, C. 2020, MNRAS, 491, 3336 [Google Scholar]
 Milgrom, M. 1983a, ApJ, 270, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983b, ApJ, 270, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983c, ApJ, 270, 371 [Google Scholar]
 Mina, M., Mota, D. F., & Winther, H. A. 2020a, A&A, 641, A107 [EDP Sciences] [Google Scholar]
 Mina, M., Mota, D. F., & Winther, H. A. 2020b, ArXiv eprints [arXiv:2007.04119] [Google Scholar]
 Mocz, P., Vogelsberger, M., Robles, V. H., et al. 2017, MNRAS, 471, 4559 [NASA ADS] [CrossRef] [Google Scholar]
 Mulder, W. A. 1983, A&A, 117, 9 [NASA ADS] [Google Scholar]
 Navez, P., & Graham, R. 2006, Phys. Rev. A, 73, 043612 [CrossRef] [Google Scholar]
 Nori, M., & Baldi, M. 2018, MNRAS, 478, 3935 [CrossRef] [Google Scholar]
 Nori, M., & Baldi, M. 2021, MNRAS, 501, 1539 [Google Scholar]
 Oh, K. S., Lin, D. N. C., & Richer, H. B. 2000, ApJ, 531, 727 [Google Scholar]
 Ostriker, E. C. 1999, ApJ, 513, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Pani, P. 2015, Phys. Rev. D, 92, 123530 [Google Scholar]
 Percival, W. J., Baugh, C. M., BlandHawthorn, J., et al. 2001, MNRAS, 327, 1297 [NASA ADS] [CrossRef] [Google Scholar]
 Pitaevskii, L. P., & Stringari, S. 2016, BoseEinstein Condensation and Superfluidity (Oxford: Oxford University Press) [CrossRef] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [Google Scholar]
 Read, J. I., Walker, M. G., & Steger, P. 2019, MNRAS, 484, 1401 [Google Scholar]
 Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Sales, L. V., Navarro, J. F., Oman, K., et al. 2016, MNRAS, 464, 2419 [Google Scholar]
 SánchezSalcedo, F. J. 2012, ApJ, 745, 135 [Google Scholar]
 SánchezSalcedo, F. J., & Brandenburg, A. 1999, ApJ, 522, L35 [NASA ADS] [CrossRef] [Google Scholar]
 SantosSantos, I. M., Brook, C. B., Stinson, G., et al. 2015, MNRAS, 455, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931 [NASA ADS] [CrossRef] [Google Scholar]
 Schive, H.Y., Chiueh, T., & Broadhurst, T. 2014, Nat. Phys., 10, 496 [Google Scholar]
 Schwabe, B., Niemeyer, J. C., & Engels, J. F. 2016, Phys. Rev. D, 94, 043513 [Google Scholar]
 Sellwood, J. A. 2014, Rev. Mod. Phys., 86, 1 [Google Scholar]
 Shadmehri, M., & Khajenabi, F. 2012, MNRAS, 424, 919 [Google Scholar]
 Shao, S., Gao, L., Theuns, T., & Frenk, C. S. 2013, MNRAS, 430, 2346 [Google Scholar]
 Sharma, A., Khoury, J., & Lubensky, T. 2019, JCAP, 2019, 054 [Google Scholar]
 Sikivie, P., & Yang, Q. 2009, Phys. Rev. Lett., 103, 111301 [Google Scholar]
 Skrbek, L. 2011, J. Phys.: Conf. Ser., 318, 012004 [Google Scholar]
 Skrbek, L., & Sreenivasan, K. R. 2012, Phys. Fluids, 24, 011301 [CrossRef] [Google Scholar]
 Slepian, Z., & Goodman, J. 2012, MNRAS, 427, 839 [Google Scholar]
 Soulaine, C., Quintard, M., Baudouy, B., & Van Weelderen, R. 2017, Phys. Rev. Lett., 118, 074506 [Google Scholar]
 Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [Google Scholar]
 Strigari, L. E. 2018, Rep. Progr. Phys., 81, 056901 [Google Scholar]
 Tamfal, T., Mayer, L., Quinn, T. R., et al. 2020, ApJ, submitted [arXiv:2007.13763] [Google Scholar]
 Taylor, E., & Griffin, A. 2005, Phys. Rev. A, 72, 8739 [Google Scholar]
 Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, ApJ, 606, 702 [NASA ADS] [CrossRef] [Google Scholar]
 Thun, D., Kuiper, R., Schmidt, F., & Kley, W. 2016, A&A, 589, A10 [EDP Sciences] [Google Scholar]
 Toro, E. 2006, Appl. Numer. Math., 56, 1464 [CrossRef] [Google Scholar]
 TrujilloGomez, S., Klypin, A., Primack, J., & Romanowsky, A. J. 2011, ApJ, 742, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Tulin, S., & Yu, H.B. 2018, Phys. Rep., 730, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, M. 2013, in Planets, Stars and Stellar Systems: Volume 5: Galactic Structure and Stellar Populations, eds. T. D. Oswalt, & G. Gilmore (Dordrecht: Springer), 1039 [Google Scholar]
 Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2009, ApJ, 704, 1274 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, M. Y., Boer, T. D., Pieres, A., et al. 2019a, ApJ, 881, 118 [Google Scholar]
 Wang, M. Y., Koposov, S., DrlicaWagner, A., et al. 2019b, ApJ, 875, L13 [Google Scholar]
 Weinberg, M. D. 1985, MNRAS, 213, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, D. H., Bullock, J. S., Governato, F., Naray, R. K. D., & Peter, A. H. G. 2015, Proc. Natl. Acad. Sci., 112, 12249 [Google Scholar]
 Zhao, H. 1996, MNRAS, 278, 488 [Google Scholar]
 Zhu, Q., Marinacci, F., Maji, M., et al. 2016, MNRAS, 458, 1559 [Google Scholar]
All Tables
Halo mass profile parameters, using Eq. (80), with values from Cole et al. (2012).
Projected radial distances, masses, and core radii of the GCs (not including the sixth found by Wang et al. 2019b) in the Fornax dSph, taken from Mackey & Gilmore (2003).
All Figures
Fig. 1. Dynamical friction from linear perturbation theory using the finitetime approach (solid lines) and the steadystate approach (dotted lines) as a function of V. As time passes, the finitetime result, Eq. (69), approaches the steadystate result, Eq. (58). In the zerotemperature limit, we have T = 0, while in the normal fluid case we have ρ_{s} = 0. (a) t = 0.1R_{max}/V. (b) t = R_{max}/V. (c) t = 10R_{max}/V. 

In the text 
Fig. 2. Dynamical friction against velocity for the superfluid with varying , and for the zero temperature and the fully normal fluid limits. The results from finitetime linear perturbation theory are shown with dotted lines of the same colors. Even with the critical velocity included, the superfluid case gives a dynamical friction force of the same magnitude as the zero temperature limit. When is decreased, the superfluid approaches the fully normal fluid limit as thermal counterflow is increasingly limited. The sound speeds c_{T = 0} and c_{n} are indicated by the vertical dotted lines. 

In the text 
Fig. 3. Density profiles and streamlines for V = 1.5c_{T = 0}. The mass density profiles are superimposed by the net mass density velocity, v = j/ρ, while the entropy density is superimposed by the relative velocity w = v_{s} − v_{n}. The perturber has mass M = 5 M_{⊙}, the simulation volume is L = 75 pc, and the time is t = 50 pc/V. (a) T = 0 limit. (b) ρ_{s} = 0 limit. (c) . (d) . (e) . 

In the text 
Fig. 4. Dynamical friction against the effective critical velocity , for V = 1.5c_{T = 0}, with the results from finitetime linear perturbation theory included with dotted lines of the same colors. As is lowered, the dynamical friction goes from about the same value as the zero temperature limit to the value in the fully normal fluid limit, changing by about two orders of magnitude. 

In the text 
Fig. 5. Dynamical friction against the perturber mass M, for V = 1.5c_{T = 0}, with the results from finitetime linear perturbation theory included with dotted lines of the same colors. The departure from perturbation theory for the zerotemperature case at high M is due to nonlinear effects. Increasing the mass of the perturber causes the superfluid to behave increasingly like a normal fluid, similarly to the effect of decreasing . 

In the text 
Fig. 6. Change in the orbital decay time as parameters related to the modeling of the Fornax dSph, the GCs, and the dynamical friction are varied. The reference values, which are for GC3 in the LC model, are labeled with the subscript “ref”. 

In the text 
Fig. 7. Criteria listed in the text, and the orbital decay timescale for GC3 at T/T_{c} = 10^{−4} in the LC model for reference (solid black line). Permitted parameter space; the left side is from the constraint on the halo core radius in hydrostatic equilibrium, Eq. (82); the upper right side from the constraint from galaxy cluster collisions, Eq. (84); and the lower right side from the minimum relaxation rate needed to thermalize the fluid across the halo, Eq. (86) (solid blue line). V = c_{−}, with V < c_{−} on the left side (dotted blue line). Criterion for linear perturbation theory to be properly valid, with δρ/ρ_{0} < 1 satisfied on the left side (dashed blue line). Supplementary criterion for the T = 0 treatment of the hydrostatic halo size, with Eq. (83) satisfied on the left side (solid red line). , where the superfluid dynamical friction transitions from superfluid on the left side, to normal fluid on the right (dashed red line). Constraint on the halo core radius in hydrostatic equilibrium with thermal pressure included, with halo cores smaller than the core as modeled by Eq. (80) on the right side. 

In the text 
Fig. 8. Decay time of GC3 and GC4, as listed in Table 2, in the LC and WC models for the Fornax dSph density profile from Table 1 (solid line). Permitted parameter space; the left side is from the constraint on the halo core radius in hydrostatic equilibrium, Eq. (82); the upper right side from the constraint from galaxy cluster collisions, Eq. (84); and the lower right side from the minimum relaxation rate needed to thermalize the fluid across the halo, Eq. (86) (dotted line). Criterion for linear perturbation theory to be properly valid, with δρ/ρ_{0} < 1 satisfied on left side (dashed line). Limit due to a hydrostatic halo with thermal pressure included, with resulting core radii smaller than the core of the Fornax dSph as modeled by Eq. (80) to the right. Changing the temperature only changes the decay time of the normal fluid phase, as well as the crossover from superfluid to normal fluid. However, for the temperatures shown and lower, the normal phase is well outside the parameter space where perturbation theory is valid. (a) GC3 LC, T/T_{c} = 10^{−2}. (b) GC3 LC, T/T_{c} = 10^{−4}. (c) GC3 LC, T/T_{c} = 10^{−6}. (d) GC3 WC, T/T_{c} = 10^{−2}. (e) GC3 WC, T/T_{c} = 10^{−4}. (f) GC3 WC, T/T_{c} = 10^{−6}. (g) GC4 LC, T/T_{c} = 10^{−2}. (h) GC4 LC, T/T_{c} = 10^{−4}. (i) GC4 LC, T/T_{c} = 10^{−6}. (j) GC4 WC, T/T_{c} = 10^{−2}. (k) GC4 WC, T/T_{c} = 10^{−4}. (l) GC4 WC, T/T_{c} = 10^{−6}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.