© 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; Boylan-Kolchin 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 Maxwell-Boltzmann 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ánchez-Salcedo & Brandenburg 1999; Lee & Stahler 2011, 2014; Thun et al. 2016), relativistic (Barausse 2007; Katz et al. 2019), and magnetized gases (Sánchez-Salcedo 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, large-scale structure, the expansion history of the Universe, and important properties of galaxies (Davis et al. 1985; Percival et al. 2001; Tegmark et al. 2004; Trujillo-Gomez 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 & Boylan-Kolchin 2017). These discrepancies may have their solution within ΛCDM by including more realistic models of baryonic physics in simulations (Santos-Santos 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ía-Bellido 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; Bar-Or et al. 2019; Lancaster et al. 2020), and self-interacting Bose-Einstein condensed (SIBEC) DM (Berezhiani et al. 2019), also known as superfluid DM. A number of studies have considered finite-temperature 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 large-scale success of CDM. For the fifth force to be MONDian, the DM particles need exotic three-body self-interactions, and the DM fluid has to be above a certain temperature to be well-behaved. Finite-temperature 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 well-suited 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; Arca-Sedda & Capuzzo-Dolcetta 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 fine-tuning 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 so-called timing-problem, 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 zero-temperature 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 self-interactions, 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 SIBEC-DM, 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 finite-temperature 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 steady-state and a finite-time 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 SIBEC-DM are inferred. In Sect. 8 a summary of this work and the main results are presented. Natural units are used throughout.

2. Hydrodynamics of finite-temperature superfluids

In the standard treatment, superfluids are often related to Bose-Einstein 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 single-particle wave-function. This wave-function is usually associated with the superfluid, and can therefore be regarded as a quantum mechanical effect at macroscopic scales. The wave-function ψ at the mean-field level is governed by the Gross-Pitaevskii equation, a non-linear Schrödinger equation with effective contact interactions parameterized by g;

$$\begin{array}{c}\hfill i\frac{\partial \psi }{\partial t}=[\frac{-{\mathrm{\nabla }}^2}{2m}+g{|\psi |}^2+m{V}_{\mathrm{ext}}]\psi .\end{array}$$(1)

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 = |ψ|², and mass density ρ = m|ψ|².

By inserting for the wave-function

$$\begin{array}{c}\hfill \psi =\sqrt{n}{e}^{\mathrm{iS}}=\sqrt{\frac{\rho }{m}}{e}^{\mathrm{iS}},\end{array}$$(2)

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

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{v}}{\partial t}+(\mathit{v}·\mathbf{\nabla })\mathit{v}+\mathrm{\nabla }(\frac{g\rho }{m^2}+Q+V_{\mathrm{ext}})=0.\hfill \end{array}$$(4)

These are the so-called 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

$$\begin{array}{c}\hfill Q=-\frac{1}{2m^2}\frac{{\mathrm{\nabla }}^2\sqrt{\rho }}{\sqrt{\rho }},\end{array}$$(5)

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

$$\begin{array}{c}\hfill m{\displaystyle \oint }\mathit{v}·\mathrm{d}\mathit{l}=2\pi N,N\in \mathbb{Z},\end{array}$$(6)

because the complex wave-function must be single-valued. 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 small-scale structure that is predicted in N-body 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 two-fluid picture of superfluids. The hydrodynamic equations for a finite-temperature superfluid are (neglecting the quantum potential; Taylor & Griffin 2005; Chapman et al. 2014):

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·\mathit{j}=0,\hfill \end{array}$$(7)

$$\begin{array}{cc}& \frac{\partial S}{\partial t}+\mathbf{\nabla }·(S{\mathit{u}}_{\mathrm{n}})=0,\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{\partial {\mathit{u}}_{\mathrm{s}}}{\partial t}+\mathbf{\nabla }(\mu +\frac{1}{2}{\mathit{u}}_{\mathrm{s}}^2)=-\mathbf{\nabla }\mathrm{\Phi },\hfill \end{array}$$(9)

$$\begin{array}{cc}& \frac{\partial \mathit{j}}{\partial t}+\mathbf{\nabla }P+{\rho }_{\mathrm{s}}({\mathit{u}}_{\mathrm{s}}·\mathbf{\nabla }){\mathit{u}}_{\mathrm{s}}+{\rho }_{\mathrm{n}}({\mathit{u}}_{\mathrm{n}}·\mathbf{\nabla }){\mathit{u}}_{\mathrm{n}}\hfill \\ \hfill & +{\mathit{u}}_{\mathrm{s}}[\mathbf{\nabla }·({\rho }_{\mathrm{s}}{\mathit{u}}_{\mathrm{s}})]+{\mathit{u}}_{\mathrm{n}}[\mathbf{\nabla }·({\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}})]=-\rho \mathbf{\nabla }\mathrm{\Phi }.\hfill \end{array}$$(10)

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 = ρ_nu_n + ρ_su_s. The fluid pressure is P, the entropy density S, temperature T, and $\mu =[P+U-ST-\frac{1}{2}{\rho }_{\mathrm{n}}{({\mathit{u}}_{\mathrm{s}}-{\mathit{u}}_{\mathrm{n}})}^2]/\rho$.

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 co-called 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,

