Dark matter superfluidity and galactic dynamics

Abstract

We propose a unified framework that reconciles the stunning success of MOND on galactic scales with the triumph of the ΛCDM model on cosmological scales. This is achieved through the physics of superfluidity. Dark matter consists of self-interacting axion-like particles that thermalize and condense to form a superfluid in galaxies, with about milliKelvin critical temperature. The superfluid phonons mediate a MOND acceleration on baryonic matter. Our framework naturally distinguishes between galaxies (where MOND is successful) and galaxy clusters (where MOND is not): dark matter has a higher temperature in clusters, and hence is in a mixture of superfluid and normal phase. The rich and well-studied physics of superfluidity leads to a number of striking observational signatures.

Introduction

The standard Λ Cold Dark Matter (ΛCDM) model does very well at fitting large scale observables. On galactic scales, however, a number of challenges have emerged. Disc galaxies display a tight correlation between total baryonic mass and asymptotic velocity, the baryonic mass going as the fourth power of the circular velocity, known as the Baryonic Tully–Fisher Relation (BTFR). Hydrodynamical simulations can reproduce the BTFR by tuning baryonic feedback processes, but their stochastic nature naturally results in a much larger scatter. Furthermore, the mass and phase-space distributions of dwarf satellites in the Local Group are puzzling.

A radical alternative is MOdified Newtonian Dynamics (MOND), which replaces dark matter (DM) with a modification of gravity at low acceleration: the acceleration equals the Newtonian value when that value is large compared with a critical scale, and equals the geometric mean of the Newtonian value and the critical scale when the Newtonian value is small compared with it, with best-fit value for the critical acceleration of about 1.2 × 10⁻⁸ cm per second squared. This empirical force law has been remarkably successful at explaining a wide range of galactic phenomena. In the MOND regime, a test particle orbits an isolated source in such a way that the square of its velocity divided by the radius equals the square root of Newton’s constant times the baryonic mass times the critical acceleration, divided by the radius squared. This gives a constant asymptotic velocity, with the square of the circular velocity equal to the square root of Newton’s constant times the baryonic mass times the critical acceleration, which in turn implies the BTFR.

The empirical success of MOND, however, is limited to galaxies. The predicted temperature profile in galaxy clusters conflicts with observations. The Tensor–Vector–Scalar (TeVeS) relativistic extension fails to reproduce the CMB and matter spectra. The lensing features of merging clusters are problematic. This has motivated various hybrid proposals that include both DM and MOND.

In this Letter, together with a longer companion paper, we propose a novel framework that unifies the DM and MOND phenomena through the physics of superfluidity. There are two central ideas underlying our work, which must be carefully distinguished. The first is the very general idea that DM forms a superfluid inside galaxies, with a coherence length of order the size of galaxies. The critical temperature is about a milliKelvin, which intriguingly is comparable to Bose–Einstein condensation (BEC) critical temperatures for cold atom gases. Indeed, in many ways our DM behaves like cold dark atoms. The generic idea of DM superfluidity leads to a number of remarkable observational consequences. The superfluid nature of DM dramatically changes its macroscopic behavior in galaxies. Instead of evolving as independent particles, DM is more aptly described as collective excitations. The second central idea is the postulate that superfluid phonons mediate a MOND-like force between baryons. Since superfluidity only occurs at low temperature, our framework naturally distinguishes between galaxies (where MOND is successful) and galaxy clusters (where MOND is not). Due to the larger velocity dispersion in clusters, DM has a higher temperature and hence is in a mixture of superfluid and normal phases.

The superfluid interpretation makes the non-analytic nature of the MOND scalar action more palatable. The Unitary Fermi Gas, which has attracted much excitement in cold atom physics, is also governed by a non-analytic kinetic term. Our equation of state, with pressure proportional to the cube of the density, suggests that the DM superfluid arises through three-body interactions. It would be fascinating to find precise cold atom systems with the same equation of state, as this would give important insights on the microphysical interactions underlying our superfluid. Tantalizingly, this might allow laboratory simulations of galactic dynamics.

The idea of DM BEC has been studied before, with important differences from our work. In BEC DM galactic dynamics are caused by the condensate density profile; in our case phonons play a key role in explaining the BTFR. Moreover, BEC DM has pressure proportional to the square of the density instead of the cube. This implies a much lower sound speed, which puts BEC DM in tension with observations.

DM condensation

