© ESO 2020

1. Introduction

A universe with cold dark matter (CDM), a cosmological constant (Λ), and inflationary initial conditions forms the foundation of the standard ΛCDM paradigm that has proven successful at explaining a wide range of observables, such as the expansion history of the universe, the cosmic microwave background, formation of large-scale structure, the matter power spectra, and the abundance of light elements (Tegmark et al. 2004; Planck Collaboration XIII 2016; Cyburt et al. 2016). Nonetheless, it is a phenomenological model that is ignorant of the origin of the cosmological constant and the identity of dark matter (DM), which remain two of the greatest mysteries in fundamental physics today.

A number of challenges to ΛCDM have emerged as both observations and numerical simulations become increasingly more precise, especially on small scales. The cores of DM halos predicted from N-body simulations are denser and more cuspy than observed, and the number of dwarf galaxies in the Local Group is far smaller than expected from pure ΛCDM simulations. These issues are known as the too-big-to-fail, cusp-core, and missing satellite problems (see e.g., Del Popolo & Le Delliou 2017; Bullock & Boylan-Kolchin 2017 and references therein). Another puzzling phenomenology on the scale of galaxies is the empirical baryonic Tully-Fisher relation (BTFR; McGaugh et al. 2000; McGaugh 2005; Lelli et al. 2015). This relates the baryonic mass of galaxies M_b with the asymptotic circular velocity v_c through $M_{\text{b}}\sim v_{\text{c}}^4$ and holds for many orders of magnitude with remarkably small scatter. The ΛCDM prediction for this relation is $M\sim v_{\text{c}}^3$ with the total mass M from both baryons and DM (McGaugh 2012). It is the latter that dominates the gravitational pull in galaxies, which only adds to the strangeness of the BTFR.

Solutions to these problems within the framework of ΛCDM have been proposed by including baryonic physics (Santos-Santos et al. 2015; Sales et al. 2016; Zhu et al. 2016; Sawala et al. 2016), but it is unclear if they can completely cure the ails of ΛCDM. These processes are not yet fully understood and are difficult to model in simulations of galaxy formation, and their stochastic nature makes it even more puzzling as to how they can be responsible for the tight correlation in scaling relations, such as the BTFR.

An alternative possibility is that the mismatch between observations and simulations is an indication of physics beyond the standard model, either through modified theories of gravity, the particle nature of DM, or both. An example of such a model is modified Newtonian dynamics (MOND; Milgrom 1983a,b,c; Famaey & McGaugh 2012), in which the Newtonian law of gravity in low-acceleration regions is modified to explain the rotation curves of galaxies without the need of resorting to DM. One of its most appealing features is that the BTFR and its small scatter is a direct consequence of it. However, MOND and its relativistic extensions face challenges of their own on extragalactic scales where the CDM paradigm is successful (Zuntz et al. 2010; Dodelson 2011; Angus et al. 2013, 2014). This has, somewhat ironically, motivated extended models of DM where MOND is an emergent fifth force on small scales (Berezhiani & Khoury 2015; Khoury 2016). This is achieved by DM undergoing Bose–Einstein condensation on galactic scales and adding a coupling designed to give a MONDian long-range force between baryons mediated by phonons in the superfluid cores of galaxies. Outside galaxies, the DM fluid ceases to be superfluid, and the extra force disappears, preserving the success of CDM on large scales.

Superfluid dark matter (SFDM) models are also interesting on more general grounds. From condensed matter physics, we know that self-interacting boson gases can become superfluid given sufficiently high densities and low temperatures. In the weakly interacting Bose gas, the critical temperature that marks the onset of superfluidity depends almost solely on the particle mass and number density. We can therefore expect boson DM candidates with self-interactions to exhibit superfluid behavior in certain mass ranges.

Observations of the large-scale structure of the universe strongly favor cold and collisionless DM, but for SFDM this is no longer the case since the transition in and out of the superfluid phase requires both self-interactions and finite temperatures. We must therefore be wary of how structure forms in SFDM. Studies of other DM models with pressure-like terms, such as fuzzy dark matter (Hu et al. 2000; Schive et al. 2014; Schwabe et al. 2016; Mocz et al. 2017) and self-interacting dark matter (Spergel & Steinhardt 2000; Elbert et al. 2015; Tulin & Yu 2018), find they can help remove the surplus of small-scale structure in ΛCDM. So far, there has been little work done on structure formation in SFDM and how it differs from conventional DM fluids. In this paper, we aim to provide preliminary answers to these questions by considering the spherical collapse of SFDM.

The paper is organized as follows: in Sect. 2, the equations for superfluid hydrodynamics used to describe the collapse of SFDM are introduced, as well as the critical temperature and the critical velocity, which are important for the superfluid phenomenology. The linear expansion of the superfluid equations was derived to better understand how superfluidity changes the behavior of the DM fluid. In Sect. 3, the results are presented and discussed, and we draw our conclusions in Sect. 4.

2. Method

2.1. Superfluid hydrodynamics

To describe a finite-temperature superfluid, we employed the superfluid hydrodynamic equations (Taylor & Griffin 2005; Chapman et al. 2014), which in proper coordinates and physical variables are;

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

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

$$\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}$$(3)

$$\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}$$(4)

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

This set of equations describes the evolution of the fluid mass density ρ, entropy density S, superfluid velocity u_s, momentum density j, and energy density E under the influence of the gravitational potential Φ sourced by matter and a cosmological constant,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G(\rho -2{\rho }_{\mathrm{\Lambda }}).\end{array}$$(6)

Equations (2) and (5) are degenerate in our set of equations if the solution is free of shocks, otherwise entropy is generated. The former is used in this work, but both are given for completeness.

A superfluid differs from a classical fluid in that it consists of two fluid components; the “superfluid” with density ρ_s and velocity u_s, and the “normal fluid” with density ρ_n and velocity u_n. The sum of the two component densities gives the total fluid density ρ = ρ_n + ρ_s, and likewise for momentum, j = ρ_nu_n + ρ_su_s. However, only the normal fluid transports entropy and thermal energy, as can be seen from Eqs. (2) and (5), and the superfluid velocity evolves according to its own potential given in Eq. (3), where the chemical potential is $\mu =[P+U-ST-\frac{1}{2}{\rho }_n{({\mathit{u}}_{\mathrm{s}}-{\mathit{u}}_{\mathrm{n}})}^2]/\rho$. Since there are two fluid components with separate velocity fields a superfluid can have two sound modes. One is called first sound and is associated with density perturbations, which we are familiar with from regular hydrodynamics. The other is called second sound and is associated with temperature perturbations. This is made possible by the fact that only the normal component carries entropy, hence the normal and superfluid components can oscillate in such a way that perturbations in temperature, and not density, are propagated through the fluid. As we will see it is this property that is responsible for the difference in collapse of superfluid and non-superfluid DM.

The remaining variables in the above set of equations are pressure P, internal energy density U, and temperature T. In the limit ρ_s = 0, they reduce to the Euler equations of fluid dynamics.

2.2. Critical temperature and velocity, and equation of state