$$\begin{array}{c}\hfill v<v_{\mathrm{c}}=\underset{\mathit{p}}{min}\frac{ϵ(\mathit{p})}{p},\end{array}$$(11)

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²/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_sp. 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 self-interactions (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 finite-temperature 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;

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(12)

$$\begin{array}{cc}& \frac{\partial S}{\partial t}+\mathbf{\nabla }·(S\mathit{v})-\mathbf{\nabla }·(\frac{S{\rho }_{\mathrm{s}}}{\rho }\mathit{w})=0.\hfill \end{array}$$(13)

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 three-dimensional, homogeneous, ideal Bose gas, this critical temperature is

$$\begin{array}{c}\hfill T_{\mathrm{c}}=\frac{2\pi {ħ}^2}{m^{5/3}}{\left(\frac{\rho }{\zeta (3/2)}\right)}^{2/3},\end{array}$$(14)

where ζ(x) is the Riemann Zeta-function, 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 two-body interactions,

$$\begin{array}{cc}& P=\frac{1}{2}\frac{g}{m^2}{\rho }^2+\zeta (5/2){\left(\frac{m}{2\pi }\right)}^{3/2}{T}^{5/2},\hfill \end{array}$$(15)

$$\begin{array}{cc}& S=\frac{5}{2}\zeta (5/2){\left(\frac{m}{2\pi }\right)}^{3/2}{T}^{3/2},\hfill \end{array}$$(16)

valid only for T < T_c. The fraction of particles in the condensate f₀ and the superfluid f_s = ρ_s/ρ are both taken to be equal to the condensate fraction in the ideal case;

$$\begin{array}{c}\hfill f_{\mathrm{s}}=f_0=1-(\frac{T}{T_{\mathrm{c}}}{)}^{3/2}·\end{array}$$(17)

The critical velocity is approximated as

$$\begin{array}{c}\hfill v_{\mathrm{c}}=\sqrt{\frac{gn{f}_0}{m}}·\end{array}$$(18)

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 steady-state 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:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G[\rho +{\rho }_{\mathrm{pert}}].\end{array}$$(19)

The superfluid is assumed to be homogeneous, and so the fluid variables are expanded to linear order, ρ = ρ₀ + δρ, S = S₀ + δS, u_s = δu_s, and so on. The linear equations are

$$\begin{array}{cc}& \frac{\partial \delta \rho }{\partial t}+\mathbf{\nabla }·\delta \mathit{j}=0,\hfill \end{array}$$(20)

$$\begin{array}{cc}& \frac{\partial \delta S}{\partial t}+S_0\mathbf{\nabla }·\delta {\mathit{u}}_{\mathrm{n}}=0,\hfill \end{array}$$(21)

$$\begin{array}{cc}& \frac{\partial \delta {\mathit{u}}_{\mathrm{s}}}{\partial t}+\frac{1}{{\rho }_0}\mathbf{\nabla }\delta P-\frac{S_0}{{\rho }_0}\mathbf{\nabla }\delta T=-\mathbf{\nabla }\delta \mathrm{\Phi },\hfill \end{array}$$(22)

$$\begin{array}{cc}& \frac{\partial \delta \mathit{j}}{\partial t}+\mathbf{\nabla }\delta P=-{\rho }_0\mathbf{\nabla }\delta \mathrm{\Phi },\hfill \end{array}$$(23)

$$\begin{array}{cc}& \delta {\mathit{u}}_{\mathrm{n}}=\frac{1}{{\rho }_0}\delta \mathit{j}-\frac{{\rho }_{s0}}{{\rho }_0}\delta {\mathit{u}}_{\mathrm{s}},\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\mathrm{\nabla }}^2\delta \mathrm{\Phi }=4\pi G[\delta \rho +{\rho }_{\mathrm{pert}}].\hfill \end{array}$$(25)

These can be combined into two coupled equations for δρ and δS;

$$\begin{array}{cc}\hfill \frac{{\partial }^2\delta \rho }{\partial t^2}& -[{\left(\frac{\partial P}{\partial \rho }\right)}_0{\mathrm{\nabla }}^2+4\pi G{\rho }_0]\delta \rho -{\left(\frac{\partial P}{\partial S}\right)}_0{\mathrm{\nabla }}^2\delta S=4\pi G{\rho }_0{\rho }_{\mathrm{pert}},\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill \frac{{\partial }^2\delta S}{\partial t^2}& -\frac{S_0}{{\rho }_0}[{\left(\frac{\partial P}{\partial \rho }\right)}_0{\mathrm{\nabla }}^2+S_0\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial \rho }\right)}_0{\mathrm{\nabla }}^2+4\pi G{\rho }_0]\delta \rho \hfill \\ \hfill & -\frac{S_0}{{\rho }_0}[{\left(\frac{\partial P}{\partial S}\right)}_0{\mathrm{\nabla }}^2+S_0\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0{\mathrm{\nabla }}^2]\delta S=4\pi G{S}_0{\rho }_{\mathrm{pert}}.\hfill \end{array}$$(27)

