© T. Mistele et al. 2022

1. Introduction

In 2015, Berezhiani & Khoury (2015) proposed a new hypothesis that combines features of cold dark matter (CDM) and modified Newtonian dynamics (MOND; Milgrom 1983a,b,c; Bekenstein & Milgrom 1984): superfluid dark matter (SFDM). In SFDM, dark matter is composed of a light (on the order of eV) scalar field that can condense to a superfluid. In the superfluid phase, phonons mediate a force that is similar to the force of MOND. This hypothesis has since passed several observational tests (Berezhiani et al. 2018; Hossenfelder & Mistele 2019, 2020).

However, it was recently found that SFDM needs about 20% less baryonic mass than MOND to fit the Milky Way rotation curve at R ≲ 25 kpc (Hossenfelder & Mistele 2020). Though a modest effect, this underestimates the stellar mass required by microlensing (Wegg et al. 2017). It also underestimates the amplitude of the spiral structure required to reconcile the Galactic rotation curve measured independently by stars and gas (McGaugh 2019). This offset is similar to that found for emergent gravity (Verlinde 2017) by Lelli et al. (2017a), which shares some properties of SFDM, thus raising the prospect that it might be a general trend. Consequently, SFDM may require a systematically smaller stellar mass-to-light ratio (M/L_*) than MOND. Since MOND generally agrees with the M/L_* expected from stellar population synthesis (SPS) models (McGaugh 2004), such a systematic trend can be problematic for SFDM. To investigate this, we fitted SFDM to the Spitzer Photometry and Accurate Rotation Curves (SPARC) data (Lelli et al. 2016) with M/L_* as a fitting parameter.

2. Models

Four parameters are required for SFDM; we used the fiducial values from Berezhiani et al. (2018), m = 1 eV, Λ = 0.05 meV, α = 5.7, and β = 2. We kept those parameters fixed during our analysis. In Appendix D.2.6, we argue that our conclusions are generally robust against variations in these parameters.

The total acceleration inside the superfluid core of a galaxy is a_tot = a_θ + a_b + a_SF, where a_θ is the acceleration created by the phonon force, a_SF the acceleration stemming from the normal gravitational attraction of the superfluid, and a_b that stemming from the mass of the baryons. The position dependence of those accelerations is determined by the SFDM equations of motion and the distribution of baryonic mass. At a transition radius where the superfluid condensate is estimated to break down, one matches the superfluid core to a Navarro-Frenk-White (NFW) halo (Berezhiani et al. 2018).

From integrating the standard Poisson equation including the superfluid’s energy density ρ_SF as a source term, one obtains , where μ_nr is the chemical potential and ϕ_N(x) is the Newtonian gravitational potential. The gradient of ϕ_N(x) gives a_b + a_SF. In the so-called no-curl approximation, one obtains the phonon force, a_θ, as an algebraic function of a_b and ε_*(x) (see Appendix A.1),

(1)

where M_Pl is the Planck mass (it enters through Newton’s constant). The quantity ε_*(x) controls how closely SFDM resembles MOND. We refer to |ε_*|≪1 as the MOND limit and to |ε_∗| = 𝒪(1) as the pseudo-MOND limit. In the MOND limit and assuming the no-curl approximation, the phonon force points into the same direction as a_b with magnitude . Here, a₀ = α³Λ²/M_Pl is the acceleration scale below which the phonon force becomes important compared to a_b. At larger accelerations, it is subdominant. This gives a typical MOND-like total acceleration a_tot, at least as long as a_SF stays negligible. Usually, a_SF is indeed negligible in the proper MOND limit |ε_*|≪1 but less so in the pseudo-MOND limit |ε_*|=𝒪(1). More details on the definition and rationale behind these limits are in Appendix A (see also Mistele 2021). The actual value of ε_* depends on the baryonic mass distribution and a boundary condition needed to solve the equations of motion.

A big advantage of MOND is that galactic scaling relations such as the radial acceleration relation (RAR; Lelli et al. 2017b) arise automatically with no intrinsic scatter. The same goes for SFDM in the MOND limit |ε_*|≪1. In this limit, SFDM predicts a tight RAR irrespective of the precise value of the boundary condition. This is different outside the MOND limit where the total acceleration a_tot depends sensitively on the choice of the boundary condition. Thus, outside the MOND limit, scaling relations such as the RAR can arise only by adjusting this boundary condition separately for each galaxy. Otherwise, increased scatter and systematic deviations are likely.

In principle, it might be possible that galaxy formation selects the right boundary condition for each galaxy to produce a tight RAR even outside the MOND limit. However, then SFDM loses one of its main advantages over ΛCDM and one might as well stick with ΛCDM.

We compared SFDM to MOND using one of the standard interpolation functions (Lelli et al. 2017b),

(2)

where y = |a_b|/a₀ and a₀ is the one free parameter in MOND. In SFDM the interpolation function is slower to reach its limits for large and small y. Also, usually a₀ is chosen smaller in SFDM compared to MOND to account for the presence of a_SF (Berezhiani et al. 2018). For MOND, we adopt from Lelli et al. (2017b). For SFDM, the fiducial parameters from Berezhiani et al. (2018) give .

To check how sensitive our results are to the particular theoretical realization of SFDM, we included the two-field model from Mistele (2021). In this two-field model, the phenomenology on galactic scales is similar to standard SFDM, but it has the advantages that (a) it does not require ad-hoc finite-temperature corrections for stability, (b) its phonon force is always close to its MOND-limit, and (c) the superfluid can remain in equilibrium much longer than galactic timescales. Both models are described in more detail in Appendix A.

3. Data

We took the observed rotation velocity V_obs directly from SPARC (Lelli et al. 2016). To find the best SFDM fit, we then needed the baryonic energy density ρ_b(R, z) because it is a source for the equation of motion of the superfluid. For this, we used updated high-resolution mass models including resolved gas surface density profiles for 169 of the 175 SPARC galaxies (Lelli 2021, priv. comm.). We excluded the six galaxies that lack radial profiles for the gas distribution.

These mass models provide surface densities Σ for the bulge, the stellar disk, and the HI disk of each galaxy for a discrete set of positions starting at R = 0. We linearly interpolated the data points and assumed zero surface density outside the outermost surface density data point. This gives a simple, data-compatible approximation for the density distribution at all radii.

For the bulge, we assumed spherical symmetry and extracted its energy density from its surface density using an Abel transform,

(3)

For the stellar disk, we assumed a scale height, h_*, of (Lelli et al. 2016)

(4)

where R_disk is the disk scale length from SPARC. Again, we used a linear interpolation of the SPARC surface brightness data points.

For the gas disk we did the same as for the stellar disk, except that in this case we assumed a fixed scale height, h_g = 0.130 kpc. This is the same scale height used in Hossenfelder & Mistele (2020). We do not expect this choice of scale height to significantly affect the results. To account for the non-HI gas, we multiplied the HI surface density by 1.4 (McGaugh et al. 2020).

4. Method

Solutions of the equations of motion can be parameterized by one boundary condition, ε := ε_*(R_mid), where R_mid := (R_min + R_max)/2, and R_min (R_max) is the smallest (largest) radius with a rotation curve data point. The value of ε quantifies how closely the phonon force resembles a MOND force in the middle of the observed rotation curve.

In our fitting procedure, we kept V_obs and the fiducial model parameters of SFDM fixed, but we allowed a common factor, Q_*, to adjust the stellar disk and bulge M/L_* relative to the nominal stellar population values in the Spitzer [3.6] band (M/L_*)_disk = 0.5 and (M/L_*)_bulge = 0.7 (Lelli et al. 2017b),

(5)

Using this total baryonic energy density, we solved the SFDM equations of motion for different boundary conditions. From that we then obtained the expected rotation curve.

We assume that all rotation curve data points are within the superfluid core; otherwise, rotation curves cannot be automatically MOND-like since the MOND-like phonon force is active only inside the superfluid. In our fits, we took the superfluid phase to end only when its energy density reaches zero. That is, we only required that ε_*(x) (see Eq. (1)) is larger than an algebraic minimum value ε_min everywhere within the superfluid. This minimum value is reached when ρ_SF vanishes and (for the case of β = 2) is given by .

For the best fits, we then checked whether all data points lie within the superfluid core according to a different criterion based on thermal equilibrium. It turns out that for 31 of the 169 galaxies this is not the case. However, this criterion for the value of the transition radius to the NFW halo is quite ad hoc. We therefore do not discard these solutions, though we checked that they do not alter the main conclusions (see also Appendix D.2.7).

Then we compared how well this rotation curve matches with the observed velocities, V_obs, from SPARC. For this, we defined the best fit for each galaxy as that with the smallest χ²,

(6)

Here, N is the number of data points in the galaxy, f = 2 is the number of fit parameters (Q_* and ε), σ_Vobs is the uncertainty on the velocity V_obs from SPARC, V_c(R) is the calculated rotation curve in SFDM, and the sum is over the data points at radius R.

We minimized χ² for

(7)

(8)

In our fit code, we scanned values of log₁₀(Q_*) and log₁₀(ε − ε_min).

In the SPARC data, the Newtonian acceleration due to gas sometimes points outward from the galactic center, not toward it, because of a hole in the HI data, possibly due to a transition from atomic to molecular gas. Usually, such a negative gas contribution is countered by the positive contributions from the stellar disk and the bulge and does not pose a problem. When this is not the case, there is technically no stable circular orbit so we cannot calculate a rotation curve. When this happened, we omitted those data points when calculating χ².

As a cross-check and as a comparison for SFDM, we also fitted the RAR to the SPARC data, that is, we fitted the SPARC data with MOND assuming no curl term and the exponential interpolation function ν_e (Lelli et al. 2017b). In this case, we have only one free fit parameter, Q_*, and consequently, when calculating χ², we set f = 1. We describe our fitting and calculation methods in more detail in Appendix C.

5. Results

The result of our MOND fit is similar to that of Li et al. (2018), which also fitted the RAR to SPARC galaxies. The major difference is that Li et al. (2018) used a Markov chain Monte Carlo (MCMC) procedure with Gaussian priors, while we used a simple parameter scan to minimize χ². We also did not vary distance and inclination and did not separately vary the mass-to-light ratio of the stellar disk and the bulge. As a consequence of this simplified fitting procedure, our distribution of best-fit M/L_* has more outliers and looks less Gaussian than that of Li et al. (2018).

Still, our median best-fit stellar mass-to-light ratios and the best-fit χ² values are similar to those from Li et al. (2018). The median stellar disk M/L_* is 0.39. When we restrict ourselves to galaxies with high quality data (q = 1), this becomes 0.47, very similar to the 0.50 from Li et al. (2018). We show the χ² cumulative distribution function (CDF) in Fig. 1, which is also in reasonable agreement with Li et al. (2018).

In Fig. 2 one sees that some galaxies end up at the minimum stellar mass-to-light ratio allowed in our fitting method, corresponding to Q_* ≈ 0.01. If we do not restrict ourselves to q = 1, this peak at Q_* ≈ 0.01 is even more pronounced. As discussed in Appendix D.1, this is an artifact of our fitting procedure and can be ignored in what follows.

5.1. MOND versus SFDM

Figure 2 shows the best-fit Q_* for the 97 galaxies with q = 1. Contrary to what one might naively expect from the Milky Way result (Hossenfelder & Mistele 2019), the SFDM fits do not have significantly smaller Q_* than the MOND fits. Indeed, the median Q_* for the q = 1 galaxies is about 4% larger than for MOND.

One reason for this is that for many galaxies the superfluid is not in the MOND limit |ε|≪1, as one sees from Fig. 3. We theoretically explain why going outside the MOND limit allows for larger M/L_* in Appendix B.2. To confirm this, we did the fits again but required that the galaxies are in the MOND limit, |ε|< 0.4. For the rationale behind the precise value 0.4, please refer to Appendix D.2.3. The resulting Q_* values are shown in Fig. 4.

As one can see from Fig. 1, the fits with the requirement |ε|< 0.4 are not much worse than those without. The averaged Q_* is now smaller than in MOND; for the q = 1 galaxies, the median stellar disk M/L_* is about 10% smaller than for MOND. This confirms superfluids outside the MOND limit as one reason for the large Q_* values in SFDM (see also Appendix D.2.2).