When a boson gas is cooled below a critical temperature T_c, the particles begin accumulating in the quantum ground state of the system and form a Bose–Einstein condensate (BEC). In the three-dimensional homogeneous and ideal Bose gas this critical temperature is

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

where ζ(x) is the Riemann Zeta-function. This result holds approximately for weakly interacting gases as well (Sharma et al. 2019), apart from a small interaction-dependent shift (Andersen 2004) that we neglect.

The formation of a BEC does not automatically imply a superfluid. A further criterion must be satisfied as realized by Landau (1941). He assumed that if dissipation and heating happens through the creation of elementary excitations in the fluid, and if these excitations can no longer spontaneously appear the fluid will become superfluid. This gives the so-called Landau criterion and requires the relative motion w = u_s − u_n to be smaller than the critical velocity v_c,

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

where ϵ(p) is the energy of an elementary excitation with momentum p. Clearly, we must have v_c >  0, otherwise any motion will destroy the superflow, and it no longer makes sense to refer to it as a superfluid. An ideal Bose gas can therefore not be superfluid since the elementary excitations are ϵ(p) = p²/2m so that v_c = 0. In an interacting Bose gas, on the other hand, the condensation of the gas makes the energy spectrum phonon-like at small momenta, ϵ(p) = c_sp. The critical velocity is in this case finite, v_c = c_s, and we get superfluidity.

As w approaches and exceeds the critical velocity, the superfluid flow begins to decay. This happens because a tangle of superfluid vortices, so-called quantum turbulence, forms and interacts with the excitations that make up the normal fluid, resulting in a dissipative mutual friction between the normal and superfluid components (Skrbek 2011; Skrbek & Sreenivasan 2012; Barenghi et al. 2014). This effect is not included in the equations for superfluid hydrodynamics and must be added through additional terms. However, this would require us to assume the dependence of this force on the fluid variables and specify the extra parameters introduced to our model (for examples of this in numerical studies of superfluid helium, see Doi et al. 2008; Darve et al. 2012; Soulaine et al. 2017). To capture the basic consequence of Landau’s criterion relevant for this work, which is that the counterflow w is limited by the critical velocity, we instead assume the mutual friction only takes place once the critical velocity is exceeded, and that the complicated processes taking place happen on time and length scales much shorter than we are considering. The mutual friction is therefore effectively instantaneous, and since it is dissipative, there is a conversion of kinetic energy into internal energy, heating the fluid and generating entropy. Stated more precisely, we enforce the superfluid critical velocity at every position in our numerical scheme by converting kinetic energy of the two fluid components (while conserving the total momentum) into internal energy and generated entropy so that w <  v_c is always satisfied. See Appendix B.3 for further details.

We must also specify the equation of state (EOS) that defines how the thermodynamic quantities depend on the temperature and particle density. In superfluids, the EOS is also a function of the counterflow w (Landau & Lifshitz 1987; Khalatnikov 2000), but we neglected this dependence and used the EOS corresponding to the w = 0 limit. While this work is motivated by the superfluid DM model presented by Berezhiani & Khoury (2015), it lacks a complete EOS at finite temperatures. We therefore employed the weakly interacting Bose gas with effective repulsive two- and three-body contact interactions as described by Sharma et al. (2019), where the three-body case corresponds most closely to the model by Berezhiani & Khoury (2015). Effective contact interactions can describe the s-wave scattering limit of more complicated interactions through the Born approximation, which makes this class of models a more general description of superfluids (Pethick & Smith 2008). The coupling term between the DM fluid and baryons that gives rise to the emergent fifth force is not included in this work. For computational speed, we approximated the EOS in the sub-T_c regime by an ideal Bose gas with contributions from interactions at zero temperature. Notably, the superfluid fraction is approximated as the fraction of condensed particles in an ideal BEC, f_s = ρ_s/ρ = 1 − (T/T_c)^3/2. This might appear paradoxical since we already stated that an ideal Bose gas cannot be superfluid, but in the weakly interacting gas these quantities can be seen from Fig. 1 to be closely related. For strong interactions, the superfluid fraction can approach unity while the condensate fraction remains small, which is the case in superfluid ⁴He (Glyde 2013), but this scenario is outside the scope of this paper. See Appendix A for further details on the EOS.

Fig. 1. Superfluid fraction fs for the two-body interacting Bose gas calculated by Sharma et al. (2019) with ρ = 108 ρc0 and m = 1 eV compared to the approximation fs, ap. = 1 − (T/Tc)3/2. For sufficiently weak interactions, fs can be approximated by the condensate fraction in an ideal Bose gas.

2.3. Super-comoving variables

Since we are interested in the evolution of the superfluid in an expanding space, we introduce the peculiar velocity v = u − Hr and super-comoving variables (Martel & Shapiro 1998), denoted by a tilde-sign, to rewrite the hydrodynamic equations in a more convenient form¹:

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

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

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

$$\begin{array}{cc}& \frac{\partial \mathit{j}}{\partial t}+\mathbf{\nabla }P+{\rho }_{\mathrm{s}}({\mathit{v}}_{\mathrm{s}}·\mathbf{\nabla }){\mathit{v}}_{\mathrm{s}}+{\rho }_{\mathrm{n}}({\mathit{v}}_{\mathrm{n}}·\mathbf{\nabla }){\mathit{v}}_{\mathrm{n}}\hfill \\ \hfill & +{\mathit{v}}_{\mathrm{s}}[\mathbf{\nabla }·({\rho }_{\mathrm{s}}{\mathit{v}}_{\mathrm{s}})]+{\mathit{v}}_{\mathrm{n}}[\mathbf{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}})]=-\rho \mathbf{\nabla }\varphi ,\hfill \end{array}$$(12)

$$\begin{array}{cc}& \frac{\partial E}{\partial t}+\mathbf{\nabla }·[(U+\frac{1}{2}{\rho }_{\mathrm{n}}{v}_{\mathrm{n}}^2+P){\mathit{v}}_{\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{s}}{v}_{\mathrm{s}}^2{\mathit{v}}_{\mathrm{s}}\hfill \\ \hfill & +\mu {\rho }_{\mathrm{s}}({\mathit{v}}_{\mathrm{s}}-{\mathit{v}}_{\mathrm{n}})]=-H(3P-2U)-\mathit{j}·\mathbf{\nabla }\varphi .\hfill \end{array}$$(13)

The super-comoving quantities are re-scaled to reduce the dependence on the scale factor a, with the variables defined as before: $\stackrel{\mathbf{\sim }}{\mathit{j}}={\stackrel{\sim }{\rho }}_{\mathrm{n}}{\stackrel{\mathbf{\sim }}{\mathit{v}}}_{\mathrm{n}}+{\stackrel{\sim }{\rho }}_{\mathrm{s}}{\stackrel{\mathbf{\sim }}{\mathit{v}}}_{\mathrm{s}}$ and $\stackrel{\sim }{E}=\stackrel{\sim }{U}+\frac{1}{2}{\stackrel{\sim }{\rho }}_{\mathrm{n}}{\stackrel{\sim }{v}}_{\mathrm{n}}^2+\frac{1}{2}{\stackrel{\sim }{\rho }}_{\mathrm{s}}{\stackrel{\sim }{v}}_{\mathrm{s}}^2$. The only real difference is the peculiar gravitational potential $\stackrel{\sim }{\varphi }$ that is now given by (in a flat universe with matter and a cosmological constant)

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathrm{\nabla }}}^2\stackrel{\sim }{\varphi }=6a(\stackrel{\sim }{\rho }-1).\end{array}$$(14)