In order for DM particles to condense in galaxies, their de Broglie wavelength, of order the inverse of the particle mass times its velocity, must be larger than the interparticle separation, of order the cube root of the mass divided by the virial density. From standard collapse theory, the density at virialization is about 5.4 × 10⁻²⁸ g/cm³ times the cube of one plus the virialization redshift, while the virial velocity is 113 km/s times the cube root of the halo mass in units of 10¹² solar masses times the square root of one plus the virialization redshift.

Equation 1. The de Broglie condition then implies that the particle mass must be smaller than about 2.3 eV times the three-eighths power of one plus the virialization redshift, divided by the fourth root of the halo mass in units of 10¹² solar masses.

We work in units where h-bar equals one. The second condition is that DM thermalizes, with temperature set by the virial velocity. The interaction rate is the number density times the velocity times the cross section per unit mass times the Bose enhancement factor, the latter being the phase-space density expressed through the virial density, the mass and the velocity. The rate should be larger than the inverse dynamical time, so that the coherence length will span the halo. This translates into a bound on the cross section.

Equation 2. The cross section per unit mass must be larger than about 52 cm²/g times the fourth power of the mass in electronvolts times the two-thirds power of the halo mass in units of 10¹² solar masses, divided by the seven-halves power of one plus the virialization redshift.

Later on, we will adopt m = 0.6 eV as a fiducial value. For a halo of 10¹² solar masses virialized at redshift 2, the inequality becomes a cross section per unit mass larger than about 0.1 cm²/g. The lower end is consistent with current constraints on the cross section per unit mass for self-interacting dark matter (SIDM), though these constraints must be carefully revisited in the superfluid context.

Equation 3. The critical temperature, obtained by equipartition — Boltzmann’s constant times the critical temperature equals one third of the mass times the square of the critical velocity — is in the milliKelvin range: about 6.5 mK times the square of one plus the virialization redshift, divided by the five-thirds power of the mass in electronvolts.

For temperatures between zero and the critical temperature, the system is a mixture of condensate and normal components. The fraction of condensed particles, one minus the three-halves power of the ratio of the temperature to the critical temperature, is shown in Fig. 1 as a function of halo mass assuming virialization at redshift zero. For a mass of about an electronvolt, galaxies are almost completely condensed while massive clusters have a significant normal component. Also, since the critical temperature depends on redshift, halos at higher redshift tend to have more superfluid than the ones formed more recently.

Fig. 1. Fraction of DM particles in the condensate.

Superfluid phase

The relevant low-energy degrees of freedom of a superfluid are phonons, described by a scalar field θ. In the presence of a gravitational potential, the non-relativistic effective action is the pressure evaluated on the kinetic variable X, where X is the time derivative of θ minus the mass times the potential, minus the square of the gradient of θ divided by twice the mass. The nature of the superfluid — that is, its equation of state — is encoded in the choice of the pressure function. Up to this point, the discussion has been very general. However, in order to endow our superfluid with MOND-like phenomenology, we conjecture that DM superfluid phonons are governed by the MOND action.

Equation 4. The Lagrangian is two thirds of Λ times the three-halves power of twice the mass, times X multiplied by the square root of the magnitude of X, minus α times θ times the baryonic density divided by the Planck mass. Here Λ is a mass scale, and α is a dimensionless constant.

This action should only be trusted away from X = 0, as we will see later. The matter coupling breaks the shift symmetry at the level of one over the Planck mass and is thus technically natural. Remarkably, Equation 4 is strikingly reminiscent of the Unitary Fermi Gas, whose Lagrangian goes as the five-halves power of X, which is also non-analytic.

The phonon action uniquely fixes the properties of the condensate through standard thermodynamics. At finite chemical potential, and ignoring the gravitational potential, the pressure is given by the Lagrangian density.

Equation 5. The pressure as a function of the chemical potential is two thirds of Λ times the three-halves power of twice the mass times the chemical potential.

This is the grand canonical equation of state for the condensate. The number density, obtained as the derivative of the pressure with respect to the chemical potential, is Λ times the three-halves power of twice the mass times the square root of the chemical potential.

Equation 7. Combining these expressions with the mass density equal to the mass times the number density, we obtain a pressure equal to the cube of the density divided by twelve times the square of Λ times the sixth power of the mass.

This is a polytropic equation of state with index one half. In comparison, BEC DM has pressure proportional to the square of the density.

Including phonon excitations, the quadratic action for the fluctuation has a sound speed that can be immediately read off.