Another reason why SFDM does not universally give smaller Q_* than MOND is that the best-fit Q_* depends on the type of galaxy. In SFDM, Q_* is systematically smaller for galaxies with relatively large accelerations a_b, but not for those with small accelerations. This can be seen, for example, in the right panel of Fig. 5, which shows the best-fit Q_* of each galaxy in SFDM relative to the best-fit Q_* for MOND as a function of the observed asymptotic rotation velocity, V_flat. A larger V_flat is associated with larger accelerations – this is, where SFDM systematically gives smaller Q_* than MOND. There are similar trends for surface brightness and the gas fraction; both also correlate with the accelerations a_b (see Appendix D.2.4 for more details).

The reason for this trend is that the smaller a₀ value of SFDM makes the acceleration smaller than in MOND. This acceleration is dominant at small a_b, so that SFDM needs more baryonic mass than MOND to get the same total acceleration (at least if a_SF is negligible). This is explained in more detail in Appendix B.1.

This trend in shows not only how SFDM is different from MOND but also how SFDM does not comply with expectations from SPS models. The general idea is that MOND is known to be in good agreement with the expectations of SPS (Sanders 1996; Sanders & Verheijen 1998; McGaugh 2004, 2020), so any systematic trend in is potentially problematic.

In our case, both the MOND and the SFDM fits show increased scatter in Q_* for small galaxies. This is shown in the left panel of Fig. 5. One reason for the increased scatter is that the data for smaller galaxies is generally of lower quality. Another reason is that these galaxies tend to be gas-dominated, in which case adjusting Q_* has only a small effect on the overall fit. Consequently, larger changes in Q_* are needed to impact the fit quality. Indeed, the scatter increases dramatically at precisely the scale where gas typically begins to dominate the mass budget (McGaugh 2011; Lelli 2022).

The SFDM fits show an upward trend in Q_* for small galaxies. This corresponds to the upward trend in shown in the right panel of Fig. 5. The problem is that such trends of the stellar M/L_* with galaxy properties are not expected from SPS models for the late type galaxies that compose the SPARC sample. If anything, we expect the stellar M/L_* to increase with mass (e.g., Bell & de Jong 2001), opposite the sense of the trend found here. Indeed, we utilize the near-infrared Spitzer [3.6] band specifically to minimize variations in the mass-to-light ratio. In the most recent stellar population models of late type galaxies (e.g., Schombert et al. 2019, 2022), accounting for the shape of the stellar metallicity distribution tends to counteract the modest effect of stellar age in the near-infrared, leading to the expectation of a nearly constant M/L_*.

Given our simplistic fitting procedure, Fig. 5, left, alone may or may not be convincing evidence for a systematic trend in Q_*. Still, our fitting procedure is well suited to identify relative differences between MOND and SFDM, and we have a good theoretical understanding of this difference. Thus, we expect the trends in our best-fit to be robust. Since is known to be in good agreement with SPS expectations, we interpret the systematic trend seen in the right panel of Fig. 5 as a good indicator of trends in absolute Q_* revealing a tension between SFDM and SPS.

In Appendix D.3 we show how our fit results illustrate that only the MOND limit of SFDM can reproduce MOND-like galactic scaling relations such as the RAR without having to adjust the boundary condition ε separately for each galaxy.

5.2. Tension with strong lensing

Irrespective of the resulting M/L_* values, there is a price to pay for enforcing the MOND limit in SFDM. A MOND-like rotation curve in the MOND limit |ε_*|≪1 can only be achieved by reducing the acceleration created by the gravitational pull of the superfluid. As a result, the total dark matter mass in those galaxies, , comes out to be quite small. Here, is the dark matter mass within the radius r₂₀₀ where the mean dark matter density drops below ρ₂₀₀ = 200 × 3H²/(8πG) with the Hubble constant H. We adopt H = 67.3 km s⁻¹ Mpc⁻¹.

A small is not a problem for fitting SFDM to the observed rotation curves, but it is a problem if one also wants to fit strong lensing data. Indeed, Hossenfelder & Mistele (2019) find that SFDM requires ratios to fit strong lensing constraints, where M_b is the total baryonic mass. Requiring a rotation curve in the MOND limit |ε_*|≪1 for the SPARC galaxies produces average masses at least an order of magnitude smaller.

To illustrate the problem with strong lensing, we have in Fig. 6 plotted the (logarithm of) the total baryonic and the maximum possible total dark matter mass given our requirement |ε|< 0.4 in comparison to the values found in Hossenfelder & Mistele (2019). For this, we used Q_* = 1 for all SPARC galaxies because the precise stellar mass-to-light ratio is irrelevant here. “Maximum possible” here refers not only to the requirement |ε|< 0.4 but also to uncertainties in how to determine the radius where the superfluid core is matched to an NFW halo: We used the transition radius that gives the largest total dark matter mass (see Appendix D.4 for details).

The best SFDM fits to strong lensing data tend to have M_b ≳ 10¹¹ M_⊙ and . In contrast, despite our generous NFW matching procedure, the SPARC galaxies with M_b > 10¹¹ M_⊙ have when restricted to have rotation curves in the MOND limit |ε|< 0.4. This is a stark contrast.

The SPARC galaxies do not reach baryonic masses quite as large as the lensing galaxies from Hossenfelder & Mistele (2019). But from Fig. 6 it seems clear that the trend goes into the wrong direction: The larger the galaxy, the smaller the maximum possible (given |ε|< 0.4).

The quoted values of for the strong lensing fits may seem high. But at least from a ΛCDM abundance matching perspective, these are actually expected due to the large baryonic masses of the lensing galaxies (Hossenfelder & Mistele 2019). Nevertheless, somewhat smaller ratios may be possible. The fitting procedure of Hossenfelder & Mistele (2019) did not aim to produce small values. It only aimed to simultaneously fit the observed Einstein radii and velocity dispersions of the lensing galaxies. Probably somewhat smaller ratios are possible at the cost of somewhat worse fits of the Einstein radii and velocity dispersions. However, given the size of the discrepancy in Fig. 6, we do not expect that the proper MOND limit of SFDM can reasonably fit these data.

To study this closer, we did another calculation in which we allowed galaxies into the pseudo-MOND limit. Concretely, we redid the maximum calculation with the requirement |ε|< 5. The precise value 5 is again somewhat arbitrary. We explain why this is a pragmatic choice in Appendix D.2.3. We see from Fig. 6 that in the pseudo-MOND-limit galaxies with M_b > 10¹¹ M_⊙ still have smaller total dark matter masses than what is required for strong lensing, although the problem is less severe than in the proper MOND limit. Somewhat worse but still acceptable fits to the strong lensing data might be able to ameliorate this.

The pseudo-MOND limit, however, is unsatisfactory for two reasons. First, it relies sensitively on ad hoc finite-temperature corrections of SFDM that may be unphysical. Second, the pseudo-MOND limit has the disadvantage that the acceleration from the superfluid, a_SF, can be significant. If a_SF is significant, we do not automatically get the MOND-type galactic scaling relations, since then the superfluid boundary condition must be adjusted for each galaxy to get the correct total acceleration. In this case, SFDM loses its advantage over CDM despite the phonon force being close to .

Figure 7 shows the size of a_SF relative to a_b + a_θ at the last rotation curve data point at R = R_max, assuming the maximum total dark matter masses from Fig. 6. Indeed, for the pseudo-MOND limit, a_SF is significant for the galaxies with M_b > 10¹¹ M_⊙ relevant for strong lensing. This is despite SFDM having a very cored dark matter profile. Thus, also with the pseudo-MOND limit, we cannot get MOND-like rotation curves and strong lensing at the same time.

5.3. Two-field SFDM

For two-field SFDM, the Q_*-distribution (Fig. 8) and the corresponding CDF (Fig. 1) are similar to those of standard SFDM. However, the two-field model is constructed so that it is easier for the phonon force to be close to the MOND-like value . For this reason, the best fits for two-field SFDM all have |ε_*|≪1, as expected (Fig. 3). Only for two galaxies (NGC 6789, UGC 0732) does ε_* become larger than 0.1. Its largest value is 0.36 for NGC 6789. That is, the acceleration a_b + a_θ is almost always close to the MOND-like value (see Appendix D.5 for more details).

Thus, two-field SFDM can easily have large dark matter masses and |ε_*|≪1 at the same time. It does not have the same problem with strong lensing as the proper MOND limit |ε_*|≪1 of standard SFDM. Two-field SFDM does, however, still have a problem with strong lensing similar to the pseudo-MOND limit of standard SFDM. Large total dark matter masses imply that the rotation curve receives significant corrections from the superfluid’s gravitational pull a_SF. This is despite two-field SFDM having, like standard SFDM, a very cored density profile. For this reason, large total dark matter masses imply systematically higher rotation curve velocities than MOND.

To illustrate this problem we depict in Fig. 9 the maximum possible total dark matter mass for the two-field model, given the requirement that a_SF is at most 30% as large as a_b  +  a_θ at the last rotation curve data point at R  =  R_max (see Appendix D.6). The scatter in the distribution is smaller in the two-field model because it depends less on the details of the baryonic matter distribution (see Appendix D.6). As one can see, in the two-field model the discrepancy with the lensing data is weaker than for the proper MOND limit |ε_*|  ≪  1 of standard SFDM, but still present. Avoiding this tension with the lensing data would require rotation curves that are even less MOND-like. Whether or not somewhat worse but still acceptable fits to the strong lensing data could ameliorate this problem requires further investigation.

6. Conclusion

We have found that it is difficult to reproduce the achievements of MOND with the models that have so far been proposed for SFDM.

---

Acknowledgments

This work was supported by the DFG (German Research Foundation) under grant number HO 2601/8-1 together with the joint NSF grant PHY-1911909.

References

Appendix A: The models

Here, we introduce both the original SFDM model from Berezhiani & Khoury (2015) and the two-field model from Mistele (2021) in more detail.

A.1. Standard SFDM

In standard SFDM, in an equilibrium superfluid core of a galaxy, the phonon field, θ, is determined by the equation (Berezhiani et al. 2018)

(A.1)

and the field is determined by the Poisson equation

(A.2)

with the superfluid energy density, ρ_SF,

(A.3)

Here, m, Λ, and α are model parameters. The quantity is a combination of the (constant) non-relativistic chemical potential μ_nr and the Newtonian gravitational potential ϕ_N(x). It controls how much the superfluid weighs, depending on a boundary condition (see Appendix C.2). The parameter β parametrizes finite-temperature corrections, which are needed to avoid an instability (Berezhiani & Khoury 2015). The phonon force a_θ is given by

(A.4)

We mainly used the no-curl approximation for the θ equation of motion. That is, for the solution of this equation, which is of the form ∇(g ⋅ ∇θ) = a_b for some g, we assumed g∇θ = a_b. This is a standard approximation in MOND and it works well also for SFDM (Hossenfelder & Mistele 2020).

In the no-curl approximation, the quantity ε_*(x) (see Eq. (1)) is useful. As we will see, it controls how closely SFDM resembles MOND. As discussed in Mistele (2021), we have

(A.5)

where

(A.6)

and where x_β(ε_*) is determined by the cubic equation

(A.7)

That is, a_θ is an algebraic function of a_b and ε_*. This also allows us to write ρ_SF as a function of |a_b| and ε_*,

(A.8)

where

(A.9)

For both a_θ and ρ_SF, we have a prefactor proportional to multiplied by a function that depends on ε_* and β only. For later use, we record the expansion of this second function for small and large values of ε_*. For |ε_*|≪1, we have

(A.10a)

(A.10b)

For ε_* ≫ 1,

(A.11a)

(A.11b)

In general, f_β is a monotonically increasing, concave function of ε_* (see Fig. A.2). The function is not monotonic (see Fig. A.1).

To avoid a negative or imaginary ρ_SF as well as an instability, ε_* must be larger than some minimum value ε_*min. Berezhiani & Khoury (2015) assumed , corresponding to ε_* > 0, but this is not required from their Lagrangian. It is an assumption with unclear justification. Here, we were more generous to the model and allowed to become negative as long as ρ_SF stays positive. The corresponding minimum value of ε_* is determined by ρ_SF = 0, which is equivalent to f_β(ε_*min) = 0. For example,

(A.12)

With Eq. (1), this translates into a minimum value for ,

(A.13)

A.1.1. MOND limit

In SFDM, the total acceleration inside the superfluid core of a galaxy can be written as

(A.14)

where a_θ is the acceleration created by the phonon force, a_SF the acceleration stemming from the normal gravitational attraction of the superfluid, and a_b that stemming from the mass of the baryons.