$\stackrel{\sim }{H}$ is the super-comoving Hubble parameter.

2.4. Linear perturbation expansion

The superfluid hydrodynamic equations at linear order can tell us a lot about the collapse of a superfluid, in particular how it will differ from CDM and non-superfluid thermal DM. The fluid variables are expanded around their background values, $\stackrel{\sim }{\rho }={\stackrel{\sim }{\rho }}_0+\delta \stackrel{\sim }{\rho }$, $\stackrel{\sim }{S}={\stackrel{\sim }{S}}_0+\delta \stackrel{\sim }{S,}$ etc. The peculiar background velocities are zero, so ${\stackrel{\sim }{\mathit{v}}}_{\mathrm{s}}=\delta {\stackrel{\sim }{\mathit{v}}}_{\mathrm{s}}$, ${\stackrel{\sim }{\mathit{v}}}_{\mathrm{n}}=\delta {\stackrel{\sim }{\mathit{v}}}_{\mathrm{n}}$, and $\stackrel{\sim }{\mathit{j}}=\delta \stackrel{\sim }{\mathit{j}}$. We also have ${\stackrel{\sim }{\mathrm{\nabla }}}^2\delta \stackrel{\sim }{\varphi }=6a\delta \stackrel{\sim }{\rho }$. This gives the following linear equations;

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

$$\begin{array}{cc}& \frac{\partial \delta \stackrel{\sim }{S}}{\partial \stackrel{\sim }{t}}+{\stackrel{\sim }{S}}_0\stackrel{\sim }{\mathbf{\nabla }}·\delta {\stackrel{\sim }{\mathit{v}}}_{\mathrm{n}}=0,\hfill \end{array}$$(16)

$$\begin{array}{cc}& \frac{\partial \delta {\stackrel{\sim }{\mathit{v}}}_{\mathrm{s}}}{\partial \stackrel{\sim }{t}}+\stackrel{\sim }{\mathbf{\nabla }}(\delta \stackrel{\sim }{\mu }+\delta \stackrel{\sim }{\varphi })=\mathbf{0},\hfill \end{array}$$(17)

$$\begin{array}{cc}& \frac{\partial \delta \stackrel{\sim }{\mathit{j}}}{\partial \stackrel{\sim }{t}}+\stackrel{\sim }{\mathbf{\nabla }}\delta \stackrel{\sim }{P}+\stackrel{\sim }{{\rho }_0}\stackrel{\sim }{\mathbf{\nabla }}\delta \stackrel{\sim }{\varphi }=\mathbf{0}.\hfill \end{array}$$(18)

These can be combined into two coupled equations for $\delta \stackrel{\sim }{\rho }$ and $\delta \stackrel{\sim }{S}$ in $\stackrel{\sim }{k}$-space;

$$\begin{array}{c}\hfill \frac{{\partial }^2\delta {\stackrel{\sim }{\rho }}_{\stackrel{\sim }{k}}}{\partial {\stackrel{\sim }{t}}^2}+[{\left(\frac{\partial \stackrel{\sim }{P}}{\partial \stackrel{\sim }{\rho }}\right)}_0{\stackrel{\sim }{k}}^2-6a{\stackrel{\sim }{\rho }}_0]\delta {\stackrel{\sim }{\rho }}_{\stackrel{\sim }{k}}+{\left(\frac{\partial \stackrel{\sim }{P}}{\partial \stackrel{\sim }{S}}\right)}_0{\stackrel{\sim }{k}}^2\delta {\stackrel{\sim }{S}}_{\stackrel{\sim }{k}}=0,\end{array}$$(19)

$$\begin{array}{cc}\hfill \frac{{\partial }^2\delta {\stackrel{\sim }{S}}_{\stackrel{\sim }{k}}}{\partial {\stackrel{\sim }{t}}^2}& +{\stackrel{\sim }{S}}_0[\{\frac{1}{{\stackrel{\sim }{\rho }}_0}{\left(\frac{\partial \stackrel{\sim }{P}}{\partial \stackrel{\sim }{\rho }}\right)}_0+\frac{{\stackrel{\sim }{S}}_0{\stackrel{\sim }{\rho }}_{s0}}{{\stackrel{\sim }{\rho }}_0{\stackrel{\sim }{\rho }}_{n0}}{\left(\frac{\partial \stackrel{\sim }{T}}{\partial \stackrel{\sim }{\rho }}\right)}_0\}{\stackrel{\sim }{k}}^2-6a]\delta {\stackrel{\sim }{\rho }}_{\stackrel{\sim }{k}}\hfill \\ \hfill & +{\stackrel{\sim }{S}}_0[\frac{1}{{\stackrel{\sim }{\rho }}_0}{\left(\frac{\partial \stackrel{\sim }{P}}{\partial \stackrel{\sim }{S}}\right)}_0+\frac{{\stackrel{\sim }{S}}_0{\stackrel{\sim }{\rho }}_{s0}}{{\stackrel{\sim }{\rho }}_0{\stackrel{\sim }{\rho }}_{n0}}{\left(\frac{\partial \stackrel{\sim }{T}}{\partial \stackrel{\sim }{S}}\right)}_0]{\stackrel{\sim }{k}}^2\delta {\stackrel{\sim }{S}}_{\stackrel{\sim }{k}}=0,\hfill \end{array}$$(20)

where the subscript “0” indicates the background values.

We would like to enforce the critical velocity in the linear approach, though we cannot do it in the same way as for the full hydrodynamic equations. Since the effect of the critical velocity is to essentially restrict the two-fluid nature of the superfluid, forcing the whole fluid to evolve like a normal fluid, we can, as a rough approximation, set ρ_s = 0 and ρ_n = ρ once ${\stackrel{\sim }{w}}_{\stackrel{\sim }{k}}\ge {\stackrel{\sim }{v}}_{\mathrm{c}}$, where ${\stackrel{\sim }{w}}_{\stackrel{\sim }{k}}$ is the relative velocity of mode $\stackrel{\sim }{k}$ and evolves at linear order according to

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{w}}_{\stackrel{\sim }{k}}}{\partial \stackrel{\sim }{t}}=\stackrel{\sim }{k}\frac{{\stackrel{\sim }{S}}_0}{{\stackrel{\sim }{\rho }}_{n0}}[{\left(\frac{\partial \stackrel{\sim }{T}}{\partial \stackrel{\sim }{\rho }}\right)}_0\delta {\stackrel{\sim }{\rho }}_{\stackrel{\sim }{k}}+{\left(\frac{\partial \stackrel{\sim }{T}}{\partial \stackrel{\sim }{S}}\right)}_0\delta {\stackrel{\sim }{S}}_{\stackrel{\sim }{k}}].\end{array}$$(21)

