Lasha Berezhiani

Benoit Famaey

Justin Khoury

Phenomenological consequences of superfluid dark matter with baryon-phonon coupling

Keywords: numbers: 98.10.+z, 98.62.Dm, 95.35.+d, 95.30.Sf

published or not. The documents may come

Phenomenological consequences of superfluid dark matter with baryon-phonon coupling

Lasha Berezhiani, Benoit Famaey, Justin

I. INTRODUCTION

The standard Λ Cold Dark Matter (ΛCDM) cosmological model provides an excellent fit to the Cosmic Microwave Background (CMB) angular power spectrum, to the large-scale structure matter power spectrum, to the Hubble diagram of Type Ia supernovae, and to the detailed scale of baryonic acoustic oscillations. In the context of simulations of structure formation, it has proven successful down to the scale of galaxy clusters and groups.

On galactic scales, however, the ΛCDM model faces some challenges [START_REF] Bullock | Small-Scale Challenges to the ΛCDM Paradigm[END_REF]. These include the internal structure of dwarf galaxies in the Local Group [START_REF] Boylan-Kolchin | Too big to fail? The puzzling darkness of massive Milky Way subhaloes[END_REF][START_REF] Boylan-Kolchin | The Milky Way's bright satellites as an apparent failure of LCDM[END_REF] and possibly beyond [START_REF] Papastergis | Is there a ?too big to fail? problem in the field?[END_REF], and the vast planar structures seen around the Milky Way and Andromeda [START_REF] Pawlowski | The VPOS: a vast polar structure of satellite galaxies, globular clusters and streams around the Milky Way[END_REF][START_REF] Pawlowski | Dwarf Galaxy Planes: the discovery of symmetric structures in the Local Group[END_REF][START_REF] Pawlowski | The rotationally stabilized VPOS and predicted proper motions of the Milky Way satellite galaxies[END_REF][START_REF] Ibata | A Vast Thin Plane of Co-rotating Dwarf Galaxies Orbiting the Andromeda Galaxy[END_REF]. Some other issues, notably those related to the frequency, size and pattern speeds of galactic bars, are related to dynamical friction between baryons and a live dark matter (DM) halo, e.g., [START_REF] Debattista | Dynamical friction and the distribution of dark matter in barred galaxies[END_REF][START_REF] Sellwood | Bar Instability in Disk-Halo Systems[END_REF].

Chief among these challenges is arguably the tight correlation between the distribution of baryons and the gravitational acceleration in galaxies [START_REF] Mcgaugh | Radial Acceleration Relation in Rotationally Supported Galaxies[END_REF][START_REF] Lelli | One Law to Rule Them All: The Radial Acceleration Relation of Galaxies[END_REF]. This apparent conspiracy, known as the mass discrepancy acceleration relation (MDAR), states that the acceleration expe-rienced by a baryonic particle (irrespective of whether it is actually due to DM, modified gravity, or both) can be uniquely predicted from the baryon density profile. At large distances, the MDAR implies the baryonic Tully-Fisher relation (BTFR) [START_REF] Mcgaugh | The Baryonic Tully-Fisher Relation of Gas Rich Galaxies as a Test of LCDM and MOND[END_REF], which relates the total baryonic mass to the asymptotic/flat rotation velocity as M b ∼ V 4 f . This relation holds over five decades in mass, with remarkably small scatter [START_REF] Mcgaugh | The Small Scatter of the Baryonic Tully-Fisher Relation[END_REF]. Another scaling relation is the correlation between the central stellar and dynamical surface densities in disk galaxies [START_REF] Lelli | The Relation between Stellar and Dynamical Surface Densities in the Central Regions of Disk Galaxies[END_REF].

All of these relations involve a universal acceleration constant a 0 ∼ 10 -8 cm/s 2 . This constant independently defines i) the zero-point of the BTFR; ii) the transition of the acceleration at which the mass discrepancy between baryonic and total mass appears; iii) the transition surface density between DM-dominated Low Surface Brightness (LSB) galaxies and baryon-dominated High Surface Brightness (HSB) galaxies; iv) the observed maximal baryonic surface density for disk galaxies [START_REF] Mcgaugh | Galaxy selection and the surface brightness distribution[END_REF]; and v) the DM halo central surface density [START_REF] Donato | A constant dark matter halo surface density in galaxies[END_REF].

A. Milgrom's law These independent occurrences of a 0 in observations can all be summarized by a simple empirical law proposed by Milgrom over thirty years ago [START_REF] Milgrom | A Modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis[END_REF][START_REF] Milgrom | A Modification of the Newtonian dynamics: Implications for galaxies[END_REF][START_REF] Milgrom | A modification of the Newtonian dynamics: implications for galaxy systems[END_REF]. Milgrom's law states that the total gravitational attraction a is approximately a N in the regime a N a 0 , and approaches the geometric mean √ a N a 0 whenever a N a 0 , where a N is the usual Newtonian gravitational field generated from the observed distribution of baryonic matter alone. This simple recipe successfully accounts for a wide range of galactic phenomena. It is remarkable that it was suggested by Milgrom well before being corroborated by detailed data: hence, Milgrom made succesful predictions for the behavior of, for instance, LSB galaxies which were not even known to exist at the time. Today, the detailed shape of rotation curves of disk galaxies with vastly different structural parameters is reproduced with uncanny precision by this empirical law. See [START_REF] Sanders | Modified Newtonian dynamics as an alternative to dark matter[END_REF][START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF] for comprehensive reviews.

Although historically proposed as an alternative to DM, hence the name Modified Newtonian Dynamics (MOND), the empirical power of Milgrom's formula at fitting galaxy properties is unequivocal, especially in rotationally supported systems. A debate exists, however, as to the possible origin of this law:

1. The conservative viewpoint is that Milgrom's law is an emergent phenomenon, linked to efficient feedback effects on galaxy scales. This possibility faces a number of challenges. Firstly, semi-empirical models with cored DM halos can reproduce the observed slope and normalization of the BTFR,

but not yet the small scatter [START_REF] Mcgaugh | The Small Scatter of the Baryonic Tully-Fisher Relation[END_REF][START_REF] Di Cintio | The mass discrepancy acceleration relation in a ΛCDM context[END_REF][START_REF] Desmond | The scatter, residual correlations and curvature of the SPARC baryonic Tully-Fisher relation[END_REF]. Secondly, the diversity of shapes of rotation curves at a given maximal circular velocity remains a puzzle in ΛCDM [START_REF] Oman | The unexpected diversity of dwarf galaxy rotation curves[END_REF] (though see also [START_REF] Read | Understanding the shape and diversity of dwarf galaxy rotation curves in ΛCDM[END_REF]). Finally, the observation that the central slope of the rotation curve correlates instead with the baryonic surface density means that feedback processes should be more efficient at removing DM from the central regions of galaxies with lower surface density, which are often more gas-dominated. Hence, the feedback efficiency would have to increase with decreasing star formation rate [START_REF] Mcgaugh | The Baryonic Tully-Fisher Relation of Gas Rich Galaxies as a Test of LCDM and MOND[END_REF], and be tightly anticorrelated with baryonic surface density (or correlated with scale-length at a given mass scale).

More generally, the feedback explanation could present a fine-tuning problem for galaxies that have had a complex formation history, as all effects not taken into account when dealing with internal feedback processes [START_REF] Read | Understanding the shape and diversity of dwarf galaxy rotation curves in ΛCDM[END_REF], such as the effect of environment and the merger history of each individual galaxy, must then be precisely tuned with the internal feedback in order to give precisely the observed relation [START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF].

2. At the other end of the spectrum is the viewpoint that Milgrom's law represents a fundamental modification of gravity (e.g., [START_REF] Bekenstein | Does the missing mass problem signal the breakdown of Newtonian gravity?[END_REF][START_REF] Bekenstein | Relativistic gravitation theory for the MOND paradigm[END_REF][START_REF] Milgrom | Bimetric MOND gravity[END_REF]). See [START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF][START_REF] Bruneton | Fieldtheoretical formulations of MOND-like gravity[END_REF] for reviews of such theories. By now this possibility seems however rather unlikely, given the observational evidence for DM behaving as a collisionless fluid on cosmological and cluster scales. Indeed, Milgrom's law does not work on galaxy cluster scales [START_REF] Gerbal | Physical and dynamical properties of two clusters of galaxies derived from X-ray data[END_REF][START_REF] Sanders | The Virial Discrepancy in Clusters of Galaxies in the Context of Modified Newtonian Dynamics[END_REF][START_REF] Clowe | Weak lensing mass reconstruction of the interacting cluster 1E0657-558: Direct evidence for the existence of dark matter[END_REF][START_REF] Clowe | A direct empirical proof of the existence of dark matter[END_REF][START_REF] Angus | On the Law of Gravity, the Mass of Neutrinos and the Proof of Dark Matter[END_REF][START_REF] Angus | X-ray group and cluster mass profiles in MOND: unexplained mass on the group scale[END_REF], and modified gravity theories cannot a priori explain the angular power spectrum of the CMB without additional DM or unrealistically massive ordinary neutrinos, with a mass above the terrestrial experimental constraints [START_REF] Skordis | Large Scale Structure in Bekenstein's theory of relativistic Modified Newtonian Dynamics[END_REF]. Such theories are also highly constrained in the Solar System [START_REF] Hees | Constraints on modified Newtonian dynamics theories from radio tracking data of the Cassini spacecraft[END_REF][START_REF] Hees | Combined Solar System and rotation curve constraints on MOND[END_REF] (though see [START_REF] Babichev | Improving relativistic MOND with Galileon k-mouflage[END_REF]). Moreover, Milgrom's law has some problems at subgalactic scales (e.g., [START_REF] Ibata | The globular cluster NGC 2419: a crucible for theories of gravity[END_REF]).

3. A middle-ground interpretation is that Milgrom's law is telling us something about the fundamental nature of DM. Indeed, while the existence of a new degree of freedom acting as a pressureless and dissipationless fluid is on firm grounds within the linear regime in cosmology, the exact nature of this dark fluid is far from known.

Already in the context of self-interacting DM (SIDM) [START_REF] Spergel | Observational evi-dence for selfinteracting cold dark matter[END_REF], some recent encouraging results have shown how underdense halos could indeed be associated with extended baryonic disks [START_REF] Kamada | How the Self-Interacting Dark Matter Model Explains the Diverse Galactic Rotation Curves[END_REF][START_REF] Creasey | Spreading out and staying sharp -creating diverse rotation curves via baryonic and selfinteraction effects[END_REF], in line with observations. However, getting a BTFR as tight as observed, as well as a MDAR at all radii remains challenging. A more radical possibility is Modified DM, inspired by gravitational thermodynamics [START_REF] Ho | Cold Dark Matter with MOND Scaling[END_REF][START_REF] Ho | Quantum Gravity and Dark Matter[END_REF][START_REF] Ho | Dark Matter, Infinite Statistics and Quantum Gravity[END_REF][START_REF] Edmonds | Modified Dark Matter: Relating Dark Energy, Dark Matter and Baryonic Matter[END_REF].