For SFDM to make sense, one needs the superfluid to at least approximately reproduce MOND rotation curves without being sensitive to the choice of the boundary condition ε. Otherwise, one does not get the observed MOND-like scaling relations without carefully adjusting the boundary condition separately for each galaxy. That is to say, without the MOND limit of SFDM, one might as well use CDM.

Rotation curves in SFDM approximate those in MOND when the total acceleration a_tot approximately has the form . This corresponds to two conditions. First, the phonon force a_θ must be close to . Second, the superfluid’s gravitational pull a_SF must be negligible.

The numerical values of the model parameters and the boundary condition of the Poisson equation for determine in which coordinate-range SFDM approximates MOND for a given baryonic mass distribution. Specifically, the MOND limit corresponds to |ε_*(x)| ≪ 1. In this |ε_*|≪1 limit, both conditions to reproduce MOND rotation curves are automatically fulfilled: The phonon force is close to and the superfluid’s gravitational pull a_SF is negligible. The phonon force is close to because, for |ε_*|≪1, the (no-curl version of) the phonon field equation Eq. (A.1) has the MOND-like form |a_θ|a_θ = a₀a_b. This corresponds to the small-|ε_*| expansion from Eq. (A.10a). We explicitly show that a_SF is negligible (i.e., that the second condition is fulfilled) at the end of this subsection.

However, even when ε_* is of order one, deviations of the phonon force from the MOND form remain within the percent range, at least for β = 2 (see Fig. A.1). It will therefore in the following be handy to define a “pseudo-MOND limit,” ε_* = 𝒪(1). If this condition is fulfilled, the phonon field no longer satisfies a MOND-like equation, but the acceleration a_θ of an isolated¹ galaxy is numerically still relatively close to . One difference to the proper MOND limit |ε_*|≪1 is that now the second condition for having MOND-like rotation curves is not automatically fulfilled. The superfluid’s gravitational pull a_SF can be significant. So the observed scaling relations are fulfilled automatically only if a_SF stays sufficiently small, which needs to be checked separately for each solution.

A different problem with the pseudo-MOND limit is that it depends sensitively on the details of the ad hoc finite-temperature corrections introduced in Berezhiani & Khoury (2015) to avoid an instability. For example, the pseudo-MOND limit works only for β close to 2, as can be seen from Fig. A.1, left. Just as these ad hoc finite-temperature corrections, the pseudo-MOND limit may turn out to be unphysical.

It now remains to show that the superfluid’s gravitational pull a_SF is negligible in the proper MOND limit |ε_*|≪1. For simplicity, we assume a point mass baryonic energy density, ρ_b(x) = M_bδ(x), which gives a_b = GM_b/r². Then, for |ε_*|≪1, we have ρ_SF ∝ 1/r (see Eq. (A.10b)). The superfluid’s mass is then

(A.15)

where

(A.16)

We can now estimate the superfluid’s gravitational pull a_SF compared to a_b + a_θ. Roughly,

(A.17)

with . This ratio a_SF/(a_b + a_θ) can be larger than a fraction δ only at a radius r_δ that satisfies

(A.18)

where we used the fiducial numerical parameters from Berezhiani et al. (2018) for the last equality. That is, assuming the proper MOND limit |ε_*|≪1, the superfluid’s mass becomes important only at radii larger than where rotation curves are measured.

A.1.2. Reaching the proper MOND limit

As already mentioned in Mistele (2021), reaching the proper MOND limit |ε_*|≪1 is not always possible. To avoid a negative ρ_SF there is a minimum value for (see Eq. (A.10a)). Typically, is a decreasing function of galactocentric radius and the Poisson equation Eq. (A.13) tells us that has a derivative of about −GM/r² where M includes both the baryonic and superfluid mass. Using the baryonic mass M_b as a lower bound on M then gives a lower bound on . Roughly, . This translates into a rough lower bound on ε_*,

(A.19)

where we used a_b = GM_b/r² and the fiducial parameter values from Berezhiani et al. (2018).

Thus, small galaxies can easily reach the proper MOND limit |ε_*|≪1 over the whole range where their rotation curve is measured. One just needs to ensure that the superfluid mass is not too large, which is usually possible.

In contrast, larger galaxies sometimes struggle to satisfy the MOND limit condition |ε_*|≪1, even when the superfluid mass is as small as possible.

A.1.3. Naive upper bound on MOND limit dark matter mass

In Appendix A.1.1, we saw that being in the proper MOND limit |ε_*|≪1 restricts the superfluid’s gravitational pull a_SF to be relatively small. Similarly, the MOND limit restricts the total dark matter mass to be relatively small, even if we include the non-condensed phase outside the superfluid core.

To see this, consider a galaxy with a superfluid core in the MOND limit |ε_*|≪1 and, for simplicity, assume a point mass baryonic mass distribution ρ_b = M_bδ(x). Then, the superfluid’s mass is M_SF = M_b ⋅ (r/r_c)² (see Eq. (A.15)). In SFDM one usually assumes that the superfluid ends at some finite radius r_NFW where the superfluid’s density is matched to that of an NFW halo. The total dark matter mass can be calculated from

(A.20)

Here, M_NFW denotes the mass of the NFW halo between the radii r_NFW and r₂₀₀. The NFW halo energy density falls off faster than the superfluid energy density (i.e., faster than 1/r). Thus, M_NFW grows slower than quadratically in r and we have the inequality

(A.21)

That is,

(A.22)

which is equivalent to

(A.23)

Numerically, with the fiducial parameter values from Berezhiani et al. (2018) and H = 67.3 km/(s ⋅ Mpc), this is

(A.24)

This is too little for strong lensing even for very massive galaxies (see Appendix D.4). This upper bound is independent of the matching procedure to the NFW halo.

A.2. Two-field SFDM

Two-field SFDM contains two fields θ₊ and θ₋ instead of just one field θ like standard SFDM (Mistele 2021). Still, in equilibrium only two nontrivial equations must be solved. One for θ₊ that carries the MOND-like force and one for the Newtonian gravitational potential ϕ_N. As in standard SFDM, we write the equations in terms of where , μ_nr is the non-relativistic chemical potential, and m is the mass of the superfluid’s constituent particles. Also, as in standard SFDM, the MOND limit of the phonon force is controlled by a quantity . Thus, we use the same notation ε_* in both models.

This model has two contributions to the superfluid energy density, ρ_SF+ and ρ_SF−. As discussed in Mistele (2021), usually ρ_SF− dominates. For our calculation below, we assumed that this is the case and neglected ρ_SF+. We verified that ρ_SF+ is always at most 5% as large as ρ_SF− at R = R_mid for the best fits.

As in standard SFDM, we can get the phonon force from a no-curl approximation as a function of a_b and ε_*. We use the same notation a_θ as in standard SFDM. The quantity is determined as in standard SFDM just with a different ρ_SF. Concretely, we have in the no-curl approximation for an equilibrium superfluid

(A.25)

and

(A.26)

where r₀ is a parameter of the model.

We used the numerical parameter values from Mistele (2021). That is, r₀ = 50 kpc, a₀ = 0.87 ⋅ 10⁻¹⁰ m/s², and, unless stated otherwise, . The quantity m²/α that enters ε_* is a combination of these, namely .

One difference to standard SFDM is that ε_* is almost always small so that the phonon force a_θ almost always has the MOND-like form . But, in contrast to standard SFDM, a small |ε_*| implies only that , not that a_SF is small. Thus, even for |ε_*|≪1 one must check that a_SF is small in order to get MOND-like rotation curves. Otherwise, the total acceleration will be systematically larger than in MOND.

The energy density ρ_SF reaches zero for in two-field SFDM (i.e. ε_*min = 0).

Appendix B: Comparison to MOND

B.1. Assuming the MOND limit of SFDM

For SFDM in the MOND limit we approximately have , which, in MOND, would correspond to the interpolation function

(B.1)

with y = a_b/a₀.

At baryonic accelerations not much smaller or much larger than a₀ (i.e., y = 𝒪(1)), the additional acceleration from SFDM is significantly larger than what one obtains from standard MOND interpolation functions such as (Lelli et al. 2017b)

(B.2)

It is because of this difference in the interpolation functions that one may naively expect SFDM to require less baryonic mass than standard MOND models, at least in the MOND limit.

This is illustrated in Fig. B.1, top. The total acceleration in SFDM is always larger than in MOND, if both use the same baryonic a_b. At intermediate accelerations (a_b ∼ a₀ ∼ 10⁻¹⁰ m/s²) the difference between MOND and SFDM is significant. This can be countered by making a_b in SFDM smaller, that is to say, by choosing a smaller mass-to-light ratio in SFDM than in MOND.

This discussion so far assumes the same a₀ for both SFDM and MOND. However, in practice, one usually chooses a somewhat smaller value for a₀ in SFDM. Indeed, Berezhiani et al. (2018) chose , while MOND typically requires (Lelli et al. 2017b). The motivation of Berezhiani et al. (2018) to choose a lower value is to take into account a possible effect of the superfluid’s gravitational pull a_SF. Indeed, at small accelerations a_b, the total acceleration in MOND is close to

(B.3)

while in SFDM we have

(B.4)

The smaller value of a₀ in SFDM allows us to get the same total acceleration even with a nonzero a_SF. Numerically, the smaller a₀ value is compensated for when a_SF is about .

Neglecting a_SF, this smaller value of a₀ makes the total acceleration smaller, so it counters the need for less baryonic mass in SFDM. This is illustrated in Fig. B.1, bottom. The smaller a₀ value has the biggest impact at small accelerations a_b. At small accelerations, a_b ≪ a₀, SFDM may even require more baryonic mass than MOND, at least if we neglect a_SF. Indeed, a_SF is usually negligible in the proper MOND limit |ε_*|≪1, as discussed above. Thus, assuming the proper MOND limit, we expect to find systematically smaller M/L_* in SFDM than in MOND for galaxies with large a_b but not for galaxies with small a_b. This is roughly what we find in our fits below (see Appendix D.2.4).

B.2. Caveat: MOND limit

The above discussion applies in the MOND limit of SFDM. Outside this MOND limit, the phonon force does not necessarily have the form and the superfluid’s gravitational pull may not be negligible. For example, at ε_* → ∞, we find that (see Fig. A.1, right). That is, having a large ε_* makes a_θ small. A smaller acceleration may allow for larger baryonic masses. Thus, having galaxies end up at ε_* ≫ 1 is a way to allow for relatively large mass-to-light ratios in our fits.

One might be skeptical of this argument for the following reason. The argument relies on the total acceleration a_tot becoming smaller for ε_* ≫ 1. But this is not necessarily the case. A large ε_* does make a_θ smaller. But it is possible that the decrease in a_θ is compensated for by an increase in a_SF. Indeed, at large ε_*, the superfluid’s energy density scales as

(B.5)

Thus, at fixed a_b, a larger ε_* makes the superfluid heavier and thus a_SF larger. For ε_* → ∞, the acceleration a_SF can become arbitrarily large. Thus, the total acceleration does not become smaller for ε_* → ∞, despite the smaller phonon acceleration a_θ.

Still, in practice there is a significant window of large values of ε_* where the total acceleration does become smaller. To see this explicitly, expand for large large ε_*. Then, roughly, a_tot scales with ε_* as

(B.6)

where we treated ε_* as a constant that we can pull out of a_SF (see Appendix A.1). Thus, at fixed a_b, the total acceleration a_tot decreases as a function of large ε_* as long as

(B.7)

with the MOND radius and r_c as defined in Appendix A.1.1. Numerically, for the fiducial parameter values from Berezhiani et al. (2018),

(B.8)

Thus, the total acceleration is a decreasing function of ε_* for a significant range of large ε_* values so that going to large ε_* is one way to allow for relatively large baryonic masses.

Appendix C: Method

C.1. Data

As already mentioned in Sect. 3, we used the observed rotation velocity V_obs directly from SPARC. We did not allow distance or inclination as a fit parameter, so we did not vary V_obs in our fitting procedure. As also described in Sect. 3, we obtained the baryonic energy density ρ_b(R, z) from the surface densities provided by SPARC,

(C.1)

The fit parameter Q_* parametrizes the stellar mass-to-light ratio. For later use, we numerically solved the Poisson equation

(C.2)