This approximation is further justified by the fact that the critical velocity decreases with the DM density. Once w reaches v_c, it only becomes smaller in the linear regime, forcing the superfluid to behave even more like a normal fluid.

A few qualitative statements can be made from Eqs. (19) and (20). Both mass density and entropy perturbations grow due to gravity, but this growth is slowed by pressure terms that are scale dependent through the ${\stackrel{\sim }{k}}^2$ factor, as expected in a self-gravitating fluid with nonzero pressure. In a superfluid, however, there are additional effective pressure terms that suppress the growth of entropy perturbations, and hence thermal pressure, that are absent in conventional fluids. This in turn allows the mass density perturbations to collapse more efficiently, even though the DM fluid may have relatively high temperatures. The reason for this behavior is the superfluid component’s attraction to higher temperatures. The normal component tends to transport mass and entropy from hot to cold regions, while the superfluid tends to flow in the opposite direction and balance the mass-loss due to the normal component, resulting in a thermal flux that can be large compared to the net mass flux. This effect, called thermal counterflow, makes superfluids very efficient at conducting heat.

3. Results and discussion

The hydrodynamic equations were integrated numerically using a modified first-order FORCE scheme (see Toro 2006 and Appendix B for further details) for a spherically symmetric system with an initial density contrast of the form $\delta \stackrel{\sim }{\rho }/{\stackrel{\sim }{\rho }}_0={\mathrm{\Delta }}_0{e}^{-{(\stackrel{\sim }{r}/\stackrel{\sim }{L})}^2}$ and $\delta \stackrel{\sim }{S}=0$, where $\stackrel{\sim }{L}$ is the size of the overdensity. The initial state is at approximately the same T/T_c, and hence the same mixing fraction of the normal and superfluid components, throughout the system. A flat ΛCDM background cosmology with Ω_m0 = 0.3, Ω_Λ0 = 0.7, and h = 0.7 was used, and the integration started at redshift z = 1000 with Δ₀ = 5 × 10⁻³. An example of a collapsing SFDM halo at various redshifts can be seen in Fig. 2, illustrating that as the halo collapses, a thermal counterflow carrying entropy away from the halo center develops, slowing down the growth of entropy until the critical velocity is reached and the fluid starts heating up.

Fig. 2. Profiles of a collapsing SFDM halo with an initial Gaussian density contrast, m = 30 eV, g = 10−5 eV−2, L = 100 kpc, and T/Tc = 0.1. A thermal counterflow develops and the growth of entropy perturbations is at first suppressed. This also gives a slight decrease in the ratio T/Tc, and hence the superfluid fraction, since fs = ρs/ρ = 1 − (T/Tc)3/2. As the critical velocity is reached, entropy is generated, and T/Tc increases.

3.1. Growth of structure

In Figs. 3–5, the redshift when the central density contrast reaches 200 is shown for various parameters for both the superfluid and non-superfluid (a conventional fluid with ρ_s = 0, ρ_n = ρ, and the same EOS) cases. While the growth of structure is slower compared to CDM, the SFDM halos collapse more efficiently than their non-superfluid counterparts as the interaction strength is increased until a maximum is reached, after which the growth of structure in both super-and non-superfluid DM is suppressed. This is counter to what one would intuitively expect, since an increase in interactions also means an increase in pressure. It can, however, be understood as follows: for small interactions, the superfluid behaves nearly the same as a normal fluid because the critical velocity, which scales as $v_{\mathrm{c}}\sim \sqrt{g_i}$, is reached very early. When this happens, the flow of the normal and superfluid components become “locked” to one another, unable to efficiently conduct heat away from the halo core. As the interaction increases, the thermal counterflow can both be larger and last longer, resulting in an increased suppression of thermal gradients and thus allows for a faster collapse. For sufficiently large interactions, the collapse is instead suppressed due to large zero-temperature pressure gradients that the superfluid is unable to wash out.

Fig. 3. Comparison of redshifts when the central density contrast reaches 200 as function of the interaction strength for various particle masses, with T/Tc = 0.1 and L = 100 kpc. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown. For constant T/Tc, the temperature is increased for decreasing mass, since Tc ∼ m−5/3. The comparison of the collapse for various masses is therefore not done at the same temperature, but instead at a similar place in the superfluid phase.

Fig. 4. Comparison of redshifts when the central density contrast reaches 200 as function of the interaction strength for various scales, with T/Tc = 0.1 and m = 50 eV. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown.

Fig. 5. Comparison of redshifts when the central density contrast reaches 200 as function of the interaction strength for various temperatures, with m = 50 eV and L = 100 kpc. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown.

Most production of entropy due to mutual friction as the Landau criterion is broken takes place away from the center of the halo. The resulting extra thermal pressure acting on the interior causes the central density contrast to grow slightly faster and can best be seen by the gap between collapse times of the superfluid and non-superfluid cases at low g₂. If entropy was not produced, this gap would vanish.

3.2. Dependence on equation of state

The Bose gas with two-body interactions is compared with three-body interactions in Fig. 6. The same qualitative behavior is present in both cases and is expected to be a general feature regardless of the EOS used, as long as there is superfluidity. In the linear expansion of the superfluid equations, Eqs. (19) and (20), the additional effective pressure terms due to a superfluid component require only the temperature to be dependent on mass density or entropy. Indeed, the approximated two-body and three-body EOS used in this work both have a temperature profile that is independent of mass density for T <  T_c, so that one of the effective pressure terms in Eq. (20) is absent. For EOS where the temperature is dependent on both the mass density and entropy, the collapse of SFDM may be even more efficient.

Fig. 6. Comparison of redshifts when the central density contrast reaches 200 for two-body and three-body interactions as function of the interaction strengths g2 and g3, respectively, with m = 30 eV, L = 100 kpc, and T/Tc = 0.1. The three-body interaction is multiplied by $\sqrt{10}\times {10}^{-4}$ to make the comparison clearer. Both the superfluid case (solid lines) and the non-superfluid (striped lines) are shown.

3.3. Effect of small-scale and nonradial motion

In this work, we assumed perfect radial infall of DM. The relative velocity w is simply the difference between the radial velocities of the two fluid components. In a real system, there is expected to be additional small-scale motion in all directions, such as turbulence that our simplified model averages over. The superfluid critical velocity may therefore be exceeded on small scales, while the large-scale radial average only appears to have w <  v_c. In this case, the superfluid would behave like a conventional fluid at much smaller w. In other words, there is an effective superfluid critical velocity $v_{\mathrm{c}}^{\text{eff}}<v_{\mathrm{c}}$ that is a decreasing function of the local turbulence. This leads to a difference in collapse times of halos with different amounts of turbulence, the turbulent ones collapsing at a slower rate, as seen in Fig. 7.

Fig. 7. Comparison of redshifts when the central density contrast reaches 200 for various effective critical velocities as function of the interaction strength, with m = 30 eV, T/Tc = 0.1, and L = 100 kpc. Both the superfluid case (solid lines) and the non-superfluid (striped line) are shown.