Our viewpoint here is that the correlation between DM and visible matter embodied in Milgrom's law could be the result of a novel interaction between DM and baryons [START_REF] Blanchet | Gravitational polarization and the phenomenology of MOND[END_REF][START_REF] Blanchet | Model of Dark Matter and Dark Energy Based on Gravitational Polarization[END_REF][START_REF] Khoury | An Alternative to Particle Dark Matter[END_REF][START_REF] Blanchet | Dark Matter via Massive (bi-)Gravity[END_REF][START_REF] Blanchet | Dipolar Dark Matter with Massive Bigravity[END_REF][START_REF] Salucci | Evidences for Collisional Dark Matter In Galaxies?[END_REF]. Indeed, this could naturally avoid the problems of pure-CDM and pure-MOND, while in principle also reproducing the observed phenomenology on galactic scales.

B. Superfluid dark matter

Recently, in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF], we proposed a novel theory of DM superfluidity, inspired by recent developments in cold atom physics. See also [START_REF] Khoury | Another Path for the Emergence of Modified Galactic Dynamics from Dark Matter Superfluidity[END_REF][START_REF] Khoury | Dark Matter Superfluidity[END_REF][START_REF] Hodson | Galaxy Clusters in the Context of Superfluid Dark Matter[END_REF]. The framework marries the phenomenological success of the ΛCDM model on cosmological scales with that of Milgrom's empirical law [START_REF] Milgrom | A Modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis[END_REF] at fitting galaxy properties. The DM and MOND components have a common origin, representing different phases of a single underlying substance. This is achieved through superfluidity.

The DM candidate proposed in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] consists of axion-like particles with sufficiently strong selfinteractions (σ/m ∼ > 0.1 cm 2 /g, where σ is the selfinteraction cross-section and m the particle mass) such As in General Relativity that they thermalize in galaxies. With m ∼ eV, their de Broglie wavelengths overlap in the (cold and dense enough) central region of galaxies, resulting in Bose-Einstein condensation into a superfluid phase.

The superfluid nature of DM dramatically changes its macroscopic behavior in galaxies. Instead of behaving as individual collisionless particles, the DM is more aptly described as collective excitations, in particular phonons. The DM phonons couple to ordinary matter, thereby mediating a long-range force (beyond Newtonian gravity) between baryons. In contrast with theories that propose to fundamentally modify Newtonian gravity, in this case the new long-range force mediated by phonons is an emergent property of the DM superfluid medium.

For a particular choice of the superfluid equation of state, the resulting phonon effective Lagrangian is similar to the Bekenstein-Milgrom MOND theory [START_REF] Bekenstein | Does the missing mass problem signal the breakdown of Newtonian gravity?[END_REF]. Remarkably, as reviewed below in Sec. II, the desired superfluid phonon effective theory is strikingly similar to that of the Unitary Fermi Gas [START_REF] Zwerger | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF][START_REF] Giorgini | Theory of ultracold atomic Fermi gases[END_REF], which has generated much excitement in the cold atom community in recent years.

The aim of this paper is to develop novel distinctive predictions of such a model, especially on galaxy scales. For this purpose we will derive a realistic superfluid halo density profile, improving on the original derivation of [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] by carefully modeling the superfluid at finite temperature. The density profile of [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] used the mean-field Gross-Pitaevskii approach [START_REF] Gross | Structure of a quantized vortex in boson systems[END_REF][START_REF] Pitaevskii | Vortex lines in an imperfect Bose gas[END_REF], which assumes that DM is entirely in the condensed phase. In reality, owing to their velocity dispersion DM particles have a small non-zero temperature in galaxies. As is familiar from liquid helium, a superfluid at finite (sub-critical) temperature is best described phenomenologically as a mixture of two fluids [64][65][66]: i) the superfluid, which by definition has vanishing viscosity; ii) the "normal" component, comprised of excitations, which is viscous and carries entropy.

The DM density profile should consist of a superfluid core, with approximately homogeneous density, sur-rounded by an envelope of DM particles in the normal phase following a collisionless Navarro-Frenk-White (NFW) density profile [START_REF] Navarro | A Universal density profile from hierarchical clustering[END_REF]. In the present contribution, the transition radius will be estimated in two ways. The first estimate, denoted by R T , is determined by the density dropping to a value where the interaction rate is too low to maintain local thermal equilibrium. The second estimate, denoted by R NFW , is the radius at which both density and pressure (or equivalently, velocity dispersion) can be matched continuously to an NFW profile. As we shall demonstrate hereafter, for the range of theory parameters that give reasonable fit to rotation curves, R T and R NFW coincide to within 30%.

A key difference with the original papers [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] is that we now choose parameters such that the superfluid core makes up only a fraction of the total mass/volume of the halo. For the representative massive disk galaxy described in Sec. VII, for instance, the superfluid core makes up 22% of the total DM mass. The superfluid core must be large enough to encompass the observed rotation curve, since we rely on the phonon-mediated force to reproduce the rotation curve phenomenology. But it is small enough that most of the DM halo mass is approximately in collisionless form and hence should be triaxial exactly as ΛCDM halos and consistent with observations [START_REF] Clampitt | Lensing measurements of the ellipticity of luminous red galaxies dark matter haloes[END_REF].

A by-product of having a relatively small superfluid core is that most of the gravitational lensing signal will come from the NFW envelope. As a corollary, to reproduce galaxy-galaxy lensing statistics we need not assume any non-minimal coupling between DM superfluid phonons and photons. Both photons and gravitons travel at the speed of light along the same geodesics, consistent with the tight constraint recently established by the neutron star merger GW170817 [START_REF] Abbott | GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral[END_REF]. In contrast, some DMfree relativistic completions of MOND, such as [START_REF] Bekenstein | Relativistic gravitation theory for the MOND paradigm[END_REF], rely on non-minimal photon couplings and are now ruled out by GW170817 [START_REF] Boran | GW170817 Falsifies Dark Matter Emulators[END_REF]. Potentially measurable effects of nonvanishing phonon sound-speed in the context of future gravitational wave experiments were discussed in [START_REF] Cai | Gravitational wave as probe of superfluid dark matter[END_REF].

After a few general constraints on the theory parameters in Secs. III and IV, the main focus of this paper will be the explicit rotation curve fitting of disk galaxies in Secs. V-VII. For this purpose, we will illustrate the method with two representative LSB and HSB disk galaxies (Sec. VII). Given the lack of symmetry, an exact solution for the superfluid/phonon profile would be a daunting task. Instead we will adopt the following hybrid approximation method.

In general the acceleration a on a test baryon particle receives three contributions within the superfluid region:

a = a b + a DM + a phonon . (1)
We calculate the actual Newtonian acceleration generated by the baryons, a actual b , exactly from the actual nonspherical baryon density. The second term is the gravitational acceleration from the superfluid core, which we derive from a spherically symmetric approximation. The phonon-mediated acceleration a phonon is sourced by the actual non-spherical baryon density distribution, but also depends indirectly on the Newtonian gravitational potential Φ, which we approximate as spherically symmetric for this derivation. As we will show in Sec. VII, the superfluid model offers a satisfactory fit of the rotation curve in both the LSB and HSB cases.

In Secs. VIII-X, we will describe some of the observational consequences of DM superfluidity for other systems, including galaxy clusters, dwarf satellite galaxies, globular clusters and ultra-diffuse galaxies residing in clusters. A key distinction compared to the MOND predictions for satellites, in particular, is that the MON-Dian force in our case is mediated by superfluid phonons and therefore requires superfluidity. In particular, the phonon-mediated force between two bodies only applies if both bodies reside within the superfluid region. If one body is inside while the other is outside, the only force acting on them is gravity.

The same applies to the so-called external field effect (EFE), an essential aspect of MOND phenomenology. Consider a subsystem with low internal acceleration a int a 0 in the presence of a large homogeneous external acceleration a ext a 0 . In General Relativity (GR), we know from the equivalence principle that a ext has no physical consequence and can be removed by moving to the freely-falling elevator. In MOND, however, the background acceleration is physical and renders the subsystem Newtonian. The EFE is an example of a more general phenomenon in scalar field theories known as kinetic screening [START_REF] Babichev | k-Mouflage gravity[END_REF][START_REF] Brax | Screening fifth forces in k-essence and DBI models[END_REF][START_REF] Burrage | Screening of scalar fields in Dirac-Born-Infeld theory[END_REF]. In theories with gradient interactions, non-linearities in the scalar field gradient -the scalar acceleration -can result in the suppression of the scalar field effects and the local recovery of standard gravity. See [START_REF] Joyce | Beyond the Cosmological Standard Model[END_REF] for a review.

The key difference in the superfluid context is that the EFE is the result of phonon non-linearities, and therefore only applies within the superfluid core. In particular, close-by globular clusters, such as Pal 5 [START_REF] Thomas | Stellar streams as gravitational experiments II. Asymmetric tails of globular cluster streams[END_REF], and satellite galaxies residing within the superfluid core of the Milky Way should follow MOND predictions with the EFE. On the other hand, globular clusters at large distances from the Milky Way, such as NGC 2419 [START_REF] Ibata | The globular cluster NGC 2419: a crucible for theories of gravity[END_REF] or Pal 14 [START_REF] Jordi | Testing Fundamental Physics with Distant Star Clusters: Analysis of Observational Data on Palomar 14[END_REF], are expected to be Newtonian as long as they do not harbor their own DM halo. The same applies to tidal dwarf galaxies resulting from the interaction of massive spiral galaxies. Those are not expected to harbor a significant DM halo, and thus should be Newtonian [START_REF] Lelli | Gas dynamics in tidal dwarf galaxies: Disc formation at z = 0[END_REF][START_REF] Lelli | Tidal Dwarf Galaxies: Disc Formation at z = 0[END_REF] as long as they are located outside the superfluid core of their host.

As a preview of the main results of this paper, Table I summarizes the observational consequences of superfluid DM for various systems.

II. EFFECTIVE DESCRIPTION OF SUPERFLUID DM

In field theory language, an (abelian) superfluid is described by the theory of a spontaneously broken global U (1) symmetry, in a state of finite U (1) charge density. At low energy the relevant degree of freedom is the Goldstone boson for the broken symmetry -the phonon field φ. The U (1) symmetry acts non-linearly on φ as a shift symmetry, φ → φ + c. Furthermore, in the non-relativistic regime the theory should be Galilean invariant (ignoring gravity). At finite chemical potential µ, the most general effective theory at leading order in derivatives consistent with these symmetries is [START_REF] Greiter | Hydrodynamic Relations in Superconductivity[END_REF][START_REF] Son | Low-energy quantum effective action for relativistic superfluids[END_REF]