The “0” subscript indicates that the quantity is evaluated at the background level. As expected, there are scale-dependent 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 δρ = αρ₀, and Fourier transforming into momentum (k) and frequency (k₀) space, the solutions of the k-modes α_k are found:

$$\begin{array}{c}\hfill {\alpha }_k=-4\pi G{\rho }_{\mathrm{pert},k}\frac{k_0^2-A{k}^2}{(k_0^2-{\omega }_{k+}^2)(k_0^2-{\omega }_{k-}^2)},\end{array}$$(28)

where the dispersion relation is

$$\begin{array}{c}\hfill {\omega }_{k\pm }^2=C_4{k}^2-C_2\pm \sqrt{C_3{k}^4-2C_1{C}_2{k}^2+C_2^2},\end{array}$$(29)

and

$$\begin{array}{cc}& A=\frac{S_0^2}{{\rho }_0}\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0,\hfill \end{array}$$(30)

$$\begin{array}{cc}& C_1=\frac{1}{2}{\left(\frac{\partial P}{\partial \rho }\right)}_0+\frac{S_0}{2{\rho }_0}{\left(\frac{\partial P}{\partial S}\right)}_0-\frac{S_0^2}{2{\rho }_0}\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0,\hfill \end{array}$$(31)

$$\begin{array}{cc}& C_2=2\pi G{\rho }_0,\hfill \end{array}$$(32)

$$\begin{array}{cc}& C_3=C_4^2+\frac{S_0^2}{{\rho }_0}\frac{{\rho }_{s0}}{{\rho }_{n0}}[{\left(\frac{\partial P}{\partial S}\right)}_0{\left(\frac{\partial T}{\partial \rho }\right)}_0-{\left(\frac{\partial P}{\partial \rho }\right)}_0{\left(\frac{\partial T}{\partial S}\right)}_0],\hfill \end{array}$$(33)

$$\begin{array}{cc}& C_4=\frac{1}{2}{\left(\frac{\partial P}{\partial \rho }\right)}_0+\frac{S_0}{2{\rho }_0}{\left(\frac{\partial P}{\partial S}\right)}_0+\frac{S_0^2}{2{\rho }_0}\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0·\hfill \end{array}$$(34)

The dynamical friction is given by the change in the energy of the perturber,

$$\begin{array}{c}\hfill F_{\mathrm{DF}}=-\frac{M}{V}\frac{\partial {\mathrm{\Phi }}_{\alpha }}{\partial t},\end{array}$$(35)

where M and V are the mass and velocity of the perturber, and Φ_α is the gravitational potential of the background fluid,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2{\mathrm{\Phi }}_{\alpha }=4\pi G{\rho }_0\alpha .\end{array}$$(36)

This is readily found in k-space,

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\alpha ,k}=-\frac{4\pi G{\rho }_0{\alpha }_k}{k^2},\end{array}$$(37)

which can be Fourier transformed back into position-space to give the dynamical friction,

$$\begin{array}{cc}\hfill F_{\mathrm{DF}}& =\frac{M}{V}\frac{\partial }{\partial t}{\displaystyle \int }\frac{\mathrm{d}{k}^4}{{(2\pi )}^4}{e}^{i{k}_{0}t-i\mathit{k}·\mathit{x}}{\mathrm{\Phi }}_{\alpha ,k}\hfill \\ \hfill & =-\frac{4\pi G{M}^2{\rho }_0}{V}{\displaystyle \int }\frac{\mathrm{d}{k}^4}{{(2\pi )}^4}\frac{i{k}_0}{k^2}{e}^{i{k}_{0}t-i\mathit{k}·\mathit{x}}{\alpha }_k.\hfill \end{array}$$(38)

Approximating the perturber as a point particle moving along the z-axis with constant velocity V,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{pert}}(\mathit{x},t)=M\delta (x)\delta (y)\delta (z-Vt),\end{array}$$(39)

or in k-space

$$\begin{array}{c}\hfill {\rho }_{\mathrm{pert},k}=2\pi M\delta (k_0-V{k}_z),\end{array}$$(40)

yields the expression for the dynamical friction as

$$\begin{array}{cc}\hfill F_{\mathrm{DF}}=& \frac{32{\pi }^3{G}^2{M}^2\rho }{V}{\displaystyle \int }\frac{\mathrm{d}{k}^4}{{(2\pi )}^4}\frac{i{k}_0}{k^2}{e}^{i{k}_{0}t-i\mathit{k}·\mathit{x}}\hfill \\ \hfill & \times \frac{(k_0^2-A{k}^2)\delta (k_0-V{k}_z)}{(k_0^2-{\omega }_{k+}^2)(k_0^2-{\omega }_{k-}^2)}·\hfill \end{array}$$(41)

Equation (41) can be tackled by extending the k₀-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