3.4. Evolution of superfluid fraction

In a conventional fluid, the entropy and mass density collapses at the same rate so that the ratio T/T_c is constant. A fluid that is initially in the normal phase will therefore remain so. A collapsing superfluid, on the other hand, experiences an increase in the superfluid fraction due to thermal counterflow until the critical velocity is reached. At this point, entropy is generated causing T/T_c to rise, and thus the superfluid fraction to fall; though it takes time for the full effect of this to propagate to the center of the halo, as shown in Figs. 2 and 8. It may be, however, that Eqs. (1)–(5) do not properly describe super-critical flow, and too much entropy is generated in our numerical scheme for enforcing the critical velocity. The evolution of T/T_c when no entropy is generated is therefore also shown in Fig. 8 as the opposite extreme. This case behaves similarly until near the end of the collapse, where T/T_c rises only modestly. Profiles are shown in Fig. 9, which corresponds to Fig. 2 with no production of entropy.

Fig. 8. Evolution of T/Tc in the halo center during collapse for various masses and initial temperatures with g = 10−5 eV−2 and L = 100 kpc. Both the evolution with entropy production (solid lines) and without (striped lines) are shown until the overdensity reach 105. The two cases differ only in the end stage of the collapse, well after the critical velocity is first reached, indicated by the colored vertical lines.

Fig. 9. Profiles of a collapsing SFDM halo where no entropy is produced as w = vv with an initial Gaussian density contrast, m = 30 eV, g = 10−5 eV−2, L = 100 kpc, and T/Tc = 0.1.

The decrease in T/T_c during collapse becomes smaller as the temperature approaches T_c, where the superfluid fraction goes to zero and thermal counterflow becomes inefficient. The formation of DM halos with much higher superfluid fractions than the background, as required in the emergent MOND scenario of Berezhiani & Khoury (2015), therefore appears unlikely through collapse alone. Additional cooling mechanisms during or after collapse are necessary.

3.5. Dark matter self-interaction constraints

The distribution of DM, gas, and stellar mass in cluster collisions provides constraints on the cross-section of DM self-interactions, σ/m <  0.5 cm² g⁻¹ (Harvey et al. 2015). In terms of the two-body interaction strength, this corresponds to (Pitaevskii & Stringari 2016)

$$\begin{array}{c}\hfill g_2=\sqrt{4\pi \sigma }\frac{{ħ}^2}{m}<5\times {10}^{-12}{\left(\frac{1\text{eV}}{m}\right)}^{1/2}{\text{eV}}^{-2}.\end{array}$$(22)

The values of g₂ in the above results do not generally satisfy this constraint, but we chose to relax it since we do not know how it translates to SFDM. In any case, the above features were also found for smaller g₂ using perturbation theory (while simultaneously lowering m and T/T_c) that do satisfy the constraints, as is exemplified in Fig. 10.

Fig. 10. Redshifts when the linear density contrast for the mode k = 2/100 kpc−1 with T/Tc = 2 × 10−6 reaches unity for various masses and interaction strengths. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown, illustrating that the same features can be found for a choice of parameters that satisfy the constraint from cluster collisions on DM mass and self-interaction.

4. Conclusions

When superfluid behavior is included in a finite-temperature DM fluid, the formation of structure is found to be much more efficient in certain regions of parameter space than one would, naively, expect, through it is still slower compared to CDM. The effect of thermal counterflow is most prominent when the thermal suppression is large, such as at small scales and relatively high temperatures. The increased collapse efficiency is also expected to be a general feature of SFDM regardless of the EOS used, though the specific model in question will certainly affect the finer details through the dependence of entropy, pressure, and critical velocity on temperature, mass density, and the model parameters. The toy models used in this work were motivated by condensed matter physics, but suffer some severe limitations at high redshifts. Both are derived under the assumption that the interactions are weak and the number density is not too large, which is invalid at very early times. Furthermore, the zero-temperature pressure depends on the number density through n² and n³, resulting in very high pressures at high redshifts that might wash out the initial perturbations set up by inflation. The generalization of this work to more exotic DM fluids and adding interactions between DM and baryons, which has recently been considered in the literature, is therefore of interest in the further study of SFDM models. It may also be of interest to study the case when thermal equilibrium is not always assumed so that the DM fluid can fall in and out of equilibrium and superfluidity can vanish and reappear.

Superfluid models of DM involve processes that require the superfluid hydrodynamic equations to be properly described. Throughout this work, spherical symmetry was assumed, but nonradial and turbulent motion is expected to have a significant impact on the superfluid dynamics, especially through the critical velocity, which is broken at smaller radial thermal counterflows. It is also important to understand the effect of mergers in SFDM. Large-scale and high-resolution simulations will therefore be essential for the further study of structure formation. The main challenge in this line of inquiry may be developing numerical schemes that are faster and more accurate than the modified first-order FORCE scheme used in this work that can capture the small-scale motion of the superfluid and its effect on structure formation.

Acknowledgments

We thank the Research Council of Norway for their support. We are grateful for the helpful comments and suggestions made by the anonymous referee.

References

Appendix A: equation of state

An EOS for a weakly interacting Bose gas valid at all temperatures was recently proposed by Sharma et al. (2019) for two-body and three-body interactions. Since we do not know the true EOS of DM and must resort to toy models, we instead approximated the EOS by using an ideal Bose gas with zero-temperature contributions from weak interactions. At very low temperatures, this approximation breaks down as the interactions become increasingly important, but we generally remain well above this regime.

An important quantity is the critical temperature T_c, above which the fluid behaves as a normal fluid, while below the fluid condenses into a BEC and becomes superfluid;

$$\begin{array}{c}\hfill T_{\mathrm{c}}=\frac{2\pi {ħ}^2}{m{k}_{\text{B}}}{\left(\frac{n}{\zeta (3/2)}\right)}^{2/3}=\frac{2\pi {ħ}^2}{m^{5/3}{k}_{\text{B}}}{\left(\frac{\rho }{\zeta (3/2)}\right)}^{2/3},\end{array}$$(A.1)

where n = ρ/m is the particle number density.

As an estimate for the superfluid fraction f_s = ρ_s/ρ we use the fraction of particles in the BEC in an ideal Bose gas;