L T =0 = P (X) , (2)
where X = µ-mΦ+ φ-( ∇φ) 2 /2m. Here m is the particle mass and Φ the Newtonian gravitational potential. The conjecture of [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] is that the DM superfluid phonons are governed by the MOND Lagrangian (see [START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF][START_REF] Bruneton | Fieldtheoretical formulations of MOND-like gravity[END_REF])

L DM, T =0 = 2Λ(2m) 3/2 3 X |X| . (3)
Remarkably, (3) is strikingly similar to that of the Unitary Fermi Gas (UFG) [START_REF] Zwerger | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF][START_REF] Giorgini | Theory of ultracold atomic Fermi gases[END_REF], which has generated much excitement in the cold atom community in recent years. Indeed, the fractional power of X would be strange if (3) described a fundamental scalar field. As a theory of phonons, however, the power determines the superfluid equation of state, and fractional powers are not uncommon. Indeed, the effective field theory for the UFG superfluid is L ∼ X n , where n = 5/2 in 3+1 dimensions and 3/2 in 2+1 dimensions, and is therefore also nonanalytic [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF].

To mediate a force between baryons, DM phonons must couple to the baryon density as

L int = αΛ φ M Pl ρ b , (4)
with α being a constant, ρ b the baryonic density, and M Pl the Planck mass. At zero temperature, the effective theory thus has three parameters: the particle mass m, a parameter Λ related to the self-interaction strength, and the coupling constant α between phonons and baryons.

A fourth parameter of the particles themselves is their self-interaction cross-section σ setting the conditions for their thermalization, while a fifth parameter β will later be introduced to accommodate for finite-temperature effects. 1From the superfluid perspective, the form of the coupling term ( 4) is unusual: it breaks the U (1) shift symmetry explicitly (albeit softly), making the phonon φ a pseudo-Goldstone boson. The origin of this particle non-conserving term could be non-perturbative; or more simply, there could be a soft fundamental coupling between DM and baryons which explicitly breaks the relevant symmetries. One way or the other, this makes our quantum liquid a pseudo-superfluid. In particular, the phonon will acquire a mass via radiative corrections, though it is easy to check that the explicit breaking is soft enough that this has no observable effects on galactic scales.

As mentioned already, we do not assume any nonminimal coupling to photons and/or gravitons. Both electromagnetic and gravitational waves travel at the speed of light, consistent with the recent observations of GW170817 [START_REF] Abbott | GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral[END_REF].

Let us briefly review why the action, composed of the terms ( 3) and ( 4), gives rise to the MOND force law. Assuming a static profile, the phonon equation of motion then is

∇ •   ( ∇φ) 2 -2mμ ( ∇φ) 2 -2mμ ∇φ   = αρ b 2M Pl , (5)
where μ ≡ µ -mΦ. In the limit ( ∇φ)2 2mμ, the solution is, ignoring a homogeneous curl term,

| ∇φ| ∇φ αM Pl a b , (6)
where a b is the Newtonian acceleration due to baryons only. The φ-mediated acceleration that derives from ( 4) is

a φ = α Λ M Pl ∇φ . (7)
Thus ( 6) implies

a φ = α 3 Λ 2 M Pl a b . (8)
As advocated, this matches the deep-MOND form with critical acceleration

a 0 = α 3 Λ 2 M Pl . (9)
In the context of MOND, it has long been established that the best-fit value for a 0 from galactic rotation curves is

a MOND 0 1.2 × 10 -8 cm/s 2 . ( 10
)
In our case, the situation is different for two reasons: i) the deep-MOND form ( 8) is not exact and only applies in the regime ( ∇φ) 2 2mμ; ii) the total acceleration experienced by baryons includes not only a b and a φ , but also a DM -the Newtonian acceleration from the DM halo itself. Because the superfluid core is pressure supported, the DM density profile is approximately homogeneous in the central region of galaxies. Its density is moreover rather low. This should be contrasted with the CDM scenario, where DM-only simulations display a cuspy profile [START_REF] Navarro | A Universal density profile from hierarchical clustering[END_REF].

The story is further complicated by the fact that perturbations around this zero-temperature, static background are unstable (ghost-like). However it was argued in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] that this instability can naturally be cured by finite-temperature effects. Indeed, our DM particles in galactic halos have non-zero temperature, owing to their velocity dispersion, hence we expect the zerotemperature Lagrangian (3) to receive finite-temperature corrections in galaxies. One way to understand the instability physically is that, if a galactic halo had zero temperature, then the emergent MOND force would destabilize it.

As it was originally suggested by Landau, a superfluid at finite temperature phenomenologically behaves as a mixture of two fluids. Only a fraction of the total mass forms a degenerate quantum liquid, hydrodynamically described as a potential flow carrying no entropy. The rest of the mass is stored in thermally excited states, behaving as a regular fluid which carries the entropy of the entire system. The general, finite-temperature effective theory, once again at leading order in derivatives, is [START_REF] Nicolis | Low-energy effective field theory for finitetemperature relativistic superfluids[END_REF]

L T =0 = F (X, B, Y ) . ( 11
)
This generalizes (2) through the dependence on two additional scalar quantities B and Y defined by

B ≡ det ∂ µ ψ I ∂ µ ψ J ; Y ≡ μ + φ + v • ∇φ , (12)
where ψ I ( x, t), I = 1, 2, 3 and v are respectively the Lagrangian coordinates and velocity vector of the normal fluid component. Physically, B measures the density of the normal fluid, while Y is the scalar product of the fluid velocities.

The general finite-temperature effective theory ( 11) is only mildly constrained, and further specification of its form requires first-principle knowledge of the superfluid constituents. Since a fundamental underlying description of our DM superfluid is still lacking, [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] proceeded empirically and considered various finite-temperature actions leading to MOND-like static profiles with stable perturbations. One such form, which in this paper we will take as our working model, is

L = 2Λ(2m) 3/2 3 X |X -βY | -α Λ M Pl φρ b , (13)
where for completeness we have included the coupling term ( 4). This fiducial action may appear ad hoc, and simpler finite-temperatures corrections were considered in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF], but it has the practical advantage of making the algebra straightforward. Here β is a dimensionless constant that parametrizes finite-temperature effects. In particular the limit β → 0 recovers the original action [START_REF] Boylan-Kolchin | The Milky Way's bright satellites as an apparent failure of LCDM[END_REF].

In what follows we will fix β = 2 as the fiducial value.

The temperature dependence of the theory is also welcome for another reason. To obtain an acceptable background cosmology and linear perturbation growth, both Λ and α must assume different values cosmologically than in galaxies. This is not unreasonable, since parameters of the superfluid effective field theory are expected to depend on temperature, and thus on velocity, most naturally through the ratio of the temperature over the critical temperature for phase-transition T /T c . In particular they can assume different values on cosmological scales (where T /T c ∼ 10 -28 ) than in galaxies (where T /T c ∼ 10 -6 -10 -2 ). In [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] it was shown that α must be ∼ 10 -4 smaller cosmologically, while Λ must be ∼ 10 4 larger, in order to obtain an acceptable cosmology. Hence, the phonon-baryon coupling (4) can be taken as constant.

In the present work, we perform a more detailed analysis of the rotation curves for a more realistic distribution of baryons in two different types of galaxies. As a result we find that, although subdominant in the central parts of the galaxy, the gravitational force due to the superfluid is not always negligible and can be as large as ∼ 30%. To compensate for this effect and reproduce the observed rotation curves, we will use a somewhat lower best-fit value, a 0 0.87 × 10 -8 cm/s 2 .

Another important point to bear in mind concerns the DM mass distribution throughout the halo. It is well known that galaxy clusters pose a problem for MOND. Therefore in order to avoid the tension between our model and observations on cluster scales, it would be most desirable to avoid having a superfluid halo in clusters altogether. It was argued in [START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] that this is possible if the mass of the DM particles is appropriately chosen, and then partly tested in [START_REF] Hodson | Galaxy Clusters in the Context of Superfluid Dark Matter[END_REF], albeit with slightly more simplified assumptions than in the following. Our more careful analysis will yield somewhat altered bounds on the parameters m, Λ and α of the theory.

III. BACK-OF-THE-ENVELOPE ESTIMATES

In this Section we revisit the derivation of the DM-only superfluid halo profile originally given in [START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF], and refine various steps in the argument. As in standard ΛCDM, we expect virialization of halos to proceed through violent relaxation, a manifestly out-of-equilibrium process that should take DM particles out of any pre-existing superfluid state. Therefore, at virialization the halo is expected to have an NFW profile [START_REF] Navarro | A Universal density profile from hierarchical clustering[END_REF]. This initial profile will then be altered in the central regions of the halo, where DM thermalization and condensation occur.

The final equilibrium profile consists of a superfluid core, within which DM particles interact sufficiently frequently to thermalize and form a Bose Einstein condensate, surrounded by an envelope describing approximately collisionless particles, which we assume for concreteness to be of the NFW form. For the purpose of the simple analytical estimates of this Section, we will approximate the collisionless envelope as a ρ ∼ r -3 tail of the NFW profile. In the next Section, we will come back and describe a more realistic matching procedure to the full NFW profile, to be applied using numerical analysis.

The conditions for thermal equilibrium and Bose Einstein condensation both depend on the density and velocity profiles, and therefore on the distance r from the center. The radii within which thermalization and condensation take place will be denoted respectively by R T and R degen . It is easy to argue, however, that R T must be smaller than R degen . Indeed, otherwise there would be a region R degen < r < R T where thermalization is reached without degeneracy, which would in turn require an unacceptably large self-interaction cross section.

On the other hand, since thermalization is necessary to achieve condensation, it would be inconsistent to assume that DM remains in a coherent BEC state for R T < r < R degen . Therefore, we conclude there is only one radius of interest, R T . For r < R T , our DM is degenerate and thermalizes through self-interactions. Within this region our density profile consists of a superfluid core, with approximately homogeneous density. As mentioned above, for the purpose of making analytical estimates we assume that at r = R T the profile switches to the NFW fall-off ρ ∼ r -3 , which extends up to the virial radius.

Before deriving an expression for R T , we should be more precise about how large the self-interaction cross section σ (assumed velocity-independent) should be to achieve thermal equilibrium. Thermal equilibrium requires the interaction rate of DM particles to be greater than the inverse of the dynamical time, which we take here as

t dyn ≡ r v(r) , ( 14
)
where v is the characteristic velocity dispersion at r. It is natural to compare DM interactions to t dyn , since the latter sets the characteristic time scale for the motion of perturbers (such as DM subhalos or halo stars) which can disturb the equilibrium state.

The local self-interaction rate in a region with density ρ and velocity v is given by

Γ = σ m N vρ , (15)
where N ≥ 1 denotes the Bose degeneracy factor,