$$\begin{array}{cc}\hfill F_{\mathrm{DF}}=& -\frac{16{\pi }^3{G}^2{M}^2{\rho }_0}{V}{\displaystyle \int }\frac{\mathrm{d}{k}^3}{{(2\pi )}^3}\frac{1}{k^2}\frac{1}{{\omega }_{k+}^2-{\omega }_{k-}^2}\hfill \\ \hfill & \times [e^{i{\omega }_{k+}t-i\mathit{k}·\mathit{x}}({\omega }_{k+}^2-A{k}^2)\delta ({\omega }_{k+}-V{k}_z)\hfill \\ \hfill & -e^{i{\omega }_{k-}t-i\mathit{k}·\mathit{x}}({\omega }_{k-}^2-A{k}^2)\delta ({\omega }_{k-}-V{k}_z)].\hfill \end{array}$$(42)

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 z-axis, and the force is evaluated at the position of the perturber, $\mathit{x}=Vt\widehat{\mathit{z}}$. 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, $k<k_{\max }^{\ast ,\pm }$, where $k_{\max }^{\ast ,\pm }$ satisfies kV = ω_k±. Further constraints are placed on k: The remaining k-integral 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 $k>k_{\min }^{\ast ,\pm }$, where $k_{\min }^{\ast ,\pm }$ 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 $k_{\max }^{\pm }$, and the lower limits by $k_{\min }^{\pm }$. Inserting the expression for ω_k± and using that C₄ − A = C₁, the dynamical friction becomes

$$\begin{array}{cc}\hfill F_{\mathrm{DF}}=& -\frac{4\pi G^2{M}^2{\rho }_0}{V^2}[{\displaystyle {\int }_{k_{\min }^{+}}^{k_{\max }^{+}}\frac{\mathrm{d}k}{2k}(\frac{C_1{k}^2-C_2}{\sqrt{C_3{k}^4-2C_1{C}_2{k}^2+C_2^2}}+1)}\hfill \\ \hfill & -{\displaystyle {\int }_{k_{\min }^{-}}^{k_{\max }^{-}}\frac{\mathrm{d}k}{2k}(\frac{C_1{k}^2-C_2}{\sqrt{C_3{k}^4-2C_1{C}_2{k}^2+C_2^2}}-1)]·}\hfill \end{array}$$(43)

There is an implicit criterion that $k_{\max }^{\pm }>k_{\min }^{\pm }>0$, 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 self-gravitation.

3.1. Zero-temperature limit

Taking the limit T → 0 (under the assumption that terms such as S²ρ_s/ρ_n go to zero as well) yields one band for the dispersion relation,

$$\begin{array}{c}\hfill {\omega }_k^2=c_{T=0}^2{k}^2-4\pi G{\rho }_0,\end{array}$$(44)

where the sound speed at zero temperature is

$$\begin{array}{c}\hfill c_{T=0}^2={\left(\frac{\partial P}{\partial \rho }\right)}_0·\end{array}$$(45)

The dynamical friction becomes

$$\begin{array}{c}\hfill F_{\mathrm{DF}}=-\frac{4\pi G^2{M}^2{\rho }_0}{V^2}ln(\frac{k_{\max }}{k_{\min }}),\end{array}$$(46)

with

$$\begin{array}{cc}& k_{\max }=min(2\pi R_{\min }^{-1},\sqrt{\frac{4\pi G{\rho }_0}{c_{T=0}^2-V^2}}),\hfill \end{array}$$(47)

$$\begin{array}{cc}& k_{\min }=max(2\pi R_{\max }^{-1},\sqrt{\frac{4\pi G{\rho }_0}{c_{T=0}^2}})·\hfill \end{array}$$(48)

3.2. Normal fluid limit

Taking the fully normal fluid limit ρ_s → 0 also gives one band for the dispersion relation,

$$\begin{array}{c}\hfill {\omega }_k^2=c_{\mathrm{n}}^2{k}^2-4\pi G{\rho }_0,\end{array}$$(49)

with the sound speed in the fully normal fluid

$$\begin{array}{c}\hfill c_{\mathrm{n}}^2={\left(\frac{\partial P}{\partial \rho }\right)}_0+\frac{S_0}{{\rho }_0}{\left(\frac{\partial P}{\partial S}\right)}_0·\end{array}$$(50)

The dynamical friction is again given by Eq. (46), but with

$$\begin{array}{cc}& k_{\max }=min(2\pi R_{\max }^{-1},\sqrt{\frac{4\pi G{\rho }_0}{c_{\mathrm{n}}^2-V^2}}),\hfill \end{array}$$(51)

$$\begin{array}{cc}& k_{\min }=max(2\pi R_{\min }^{-1},\sqrt{\frac{4\pi G{\rho }_0}{c_{\mathrm{n}}^2}})·\hfill \end{array}$$(52)

This is the same as the zero-temperature case, but with a different sound speed.

3.3. Small-velocity limit

At sufficiently small velocities, $V^2\ll C_4-\sqrt{C_3}$, assuming that the finite sizes of the background cloud and perturber do not set the integral limits in Eq. (43), the dynamical friction becomes

$$\begin{array}{c}\hfill F_{\mathrm{DF}}=-\frac{2\pi G^2{M}^2{\rho }_0}{c_{T=0}^2}·\end{array}$$(53)