Equation 8. The sound speed is the square root of twice the chemical potential divided by the mass.

Using the equation of state, we compute the static, spherically-symmetric density profile of the DM condensate halo. Introducing dimensionless variables for the density and the radius, with a central density, hydrostatic equilibrium implies the Lane–Emden equation with index one half, with boundary conditions that the dimensionless density is one at the origin and its derivative vanishes there. The numerical solution is shown in Fig. 2. It vanishes at a dimensionless radius of about 2.75, which defines the halo size. Meanwhile the central density is related to the halo mass in the standard Lane–Emden way. Combining these results, we obtain the following.

Equation 9. The central density is about 7 × 10⁻²⁵ g/cm³ times the two-fifths power of the halo mass in units of 10¹² solar masses, times the eighteen-fifths power of the mass in electronvolts, times the six-fifths power of Λ in milli-electronvolts; and the halo radius is about 36 kpc times the one-fifth power of the halo mass, divided by the six-fifths power of the mass in electronvolts and the two-fifths power of Λ in milli-electronvolts.

Remarkably, with a mass of order an electronvolt and Λ of order a milli-electronvolt we obtain DM halos of realistic size. Concretely, as fiducial values we will fix m = 0.6 eV and Λ = 0.2 meV. This implies a halo radius of about 125 kpc for a dark matter mass of 10¹² solar masses.

Fig. 2. Numerical solution of the Lane–Emden equation.

Figure 2 represents the density profile inside the superfluid core of a galaxy, to the extent that baryons can be neglected. The profile is cored — the density varies slowly, changing by order unity through most of the volume. In contrast, the density in a Navarro–Frenk–White (NFW) halo would vary more rapidly (as an approximate power law) over the same volume. In general, we expect the superfluid core to be surrounded by a cloud of DM particles in the normal phase, likely described by an NFW profile.

Including baryons

We now derive the phonon profile in galaxies, including baryons. For this purpose, we model baryons as a static, spherically-symmetric source. The assumption of spherical symmetry simplifies the treatment and becomes increasingly valid far from the source. The equation of motion can be readily integrated.

Equation 11. The divergence of the square root of twice the mass times the magnitude of X, multiplied by the gradient of the phonon field, equals α times the enclosed baryonic mass divided by eight pi times the Planck mass times the radius squared — a quantity we call κ(r).

There are two branches of solutions, depending on the sign of X. We focus on the MOND branch, the one with X negative, where the gradient of the phonon field is fixed by κ and by the shifted chemical potential. In the limit where κ divided by the mass is much larger than the shifted chemical potential, the gradient of the phonon field approaches the square root of κ divided by the mass. In this limit the scalar acceleration on a baryonic particle is α times that gradient divided by the Planck mass. Matching to the MOND result — the acceleration being the geometric mean of the critical acceleration and the Newtonian acceleration — fixes α in terms of Λ and the critical acceleration.

Equation 14. The three-halves power of α times Λ equals the critical acceleration times the Planck mass, which is about 0.8 meV. For the fiducial value Λ = 0.2 meV, we get α about 2.5.

As it stands, the solution with X negative is unstable. It leads to unphysical halos, with growing DM density profiles. The instability can be seen by expanding the action in the fluctuation about the background profile: the kinetic term is ghostly when the background X is negative. (The branch with X positive is stable but does not admit a MOND regime.)

Since DM in actual galactic halos is not at zero temperature, however, we expect the action to receive finite-temperature corrections in galaxies. At finite sub-critical temperature, the Lagrangian depends on X and two additional scalars: the square root of the determinant built from the gradients of the Lagrangian coordinates of the normal fluid, and a second scalar built from the shifted chemical potential, the time derivative of the phonon field, and the normal-fluid velocity contracted with the phonon gradient.

Equation 15. For instance, consider the two-derivative operator equal to a temperature-dependent mass squared times the square of that second scalar, which in the rest frame of the normal fluid reduces to the mass squared times the square of the shifted chemical potential plus the time derivative of the phonon field.

This leaves the static profile unchanged, but modifies the quadratic Lagrangian, restoring stability for sufficiently large values of that mass. Specifically this is the case for a mass of order an electronvolt or greater, remarkably of the same order as the particle mass. By the same token it corrects the condensate pressure by an amount equal to the mass squared times the chemical potential squared, which obliterates the unwanted growth in the DM density profile, resulting in localized, finite-mass halos.