N = ρ m 2π mv 3 . (16)
We will then require Γ > t -1 dyn for thermal equilibrium. In the above, both v and ρ depend on the distance r from the center of the halo. This is a key difference from [START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF], where v and ρ were taken as the virial velocity and density respectively. In what follows we will assume N 1, consistent with the degeneracy discussion above. It is easy to show a posteriori that this is justified for galaxies if m ∼ eV.

We will assume for concreteness that thermalization takes place in the outskirts of the halo where the density distribution is given by the tail of the NFW profile

ρ(r) = ρ(R 200 ) R 200 r 3 . (17)
In this region, up to slowly varying factors (related to the logarithmic growth of the enclosed mass and the varying ratio of velocity dispersion over circular velocity), the velocity profile roughly falls off in a Keplerian fashion as

v ∼ v 200 R 200 r . (18)
In the above, R 200 and v 200 are as usual defined as the radius and velocity at which the mean density is 200 times the present critical density:

ρ 200 = 200 3H 2 0 8πG N 1.95 × 10 -27 g/cm 3 ; R 200 = 3M 4πρ 200 1/3 203 M 10 12 M 1/3 kpc ; v 200 = 1 3 G N M R 200 85 M 10 12 M 1/3 km/s , ( 19
)
where we have assumed H 0 = 70 kms -1 Mpc -1 for concreteness.

Meanwhile, the quantity ρ(R 200 ), on the other hand, is the local density at R 200 , which is of course less than the mean density ρ 200 within that radius. For instance, for a NFW profile with concentration in the range 5 < c < 15, one finds 4 < ∼ ρ200 ρ(R200) < ∼ 6. For the purposes of our back-of-the-envelope estimate, we will assume the central value

ρ 200 ρ(R 200 ) = 5 . ( 20
)
Substituting these expressions, the requirement Γ > t -1 dyn translates to a bound on the cross section:

σ m > ∼ m eV 4 M 10 12 M 2/3 r R 200 7/2 0.2 cm 2 g . (21)
It is thus easier to reach equilibrium in the inner regions where the density is high.

Conversely, for a given σ/m, (21) translates to a threshold radius R T within which thermal equilibrium is achieved. The result is

R T 310 kpc m eV -8/7 M 10 12 M 1/7 σ/m cm 2 /g 2/7 . ( 22
)
Note that this bound is a rather insensitive function of M , and depends more sensitively on m.

Since we hope to explain galaxy rotation curves using DM superfluidity, the superfluid core should extend at least as far as the last observed point on the rotation curve in each galaxy. In other words, the viability of our scenario requires R T to be larger than the radius within which the circular motion of stars and gas is observed. For instance, consider a Milky Way-like galaxy with M 10 12 M , whose rotation curve (or the equivalent circular velocity curve derived from stellar kinematics, see [START_REF] Xue | The Milky Way's Circular Velocity Curve to 60 kpc and an Estimate of the Dark Matter Halo Mass from the Kinematics of 2400 SDSS Blue Horizontal-Branch Stars[END_REF]) is typically measured out to 60 kpc. For such a galaxy, we impose that R T > 60 kpc ;

for M = 10 12 M .

Using [START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF] this translates to an upper bound on the mass of the DM particle:

m 4.2 σ/m cm 2 /g 1/4
eV .

Smaller and less massive galaxies result in a somewhat weaker bound.

The bound (24) on the DM particle mass is the main result of this Section. It shows that for values of σ/m satisfying the merging-cluster bound ∼ 1 cm 2 /g [85-88], m must be somewhat below 4 eV. The dependence on the cross section is rather weak, however, scaling as the 1/4 power. It should be mentioned that the upper bound [START_REF] Desmond | The scatter, residual correlations and curvature of the SPARC baryonic Tully-Fisher relation[END_REF] would be somewhat tighter had we assumed a ρ ∝ r -2 transition density profile outside the superfluid core, instead of ρ ∝ r -3 .

In practice, for the fiducial parameter values that we will use later to fit rotation curves, we will find that the superfluid cores make up a modest fraction of the total DM mass, ranging between ∼ 20% and ∼ 50%. This has important consequences for observations. Firstly, since most of the mass is in the approximately collisionless NFW envelope, our halos should be triaxial near the virial radius, exactly as in ΛCDM halos and consistent with observations [START_REF] Clampitt | Lensing measurements of the ellipticity of luminous red galaxies dark matter haloes[END_REF].

Secondly, most of the gravitational lensing signal will come from the NFW envelope, particularly for HSB galaxies, hence we need not assume any non-minimal coupling between DM superfluid phonons and photons. If we assume that the relation between stellar mass and halo mass is mostly unaltered with respect to ΛCDM (we will explicitly show that this leads to acceptable rotation curves in Sec. VII), then galaxy-galaxy lensing statistics would be mostly unaletered too. In particular, if Milgrom's law were valid as a fundamental law out to large distances from the baryonic bulk of galaxies, then one would expect a fine relation between the baryonic mass and the weak lensing signal, of the form r E ∝ √ M b , where r E is the Einstein radius [START_REF] Tian | The relation between stellar mass and weak lensing signal around galaxies: implications for modified Newtonian dynamics[END_REF]. Such a relation is not expected to hold in the superfluid DM context.

IV. LOWER BOUND ON HALO MASS

Simulations indicate that virialization takes place when the average density of a perturbation reaches about two hundred times the mean density. Simply put, the over-density is then sufficient for the perturbation to withstand the Hubble expansion. This fact, combined with the observation that pressure supported halos tend to have lower central densities, raises the caution flag regarding the lightest possible haloes one could have in the model in question (see [START_REF] Hui | Ultralight scalars as cosmological dark matter[END_REF] and references therein).

As argued in Sec. III, the galactic halo consists of a superfluid core surrounded by the NFW envelope. However, the picture simplifies for low mass galaxies, for which ρ T ρ 200 . Here, ρ T denotes the DM density at the thermalization radius given by [START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF]. In this case most of the mass resides in the superfluid phase. Ignoring the contribution from phonon gradients, the equation of state implied by ( 3) is

P = ρ 3 12Λ 2 m 6 . (25)
In [START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] it was shown that this equation of state, together with the condition of hydrostatic equilibrium, results in a cored density profile with central density

ρ 0 M 10 12 M 2/5
Λ meV

m 3 eV 3 6/5 6.2 × 10 -25 g/cm 3 . ( 26
)
For consistency, the central density ρ 0 should be larger than ρ 200 , given in [START_REF] Milgrom | A Modification of the Newtonian dynamics: Implications for galaxies[END_REF]. (Instead of ρ 0 , [START_REF] Hui | Ultralight scalars as cosmological dark matter[END_REF] uses the average density within r 1/2 , resulting in a tighter constraint.) This translates to a minimum halo mass M min for the superfluid scenario

M min 12 Λ meV m 3 eV 3 -3 10 9 M . (27)
The superfluid pressure suppresses the formation of halos of mass less than M min .

To be consistent with observations, M min should be sufficiently small to accommodate the least massive galaxies observed. Ultra-faint dwarf satellites around the Milky Way are estimated to have virial masses between 10 8 and 10 9 M (before tidal stripping). Demanding that M min < ∼ 10 9 M to be generous translates to Λm 3 0.08 meV × eV 3 .

(

) 28
The combination Λm 3 , which will be a recurring theme in our analysis, traces back to the equation of state [START_REF] Oman | The unexpected diversity of dwarf galaxy rotation curves[END_REF] depending on this combination of parameters. In particular, we will see in the next Section that rotation curves are controlled precisely by this combination and tend to favor lower values of Λm 3 . There is therefore an interesting tension between the ability to form light enough halos, and obtaining a reasonable fit to rotation curves.

In practice we will find that values of Λm 3 that barely satisfy [START_REF] Bekenstein | Relativistic gravitation theory for the MOND paradigm[END_REF] result in acceptable fits.

V. HALO PROFILE: ALGORITHM

In this Section we show how to derive the superfluid DM density profile in a galaxy, and the predicted rotation curve, for a given (observed) baryonic density profile. As already discussed, the central parts of the halo are expected to be in the superfluid phase while the outskirts of the halo are in the normal phase. Apart from theory parameters, our density profile will depend ultimately on a single parameter -the total halo mass. This parameter will in turn fix the central value of the potential and DM density. To go further, we thus need a prescription relating the baryonic mass to the halo mass.

Since cosmological simulations of DM superfluid are yet to be devised, we must speculate somewhat on a realistic scenario, which will have to be backed by simulations. To cover the range of possibilities, we will focus on two radically different scenarios. The first scenario assumes that the DM halo mass function is the same as in ΛCDM (at least at halo masses above the cutoff), and that the same abundance matching prescriptions can be applied (e.g., [START_REF] Behroozi | The Average Star Formation Histories of Galaxies in Dark Matter Halos from z =0-8[END_REF]). This scenario implies very different baryon fractions in galaxies of different stellar masses. The second scenario assumes that a given DM halo can, to a first approximation, be thought of as a microcosm of the whole Universe, and can a priori be expected to host a fraction of baryons representative of the cosmological ratio. On galaxy scales, we could expect a constant fraction (typically of the order of one half) of those baryons to condense towards the center of the halo to form the visible galaxy, thus assuming a negligible effect of feedback. This second scenario thus implies assuming a constant DM over visible baryons ratio M DM /M b = 10. Remarkably, we will show in Sec. VII that those two extreme assumptions have a relatively small influence on the rotation curve predictions.

Let us begin with the superfluid region. The Poisson equation is as usual given by

∇ 2 Φ = 4πG N (ρ SF + ρ b ) . (29)
The baryonic mass density ρ b is prescribed by observations. The superfluid DM mass density ρ SF can readily obtained by differentiating the action [START_REF] Mcgaugh | The Baryonic Tully-Fisher Relation of Gas Rich Galaxies as a Test of LCDM and MOND[END_REF] with respect to the Newtonian gravitational potential Φ:

ρ SF = 2 √ 2m 5/2 Λ 3(β -1)μ + (3 -β) ( ∇φ)2 2m 3 (β -1)μ + ( ∇φ)2 2m . (30)
We will solve Poisson's equation together with the phonon equation

∇ •   ( ∇φ) 2 + 2m 2β 3 -1 μ ( ∇φ) 2 + 2m(β -1)μ ∇φ   = αρ b 2M Pl . (31)
The baryon density is of course not spherically symmetric in general. Nevertheless, if we ignore a homogeneous "pure curl" solution, then this equation can be readily integrated:

( ∇φ) 2 + 2m 2β 3 -1 μ ( ∇φ) 2 + 2m(β -1)μ ∇φ = αM Pl a b , (32)
where a b is the Newtonian acceleration due to baryons only. The approximation of neglecting a curl term has been shown in MOND to be accurate to within ∼ 20% [START_REF] Brada | Exact solutions and approximations of MOND fields of disk galaxies[END_REF]. The upshot of this approximation is that (32) implies an algebraic equation for the gradient of the phonon field. Specifically, it can be massaged into a cubic equation for ( ∇φ) 2 . Only one root of this cubic equation is physical -it is easy to see that the other two roots are complex at large distances where ( ∇φ) 2 → 0. The explicit expression for the physical solution is complicated and not particularly illuminating. Suffice to say the result is we have now determined the phonon gradient as a function of Φ and a b :

∇φ = F (Φ, a b ) a b . (33)
Substituting the resulting ( ∇φ) 2 into (30) then gives the superfluid mass density in terms of Φ and a b :

ρ SF = ρ SF (Φ, a b ) . (34)
Let us stress again that this result holds for an arbitrary (non spherically-symmetric) baryon distribution. Our only approximation was neglecting a pure-curl homogeneous solution in [START_REF] Sanders | The Virial Discrepancy in Clusters of Galaxies in the Context of Modified Newtonian Dynamics[END_REF]. Given a general ρ b , the DM-free Poisson equation ∇ • a b = -4πG N ρ b can be solved for a b . In turn, the resulting a b can be substituted in [START_REF] Clowe | A direct empirical proof of the existence of dark matter[END_REF] to obtain ρ SF as a function of Φ, which remains as the last unknown variable.

The gravitational potential is of course determined by Poisson's equation [START_REF] Milgrom | Bimetric MOND gravity[END_REF]. Upon substituting ρ SF (Φ, a b ), this becomes a highly non-linear equation, whose general solution absent any symmetries is complicated to find. To simplify the analysis in the present paper, we will replace the real baryon distribution by a spherical density profile ρb (r) which gives the same purely-baryonic Newtonian acceleration a b in the disk. The simplified Poisson equation becomes

1 r 2 d dr r 2 dΦ dr = 4πG N ρb (r) + ρ SF (Φ(r), a b ) , (35)
where ãb = G N Mb (r)/r 2 is the purely-baryonic Newtonian acceleration derived from the spherical distribution ρb (r).

We integrate this equation numerically with two boundary conditions: i) imposing regularity at the origin, dΦ/dr| r=0 = 0; ii) specifying the value of Φ at the origin, Φ(r = 0). The latter determines the size and the mass of the superfluid core. In other words, halo cores of different sizes have different central values of the gravitational potential.