separately for ρ_x ∈ {ρ_gas, ρ_*, ρ_bulge} using the Mathematica code used in Hossenfelder & Mistele (2020). This allows us to quickly get a solution to the Poisson equation sourced by the full ρ_b with arbitrary Q_* by the rescaling

(C.3)

where the quantity is minus the standard Newtonian gravitational potential up to an additive constant (see also the next subsection). The numerical procedure solves the Poisson equation within a sphere with radius r_∞ assuming a z → −z symmetry. We assumed spherically symmetric boundary conditions for . Specifically,

(C.4)

This is reasonable for sufficiently large r_∞. We used r_∞ = 100 kpc except when the SPARC V_obs data extend to radii larger than 100 kpc. Then, we increased r_∞ in steps of 5 kpc until r_∞ was larger than the maximum radius of the V_obs data points.

C.2. SFDM calculation

We assume that each galaxy’s V_obs data points lie within its superfluid core. This is discussed in more detail in Appendix D.2.7. Then, in SFDM, there are two equations for a galaxy in equilibrium inside the superfluid core. One for the phonon acceleration, a_θ, and one for the quantity , which contains the Newtonian gravitational potential (see Appendix A.1).

Even in a fully axisymmetric calculation, one can impose spherically symmetric boundary conditions for the fields and θ at some large radius r_∞ (Hossenfelder & Mistele 2020). The value of θ at r_∞ is inconsequential, so one can choose θ(r_∞) = 0. For , its value μ_∞ at r_∞ is important. It determines the size of the superfluid halo and is a free parameter in the boundary conditions. We used a parameter similar to μ_∞ as a free fit parameter in our fitting procedure.

It is useful to split into a part called sourced only by ρ_b and the rest called . That is, with

(C.5)

(C.6)

We used boundary conditions and . We calculated as described in the previous subsection.

C.2.1. A simple approximation

For our fits, we did not do a fully axisymmetric calculation. Instead, we used an approximation that is much faster to compute. Our approximation mainly consists of using a no-curl approximation for a_θ and assuming spherical symmetry for . As discussed in Appendix A.1, the no-curl approximation means that we get a_θ as an algebraic function of a_b and .

The second part of our approximation is assuming spherical symmetry for in . That is, the part of due to the superfluid’s self-gravity is spherically symmetric. Only the baryonic part produces an axisymmetric . This is a reasonable approximation for the following reason. A fully axisymmetric calculation gives a that is not spherically symmetric only at relatively small radii. At these radii, dominates. At larger radii, even a fully axisymmetric calculation gives a spherically symmetric (Hossenfelder & Mistele 2020). Indeed, we imposed spherically symmetric boundary conditions at larger radii. Only at these larger radii does usually become important.

We calculated from Eq. (C.6), which contains the function . To solve this equation in spherical symmetry, we need to make a choice of which R and z to use in evaluating for each r. The same applies to a_b(R, z), which enters indirectly through a_θ. We chose R = r and z = 0. Different choices may give slightly different results.

For ρ_SF, we used the expression Eq. (A.8) valid in the no-curl approximation. The function f_β(ε_*) in this expression for ρ_SF is known analytically but is relatively slow to evaluate numerically. To speed up the calculation, for a given β, we evaluated f_β as a function of log₁₀(ε_* − ε_*min) on an evenly spaced grid with grid spacing 0.01 and linearly interpolated between the grid points. We used the resulting linear interpolation in our calculation since it is faster to evaluate numerically than the analytical form of f_β.

Below, we refer to this approximation as the “simple” approximation. In Appendix C.2.3, we explicitly demonstrate that it works well using a few example galaxies.

C.2.2. Calculation using the “simple” approximation

We calculated as described above in Appendix C.1. From this, we got a_b as . In accordance with our simple approximation, we used the no-curl approximation for a_θ so that we got a_θ as an algebraic function of a_b and (see Appendix A.1). The remaining part was to calculate . This then also gave a_SF as .

For we assumed spherical symmetry and we used the form Eq. (A.8) for ρ_SF valid in the no-curl approximation. Then, Eq. (C.6) becomes

(C.7)

where

(C.8)

The non-spherically symmetric functions a_b(R, z) and are evaluated at R = r, z = 0.

This is a second-order ODE for . As boundary conditions we chose

(C.9)

(C.10)

With dr = 0, the first boundary condition is a standard regularity condition at the origin. To avoid numerical issues, we chose a small nonzero value for dr, usually dr = 10⁻⁸ kpc. The second condition corresponds to a choice of ε = ε_*(R_mid) (see Eq. (1)). It parametrizes how close to the MOND limit |ε_*|≪1 a galaxy is in the middle of the V_obs data points.

Solutions for are such that they typically reach their minimum allowed value (see Eq. (A.13)) at some finite radius. Beyond this radius, assuming a superfluid core makes no sense. Thus, whenever solutions ended up with at some radius dr ≤ r ≤ R_max, we discarded them, since we assumed all data points lie within the superfluid core.

Sometimes, Mathematica fails to solve the equations for numerical reasons. This is indicated by its “NDSolve” producing a “FindRoot::sszero,” a “NDSolveValue::berr,” or a “NDSolveValue::evcvmit” error that we can check for. In this case, we automatically decreased dr by factor of 100 and retried.

C.2.3. Validating the simple approximation

Here, we explicitly compare the simple approximation described in Appendix C.2.1 against a fully axisymmetric calculation. For definiteness, we used the best-fit Q_* and f_ε* values for SFDM (see Appendix C.3).

For the fully axisymmetric calculation we used the Mathematica code from Hossenfelder & Mistele (2020). This code expects boundary conditions in the form for some r_∞ and μ_∞. We chose r_∞ = 100 kpc unless stated otherwise. Our simple approximation instead uses a value of as a boundary condition. To compare our simple calculation and the fully axisymmetric calculation for the same physical boundary conditions, we first did the simple calculation with the best-fit values for Q_* and f_ε*. We then evaluated the solution from this simple calculation at r = r_∞ and used the resulting value as the boundary condition μ_∞ for the fully axisymmetric calculation.

Our simple calculation makes two main approximations. First, we used the no-curl approximation for the phonon force. Second, we assumed spherical symmetry for . When our simple calculation disagrees with the fully axisymmetric calculation we want to know which of these two parts is responsible for the deviation. To this end, we did a third calculation where we used the fully axisymmetric calculation for , but then used the no-curl approximation when calculating a_θ for the rotation curve. We refer to this as the “full+nocurl” calculation.

In Fig. C.1, left, we show the rotation curve and for NGC 2403 for the different types of calculation described above. The calculations differ by a few percent at intermediate radii. The full+nocurl rotation curve lies pretty much on top of the simple rotation curve, while the “full” rotation curve differs from the two others at intermediate radii. Thus, the no-curl approximation is the source of the this difference between the full and the simple calculations. For , all calculations agree almost perfectly with each other (see Fig. C.1, right).

This same qualitative result holds for DDO 064 shown in Fig. C.2. This is an example of a galaxy that is in the proper MOND limit |ε_*|≪1 almost everywhere at R_min ≤ R ≤ R_max. The no-curl approximation does not always lead to visible deviations between the full and the simple calculations. An example is IC 2574 where the full and simple calculations agree almost perfectly with each other (see Fig. C.3).

Thus, our simple approximation works well, with the main error being due to the no-curl approximation.

C.3. Fitting method

For SFDM, we used the two parameters

(C.11)

for the stellar mass-to-light ratio (see Eq. (C.1)), and

(C.12)

for the superfluid halo (see Eq. (C.12)), as fit parameters. Here, ε_*min is the minimum possible value of ε_* where ρ_SF vanishes (see Eq. (A.12)). We did not vary the model parameters m, α, Λ, and β. We used Mathematica’s “NMinimize” with the “NelderMead” method to find the smallest χ² for each galaxy,

(C.13)

Here, N is the number of data points in the galaxy, f = 2 is the number of fit parameters, σ_Vobs is the uncertainty on the velocity V_obs from SPARC, V_c(R) is the calculated rotation curve in SFDM, and the sum is over the data points at radii R.

We minimized χ² for f_Y and f_ε* in the range Eqs. (7) and (8). When f_ε* is too small, it can happen that does not exist with the desired parameters, as discussed in Appendix C.2.2. In this case, we artificially set χ² = 10¹⁰. Then, NMinimize continued searching elsewhere.

The NelderMead search method is faster than a simple grid search but can get stuck in local minima. To avoid this, we ran NMinimize three times with different starting points. Of the three results, we used that with the smallest χ². The first run is with the “RandomSeed” option set to 0, the second with the RandomSeed option set to 1, and the third run is with the starting points (0, 0), ( − 0.5, 0), and (0, −0.5). The third run is to guarantee that the point f_Y = 0 is visited at least once, since this point corresponds to the M/L_* expected from SPS models.

To further reduce the needed computation time we rounded f_Y and f_ε* to 0.01 before any calculation. For (un-rounded) f_Y and f_ε* that give the same rounded values as a previous calculation, we reused the previous results without a new computation.

This fitting method is much simpler than the MCMC method used in Li et al. (2018). Still, as we will see in Appendix D.1, we found similar results for the stellar M/L_* as Li et al. (2018) for a standard MOND model. In addition, for SFDM we could not set up informative priors on f_ε* anyway since there is so far no cosmology from which to infer such a prior.

In the SPARC data, the Newtonian acceleration due to gas sometimes points outward from the galactic center, not toward it. Usually, such a negative gas contribution is countered by the positive contributions from the stellar disk and bulge such that the total a_b points to the galactic center. But sometimes this is not the case, especially for small f_Y. When this happened, we simply ignored the data points where a_b is negative when calculating χ².

As a cross-check and as a comparison for SFDM, we also fitted the RAR to the SPARC data, that is, we fitted the SPARC data with MOND assuming no curl term and the exponential interpolation function ν_e (Lelli et al. 2017b). We call this the “MOND” model. In this case, we have only one free fit parameter, f_Y. Thus, when calculating χ², we set f = 1. Also, we used the “SimulatedAnnealing” method of Mathematica’s NMinimize function with one run instead of the NelderMead method with three runs. We did not round f_Y to 0.01 for these MOND fits.

Below we consider modifications of both the SFDM model and the MOND model. The SFDM-based models will be fitted as the “SFDM” model. The MOND-based models will be fitted as the MOND model. We will discuss the details of these modifications below.

For the SFDM models, we parametrize the total dark matter within the last rotation curve data point R_max by a parameter, f_MDM,

(C.14)

where is defined by

(C.15)

with and with the SPS M/L_* values for a_b (i.e. M/L_* = 0.5 for the disk and M/L_* = 0.7 for the bulge). The parameter f_MDM measures how far the dark matter mass at R_max is away from the reference value . This reference value is defined such that the associated dark matter acceleration a_SF counters the difference in between MOND and SFDM due to the different choice for a₀ (see Appendix B.1). Thus, f_MDM parametrizes how large the dark matter mass is relative to the mass that cancels this a₀ difference.

C.4. Two-field SFDM calculation

We can use the same calculation and fitting procedure as for standard SFDM. We simply have to adjust the expression for the superfluid energy density and the algebraic no-curl solution of the phonon force. Apart from that, we adjusted the calculation only in two ways that we now explain.

The superfluid energy density, ρ_SF, of two-field SFDM is linear in and depends on no other fields. This allows the calculation to be sped up. For a given galaxy, we first calculated one particular solution, , of the full, inhomogeneous equation as previously described,

(C.16)

To get solutions for the same galaxy with different boundary conditions, we can add solutions of the homogeneous equation to the so-obtained . Since we assume spherical symmetry, the solutions to the homogeneous equation are Asin(r/r₀)/r for arbitrary A. To get a solution for some desired boundary condition, we just needed to choose an appropriate A.

For standard SFDM, we used ε_*(R_mid) in the range ε_*min + 10⁻² to ε_*min + 10⁴ as the boundary condition for . For two-field SFDM, we must adjust this range. This is because, in two-field SFDM, the phonon force can more easily reach the MOND limit |ε_*|≪1 (i.e., typical values of |ε_*| are much smaller). Thus, we changed the range of ε_*(R_mid) values scanned by our fit code to be

(C.17)

We note that ε_*min = 0 in two-field SFDM. We will later see that no galaxies end up at the boundaries of this range, so it seems to be reasonable.

Appendix D: Results

D.1. M/L_* in MOND