$$\begin{array}{c}\hfill f_{\mathrm{s}}=\{\begin{array}{c}1-{\left(\frac{T}{T_{\mathrm{c}}}\right)}^{3/2},T\le T_{\mathrm{c}}\hfill \\ 0,T>T_{\mathrm{c}}.\hfill \end{array}\end{array}$$(A.2)

For the other thermodynamic quantities, such as pressure, entropy, etc., we must consider them above and below T_c separately. Both two-body and three-body interactions are given, parameterized by g₂ and g₃, respectively.

A.1. T  >   T_c

The pressure is given by

$$\begin{array}{c}\hfill P=\{\begin{array}{cc}{g}_2{n}^2+\frac{\sqrt{2}\mathrm{\Gamma }(5/2){(k_{\text{B}}T)}^{5/2}{m}^{3/2}}{3{\pi }^2{ħ}^3}{\text{Li}}_{5/2}\left(e^{\beta (\mu -2g_{2}n)}\right),\hfill & \text{two-body}\hfill \\ 4g_3{n}^3+\frac{\sqrt{2}\mathrm{\Gamma }(5/2){(k_{\text{B}}T)}^{5/2}{m}^{3/2}}{3{\pi }^2{ħ}^3}{\text{Li}}_{5/2}\left(e^{\beta (\mu -6g_3{n}^2)}\right),\hfill & \text{three-body},\hfill \end{array}\end{array}$$(A.3)

where Γ(x) is the gamma function, Li_z(x) is the polylogarithmic function, β = 1/k_BT, and the chemical potential μ is determined by the equation for the number density

$$\begin{array}{c}\hfill n=\{\begin{array}{cc}\frac{\sqrt{2}\mathrm{\Gamma }(5/2)}{3{\pi }^2{ħ}^3}{(k_{\text{B}}T)}^{3/2}{m}^{3/2}{\text{Li}}_{3/2}\left(e^{\beta (\mu -2g_{2}n)}\right),\hfill & \text{two-body}\hfill \\ \frac{\sqrt{2}\mathrm{\Gamma }(5/2)}{3{\pi }^2{ħ}^3}{(k_{\text{B}}T)}^{3/2}{m}^{3/2}{\text{Li}}_{3/2}\left(e^{\beta (\mu -6g_3{n}^2)}\right),\hfill & \text{three-body}.\hfill \end{array}\end{array}$$(A.4)

The entropy is

$$\begin{array}{c}\hfill S=\{\begin{array}{cc}\frac{5}{2}\frac{P}{T}-n\beta (\mu -2g_{2}n),\hfill & \text{two-body}\hfill \\ \frac{5}{2}\frac{P}{T}-n\beta (\mu -6g_3{n}^2),\hfill & \text{three-body},\hfill \end{array}\end{array}$$(A.5)

and the internal energy is

$$\begin{array}{c}\hfill U=ST-P+\mu n.\end{array}$$(A.6)

The sound speed used when determining the time-stepping in the numerical scheme was

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\sqrt{\frac{5}{3}\frac{k_{\text{B}}T}{m}}·\end{array}$$(A.7)

In the limit of very high temperature, these reduce to the classical ideal gas.

A.2. T ≤ T_c

The EOS below the critical temperature is given by the ideal Bose gas plus some zero-temperature contributions due to interactions;

$$\begin{array}{cc}\hfill P& =\{\begin{array}{cc}\frac{1}{2}{g}_2{n}^2+\zeta (5/2){\left(\frac{m}{2\pi {ħ}^2}\right)}^{3/2}{\left(k_{\text{B}}T\right)}^{5/2},\hfill & \text{two-body}\hfill \\ \frac{2}{3}{g}_3{n}^3+\zeta (5/2){\left(\frac{m}{2\pi {ħ}^2}\right)}^{3/2}{\left(k_{\text{B}}T\right)}^{5/2},\hfill & \text{three-body},\hfill \end{array}\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill S& =\frac{5}{2}\zeta (5/2){\left(\frac{m}{2\pi {ħ}^2}\right)}^{3/2}{k}_{\text{B}}^{5/2}{T}^{3/2},\text{two-and three-body},\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill \mu & =\{\begin{array}{cc}{g}_{2}n,\hfill & \text{two-body}\hfill \\ g_3{n}^2,\hfill & \text{three-body},\hfill \end{array}\hfill \end{array}$$(A.10)

and the internal energy is again given by Eq. (A.6). The fastest sound speed was approximated using

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\sqrt{\frac{\zeta (5/2)}{\zeta (3/2)}\frac{5}{3}\frac{k_{\text{B}}T}{m}},\end{array}$$(A.11)

and the critical velocity given by

$$\begin{array}{c}\hfill v_{\mathrm{c}}=\{\begin{array}{cc}\sqrt{\frac{g_{2}n}{m}[1-{(T/T_{\mathrm{c}})}^{3/2}]},\hfill & \text{two-body}\hfill \\ \sqrt{\frac{2g_3{n}^2}{m}[1-{(T/T_{\mathrm{c}})}^{3/2}][1+2{(T/T_{\mathrm{c}})}^{3/2}]},\hfill & \text{three-body.}\hfill \end{array}\end{array}$$(A.12)

There is a small discontinuity at the critical temperature, with μ = 2ng₂ above and μ = ng₂ below for the two-body interaction (and a similar jump in zero-temperature pressure and internal energy). There should be a crossover region as the condensate fraction increases, but during this crossover the thermal contributions dominates and the discontinuity is negligible.

Appendix B: Numerical scheme

In this work, we employed a modified first-order FORCE scheme (Toro 2006) – an incomplete Riemann solver – for the superfluid hydrodynamic equations with source terms due to gravity and from using spherical coordinates. The source terms were evaluated at two stages during each time-step: once before the advection step, and once after, at which point the average of the two evaluations was added to the solution. Gravity was also evaluated with half a time-step when computing fluxes during the advection step. Finally, we enforced the critical velocity, which was done in three stages; once when computing fluxes, once after the fluxes from the advection step were applied, and a final time after the source terms were applied. Further details are presented below.

For spherical collapse, this scheme was found to be sufficient since the solutions are mostly smooth, evolve slowly, and are one-dimensional. For more complex and higher dimensional cases where shock fronts arise and the solutions undergo fast changes, this scheme is expected to perform suboptimally, primarily because it is first-order. There is a well-known way to increase the order and thus accuracy of the scheme through slope reconstruction and slope limiters. However, instabilities arose when the superfluid component was included, and adding further restrictions to the reconstructed slopes with modified slope limiters failed to fix this. Slope reconstruction was therefore not used.

B.1. First-order FORCE scheme

The FORCE scheme is a variant of Godunov’s method for solving partial differential equations. In this method, the domain is divided into finite-volume elements, or cells, and the Riemann problem at each cell interface is solved. The Riemann problem is the initial value problem with two piece-wise constant initial regions connected by a discontinuity, then asking how this evolves in time and what the net flux across the interface is. The scheme for computing or approximating this flux is called a Riemann solver and is what characterizes the different ways of implementing Godunov’s method.

To see how this works, one can consider the m-component state vector U that obeys the one-dimensional conservative equation

$$\begin{array}{c}\hfill {\partial }_t\mathit{U}+{\partial }_x\mathit{F}(\mathit{U})=\mathbf{0},\end{array}$$(B.1)

where F is the flux. By integrating over the time interval [tⁿ, t^n + 1] and cell-volume [x_{i − 1/2}, x_i + 1/2], we get

$$\begin{array}{c}\hfill {\mathit{U}}_i^{n+1}={\mathit{U}}_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}[{\mathit{F}}_{i-1/2}-{\mathit{F}}_{i+1/2}],\end{array}$$(B.2)

where

$$\begin{array}{cc}& {\mathit{U}}_i^n=\frac{1}{\mathrm{\Delta }x}{\displaystyle {\int }_{x_{i-1/2}}^{x_{i+1/2}}}\mathit{U}(x,t^n)\text{d}x,\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& {\mathit{F}}_{i+1/2}=\frac{1}{\mathrm{\Delta }t}{\displaystyle {\int }_{t^n}^{t^{n+1}}}\mathit{F}(x_{i+1/2},t)\text{d}t.\hfill \end{array}$$(B.4)

In the first-order Godunov scheme, the state U is assumed to be piece-wise constant in each cell, given by the cell average ${\mathit{U}}_i^n$. To compute F_i + 1/2, the states on the left and right sides of the interface is used, ${\mathit{U}}_{i+1/2,\text{L}}={\mathit{U}}_i^n$ and ${\mathit{U}}_{i+1/2,\text{R}}={\mathit{U}}_{i+1}^n$, and the corresponding Riemann problem is solved or approximated. The time-step is chosen so that no signal in the domain travels further than one cell length Δx. This is given by a Courant-Friedrich-Lewy (CFL) type condition

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\text{s}}=C_{\text{CFL}}^{\text{s}}\frac{\mathrm{\Delta }x}{v_{\text{max}}},\end{array}$$(B.5)

where v_max is the maximum signal speed in the domain, and $C_{\text{CFL}}^{\text{s}}$ is a number less than one that controls how far across a cell the fastest signal is allowed to move during each time-step. In simulations with gravity and expansion, additional constraints need to be added to the time-stepping. For gravity, the free-fall distance in each cell, with acceleration g, must be smaller than the cell lengths,

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\text{ff}}=C_{\text{CFL}}^{\text{ff}}\sqrt{\frac{2\mathrm{\Delta }x}{g}},\end{array}$$(B.6)

and for expansion, the relative change in the scale factor is restricted:

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\text{exp}}=C_{\text{CFL}}^{\text{exp}}\frac{1}{H}·\end{array}$$(B.7)

Here, $C_{\text{CFL}}^{\text{ff}}$ and $C_{\text{CFL}}^{\text{exp}}$ are also numbers less than one. In this work, we used $C_{\text{CFL}}^{\text{s}}=0.5$, $C_{\text{CFL}}^{\text{ff}}=0.5$, and $C_{\text{CFL}}^{\text{exp}}=0.01$. The final value for the time-step is the smallest of the above,

$$\begin{array}{c}\hfill \mathrm{\Delta }t=min[\mathrm{\Delta }{t}_{\text{s}},\mathrm{\Delta }{t}_{\text{ff}},\mathrm{\Delta }{t}_{\text{exp}}].\end{array}$$(B.8)

The FORCE scheme approximates the interface flux F (given the left and right states U_L and U_R) as the average of the Lax–Friedrichs flux and the two-step Lax–Wendroff flux;

$$\begin{array}{cc}& {\mathit{F}}^{\text{FORCE}}=\frac{1}{2}[{\mathit{F}}^{\text{LF}}+{\mathit{F}}^{\text{LW}}],\hfill \\ \hfill & {\mathit{F}}^{\text{LF}}=\frac{1}{2}[\mathit{F}({\mathit{U}}_{\mathrm{L}})+\mathit{F}({\mathit{U}}_{\mathrm{R}})]-\frac{1}{2}\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}[{\mathit{U}}_{\mathrm{R}}-{\mathit{U}}_{\mathrm{L}}],\hfill \\ \hfill & {\mathit{F}}^{\text{LW}}=\mathit{F}({\mathit{U}}^{\text{LW}}),\hfill \\ \hfill & {\mathit{U}}^{\text{LW}}=\frac{1}{2}[{\mathit{U}}_{\mathrm{L}}+{\mathit{U}}_{\mathrm{R}}]-\frac{1}{2}\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}[\mathit{F}({\mathit{U}}_{\mathrm{R}})-\mathit{F}({\mathit{U}}_{\mathrm{L}})].\hfill \end{array}$$(B.9)

We modified this by enforcing the critical velocity on the intermediate state U^LW before computing the flux F^LW.

B.2. Sources

Gravity and extra terms when using spherical coordinates and super-comoving variables appear as source terms S in the superfluid equations. Continuing with the above example, we have

$$\begin{array}{c}\hfill {\partial }_t\mathit{U}+{\partial }_x\mathit{F}(\mathit{U})=\mathit{S}.\end{array}$$(B.10)

To modify our Godunov scheme to incorporate the sources in the flux, we did the following: at the beginning of each time-step, we had the states ${\mathit{U}}_i^n$. To do the advection (the Godunov step), we input the left and right states at each boundary i + 1/2; U_{i + 1/2, L} = U_i, U_{i + 1/2, R} = U_i + 1. But before we computed the interface flux, we applied half a time-step of the source due to gravity,

$$\begin{array}{cc}& {\mathit{U}}_{i+1/2,\mathrm{L}}^{\ast }={\mathit{U}}_{i+1/2,\mathrm{L}}+\frac{1}{2}\mathrm{\Delta }t{\left({\mathit{S}}_{i+1/2,\mathrm{L}}^n\right)}_{\text{grav}},\hfill \\ \hfill & {\mathit{U}}_{i+1/2,\mathrm{L}}^{\ast }={\mathit{U}}_{i+1/2,\mathrm{R}}+\frac{1}{2}\mathrm{\Delta }t{\left({\mathit{S}}_{i+1/2,\mathrm{R}}^n\right)}_{\text{grav}},\hfill \end{array}$$(B.11)

where ${\mathit{S}}_{i+1/2,\text{L}}^n$ and ${\mathit{S}}_{i+1/2,\text{R}}^n$ are the left and right values for the sources. In this work, we computed these using the average gravitational acceleration ${(\mathbf{\nabla }\mathrm{\Phi })}_{i+1/2,\mathrm{L}/\mathrm{R}}^n=\frac{1}{2}[{(\mathbf{\nabla }\mathrm{\Phi })}_i^n+{(\mathbf{\nabla }\mathrm{\Phi })}_{i+1}^n]$, and the left and right states U_{i + 1/2, L/R}. We then used ${\mathit{U}}_{i+1/2,\mathrm{L}}^{\ast }$ and ${\mathit{U}}_{i+1/2,\mathrm{R}}^{\ast }$ as the input states in the Godunov scheme to get F_i + 1/2, and updated the state vectors from the previous time-step:

$$\begin{array}{c}\hfill {\mathit{U}}_i^{n+1,\ast }={\mathit{U}}_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}[{\mathit{F}}_{i-1/2}-{\mathit{F}}_{i+1/2}].\end{array}$$(B.12)