Upon matching to the normal-phase NFW envelope, discussed below, the mass of the superfluid core in turn fixes the total mass of the halo. In practice, when fitting rotation curves (Sec. VII) we will proceed through reverse engineering: we first prescribe the total mass M DM of the halo, either as a fixed factor of the baryonic mass (e.g., M DM = 10 M b ) or to its value determined by the abundance matching prescription in ΛCDM; and then find numerically the value of Φ(r = 0) which gives the desired M DM .

A. Matching to NFW

Let us discuss in more detail our matching procedure to the NFW envelope. For this purpose we will improve on the basic analytical estimates of Sec. III, where we assumed for simplicity a sharp transition at the thermalization radius R T to a 1/r 3 NFW tail. The first improvement is that we now match to a full NFW profile

ρ NFW (r) = ρ c x(1 + x) 2 ; x ≡ r r s . ( 36
)
The density ρ c and concentration radius r s (equivalently, the concentration parameter c = R 200 /r s ) are determined by our matching conditions. The second improvement pertains to the matching radius. In reality, we do not expect a sharp transition at R T . Instead, there should be a transition region around R T within which the DM is nearly in equilibrium but is incapable of maintaining long range coherence. In the simpler context of self-interacting DM, for instance, recent simulations have shown that the density profile starts to deviate from the NFW profile at a radius within which DM particles had a chance to re-scatter at least once over the lifetime of the halo [START_REF] Kaplinghat | Dark Matter Halos as Particle Colliders: Unified Solution to Small-Scale Structure Puzzles from Dwarfs to Clusters[END_REF].

To shed light on the nature and spatial extent of the transition regime, we introduce a second characteristic radius R NFW as the location where the profile can be consistently matched to an NFW profile. The matching conditions at R NFW are:

1. Continuity of density: ρ SF = ρ NFW .

2. Continuity of pressure: P SF = P NFW .

To evaluate the second condition, we compute the radial superfluid pressure using

P SF = L + ∂L ∂X (∂rφ) 2
m , reducing to the following expression for the Lagrangian density (13) (with β = 2)

P SF = 2(2m) 3/2 Λ 3 μ μ + ( ∇φ) 2 2m + 2 ( ∇φ) 2 2m 2 ( ∇φ) 2 2m + μ . (37)
The NFW pressure simply follows from integrating the hydrostatic equilibrium condition:

P NFW (R NFW ) = ∞ RNFW dr ρ NFW (r) G N M (r) r 2 , (38)
where M (r) is the total (DM + baryonic) mass enclosed.

Together with the constraint of fixed total mass and these matching conditions we need to impose an additional requirement in order to uniquely fix the NFW profile and determine R NFW . In particular, we further require that the concentration parameter thus determined is close to its ΛCDM value.

Meanwhile, as before the thermalization radius R T is defined as the location where the interaction rate Γ given by ( 15) drops to t -1 dyn = √ G N ρ. As another refinement on our earlier estimates, we will now evaluate this selfconsistently by keeping tracking of Γ and t dyn as we integrate the equations of motion numerically from the center. (Note that for this purpose Γ correctly includes the degeneracy factor N appropriate for the superfluid region, but is otherwise a "vacuum" scattering rate, i.e., it ignores possible dressing effects from the superfluid medium.) For the range of theory parameters that give a good fit to rotation curves, we found that R NFW and R T agree to within 30%. A calculation of the exact profile between these radii is beyond the scope of the present work and is left for the future. But these estimates at least instruct us on the likely spatial extent of the transition region.

B. Fiducial theory parameters

The DM density profile is determined by two combinations of theory input parameters. The first combination

is m σ/m cm 2 /g -1/4
, which enters in Γt dyn and hence determines the size of the superfluid region. As fiducial values we will set

m = 1 eV ; σ m = 0.01 cm 2 g . ( 39
)
These values are optimal for ensuring that the superfluid core is large enough in galaxies to encompass the baryonic disk, while being small enough in clusters not to contradict observations. The second combination is Λm 3 . This controls the superfluid pressure and correspondingly the density profile of the core. Furthermore, as argued in Sec. IV, this quantity also determines the lowest possible halo mass that can form in our scenario. As a fiducial value we will choose Λm 3 = 0.05 meV × eV 3 .

(

) 40
Although strictly speaking this violates the bound [START_REF] Bekenstein | Relativistic gravitation theory for the MOND paradigm[END_REF] on the lowest mass halo, we will see in Sec. VII that values of Λm 3 consistent with ( 28) also give reasonable fits to rotation curves. Finally, the coupling constant will be set to α = 5.7 to allow for our best-fit value of a 0 = 0.87 × 10 -8 cm/s 2 , while we have set from the start β = 2 in the action [START_REF] Mcgaugh | The Baryonic Tully-Fisher Relation of Gas Rich Galaxies as a Test of LCDM and MOND[END_REF].

VI. HALO PROFILE FOR A TOY BARYON DISTRIBUTION

Before comparing to actual data, it is instructive to calculate our predicted profile and rotation curve for a toy baryonic distribution. For this purpose we assume a "spherical exponential" profile

ρ toy b (r) = M b 8πL 3 e -r/L , (41)
where the characteristic scale L plays the role of a radial "scale length". Following the algorithm outlined in the previous Section, we can calculate the DM density profile for this exponential baryon distribution once we fix the DM halo mass. Figure 1

a difference of 13%. Relative to R 200 97 kpc for this galaxy, (42) implies a significant superfluid core, encompassing 39% of the total DM mass.

The matching conditions then fix the concentration parameter of the NFW envelope c = 6.5. For comparison, also shown in the Figure are the "pure" superfluid profile (T = 0, no baryons) derived in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] (blue dasheddotted curve), together with the NFW profile (red dashed curve) for the specified mass and concentration. Note that the small feature within r < ∼ 2 kpc is due to the energy density in phonon gradients dominating the superfluid DM mass density.

Using this profile, we can also calculate the radial acceleration on a test baryonic particle. Within the superfluid core this is given by the sum of three contributions:

a = a b (r) + a SF (r) + a phonon (r) . ( 43
)
The first term is just the baryonic Newtonian acceleration. The second term is the gravitational acceleration from the superfluid core, a SF = G N M SF (r)/r 2 . The third term is the phonon-mediated acceleration [START_REF] Pawlowski | The rotationally stabilized VPOS and predicted proper motions of the Milky Way satellite galaxies[END_REF]:

a phonon (r) = α Λ M Pl dφ dr . ( 44
)
The strength of this term is set by α and Λ, or equivalently the critical acceleration a 0 given by ( 9). It is instructive to compare our predicted rotation curve with the MOND prediction for a fixed toy baryonic profile. For the MOND calculation, we assume the 'simple' interpolating function:

a MOND (r) = a b (r) 2   1 + 1 + 4a MOND 0 a b (r)   , ( 45
)
where a b (r) = G N M b (r)/r 2 . This is well known to provide good fits to galaxy rotation curves [START_REF] Famaey | Modified Newtonian dynamics in the Milky Way[END_REF][START_REF] Gentile | THINGS about MOND[END_REF], as well as elliptical galaxies [START_REF] Richtler | Remarks on the properties of elliptical galaxies in modified Newtonian dynamics[END_REF][START_REF] Dabringhausen | Understanding the internal dynamics of elliptical galaxies without non-baryonic dark matter[END_REF][START_REF] Chae | Radial Acceleration Relation by Dark Matter and Baryons in the Intermediate Acceleration Regime[END_REF]. Comparison of our predicted radial acceleration curve (solid black) for the toy baryon exponential profile [START_REF] Ibata | The globular cluster NGC 2419: a crucible for theories of gravity[END_REF]. The parameters are the same as in Fig. 1. The critical acceleration is a0 = 0.87 × 10 -8 cm/s 2 . The MOND prediction (red dashed curve) with the simple interpolating function ( 45) and a MOND 0 1.2 × 10 -8 cm/s 2 for the same baryonic profile is found to agree well (within 25%) with our prediction.

same M b = 10 10 M and L = 2 kpc galaxy as before. All parameters are identical to Fig. 1. The critical acceleration a 0 , which normalizes the phonon-mediated force is assumed to be

a 0 = 0.87 × 10 -8 cm/s 2 , (46)
which, as we will show in Sec. VII, offers an excellent fit to actual galaxy rotation curves. With Λ = 0.05 meV, this corresponds to α 5.7.

As can be seen from the Figure, we find good agreement with the MOND prediction, to within < ∼ 25%. When fitting galaxy rotation curves in Sec. VII, we will keep our parameters fixed to the values quoted above.

VII. GALAXY ROTATION CURVES