Our MOND fit should give results roughly comparable to Li et al. (2018), which also fitted the RAR to SPARC galaxies. A difference is that Li et al. (2018) used an MCMC procedure with Gaussian priors, while we used a simple parameter scan to minimize χ². We also did not vary distance and inclination and did not separately vary the mass-to-light ratio of the stellar disk and the bulge. As a consequence of this simplified fitting procedure, the distribution of best-fit M/L_* has more outliers and looks less like a Gaussian in our case compared to Li et al. (2018). This can be seen for example in Fig. D.2, which shows the histograms for the best-fit f_Y for the galaxies with the SPARC quality flag q = 1.

Still, the median best-fit stellar mass-to-light ratios and the best-fit χ² values are similar to those from Li et al. (2018). The median stellar disk M/L is 0.39. When we restrict the ourselves to q = 1 galaxies, this becomes 0.47. This is shown in Table D.1. This is in reasonable agreement with Li et al. (2018), who obtained 0.50. We show the cumulative χ² distribution in Fig. D.1, which is also in reasonable agreement with Li et al. (2018).

In Fig. D.2, one sees that some galaxies end up at the minimum stellar mass-to-light ratio allowed in our fitting method, corresponding to Q_* ≈ 0.01. If we do not restrict ourselves to q = 1, this peak at Q_* ≈ 0.01 is even more pronounced. These galaxies with Q_* ≈ 0.01 come about as follows. Consider a galaxy where the observed V_obs is smaller than that computed in MOND. The computed rotation curve can be brought closer to V_obs by decreasing M/L_*. It can happen that M/L_* must be reduced so much that the gas component, which is unaffected by M/L_*, dominates. When this happens and when V_obs is still smaller than the computed rotation velocity, the fitting code will continue to decrease M/L_* to improve χ². But of course this will only barely change χ² since the gas component dominates, and the galaxy ends up at the minimum allowed mass-to-light ratio corresponding to Q_* = 0.01. This does not happen in Li et al. (2018), since varying the distance and inclination can avoid such situations and also because the Gaussian priors discourage going to the minimum allowed M/L_*.

Thus, the best-fit mass-to-light ratios of the galaxies at Q_* = 0.01 should not be taken seriously. They are an artifact of our simplified fitting procedure. They have a comparably good fit also with larger Q_*. We verified that the galaxies at Q_* = 0.01 are gas-dominated at their best-fit Q_*. To check that our results do not depend on these outlier galaxies, we also include a column in Table D.1 that averages only over the galaxies where the MOND fit gives log₁₀ Q_* > −1.5. This gives a median stellar disk M/L of 0.44. This lies between the result for the q = 1 galaxies and the one we got when not restricting the galaxies.

In Table D.1 and Fig. D.1, we also show the results for a fourth quality cut we call “thermal ok.” This refers to the galaxies where, in our SFDM fit discussed below, a simple estimate shows that all SPARC data points lie within the superfluid core of the galaxy (see Appendix D.2.7 for more details). Here, we just note that this quality cut does not qualitatively change our results.

We show the mean stellar disk mass-to-light ratios in Table D.2. These differ from the median values for all quality cuts because the resulting f_Y distributions are not Gaussian as already discussed.

D.2. M/L_* in SFDM

For SFDM, we show the χ² CDF in Fig. D.1 and the Q_* and f_MDM histograms for the q = 1 galaxies in Fig. D.2 and Fig. D.3. The χ² CDF and the Q_* histogram look qualitatively similar to those from the MOND fit, just with some numerical differences. For example, as for the MOND fit, there are some galaxies at the minimum value Q_* = 0.01. These are the galaxies that become gas-dominated during the fitting procedure as explained in Appendix D.1. The precise distribution of best-fit f_MDM values should not be taken too seriously, especially at smaller values. This is because the superfluid halo’s Newtonian gravitational pull is often subdominant in SFDM, so that our fitting method cannot really distinguish different f_MDM values, as long as a_SF stays subdominant.

D.2.1. Stellar mass-to-light ratio

Our initial question was whether or not SFDM needs a smaller M/L_* than standard MOND models. We find that this is not necessarily the case. In SFDM, the median stellar disk M/L is 0.49 for the q = 1 galaxies. This is not much smaller than MOND. The numerical details depend on whether one considers the mean or the median and on the chosen galaxy cuts (see Tables D.1 and D.2). Still, a robust finding across all of these choices is that SFDM does not give a significantly smaller Q_* than MOND.

There are two reasons for this. First, contrary to what one would hope for in SFDM, many galaxies do not end up in the MOND limit |ε_*|≪1. If the phonon force a_θ is close to its MOND-limit value , SFDM does give smaller averaged Q_* than MOND. Second, even in the MOND limit, the best-fit Q_* values in SFDM are systematically smaller than in MOND only for certain galaxy types. Such trends are not expected from SPS models. It also means we can never say that SFDM universally requires a smaller or larger Q_* than MOND. We can make such statements only for a given galaxy sample. We now discuss these two points in more detail.

D.2.2. Effect of going outside the MOND limit on M/L_*

As discussed in Appendix B.2, going to ε_* ≫ 1 allows us to make a_θ smaller so that a larger M/L_* is possible. This could be one reason why the averaged M/L_* is relatively large in SFDM.

As a first check, we show a scatter plot of Q_* versus ε_*(R_mid) (see Fig. 3). Indeed, many galaxies have ε_* ≫ 1. In addition, Fig. 3 shows a correlation between ε_* and Q_*. Galaxies with ε_* ≫ 1 tend to have a larger Q_*. This fits with the idea that we do not find a smaller M/L_* for SFDM because many galaxies are not in the MOND limit.

To confirm this, we redid the SFDM fit, but with the phonon acceleration a_θ replaced by its MOND limit value when calculating rotation curves. The calculation of was left untouched. This is the model shown as “SFDM ” in, for example, Table D.1 and Fig. D.2. With this model, the trick of going to ε_* ≫ 1 to make the phonon acceleration a_θ small does not work. Indeed, the averaged M/L_* is now significantly smaller than for the original SFDM fit. This result is again robust against different choices for the galaxies we consider and different choices for the averaging function. We can also see explicitly in Fig. D.4 that the distribution of best-fit ε values has migrated to smaller values compared to the original SFDM fit.

As a third check, we redid the SFDM fit, but with the model parameter β set to 1.55 instead of 2. This choice makes it much harder to make the phonon acceleration a_θ small by going to large ε_* ≫ 1, as can be seen from Fig. A.1. This is the model shown as “SFDM β = 1.55” in, for example, Table D.1 and Fig. D.2. If our explanation for the large M/L_* in SFDM is correct, this modified model should again have significantly smaller M/L_*. Indeed, this SFDM β = 1.55 model gives results that are comparable to those from the SFDM model. That is, M/L_* is now significantly smaller than for the SFDM fit. Similarly, the resulting ε values are much smaller than in the original SFDM fit (see Fig. D.4).

Thus, one reason for the relatively large M/L_* in SFDM is indeed that many galaxies are not actually in the MOND limit.

D.2.3. Enforcing the MOND limit

Its MOND limit is one of the main motivations of SFDM, because then rotation curves are automatically MOND-like. This is not the case outside the MOND limit (i.e., when ε_* is not small). Then MOND-like rotation curves are not possible without adjusting the boundary condition separately for each galaxy. Thus, our fit results for SFDM go against the original motivation behind SFDM.

We may wonder if large ε_* values are really necessary for SFDM to get reasonable fits of the SPARC data. It is possible that our fit code went to ε_* ≫ 1 for little gain in χ². To check this, we redid the SFDM fit, but with ε = ε_*(R_mid) restricted to |ε|< 0.1. Whenever we solved the SFDM equations and found |ε|≥0.1, we manually set χ² = 10⁷(1 + |ε|) so that the fitting code went elsewhere. The |ε| in 10⁷(1 + |ε|) is to help the χ²-minimizing fit algorithm to find small |ε|. In this fit, all galaxies are restricted to stay in the proper MOND limit |ε|≪1 of SFDM.

This works, but only for galaxies that are not too large. The restriction |ε|< 0.1 is impossible to satisfy for many larger galaxies. Our fitting procedure did not find any fit for 38 out of the 169 SPARC galaxies (i.e., 38 galaxies end up with χ² > 10⁷). This is not unexpected since our estimate ε_* ≳ 0.1 + 0.4 ⋅ (r/18 kpc − 1) from Eq. A.19 rules out |ε|< 0.1 for many larger galaxies. Indeed, Fig. D.5, top, shows that the galaxies that cannot be fit with |ε|< 0.1 tend to be those with R_mid ≳ 15 kpc.

Since we set χ² = 10⁷(1 + |ε|) when the condition |ε|< 0.1 was not satisfied, our fitting algorithm actually minimized |ε| until it satisfied |ε|< 0.1. Thus, we can get the minimum possible |ε| for each galaxy where |ε|< 0.1 could not be satisfied. The results are shown in Fig. D.5, bottom. We see that many galaxies only barely fail to satisfy the condition |ε|< 0.1. Indeed, if we allowed |ε| up to 0.4, almost all galaxies could be fit.

Of course, 0.4 is not that small, so it is debatable whether or not a value |ε| = 0.4 still counts as the proper MOND limit |ε_*|≪1. Here, we do not dwell on this point. However, we did redo our fit with the condition |ε|< 0.1 replaced by the condition |ε|< 0.4. As expected, we then obtained a fit for almost all galaxies. We did not find a fit for only four galaxies. Thus, if by “proper MOND limit” we mean |ε|< 0.1, SFDM’s proper MOND limit does not work for lager galaxies. But if we allow |ε| up to 0.4, it might.

The resulting best-fit χ² CDF is shown in Fig. D.1. In Fig. D.6, we show the changes in χ² between the “SFDM |ε|< 0.4” fit and the unrestricted SFDM fit. The resulting χ² values tend to be worse than for the unrestricted SFDM fit, but generally still acceptable. Indeed, they are quite similar to those of the MOND fit (see the CDF in Fig. D.1).

For some galaxies, the condition |ε|< 0.4 can barely be satisfied. After satisfying this condition they have basically no freedom left to actually fit the observed rotation curve data and they end up with very bad χ². Specifically, there are seven galaxies with ε > 0.38 and χ² > 100. Since these are hardly useful in assessing the M/L_* required in SFDM, we separately list the M/L_* results with these galaxies excluded in Tables D.2 and D.1. We also exclude them in Fig. D.6.

Somewhat surprisingly, some galaxies even have a better best-fit χ² with the |ε|< 0.4 restriction than without (see Fig. D.6). Some of these are just very slightly better than the previous best fit. For two galaxies it improves by more than 15%, specifically by 46% for NGC 1090 and by 19% for NGC 2683. For all galaxies with an improved χ², the corresponding best-fit M/L_* changes by less than 10%. These differences are insignificant for our purposes. They just show that our fitting algorithm is not perfect and does not always find the very best χ².

If we exclude the galaxies with a bad χ² because they can only barely satisfy |ε|< 0.4 as described above, the resulting averaged stellar mass-to-light ratios are between 4% and 11% smaller than for the MOND fit. The numerical details depend on the averaging procedure and the cut of galaxies. We discuss the resulting M/L_* values in more detail in Appendix D.2.4.

An alternative to the proper MOND limit |ε_*|≪1 is the pseudo-MOND limit discussed in Appendix B.2. At ε_* = 𝒪(1) (for β = 2), the phonon acceleration a_θ is still numerically close to its MOND limit value although it does not satisfy a MOND-like equation (see Fig. A.1, left). To test this regime, we redid the SFDM fit with ε restricted to |ε|< 5. This is the model shown as “SFDM |ε|< 5” in, for example, Table D.1 and Fig. D.2. Allowing values of |ε| up to 5 allows the phonon acceleration a_θ to deviate by up to about 5% from its MOND limit value (at R = R_mid) (see Fig. A.1). The resulting χ² values and stellar mass-to-light ratios are roughly comparable to those of the |ε|< 0.4 fits. As always, the numerical details depend on the choice of galaxies and on whether we average using the median or the mean.

Thus, fitting the SPARC data does not require ε_* ≫ 1 (i.e., it does not require going outside the MOND limit). Both the proper MOND limit and the pseudo-MOND limit also give reasonable χ² values. In this case, the averaged M/L_* is a bit smaller than in standard MOND models. In Appendix D.2.4, we discuss the M/L_* of these fits in more detail.

D.2.4. Trends of M/L_* with galaxy type

