Issue 
A&A
Volume 639, July 2020



Article Number  A90  
Number of page(s)  10  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201937263  
Published online  14 July 2020 
Collapse of spherical overdensities in superfluid models of dark matter
Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern 0315 Oslo, Norway
email: stian.hartman@gmail.com
Received:
5
December
2019
Accepted:
7
May
2020
Aims. We intend to understand cosmological structure formation within the framework of superfluid models of dark matter with finite temperatures. Of particular interest is the evolution of smallscale structures where the pressure and superfluid properties of the dark matter fluid are prominent. We compare the growth of structures in these models with the standard cold dark matter paradigm and nonsuperfluid dark matter.
Methods. The equations for superfluid hydrodynamics were computed numerically in an expanding ΛCDM background with spherical symmetry; the effect of various superfluid fractions, temperatures, interactions, and masses on the collapse of structures was taken into consideration. We derived the linear perturbation of the superfluid equations, giving further insights into the dynamics of the superfluid collapse.
Results. We found that while a conventional dark matter fluid with selfinteractions and finite temperatures experiences a suppression in the growth of structures on smaller scales, as expected due to the presence of pressure terms, a superfluid can collapse much more efficiently than was naively expected due to its ability to suppress the growth of entropy perturbations and thus gradients in the thermal pressure. We also found that the cores of the dark matter halos initially become more superfluid during the collapse, but eventually reach a point where the superfluid fraction falls sharply. The formation of superfluid dark matter halos surrounded by a normal fluid dark matter background is therefore disfavored by the present work.
Key words: cosmology: theory / dark matter / largescale structure of Universe
© 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 largescale 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 Nbody 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 toobigtofail, cuspcore, and missing satellite problems (see e.g., Del Popolo & Le Delliou 2017; Bullock & BoylanKolchin 2017 and references therein). Another puzzling phenomenology on the scale of galaxies is the empirical baryonic TullyFisher 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 and holds for many orders of magnitude with remarkably small scatter. The ΛCDM prediction for this relation is 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 (SantosSantos 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 lowacceleration 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 longrange 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 selfinteracting 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 selfinteractions to exhibit superfluid behavior in certain mass ranges.
Observations of the largescale 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 selfinteractions and finite temperatures. We must therefore be wary of how structure forms in SFDM. Studies of other DM models with pressurelike terms, such as fuzzy dark matter (Hu et al. 2000; Schive et al. 2014; Schwabe et al. 2016; Mocz et al. 2017) and selfinteracting dark matter (Spergel & Steinhardt 2000; Elbert et al. 2015; Tulin & Yu 2018), find they can help remove the surplus of smallscale 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 finitetemperature superfluid, we employed the superfluid hydrodynamic equations (Taylor & Griffin 2005; Chapman et al. 2014), which in proper coordinates and physical variables are;
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,
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 = ρ_{n}u_{n} + ρ_{s}u_{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 . 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 nonsuperfluid 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 threedimensional homogeneous and ideal Bose gas this critical temperature is
where ζ(x) is the Riemann Zetafunction. This result holds approximately for weakly interacting gases as well (Sharma et al. 2019), apart from a small interactiondependent 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 socalled Landau criterion and requires the relative motion w = u_{s} − u_{n} to be smaller than the critical velocity v_{c},
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^{2}/2m so that v_{c} = 0. In an interacting Bose gas, on the other hand, the condensation of the gas makes the energy spectrum phononlike at small momenta, ϵ(p) = c_{s}p. 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, socalled 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 threebody contact interactions as described by Sharma et al. (2019), where the threebody case corresponds most closely to the model by Berezhiani & Khoury (2015). Effective contact interactions can describe the swave 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 subT_{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 ^{4}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 f_{s} for the twobody interacting Bose gas calculated by Sharma et al. (2019) with ρ = 10^{8} ρ_{c0} and m = 1 eV compared to the approximation f_{s, ap.} = 1 − (T/T_{c})^{3/2}. For sufficiently weak interactions, f_{s} can be approximated by the condensate fraction in an ideal Bose gas. 

Open with DEXTER 
2.3. Supercomoving variables
Since we are interested in the evolution of the superfluid in an expanding space, we introduce the peculiar velocity v = u − Hr and supercomoving variables (Martel & Shapiro 1998), denoted by a tildesign, to rewrite the hydrodynamic equations in a more convenient form^{1}:
The supercomoving quantities are rescaled to reduce the dependence on the scale factor a, with the variables defined as before: and . The only real difference is the peculiar gravitational potential that is now given by (in a flat universe with matter and a cosmological constant)
is the supercomoving 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 nonsuperfluid thermal DM. The fluid variables are expanded around their background values, , etc. The peculiar background velocities are zero, so , , and . We also have . This gives the following linear equations;
These can be combined into two coupled equations for and in space;
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 twofluid 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 , where is the relative velocity of mode and evolves at linear order according to
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 factor, as expected in a selfgravitating 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 massloss 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 firstorder FORCE scheme (see Toro 2006 and Appendix B for further details) for a spherically symmetric system with an initial density contrast of the form and , where 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 Δ_{0} = 5 × 10^{−3}. 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/T_{c} = 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/T_{c}, and hence the superfluid fraction, since f_{s} = ρ_{s}/ρ = 1 − (T/T_{c})^{3/2}. As the critical velocity is reached, entropy is generated, and T/T_{c} increases. 

Open with DEXTER 
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 nonsuperfluid (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 nonsuperfluid counterparts as the interaction strength is increased until a maximum is reached, after which the growth of structure in both superand nonsuperfluid 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 , 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 zerotemperature 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/T_{c} = 0.1 and L = 100 kpc. Both the superfluid case (solid lines) and the corresponding nonsuperfluid case (striped lines) are shown. For constant T/T_{c}, the temperature is increased for decreasing mass, since T_{c} ∼ 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. 

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

Open with DEXTER 
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 nonsuperfluid case (striped lines) are shown. 

Open with DEXTER 
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 nonsuperfluid cases at low g_{2}. If entropy was not produced, this gap would vanish.
3.2. Dependence on equation of state
The Bose gas with twobody interactions is compared with threebody 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 twobody and threebody 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 twobody and threebody interactions as function of the interaction strengths g_{2} and g_{3}, respectively, with m = 30 eV, L = 100 kpc, and T/T_{c} = 0.1. The threebody interaction is multiplied by to make the comparison clearer. Both the superfluid case (solid lines) and the nonsuperfluid (striped lines) are shown. 

Open with DEXTER 
3.3. Effect of smallscale 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 smallscale 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 largescale 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 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/T_{c} = 0.1, and L = 100 kpc. Both the superfluid case (solid lines) and the nonsuperfluid (striped line) are shown. 

Open with DEXTER 
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 supercritical 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/T_{c} 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 10^{5}. 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. 

Open with DEXTER 
Fig. 9.
Profiles of a collapsing SFDM halo where no entropy is produced as w = v_{v} with an initial Gaussian density contrast, m = 30 eV, g = 10^{−5} eV^{−2}, L = 100 kpc, and T/T_{c} = 0.1. 

Open with DEXTER 
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 selfinteraction constraints
The distribution of DM, gas, and stellar mass in cluster collisions provides constraints on the crosssection of DM selfinteractions, σ/m < 0.5 cm^{2} g^{−1} (Harvey et al. 2015). In terms of the twobody interaction strength, this corresponds to (Pitaevskii & Stringari 2016)
The values of g_{2} 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_{2} 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/T_{c} = 2 × 10^{−6} reaches unity for various masses and interaction strengths. Both the superfluid case (solid lines) and the corresponding nonsuperfluid 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 selfinteraction. 

Open with DEXTER 
4. Conclusions
When superfluid behavior is included in a finitetemperature 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 zerotemperature pressure depends on the number density through n^{2} and n^{3}, 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. Largescale and highresolution 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 firstorder FORCE scheme used in this work that can capture the smallscale motion of the superfluid and its effect on structure formation.
The temperature and entropy in supercomoving variables are not given in Martel & Shapiro (1998, MS). We define them here as and , where T_{*} is a free parameter, , with ρ_{*} and v_{*} given in MS.
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
 Andersen, J. O. 2004, Rev. Mod. Phys., 76, 599 [CrossRef] [Google Scholar]
 Angus, G., Diaferio, A., Famaey, B., Gentile, G., & van der Heyden, K. 2014, J. Cosmol. Astropart. Phys., 2014, 079 [CrossRef] [Google Scholar]
 Angus, G. W., Diaferio, A., Famaey, B., & van der Heyden, K. J. 2013, MNRAS, 436, 202 [NASA ADS] [CrossRef] [Google Scholar]
 Barenghi, C. F., Skrbek, L., & Sreenivasan, K. R. 2014, Proc. Nat. Acad. Sci., 111, 4647 [CrossRef] [Google Scholar]
 Berezhiani, L., & Khoury, J. 2015, Phys. Rev. D, 92, 103510 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., & BoylanKolchin, M. 2017, ARA&A, 55, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S., Hoyos, C., & Oz, Y. 2014, J. High Energy Phys., 2014, 27 [CrossRef] [Google Scholar]
 Cyburt, R. H., Fields, B. D., Olive, K. A., & Yeh, T.H. 2016, Rev. Mod. Phys., 88, 015004 [NASA ADS] [CrossRef] [Google Scholar]
 Darve, C., Bottura, L., Patankar, N. A., & Van Sciver, S. 2012, AIP Conf. Proc., 1434, 247 [CrossRef] [Google Scholar]
 Del Popolo, A., & Le Delliou, M. 2017, Galaxies, 5, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Dodelson, S. 2011, Int. J. Mod. Phys. D, 20, 2749 [NASA ADS] [CrossRef] [Google Scholar]
 Doi, D., Shirai, Y., & Shiotsu, M. 2008, AIP Conf. Proc., 985, 648 [CrossRef] [Google Scholar]
 Elbert, O. D., Bullock, J. S., GarrisonKimmel, S., et al. 2015, MNRAS, 453, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Living Rev. Relativ., 15, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Glyde, H. R. 2013, J. Low Temp. Phys., 172, 364 [CrossRef] [Google Scholar]
 Harvey, D., Massey, R., Kitching, T., Taylor, A., & Tittley, E. 2015, Science, 347, 1462 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Khalatnikov, I. M. 2000, An Introduction to the Theory Of Superfluidity, 1st edn. (Boulder, CO: Westview Press) [Google Scholar]
 Khoury, J. 2016, Phys. Rev. D, 93, 103533 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L. 1941, Phys. Rev., 60, 356 [CrossRef] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1987, Course of Theoretical Physics, 2nd edn. (Oxford, UK: ButterworthHeinemann), 6 [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 816, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Martel, H., & Shapiro, P. R. 1998, MNRAS, 297, 467 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2005, ApJ, 632, 859 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2012, ApJ, 143, 40 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Milgrom, M. 1983a, ApJ, 270, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983b, ApJ, 270, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983c, ApJ, 270, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Mocz, P., Vogelsberger, M., Robles, V. H., et al. 2017, MNRAS, 471, 4559 [NASA ADS] [CrossRef] [Google Scholar]
 Pethick, C. J., & Smith, H. 2008, BoseEinstein Condensation in Dilute Gases, 2nd edn. (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Pitaevskii, L. P., & Stringari, S. 2016, BoseEinstein Condensation and Superfluidity (Oxford, UK: Oxford University Press) [CrossRef] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sales, L. V., Navarro, J. F., Oman, K., et al. 2016, MNRAS, 464, 2419 [NASA ADS] [CrossRef] [Google Scholar]
 SantosSantos, I. M., Brook, C. B., Stinson, G., et al. 2015, MNRAS, 455, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931 [NASA ADS] [CrossRef] [Google Scholar]
 Schive, H. Y., Chiueh, T., & Broadhurst, T. 2014, Nat. Phys., 10, 496 [CrossRef] [Google Scholar]
 Schwabe, B., Niemeyer, J. C., & Engels, J. F. 2016, Phys. Rev. D, 94, 043513 [CrossRef] [Google Scholar]
 Sharma, A., Khoury, J., & Lubensky, T. 2019, J. Cosmol. Astropart. Phys., 2019, 054 [CrossRef] [Google Scholar]
 Skrbek, L. 2011, J. Phys. Conf. Ser., 318, 012004 [CrossRef] [Google Scholar]
 Skrbek, L., & Sreenivasan, K. R. 2012, Phys. Fluids, 24, 011301 [CrossRef] [Google Scholar]
 Soulaine, C., Quintard, M., Baudouy, B., & Van Weelderen, R. 2017, Phys. Rev. Lett., 118, 074506 [CrossRef] [Google Scholar]
 Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Taylor, E., & Griffin, A. 2005, Phys. Rev. A, 72, 8739 [CrossRef] [Google Scholar]
 Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, ApJ, 606, 702 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. 2006, Appl. Numer. Math., 56, 1464 [CrossRef] [Google Scholar]
 Tulin, S., & Yu, H.B. 2018, Phys. Rep., 730, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Q., Marinacci, F., Maji, M., et al. 2016, MNRAS, 458, 1559 [NASA ADS] [CrossRef] [Google Scholar]
 Zuntz, J., Zlosnik, T. G., Bourliot, F., Ferreira, P. G., & Starkman, G. D. 2010, Phys. Rev. D, 81, 104015 [NASA ADS] [CrossRef] [Google Scholar]
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 twobody and threebody 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 zerotemperature 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;
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;
For the other thermodynamic quantities, such as pressure, entropy, etc., we must consider them above and below T_{c} separately. Both twobody and threebody interactions are given, parameterized by g_{2} and g_{3}, respectively.
A.1. T > T_{c}
The pressure is given by
where Γ(x) is the gamma function, Li_{z}(x) is the polylogarithmic function, β = 1/k_{B}T, and the chemical potential μ is determined by the equation for the number density
The entropy is
and the internal energy is
The sound speed used when determining the timestepping in the numerical scheme was
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 zerotemperature contributions due to interactions;
and the internal energy is again given by Eq. (A.6). The fastest sound speed was approximated using
and the critical velocity given by
There is a small discontinuity at the critical temperature, with μ = 2ng_{2} above and μ = ng_{2} below for the twobody interaction (and a similar jump in zerotemperature 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 firstorder 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 timestep: 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 timestep 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 onedimensional. 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 firstorder. There is a wellknown 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. Firstorder 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 finitevolume elements, or cells, and the Riemann problem at each cell interface is solved. The Riemann problem is the initial value problem with two piecewise 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 mcomponent state vector U that obeys the onedimensional conservative equation
where F is the flux. By integrating over the time interval [t^{n}, t^{n + 1}] and cellvolume [x_{i − 1/2}, x_{i + 1/2}], we get
where
In the firstorder Godunov scheme, the state U is assumed to be piecewise constant in each cell, given by the cell average . To compute F_{i + 1/2,} the states on the left and right sides of the interface is used, and , and the corresponding Riemann problem is solved or approximated. The timestep is chosen so that no signal in the domain travels further than one cell length Δx. This is given by a CourantFriedrichLewy (CFL) type condition
where v_{max} is the maximum signal speed in the domain, and is a number less than one that controls how far across a cell the fastest signal is allowed to move during each timestep. In simulations with gravity and expansion, additional constraints need to be added to the timestepping. For gravity, the freefall distance in each cell, with acceleration g, must be smaller than the cell lengths,
and for expansion, the relative change in the scale factor is restricted:
Here, and are also numbers less than one. In this work, we used , , and . The final value for the timestep is the smallest of the above,
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 twostep Lax–Wendroff flux;
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 supercomoving variables appear as source terms S in the superfluid equations. Continuing with the above example, we have
To modify our Godunov scheme to incorporate the sources in the flux, we did the following: at the beginning of each timestep, we had the states . 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 timestep of the source due to gravity,
where and are the left and right values for the sources. In this work, we computed these using the average gravitational acceleration , and the left and right states U_{i + 1/2, L/R}. We then used and as the input states in the Godunov scheme to get F_{i + 1/2}, and updated the state vectors from the previous timestep:
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;
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, , , S^{l}, etc. If , the Landau criterion is satisfied and we do nothing. If instead , we apply a small change to to update it to l + 1 and decrease w,
By keeping j constant and assuming that the change in the superfluid fraction is negligible compared to the change in velocity, we get
Using conservation of energy the change in internal energy is equal to the change in kinetic energy;
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:
The updated entropy is
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 , which was chosen as
where
The numerical factor in Eq. (B.20) was tuned to give as smooth wprofile as possible while keeping the scheme from becoming too slow.
All Figures
Fig. 1.
Superfluid fraction f_{s} for the twobody interacting Bose gas calculated by Sharma et al. (2019) with ρ = 10^{8} ρ_{c0} and m = 1 eV compared to the approximation f_{s, ap.} = 1 − (T/T_{c})^{3/2}. For sufficiently weak interactions, f_{s} can be approximated by the condensate fraction in an ideal Bose gas. 

Open with DEXTER  
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/T_{c} = 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/T_{c}, and hence the superfluid fraction, since f_{s} = ρ_{s}/ρ = 1 − (T/T_{c})^{3/2}. As the critical velocity is reached, entropy is generated, and T/T_{c} increases. 

Open with DEXTER  
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/T_{c} = 0.1 and L = 100 kpc. Both the superfluid case (solid lines) and the corresponding nonsuperfluid case (striped lines) are shown. For constant T/T_{c}, the temperature is increased for decreasing mass, since T_{c} ∼ 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. 

Open with DEXTER  
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/T_{c} = 0.1 and m = 50 eV. Both the superfluid case (solid lines) and the corresponding nonsuperfluid case (striped lines) are shown. 

Open with DEXTER  
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 nonsuperfluid case (striped lines) are shown. 

Open with DEXTER  
In the text 
Fig. 6.
Comparison of redshifts when the central density contrast reaches 200 for twobody and threebody interactions as function of the interaction strengths g_{2} and g_{3}, respectively, with m = 30 eV, L = 100 kpc, and T/T_{c} = 0.1. The threebody interaction is multiplied by to make the comparison clearer. Both the superfluid case (solid lines) and the nonsuperfluid (striped lines) are shown. 

Open with DEXTER  
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/T_{c} = 0.1, and L = 100 kpc. Both the superfluid case (solid lines) and the nonsuperfluid (striped line) are shown. 

Open with DEXTER  
In the text 
Fig. 8.
Evolution of T/T_{c} 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 10^{5}. 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. 

Open with DEXTER  
In the text 
Fig. 9.
Profiles of a collapsing SFDM halo where no entropy is produced as w = v_{v} with an initial Gaussian density contrast, m = 30 eV, g = 10^{−5} eV^{−2}, L = 100 kpc, and T/T_{c} = 0.1. 

Open with DEXTER  
In the text 
Fig. 10.
Redshifts when the linear density contrast for the mode k = 2/100 kpc^{−1} with T/T_{c} = 2 × 10^{−6} reaches unity for various masses and interaction strengths. Both the superfluid case (solid lines) and the corresponding nonsuperfluid 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 selfinteraction. 

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