In this Section we apply the method outlined above to fit actual galaxy rotation curves. As it would represent an extremely long exercise to fit a large sample of galaxy rotation curves, we delay such a systematic analysis to future work, and only analyze here, for illustrative purposes, the rotation curves of two disk galaxies, characteristic of the LSB and HSB regime. The LSB galaxy is IC 2574, with total baryonic content M b ∼ 2 × 10 9 M . The HSB galaxy is UGC 2953 with M b ∼ 1.6 × 10 11 M .

Our numerical calculation of the superfluid profile and phonon-mediated force in Sec. VI was performed for a spherically symmetric baryon distribution, whereas the real baryonic distribution is of course far from spherical. To proceed we adopt a hybrid method. As before, the acceleration consists of three terms:

a hybrid = a actual b + a DM + a phonon . (47)
The different terms are computed as follows: 3. The phonon-mediated acceleration a phonon is then computed using [START_REF] Sanders | The Virial Discrepancy in Clusters of Galaxies in the Context of Modified Newtonian Dynamics[END_REF] sourced by the actual baryon distribution a actual b , but with the Newtonian potential Φ(r) taken from our numerical solution to the spherically symmetric problem. Interestingly, despite the non-standard form for the phonon equation [START_REF] Pawlowski | The VPOS: a vast polar structure of satellite galaxies, globular clusters and streams around the Milky Way[END_REF], for instance compared to its Bekenstein-Milgrom counterpart [START_REF] Bekenstein | Does the missing mass problem signal the breakdown of Newtonian gravity?[END_REF], we nevertheless found that a phonon √ a 0 a b to within a couple of percent.

In other words, the phonon force closely matches the deep-MOND acceleration.

A. LSB galaxy (IC 2574)

Our example of an LSB galaxy is IC 2574 (M b ∼ 2×10 9 M ). The simplified spherical baryon distribution ρ spherical b , necessary to evaluate the DM density profile, was modeled as a constant-density core out to 7.85 kpc, matched to a 1/r 2 profile out to 9.22 kpc, such that the total mass is M b 2 × 10 9 M , and that the slope of the baryonic rotation curve is roughly reproduced.

Figure 3 shows the DM profile for the fiducial theory parameters [START_REF] Hees | Combined Solar System and rotation curve constraints on MOND[END_REF] and [START_REF] Babichev | Improving relativistic MOND with Galileon k-mouflage[END_REF], assuming a constant baryon fraction M DM = 10 M b . The resulting transition radii (dashed vertical lines) in this case are FIG. 4: Predicted rotation curve of IC 2574, evaluated using the hybrid method outlined in the text. The orange points are the data from [START_REF] Lelli | SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves[END_REF] assuming a distance of 3 Mpc [START_REF] Tully | Our Peculiar Motion Away from the Local Void[END_REF]. The black curve is for MDM = 10 M b , and assumes the same theory parameter values as Fig. 3 and a0 = 0.87×10 -8 cm/s 2 . The red dashed curve assumes the same parameter values except for MDM 50 M b , which is the upper total mass obtained from ΛCDM abundance matching. The gray band corresponds to a0 ∈ (0.6, 1.2) × 10 -8 cm/s 2 , while the light blue band corresponds to Λ ∈ (0.02, 0.1) meV, with MDM = 10 M b and other parameters fixed to their fiducial values. The inset zooms in to show the width of the light blue band.

R NFW = 40 kpc ; R T = 49 kpc , (48)
a difference of 18%. The matching conditions then fix the concentration radius to r s = 10 kpc. Relative to R 200 57 kpc for this galaxy, (48) implies a relatively large superfluid core, encompassing 55% of the total DM mass.

In the abundance matching case [START_REF] Behroozi | The Average Star Formation Histories of Galaxies in Dark Matter Halos from z =0-8[END_REF], where M DM 51 M b and R 200 96 kpc, we instead get R NFW = 47 kpc and R T = 67 kpc; corresponding to 30% difference. In this case the superfluid core containing only 33% of the total DM mass.

Figure 4 compares the predicted rotation curve, calculated following the hybrid method described above, with the actual data for IC 2574 [START_REF] Lelli | SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves[END_REF].

• The black solid curve corresponds to M DM = 10 M b , with theory parameters set to their fiducial values.

• The red dashed curve (which in this case is barely distinguishable from the black curve) instead assumes M DM 51 M b , which is consistent with ΛCDM abundance matching [START_REF] Behroozi | The Average Star Formation Histories of Galaxies in Dark Matter Halos from z =0-8[END_REF].

• The shaded bands illustrate the sensitivity of our result to variations in Λ (light blue band, spanning 0.02 to 0.1 meV) and a 0 (gray band, spanning 0.6 to 1.2 × 10 -8 cm/s 2 ), assuming M DM = 10 M b and keeping other parameters fixed.

The superfluid model recovers the observed rotation curve (orange points) rather well over the specified range of parameters. This profile approximates a deprojected Sérsic profile [START_REF] Mellier | Deprojection of the de Vaucouleurs R exp 1/4 brightness profile[END_REF]. The corresponding mass of the bulge is M bulge = 3.5 × 10 10 M .

A sphericized exponential disk Σ disk

r e -r/L disk , with scale length L disk = 4 kpc and central surface density Σ disk = 6.25 × 10 8 M /kpc 2 . This profile is cut off at r = 60 kpc, corresponding to the last measured data point on the rotation curve. The total enclosed mass in the disk is M disk = 1.25×10 11 M .

Figure 5 shows the superfluid DM profile for the fiducial theory parameters ( 39) and ( 40), assuming a constant baryon fraction M DM = 10 M b . The resulting transition radii (dashed vertical lines) are

R NFW = 76 kpc ; R T = 82 kpc , (49)
a difference of only 7%. The matching conditions then fix the concentration radius to r s = 45 kpc. Relative to R 200 245 kpc for this galaxy, (49) implies a modest superfluid core, encompassing 22% of the total DM mass. Figure 6 compares the predicted rotation curve with data (orange points) obtained from [START_REF] Noordermeer | The mass distribution in early-type disk galaxies: Declining rotation curves and correlations with optical properties[END_REF]. The parameter values and conventions are the same as for IC 2574. The main difference is the red dashed curve, where the total DM mass in this case is set to the ΛCDM abundance matching value of [START_REF] Behroozi | The Average Star Formation Histories of Galaxies in Dark Matter Halos from z =0-8[END_REF] M DM = 65 M b . For this case, we get R NFW = 95 kpc and R T = 129 kpc, corresponding to a 26% difference. Therefore, the superfluid core, encompassing about 13% of the total DM mass, is significantly smaller than R 200 446 kpc. Interestingly, we see a slight increase in the rotation curve at large distances. In contrast with MOND, where the rotation curve becomes flat (or slightly decreases, because of the external field effect [START_REF] Hees | Combined Solar System and rotation curve constraints on MOND[END_REF][START_REF] Haghi | Declining rotation curves of galaxies as a test of gravitational theory[END_REF]), the slight increase in our case is due to the gravitational force of the superfluid DM itself, on top of the MONDian phonon-mediated force. A precise assessment of the rise of observed rotation curves of HSB galaxies towards the edge of their disk would be a subtle but clean way to distinguish our model from traditional MOND.

As a final note, we point out that the BTFR in the superfluid DM context is thus tight and MONDian by construction, but that statistical galaxy-galaxy lensing observations, probing larger distances in the halo, should in principle not give the same "BTFR", but rather give results closer to those expected in ΛCDM, especially if the stellar to halo mass relation remains the same.

To summarize, the above fiducial galaxies illustrate that the predicted rotation curves of spiral galaxies in the superfluid framework closely resemble those of MOND, both for LSB and HSB galaxies. The main difference is limited to the rise of HSB galaxy rotation curves towards the edge of their observed disks.

VIII. GALAXY CLUSTERS

In this Section, we make a few basic observations about superfluid DM in galaxy clusters [START_REF] Hodson | Galaxy Clusters in the Context of Superfluid Dark Matter[END_REF]. For this purpose it suffices to revert to the back-of-the-envelope approach of Sec. III.

Because of the larger velocity dispersion in galaxy clusters, it is easy to show that for the parameter values of interest the Bose enhancement factor N given in ( 15) is small at distances comparable to R 200 . Thermalization cannot be achieved at such large distances, thus galaxy clusters are mostly in the normal (collisionless) phase. Nevertheless, as we go towards the central region and density increases, DM particles can become degenerate and thermalize efficiently.

Unlike galaxies, clusters are better off, observationally speaking, without a significant superfluid component. Taking into account that the thermalized part of the halo is expected to be highly spherical, it is reasonable to demand

R T R 200 < ∼ 0.1 . (50)
Equivalently, substituting for R 200 using [START_REF] Milgrom | A Modification of the Newtonian dynamics: Implications for galaxies[END_REF] we can write this condition as

R T < ∼ 200 M 10 15 M 1/3 kpc . (51)
As in Sec. III, we use the criterion Γ t dyn ∼ 1 as an estimator for the size of the superfluid core. For such relatively small R T , however, we cannot assume that the phase transition happens in the 1/r 3 tail of the NFW profile. We must instead use an actual NFW profile. On the other hand, because the superfluid region is constrained to be small, we can safely assume that the NFW profile matches the ΛCDM expectation. Specifically, for a M = 10 15 M cluster (corresponding to R 200 2 Mpc), we assume c = 4 as concentration parameter, consistent with the mass-concentration relation in ΛCDM cosmology [START_REF] Dutton | Cold dark matter haloes in the Planck era: evolution of structural parameters for Einasto and NFW profiles[END_REF].

Using this NFW density profile, and the corresponding velocity profile, we calculate Γ t dyn as a function of r, starting from R 200 and moving inwards. The thermal radius is estimated as in Sec. III to be the location when Γ t dyn = 1. Substituting the result into (50), we obtain a lower bound on the DM mass

m 2.7 σ/m cm 2 /g 1/4
eV .

Combined with the upper bound ( 24) derived from galaxies, we obtain the allowed range:

2.7 eV < ∼ m σ/m cm 2 /g -1/4 < ∼ 4.2 eV . (53)
Thus with σ/m around the upper limit cm2 /g from merging clusters [START_REF] Randall | Constraints on the Self-Interaction Cross-Section of Dark Matter from Numerical Simulations of the Merging Galaxy Cluster 1E 0657-56[END_REF][START_REF] Massey | Cluster Bulleticity[END_REF][START_REF] Harvey | The non-gravitational interactions of dark matter in colliding galaxy clusters[END_REF][START_REF] Wittman | The Mismeasure of Mergers: Revised Limits on Selfinteracting Dark Matter in Merging Galaxy Clusters[END_REF], m must be in the few eV range. 2 The fiducial values m = eV and σ/m = 0.01 cm 2 /g adopted in fitting the galaxy rotation curves fall within this range. For instance, they imply R T 163 kpc for the aforementioned cluster, which satisfies [START_REF] Khoury | An Alternative to Particle Dark Matter[END_REF].