Equation 16. Another possibility is the finite-temperature Lagrangian: two thirds of Λ times the three-halves power of twice the mass, times the square root of X, times the magnitude of X minus a temperature-dependent coefficient β times the second scalar, where β must be greater than three halves, as required for stability.

The DM condensate pressure is then identical to the earlier result modulo the replacement of Λ by Λ times the square root of β minus one. The halo density profile is therefore identical. Including baryons, it is easy to show that the integrated equation of motion picks up corresponding β-dependent factors.

At small distances the solution approximates the MOND profile, with the phonon gradient going as the inverse radius, and this gives the dominant force on a test baryonic particle. At large distances the solution tends to a different regime, but there it is subdominant to the gravitational acceleration due to the DM halo. The transition radius delineating the MOND regime from the DM-halo regime occurs when κ becomes comparable to the mass times the shifted chemical potential.

Equation 18. Substituting and approximating the shifted chemical potential by its central value, the transition radius is about 28 kpc times the one-tenth power of the baryonic mass in units of 10¹¹ solar masses, times the two-fifths power of the baryon-to-dark-matter mass ratio, divided by the eight-fifths power of the mass in electronvolts and the eight-fifteenths power of Λ in milli-electronvolts.

For a galaxy with a baryonic mass of 3 × 10¹¹ solar masses and the cosmic dark-matter-to-baryon ratio of about 6, the MOND regime extends to about 70 kpc.

Complex scalar field description

It is well-known that a superfluid is described in the weak-coupling regime as a self-interacting complex scalar field with global U(1) symmetry. Consider the relativistic field theory whose Lagrangian is minus the squared gradient of the complex field, minus the mass squared times the squared magnitude of the field, minus a self-interaction term built from the squared gradient and the squared magnitude divided by a denominator containing a scale Λ_c, which makes the field vanishing point well-defined. Substituting the field written in modulus-and-phase form and taking the non-relativistic limit, we obtain the corresponding non-relativistic Lagrangian.

To leading order in gradients we ignore the squared-gradient contributions of the modulus and integrate out the modulus. In the limit where the modulus is much larger than Λ_c this gives the squared modulus equal to Λ times the square root of twice the mass times the magnitude of X. Substituting this back, we recover the MOND action — given by the kinetic part of Equation 4 — in that limit.

A finite Λ_c implies that MOND is restricted to values of the modulus above Λ_c, that is, to scalar accelerations above a floor set by the ratio of Λ_c to Λ times the critical acceleration divided by α squared. Thus the MOND regime does not apply to arbitrarily small accelerations. By choosing Λ_c a factor of a few smaller than Λ, the predicted departure from MOND can occur around the acceleration scale of the Milky Way dwarf spheroidals, which are well-known to pose a challenge for MOND.

Cosmology

Given their mass, our DM particles are axion-like. The simplest genesis scenario is through vacuum displacement, with DM being generated at a time when the Hubble rate is of order the mass, of order an electronvolt. This corresponds to a photon temperature of order a TeV, that is, around the weak scale. The DM rapidly reaches thermal equilibrium with itself and becomes superfluid, but is decoupled from ordinary matter.

Naturally DM is much colder cosmologically than in collapsed structures. The ratio of the temperature to the critical temperature is constant cosmologically and can be evaluated at matter-radiation equality; the result is about 10⁻²⁸, which is much colder than the typical range of 10⁻⁶ to 10⁻² found in galaxies. We will see shortly that, in order to obtain an acceptable cosmology, we need Λ and α to assume different values cosmologically. To avoid confusion, in what follows we will denote their cosmological values with a nought subscript.

Equation 22. On a cosmological background, the phonon equation of motion derived from Equation 4 sets the time derivative of the cosmological Λ times the three-halves power of twice the mass times the cube of the scale factor times the square root of the phonon time derivative equal to minus the cosmological α times the cube of the scale factor times the baryonic density, divided by the Planck mass.

Since the comoving baryonic density is constant, this can be integrated straightforwardly.

Equation 23. The resulting non-relativistic energy density is minus the cosmological α times the cosmological Λ times the baryon density times the mass times the age, divided by the Planck mass, plus a dust term that redshifts as the inverse cube of the scale factor.

In the matter-dominated era, the baryonic contribution redshifts more slowly than the dust term. In order for the superfluid to behave as ordinary dust, the second term should dominate over the first all the way to the present time. Substituting the age of the universe, 13.9 billion years, and assuming a dark-matter-to-baryon ratio of 6, we obtain a bound on the cosmological α: it must be smaller than about 2.4 × 10⁻⁵ times the cosmological Λ in electronvolts, for a mass of order an electronvolt.