This is equal to the friction force at T = 0 in the same limit, as opposed to when ρ_s = 0;

$$\begin{array}{c}\hfill F_{\mathrm{DF}}=-\frac{2\pi G^2{M}^2{\rho }_0}{c_{\mathrm{n}}^2}·\end{array}$$(54)

The dynamical friction of a superfluid therefore approaches the zero-temperature limit even when there is a significant thermal contribution. This happens because counterflow in the superfluid conspires against thermal perturbations, allowing the mass over-density to behave similarly to a zero-temperature 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 self-gravitation

The numerical results presented in Sect. 6, as well as the decay times of globular clusters in Sect. 7, are obtained when self-gravitation is neglected. It is therefore of interest to see what the steady-state linear theory predicts in this case as well.

Neglecting self-gravitation amounts to setting C₂ = 0. The dispersion relation becomes

$$\begin{array}{c}\hfill {\omega }_{k\pm }^2=(C_4\pm \sqrt{C_3})k^2=c_{\pm }^2{k}^2.\end{array}$$(55)

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

$$\begin{array}{c}\hfill c_{+}=\sqrt{\frac{c_{\mathrm{n}}^2-c_{T=0}^2}{f_{\mathrm{n}}}},\end{array}$$(56)

$$\begin{array}{c}\hfill c_{-}=c_{T=0}.\end{array}$$(57)

We note that for c_n ≫ c_T = 0, we have $c_{+}\approx c_{\mathrm{n}}/\sqrt{f_{\mathrm{n}}}\gg c_{\mathrm{n}}$. The dynamical friction takes the form

$$\begin{array}{cc}\hfill F_{\mathrm{DF}}=& -\frac{4\pi G^2{M}^2{\rho }_0}{V^2}ln(\frac{R_{\max }}{R_{\min }})\hfill \\ \hfill & \times \frac{1}{2}[(1-\frac{C_1}{\sqrt{C_3}})\theta (V-c_{-})+(1+\frac{C_1}{\sqrt{C_3}})\theta (V-c_{+})].\hfill \end{array}$$(58)

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 steady-state case, as considered in this section, the linear over-density 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 first-order perturbation of the background becomes symmetric. In order to overcome this shortcoming of steady-state linear perturbation theory, other studies have broken this symmetry by switching on the perturber for a finite time (Ostriker 1999; Sánchez-Salcedo 2012), or by going to second-order perturbations (Lee & Stahler 2011; Shadmehri & Khajenabi 2012). In the following section, the finite-time approach is employed for a superfluid.

4. Dynamical friction from finite-time linear perturbation theory

For the finite-time calculation, Eqs. (26) and (27) are used without self-gravitation, and an approach very similar to the one used by Ostriker (1999) is followed.

The equations can be written in matrix form as

$$\begin{array}{c}\hfill \frac{{\partial }^{2}Y}{\partial t^2}+A{\mathrm{\nabla }}^{2}Y=F{\rho }_{\mathrm{pert}},\end{array}$$(59)

where

$$\begin{array}{cc}& Y=\left(\begin{array}{c}\delta \rho \\ \delta S\end{array}\right),\hfill \end{array}$$(60)

$$\begin{array}{cc}& A=\left(\begin{array}{cc}{\left(\frac{\partial P}{\partial \rho }\right)}_0& {\left(\frac{\partial P}{\partial S}\right)}_0\\ \frac{S_0}{{\rho }_0}{\left(\frac{\partial P}{\partial \rho }\right)}_0+\frac{S_0^2}{{\rho }_0}\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial \rho }\right)}_0& \frac{S_0}{{\rho }_0}{\left(\frac{\partial P}{\partial S}\right)}_0+\frac{S_0^2}{{\rho }_0}\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0\end{array}\right),\hfill \end{array}$$(61)

and

$$\begin{array}{c}\hfill F=\left(\begin{array}{c}4\pi G{\rho }_0\\ 4\pi G{S}_0\end{array}\right).\end{array}$$(62)

By diagonalizing matrix A, the coupled set of equations can be transformed into two decoupled wave equations of the form

$$\begin{array}{c}\hfill \frac{{\partial }^2{\chi }_i}{\partial t^2}-c_i^2{\mathrm{\nabla }}^2{\chi }_i=f_i,\end{array}$$(63)

which are solved using the retarded Green’s function for the wave equation in three dimensions:

$$\begin{array}{c}\hfill {\chi }_i(\mathit{x},t)={\displaystyle \int {\mathrm{d}}^{3}x\prime \int \mathrm{d}t\prime \frac{\delta (t\prime -(t-|\mathit{x}-\mathit{x}\prime |/c_i))f_i(\mathit{x}\prime ,t\prime )}{4\pi c_i^2|\mathit{x}-\mathit{x}\prime |}·}\end{array}$$(64)

For a point source switched on at the origin at t = 0 and moving at speed $\mathit{V}=V\widehat{\mathit{z}}$,

$$\begin{array}{c}\hfill f_i(\mathit{x},t)=K_i\delta (x)\delta (y)\delta (z-Vt)H(t),\end{array}$$(65)