We now come back to the question of why SFDM does not necessarily need a smaller averaged M/L_* compared to MOND. Above, we already identified one reason, namely that many galaxies are not in the MOND limit |ε_*|≪1. But this is not the whole story, as we will now explain. To this end, we consider the fits with the MOND limit enforced as introduced in the previous subsection. This excludes effects from going outside the MOND limit.

In Appendix B.1, we argued that the MOND limit of SFDM likely requires a systematically smaller M/L_* than MOND only for high-acceleration galaxies. Galaxies with smaller accelerations may even require a larger M/L_* than MOND. The main reason is that SFDM has a smaller value of a₀ than MOND, which becomes important at small accelerations. If this is true, having a smaller or larger Q_* in SFDM is not just a property of the model but also a property of the galaxy sample. An example of this is Fig. 5, which shows the best-fit Q_* of each galaxy in the SFDM |ε|< 0.4 model relative to the best-fit Q_* for MOND as a function of the observed flat rotation velocity V_flat. A larger V_flat is associated with larger accelerations. Indeed, large V_flat values are where SFDM systematically gives smaller Q_* than MOND. Similarly, a smaller gas fraction and a larger surface brightness are associated with larger accelerations. The effective surface brightness Σ_eff and the ratio M_HI/L_[3.6], a proxy for gas fraction, of each galaxy are also part of SPARC. And indeed, Fig. D.7 shows that SFDM has a systematically smaller M/L_* than MOND for galaxies with a large Σ_eff and a small M_HI/L_[3.6].

To further test our understanding, we redid the MOND fit but using both the interpolation function ν_θ instead of ν_e (see Appendix B.1) and the smaller value of a₀. This should give fits qualitatively similar to those of the MOND limit of SFDM. Indeed, we verified that the resulting best-fit M/L_* show similar trends with, for example, V_flat as SFDM. As we explain in Appendix B.1, the different shape of the interpolation function is responsible for the systematically smaller M/L_* at large accelerations. The smaller a₀ value is responsible for the fact that this is not true at smaller accelerations. Thus, in a MOND model with the SFDM-like interpolation function ν_θ but with the larger a₀ value , we would expect to see a smaller M/L_* consistently across all accelerations. To test this, we redid the MOND fit with the interpolation function ν_θ but keeping the larger value of a₀. And indeed, this gives a consistently smaller M/L_* than in MOND. There is no clear trend with, for example, V_flat (see Fig. D.8).

Trends of the stellar M/L in the [3.6] micron band with galaxy properties are not expected from SPS models (Schombert et al. 2019). This disfavors SFDM, especially compared to MOND, which does not show such trends. If anything, one expects the opposite trend: more massive galaxies should have higher mass-to-light ratios than dwarfs, especially in optical bands.

More precisely, when we do not normalize Q_* to , our fitting procedure reproduces the SPS expectations for neither MOND nor SFDM. But this is a consequence of our simplistic fitting procedure. More sophisticated fits do reproduce the SPS expectations for MOND (McGaugh 2004). Moreover, since our simple fitting procedure is suited to identify relative differences between MOND and SFDM, we expect the trends in to be robust. It is unlikely that these would be mitigated by a more sophisticated fitting procedure. That is, we expect that this conflict between SFDM and SPS expectations is real.

Figures 5 and D.7 not only show that SFDM has a systematically small Q_* for large accelerations. They also show that there is more scatter in at smaller accelerations. One reason for this scatter is that many of the small-acceleration galaxies are gas-dominated. In gas-dominated galaxies, the formal best-fit value for the stellar mass-to-light ratio (i.e., Q_*) may not mean much, since a_b is relatively independent of Q_*. This allows for more scatter in Q_*.

D.2.5. The Milky Way

Hossenfelder & Mistele (2020) find that SFDM requires about 20% less baryonic mass than standard MOND models to fit the Milky Way rotation curve at R ≲ 25 kpc. Specifically, it requires about 20% less baryonic mass than the MOND model from McGaugh (2019). This is a significantly larger difference than what we find, on average, for the SPARC galaxies (see Tables D.1 and D.2). Here, we confirm and discuss this result.

To confirm the result of Hossenfelder & Mistele (2020), we fitted the Milky Way with the same method we used for the SPARC galaxies. For this, we used ρ_b and V_obs as in Hossenfelder & Mistele (2020). The V_obs data based on Portail et al. (2017) that is used in Hossenfelder & Mistele (2020) for R ≲ 2 kpc does not come with error bars. As a simple way to still get a result, we assumed an error of 5 km/s for these V_obs data points. For easier comparison to the SPARC fits, we rescaled the stellar disk and bulge densities such that stellar mass-to-light ratios of 0.5 (for the stellar disk) and 0.7 (for the bulge) correspond to the baryonic mass model used in McGaugh (2019). That is, the factor 10^f_Y tells us how much less stellar mass SFDM uses compared to the standard MOND model from McGaugh (2019). We found a best-fit χ² of 2.69, a best-fit Q_* of 0.79, and a best-fit ε of 4.37. This confirms the estimate from Hossenfelder & Mistele (2020) of about 20% less baryonic mass compared to standard MOND. This fit stays roughly within the pseudo-MOND limit ε_* = 𝒪(1). With the best-fit parameters, ε_* stays below 15 at R < 25 kpc. Thus, the phonon acceleration cannot be much suppressed and cannot allow for an increased M/L_*, in contrast to galaxies at ε_* ≫ 1 (see Appendix D.2.2).

The Milky Way’s a_b ranges from about 10⁻⁹ m/s² to about 10^−10.8 m/s² at R < 25 kpc for the best-fit parameters. These are relatively large. So, from the discussion in Appendix D.2.4, a smaller M/L_* than in MOND is what one would expect.

D.2.6. SFDM model parameters

Above, we used the fiducial parameter values from Berezhiani et al. (2018) and kept them fixed during the fitting procedure. Here, we discuss whether our conclusions could be changed by adjusting these parameters.

Our calculation does not depend on each of the four parameters, α, Λ, m, and β, separately. We need only the combinations

(D.1)

This can be seen directly from the SFDM Lagrangian (Berezhiani & Khoury 2015) by rescaling the phonon field, θ → θ(αΛ/M_Pl)⁻¹. It also follows from (Λm³)² = a₀(m²/α)³M_Pl.

First, the acceleration scale a₀. To reproduce MOND, a₀ must be close to 10⁻¹⁰ m/s². Still, we could choose the same value as in standard MOND rather than the somewhat smaller value that Berezhiani et al. (2018) chose. This would give M/L_* values that are smaller than those for MOND for all galaxies, not just the small-acceleration ones (see Appendix D.2.4), at least as long as the superfluid’s gravitational pull stays negligible. That is, SFDM would give M/L_* values similar to our “MOND ν_θ” fit (see Appendix D.2.4). These are relatively small. To get closer to M/L_* values as expected from SPS modeling, one would have to change not only the value of a₀ but also the form of the interpolation function ν_θ. This might be possible by adjusting the Lagrangian, that is, by changing what is usually called the function P(X) in superfluid low-energy effective field theories. Exploring this is beyond the scope of this paper.

One effect of the parameter β is that it controls the phonon force outside the proper MOND limit |ε_*|≪1 (see Appendix B.2). For one example (β = 1.55), we explicitly explored the effect of a different value of β (see Appendix D.2.2). Still, it is better not to tune this parameter for better fits to the data. This is because the value of β and the form of the finite-temperature corrections it is supposed to represent are completely ad hoc. So they might turn out to be unphysical. It is better to not rely too sensitively on any specific value of β for the fits.

The combination m²/α multiplies both ε_* and the superfluid’s energy density ρ_SF. One motivation to change m²/α is to make |ε_*| small in order to allow more galaxies to reach SFDM’s proper MOND limit |ε_*|≪1 (see Appendix A.1.2). The problem with this is that then ρ_SF ∝ (m²/α)⋅f_β(ε_*) becomes small as well. Indeed, in this case the problem regarding strong lensing described in Appendix D.4 becomes even worse. Conversely, making m²/α larger in order to solve the strong lensing tension means even fewer galaxies can reach the MOND limit |ε_*|≪1. For example, increasing m²/α by a factor of 10 would give ε_* ≳ 0.1 + 0.41(r/1.8 kpc − 1) from Eq. (A.19). Then, only the smallest galaxies could reach the proper MOND limit.

An adjusted function P(X) might invalidate this argument since this function determines not only the phonon force but also the superfluid’s energy density in SFDM. But, again, this is beyond the scope of the present work.

To sum up, adjusting the SFDM model parameters might change some of our conclusions regarding the best-fit M/L_* values, but probably not in a way that is completely satisfactory from the perspective of SPS models. Moreover, we expect that the tradeoff between having MOND-like rotation curves and producing sufficient strong lensing described in Appendix D.4 remains.

D.2.7. Thermal radius check

We assumed that all SPARC rotation curve data points of a given galaxy lie within this galaxy’s superfluid core. This is necessary for one of the main motivations behind SFDM, namely to automatically reproduce MOND-like rotation curves without having to adjust the boundary condition separately for each galaxy. The superfluid phase ends at the very latest when ρ_SF reaches zero. Thus, we discarded solutions where ρ_SF vanishes inside the V_obs data points, as discussed in Appendix C.2.2.

But this may not be sufficient since the superfluid phase may end even before ρ_SF reaches zero. Consider, for example, the simplest estimate for the radius where the superfluid phase transitions to the non-superfluid phase, the so-called thermal radius R_T (Berezhiani et al. 2018). According to this estimate, the superfluid phase corresponds to

(D.2)

where Γ is the local self-interaction rate and t_dyn is the dynamical time. Here, Γ = (σ/m) 𝒩 v ρ, where σ is the self-interaction rate, 𝒩 = (ρ/m)(2π/mv)³ is the Bose-degeneracy factor, and v is the average velocity of the particles. Following Berezhiani et al. (2018), we take σ/m = 0.01 cm²/g and t_dyn = R/v.

As a simple check, we evaluated the quantity for each galaxy at the last V_obs data point at R_max for the SFDM best fits. We found that 31 of 169 galaxies violate the condition at R_max. In principle, we should discard these solutions, just as we discard solutions where ρ_SF reaches zero before R_max. Here, we do not do this. The reason is that the condition is quite ad hoc. For example, the value of σ/m is chosen ad hoc and not derived from an underlying Lagrangian. Also, the transition radius derived from can, in general, differ wildly from the transition radius derived from the so-called NFW matching procedure (Berezhiani et al. 2018; Hossenfelder & Mistele 2020) where the density and pressure are matched to those of an NFW halo for a fixed NFW concentration parameter, as discussed in Mistele (2021). This makes any particular choice for discarding solutions based on the thermal or NFW matching somewhat arbitrary.

We avoid this arbitrariness by discarding solutions only based on the criterion that ρ_SF must be positive. Still, this means we do not discard some solutions that we maybe should discard. This might affect our M/L_* results. To get a rough idea of possible effects due to this, we also show our results for M/L_* in Table D.1 and Table D.2 with the 31 galaxies violating excluded, labeled as “thermal ok.” Often, this does not significantly change the M/L_* fit results, though it does give a larger M/L_* for example for the SFDM model.

Overall, we expect that actually enforcing the rotation curve data to lie within the superfluid phase does not significantly change our M/L_* results. But keep in mind that there is a considerable theoretical uncertainty around the transition from the superfluid to the non-superfluid phase. We stress that the tension with strong lensing (see Appendix D.4) is not affected by this uncertainty, because there we anyway choose the transition radius to maximize the resulting total dark matter mass.

D.3. The RAR implied by the SFDM best fits

We show the RAR (i.e., the relation between a_b and a_tot) implied by our SFDM fits for the q = 1 galaxies in Fig. D.9. For each galaxy, we show one data point for each radius with a SPARC rotation curve data point V_obs. As expected, the fit restricted to the proper MOND limit (|ε|< 0.4) shows a tight relation with almost no scatter. The unrestricted fit deviates from this in two ways. First, it produces significantly more scatter. Second, it systematically puts more data points below the relation implied by the proper MOND limit fit (i.e., it puts more data points at larger a_b, smaller a_tot, or both). Of course, outside the MOND limit of SFDM, there is no guarantee that one ends up with a MOND-like RAR. Still, we can understand in more detail why our unrestricted fits look the way they do.