This modification to the Godunov scheme was to include the effect of gravity on the flux, but explicitly adding the sources to the solution remains to be done. For this, we used the average before and after the advection step;

$$\begin{array}{c}\hfill {\mathit{U}}_i^{n+1}={\mathit{U}}_i^{n+1,\ast }+\mathrm{\Delta }t\frac{1}{2}[{\mathit{S}}_i^n+{\mathit{S}}_i^{n+1,\ast }].\end{array}$$(B.13)

B.3. Enforcing critical velocity

The critical velocity was enforced by iteratively converting kinetic energy into internal energy and generated entropy in all cells until w <  v_c. The scheme works as follows: we consider a cell with the state vector U^l, where l denotes the current step in the iterative scheme to enforce v_c, and l = 0 is the initial state. From this we get the fluid variables of the cell, $v_{\text{s}}^l$, $v_{\text{n}}^l$, S^l, etc. If $w^l<v_{\text{c}}^l$, the Landau criterion is satisfied and we do nothing. If instead $w^l>v_{\text{c}}^l$, we apply a small change $\mathrm{\Delta }{v}_{\text{s}}^l$ to $v_{\text{s}}^l$ to update it to l + 1 and decrease w,

$$\begin{array}{c}\hfill {\mathit{v}}_s^{l+1}={\mathit{v}}_s^l+\mathrm{\Delta }{\mathit{v}}_{\mathrm{s}}^l.\end{array}$$(B.14)