where H(x) is the Heaviside function, the solution of χ becomes, upon defining s = z − Vt, ℳ_i = V/c_i, and R² = x² + y²,

$$\begin{array}{cc}& {\chi }_i(\mathit{x},t)=\frac{K_i}{4\pi c_i^2\sqrt{s^2+R^2(1-{\mathcal{M}}_i^2)}}{\mathcal{H}}_i,\hfill \end{array}$$(66)

$$\begin{array}{cc}& {\mathcal{H}}_i=\{\begin{array}{cc}1\hfill & \mathrm{for}{R}^2+z^2<{(c_{i}t)}^2,\hfill \\ 2\hfill & \mathrm{for}{\mathcal{M}}_i>1,R^2+z^2>{(c_{i}t)}^2,\hfill \\ \hfill & s/R<-\sqrt{{\mathcal{M}}_i^2-1},\mathrm{and}z>c_{i}t/{\mathcal{M}}_i,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\hfill \end{array}$$(67)

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,

$$\begin{array}{c}\hfill F_{\mathrm{DF}}=2\pi GM{\displaystyle \int \mathrm{d}s\int \mathrm{d}R\frac{sR\delta \rho }{{(s^2+R^2)}^{3/2}}·}\end{array}$$(68)

In spherical polar coordinates, s = r cos θ = rx and $R=rsin\theta =r\sqrt{1-x^2}$, we get

$$\begin{array}{cc}& F_{\mathrm{DF}}=-\frac{4\pi G^2{M}^2{\rho }_0}{V^2}({\mathcal{I}}_{+}+{\mathcal{I}}_{-}),\hfill \end{array}$$(69)

$$\begin{array}{cc}& {\mathcal{I}}_i=-D_i{\displaystyle {\int }_{R_{\min }}^{R_{\max }}\frac{\mathrm{d}r}{2r}{\int }_{-1}^1\mathrm{d}x\frac{x{\mathcal{M}}_i^2{\mathcal{H}}_i}{\sqrt{1-{\mathcal{M}}_i^2+x^2{\mathcal{M}}_i^2}},}\hfill \end{array}$$(70)

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

$$\begin{array}{cc}& D_{+}=-\frac{S_0{\left(\frac{\partial P}{\partial S}\right)}_0[S_0\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial \rho }\right)}_0+c_{+}^2]}{{\rho }_0(c_{+}^2-c_{-}^2)[{\left(\frac{\partial P}{\partial \rho }\right)}_0-c_{+}^2]},\hfill \end{array}$$(71)

$$\begin{array}{cc}& D_{-}=\frac{S_0{\left(\frac{\partial P}{\partial S}\right)}_0[S_0\frac{{\rho }_{s0}}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial \rho }\right)}_0+c_{-}^2]}{{\rho }_0(c_{+}^2-c_{-}^2)[{\left(\frac{\partial P}{\partial \rho }\right)}_0-c_{-}^2]}·\hfill \end{array}$$(72)

The dynamical friction from the finite-time calculation is compared to the steady-state 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² dependence as in the steady-state result. As time passes, the finite-time result approaches the steady-state result, as expected.

Fig. 1. Dynamical friction from linear perturbation theory using the finite-time approach (solid lines) and the steady-state approach (dotted lines) as a function of V. As time passes, the finite-time result, Eq. (69), approaches the steady-state result, Eq. (58). In the zero-temperature limit, we have T = 0, while in the normal fluid case we have ρs = 0. (a) t = 0.1Rmax/V. (b) t = Rmax/V. (c) t = 10Rmax/V.

Both approaches predict F_DF in the superfluid phase to be very close to the zero-temperature value, even when thermal pressure dominates over the contribution from self-interactions. 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,

$$\begin{array}{c}\hfill \frac{\partial \mathit{w}}{\partial t}=\frac{S_0}{{\rho }_{n0}}\mathbf{\nabla }\delta T=\frac{S_0}{{\rho }_{n0}}[{\left(\frac{\partial T}{\partial S}\right)}_0\mathbf{\nabla }\delta S+{\left(\frac{\partial T}{\partial \rho }\right)}_0\mathbf{\nabla }\delta \rho ]·\end{array}$$(73)

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 $\overline{\mathit{w}}$ of the system, there is an interpolating function $f(\overline{\mathit{w}},{\mathit{v}}_{\mathrm{c}})$ with $f(\overline{\mathit{w}}\ll {\mathit{v}}_{\mathrm{c}})\to 1$, $f(\overline{\mathit{w}}\gg {\mathit{v}}_{\mathrm{c}})\to 0$, and a transitional region around $\overline{\mathit{w}}\sim {\mathit{v}}_{\mathrm{c}}$, such that

$$\begin{array}{c}\hfill F_{\mathrm{DF}}=f(\overline{w},v_{\mathrm{c}})F_{\mathrm{DF}}^{\mathrm{sf}}+[1-f(\overline{w},v_{\mathrm{c}})]F_{\mathrm{DF}}^{\mathrm{nf}},\end{array}$$(74)