One contribution to the increased scatter is as follows. The observed rotation curves contain some scatter from observational uncertainties. Our unrestricted SFDM fits match these observed rotation curves more closely than the fits restricted to the MOND limit (see, e.g., Fig. D.6). Thus, overfitting the noise in the observational data may, at least partly, explain the increased scatter in the RAR implied by our unrestricted fits.

We can also identify a possible cause for the systematic deviation from a MOND-like RAR. Namely, some observed rotation curves can be matched more closely outside the MOND limit (i.e., with a large |ε_*|). As discussed in Appendices B.2 and D.2.2, going outside the MOND limit allows for larger M/L_*: a large |ε_*| makes a_tot smaller, which must be countered by a larger M/L_*. These larger M/L_* values explain why the RAR implied by our unrestricted SFDM fit puts significantly more data points below the tight RAR from the MOND limit than above it. Since this mechanism affects only some galaxies and since it affects different galaxies differently, this likely also contributes to the increased scatter.

When comparing the RAR implied by our fits to the actually observed RAR, one should keep in mind that, as mentioned above, observations always add scatter on top of what an underlying theory predicts. So, for example, if our unrestricted SFDM best fits are the ground truth, the observed RAR should contain even more scatter than what Fig. D.9, left, shows. At least if one ignores that some of the scatter in Fig. D.9, left, already reflects the scatter in the observed rotation curve data due to overfitting. The result may be even more scatter than what one obtains in actual observations with the simple prescription (M/L_*)_disk = 0.5 and (M/L_*)_bulge = 0.7 (Lelli et al. 2017b). Since SFDM aims to explain this observational fact, one might reject our unrestricted fit results purely on theses grounds.

Still, the objective of our fits was to match the observed rotation curve data, not to get a tight RAR. So our fit results do not directly imply that a non-tight RAR is an intrinsic property of going outside the MOND limit of SFDM. It may still be possible to get a tight RAR even outside the MOND limit. However, as discussed in Sect. 2, for this one would have to carefully adjust the boundary condition ε for each galaxy. Then, one must rely on galaxy formation to always pick the right values. This is contrary to the aim of SFDM to explain scaling relations such as the RAR without having to resort to the details of galaxy formation.

D.4. Tension with strong lensing in SFDM

The proper MOND limit |ε_*|≪1 of SFDM is useful for fitting rotation curves due to its MOND-like phonon force. Above, we saw that most SPARC rotation curves can be reasonably fit with the proper MOND limit |ε_*|≪1. At least if we count values of |ε_*| as large as 0.4 as still satisfying |ε_*|≪1. But for certain other observables, such as strong lensing, the phonon force plays no role. This is because the observation of GW170817 requires that the phonon force does not affect photons (Hossenfelder & Mistele 2019; Sanders 2018; Boran et al. 2018). Thus, the strong lensing signal is produced only by the standard gravitational pull of the mass of the baryons and the superfluid, not by the phonon force.

As a consequence, the MOND limit of SFDM has a serious problem, because the superfluid’s energy density is too small to produce significant lensing. Specifically, assuming the whole superfluid core to be in the MOND limit |ε_*|≪1, a rough upper bound is (see Appendix A.1.3)

(D.3)

where we assumed the numerical parameter values from Berezhiani et al. (2018). Producing sufficient strong lensing and a superfluid core in the MOND limit seem to be mutually exclusive. Choosing different parameter values may help, especially a larger m²/α. But this would imply that fewer galaxies can reach the MOND limit |ε_*|≪1, as discussed in Appendix D.2.6.

This conclusion can, in principle, be avoided if the superfluid core is in the MOND limit only at smaller radii (where rotation curves are measured), but not at larger radii (where part of the lensing signal comes from). But even in this case there are limits on how large can be. This is because, given that we have |ε_*|≪1 at relatively small radii, the superfluid’s energy density cannot be arbitrarily large at larger radii. Here, we check whether or not the SPARC galaxies can possibly have a sufficiently large M₂₀₀/M_b for strong lensing, if we assume the proper MOND limit for the rotation curves.

Specifically, for each galaxy, we will find the largest possible value of that is compatible with a rotation curve in the proper MOND limit. For the requirement that the rotation curve is in the proper MOND limit we impose |ε|< 0.4 as above.

In SFDM one usually assumes that the superfluid ends at a radius r_NFW beyond which the energy density is that of an NFW halo (Berezhiani & Khoury 2015; Berezhiani et al. 2018; Hossenfelder & Mistele 2020). For simplicity, we assumed the NFW profile ρ_NFW(r) to be proportional to 1/r³ at the radii of interest, that is to say, we continued the superfluid density with an NFW tail instead of a full NFW profile. We do not expect this to significantly affect our results (Hossenfelder & Mistele 2019). The usual procedure for matching the superfluid density to the NFW profile is heuristic and not derived from first principles (Berezhiani et al. 2018; Mistele 2021). To be independent of the details of this matching procedure we matched the NFW density ρ_NFW to the superfluid density ρ_SF at a radius r_NFW that is chosen to maximize . This gives the most conservative upper bound for . We restricted r_NFW only in two ways. First, we assumed a positive superfluid energy density within the superfluid core, which implies that r_NFW is smaller than the radius r_m where ρ_SF vanishes. Second, we assumed all rotation curve data points to be within the superfluid core. Thus, we also restricted r_NFW to be larger than R_max. This is illustrated in Fig. D.10.

To get the largest possible compatible with our constraints, we scanned values of ε in (ε_*min, 0.4) and values of r_NFW in (R_max, r_m) and recorded the largest as an upper bound . For a given ε and r_NFW we calculated by solving the equation

(D.4)

for r₂₀₀ and then plugging the result into . We calculated M_SF(r) as . For some galaxies, it is not possible to solve the equation with |ε|< 0.4, if we require a positive ρ_SF up to the last rotation curve data point. Such galaxies are excluded in the results shown below.

For the baryonic mass distribution we used the same values (M/L_*)_disk = 0.5 and (M/L_*)_bulge = 0.7 for all galaxies. These may not give the best fit to the measured rotation curves for all galaxies. However, here we are not interested in the fit to the rotation curve data, but in whether or not it is possible to have, at the same time, both a rotation curve in the proper MOND limit and sufficient for strong lensing. The precise value of M/L_* is irrelevant for this.

We can simplify the scanning procedure a bit. Namely, Fig. D.13 and Fig. D.14 suggest that increasing the boundary condition ε ≡ ε_*(R_mid) increases ε_*(r) at all radii, not just at R_mid. We verified this numerically for various galaxies and boundary conditions. Here, we assumed that this is true in general.² Consequences of this are that a larger boundary condition ε implies an everywhere larger superfluid energy density and a larger radius r_m where this energy density reaches zero. Therefore, for a fixed r_NFW, the quantity is a monotonically increasing function of ε. Thus, we can simplify our scanning procedure by always setting ε = 0.4 and scanning only values of r_NFW in (R_max, r_m). For this, we used Mathematica’s NMaximize function with its default options.

It can happen that the NFW radius r_NFW is bigger than the radius r_∞ up to which we numerically solved for (usually 100 kpc; see Appendix C.1). In this case, we must continue beyond r_∞ since solving the equation requires and a_b up to r_NFW. For simplicity, we continued beyond r_∞ assuming spherical symmetry and zero baryonic energy density at r > r_∞,

(D.5)

This implies . We expect errors due to this to not significantly change our results.

Below we need the total baryonic mass M_b of each galaxy. We adopt

(D.6)

where L_[3.6], L_bulge, and M_HI are taken directly from SPARC. They denote the total [3.6] luminosity, the total bulge luminosity, and the HI mass, respectively. The stellar luminosities are weighted by their respective mass-to-light ratio. The factor 1.4 in front of the HI mass is to take into account helium and molecular hydrogen (McGaugh et al. 2020).

We show the results in Fig. 6. We see that large ratios are easier to reach for galaxies with relatively small baryonic masses M_b. In part, this is due to the factor 1/M_b in . For strong lensing, relatively large baryonic masses are relevant. To illustrate this, we also show the best-fit values from the SFDM strong lensing analysis from Hossenfelder & Mistele (2019) in Fig. 6. These lensing galaxies tend to have M_b ≳ 10¹¹ M_⊙ and . In contrast, the SPARC galaxies with M_b > 10¹¹ M_⊙ almost all have when restricted to have rotation curves in the MOND limit (i.e., when restricted to |ε|< 0.4). This is a stark contrast. The SPARC galaxies do not reach baryonic masses quite as large as the lensing galaxies from Hossenfelder & Mistele (2019). But from Fig. 6, it seems clear that the trend goes into the wrong direction: The larger the galaxy, the smaller the ratio (assuming |ε|< 0.4).

Thus, it seems that a rotation curve in the MOND limit and sufficient dark matter for strong lensing are indeed mutually exclusive in standard SFDM. A caveat is that the galaxies in the SPARC sample are not ellipticals, in contrast to the lensing sample used in Hossenfelder & Mistele (2019). One might think that, for a given M_b, the maximum possible is not sensitive to the details of the baryonic mass distribution since the main contributions to come from large radii where, to a first approximation, only the total M_b plays a role. However, we impose the condition |ε|< 0.4 at relatively small radii R = R_mid where the details of the baryonic mass distribution may still matter. Indeed, the sensitivity of to these details is reflected in the scatter in Fig. 6 (see also the end of Appendix D.6). Still, even being generous with this scatter, it seems unlikely that the difference in galaxy type explains why the SPARC galaxies restricted to |ε|< 0.4 cannot reach larger values of . Also, Fig. 6 suggests that for values closer to those required for strong lensing (when restricting to |ε|< 5 instead of |ε|< 0.4) there is less scatter (i.e., less sensitivity to the baryonic mass distribution beyond the total M_b).

This suggests that strong lensing galaxies cannot have their inner parts (where rotation curves or velocity dispersions are measured) be in the proper MOND limit |ε_*|≪1. This is not in direct contradiction with measurements. Indeed, Hossenfelder & Mistele (2019) successfully fitted strong lensing data in SFDM. But one cannot easily keep the key idea of SFDM that the inner parts of galaxies are always in the proper MOND limit. Either one has to give up this key idea or one has to postulate that it does not apply to strong lensing galaxies for some reason.

D.5. M/L_* in two-field SFDM

By construction, the phonon acceleration in two-field SFDM should almost always be close to (Mistele 2021). The superfluid’s Newtonian gravitational pull can be comparable to that of standard SFDM (Mistele 2021). Thus, we expect the fit results for two-field SFDM to be close to that of the SFDM model discussed above.

This should be true at least for the stellar mass-to-light ratios and the best-fit χ². The results for f_MDM may be not as close since the superfluid halo’s shape is different. We discuss the superfluid halo in more detail below.

As expected, the best fits for two-field SFDM are almost all in the |ε_*|≪1 limit so that their phonon force a_θ is close to . This is shown in Fig. D.4. Only for two galaxies (NGC 6789, UGC 07232) does ε_* become larger than 0.1. Its largest values is 0.36 for NGC 6789.

The averaged best-fit stellar mass-to-light ratios and the best-fit χ² can be found in Table D.1, Table D.2, and Fig. D.12. As expected, these are almost identical to those of the SFDM model discussed above.

Two-field SFDM gives a systematically smaller M/L_* than MOND only for high-acceleration galaxies, not for low accelerations, just as standard SFDM (see Sect. 5.1 and Appendix D.2.4, and see Fig. D.11 for an example).

D.5.1. The minimum acceleration

In two-field SFDM, the equilibrium on galactic scales is stable only for phonon accelerations above a certain minimum acceleration. This minimum acceleration depends on the value of the field . Assuming |ε_*|≪1, as is usually the case, stability requires

(D.7)

where is one of this model’s parameters. Often, is on the order of 10⁻⁷ on galactic scales. Thus, roughly, stability requires

(D.8)

When is smaller, the instability sets in earlier. Mistele (2021) chose so that this model does not predict standard MOND-like behavior for dwarf spheroidals that may start to deviate from MOND-like behavior around a_b ∼ 10⁻¹² m/s² (Lelli et al. 2016). In our two-field SFDM fit we adopted this value of .

For galaxies violating Eq. (D.7), we should in principle model what happens beyond the minimum acceleration in two-field SFDM. Here, we did not do this for two reasons. First, this regime is not well-understood. Second, we are interested in the MOND regime inside the superfluid core. Modeling the behavior beyond the minimum acceleration will not help us understand whether or not the MOND regime of two-field SFDM requires a larger or smaller M/L_* than standard MOND models.