Since the phonon time derivative scales as the inverse cube of the scale factor, the phonon velocity increases with redshift and inevitably results in a breakdown of the non-relativistic approximation. In other words, the superfluid becomes relativistic at sufficiently high density. This is consistent with the equation of state: the equation-of-state parameter is the density squared divided by twelve times the cosmological Λ squared times the sixth power of the mass. Demanding that it be non-relativistic by matter-radiation equality puts a lower bound on the cosmological Λ of about 0.1 eV.

This is roughly four orders of magnitude larger than the fiducial value of 0.2 meV assumed in galaxies. This can be achieved, for instance, if Λ depends on temperature through a factor of one plus a large constant times the fourth root of the temperature ratio, with that constant of order 10⁴. Meanwhile, the earlier bound implies a cosmological α smaller than about 10⁻⁴, which is roughly four orders of magnitude smaller than the order-unity value obtained in galaxies. This can be achieved, for instance, by a similar temperature dependence for α. Note that the scale appearing in the phonon-baryon coupling in Equation 4 is nearly temperature-independent.

Gravitational lensing

In TeVeS the complete absence of DM requires introducing a time-like vector field, as well as a complicated coupling between the scalar, the vector and baryons in order to reproduce lensing observations. In our case, there is no need to introduce an extra vector, as the normal fluid already provides us with a time-like vector. Moreover, our DM contributes to lensing, so we are free to generalize the TeVeS coupling.

Equation 26. Suppose matter fields couple to an effective metric equal to the spacetime metric minus twice α times the phonon field divided by the Planck mass, multiplied by a combination of the metric weighted by γ and the normal-fluid four-velocity product weighted by one plus γ, with γ equal to one corresponding to TeVeS.

In the weak-field limit, this gives metric perturbations in which the time-time component carries the sum of the Newtonian potential and the phonon term, while the spatial components carry the Newtonian potential and γ times the phonon term, the Newtonian potential itself being sourced by both the baryonic and the dark matter density. Hence the lensing signal arises from a combination of the disformal coupling to the normal-fluid velocity and the DM condensate density profile. Determining the allowed range of γ will require a detailed study, which is beyond the scope of this paper.

Observational implications

We conclude with some astrophysical implications of our DM superfluid.

Vortices. When spun faster than a critical velocity, a superfluid develops vortices. The typical angular velocity of halos is well above critical, giving rise to an array of DM vortices permeating the disc. It will be interesting to see whether these vortices can be detected through substructure lensing, for example with ALMA.

Galaxy mergers. A key difference with CDM is the merger rate of galaxies. Applying Landau’s criterion, we find two possible outcomes. If the infall velocity is less than the phonon sound speed (of order the virial velocity), then halos will pass through each other with negligible dissipation, resulting in multiple encounters and a longer merger time. If the infall velocity is comparable to or greater than the sound speed, however, the encounter will excite DM particles out of the condensate, resulting in dynamical friction and rapid merger.

Bullet cluster. For merging galaxy clusters, the outcome also depends on the relative fraction of superfluid versus normal components in the clusters. For subsonic mergers, the superfluid cores should pass through each other with negligible friction (consistent with the Bullet Cluster), while the normal components should be slowed down by self interactions. Remarkably this picture is consistent with the lensing map of the Abell 520 “train wreck”, which shows lensing peaks coincident with galaxies (superfluid components), as well as peaks coincident with the X-ray luminosity peaks (normal components).

Dark-bright solitons. Galaxies in the process of merging should exhibit interference patterns (so-called dark-bright solitons) that have been observed in BECs counterflowing at super-critical velocities. This can potentially offer an alternative mechanism to generate the spectacular shells seen around elliptical galaxies.

Acknowledgements

We thank A. Arvanitaki, L. Blanchet, A. Erickcek, B. Famaey, L. Hui, B. Jain, R. Kamien, A. Kosowsky, T. Lubensky, S. McGaugh, A. Nicolis, M. Pawlowski, J. Peebles, R. Sheth, D. Spergel, P. Steinhardt and M. Zaldarriaga. J.K. is supported by NSF CAREER Award PHY-1145525 and NASA ATP grant NNX11AI95G. L.B. is supported by funds provided by the University of Pennsylvania.