where $F_{\mathrm{DF}}^{\mathrm{sf}}$ and $F_{\mathrm{DF}}^{\mathrm{nf}}$ is the dynamical friction from linear theory for the superfluid and fully normal fluid, respectively. Using Eq. (73) we can write

$$\begin{array}{c}\hfill \overline{w}=\frac{S_0}{{\rho }_{n0}}\frac{\delta T(0)}{L}\mathrm{\Delta }t=\frac{S_0}{{\rho }_{n0}}[{\left(\frac{\partial T}{\partial S}\right)}_0\frac{\delta S(0)}{L}+{\left(\frac{\partial T}{\partial \rho }\right)}_0\frac{\delta \rho (0)}{L}]\mathrm{\Delta }t.\end{array}$$(75)

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 well-defined. 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

$$\begin{array}{c}\hfill \delta S(R_{\min }/2)\approx \frac{2S_{0}GM}{c_{+}^2{R}_{\min }}·\end{array}$$(76)

The rough estimate of the counterflow is therefore

$$\begin{array}{c}\hfill \overline{w}=\frac{S_0^2}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0\frac{2GM}{c_{+}^3{R}_{\min }}·\end{array}$$(77)

Only the form of the interpolating function $f(\overline{\mathit{w}},{\mathit{v}}_{\mathrm{c}})$ remains to be specified. The simple but rather arbitrary choice,

$$\begin{array}{c}\hfill f(\overline{w},v_{\mathrm{c}})=\frac{v_{\mathrm{c}}}{v_{\mathrm{c}}+\overline{w}}={[1+\frac{S_0^2}{{\rho }_{n0}}{\left(\frac{\partial T}{\partial S}\right)}_0\frac{2G}{c_{+}^3{R}_{\min }}\frac{M}{v_{\mathrm{c}}}]}^{-1},\end{array}$$(78)

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 ${\rho }_{\text{pert}}=3M/4\pi R_{\text{pert}}^3$. 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 $\mathit{V}=-V\widehat{z}$. 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 self-gravitation 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⁷ M_⊙ kpc⁻³, T = 0.2T_c, m = 500 eV, and g = 2 × 10⁻³ eV⁻². These parameters are chosen only to illustrate the basic features of dynamical friction in superfluids while keeping the simulation run-times 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 zero-temperature and normal fluid limits, for which the numerical scheme was made second-order in time and space using a MinMod slope-limited MUSCL-Hancock scheme (Toro 2006) without stability issues, a lower resolution of 2048 × 1024 is used.

An effective critical velocity ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$, 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 self-interaction.

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 zero-temperature 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 well-defined supersonic cone that trails the pertuber, and the density contrast reaches about 4.5. The finite-temperature superfluid has a similar density contrast and supersonic cone, though not as well-defined, 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 ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$ 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 ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$, and for the zero temperature and the fully normal fluid limits. The results from finite-time 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 ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$ is decreased, the superfluid approaches the fully normal fluid limit as thermal counterflow is increasingly limited. The sound speeds cT = 0 and cn are indicated by the vertical dotted lines.

Fig. 3. Density profiles and streamlines for V = 1.5cT = 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 = vs − vn. 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) ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}={\mathit{v}}_{\mathrm{c}}$. (d) ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}={10}^{-1}{\mathit{v}}_{\mathrm{c}}$. (e) ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}={10}^{-5}{\mathit{v}}_{\mathrm{c}}$.

Fig. 4. Dynamical friction against the effective critical velocity ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$, for V = 1.5cT = 0, with the results from finite-time linear perturbation theory included with dotted lines of the same colors. As ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$ 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 $\overline{\mathit{w}}\ll {\mathit{v}}_{\mathrm{c}}$, 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.5cT = 0, with the results from finite-time linear perturbation theory included with dotted lines of the same colors. The departure from perturbation theory for the zero-temperature 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 ${\mathit{v}}_{\mathrm{c}}^{\mathrm{eff}}$.

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, self-interaction, 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 SIBEC-DM, 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;

$$\begin{array}{c}\hfill {\tau }_{\mathrm{DF}}=\frac{L}{r|F_{\mathrm{DF}}|}=\frac{\mathit{MV}}{|F_{\mathrm{DF}}|},\end{array}$$(79)

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 double-power-law profile¹ of the form

$$\begin{array}{c}\hfill \rho (r)=\overline{\rho }{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{-{\gamma }_0}{[1+{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{\eta }]}^{\frac{{\gamma }_0-{\gamma }_{\infty }}{\eta }}.\end{array}$$(80)

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 SIBEC-DM, 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 ρ₀ = ρ(r) in the computation of the SIBEC-DM 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”.

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), $r=2r_{\perp }/\sqrt{3}$ 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,

$$\begin{array}{c}\hfill V=\sqrt{\frac{G{M}_{\mathrm{encl}}}{r}}·\end{array}$$(81)

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 self-interaction 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 $\overline{\mathit{w}}\sim {\mathit{v}}_{\mathrm{c}}$, 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

$$\begin{array}{c}\hfill g\lesssim \pi G{R}_{\mathrm{c}}^2{m}^2.\end{array}$$(82)

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

$$\begin{array}{c}\hfill g\gtrsim \frac{m^2}{9\rho }\frac{R_{\mathrm{c}}^2}{\mathrm{\Delta }{t}^2},\end{array}$$(83)

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 zero-temperature estimate used to derive Eq. (82).

Fig. 7. Criteria listed in the text, and the orbital decay timescale for GC3 at T/Tc = 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). $\overline{\mathit{w}}={\mathit{v}}_{\mathrm{c}}$, 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 self-interaction 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⁻² g⁻¹. In terms of the self-interaction parameter g, this constraint reads (Pitaevskii & Stringari 2016)