IX. SUPERFLUID EXTERNAL FIELD EFFECT

An important aspect of the MOND phenomenology is the so-called external field effect (EFE), which is a consequence of the breaking of the strong equivalence principle. In GR, the equivalence principle tells us that a homogeneous acceleration has no physical consequence -it can be removed by moving to the freely-falling elevator. This is not so in MOND. The internal dynamics of a subsystem are affected by the presence of an external acceleration. See [START_REF] Sanders | Modified Newtonian dynamics as an alternative to dark matter[END_REF][START_REF] Famaey | Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions[END_REF] for reviews.

For instance, consider a subsystem with low internal acceleration (a int a 0 ) in the background of a large external acceleration (a ext a 0 ). Then the EFE implies a screening effect, resulting in approximately Newtonian dynamics in the subsystem. Consider then the situation of an external-field dominated subsystem (a int a ext ) immersed in a MONDian background acceleration (a ext a 0 ). In this case, the EFE implies that the subsystem will remain Newtonian, but with an enhanced Newton's constant G eff ∼ a0 a G N . The EFE is not special to MOND but is an example of a more general phenomenon in scalar field theories with kinetic interactions known as kinetic screening [START_REF] Babichev | Improving relativistic MOND with Galileon k-mouflage[END_REF][START_REF] Babichev | k-Mouflage gravity[END_REF][START_REF] Brax | Screening fifth forces in k-essence and DBI models[END_REF][START_REF] Burrage | Screening of scalar fields in Dirac-Born-Infeld theory[END_REF]. In P (X) theories, of which (3) is an example, non-linearities in the scalar field triggered by ∇φ (the acceleration) becoming large can result in the suppression of the scalar field effects and the local recovery of GR. See [START_REF] Joyce | Beyond the Cosmological Standard Model[END_REF] for a review.

The main difference in the superfluid context is that the MONDian effective theory (3), or more precisely its finite-temperature version [START_REF] Mcgaugh | The Baryonic Tully-Fisher Relation of Gas Rich Galaxies as a Test of LCDM and MOND[END_REF], is only valid within the superfluid core. This implies, first of all, that the phononmediated force between two bodies only applies if both bodies reside within the superfluid region. If one body is inside while the other is outside, the only force acting on them is gravity. The same holds true for the EFE -it is only active within the superfluid core of the host.

Below we apply these observations to various subsystems in galaxies and in galaxy clusters.

A. Satellite galaxies and globular clusters

Based on the above discussion, we expect that within the superfluid core of the Milky Way, both globular clusters and satellite galaxies are expected to follow MOND predictions. In particular, the tidal streams of globular clusters, such as Pal 5, orbiting within the superfluid core should be affected by the EFE just as in traditional MOND [START_REF] Thomas | Stellar streams as gravitational experiments II. Asymmetric tails of globular cluster streams[END_REF]. The case of dwarf spheroidal galaxies depends on the actual exact size of the superfluid core around the Milky Way. With a core of ∼ 60 kpc, most dwarf spheroidals would sit outside of the core, and should thus not be affected by the phonons of the superfluid core. In this case, if the dwarf spheroidals are primordial (hence not tidal dwarfs), the expected DM masses from abundance matching being very large, they would mostly be expected to display higher velocities than predicted in MOND, even without the EFE, as observed [START_REF] Spergel | Dwarf spheroidal galaxies and non-Newtonian gravity[END_REF][START_REF] Milgrom | MOND and the seven dwarfs[END_REF][START_REF] Angus | Dwarf Spheroidals in MOND[END_REF][START_REF] Lüghausen | A census of the expected properties of classical Milky Way dwarfs in Milgromian dynamics[END_REF]. Also, globular clusters at large distances from the Milky Way, such as NGC 2419 [START_REF] Ibata | The globular cluster NGC 2419: a crucible for theories of gravity[END_REF] or Pal 14 [START_REF] Jordi | Testing Fundamental Physics with Distant Star Clusters: Analysis of Observational Data on Palomar 14[END_REF], would be expected to be Newtonian as long as they do not harbor their own DM halo.

In the case of a larger superfluid core of size ∼ 120 kpc, dwarfs such as Sextans and Draco would sit well within the core, and would thus need to be explained by a combination of outliers, binaries, extreme anisotropy, and out-of-equilibrium dynamics [START_REF] Serra | Implications for dwarf spheroidal mass content from interloper removal[END_REF][START_REF] Mcgaugh | Local Group Dwarf Spheroidals: Correlated Deviations from the Baryonic Tully-Fisher Relation[END_REF]. Globular clusters such as Pal 14 and NGC 2419 would moreover become MONDian. Such a large core would have the advantage, however, to naturally explain the observed velocity dispersion of Crater II [START_REF] Mcgaugh | MOND Prediction for the Velocity Dispersion of the Feeble Giant Crater II[END_REF][START_REF] Caldwell | Crater 2: An Extremely Cold Dark Matter Halo[END_REF], which is surprisingly low in the standard context.

The same reasoning applies to tidal dwarf galaxies resulting from the interaction of massive spiral galaxies. Those are not expected to harbor a significant DM halo, and thus should be Newtonian [START_REF] Lelli | Gas dynamics in tidal dwarf galaxies: Disc formation at z = 0[END_REF][START_REF] Lelli | Tidal Dwarf Galaxies: Disc Formation at z = 0[END_REF] as long as they are located outside the superfluid core of their host.

Finally, tidal stellar streams of satellite galaxies such as the Sgr dwarf would be expected to display very peculiar features linked to the crossing of the transition radius between the superfluid and normal phase. Once reaching this transition radius, they could be ejected to much larger distances and then come back within the superfluid core on a different orbit. This could leave distinct signatures in tidal streams, such as bifurcations, which would be unexplainable in any other framework. We leave a detailed study of tidal streams to future work.

B. Ultra-diffuse galaxies in galaxy clusters

As shown in Sec. VIII, galaxy clusters have a relatively small superfluid core, of order 100 kpc for m = 1 eV. Nevertheless, galaxies within galaxy clusters are still expected to be surrounded by their own superfluid cores. This is interesting, for instance, in view of the recent observations of ultra-diffuse galaxies within clusters [START_REF] Van Dokkum | Forty-Seven Milky Way-Sized, Extremely Diffuse Galaxies in the Coma Cluster[END_REF][START_REF] Koda | Approximately A Thousand Ultra Diffuse Galaxies in the Coma cluster[END_REF], which are observed to be extremely DM dominated [START_REF] Van Dokkum | Spectroscopic Confirmation of the Existence of Large, Diffuse Galaxies in the Coma Cluster[END_REF]. This class of objects poses a problem for MOND because of the EFE [START_REF] Hodson | Galaxy Clusters in the Context of Superfluid Dark Matter[END_REF][START_REF] Milgrom | Ultra-diffuse cluster galaxies as key to the MOND cluster conundrum[END_REF]. Indeed, although ultra-diffuse galaxies have small internal accelerations, and thus in isolation should be deeply MONDian, the large external acceleration from the galaxy cluster should cause them to be Newtonian, in conflict with observations.

In the superfluid case, however, there is no EFE in those galaxies since they are clearly orbiting outside of the superfluid core of the cluster itself. Moreover, these galaxies are very extended, meaning that their baryon fraction must be smaller than the one assumed here for spiral galaxies, in order for their superfluid core to encompass their large observed size. Hence, these LSB galaxies would share some characteristics of the HSB galaxy analyzed above in terms of the mass of their DM halo. This means that they should slightly deviate from the BTFR, with velocities slightly above the BTFR expectation.

X. DYNAMICAL FRICTION

A distinctive prediction of the superfluid framework is the absence of dynamical friction for subsonic motion within the superfluid region [START_REF] Hui | Ultralight scalars as cosmological dark matter[END_REF][START_REF] Ostriker | Dynamical friction in a gaseous medium[END_REF]. This may alleviate a number of minor problems for ΛCDM. For instance, instead of being slowed down by dynamical friction, galactic bars in spiral galaxies should achieve a nearly constant velocity [START_REF] Goodman | Chance and Chandra[END_REF], as favored by observations [START_REF] Debattista | Dynamical friction and the distribution of dark matter in barred galaxies[END_REF].

It may also offer a natural explanation to the longstanding puzzle of why the five globular clusters orbiting Fornax have not merged to the center to form a stellar nucleus. Indeed, in the context of CDM, dynamical friction should have caused the globular clusters to rapidly fall towards the center of Fornax [119,120]. In reality Fornax shows no sign of such mergers. See [START_REF] Cole | The mass distribution of the Fornax dSph: constraints from its globular cluster distribution[END_REF][START_REF] Sanchez-Salcedo | An extensive study of dynamical friction in dwarf galaxies: the role of stars, dark matter, halo profiles and mond[END_REF] for possible explanations within the classical DM particles context. In our case, the globular clusters should be happily swimming within the superfluid core without dissipation.

As another speculation, the suppression of dynamical friction might shed light on the luminous red galaxy (LRG) two-point correlation function [START_REF] Masjedi | Very small-scale clustering and merger rate of luminous red galaxies[END_REF], which remains surprisingly featureless (ξ(r) ∼ 1/r 2 ) down to small distances where one would expect dynamical friction and mergers to show up. This observation in fact implies an upper bound for the LRG-LRG merger rate of < ∼ 0.6 × 10 4 Gyr -1 /Gpc 3 [START_REF] Masjedi | Very small-scale clustering and merger rate of luminous red galaxies[END_REF]. The suppression of dynamical friction within the superfluid cores of the LRG may partially explain these observations. Two other speculations related to reduced dynamical friction are that it could have interesting consequences on the long term evolution of Hickson compact groups [START_REF] Kroupa | Galaxies as simple dynamical systems[END_REF], provided that their cores are large enough, as well as on the possible history of the Local Group, allowing for a past Milky Way-Andromeda encounter [START_REF] Zhao | Local Group timing in Milgromian dynamics. A past Milky Way-Andromeda encounter at z ¿ 0.8[END_REF][START_REF] Banik | Dynamical history of the Local Group in ΛCDM[END_REF]. The latter may require a large superfluid core for the Milky Way. By reducing (39) down to 0.85 eV, which saturates the bound acquired for clusters, one can increase the superfluid core radius to ∼ 100 kpc. One could push it further to even greater values by violating (52); however the scattering cross section does not have to be a constant, in fact the mild velocity dependence would suffice to push the core radius well-above 100 kpc for a HSB galaxy, while simultaneously satisfying (52) in clusters. Moreover, as we already mentioned, if we set the total DM mass of HSB to the ΛCDM abundance matching value then the superfluid core radius is already ∼ 100 kpc without even altering [START_REF] Hees | Combined Solar System and rotation curve constraints on MOND[END_REF].