Still, 98 of 169 SPARC galaxies violate the condition Eq. (D.7) at R_max. For these galaxies, our fit is not meaningful since it relies on an unstable equilibrium. This is expected for the dwarf spheroidals with a_b around 10⁻¹² m/s². But the 98 galaxies violating Eq. (D.7) include many more galaxies, also at . We will have to deal with this one way or another.

To further explore this, we redid the two-field SFDM fit with the much smaller value . This is listed as “two-field ” in our tables and figures. The resulting M/L_* and χ² values are almost identical to those of the previous two-field SFDM fit. But still 71 galaxies violate Eq. (D.7).

Thus, Eq. (D.7) is often not violated because a_b is smaller than but because is smaller than 10⁻⁷. A small corresponds to a small superfluid mass M_DM. Indeed, Fig. D.3 shows that many galaxies have a smaller f_MDM in two-field SFDM compared to standard SFDM, for both values discussed above.

This raises two questions. The first is whether two-field SFDM really needs to go to small to fit the SPARC data, thus often violating Eq. (D.7). The second is why standard SFDM tends to end up at larger f_MDM values than two-field SFDM.

D.5.2. Origin of the small f_MDM values in two-field SFDM

The best-fit χ² and M/L_* are almost identical for the two-field SFDM fits and for the standard SFDM fit with . But the superfluid halos of two-field SFDM reach much smaller masses than those in standard SFDM (see Fig. D.3). For example, there are no galaxies at f_MDM < −2.3 in any standard SFDM fit, but many such galaxies in two-field SFDM. The only relevant difference between two-field SFDM and the SFDM model is the difference in their ρ_SF. Thus, this different superfluid energy density must be the reason for the qualitative difference in f_MDM. Here, we explain this in more detail.

As discussed above, we enforced a positive superfluid energy density at radii smaller than the last rotation curve data point (i.e., ρ_SF > 0 at R ≤ R_max). The difference between standard and two-field SFDM regarding small dark matter masses boils down to what this condition implies.

Consider first two-field SFDM. In two-field SFDM, the superfluid energy density vanishes when . Equivalently, when ε_* = 0. Typically, is a decreasing function of galactocentric radius.³ Thus, whenever the condition ρ_SF > 0 is fulfilled at the largest radius of interest, R = R_max, it is fulfilled at all radii of interest (i.e., at R ≤ R_max). In a given galaxy, a smaller superfluid mass M_SF(R) at some radius, for example at R = R_max, implies a smaller everywhere in the superfluid core (not just at R = R_max). Indeed, different superfluid masses just correspond to adding a term sin(r/r₀)/r with a different prefactor to ; see Appendix C.4. Thus, if we go to smaller and smaller f_MDM, we go to smaller and smaller . At some point, will become negative at R = R_max. Then, we are at the minimum possible f_MDM allowed by our condition ρ_SF > 0.

The condition ρ_SF > 0 at R = R_max enforces a minimum possible mass also in standard SFDM. A difference is that the superfluid energy density in standard SFDM vanishes not at ε_* = 0 but at a negative value ε_*min ≈ −0.31 (for β = 2). Using the definition of ε_* from Eq. (1), this constant negative lower bound on ε_∗(x) becomes a nonconstant lower bound on . This is illustrated in Fig. D.13 for NGC 2403. The top panel shows solutions for various boundary conditions ε. The bottom panel shows ε_* for the same boundary conditions. Both panels also show the lower bounds on ε_* and , respectively, which ensure ρ_SF > 0. For the smallest boundary condition value ε = 0.08 shown in Fig. D.13, both ε_* and reach this lower bound at R = R_max. Thus, ε = 0.08 corresponds to the smallest possible superfluid mass that is allowed by the condition ρ_SF > 0. Lower superfluid masses would need ρ_SF to reach zero before R = R_max, which we do not allow. This is similar to as in two-field SFDM.

However, in contrast to two-field SFDM, this constraint at R = R_max is not the only constraint on superfluid masses in standard SFDM. In standard SFDM, having ρ_SF > 0 at R = R_max does not imply ρ_SF > 0 at smaller radii. The main issue is at R = 0. This is illustrated in Fig. D.14 for DDO 168. For the smallest boundary condition value ε = −0.06 shown in Fig. D.14, at R = R_max, both and ε_* are not close to the minimum values that ρ_SF > 0 allows. However, at R = 0, the minimum values are reached (i.e., ρ_SF vanishes). Thus, the boundary condition ε = −0.06 corresponds to the minimum possible f_MDM allowed by our condition ρ_SF > 0. But now the constraint comes from R = 0 rather than R = R_max.

To understand this, consider small radii, R → 0. At R → 0, the field reaches a finite value . The Newtonian baryonic acceleration |a_b(R)| goes to zero, typically a_b ∝ R for R → 0. Thus, we have from the definition of ε_* Eq. (1)

(D.9)

That is, ε_* tends to ±∞ with the sign being that of . A positive ρ_SF requires ε_*(0)≥ε_*min. Thus, we must have

(D.10)

This condition must also hold in two-field SFDM. But, in two-field SFDM, follows from ρ_SF > 0 at R = R_max. This is because is a decreasing function of radius and because ρ_SF > 0 requires even at R = R_max. So this adds nothing in two-field SFDM.

This is different in standard SFDM. A positive ρ_SF at R = R_max does not necessarily imply a positive and thus a positive ρ_SF at R = 0. Since is a decreasing function of radius, violating is possible only if is negative already at R = R_max. Consider this case where is negative at R = R_max. Since decreases with radius, it can happen that grows sufficiently between R = R_max and R = 0 to become positive at R = 0. In this case the condition , corresponding to ρ_SF > 0 at R = 0, gives no additional constraint. However, this is not guaranteed to happen. When does not grow sufficiently, the condition gives an additional constraint. This is what happens for DDO 168 as illustrated in Fig. D.14.

This, then, is the reason why two-field SFDM allows smaller superfluid masses (i.e. smaller f_MDM) compared to standard SFDM. In both models, small superfluid masses correspond to ρ_SF close to zero. And in both cases one needs to be careful not to let this density become negative (or ill-defined) at R = R_max. However, in standard SFDM, being close to ρ_SF = 0 implies a negative , and in this case there can be an additional constraint at R = 0, as just discussed. This second constraint is absent in two-field SFDM where ρ_SF vanishes already at .

As mentioned above, this second constraint in standard SFDM occurs only when does not grow sufficiently between R = R_max and R = 0. We now make this more precise. The minimum allowed superfluid mass from the constraint at R = R_max corresponds to ρ_SF = 0 at R = R_max. This constraint is present in both standard and two-field SFDM. Reaching this minimum mass implies

(D.11)

A second constraint from R = 0 is avoided whenever grows between R = R_max and R = 0 by at least

(D.12)

As the example of DDO 168 (see Fig. D.14) shows, this is not always guaranteed in standard SFDM. But in some cases it is. Namely, as in Appendix C.2, we can write . Both and typically decrease with radius, with the amount they decrease being determined by the baryonic and superfluid mass, respectively. Once we know the baryonic mass distribution, we know a lower bound on how much the total grows, independently of the superfluid energy density and the boundary condition ε. As a result, if the baryonic mass alone makes grow sufficiently, we never get an additional constraint from . This is the case if

(D.13)

Thus, more quantitatively, we claim that standard SFDM cannot get masses as low as two-field SFDM because Eq. (D.13) is often not satisfied. That is, there is an additional constraint from R = 0 in standard SFDM that is not present in two-field SFDM. We expect that, without this constraint, galaxies would end up at smaller f_MDM also in standard SFDM.

As a consequence, we expect that the galaxies that end up at very small f_MDM in two-field SFDM but not in standard SFDM violate Eq. (D.13). To check this, we selected the 30 galaxies that have f_MDM < −2.3 in two-field SFDM. As mentioned above, no galaxies have such small f_MDM in standard SFDM. For the most direct comparison against two-field SFDM we use the best-fit stellar M/L_* from the SFDM fits. We find that these 30 galaxies all violate Eq. (D.13) so that they face an additional constraint in standard SFDM. This confirms our explanation why galaxies reach smaller f_MDM values in two-field SFDM compared to standard SFDM.

D.5.3. Enforcing the minimum acceleration in two-field SFDM

In Appendix D.5.1, we saw that many best fits in two-field SFDM violate the minimum acceleration condition Eq. (D.7), even with a reduced value. The reason is that many galaxies end up at small , or, equivalently, at small superfluid masses corresponding to small f_MDM (see Appendix D.5.2).

For small superfluid masses, the precise mass is often not important for fitting rotation curves since the corresponding a_SF is subdominant. Thus, it may be possible that we can find fits with sufficiently large so that all galaxies satisfy the minimum acceleration condition without getting significantly worse fits. To check this, we redid the two-field SFDM fit with but with the minimum acceleration condition Eq. (D.7) enforced. That is, whenever this condition is violated we set χ² = 10¹⁰ so that our fit code goes elsewhere. We label this model “two-field .”

As we can see from Table D.1, Table D.2, and Fig. D.12, this gives almost identical results for the best-fit χ² and M/L_* values as the previous two-field SFDM fits. This shows that it is not necessary for two-field SFDM to violate the minimum acceleration condition. The small values are not required to get a reasonable fit (see also Fig. D.3).

D.6. Tension with strong lensing in two-field SFDM

In two-field SFDM, the |ε_*|≪1 condition is almost always fulfilled so that the phonon force is almost always close to the MOND-like value . So, in contrast to standard SFDM, simultaneously being in the |ε_*|≪1 limit and producing a sufficient strong lensing signal is not a problem. However, in two-field SFDM, a small ε_* does not imply that a_SF is negligible. So two-field SFDM may still not be able to produce MOND-like rotation curves and sufficient strong lensing at the same time.

To check this, we calculated a maximum total dark matter mass in the same way as we did for standard SFDM in Appendix D.4. But instead of imposing |ε|< 0.4 we imposed

(D.14)

at R = R_max. For simplicity, we assumed , which is usually a good approximation. Then, Eq. (D.14) becomes

(D.15)

To implement this numerically, we used M_DM(R_max) instead of ε_*(R_mid) as a boundary condition. This translates to a boundary condition on instead of on . We can still solve the two-field equations as described in Appendix C.4 with this modified boundary condition.

The value 0.3 is somewhat arbitrary. It is chosen to give a non-negligible but still subdominant a_SF. There is a clear trade-off here. Larger values allow for larger total dark matter masses but also larger deviations from MOND-like rotation curves. Smaller values give more MOND-like rotation curves but also smaller total dark matter masses.

One way to find the largest possible compatible with the condition Eq. (D.15) is to just scan over all possible M_DM(R_max) satisfying the condition. But we do not have to do this here for the same reason we did not have to do it for standard SFDM (see Appendix D.4). Basically, since a larger M_DM(R_max) always gives a larger .

Above, we considered different values of the parameter , which determines the minimum acceleration. Here, we use . But the choice of does not actually affect the maximum we calculate since we keep r₀ and a₀ fixed. It only determines how small ε_* is since . That is, it determines how well our approximation works. Since, by construction, |ε_*|≪1 for any reasonable value of , this approximation works well for any reasonable .

The result is shown in Fig. 9 and discussed in Sect. 5.3. Comparing Fig. 9 for two-field SFDM and Fig. 6 for standard SFDM, there is less scatter in the derived relation between and M_b for two-field SFDM. This is likely because this relation in two-field SFDM is less sensitive to the details of the baryonic mass distribution at a given total baryonic mass M_b. There are a few reasons for this. First, the superfluid energy density ρ_SF scales as in standard SFDM (at least in the MOND limit |ε_*|≪1) but not in two-field SFDM. Second, we impose the condition |ε|< 0.4 in standard SFDM at smaller radii R = R_mid than a_SF/(a_b + a_θ) < 0.3 in two-field SFDM, which we impose at R = R_max. There is less variation in the baryonic mass distribution at larger radii. Third, these conditions depend on different quantities. The left-hand side of |ε|< 0.4 in standard SFDM scales as while the left-hand side of a_SF/(a_b + a_θ) < 0.3 in two-field SFDM scales as , at least at larger radii where dominates. The acceleration a_SF depends only on , while also depends on , which is much more sensitive to the details of the baryonic mass distribution. Similarly, is less sensitive to these details than a_b.

All Tables

All Figures