$$\begin{array}{c}\hfill g=\frac{\sqrt{4\pi \sigma }}{m}\lesssim 5\times {10}^{-12}{\left(\frac{1\mathrm{eV}}{m}\right)}^{1/2}{\mathrm{eV}}^{-2}.\end{array}$$(84)

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, ${\mathrm{\Gamma }}_{\mathrm{grav}}\sim \sqrt{G\rho }$. For two-body interactions, the relaxation rate is Γ ∼ nσδv, where σ is the scattering cross-section and δv the velocity dispersion of the particles. In terms of g, as above, the cross-section is σ = m²g²/4π. However, for a condensed Bose gas, the relaxation rate is enhanced, that is, Γ ∼ 𝒩nσδv, where

$$\begin{array}{c}\hfill \mathcal{N}=n\frac{{(2\pi )}^3}{\frac{4\pi }{3}{(m\delta v)}^3},\end{array}$$(85)

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

$$\begin{array}{c}\hfill g\gtrsim \sqrt{\frac{2}{3\pi }}\frac{m^{3/2}{G}^{1/4}V}{{\rho }^{3/4}}·\end{array}$$(86)

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 finite-time 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 finite-time formalism is the lack of self-gravitation in the background fluid, and it might seem that a better choice is to instead use the steady-state 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 large-scale 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 self-gravitation 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 timing-problem usually applies too (Cole et al. 2012; Hui et al. 2017; Arca-Sedda & Capuzzo-Dolcetta 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, $\overline{\mathit{w}}={\mathit{v}}_{\mathrm{c}}$, lies far away from the region δρ/ρ₀ < 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/Tc = 10−2. (b) GC3 LC, T/Tc = 10−4. (c) GC3 LC, T/Tc = 10−6. (d) GC3 WC, T/Tc = 10−2. (e) GC3 WC, T/Tc = 10−4. (f) GC3 WC, T/Tc = 10−6. (g) GC4 LC, T/Tc = 10−2. (h) GC4 LC, T/Tc = 10−4. (i) GC4 LC, T/Tc = 10−6. (j) GC4 WC, T/Tc = 10−2. (k) GC4 WC, T/Tc = 10−4. (l) GC4 WC, T/Tc = 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 δρ/ρ₀ < 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)

$$\begin{array}{c}\hfill F_{\mathrm{DF},\mathrm{CDM}}=-\frac{4\pi M^2{G}^2{\rho }_{0}ln\mathrm{\Lambda }}{V^2}[\mathrm{erf}(X)-\frac{2X}{\sqrt{\pi }}{e}^{-X^2}],\end{array}$$(87)

where Λ ≈ rδv²/GM, $X=V/\sqrt{2}\delta \mathit{v}$, δ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 SIBEC-DM 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 SIBEC-DM model considered here, GC3 and GC4 are currently racing towards the center of their host halo in a SIBEC-DM universe.

Let us now consider τ_DF in light of constraints on the SIBEC-DM model from the literature. By fitting rotation curves of slowly rotating SIBEC-DM 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 SIBEC-DM halos at T = 0, and found the best fit values for g/m² to be between 2.7 × 10⁻⁴ eV⁻⁴ and 5.0 × 10⁻² eV⁻⁴. For reference, the estimated limit from hydrostatic equilibrium using Eq. (82) gives g/m² of less than about 2 × 10⁻⁴ eV⁻⁴ or 10⁻³ eV⁻⁴, depending on the profile used for the dSph. As the preferred values obtained by Crǎciun & Harko (2020) for zero-temperature SIBEC-DM satisfy V < c₋ for the GCs and dSph profiles considered, leading to a vanishing dynamical friction, the T = 0 case does not have a timing-problem, 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 $R_{\mathrm{J}}\sim c_{\mathrm{s}}/\sqrt{G\rho }$ 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 finite-temperature SIBEC-DM halo–for which we expect the preferred values for g/m² 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 steady-state linear perturbation theory. The well-known issue of discontinuities in the friction force as the perturber’s velocity crosses the fluid sound speed was encountered. We therefore also employed a finite-time 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⁻² dependence as the steady-state 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 superfluid-normal 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 SIBEC-DM model for the values of g/m² obtained by Crǎciun & Harko (2020) by fitting rotation curves at T = 0. For a finite-temperature SIBEC-DM, on the other hand, for which the preferred parameter space of g/m² 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 time-stepping 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.

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