In ongoing work, we are currently analyzing quantitatively the effect of dynamical friction in the superfluid context [START_REF] Berezhiani | Dynamical friction in superfluid dark matter[END_REF].

XI. CONCLUSION

The theory of DM superfluidity proposed in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] has some commonalities with SIDM and fuzzy DM [START_REF] Hui | Ultralight scalars as cosmological dark matter[END_REF]128]. All three proposals achieve a cored DM profile in central regions of galaxies, either through interactions or quantum effects, to alleviate existing observational puzzles/tensions with ΛCDM on galactic scales. In all three proposals, the DM outside the core is in approximately collisionless form and assumes an NFW profile. Like the quantum pressure of fuzzy DM, the classical pressure of superfluid DM results in lower central densities and implies a minimal halo mass necessary for collapse and virialization.

The main difference is that DM superfluidity achieves a much larger core, encompassing the entire range of scales probed by rotation curve observations. Moreover, the superfluid collective excitations (phonons) within the core mediate a long-range force within the core, thereby affecting the dynamics of orbiting baryons and reproducing the MOND phenomenology. Unlike most attempts to modify gravity, there is no fundamental additional long-range force in the model. Instead the phonon-mediated force is an emergent phenomenon which requires the coherence of the underlying superfluid substrate.

This has important implications for the phenomenological viability of the scenario vis-à-vis solar system tests of gravity. Although typical accelerations in the solar system are large compared to a 0 , post-Newtonian tests are sensitive to small corrections to Newtonian gravity and require further suppression of the new force at short scales. The superfluid framework offers an elegant explanation. As argued in [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF], the large phonon gradient induced by an individual star results in a breakdown of superfluidity in its vicinity. Coherence is lost within the solar system. Individual DM particles still interact with baryons, but no long-range force can be mediated.

In this paper, using a combination of analytical and numerical arguments, we developed a number of novel distinctive predictions of the superfluid framework. For this purpose we derived an approximate finite-temperature density profile, consisting of a superfluid core, with approximately homogeneous density, surrounded by an envelope of normal-phase DM particles following an NFW profile. In this context, we constrained the parameters of the theory from various observational constraints. In particular, we constrained the mass m of the particles such that the superfluid core can be large enough in galaxies, but small enough in clusters. We constrained the combination Λm 3 such that massive enough halos could form.

An important difference of the present analysis compared to the original papers [START_REF] Berezhiani | Dark Matter Superfluidity and Galactic Dynamics[END_REF][START_REF] Berezhiani | Theory of dark matter superfluidity[END_REF] is that the superfluid core makes up only a modest fraction of the entire halo. It is large enough to encompass the observed rotation curves, since the phonon force is critical for reproducing MOND. But it is small enough that most of the mass lies in the approximately collisionless envelope, resulting in triaxial halos near the virial radius.

Since, as we showed, keeping unaltered the relation between stellar mass and halo mass with respect to ΛCDM abundance matching leads to acceptable rotation curves, galaxy-galaxy lensing statistics should be mostly unaltered too. In particular, contrary to MOND, one does not expect a fine relation r E ∝ √ M b between the baryonic mass and the weak lensing signal. Also, both photons and gravitons travel at the speed of light along the same geodesics, consistent with the tight constraint recently established by the neutron star merger GW170817.

The main result of this paper is the explicit rotation curve fitting of two disk galaxies: a representative LSB galaxy IC 2574, and a representative HSB galaxy UGC 2953. The superfluid model offers an excellent fit in both cases. This is particularly remarkable for IC 2574, which is notoriously difficult to recover in ΛCDM cosmology with standard feedback recipes [START_REF] Oman | The unexpected diversity of dwarf galaxy rotation curves[END_REF]. The predicted rotation curve for IC 2574 is phonon-dominated and therefore closely matches the MOND prediction.

The rotation curve for UGC 2953 shows a small but potentially interesting distinguishing feature from traditional MOND. In contrast with MOND, where the rotation curve becomes flat (or slightly decreases, because of the EFE [START_REF] Hees | Combined Solar System and rotation curve constraints on MOND[END_REF][START_REF] Haghi | Declining rotation curves of galaxies as a test of gravitational theory[END_REF]), we found a slight rise in the asymptotic velocity. This is due the gravitational pull of the superfluid DM mass itself, which adds to the phonon acceleration. A precise assessment of the rise of observed rotation curves of HSB galaxies towards the edge of their disk would be a subtle but clean way to distinguish our model from traditional MOND.

We also explored the observational consequences, albeit at a more qualitative level, for other systems, such as galaxy clusters, dwarf satellites, globular clusters and ultra-diffuse galaxies in clusters. A key difference compared to traditional MOND is that the MONDian force in our case is limited to the superfluid region. Both bodies must be in the superfluid core to experience a long-range phonon force.

Similarly, the EFE (or kinetic screening effect) of MONDian phenomenology only applies within the superfluid core. As a result, globular clusters and satellite galaxies within the superfluid core of the Milky Way should follow MOND predictions. But globular clusters at large distances from the Milky Way should be Newtonian, assuming they do not harbor their own DM halo, while dwarf spheroidals outside of the core might display larger velocities than expected from the BTFR if the usual abundance matching recipes apply. On the other hand, tidal dwarf galaxies resulting from the interaction of massive spiral galaxies should be expected to be devoid of DM and purely Newtonian if they sit outside of the core of their merged host. Finally, ultra-diffuse galaxies in galaxy clusters are immune to the EFE in our case, since they are clearly orbiting outside of the (small) superfluid core of the cluster itself. All these predictions are summarized in Table I. They will of course have to be backed by numerical simulations of superfluid DM, which are currently underway [START_REF] Elder | Numerical simulations of dark matter superfluidity[END_REF].

  illustrates the result (solid black), normalized to the central density, for M b = 10 10 M and L = 2 kpc. The value of the central density (or gravitational potential) in this case is fixed to achieve M DM = 10 M b . For the theory parameters, we assume the fiducial values (39) and (40). The resulting transition radii (shown as dashed vertical lines on the Figure) are R NFW = 52 kpc ; R T = 60 kpc ,

FIG. 1 :

 1 FIG. 1: Comparison of our DM density profile (solid black curve) for the toy baryon exponential profile (41) with M b = 10 10 M and L = 2 kpc. The parameters are m = 1 eV, Λ = 0.05 meV and MDM = 10 M b = 10 11 M . The transition radii (dashed vertical lines) are RNFW = 52 kpc and RT = 60 kpc. For comparison, also plotted are the "pure" (T = 0, no baryons) superfluid profile (blue dashed-dotted curve), and the NFW profile (red dashed curve) for the same MDM and concentration parameter c = 6.5.

Figure 2

 2 compares our predicted radial acceleration curve (solid black curve) with that of the MOND simple interpolating function[START_REF] Ho | Cold Dark Matter with MOND Scaling[END_REF] (red dashed curve) for the

  FIG. 2:Comparison of our predicted radial acceleration curve (solid black) for the toy baryon exponential profile[START_REF] Ibata | The globular cluster NGC 2419: a crucible for theories of gravity[END_REF]. The parameters are the same as in Fig.1. The critical acceleration is a0 = 0.87 × 10 -8 cm/s 2 . The MOND prediction (red dashed curve) with the simple interpolating function[START_REF] Ho | Cold Dark Matter with MOND Scaling[END_REF] and a MOND

  FIG. 3: Superfluid DM density profile (solid black curve) for IC 2574 (M b ∼ 2 × 10 9 M ). The parameters are m = 1 eV, Λ = 0.05 meV and MDM = 10 M b . The transition radii (dashed vertical lines) are RNFW = 40 kpc and RT = 49 kpc, compared to R200 = 57 kpc for this mass ratio. For the NFW envelope/profile, we have the concentration parameter c = 5.7. The conventions are the same as in Fig. 1.

2 .

 2 The DM acceleration a DM is the result of our numerical analysis for a simplified sphericallyapproximated baryon distribution ρ spherical b .

1 ×

 1 FIG. 5: Superfluid DM density profile (solid black curve) for UGC 2953 (M b 1.6 × 10 11 M ), with MDM = 10 M b . Other parameters are the same as in Fig. 3. The transition radii (dashed vertical lines) are RNFW = 76 kpc and RT = 82 kpc, compared to R200 = 245 kpc for this mass ratio. For the NFW envelope/profile, we have the concentration parameter c = 5.4. The conventions are the same as in Fig. 1.

FIG. 6 :

 6 FIG.6: The rotation curve of UGC 2953, with data (orange points) taken from[START_REF] Noordermeer | The mass distribution in early-type disk galaxies: Declining rotation curves and correlations with optical properties[END_REF]. The parameter values and conventions are the same as in Fig.4. The only difference is the red dashed curve, which in this case assumes MDM = 65 M b , as determined by ΛCDM abundance matching.

TABLE I :

 I Summary of observational consequences of superfluid DM

	System	Behavior
	Rotating Systems
	Solar system	Newtonian
	Galaxy rotation curve shapes	MOND (+ small DM component making HSB curves rise)
	Baryonic Tully-Fisher Relation	MOND for rotation curves (but particle DM for lensing)
	Bars and spiral structure in galaxies	MOND
	Interacting Galaxies
	Dynamical friction	Absent in superfluid core
	Tidal dwarf galaxies	Newtonian when outside of superfluid core
	Spheroidal Systems
	Star clusters	MOND with EFE inside galaxy host core -Newton outside of core
	Dwarf Spheroidals	MOND with EFE inside galaxy host core -MOND+DM outside of core
	Clusters of Galaxies	Mostly particle DM (for both dynamics and lensing)
	Ultra-diffuse galaxies	MOND without EFE outside of cluster core
	Galaxy-galaxy lensing	Driven by DM enveloppe =⇒ not MOND
	Gravitational wave observations

Note that we will be working mostly hereafter in natural units, where the gravitational potential (∝ c

) is dimensionless, where mass, energy, and acceleration are expressed in eV, and where length and time are in eV -1 . When dealing with rotation curves or halo masses, we will of course switch back to physical units. In natural units, m and Λ can be expressed in eV, while α and β are dimensionless.

An important caveat is that (53) assumes a constant σ for simplicity. Allowing for velocity dependence should weaken the upper bound, thereby broadening the allowed range.

Acknowledgements

We thank Lam Hui, Andrey Kravtsov, Tom Lubensky, Jerry Ostriker, Riccardo Penco and HongSheng Zhao for helpful discussions. The work of L.B. was partially supported by US Department of Energy grant de-sc0007968 at Princeton University. B.F. acknowledges the financial support from the "Programme Investissements d'Avenir" (PIA) of the IdEx from the Université de Strasbourg. J.K. is supported in part by NSF CAREER Award PHY-1145525, US Department of Energy (HEP) Award de-sc0017804, NASA ATP grant NNX11AI95G, and the Charles E. Kaufman Foundation of the Pittsburgh Foundation.