By keeping j constant and assuming that the change in the superfluid fraction is negligible compared to the change in velocity, we get

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{v}}_{\mathrm{n}}^l=-\frac{{\rho }_{\mathrm{s}}^l}{{\rho }_{\mathrm{n}}^l}\mathrm{\Delta }{\mathit{v}}_{\mathrm{s}}^l.\end{array}$$(B.15)

Using conservation of energy the change in internal energy is equal to the change in kinetic energy;

$$\begin{array}{cc}\hfill \mathrm{\Delta }{U}^l& =-\mathrm{\Delta }{E}_{\text{kin}}^l=-\mathrm{\Delta }(\frac{1}{2}{\rho }_{\mathrm{n}}^l|{\mathit{v}}_{\mathrm{n}}^l{|}^2+\frac{1}{2}{\rho }_{\mathrm{s}}^l{|{\mathit{v}}_{\mathrm{s}}^l|}^2)\hfill \\ \hfill & =-{\rho }_{\mathrm{s}}^l({\mathit{v}}_{\mathrm{s}}^l-{\mathit{v}}_{\mathrm{n}}^l)·\mathrm{\Delta }{\mathit{v}}_{\mathrm{s}}^l.\hfill \end{array}$$(B.16)

The change in entropy is ΔS^l = ΔQ^l/T^l, where ΔQ^l is the heating of the fluid, which in this case is just the change in internal energy:

$$\begin{array}{c}\hfill \mathrm{\Delta }{S}^l=\frac{\mathrm{\Delta }{Q}^l}{T^l}=\frac{\mathrm{\Delta }{U}^l}{T^l}=-\frac{{\rho }_{\mathrm{s}}^l({\mathit{v}}_{\mathrm{s}}^l-{\mathit{v}}_{\mathrm{n}}^l)·\mathrm{\Delta }{\mathit{v}}_{\mathrm{s}}^l}{T^l}·\end{array}$$(B.17)

The updated entropy is

$$\begin{array}{c}\hfill S^{l+1}=S^l+\mathrm{\Delta }{S}^l.\end{array}$$(B.18)

We arrive at the state vector U^l + 1 and repeat the above process until w <  v_c. The only part that needs to be specified is $\mathrm{\Delta }{\mathit{v}}_{\text{s}}^l$, which was chosen as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{v}}_{\mathrm{s}}^l=-C(w^l){\widehat{\mathit{w}}}^l,\end{array}$$(B.19)

where

$$\begin{array}{c}\hfill C(w)=[{10}^{-2},{10}^{-5}]w.\end{array}$$(B.20)

The numerical factor in Eq. (B.20) was tuned to give as smooth w-profile as possible while keeping the scheme from becoming too slow.

All Figures

	Fig. 1. Superfluid fraction fs for the two-body interacting Bose gas calculated by Sharma et al. (2019) with ρ = 108 ρc0 and m = 1 eV compared to the approximation fs, ap. = 1 − (T/Tc)3/2. For sufficiently weak interactions, fs can be approximated by the condensate fraction in an ideal Bose gas.
In the text

Fig. 2. Profiles of a collapsing SFDM halo with an initial Gaussian density contrast, m = 30 eV, g = 10−5 eV−2, L = 100 kpc, and T/Tc = 0.1. A thermal counterflow develops and the growth of entropy perturbations is at first suppressed. This also gives a slight decrease in the ratio T/Tc, and hence the superfluid fraction, since fs = ρs/ρ = 1 − (T/Tc)3/2. As the critical velocity is reached, entropy is generated, and T/Tc increases.

In the text

Fig. 3. Comparison of redshifts when the central density contrast reaches 200 as function of the interaction strength for various particle masses, with T/Tc = 0.1 and L = 100 kpc. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown. For constant T/Tc, the temperature is increased for decreasing mass, since Tc ∼ m−5/3. The comparison of the collapse for various masses is therefore not done at the same temperature, but instead at a similar place in the superfluid phase.

In the text

	Fig. 4. Comparison of redshifts when the central density contrast reaches 200 as function of the interaction strength for various scales, with T/Tc = 0.1 and m = 50 eV. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown.
In the text

	Fig. 5. Comparison of redshifts when the central density contrast reaches 200 as function of the interaction strength for various temperatures, with m = 50 eV and L = 100 kpc. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown.
In the text

Fig. 6. Comparison of redshifts when the central density contrast reaches 200 for two-body and three-body interactions as function of the interaction strengths g2 and g3, respectively, with m = 30 eV, L = 100 kpc, and T/Tc = 0.1. The three-body interaction is multiplied by $\sqrt{10}\times {10}^{-4}$ to make the comparison clearer. Both the superfluid case (solid lines) and the non-superfluid (striped lines) are shown.

In the text

	Fig. 7. Comparison of redshifts when the central density contrast reaches 200 for various effective critical velocities as function of the interaction strength, with m = 30 eV, T/Tc = 0.1, and L = 100 kpc. Both the superfluid case (solid lines) and the non-superfluid (striped line) are shown.
In the text

Fig. 8. Evolution of T/Tc in the halo center during collapse for various masses and initial temperatures with g = 10−5 eV−2 and L = 100 kpc. Both the evolution with entropy production (solid lines) and without (striped lines) are shown until the overdensity reach 105. The two cases differ only in the end stage of the collapse, well after the critical velocity is first reached, indicated by the colored vertical lines.

In the text

	Fig. 9. Profiles of a collapsing SFDM halo where no entropy is produced as w = vv with an initial Gaussian density contrast, m = 30 eV, g = 10−5 eV−2, L = 100 kpc, and T/Tc = 0.1.
In the text

Fig. 10. Redshifts when the linear density contrast for the mode k = 2/100 kpc−1 with T/Tc = 2 × 10−6 reaches unity for various masses and interaction strengths. Both the superfluid case (solid lines) and the corresponding non-superfluid case (striped lines) are shown, illustrating that the same features can be found for a choice of parameters that satisfy the constraint from cluster collisions on DM mass and self-interaction.

In the text