Confronting fuzzy dark matter with the rotation curves of nearby dwarf irregular galaxies

Andrés Bañares-Hernández; Andrés Castillo; Jorge Martin Camalich; Giuliano Iorio

doi:10.1051/0004-6361/202346686

Home

All issues

Volume 676 (August 2023)

A&A, 676 (2023) A63

Full HTML

Open Access

This article has an erratum: [https://doi.org/10.1051/0004-6361/202449974e]

Issue		A&A Volume 676, August 2023


Article Number		A63
Number of page(s)		21
Section		Cosmology (including clusters of galaxies)
DOI		https://doi.org/10.1051/0004-6361/202346686
Published online		07 August 2023

A&A 676, A63 (2023)

Confronting fuzzy dark matter with the rotation curves of nearby dwarf irregular galaxies

Andrés Bañares-Hernández¹^,2, Andrés Castillo¹^,2, Jorge Martin Camalich¹^,2 and Giuliano Iorio³^,4^,5

¹ Instituto de Astrofísica de Canarias, La Laguna, Tenerife 38200, Spain
e-mail: abh@iac.es
² Departamento de Astrofísica, Universidad de La Laguna, La Laguna, Tenerife 38200, Spain
³ Dipartimento di Fisica e Astronomia Galileo Galilei, Università di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy
⁴ INFN-Padova, Via Marzolo 8, 35131 Padova, Italy
⁵ INAF-Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy

Received: 18 April 2023
Accepted: 18 May 2023

Abstract

Aims. In this paper, we carry out a phenomenological investigation of the viability of fuzzy dark matter, which is composed of coherent waves of non-interacting ultralight axions with a mass of m_a ≈ 10⁻²² eV. We did so by confronting the predictions of the model, in particular, the formation of a solitonic core at the center of dark matter halos, with a homogeneous and robust sample of high-resolution rotation curves from the LITTLE THINGS in 3D catalog. This comprises a collection of isolated, dark matter-dominated dwarf-irregular galaxies that provides an optimal benchmark for cosmological studies. Our aim is to find evidence of fuzzy dark matter in the observations; alternatively, we seek to set exclusion ranges for its mass.

Methods. We used a statistical framework based on a χ² analysis of the rotation curves of the LITTLE THINGS in 3D catalog using a fuzzy dark matter profile as the theoretical model. This allows us to extract relevant parameters such as the axion mass and mass of the solitonic core, as well as the mass of the dark matter halo and its concentration parameter. We fit the data using current Markov chain Monte Carlo techniques with a rather loose set of priors, except for the implementation of a core-halo relation predicted by simulations. The results of the fits were then used to perform various diagnostics on the predictions of the model.

Results. Fuzzy dark matter provides an excellent fit to the rotation curves of the LITTLE THINGS in 3D catalog, with axion masses determined from different galaxies clustering around m_a ≈ 2 × 10⁻²³ eV. However, we find two major problems from our analysis. First, the data follow scaling relations of the properties of the core, which are not consistent with the predictions of the soliton. This problem is particularly acute in the core radius-mass relation with a tension that (at face value) has a significance of ≳5σ. The second problem is related to the strong suppression of the linear power spectrum that is predicted by fuzzy dark matter for the axion mass preferred by the data. This can be constrained very conservatively by the galaxy counts in our sample, which leads to a tension that exceeds 5σ. We estimate the effects of baryons in our analysis and discuss whether they could alleviate the tensions of the model with observations.

Key words: dark matter / galaxies: dwarf / galaxies: kinematics and dynamics / astroparticle physics

Note to the reader: In the Acknowledgements, the name of J. Garcia-Farieta has been misspelled. It was corrected on 10 August 2023.

© The Authors 2023

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

One of the most pressing open problems in contemporary physics is the fundamental nature of dark matter (DM). The prevailing hypothesis is that DM is formed by a new type of particle, or a “dark sector” that is not contained within the standard model of particle physics and behaves as a cold and collisionless fluid on large cosmological scales (Silk et al. 2010; Bertone & Hooper 2018; Chou et al. 2022). This cold DM (CDM) is a crucial ingredient of the standard cosmological model, known as ΛCDM, which has proven to be very successful in describing observations of the universe at these large scales (Planck Collaboration VI 2020; Alam et al. 2017; Workman et al. 2022; Peebles 2022). However, observations probing the formation of structure at small galactic scales might pose a challenge to ΛCDM.

Stellar and gas kinematics indeed provide one of the most direct observational probes of the internal structure of DM halos (Bosma 1978; Rubin et al. 1980; Trimble 1987). While pure N-body CDM simulations generically predict density profiles peaking as ρ(r)∝1/r at small radii (Dubinski & Carlberg 1991; Navarro et al. 1996b, 1997), the rotation curves of DM-dominated dwarf galaxies point instead to the existence of flat cores at their centers (Flores & Primack 1994; Moore 1994; de Blok & Bosma 2002; Oh et al. 2011, 2015; Read et al. 2017; Li et al. 2020; Dehghani et al. 2020). This is known as the “core-cusp” puzzle that is among the small-scale problems related to ΛCDM (see Bullock & Boylan-Kolchin 2017 and Sales et al. 2022 for recent reviews). Other possible discrepancies between the predictions of the model and observations are referred to as the “missing satellites” problem that is related to the observation of many fewer low-mass galaxies in the Local Group than otherwise expected from CDM simulations (Klypin et al. 1999; Moore et al. 1999; McConnachie 2012). Then there is the “too-big-to-fail” problem, whereby the central masses are too low compared to the most massive (sub)halos predicted in ΛCDM (Boylan-Kolchin et al. 2011; Kirby et al. 2014).

These problems have prompted a flurry of activity over the past decade theorizing and analyzing phenomenologically models beyond CDM. One prominent candidate is fuzzy DM (FDM), which is composed of ultralight axions with masses in the range of ∼10⁻²³ − 10⁻²¹ eV (Sin 1994; Hu et al. 2000; Matos & Urena-Lopez 2001; Chavanis 2011; Schive et al. 2014a,b). This type of particle has been predicted by several theoretical models, preserving the features of CDM at large cosmological scales and which can be produced “naturally” in the early universe to match the DM relic abundance (Marsh 2016; Hui et al. 2017; Niemeyer 2020; Ferreira 2021; Hui 2021). On the other hand, on distances comparable to their de Broglie wavelength, which can be in the kiloparsec (kpc) scale, these particles populate the galactic halos with large occupation numbers and behave as self-gravitating DM waves. One of the consequences of this is the appearance of a pressure-like effect on macroscopic scales leading to the formation of a flat core, or “soliton”, at the center of galaxies with a relatively marked transition to a less dense outer region that follows a CDM-like distribution (Schive et al. 2014a,b; Veltmaat et al. 2018; Mocz et al. 2018; Nori & Baldi 2021; Mina et al. 2022; Chan et al. 2022). This effect also produces a suppression of the linear power spectrum on small scales and to a decrease of the number of low-mass halos (Khlopov et al. 1985; Schive et al. 2016; Du et al. 2017; May & Springel 2021, 2023). The FDM provides, then, an elegant solution to the small-scale issues of ΛCDM, which has motivated several phenomenological studies and searches in astrophysical and cosmological probes of low-scale structure (see Ferreira 2021 for a review and Sect. 5 below for a discussion).

Other attempts to solve some of the small-scale issues of ΛCDM by modifying the fundamental properties of DM are warm DM (Sommer-Larsen & Dolgov 2001; Bode et al. 2001; Rubakov & Gorbunov 2017) including degenerate fermion DM (Tremaine & Gunn 1979; Gorbunov et al. 2008; Domcke & Urbano 2015; Bar et al. 2021; Chavanis 2022; Krut et al. 2023), self-interacting FDM Goodman 2000; Chavanis 2011; Chavanis & Delfini 2011; Delgado & Muñoz Mateo 2022 and, more generally, self-interacting DM (Spergel & Steinhardt 2000; Rocha et al. 2013; Vogelsberger et al. 2012; Tulin et al. 2013; Ren et al. 2019; Correa et al. 2022), as well as phenomenological approaches (see e.g., Sánchez Almeida et al. 2020). However, these problems could also be the result of baryonic physics at play in galaxy formation and evolution within ΛCDM (Sales et al. 2022). For instance, the missing satellites could be due to a loss of efficiency in galaxy formation within low-mass halos as a result of the suppression of gas accretion during reionization (Efstathiou 1992; Bullock et al. 2000; Sawala et al. 2016). In addition, baryonic feedback from supernovae (SN) could provide a mechanism that is powerful enough to flatten CDM profiles at the core of dwarf galaxies (Navarro et al. 1996a; Pontzen & Governato 2012; Di Cintio et al. 2014; Sawala et al. 2016; Tollet et al. 2016; Chan et al. 2015; Fitts et al. 2017; Read et al. 2016a). The extent and reach of these mechanisms is yet unclear because it depends on modeling assumptions. In the case of SN feedback, there is a certain consensus that it is particularly efficient for galaxies with stellar masses in the range of 10⁸ M_⊙ ≲ M_⋆ ≲ 10⁹ M_⊙ (Tollet et al. 2016) and it could still play a role for masses all the way down to M_⋆ ∼ 10⁶ M_⊙¹.

The aim of this paper is to compare the predictions of FDM with the dynamics of dwarf galaxies, as observed through their HI rotation curves. While recent works in this direction (e.g., in Bernal et al. 2018; Bar et al. 2018, 2022; Robles et al. 2019; Chan & Fai Yeung 2021; Street et al. 2022; Khelashvili et al. 2023; Dave & Goswami 2023) have focused on the galaxies compiled by the SPARC data base (Lelli et al. 2016), we focus here on a specific group of nearby isolated low-mass irregular galaxies from the LITTLE THINGS survey (Oh et al. 2015). Our selection is based on rotation curves published in Iorio et al. (2017), Read et al. (2016b), which were obtained using state-of-the-art 3D reconstruction techniques to reanalyze the HI data cubes. Galaxies that could have a biased dynamical analysis, dubbed “rogues”, were flagged and excluded from the core data set. This catalog, referred to as LITTLE THINGS in 3D (Iorio et al. 2017), constitutes a robust and homogeneous collection of high-resolution rotation curves of galaxies that lie on the edge of galaxy formation. It provides an excellent laboratory to investigate the small-scale issues of ΛCDM and their potential interconnections with the nature of DM.

2. FDM generalities

We start by considering a real scalar field, ϕ(t, r), with a mass, m_a, satisfying a Klein-Gordon equation coupled to the gravitational field, Φ(t, r)². In the non-relativistic limit, it is convenient to decompose it as $ϕ (t, r) = 1 / \sqrt{2} m_{a} (e^{- i m_{a} t} ψ (t, r) + c . c .)$ $\phi(t,\boldsymbol{r})=1/\sqrt{2}m_{\mathrm{a}}\left(e^{-im_{\mathrm{a}}t}\psi(t,\boldsymbol{r})+\mathrm{c.c.}\right)$ in terms of a new complex scalar field ψ(t, r) and its complex conjugate (c.c.) that we assume to be slowly varying in space and time, |∇ψ|≪m|ψ| and $| \dot{ψ} | ≪ m | ψ |$ $\vert\dot{\psi}\vert\ll m\vert\psi\vert$ . The field ψ verifies, then, the Schrödinger-Poisson (SP) system of classical wave equations,

$\begin{matrix} i \frac{\partial}{\partial t} ψ (t, r) & = - \frac{1}{2 m_{a}} \nabla^{2} ψ (t, r) + m_{a} Φ (t, r) ψ (t, r), \\ \nabla^{2} Φ (t, r) & = {4 π G | ψ (t, r) |}^{2} . \end{matrix}$ $\begin{aligned} i\frac{\partial }{\partial t}\psi (t,\boldsymbol{r})&= -\frac{1}{2m_{\rm a}}\nabla ^{2}\psi (t,\boldsymbol{r})+m_{\rm a}\Phi (t,\boldsymbol{r})\psi (t,\boldsymbol{r}), \nonumber \\ \nabla ^{2}\Phi (t,\boldsymbol{r})&=4\pi G \vert \psi (t,\boldsymbol{r})\vert ^{2}. \end{aligned}$ (1)

Space and time dependence of the field factorize in the quasi-stationary solutions of these equations,

$\begin{matrix} ψ (r, t) = (\frac{m_{a}}{\sqrt{4 π G}}) exp (- i γ m_{a} t) χ (r) . \end{matrix}$ $\begin{aligned} \psi (\boldsymbol{r},t)=\left(\frac{m_{\rm a}}{\sqrt{4\pi G}}\right)\exp (-i\gamma m_{\rm a} t)\chi (\boldsymbol{r}). \end{aligned}$ (2)

We have chosen the normalization such that the physical density of the field is

$\begin{matrix} ρ (r) = \frac{m_{a}^{2} m_{P}^{2}}{4 π} χ^{2} (r), \end{matrix}$ $\begin{aligned} \rho (\boldsymbol{r}) = \frac{m^{2}_{\rm a}m_P^2}{4\pi }\chi ^{2} (\boldsymbol{r}), \end{aligned}$ (3)

and have introduced the Planck mass m_P = G^−1/2 ≈ 10¹⁹ GeV for convenience. Assuming spherical symmetry, χ obeys a simpler version of the SP equations,

$\begin{matrix} \partial_{x}^{2} [x χ (x)] & = 2 x (Φ (x) - γ) χ (x), \\ \partial_{x}^{2} [x Φ (x)] & = x χ^{2} (x), \end{matrix}$ $\begin{aligned} \partial _x^2[x\chi (x)]&=2x(\Phi (x)-\gamma )\chi (x),\nonumber \\ \partial _x^2[x\Phi (x)]&=x\chi ^2(x), \end{aligned}$ (4)

where we have rescaled the radial extension as x = m_ar (x is dimensionless, while r is the usual radial coordinate). The solutions to these equations that are regular at x = 0 and approaching 0 at x → ∞ are obtained from solving an eigenvalue problem for the parameter γ, which encodes the energy per unit mass. In particular, by integrating the energy density that can be derived from Eqs. (1) over a virialized density distribution, it is easy to show that γ = 3E/M, where E, M are the energy and mass of the distribution respectively (Bar et al. 2018). We also note that the expressions in Eq. (4) are invariant with respect to a scale transformation,

$\begin{matrix} (χ (x), Φ (x), γ) \to λ^{2} (χ (λ x), Φ (λ x), γ), \end{matrix}$ $\begin{aligned} \Big (\,\chi (x),\,\Phi (x),\,\gamma \Big )\rightarrow \lambda ^2\,\Big (\,\chi (\lambda x),\,\Phi (\lambda x),\,\gamma \Big ), \end{aligned}$ (5)

and we can set χ(0) = 1 without loss of generality.

The gravitationally bound eigenstates require γ < 0 and they are characterized by their number of zeros, n, with the ground-state solution corresponding to n = 0 and γ ≈ −0.69. Equation (4) needs to be determined numerically, but these expressions admit an approximate analytical solution around x = 0. In particular, using Moroz et al. (1998), we have:

$\begin{matrix} χ (r) = 1 - \frac{1}{3} (γ - Φ (0)) x^{2} + O (x^{4}) . \end{matrix}$ $\begin{aligned} \chi (r)=1-\frac{1}{3}(\gamma -\Phi (0))x^2+\mathcal{O} (x^4). \end{aligned}$ (6)

Hence, the resulting density distributions are flat at the center or “cored.” That of the ground state solution is referred to as a “soliton” in the literature. For completeness, we show in Fig. 1 the solution to Eqs. (4) χ(x) for the soliton and the two first excited states.

Fig. 1.

Solutions to Eq. (4) for the ground state and first two excitations, where the scale invariance was used to set χ(0) = 1.

The FDM corresponds, then, to a realization of the solutions to the SP equations whose mass density at the center is dominated by the soliton. This is determined by the axion mass, m_a, and a rescaling parameter λ that can be fixed using a physical property of the halo, such as the central density, via Eq. (3), the size, or the mass of the core. Interestingly, this implies that we obtain “soliton scaling relations” between different observable properties of the cores for a given mass of the axion. In particular, we obtain:

$\begin{matrix} M_{c} & ≃ 5.4 \times 10^{9} (\frac{m_{a}}{10^{- 23} eV})^{- 2} (\frac{r_{c}}{kpc})^{- 1} M_{⊙}, \end{matrix}$ $\begin{aligned} M_{\rm c}&\simeq 5.4 \times 10^9\,\bigg (\frac{m_{\rm a}}{10^{-23}\,{\text{ eV}}} \bigg )^{-2} \bigg (\frac{r_{\rm c}}{{\text{ kpc}}} \bigg )^{-1}{{M}_{\odot }}, \end{aligned}$ (7)

$\begin{matrix} ρ_{c} & ≃ 1.9 (\frac{m_{a}}{10^{- 23} eV})^{- 2} (\frac{r_{c}}{kpc})^{- 4} M_{⊙} {pc}^{- 3}, \end{matrix}$ $\begin{aligned} \rho _{\rm c}&\simeq 1.9\,\bigg (\frac{m_{\rm a}}{10^{-23}\,{\text{ eV}}} \bigg )^{-2} \bigg (\frac{r_{\rm c}}{{\text{ kpc}}} \bigg )^{-4}~~{{M}_{\odot }}\,{\text{ pc}}^{-3}, \end{aligned}$ (8)

where r_c is the core radius defined as the radius where the density drops by a factor of 2, M_c is the core mass, which is the mass enclosed within r_c, and ρ_c is the central density. Therefore, for a fixed m_a the soliton is heavier and denser the smaller the corresponding core. Lastly, the soliton solution can be approximated by an analytical profile that is obtained by fitting to results of simulations as Schive et al. (2014b):

$\begin{matrix} ρ_{sol} (r) \approx \frac{ρ_{c}}{[1 + 9.1 \times 10^{- 2} {(r / r_{c})}^{2}]^{8}} . \end{matrix}$ $\begin{aligned} \rho _{\text{sol}}(r) \approx \frac{\rho _{\rm c}}{\big [1 + 9.1 \times 10^{-2} (r/r_{\rm c})^2 \big ]^8}. \end{aligned}$ (9)

The FDM can, in principle, be composed of different eigenstates of the SP equations but they will all eventually decay down to the solitonic ground state by expelling scalar field to infinity (Seidel & Suen 1994). This process would relax the FDM distribution relatively quickly inside of a certain radius within the halo (Guzman & Urena-Lopez 2006; Schwabe et al. 2016; Hui et al. 2017; Li et al. 2021), while the outer region is expected to form a virialized halo. This configuration of the FDM has been confirmed by simulations (Schive et al. 2014a,b; Veltmaat et al. 2018; Mocz et al. 2018; Nori & Baldi 2021; Chan et al. 2022). The outer halo follows the universal Navarro-Frenk-White (NFW) density profile characteristic of CDM (Navarro et al. 1997):

$\begin{matrix} ρ_{NFW} (r) = \frac{ρ_{s}}{\frac{r}{r_{s}} {(1 + \frac{r}{r_{s}})}^{2}}, \end{matrix}$ $\begin{aligned} \rho _{\rm NFW}(r) = \frac{\rho _{\rm s}}{\frac{r}{r_{\rm s}}\left(1+\frac{r}{r_{\rm s}}\right)^{2}}, \end{aligned}$ (10)

where r_s and ρ_s are the scale radius and density, respectively. Equation (10) parametrizes the mass distribution of a CDM halo formed by the gravitational collapse of an initial overdensity that is spherically symmetric and homogeneous. Following the convention of Bryan & Norman (1998), the mass of the halo, M_h, is:

$\begin{matrix} M_{h} & \equiv \frac{4}{3} π r_{vir}^{3} ζ (z) ρ_{m 0}, \end{matrix}$ $\begin{aligned} M_{\rm h}&\equiv \frac{4}{3} \pi r_{\rm vir}^3 \zeta (z) \rho _{m0}, \end{aligned}$ (11)

$\begin{matrix} ζ (z) & \equiv \frac{18 π^{2} + 82 (Ω_{m} (z) - 1) - 39 {(Ω_{m} (z) - 1)}^{2}}{Ω_{m} (z)}, \end{matrix}$ $\begin{aligned} \zeta (z)&\equiv \frac{ 18 \pi ^2 + 82(\Omega _m(z) - 1) - 39(\Omega _m(z) - 1)^2}{\Omega _m(z)}, \end{aligned}$ (12)

where ρ_m0 denotes the cosmological background matter density (ρ_m0 = Ω_m0ρ_c0, with $ρ_{c 0} = 3 H_{0}^{2} / 8 G$ $\rho_{c0} = 3 H_0^2 /8 G$ ), and r_vir is defined as the radius where the average density falls to a critical threshold ζ(z)ρ_m0, marking the boundary of the halo. Throughout this paper, we set our cosmological parameters to the Planck 2015 cosmology results in Planck Collaboration XIII (2016), with h = 0.678, Ω_m0 = 0.308 (all our calculations are at z = 0)³. An additional characteristic property of the halo is the concentration parameter which is defined as c = r_vir/r_s and can be related to ρ_s in Eq. (10) by using

$\begin{matrix} ρ_{s} = \frac{ξ (z) ρ_{m 0}}{3} \frac{c^{3}}{log (1 + c) - \frac{c}{1 + c}} . \end{matrix}$ $\begin{aligned} \rho _{\rm s}=\frac{\xi (z)\rho _{m0}}{3}\frac{c^3}{\log (1+c)-\frac{c}{1+c}}. \end{aligned}$ (13)

Another important result of the simulations has been finding a relation that links the mass of the soliton to the mass of the host halo. This “core-halo relation” was originally found by Schive et al. (2014a,b) and can be written as

$\begin{matrix} M_{c} \approx 3.1 \times 10^{9} {(1 + z)}^{1 / 2} (\frac{m_{a}}{10^{- 23} eV})^{- 1} (\frac{M_{h}}{10^{12} M_{⊙}})^{1 / 3} M_{⊙}, \end{matrix}$ $\begin{aligned} M_{\rm c} \approx 3.1\times 10^9 (1+z)^{1/2} \bigg (\frac{m_{\rm a}}{10^{-23}\,\text{ eV}}\bigg )^{-1} \bigg (\frac{M_{\rm h}}{10^{12}\,{M}_{\odot }}\bigg )^{1/3}\,{M}_{\odot }, \end{aligned}$ (14)

which depends on the redshift, z. At the moment there is no comprehensive understanding of the physical mechanism underlying the core-halo relation, although some ideas have been put forward in Schive et al. (2014b), Hui et al. (2017), Bar et al. (2018), Eggemeier & Niemeyer (2019). In fact, subsequent simulations have found that there is a significant scatter that was studied in detail in Chan et al. (2022) by using a large sample of cosmological and soliton-merger simulations. This resulted in a more general relation that at z = 0 can be parametrized as:

$\begin{matrix} M_{c} = β (\frac{m_{a}}{8 \times 10^{- 23} eV})^{- 3 / 2} + (\frac{M_{h}}{γ})^{α} (\frac{m_{a}}{8 \times 10^{- 23} eV})^{3 (α - 1) / 2} M_{⊙}, \end{matrix}$ $\begin{aligned} M_{\rm c} = \beta \,\bigg (\frac{m_{\rm a}}{8 \times 10^{-23}\,\text{ eV}} \bigg )^{-3/2} + \bigg ( \frac{M_{\rm h}}{\gamma } \bigg )^{\alpha }\bigg (\frac{m_{\rm a}}{8 \times 10^{-23}\, \text{ eV}}\bigg )^{3(\alpha - 1)/2}\,{M}_{\odot }, \end{aligned}$ (15)

where $β = 8 . 00_{- 6.00}^{+ 0.52} \times 10^{6} M_{⊙}, {log}_{10} (γ / M_{⊙}) = - 5 . 73_{- 8.38}^{+ 2.38}$ $\beta = 8.00^{+0.52}_{-6.00} \times 10^6\,{M}_{\odot}, \; \log_{10} ( \gamma/{M}_{\odot} ) = -5.73 ^{+2.38}_{-8.38} \;$ and $α = 0 . 515_{- 0.189}^{+ 0.130}$ $\alpha = 0.515^{+0.130}_{-0.189}$ are the best-fit parameters obtained in Chan et al. (2022) with their respective errors determining the scatter.

Another feature of FDM halos is that the characteristic flatness in their density profiles is not limited to the core region where the soliton dominates, but is also reflected in their NFW-like tails at greater radii. This manifests in the form of a smaller concentration parameter compared to the typical CDM halo. This leads to a modified concentration – mass relation characterized by a suppression on lower mass scales, which, following Dentler et al. (2022), can be parametrized as:

$\begin{matrix} c^{FDM} (M_{h}) = [1 + γ_{1} (\frac{M_{0}}{M_{h}})]^{- γ_{2}} c^{CDM} (M_{h}), \end{matrix}$ $\begin{aligned} c^{\text{ FDM}}(M_{\rm h}) = \bigg [1 + \gamma _1 \bigg (\frac{M_0}{M_{\rm h}}\bigg )\bigg ]^{-\gamma _2} c^{\text{ CDM}}(M_{\rm h}), \end{aligned}$ (16)

where γ₁ = 15, γ₂ = 0.3 and we have the characteristic mass scale M₀ = 1.6 × 10¹⁰(m_a/10⁻²² eV)^−4/3 M_⊙, while c^CDM(M_h) denotes the concentration parameter in CDM. As we show later in this work, the CDM relation is modeled using the relation from Dutton & Macciò (2014) as:

$\begin{matrix} {log}_{10} c^{CDM} = 1.025 - 0.097 {log}_{10} (M_{h} / [10^{12} h^{- 1} M_{⊙}]), \end{matrix}$ $\begin{aligned} \log _{10} c^{\text{ CDM}} = 1.025 - 0.097 \log _{10} \Big (M_{\rm h} / [10^{12}\,h^{-1}\,{{M}_{\odot }}] \Big ), \end{aligned}$ (17)

with a scatter of 0.16 dex⁴.

As a point of clarification, we note that the concentration parameter is defined in Eqs. (10) and (13) with respect to the NFW component of the density profile. This is irrespective of whether the point r = r_s, which can be generally defined as the point where the logarithmic slope equals −2, is located in the NFW-like region or not, meaning that the two definitions are not necessarily equivalent. Nevertheless, as we discuss in Sect. 3 for the halos that we model in this work, we find that the median values of r_s lie in the NFW-like regions and the two definitions are consistent with each other in our case.

A final property that is characteristic of FDM and is particularly useful in the context of small-scale structure issues of ΛCDM is the suppression of the primordial power spectrum that inhibits the formation of halos below a certain mass scale (Khlopov et al. 1985; Hu et al. 2000). Moreover, the subhalo population is expected to be further reduced by dynamical effects such as tidal disruption (Hui et al. 2017). This offers an alternative “exotic” solution to the missing satellites problem and, by extension, to the associated too big to fail problem. One observable statistic that is relevant to test this suppression of formation of structures at small scales is the halo mass function (HMF) for low halo masses of M_h ≲ 10¹⁰ M_⊙. This can be derived by a simple parametrization derived from numerical simulations (Schive et al. 2016; May & Springel 2021, 2023):

$\begin{matrix} {\frac{d n}{d M_{h}} |}_{FDM} = {(1 + {(\frac{M_{0}}{M_{h}})}^{- α})}^{- β} {\frac{d n}{d M_{h}} |}_{CDM}, \end{matrix}$ $\begin{aligned} \left.\frac{\mathrm{d}n}{\mathrm{d}M_{\rm h}}\right|_{\rm FDM}=\left(1+\left(\frac{M_{0}}{M_{\rm h}}\right)^{-\alpha }\right)^{-\beta } \left.\frac{\mathrm{d}n}{\mathrm{d}M_{\rm h}}\right|_{\rm CDM}, \end{aligned}$ (18)

where, as previously for the concentration – mass relation, M₀ = 1.6 × 10¹⁰(m_a/10⁻²² eV)^−4/3 M_⊙, α = 1.1 and β = 2.2. The characteristic mass M₀ defines the regime below which the HMF starts dropping substantially.

3. Fitting FDM to the LITTLE THINGS in 3D

3.1. LITTLE THINGS in 3D

Local Irregulars That Trace Luminosity Extremes, The HI Nearby Galaxy Survey (LITTLE THINGS) is a Very Large Array HI survey which has produced high-resolution rotation curves for 26 nearby dwarf irregular galaxies within 11 Mpc of the local volume (Oh et al. 2015; Hunter et al. 2012). Most of these curves are able to resolve the inner part (within ∼1 kpc of the center) of the galaxies, making this survey an optimal laboratory to investigate the core-cusp problem. We considered a representative and robust subset of these galaxies, with stellar masses in the range of 10⁵ M_⊙ ≲ M_⋆ ≲ 10⁸ M_⊙, obtained in Iorio et al. (2017), Read et al. (2017) by reanalyzing the data cubes with 3D reconstruction techniques. We will refer to this as the LITTLE THINGS in 3D (LTs) sample. The data set consists of circular velocities after the observed velocities have been corrected for the pressure support of the random gas motions⁵. A fraction of these galaxies termed “rogues” were flagged and excluded from the core LTs data set because they were nearly face-on (five inclination rogues), their measured distance was not well-determined (one distance rogue) or because they presented indications of being far from equilibrium (two disequilibrium rogues). Our fiducial cosmological analysis will include 11 of the remaining galaxies (first corresponding lines in Table 2 of Read et al. 2017). For these galaxies, we will consider the measured rotation velocities in the whole measured range, which are quite regular even if some of them manifest the effect of HI holes, which could indicate recent (fast expanding) or old (quiescent) stellar activity. The disequilibrium rogues present irregular rotation curves while the distance rogue, DDO 101, could have a stellar mass quite above M_⋆ ∼ 10⁸ M_⊙, depending on the assumed distance, see Read et al. (2016b). Therefore we did not include any of these galaxies in our analysis. Finally, we did not include the inclination rogues in our fiducial analysis either, but we do show the results obtained for DDO 133 and DDO 50 for illustration purposes later in this work⁶. We also note that some galaxies in our sample have significant errors in the inclination values. However, this is self-consistently accounted for in the rotation curve errors (Iorio et al. 2017), so that it is unnecessary (and potentially inconsistent) to use the inclination angle as a fitting parameter for LTs.

3.2. Fitting strategy

The FDM density model is a piecewise function of the soliton and NFW profiles matched at a transition radius r_t, namely, ρ_sol(r_t) = ρ_NFW(r_t). In this model, the soliton describes the inner part of the galaxy, while a NFW halo does the outer part:

$\begin{matrix} ρ_{FDM} (r) = {\begin{matrix} ρ_{sol} (r) r \leq r_{t}, \\ ρ_{NFW} (r) r > r_{t} . \end{matrix} \end{matrix}$ $\begin{aligned} \rho _{\rm FDM}(r)={\left\{ \begin{array}{ll} \rho _{\rm sol}(r) r\le r_{\rm t},\\ \rho _{\rm NFW}(r) r> r_{\rm t}. \end{array}\right.} \end{aligned}$ (19)

We use the numerical ground-state solution to the SP Eqs. (4) for ρ_sol(r) and the NFW density profile in Eq. (10) for ρ_NFW(r). We note that the numerical solution for ρ_sol(r) that we use is very similar to the phenomenological function in Eq. (9) obtained in Schive et al. (2014a,b) from simulations and which is the basis for most of the analyses of FDM.

Overall, the FDM density profile has four free parameters that we chose to be: (1) the axion mass m_a, (2) the solitonic core’s mass M_c, (3) the concentration parameter c and, (4) the halo mass, M_h, defined as in Eq. (11). Other parameters such as the core radius, r_c, the transition radius, r_t, the central density, ρ_c, or the virial radius r_vir can be obtained from these four parameters. In particular, the transition radius, r_t, can be determined by searching the first point in r, where ρ_sol(r) = ρ_NFW(r) starting from r → ∞.

There are two major approximations when using the FDM density profile in Eq. (19). Firstly, simulations have shown that the central soliton can still be affected significantly by interference with residual excited states of the SP equations. This produces 𝒪(1) oscillations in the central density, which is an effect not captured by the simple solitonic distribution (Veltmaat et al. 2018; Li et al. 2021). This also includes contributions from non spherically symmetric eigenstates of the SP equations, as described in Vicens et al. (2018), Li et al. (2021), for instance. Secondly, the SP expressions in Eq. (1) describe configurations of the FDM field neglecting the contribution of the baryons to the gravitational potential. The baryonic contributions, including baryonic SN feedback, have been recently taken into account in simulations (Chan et al. 2018; Mocz et al. 2019, 2020; Veltmaat et al. 2020) with the general result that, for the same given m_a, a soliton of a given mass is still formed but more compressed and denser than in its pure FDM configuration. This effect should be especially relevant for galaxies whose central masses are dominated by the baryons and it can be estimated semi-analytically as discussed in Bar et al. (2018, 2019a) and followed in Sect. 5.

Moreover, the structure of the boundary between the soliton and the halo at r_t can be modeled differently. This can have a non-trivial effect (Bernal et al. 2018; Vicens et al. 2018) especially if it introduces further constraints between the four parameters of the basic model.

The theoretical radial profile of the total circular velocity squared is given by the sum of the various mass contributions to the gravitational potential,

$\begin{matrix} V_{th}^{2} (r ; θ_{FDM} ; M_{⋆}) = V_{FDM}^{2} (r ; θ_{FDM}) + V_{⋆}^{2} (r ; M_{⋆}) + V_{gas}^{2} (r) . \end{matrix}$ $\begin{aligned} V_{\rm th}^{2}(r;\,\boldsymbol{\theta }_{\rm FDM};{M_{\star }})= V_{\rm FDM}^{2}(r;\,\boldsymbol{\theta }_{\rm FDM})+V_{\star }^{2}(r;\,{M_{\star }})+V_{\rm gas}^{2}(r). \end{aligned}$ (20)

In this equation $V_{FDM}^{2} (r; θ_{FDM}) = G M_{FDM} (r; θ_{FDM}) / r$ $V^2_{\rm FDM}(r;\,\boldsymbol{\theta}_{\rm FDM})=GM_{\rm FDM}(r;\,\boldsymbol{\theta}_{\rm FDM})/r$ , where M_FDM(r; θ_FDM) is the enclosed mass within radius r obtained from the density distribution in Eq. (19) and θ_FDM = (m_a, M_c, c, M_h) is the four-vector of FDM fitted parameters. To describe the baryonic components in Eq. (20), we follow the reference Read et al. (2017) and model the distributions of gas and stars as exponentially thin discs (Binney & Tremaine 2008),

$\begin{matrix} V_{⋆ / gas}^{2} = \frac{G M_{⋆ / gas}}{R_{⋆ / gas}} I (y), \end{matrix}$ $\begin{aligned} V^2_{\star /\mathrm{gas}}=\frac{GM_{\star /\mathrm{gas}}}{R_{\star /\mathrm{gas}}}\mathcal{I} ({ y}), \end{aligned}$ (21)

with M_⋆/gas and R_*/gas being the mass and the exponential scale length of the disc, respectively, and with I(y) = y²[I₀(y)K₀(y)−I₁(y)K₁(y)], where I_0, 1(y) and K_0, 1(y) are Bessel functions that depend on the dimensionless radius parameter y = r/(2R_⋆/gas). We fixed the values of R_⋆, R_gas, and M_gas to the values reported in Read et al. (2017) and we added the stellar mass M_⋆ as a new parameter in the fit.

The theoretical velocity is compared to the observed rotation curves using a log-likelihood

$\begin{matrix} ln L (θ_{FDM}, M_{⋆}) = - \frac{1}{2} \sum_{i} {(\frac{V_{obs} (r_{i}) - V_{th} (r_{i} ; θ_{FDM} ; M_{⋆})}{δ V_{obs} (r_{i})})}^{2}, \end{matrix}$ $\begin{aligned} \ln \mathcal{L} (\boldsymbol{\theta }_{\rm FDM},{M_{\star }})=-\frac{1}{2}\sum _{i}\left(\frac{V_{\rm obs}(r_{i})-V_{\mathrm{th}}(r_i;\, \boldsymbol{\theta }_{\rm FDM};\,{M_{\star }})}{\delta V_{\rm obs}(r_{i})}\right)^2, \end{aligned}$ (22)

(up to an additive constant), where V_obs(r_i) and δV_obs(r_i) are the velocity and corresponding uncertainty of the i-th data point in a galaxy’s rotation curve (the summation is over the radial data bins). We also assumed uncorrelated Gaussian errors for the data points distributions. Finally, for the inclination rogues, DDO 50 and DDO 133, we also fit the inclination as an additional parameter (see Appendix A.2 for details).

To effectively scan the parameter space of the FDM model, we used the emcee package, an affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) (Foreman-Mackey et al. 2013), based on the algorithm proposed by Goodman & Weare (2013). To deal with possible multi-modal posteriors and enhance sampling efficiency, we implemented a combination of MCMC moves, which significantly reduces autocorrelation times⁷. We also initialized MCMCs under different local optima to ensure that the posterior distributions converged to the global one. We ran the MCMCs with 48 walkers for a total of 20 000 steps, discarding the first 10 000 as a conservative burn-in period. We found these to produce consistent results under different runs, as well as meeting the criterion of exceeding 100 autocorrelation times, with acceptance fractions in the ∼0.2 − 0.4 range, in agreement with the criteria established in Foreman-Mackey et al. (2013).

The priors used for the initial analysis in this study (where the FDM mass is allowed to vary for each galaxy) are summarized in Table 1. The prior on M_⋆ is designed to yield results consistent with stellar population synthesis models, as reported in Read et al. (2016b), using the same central values and defining the errors given as the 68% CL bounds under a log-normal distribution. A similar type of prior was used for the angle i for inclination rogues, using a normal distribution instead of a log-normal one. Priors are generally designed to be very loose, covering ranges well beyond expected values and have, for the most part, little effect on the fits. The less trivial priors that are able to meaningfully constrain the parameter space of the FDM distribution, besides the priors on M_⋆, are the core-halo relation from Eq. (15) and demanding r_t ≥ r_c, (when r_t < r_max, where r_max is the maximum radius where rotation curves are measured), which are both expected from simulations. Fits with r_t ≥ r_max would yield identical soliton-dominated fits (for a given m_a and M_c). For the core-halo relation we set z = 0 using Eq. (15), where we have assumed that the uncertainty on the formation history of the halo, see Burkert (2020), is within the intrinsic scatter comprised by this formula.

Table 1.

Priors used to fit the DM and baryon contributions to the rotation curve data from the LTs catalog.

4. Results

4.1. Fits to rotation curves

The results of the fits to the rotation curves of the 11 LTs galaxies and the 2 selected inclination rogues (DDO 50 and DDO 133) are shown in Table 2. For each galaxy we list the median and 68% CL region of the posterior distributions of the five fitted parameters, M_⋆ and θ_FDM, the reduced chi-square and a subset of derived quantities characterizing the FDM distribution, including the core radius, r_c, the transition radius, r_t, and the central density, ρ_c. In Fig. 2, we show the results of the fits for two representative galaxies in our sample, WLM and DDO 154. Besides the data, we show the median, the 68% and 95% CL of the posterior distribution of the total theoretical rotation curve and the individual contributions from the FDM and baryonic density distributions for the central values of the parameters. Further details, including the corner plots and mean autocorrelation times of the MCMC for the fits of these two galaxies are shown for reference in the Appendix A. The full set of rotation curves, corner plots, and mean autocorrelation times obtained from the fits in Table 2 are publicly available under the category Variable Axion Mass in the FuzzyDM GitHub repository⁸.

Fig. 2.

Rotation curves in FDM for two representative galaxies in the LTs catalog, WLM, and DDO 154. We plot the DM, gas, and stellar components obtained from the fits (maximum posterior values) and compared to the data. The red dashed lines are the median of the posterior distribution for the total velocity and the colored band is the corresponding uncertainty at 68% and 95% CL.

Table 2.

Median and 68% confidence level (CL) of the posterior distributions of the parameters m_a, M_c, c, M_h, and M_⋆, and the reduced chi-square of the maximum posterior fit to the LTs catalog, including some of the rogues.

As shown by the reduced chi-square values obtained in the fits, the performance of the FDM model describing the rotation curve data is, in general, excellent. In two cases, DDO 52 and UGC 8508, the value of the chi-square is too low suggesting there is overfitting or overestimated errors in the data. Our fits favor large solitons and with transition radii to the halo typically overlapping with r_max. This reduces the power of the data to constrain the halo parameters c and M_h, explaining their large uncertainties. The direct relation between the scale density of the NFW distribution and c in Eq. (13) imposes an upper bound on this parameter because, given a soliton describing the core of a galaxy, the density at r ≤ r_t cannot be arbitrarily large. In case of the halo mass an increased constraining power is achieved by imposing the core-halo relation in the priors. On the other hand, the parameters of the soliton, m_a and M_c, are quite well determined, with relative uncertainties even reaching 𝒪(10%) in some cases. The sensitivity to the axion mass of these data is illustrated in Fig. 3 for the rotation curves of WLM and DDO 154. The fit is also compared to the to a pure NFW profile that is included for reference. The description of the data can be quite insufficent, and even worse than for NFW for some choices of m_a. In these cases, the fit does not converge to the one of NFW (that would result from suppressing the soliton by effectively running r_t to zero) because the soliton cannot be arbitrarily light due to the soliton-halo relation imposed in the priors.

Fig. 3.

Rotation curves in FDM for two representative galaxies in the LTs catalog: WLM and DDO 154. We plot the total circular velocity and reduced chi-squares resulting from the maximum posterior fits to the data when fixing the axion mass at different values. We compare the results also to the fit to a NFW profile.

As illustrated in Fig. 2, the rotation curves are dominated by the DM contribution in the full radial domain, although the baryonic contribution can also be significant in the central regions of DDO 210, NGC 2366, DDO 168, and NGC 6822. In fact, NGC 6822 is the only galaxy for which the FDM presents a relatively poor fit. The main tension arises from the two innermost points (at ∼70 pc and ∼140 pc), where the velocities can be extracted with a relatively high precision. For the nominal stellar mass of NGC 6822, M_⋆ = (7.6 ± 1.9) × 10⁷ M_⊙ (and fixed scale radius R_⋆), the mass distribution at the center of the galaxy is dominated by stars and the FDM model cannot fit the rotation curve without significantly reducing the stellar content down to M_⋆ ≃ 1.8 × 10⁷ M_⊙. In this context, it is important to note that there is independent evidence for the baryonic dominance of the mass at the core of NGC 6822 based on stellar kinematics (Kirby et al. 2014). The difficulty of our model to account satisfactorily for these data does not pose, a priori, a serious problem to FDM because the central densities of the soliton can be significantly affected by interference effects, as discussed above. Apart from this, a fast-expanding HI bubble was identified in this galaxy for radii below 2.5 kpc (Read et al. 2017). For these reasons, we removed NGC 6822 from our core data sample, adding it to the list of rogue galaxies in Table 2 and for the purposes of the subsequent analyses. In the next subsections, we present different consistency and diagnostic tests on the behavior of the parameters obtained in the fitting procedure.

4.2. The axion mass and the soliton scaling relations

In our initial analysis, we fit the fundamental axion mass individually for each galaxy, allowing for a broad range of values, from 10⁻²⁵ eV to 10⁻¹⁹ eV. As shown in Table 2 and in the left panel of Fig. 4, the obtained values of m_a lie within the range of ∼10⁻²³ eV − 10⁻²² eV. In particular, the values for the most constraining galaxies are different, at most, by a factor of 2 − 3, which constitutes a remarkable consistency check for FDM with the LTs data sample. Moreover, the values of m_a extracted from the selected rogues are also within this range of masses.

Fig. 4.

Axion masses obtained from the posterior distributions of the fits to the rotation curves of the LTs galaxies where we show medians and 68% and 95% CL errors. The 10 galaxies entering our fiducial analysis are represented by circles while the rogues are represented by triangles. The dot-dashed line and bands represent the mean, 68% and 95% CL regions of a weighted average to these determinations (see Eq. (23) and associated discussion). Left: Galaxies are ordered by the central value of the stellar mass reported in Read et al. (2017), which where found to be in good agreement with our posterior distributions. Stellar masses are increasing from left to right. Right: Empirical fit to a power law between m_a and M_⋆ showing evidence for a correlation between the two as extracted from the galaxies in our sample. The red solid line denotes the maximum posterior fit with the 68% and 95% CL regions.

We can obtain the preferred value of m_a by taking a weighted average of log₁₀m_a for the ten galaxies in our sample (excluding the rogues), which is the parameter actually determined in the fits and whose posterior distributions are approximately symmetric around the mean. We then obtain m_a = 1.90 ± 0.08 × 10⁻²³ eV, corresponding to an uncertainty of 0.019 dex, and a reduced chi-square $χ_{ν}^{2} = 7.68$ $\chi^2_\nu=7.68$ . This relatively poor value stems from the small uncertainties derived for m_a in some of the galaxies (see Fig. 2) and it could be due to unknown systematic uncertainties in the data or the model. Under this premise, there are several ways to account for them in the analysis (Workman et al. 2022; Hogg et al. 2010). We follow here the prescription followed by the particle data group for averages of experimental results in this case, which consists in increasing the final uncertainty by a scale factor $S = {(χ_{ν}^{2})}^{1 / 2}$ $S=(\chi^2_\nu)^{1/2}$ (Workman et al. 2022). With this method we obtain a more conservative error:

$\begin{matrix} Optimal axion mass: m_{a} = 1 . 90_{- 0.21}^{+ 0.24} \times 10^{- 23} eV, \end{matrix}$ $\begin{aligned} \text{ Optimal} \text{ axion} \text{ mass:} m_{\rm a} = 1.90^{+0.24}_{-0.21}\times 10^{-23}\,\text{ eV}, \end{aligned}$ (23)

corresponding to an uncertainty of 0.052 dex. For completeness we have also obtained an alternative determination of this optimal mass by performing a global simultaneous fit to the parameters involving the rotation curves of all ten galaxies in the sample and using the same m_a. This procedure is more robust as it accounts for possible correlations and features in the fit and posterior distributions. Nonetheless, we obtain the same central value although with a relatively smaller uncertainty. The full set of results obtained from the fits using Eq. (23) are publicly available under the category Benchmark analysis in the FuzzyDM GitHub repository⁹.

On the other hand, the values of m_a determined from the different galaxies do not seem to be randomly scattered around the optimal value but rather following a negative correlation with the stellar mass. Specifically, the lower the M_⋆ value, the higher these m_a values tend to be. This correlation is visible in the left panel of Fig. 4, where we have ordered the values of m_a according to increasing values of the stellar mass of the galaxies. Since m_a is a fundamental constant in the model, this suggests that there are astrophysical processes influencing the observed cores beyond the structure predicted by a FDM soliton with the optimal mass in Eq. (23). We can try to quantify the statistical significance of the correlation by introducing an empirical power law between the two quantities, $m_{a} = α_{⋆} M_{⋆}^{β_{⋆}}$ $m_{\rm a}=\alpha_\star\,{M_{\star}}^{\,\beta_\star}$ . We then performed a linear fit of $α_{⋆}^{'} = {log}_{10} α_{⋆}$ $\alpha^\prime_\star=\log_{10}\alpha_\star$ and β_⋆ to the 10 (log₁₀M_⋆, log₁₀m_a) data points, excluding the three rogue galaxies we have studied. For this, we use MCMCs and also include a possible source of unknown systematic uncertainties in the data points as described in Sect. 8 of Hogg et al. (2010). The results of this fit are shown in the right panel of Fig. 4. We find that $β_{⋆} = - 0 . 62_{- 0.25}^{+ 0.21}$ $\beta_\star=-0.62^{+0.21}_{-0.25}$ , corresponding to evidence at a ∼3.3σ level (with respect to β_⋆ = 0) of a correlation between the values of m_a determined from the rotation curves of the galaxies and the relevant M_⋆. We describe the implementation of these fits and the interpretation of the CL regions in Appendix A. As shown in Fig. 4 the two rogue galaxies NGC 6882 and DDO 133 seem to also follow the correlation, while DDO 50 is an outlier. We discuss the possible implications of this correlation in more detail in Sect. 5.

Another perspective on the consistency of the inferred axion masses can be obtained by testing the soliton scaling relations Eqs. (7) and (8). To quantify the agreement of the model with observations, we fit our results of ρ_c and M_c to power laws in r_c with exponents β_ρ and β_M, respectively. These fits are done again with MCMCs, including possible systematic uncertainties from Hogg et al. (2010), and they are described in Appendix A. In the left panel of Fig. 5, we present our results in the r_c − ρ_c plane. The data shows that the central densities fall with the size of the core, although at a pace much slower than predicted by FDM. Indeed, our fit to the data yields $β_{ρ} = - 1 . 18_{- 0.30}^{+ 0.22}$ $\beta_\rho=-1.18^{+0.22}_{-0.30}$ , disagreeing with the prediction of FDM, β_ρ = −4, with a significance larger than 3σ. It is remarkable that our result is compatible with β_ρ ≈ −1.3, obtained by Deng et al. (2018) using a different sample of galaxies analyzed with a Burkert profile Rodrigues et al. (2017). Our results also agree within 1σ with β_ρ ≈ −1, which is an empirical result derived from the fact that the so-called “column density,” Σ₀ = ρ_c × r_c (Donato et al. 2009), is nearly constant for a large sample of different galaxies (see also Kormendy & Freeman 2016; Burkert 2015; Karukes & Salucci 2017, and Di Paolo et al. 2019). In fact, we obtain $Σ_{0} = 51_{- 9}^{+ 12} M_{⊙} {pc}^{- 2}$ $\Sigma_0=51^{+12}_{-9}\,{{M}_{\odot}}\,\text{ pc}^{-2}$ which undershoots slightly the empirical value $Σ_{0} = 75_{- 45}^{+ 55} M_{⊙}$ $\Sigma_0=75^{+55}_{-45}\,{{M}_{\odot}}$ pc⁻² (Burkert 2020).

Fig. 5.

Soliton scaling relations compared to the results of the fits to LTs, where each point is colored according to the value of M_⋆ of the corresponding galaxy. Dot-dashed lines and gray bands represent the mean and the 68% and 95% CL regions for the correlation expected from FDM for the average axion mass in Eq. (23). The solid line and surrounding bands represent the maximum posterior and 68% and 95% CL regions of the fit to a power law of the results extracted from the individual galaxies. Left: Results in the r_c − ρ_c plane compared to Eq. (8) and where we also show the results obtained independently for a different sample of galaxies by using a Burkert profile (Rodrigues et al. 2017; Deng et al. 2018). Right: Results in the r_c − M_c plane compared to Eq. (8).

In the right panel of Fig. 5, we show our results in the r_c − M_c plane. In this case the data manifest a correlation which is almost orthogonal to the FDM prediction: The cores are heavier, the bigger (not the smaller) they are. Quantitatively we find $β_{M} = 0 . 45_{- 0.08}^{+ 0.10}$ $\beta_M=0.45^{+0.10}_{-0.08}$ which (at face value) disagrees with the FDM value β_M = −1 with a significance larger than 5σ. In addition, and as shown by the color code used in the data points of Fig. 5, r_c, and M_c are also correlated with M_⋆, as the heavier cores contain also more stars. This suggests again that there is a non-trivial interplay between the structure of the DM cores at the center of the dwarf galaxies in LTs and baryonic physics.

4.3. Properties of the halo

In addition to the properties of the solitons we can test the predictions of FDM regarding the halos they inhabit. To do this we set m_a to the optimal value found in Eq. (23), m_a ≈ 1.90 × 10⁻²³ eV, disregarding the difficulties of the model to explain the data discussed in the previous section. First, we show on the top-left panel of Fig. 6 our results in the M_h − M_c plane compared to the core-halo relation predicted from simulations in Eqs. (14) and (15). This relation is imposed as a prior in the fits and it provides the most important constraint on the halo masses for most of the galaxies. This is why the uncertainties of M_h span almost the entire allowed range for a given M_c. Next, we show on the top-right panel our results in the M_h − c plane comparing them to the predictions of CDM and FDM. The uncertainties on the concentration parameter c obtained from the fits are large although, as discussed in Sect. 4.1, this is constrained from above by the central densities reached by the soliton. This reduces systematically the values of c obtained from the data as compared to the predictions of CDM and agreeing with the expectations of FDM.

Fig. 6.

Halo properties obtained from the posteriors of the FDM fits to the rotation curves of the LTs galaxies for a fixed axion mass m_a = 1.9 × 10⁻²³ eV. Gray-colored points are galaxies for which both median and maximum posterior values of r_t are within the observed range, while gold-colored are soliton-dominated by virtue of not satisfying that condition. Uncertainties in the data points describe the credible intervals at 68% (solid) and 95% (dashed) CL. Uncertainties in the theoretical bands refer to the 1 and 2σ scatter in the predictions. Top-left: M_c − M_h relation compared to the region allowed by Eq. (15) and the prediction in Eq. (14) from Schive et al. (2014b; dotted line). Top-right: c − M_h relation compared to Eq. (16) and to the prediction in CDM, Eq. (17) from Dutton & Macciò (2014) with σ = 0.16 dex (gray dashed line). Bottom-left: Stellar-to-halo mass relation from abundance matching compared to the FDM prediction using Eq. (18), and the CDM predictions from Behroozi et al. (2019; dashed gray line) and Brook et al. (2014; dot-dashed line) with a reference σ = 0.30 dex scatter in both cases. Bottom-right: Expected number of halos with M_h ≲ 10¹¹ M_⊙ in a Local Group volume of radius 1.8 Mpc predicted in FDM as a function of m_a, using Eq. (18) and calibrated with the CDM prediction of the HMF from Brook et al. (2014). This is compared to the HMF from Schive et al. (2016; dashed line). For the CDM prediction, we take 10⁹ ≲ M_h ≲ 10¹¹ M_⊙. The dot-dashed line and gray bands represent the mean and the 68% and 95% CL regions for the axion mass in Eq. (23).

The most important test that can be performed on the properties of the FDM halos with our results is related to the abundance of dwarf galaxies. As discussed in Sect. 2, FDM predicts a suppression of structure formation at small scales that translates into a low-mass cutoff in the HMF. One classical method to test observationally the HMF is via abundance matching methods with the observed stellar mass function (SMF; Frenk et al. 1988; Vale & Ostriker 2004; Moster et al. 2010; Behroozi et al. 2013). Given the non-monotonic behavior of the HMF in Eq. (18), we follow a method inspired by Vale & Ostriker (2004), whereby we find a M_⋆ − M_h stellar-to-halo mass function by matching the cumulatives of the SMF and HMF,

$\begin{matrix} \int_{M_{⋆}}^{\infty} \frac{d n}{d M} |_{SMF} d M = \int_{M_{h}}^{\infty} \frac{d n}{d M} |_{HMF} d M . \end{matrix}$ $\begin{aligned} \int ^\infty _{{M_{\star }}} \frac{\mathrm{d}n}{\mathrm{d}M}\Big |_{\rm SMF} \mathrm{d}M=\int ^\infty _{M_{\rm h}} \frac{\mathrm{d}n}{\mathrm{d}M}\Big |_{\rm HMF} \mathrm{d}M. \end{aligned}$ (24)

In the bottom-left panel of Fig. 6 we show the results of abundance matching using the stellar-to-halo mass function reported in Behroozi et al. (2019) for CDM and corrected for FDM by Eq. (18) and the optimal axion mass in Eq. (23). We also add the prediction for CDM reported in Brook et al. (2014), which implements a HMF specific for the environment of the Local Group. The FDM prediction for the stellar-to-halo mass function features a flattening starting at values as high as M_h ≈ 5 × 10¹¹ M_⊙, which results from the strong suppression of structure formation in FDM for the axion mass value favored by the fits to the rotation curves of LTs. In particular, we observe that according to this prescription of abundance matching, galaxies with stellar masses M_⋆ ≲ 10^8.5 should be extremely rare, which is clearly in plain contradiction with the posterior distributions of the very same fits in the M_⋆ − M_h plane. On the contrary, the data seem to be more consistent with the stellar-to-halo mass function predicted by CDM (Behroozi et al. 2019; Brook et al. 2014).

A more transparent way to illustrate this tension between FDM and observations is by directly comparing the expected number of halos that are predicted to exist in a specified local volume to the corresponding number of observed galaxies. Or, inversely, we can translate the abundance of observed galaxies to a lower bound on the FDM mass that is necessary to make the existence of such galaxies statistically plausible in the model. We take as a reference the simulations of galaxy formation in the Local Group within the ΛCDM that are reported in Brook et al. (2014) using 1.8 Mpc for the radius of the Local Group Volume (LGV), which is the one we use. In the bottom-right panel of Fig. 6, we show the number of FDM halos with M_h ≲ 10¹¹ M_⊙ that are expected to be contained in this volume as a function of m_a. For this we use the HMF in Eq. (18) where we implement the Sheth-Tormen prediction for the HMF of CDM (Sheth & Tormen 1999) and we normalize it to reproduce the results of Brook et al. (2014) for the LGV (see Appendix B). For reference, we also show the abundance that is obtained using the HMF originally reported in Schive et al. (2016).

The predicted abundances are to be compared with the three galaxies contained in the LGV that we have found in our analysis. More generally, we should also compare these predictions with the properties of 116 such galaxies that can be extracted from the most recent census Putman et al. (2021)¹⁰. Using the former, we find a very conservative lower bound of m_a ≳ 4.3 × 10⁻²³ eV. Nonetheless, this alone implies that our benchmark value for m_a in Eq. (23) is in tension with observations with a significance greater than 5σ. Using the 116 galaxies from Putman et al. (2021), we can obtain a much more severe lower limit on the mass, namely, m_a ≳ 5.5 × 10⁻²² eV.

Finally, we can also study the population of dwarf galaxies in the Milky-Way’s halo using the subhalo mass function for FDM recently reported in Du (2018), Du et al. (2017). Using this, and focusing on the 59 galaxies from Putman et al. (2021) that lie within 400 kpc of the Milky Way (following Marsh & Niemeyer 2019 and which is roughly twice its virial radius $r_{vir}^{MW} \approx 200$ $r_{\mathrm{vir}}^{\mathrm{MW}}\approx200$ kpc (Bland-Hawthorn & Gerhard 2016)), we obtain the stronger lower bound m_a ≳ 7.3 × 10⁻²² eV.

4.4. Baryonic effects

Our analysis so far has ignored the effects of baryons and we may consider if a more realistic approach including them could address the difficulties of FDM to explain observational data that were discussed above. Intuitively, one might expect that in DM-dominated galaxies such as those in the LTs sample these effects should not lead to major changes in the analysis. However, as discussed in Sect. 1, the consequences of baryonic feedback from SN can be significant for CDM specifically for dwarf galaxies in the mass range we are investigating. In order to realistically account for these effects, we would require information outside the SP equation related to the baryonic components and including the modeling of star formation and evolution. Interestingly, FDM simulations of structure formation with baryonic feedback have recently been performed by Mocz et al. (2019, 2020), Veltmaat et al. (2020), finding that the inclusion of baryons in FDM halos does not impede their thermalization and formation of quasi-stationary (modified) solitons at their centers. This motivates the study of this solution even in the presence of baryonic effects.

A simple and minimal approach is to solve for the modified soliton solution under the assumption of a static background baryonic potential, which can be obtained in a very similar way to that of the original soliton. The resulting modified soliton solution has been studied in Bar et al. (2018; and subsequently in Bar et al. 2019b), which was found to accurately match simulations incorporating gravitational interactions with baryons in Chan et al. (2018). In order to proceed, we closely follow the approach of Bar et al. (2018). The methodology is very similar to that of the ordinary soliton in Eq. (1), except that in this case we introduce the baryonic density, ρ_b, solving for the total (baryonic + DM) potential Φ_tot and solve for a given central density ρ_c. The resulting modified SP equation is given by

$\begin{matrix} i \frac{\partial}{\partial t} ψ (t, r) = - \frac{1}{2 m_{a}} \nabla^{2} ψ (t, r) + m_{a} Φ_{tot} (t, r) ψ (t, r), \\ \nabla^{2} Φ_{tot} (t, r) = 4 π G ({| ψ (t, r) |}^{2} + ρ_{b} (r)) . \end{matrix}$ $\begin{aligned}&i\frac{\partial }{\partial t}\psi (t,\boldsymbol{r}) = -\frac{1}{2m_{\rm a}}\nabla ^{2}\psi (t,\boldsymbol{r}) +m_{\rm a}\Phi _{\text{tot}}(t,\boldsymbol{r})\psi (t,\boldsymbol{r}), \nonumber \\&\nabla ^{2}\Phi _{\text{tot}}(t,\boldsymbol{r})= 4\pi G \big (\vert \psi (t,\boldsymbol{r})\vert ^{2} + \rho _{\rm b}(\boldsymbol{r})\big ). \end{aligned}$ (25)

We can then apply the same quasi-stationary ansatz as in Eq. (2), along with the boundary condition χ′(0) = 0, and proceed accordingly with the solution to the eigenvalue problem to obtain the modified soliton solution; that is, we solve for the unique, constantly decreasing solution with χ(r → ∞) = 0. We note that ρ_b is taken to be static in time and spherically symmetric. Nonetheless, the validity of the spherical approximation was corroborated for DM-dominated galaxies in Bar et al. (2019b), where axially-symmetric baryonic distributions were implemented in the SP equations. The spherically averaged baryonic density is computed from the baryonic mass profile

$\begin{matrix} M_{b} (r) = \frac{r}{G} (V_{*}^{2} (r) + V_{gas}^{2} (r)), \end{matrix}$ $\begin{aligned} M_{\rm b}(r) = \frac{r}{G}\big (V_*^2(r) + V_{\text{gas}}^2(r) \big ), \end{aligned}$ (26)

using Eq. (21) so that

$\begin{matrix} ρ_{b} (r) = \frac{1}{4 π r^{2}} \frac{d M_{b}}{d r}, \end{matrix}$ $\begin{aligned} \rho _{\rm b}(r) = \frac{1}{4\pi r^2} \frac{\text{ d} M_{\rm b}}{\text{ d} r}, \end{aligned}$ (27)

where we compute the derivative analytically.

Our initial approach to assess the self-consistency of the model with baryons is to solve for the modified soliton for each galaxy fixing the axion mass and central density to the corresponding central values in Table 2. This leads to solitons more compact than the original ones, in agreement with Bar et al. (2018) and the results found in simulations (Mocz et al. 2019; Veltmaat et al. 2020). However, the differences between these solitons and their DM-only counterparts are less pronounced than in these simulations mainly because the relative proportion of baryons compared to DM is significantly lower in the LTs galaxies. This is illustrated in Fig. 7, where we show the effect of the gas and stars in the density distribution of the soliton for DDO 154 and NGC 2366.

Fig. 7.

Soliton density profiles for DDO 154 (left) and NGC 2366 (right) including the baryonic contributions to the gravitational potential. The red solid line denotes the original soliton matching the median values of the FDM mass and central density (with the latter being fixed in all cases). The blue dot-dashed line corresponds to the same FDM mass with the inclusion of baryons, while the green dashed line to the mass needed to replicate the same density profile as the original soliton. Lastly, the orange dotted line corresponds to a soliton without baryons fixed at the value of the global optimal mass from Eq. (23).

We confirm the findings of Bar et al. (2018) that solitons formed including the baryonic effects closely follow the same universal distribution as the DM-only soliton, for instance, the one approximated by Eq. (9). This does not imply necessarily that the two solitons (DM-only and with baryons) are the same but rather that their difference is captured by a modification of the relation between ρ_c and r_c for a given m_a (or of any other scaling relations such as r_c and M_c); or, conversely, that ρ_c and r_c parametrize the same density profile regardless of whether baryons are included or not but where each case corresponds to a different value of m_a. In particular, if baryons are included for a given empirical density profile, a lower m_a should be obtained to compensate for the deepening of the gravitational potential.

Therefore, we can implement the inclusion of baryons in our analysis of the rotation curves of LTs by including a systematic negative shift in the FDM mass needed to replicate the original profile favored by the fits. In Table 3 we list the corrections to the value of m_a (with respect to the results listed in Table 2) that is required to fit the rotation curves including the effects of the baryons. We find that the correction needed is typically within the original 68% CL errors, with the only exception of DDO 154 and NGC 2366 (both of which appear in Fig. 7) and NGC 6822 in case of the rogues. Based on these results, we conclude that the effect of including the background baryonic potential introduces a systematic drop of ≲15% in the axion masses which should be included as a correction in the determination of m_a in Eq. (23). We note that this effect does not explain the empirical dependence of m_a with M_⋆ that we found in our sample and does not help understanding the discrepancies of the data with the soliton scaling relations discussed in Sect. 4.2. Furthermore, this effect of the baryons only exacerbates the tension found in the previous Sect. 4.3 in relation to the excessive suppression of the HMF at small scales.

Table 3.

Systematic correction in FDM mass central value as a result of introducing the modified soliton due to a background baryonic potential.

In connection to the baryonic effects, we also investigate the inner-slope parameter of the DM distribution of the LTs galaxies as a function of their mass. This quantity characterizes the flatness of the central DM distribution and it is a figure of merit for many studies of baryonic feedback within ΛCDM. We follow Di Cintio et al. (2014), Tollet et al. (2016) and define this parameter as the logarithmic slope of the density ρ(r)∝r^α fitted to the data between 1 − 2% of the virial radius. In Fig. 8 we show the values of α for the LTs sample as a function of M_⋆ and compared to the prediction from simulations in ΛCDM (Di Cintio et al. 2014; Tollet et al. 2016). The M_⋆ range covered by the LTs galaxies overlaps with the region where baryonic feedback starts becoming relevant, M_⋆ ≳ 10⁷ M_⊙. The most massive galaxies in the sample present cores which are consistent with those expected from baryonic feedback while the less massive ones, with M_⋆ ≈ 10⁷ M_⊙, present cores generally more extended (or flatter) than predicted by these processes. We note, however, that we are comparing results from two different DM models.

Fig. 8.

Inner slope vs stellar mass for simulations running with the optimal FDM mass of 1.9 × 10⁻²³ eV. In green appears the prediction of baryonic feedback simulations from Di Cintio et al. (2014), Tollet et al. (2016). We use a 1σ scatter of 0.3, matching closely with the results of Di Cintio et al. (2014), Tollet et al. (2016), while also adding an additional 2σ scatter band. The gray band represents a conservative scatter derived analytically for NFW profiles.

5. Comparison with previous works

Various works have recently investigated the solitonic profiles using galactic rotation curves. These papers use data from the 175 galaxies compiled by the SPARC database (Bar et al. 2018, 2022; Robles et al. 2019; Chan & Fai Yeung 2021; Street et al. 2022; Khelashvili et al. 2023; Dave & Goswami 2023), while Bernal et al. (2018) also studied a set of low-surface brightness (LSB) galaxies. In terms of methodology, our analysis presents some similarities with some of these works:

In their work, Bernal et al. (2018) focuses on χ² fits to the selected 18 LSB galaxies and 6 from the SPARC data set. A large scatter of preferred axion masses is obtained although a good global fit to the LSB set is obtained for m_a ≈ 0.5 × 10⁻²³ eV. The soliton or soliton-halo scaling relations are not investigated nor the cosmological implications of the model for this mass.

In their work, Bar et al. (2018) focuses on the physical interpretation and phenomenological implications of the core-halo relation. In particular, the appearance of a “double-bump” structure in the rotation curves is predicted and then used to set an exclusion of FDM in the range of m_a ∼ 10⁻²² − 10⁻²¹ eV range. This structure is not a generic feature appearing in our analyses due to the broader range of core-halo relations that we use, as described in detail in Appendix C. In their work, Bar et al. (2022) derived, for each m_a, upper bounds on the mass of the soliton using a conservative method and the 175 galaxies in SPARC. However, this reference does not perform a full-fledged fit to the data and uses exclusively the core-halo relation as a exclusion criterion for FDM in the range m_a ∼ 10⁻²³ − 10⁻²⁰ eV. In their work, Khelashvili et al. (2023) performs a Bayesian analysis of the 175 galaxies in SPARC and studies the statistical comparison between FDM and different CDM models. This yields that FDM is preferred by the data although no single value of m_a gives a good fit to all galaxies. Moreover, the core-halo and soliton scaling relations are tested.

FDM has been studied in other different cosmological and astrophysical contexts and the derived limits on the axion mass have been recently compiled in Ferreira (2021). In the following paragraphs we include some of the cosmological bounds available in the literature.

The small-scale suppression given by FDM would impact the angular scale of the CMB acoustic peaks and anisotropies, as proposed in Hlozek et al. (2015, 2018). Using data from Planck 2015 and WiggleZ, these studies derive the lower bound m_a ≳ 10⁻²⁴ eV.

The Lyman-α forest is a powerful tracer of the linear matter power spectrum for sub-Mpc scales, where the FDM suppression arises (Niemeyer 2020). With high-resolution observations of the spectra on the Lyman-α forest, Nori et al. (2019), Kobayashi et al. (2017), Armengaud et al. (2017), Rogers & Peiris (2021) have determined lower bounds in the range m_a > (2.0 − 20)×10⁻²¹ eV.

At high redshift, the 21 cm absorption signal of neutral hydrogen at z ∼ 15 − 20, measured by EDGES (Bowman et al. 2018) is also sensitive to a suppression of the power spectrum at small scales (Lidz & Hui 2018) as it would delay galaxy formation (Schneider 2018). This would lead to the lower bound m_a > 8 × 10⁻²¹ eV (Schneider 2018).

The shear correlation of galaxies can be used to study the impact of the FDM scenario. With the Dark Energy Survey Year 1 (DES-Y1) data and the Planck cosmic microwave background anisotropies information, Dentler et al. (2022) has reported a search of the effects caused by FDM in the cosmic shear. This sets m_a > 10⁻²³ eV at 95% CL.

The abundances of lensed ultra-faint galaxies at high redshift provides strong constraints on the suppression of the HMF of FDM. With the measurements from Hubble Frontier Field (HFF) of the number density of z ≈ 6 galaxies, Menci et al. (2017) found the lower limit m_a ≥ 8 × 10⁻²² at 3σ CL.

We now include results from astrophysical searches and bounds. For instance, in the regime of zero self-coupling, the non-detection of black-hole superradiance leads to two distinctive ranges of exclusion of axion masses (Arvanitaki & Dubovsky 2011; Baryakhtar et al. 2017; Brito et al. 2017): Supermassive black holes exclude 7 × 10⁻²⁰ eV < m_a < 10⁻¹⁶ eV at the 95% CL, while stellar mass black holes yield an exclusion region 7 × 10⁻¹⁴ eV < m_a < 2 × 10⁻¹¹ eV at the 95% CL. With superradiance analyses, the Event Horizon Telescope measurement of the spin and mass of M87^* also disfavors the range 2.9 × 10⁻²¹ eV < m_a < 4.6 × 10⁻²¹ eV (Davoudiasl & Denton 2019).

The FDM’s wave behavior would heat up stellar objects in galaxies (Bar-Or et al. 2019). In case of the Milky Way this effect can be search for through the velocity dispersion of stars in the solar system vicinity, leading to the lower bound m_a > 6 × 10⁻²³ eV (Church et al. 2019). A recent analysis (Chiang et al. 2023) for the Milky Way disk heating leads to the bound m_a > 0.4 × 10⁻²² eV. By accounting for different sources of uncertainties, Chiang et al. (2023) finds that the range of m_a ∼ 0.5 − 0.7 × 10⁻²² eV is favored by the observed Milky-Way disc kinematics.

FDM affects dynamical friction and the effects can be searched for with the analysis of orbits of globular clusters in galaxies (Hui et al. 2017; Bar et al. 2021). In fact, FDM can help solving the so-called Fornax timing puzzle Tremaine (1976) as long as m_a > 10⁻²¹ eV (Lancaster et al. 2020).

Stellar streams are excellent tracers of the gravitational potential (Carballo-Bello et al. 2014) and can be used to test predictions of the sub-halo mass function in different DM models. The fluctuations in the stellar streams lead to imposing the lower limit m_a > 2 × 10⁻²¹ eV (Schutz 2020).

Orbital motions of stars in the galactic center can also probe the matter distribution through its gravitational potential. Using mock catalogs of astrometric and spectroscopic observations of the star S2 orbiting near Sagittarius A^*, Della Monica & de Martino (2023) obtained an upper limit of m_a < 10⁻²¹ eV at a 95% CL.

The suppression of the cosmic growth in FDM reduces the number of low-mass halos, which impinges on the dwarf-spheroidal galaxies’ (dSph) typical distance scales. From the observed population of satellites in the Milky Way the mass of the FDM is required to be m_a > 2.9 × 10⁻²¹ eV (Nadler et al. 2019).

By applying a Jeans analysis for Milky Way-satellite dSphs, it is possible to study the density profile predicted by FDM. By considering a soliton core and kinematic data of the classical dSphs, Chen et al. (2017) obtained the value for $m_{a} = 1 . 79_{- 0.33}^{+ 0.35} \times 10^{- 22}$ $m_{\mathrm{a}}=1.79^{+0.35}_{-0.33}\times10^{-22}$ eV (at 2σ). Comparing to simulated data on the line-of-sight velocity dispersion (σ_LOS) of Fornax and Sculptor dSphs, González-Morales et al. (2017) found a tighter upper bound m_a < 4 × 10⁻²³ eV at 97.5% CL.

The existence of the sub-halo to host Eridanus II low-surface brightness dwarf galaxy and the survival of its star cluster can also constrain the FDM mass (Bar-Or et al. 2019; Brandt 2016). As shown in Marsh & Niemeyer (2019), the sub-halo mass function in the Milky Way and the existence of Eridanus II implies m_a ≳ 8 × 10⁻²² eV. Gravitational heating can also be applied in the case of Eridanus II and its star cluster, which will be destroyed by gravitational heating from the fluctuations of soliton peak density for m_a ≲ 10⁻¹⁹ eV with the possible exception of the range 10⁻²¹ eV ≤ m_a ≤ 10⁻¹⁹ eV which is still allowed by observations (Marsh & Niemeyer 2019). We note, however, that this lower limit was recently disputed by Chiang et al. (2021).

Sizes and stellar kinematics of ultra-faint dwarf (UFD) galaxies can be used to estimate the effects of wave interference produced by FDM. In particular, in order for the FDM’s gravitational heating to be weak enough to explain the observed velocity dispersions of 2.5 − 3 km s⁻¹ and sizes of 24 − 40 pc in Segue 1 and 2 (Simon et al. 2011; Kirby et al. 2013), Dalal & Kravtsov (2022) derived a lower bound of m_a > 3 × 10⁻¹⁹ eV at 99% CL. On the other hand, a Jeans analysis of 18 galactic UFD galaxies finds a preference for higher values of FDM mass in Hayashi et al. (2021), being the strongest one for Segue 1 with $m_{a} = 1 . 1_{- 0.7}^{+ 8.3} \times 10^{- 19}$ $m_{\mathrm{a}}= 1.1^{+8.3}_{-0.7} \times 10^{-19}$ eV, which is consistent with the results obtained directly from the MUSE-Faint survey (Zoutendijk et al. 2021). Finally, Calabrese & Spergel (2016) found that FDM particle with a mass of m_a ∼ 3.7 − 5.6 × 10⁻²² eV is in agreement with the data of half-light mass of Draco II and Triangulum II; although Safarzadeh & Spergel (2020) later found an upper limit m_a ≲ 10⁻²¹ eV using a kinematic analysis of the Milky Way’s satellites.

In Fig. 9 we show a selection of the existing constraints on the axion mass along with the ones derived in this work from the rotation curves of the isolated dwarf irregular galaxies in the LTs sample. In particular, the 2σ interval stemming directly from the average of the fits to the rotation curves, Eq. (23), and the lower bounds derived from the comparison of the HMF of FDM with dwarf-galaxy counts in the LGV. Nonetheless, before closing this section, we would like to point out that many of the analyses leading to these constraints depend on astrophysical modeling assumptions and are, thus, subject to significant uncertainties. We refer, for instance, to Ferreira (2021) and, in particular, to Sect. 4.3 of Chiang et al. (2023) for a recent critical reappraisal of some of these constraints.

Fig. 9.

Bounds from cosmology and astrophysics on the axion mass. The constraint from HMF in the local group with the halos of nearby dwarf irregulars (dIrrs) from the LTs catalog is displayed with dark turquoise bars. On the dIrrs bar, we also show the optimal mass $m_{a} = 1 . 9_{- 0.4}^{+ 0.5} \times 10^{- 23}$ $m_{\mathrm{a}}=1.9^{+0.5}_{-0.4}\times10^{-23}$ eV (uncertainties at 2σ CL). Our constraints are compared to other cosmological and astrophysical limits (see main text for details).

6. Conclusions

In this work we have investigated and tested the predictions of FDM at galactic scales using a homogeneous and robust sample of high-resolution rotation curves from the LITTLE THINGS in 3D (or LTs for short in our paper) catalog. This comprises a collection of isolated, dark matter-dominated dwarf-irregular galaxies that provides an optimal benchmark for cosmological studies. The methodological basis of our study is a statistical framework based on the χ² analysis of the rotation curves using a soliton+NFW density profile as a theoretical model. This depends on four parameters: two for the soliton that we chose to be the axion mass, m_a, and the core mass, M_c, and two for the NFW halo that we choose to be the concentration parameter, c, and the halo mass M_h. We fit the data using current MCMC techniques and a rather loose set of priors, except for the core-halo relation linking M_c and M_h for which we used a broad set of predictions stemming from the FDM simulations compiled in Chan et al. (2022). From the results of the fits, we were able to perform various diagnostics on the predictions of FDM that allow us to draw the main conclusions of our work:

The soliton+NFW model provides an excellent fit to the rotation curves of the LTs sample with the inferred axion masses clustering around a relatively narrow range of values m_a ∼ (1 − 5)×10⁻²³ eV. Some of the galaxies lead to very stringent constraints on m_a (see Fig. 4) so that combining them yields a relatively poor fit. If we conservatively attribute this to possible systematic effects or presence of outliers we obtain $m_{a} = 1 . 90_{- 0.21}^{+ 0.24} \times 10^{- 23}$ $m_{\mathrm{a}}=1.90^{+0.24}_{-0.21}\times 10^{-23}$ eV.
However, we find that the individual determinations of m_a are not scattered randomly around the average, but follow instead a mild correlation with the stellar mass of the galaxy with a significance of ∼3.3σ. Moreover, displaying the results of our fits in the r_c − ρ_c and r_c − M_c planes, we find scaling relations in the data that are at odds with the predictions of FDM. In case of the r_c − M_c relation, the data show a scaling that is almost orthogonal to the one predicted by FDM with a significance that (at face value) is ≳5σ. From this, we conclude that the cores we find in the LTs catalog do not have the properties of the standard FDM solitons.
In parallel, we investigated various cosmological implications related to the FDM halos for an axion mass m_a ≈ 2 × 10⁻²³ eV. The most important result is that the mere existence of the very same galaxies we are studying would be excluded by the strong suppression in the linear power spectrum that is predicted for FDM with this mass. Combining the HMF of FDM obtained in simulations with the one of CDM obtained in simulations of the Local Group, we derive a conservative bound of m_a ≳ 4.3 × 10⁻²³ eV. This represents a tension with the average axion mass determined from the rotation curves that (at face value) again has a significance of ≳5σ. Including all galaxies from a recent census yields a much more severe lower bound m_a ≳ 5.5 × 10⁻²² eV.
Finally, we investigated the effects of baryons in our analysis. In particular, we estimated the contribution of the mass distribution of stars and gas to the structure of the soliton and found that this can be parametrized by a galaxy-dependent negative correction to the axion mass. However, these corrections are typically too small to alleviate any of the problems of FDM with data that we found. We also briefly discuss the possible role of baryonic feedback and pointed out that our galaxy sample is in a region of stellar masses that overlaps with the one where these effects can be significant in CDM.

In summary, our analysis poses a serious challenge to FDM as a contending hypothesis for solving the core-cusp problem or any of the other small-scale issues related to ΛCDM. The tension that we have found between the favored FDM mass and the abundance of halos is yet another example of the problems of the model with cosmology, which is reminiscent of the so-called “catch-22 problem” found in warm DM (Macciò et al. 2012) and renewed recently in the context of FDM (Mocz et al. 2023). This arises when one wishes to replicate the core-like structures in dwarf galaxies by invoking a particle mass that would imply a suppression of small-scale structure that is much too strong to be consistent with observations. Baryonic physics within FDM could help relieving some of these tensions, although it seems unlikely in light of our findings and the results from pioneering FDM simulations with realistic baryonic feedback. An alternative approach would be extending the dark sector with, for instance, self-interactions, as recently explored in Mocz et al. (2023) as a possible mechanism to enhance small-scale structure formation, or by demanding that only a fraction of DM is FDM, as in, for instance, Kobayashi et al. (2017). Future work is expected to shed more light on these topics. On the observational side, a reanalysis of HI data from other isolated dwarf irregular galaxies, using the same methods as in the LTs catalog, would increase the statistical power of cosmological analyses, such as the one presented in this paper.

¹

The lack of a first-principles treatment of baryonic feedback makes the assessment of their precise effects uncertain and somewhat controversial. For instance, simulations from Sawala et al. (2016) and Zhu et al. (2016) do not find baryonic-feedback induced cores for dwarf galaxies while those from Read et al. (2016a) find the cores for all of them (see Bullock & Boylan-Kolchin 2017 for a brief discussion). Note: also the systematic errors in the modeling required for extracting the rotation curves from data have been suggested to introduce considerable uncertainties, explaining the observed diversity (Roper et al. 2023).

²

The theoretical framework presented in this section is obtained specifically from Bar et al. (2018), Guzman & Urena-Lopez (2004), Moroz et al. (1998), Hui et al. (2017), Ruffini & Bonazzola (1969). In our formulae, we implicitly use natural units ℏ = c = 1.

³

Some of the authors that we cite and compare with used somewhat different cosmologies and related conventions. While this is a technical consideration one should be aware of, such modifications have little effect on the results, allowing for adequate comparisons to suit our purposes.

⁴

This was not the scatter originally reported in Dutton & Macciò (2014; 0.11 dex), but we instead decided to use the more conservative one suggested by Diemer & Kravtsov (2015). We also note that the CDM relation used in Dentler et al. (2022) was somewhat different from the one used here. However, we checked that, when applied to Eq. (16), the resulting FDM relation was very similar.

⁵

The correction takes into account that both the velocity dispersion and the velocity rotation are important to trace the gravitational potential. This effect is appreciable when the dispersion velocities are comparable with the maximum of the velocities in the rotation curves, which is the typical behavior in many LTs dwarf galaxies. The correction is equivalent to the asymmetric drift correction for the stars.

⁶

The galaxies DDO 47 and IC 1613 in the LTs sample present, in addition, significant irregularities in the rotation curves correlated with fast expanding HI bubbles, while DDO 53 is almost face on.

⁷

See https://emcee.readthedocs.io/en/stable/tutorials/moves/, where this is explained further using the same combination of moves.

⁸

https://github.com/acastillodm/FuzzyDM

⁹

https://github.com/acastillodm/FuzzyDM

¹⁰

This bound can be applied using galaxies beyond those in our sample (for which we have extracted M_h) because we are conservatively integrating from 10¹¹ M_⊙ down to negligible masses.

¹¹

We did this as such data were not explicitly available in Behroozi et al. (2019) due to their use of a non-parametric algorithm.

¹²

In fact, this difference was also noticed by Bar et al. (2018) and discussed in its Appendix A. There, the difference in the computations of E/M|_halo is absorbed into a redefinition of the halo mass which is not equivalent to the simple reparametrization of the NFW profile mentioned above.

Acknowledgments

A.B.-H., A.C. and J.M.C. would like to thank Diego Blas, Jorge Sánchez Almeida and Ignacio Trujillo Cabrera for introducing us to the small-scale issues of ΛCDM and for continuous guidance and encouragement. In addition, we thank Giuseppina Battaglia, Kfir Blum, Chris Brook and Arianna Di Cintio for guidance at different stages of the project. Finally, we would also like to thank Filippo Fraternali, Jorge García-Farieta, Théo Gayoux, Se-Heon Oh, Marc Huertas-Company, Javier Olivares and Hsi-Yu Schive for useful discussions and/or for sharing data and information with us. A.B.-H. A.C. and J.M.C. acknowledge support from the grant PGC2018-102016-A-I00 and J.M.C. is also funded by the “Ramón y Cajal” program RYC-2016-20672. G.I. acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017.

References

Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617 [Google Scholar]
Armengaud, E., Palanque-Delabrouille, N., Yèche, C., Marsh, D. J. E., & Baur, J. 2017, MNRAS, 471, 4606 [NASA ADS] [CrossRef] [Google Scholar]
Arvanitaki, A., & Dubovsky, S. 2011, Phys. Rev. D, 83, 044026 [NASA ADS] [CrossRef] [Google Scholar]
Bar, N., Blas, D., Blum, K., & Sibiryakov, S. 2018, Phys. Rev. D, 98, 083027 [NASA ADS] [CrossRef] [Google Scholar]
Bar, N., Blum, K., Eby, J., & Sato, R. 2019a, Phys. Rev. D, 99, 103020 [NASA ADS] [CrossRef] [Google Scholar]
Bar, N., Blum, K., Lacroix, T., & Panci, P. 2019b, JCAP, 07, 045 [CrossRef] [Google Scholar]
Bar, N., Blas, D., Blum, K., & Kim, H. 2021, Phys. Rev. D, 104, 043021 [NASA ADS] [CrossRef] [Google Scholar]
Bar, N., Blum, K., & Sun, C. 2022, Phys. Rev. D, 105, 083015 [NASA ADS] [CrossRef] [Google Scholar]
Bar-Or, B., Fouvry, J.-B., & Tremaine, S. 2019, ApJ, 871, 28 [Google Scholar]
Baryakhtar, M., Lasenby, R., & Teo, M. 2017, Phys. Rev. D, 96, 035019 [NASA ADS] [CrossRef] [Google Scholar]
Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143 [NASA ADS] [CrossRef] [Google Scholar]
Bernal, T., Fernández-Hernández, L. M., Matos, T., & Rodríguez-Meza, M. A. 2018, MNRAS, 475, 1447 [NASA ADS] [CrossRef] [Google Scholar]
Bertone, G., & Hooper, D. 2018, Rev. Mod. Phys., 90, 064050 [NASA ADS] [CrossRef] [Google Scholar]
Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93 [Google Scholar]
Bosma, A. 1978, PhD thesis, University of Groningen, The Netherlands [Google Scholar]
Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67 [Google Scholar]
Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40 [NASA ADS] [CrossRef] [Google Scholar]
Brandt, T. D. 2016, ApJ, 824, L31 [NASA ADS] [CrossRef] [Google Scholar]
Brito, R., Ghosh, S., Barausse, E., et al. 2017, Phys. Rev. D, 96, 064050 [NASA ADS] [CrossRef] [Google Scholar]
Brook, C. B., Di Cintio, A., Knebe, A., et al. 2014, ApJ, 784, L14 [NASA ADS] [CrossRef] [Google Scholar]
Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343 [Google Scholar]
Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 539, 517 [CrossRef] [Google Scholar]
Burkert, A. 2015, ApJ, 808, 158 [NASA ADS] [CrossRef] [Google Scholar]
Burkert, A. 2020, ApJ, 904, 161 [NASA ADS] [CrossRef] [Google Scholar]
Calabrese, E., & Spergel, D. N. 2016, MNRAS, 460, 4397 [NASA ADS] [CrossRef] [Google Scholar]
Carballo-Bello, J. A., Sollima, A., Martinez-Delgado, D., et al. 2014, MNRAS, 445, 2971 [NASA ADS] [CrossRef] [Google Scholar]
Chan, M. H., & Fai Yeung, C. 2021, ApJ, 913, 25 [CrossRef] [Google Scholar]
Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981 [NASA ADS] [CrossRef] [Google Scholar]
Chan, J. H. H., Schive, H.-Y., Woo, T.-P., & Chiueh, T. 2018, MNRAS, 478, 2686 [NASA ADS] [CrossRef] [Google Scholar]
Chan, H. Y. J., Ferreira, E. G. M., May, S., Hayashi, K., & Chiba, M. 2022, MNRAS, 511, 943 [NASA ADS] [CrossRef] [Google Scholar]
Chavanis, P.-H. 2011, Phys. Rev. D, 84, 043531 [NASA ADS] [CrossRef] [Google Scholar]
Chavanis, P.-H. 2022, Phys. Rev. D, 106, 043538 [NASA ADS] [CrossRef] [Google Scholar]
Chavanis, P. H., & Delfini, L. 2011, Phys. Rev. D, 84, 043532 [NASA ADS] [CrossRef] [Google Scholar]
Chen, S.-R., Schive, H.-Y., & Chiueh, T. 2017, MNRAS, 468, 1338 [NASA ADS] [CrossRef] [Google Scholar]
Chiang, B. T., Schive, H.-Y., & Chiueh, T. 2021, Phys. Rev. D, 103, 103019 [NASA ADS] [CrossRef] [Google Scholar]
Chiang, B. T., Ostriker, J. P., & Schive, H.-Y. 2023, MNRAS, 518, 4045 [Google Scholar]
Chou, A. S., Soares-Santos, M., Tait, T. M. P., et al. 2022, ArXiv e-prints [arXiv:2211.09978] [Google Scholar]
Church, B. V., Mocz, P., & Ostriker, J. P. 2019, MNRAS, 485, 2861 [NASA ADS] [CrossRef] [Google Scholar]
Correa, C. A., Schaller, M., Ploeckinger, S., et al. 2022, MNRAS, 517, 3045 [NASA ADS] [CrossRef] [Google Scholar]
Dalal, N., & Kravtsov, A. 2022, Phys. Rev. D, 106, 063517 [NASA ADS] [CrossRef] [Google Scholar]
Dave, B., & Goswami, G. 2023, JCAP, 2023, 015 [CrossRef] [Google Scholar]
Davoudiasl, H., & Denton, P. B. 2019, Phys. Rev. Lett., 123, 021102 [Google Scholar]
de Blok, W. J. G., & Bosma, A. 2002, A&A, 385, 816 [CrossRef] [EDP Sciences] [Google Scholar]
Dehghani, R., Salucci, P., & Ghaffarnejad, H. 2020, A&A, 643, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Delgado, V., & Muñoz Mateo, A. 2022, MNRAS, 518, 4064 [CrossRef] [Google Scholar]
Della Monica, R., & de Martino, I. 2023, A&A, 670, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Deng, H., Hertzberg, M. P., Namjoo, M. H., & Masoumi, A. 2018, Phys. Rev. D, 98, 023513 [NASA ADS] [CrossRef] [Google Scholar]
Dentler, M., Marsh, D. J. E., Hložek, R., et al. 2022, MNRAS, 515, 5646 [NASA ADS] [CrossRef] [Google Scholar]
Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415 [Google Scholar]
Di Paolo, C., Salucci, P., & Erkurt, A. 2019, MNRAS, 490, 5451 [NASA ADS] [CrossRef] [Google Scholar]
Diemer, B. 2018, ApJS, 239, 35 [NASA ADS] [CrossRef] [Google Scholar]
Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108 [Google Scholar]
Domcke, V., & Urbano, A. 2015, JCAP, 01, 002 [CrossRef] [Google Scholar]
Donato, F., Gentile, G., Salucci, P., et al. 2009, MNRAS, 397, 1169 [Google Scholar]
Du, X. 2018, Ph.D. Thesis, Georg-August-Universität Göttingen, Germany [Google Scholar]
Du, X., Behrens, C., & Niemeyer, J. C. 2017, MNRAS, 465, 941 [NASA ADS] [CrossRef] [Google Scholar]
Dubinski, J., & Carlberg, R. G. 1991, ApJ, 378, 496 [Google Scholar]
Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
Efstathiou, G. 1992, MNRAS, 256, 43 [Google Scholar]
Eggemeier, B., & Niemeyer, J. C. 2019, Phys. Rev. D, 100, 063528 [NASA ADS] [CrossRef] [Google Scholar]
Ferreira, E. G. M. 2021, A&ARv, 29, 7 [NASA ADS] [CrossRef] [Google Scholar]
Fitts, A., Michael, B.-K., Elbert, O. D., et al. 2017, MNRAS, 471, 3547 [NASA ADS] [CrossRef] [Google Scholar]
Flores, R. A., & Primack, J. R. 1994, ApJ, 427, L1 [NASA ADS] [CrossRef] [Google Scholar]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
Frenk, C. S., White, S. D. M., Davis, M., & Efstathiou, G. 1988, ApJ, 327, 507 [NASA ADS] [CrossRef] [Google Scholar]
González-Morales, A. X., Marsh, D. J. E., Peñarrubia, J., & Ureña López, L. A. 2017, MNRAS, 472, 1346 [CrossRef] [Google Scholar]
Goodman, J. 2000, New Astron., 5, 103 [NASA ADS] [CrossRef] [Google Scholar]
Goodman, J., & Weare, J. 2013, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
Gorbunov, D., Khmelnitsky, A., & Rubakov, V. 2008, JCAP, 10, 041 [CrossRef] [Google Scholar]
Guzman, F. S., & Urena-Lopez, L. A. 2004, Phys. Rev. D, 69, 124033 [NASA ADS] [CrossRef] [Google Scholar]
Guzman, F. S., & Urena-Lopez, L. A. 2006, ApJ, 645, 814 [CrossRef] [Google Scholar]
Hayashi, K., Ferreira, E. G. M., & Chan, H. Y. J. 2021, ApJ, 912, L3 [NASA ADS] [CrossRef] [Google Scholar]
Hlozek, R., Grin, D., Marsh, D. J. E., & Ferreira, P. G. 2015, Phys. Rev. D, 91, 103512 [NASA ADS] [CrossRef] [Google Scholar]
Hlozek, R., Marsh, D. J. E., & Grin, D. 2018, MNRAS, 476, 3063 [CrossRef] [Google Scholar]
Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv e-prints [arXiv:1008.4686] [Google Scholar]
Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [NASA ADS] [CrossRef] [Google Scholar]
Hui, L. 2021, ARA&A, 59, 247 [NASA ADS] [CrossRef] [Google Scholar]
Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541 [NASA ADS] [CrossRef] [Google Scholar]
Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, AJ, 144, 134 [CrossRef] [Google Scholar]
Iorio, G., Fraternali, F., Nipoti, C., et al. 2017, MNRAS, 466, 4159 [NASA ADS] [Google Scholar]
Karukes, E. V., & Salucci, P. 2017, MNRAS, 465, 4703 [Google Scholar]
Khelashvili, M., Rudakovskyi, A., & Hossenfelder, S. 2023, MNRAS, 523, 3393 [NASA ADS] [CrossRef] [Google Scholar]
Khlopov, M. I., Malomed, B. A., & Zeldovich, I. B. 1985, MNRAS, 215, 575 [NASA ADS] [CrossRef] [Google Scholar]
Kirby, E. N., Boylan-Kolchin, M., Cohen, J. G., et al. 2013, ApJ, 770, 16 [NASA ADS] [CrossRef] [Google Scholar]
Kirby, E. N., Bullock, J. S., Boylan-Kolchin, M., Kaplinghat, M., & Cohen, J. G. 2014, MNRAS, 439, 1015 [NASA ADS] [CrossRef] [Google Scholar]
Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82 [Google Scholar]
Kobayashi, T., Murgia, R., De Simone, A., Iršič, V., & Viel, M. 2017, Phys. Rev. D, 96, 123514 [NASA ADS] [CrossRef] [Google Scholar]
Kormendy, J., & Freeman, K. C. 2016, ApJ, 817, 84 [NASA ADS] [CrossRef] [Google Scholar]
Krut, A., Argüelles, C. R., Chavanis, P. H., Rueda, J. A., & Ruffini, R. 2023, ApJ, 945, 1 [NASA ADS] [CrossRef] [Google Scholar]
Lancaster, L., Giovanetti, C., Mocz, P., et al. 2020, JCAP, 01, 001 [NASA ADS] [Google Scholar]
Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157 [Google Scholar]
Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2020, ApJS, 247, 31 [Google Scholar]
Li, X., Hui, L., & Yavetz, T. D. 2021, Phys. Rev. D, 103, 023508 [NASA ADS] [CrossRef] [Google Scholar]
Lidz, A., & Hui, L. 2018, Phys. Rev. D, 98, 023011 [NASA ADS] [CrossRef] [Google Scholar]
Macciò, A. V., Paduroiu, S., Anderhalden, D., Schneider, A., & Moore, B. 2012, MNRAS, 424, 1105 [CrossRef] [Google Scholar]
Marsh, D. J. E. 2016, Phys. Rept., 643, 1 [NASA ADS] [CrossRef] [Google Scholar]
Marsh, D. J. E., & Niemeyer, J. C. 2019, Phys. Rev. Lett., 123, 051103 [CrossRef] [Google Scholar]
Matos, T., & Urena-Lopez, L. A. 2001, Phys. Rev. D, 63, 063506 [NASA ADS] [CrossRef] [Google Scholar]
May, S., & Springel, V. 2021, MNRAS, 506, 2603 [NASA ADS] [CrossRef] [Google Scholar]
May, S., & Springel, V. 2023, MNRAS, 524, 4256 [NASA ADS] [CrossRef] [Google Scholar]
McConnachie, A. W. 2012, AJ, 144, 4 [Google Scholar]
Menci, N., Merle, A., Totzauer, M., et al. 2017, ApJ, 836, 61 [NASA ADS] [CrossRef] [Google Scholar]
Mina, M., Mota, D. F., & Winther, H. A. 2022, A&A, 662, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mocz, P., Lancaster, L., Fialkov, A., Becerra, F., & Chavanis, P.-H. 2018, Phys. Rev. D, 97, 083519 [NASA ADS] [CrossRef] [Google Scholar]
Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2019, Phys. Rev. Lett., 123, 141301 [CrossRef] [Google Scholar]
Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2020, MNRAS, 494, 2027 [NASA ADS] [CrossRef] [Google Scholar]
Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2023, MNRAS, 521, 2608 [NASA ADS] [CrossRef] [Google Scholar]
Moore, B. 1994, Nature, 370, 629 [NASA ADS] [CrossRef] [Google Scholar]
Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19 [Google Scholar]
Moroz, I. M., Penrose, R., & Tod, P. 1998, Class. Quant. Grav., 15, 2733 [NASA ADS] [CrossRef] [Google Scholar]
Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903 [Google Scholar]
Nadler, E. O., Gluscevic, V., Boddy, K. K., & Wechsler, R. H. 2019, ApJ, 878, L32, [Erratum: ApJ 897, L46 (2020), Erratum: ApJ 897, L46 (2020)] [NASA ADS] [CrossRef] [Google Scholar]
Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996a, MNRAS, 283, L72 [NASA ADS] [Google Scholar]
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996b, ApJ, 462, 563 [Google Scholar]
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
Niemeyer, J. C. 2020, Prog. Part. Nucl. Phys., 113, 103787 [NASA ADS] [CrossRef] [Google Scholar]
Nori, M., & Baldi, M. 2021, MNRAS, 501, 1539 [Google Scholar]
Nori, M., Murgia, R., Iršič, V., Baldi, M., & Viel, M. 2019, MNRAS, 482, 3227 [NASA ADS] [CrossRef] [Google Scholar]
Oh, S.-H., de Blok, W. J. G., Brinks, E., Walter, F., & Kennicutt, R. C., Jr 2011, AJ, 141, 193 [NASA ADS] [CrossRef] [Google Scholar]
Oh, S.-H., Hunter, D. A., Brinks, E., et al. 2015, AJ, 149, 180 [CrossRef] [Google Scholar]
Peebles, P. J. E. 2022, Ann. Phys., 447, 169159 [NASA ADS] [CrossRef] [Google Scholar]
Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [NASA ADS] [CrossRef] [Google Scholar]
Putman, M. E., Zheng, Y., Price-Whelan, A. M., et al. 2021, ApJ, 913, 53 [NASA ADS] [CrossRef] [Google Scholar]
Read, J. I., Agertz, O., & Collins, M. L. M. 2016a, MNRAS, 459, 2573 [NASA ADS] [CrossRef] [Google Scholar]
Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2016b, MNRAS, 462, 3628 [NASA ADS] [CrossRef] [Google Scholar]
Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2017, MNRAS, 467, 2019 [NASA ADS] [Google Scholar]
Ren, T., Kwa, A., Kaplinghat, M., & Yu, H.-B. 2019, Phys. Rev. X, 9, 031020 [NASA ADS] [Google Scholar]
Robles, V. H., Bullock, J. S., & Boylan-Kolchin, M. 2019, MNRAS, 483, 289 [NASA ADS] [CrossRef] [Google Scholar]
Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81 [NASA ADS] [CrossRef] [Google Scholar]
Rodrigues, D. C., del Popolo, A., Marra, V., & de Oliveira, P. L. C. 2017, MNRAS, 470, 2410 [NASA ADS] [CrossRef] [Google Scholar]
Rogers, K. K., & Peiris, H. V. 2021, Phys. Rev. Lett., 126, 071302 [NASA ADS] [CrossRef] [Google Scholar]
Roper, F. A., Oman, K. A., Frenk, C. S., et al. 2023, MNRAS, 521, 1316 [NASA ADS] [CrossRef] [Google Scholar]
Rubakov, V. A., & Gorbunov, D. S. 2017, Introduction to the Theory of the Early Universe: Hot big bang theory (Singapore: World Scientific) [Google Scholar]
Rubin, V. C., Ford, Jr., W. K., & Thonnard, N. 1980, ApJ, 238, 471 [NASA ADS] [CrossRef] [Google Scholar]
Ruffini, R., & Bonazzola, S. 1969, Phys. Rev., 187, 1767 [NASA ADS] [CrossRef] [Google Scholar]
Safarzadeh, M., & Spergel, D. N. 2020, ApJ, 893, 21 [NASA ADS] [CrossRef] [Google Scholar]
Sales, L. V., Wetzel, A., & Fattahi, A. 2022, Nat. Astron., 6, 897 [NASA ADS] [CrossRef] [Google Scholar]
Sánchez Almeida, J., Trujillo, I., & Plastino, A. R. 2020, A&A, 642, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931 [Google Scholar]
Schive, H.-Y., Chiueh, T., & Broadhurst, T. 2014a, Nat. Phys., 10, 496 [NASA ADS] [CrossRef] [Google Scholar]
Schive, H.-Y., Liao, M.-H., Woo, T.-P., et al. 2014b, Phys. Rev. Lett., 113, 261302 [NASA ADS] [CrossRef] [Google Scholar]
Schive, H.-Y., Chiueh, T., Broadhurst, T., & Huang, K.-W. 2016, ApJ, 818, 89 [NASA ADS] [CrossRef] [Google Scholar]
Schneider, A. 2018, Phys. Rev. D, 98, 063021 [NASA ADS] [CrossRef] [Google Scholar]
Schutz, K. 2020, Phys. Rev. D, 101, 123026 [NASA ADS] [CrossRef] [Google Scholar]
Schwabe, B., Niemeyer, J. C., & Engels, J. F. 2016, Phys. Rev. D, 94, 043513 [Google Scholar]
Seidel, E., & Suen, W.-M. 1994, Phys. Rev. Lett., 72, 2516 [NASA ADS] [CrossRef] [Google Scholar]
Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [Google Scholar]
Silk, J., Moore, B., & Diemand, J. 2010, in Particle Dark Matter: Observations, Models and Searches, ed. G. Bertone (Cambridge: Cambridge Univ. Press) [Google Scholar]
Simon, J. D., Geha, M., Minor, Q. E., et al. 2011, ApJ, 733, 46 [NASA ADS] [CrossRef] [Google Scholar]
Sin, S.-J. 1994, Phys. Rev. D, 50, 3650 [NASA ADS] [CrossRef] [Google Scholar]
Sommer-Larsen, J., & Dolgov, A. 2001, ApJ, 551, 608 [NASA ADS] [CrossRef] [Google Scholar]
Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [NASA ADS] [CrossRef] [Google Scholar]
Street, L., Gnedin, N. Y., & Wijewardhana, L. C. R. 2022, Phys. Rev. D, 106, 043007 [NASA ADS] [CrossRef] [Google Scholar]
Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542 [CrossRef] [Google Scholar]
Tremaine, S. D. 1976, ApJ, 203, 345 [NASA ADS] [CrossRef] [Google Scholar]
Tremaine, S., & Gunn, J. E. 1979, Phys. Rev. Lett., 42, 407 [NASA ADS] [CrossRef] [Google Scholar]
Trimble, V. 1987, ARA&A, 25, 425 [NASA ADS] [CrossRef] [Google Scholar]
Tulin, S., Yu, H.-B., & Zurek, K. M. 2013, Phys. Rev. D, 87, 115007 [NASA ADS] [CrossRef] [Google Scholar]
Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189 [Google Scholar]
Veltmaat, J., Niemeyer, J. C., & Schwabe, B. 2018, Phys. Rev. D, 98, 043509 [NASA ADS] [CrossRef] [Google Scholar]
Veltmaat, J., Schwabe, B., & Niemeyer, J. C. 2020, Phys. Rev. D, 101, 083518 [NASA ADS] [CrossRef] [Google Scholar]
Vicens, J., Salvado, J., & Miralda-Escudé, J. 2018, ArXiv e-prints [arXiv:1802.10513] [Google Scholar]
Vogelsberger, M., Zavala, J., & Loeb, A. 2012, MNRAS, 423, 3740 [NASA ADS] [CrossRef] [Google Scholar]
Workman, R. L., Burkert, V. D., Crede, V., et al. 2022, PTEP, 2022, 083C01 [Google Scholar]
Zhu, Q., Marinacci, F., Maji, M., et al. 2016, MNRAS, 458, 1559 [Google Scholar]
Zoutendijk, S. L., Brinchmann, J., Bouché, N. F., et al. 2021, A&A, 651, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: Further details on the fits

In this appendix, we provide additional information about the fits and MCMCs performed in this paper. Various plots describing the relevant information are available in the https://github.com/acastillodm/FuzzyDM GitHub repository. There are three main directories. Variable mass refers to the analysis discussed in Sect. 4.1, where the axion mass m_a is fitted independently for each galaxy. The Benchmark Analysis contains the information regarding the fits where the axion mass in Eq. (23) was used for all the galaxies, as discussed in Sect. 4.2. Finally, NFW contains the details regarding the fits to a NFW profile discussed in App. A.3. The information in each class includes rotation curves, the corner plots and 1D posterior distributions of the parameters, and the MCMC auto-correlation times. In case of the fits to FDM we also include the density profiles for the DM component.

A.1. Results for the FDM profile

In Fig. A.1, we show the corner plots and the 1D distributions of the marginalized posteriors of the parameters for the fits to WLM and DDO 154, which are the two galaxies used as a reference in the main text and in Figs. 2 and 3. The posteriors in the m_a − M_c plane have elliptical shapes (with a negative correlation for both galaxies) that allow us to precisely determine both parameters, with relative uncertainties of ≲20% (see Table 2). Consequently, other derived core properties such as its radius, r_c, or central density, ρ_c, are also precisely determined from the fits. On the other hand, the halo parameters (concentration factor and virial mass) are not well constrained due to the lack of data at large radii and the freedom that gives the continuity condition over transition radius, r_t. However, r_t itself is determined with good precision, with relative errors also in the ballpark of ≲20%. The rotation curves corresponding to the fits displayed in Table 2 for the other galaxies in the core LTs catalog, besides WLM and DDO 154 and including NGC 6822, are shown in Fig. A.2.

Fig. A.1.

Corner and 1D posterior marginalized distributions for the fitting parameters θ_FDM and M_⋆ in the FDM plus baryons model for the two benchmark galaxies: WLM (top panel) and DDO 154 (bottom panel). Posteriors in the 1D plots show the 16%, 50%, and 84% CL ranges. The priors are shown in Table 1 and the blue lines and dots indicate the maximum posterior values of the parameters from the MCMC samples. Contours in the 2D plots display the ∼39.3%, 67.5%, and 86.4% CL regions, corresponding to the 1, 1.5, and 2σ regions of a 2D normal distribution. Cases where an additional inner contour is present denote an ∼11.8% CL or 0.5σ region.

Fig. A.2.

Rotation curves in FDM for the core galaxies in the LTs catalog including NGC 6822 and excluding WLM and DDO 154, shown in Fig. 2. We plot the DM, gas, and stellar components compared to data. The red dashed lines are the median of the posterior distribution for the total velocity and the colored band is the corresponding uncertainty at 68% and 95% CL. High-resolution figures can be found in the https://github.com/acastillodm/FuzzyDM GitHub repository.

A.2. Inclination rogues

We describe in detail the method used to fit the mass models to the galaxies with an inclination of i < 40° (inclination rogues). The inclination i is defined as the angle between the line-of-sight and the vector normal to the plane of the galaxy (Read et al. 2016b, 2017). Hence, a galaxy is edge-on when i = 90° and face-on when i = 0°. To account for this, we use the same method to fit the rotation curves (Sect. 3) but adding now the inclination as a new parameter. The dependence on i of the total rotation velocity V_R (and its uncertainty) is:

$\begin{matrix} V_{R}^{'} = V_{R} \frac{sin (i_{0})}{sin (i)}, and δ V_{R}^{'} = δ V_{R} \frac{sin (i_{0})}{sin (i)}, \end{matrix}$ $\begin{aligned} V_{R}^{\prime }=V_{R}\frac{\sin (i_{0})}{\sin ({i})},\;\;\;\; \text{ and} \;\;\;\; \delta V_{R}^{\prime }=\delta V_{R}\frac{\sin (i_{0})}{\sin ({i})}, \end{aligned}$ (A.1)

where i₀ is a “reference” angle. Therefore, for the regime of i < 40° (nearly face-on galaxies), small errors on i can have a strong impact on the reconstruction of the observed rotational velocity. The i₀ values are taken (with uncertainties) from the original extraction of the rotation curves reported in Read et al. (2017). To compute the circular velocity, we add in quadrature the gravitational components and the asymmetric drift correction. Correspondingly, the final errors are obtained from the transformed uncertainties and the errors in the asymmetric drift correction V_a,

$\begin{matrix} δ V_{c} = \frac{\sqrt{V_{R}^{^{'} 2} δ V_{R}^{^{'} 2} + V_{a}^{2} δ V_{a}^{2}}}{V_{c}} . \end{matrix}$ $\begin{aligned} \delta V_{c} = \frac{\sqrt{V_{R}^{^{\prime }2}\delta V_{R}^{^{\prime }2} + V_{\rm a}^{2}\delta V_{\rm a}^{2}}}{V_{c}}. \end{aligned}$ (A.2)

For the fits to the total velocity of Eq. (20) with the mass model in Eq. (19), we assume a normal distribution prior centered around i₀ and error determined by its uncertainty. The remaining parameters follow the standard priors of Table 1, and the fitting algorithm, based on MCMC, uses the same attributes as those described in Sect. 3.2.

An additional correction for the log-likelihood is necessary since now the errors in eq. (A.1) are i-dependent,

$\begin{matrix} ln L (θ_{FDM}, i, M_{⋆}) = ln L^{'} (θ_{FDM}, M_{⋆}) - \frac{1}{2} \sum_{n} ln (δ V_{R_{n}}^{^{'} 2}), \end{matrix}$ $\begin{aligned} \ln \mathcal{L} (\boldsymbol{\theta }_{\rm FDM},i,{M_{\star }})=\ln \mathcal{L} ^{\prime }(\boldsymbol{\theta }_{\rm FDM},{M_{\star }})-\frac{1}{2}\sum _{n}\ln (\delta V_{R_{n}}^{^{\prime }2}), \end{aligned}$ (A.3)

where ℒ′(θ_FDM,M_⋆) is the log-likelihood function in Eq. (22) changing the data points for velocities and their errors using transformation of Eq. (A.1).

The inclination rogues analyzed in our work are DDO 50 and DDO 133, which have stellar masses in the selected range of M_⋆ ≲ 10⁸ M_⊙. In the case of FDM with variable mass, the parameters are reported in the last rows of Table 2. We show the corresponding rotation curves compared to the data in Fig. A.3. We also obtained the inclination angles $i = (40 . 7_{- 2.4}^{+ 2.3}) °$ $i=(40.7^{+2.3}_{-2.4})\circ$ (DDO 50) and $i = (38 . 1_{- 1.7}^{+ 1.6}) °$ $i=(38.1^{+1.6}_{-1.7})\circ$ (DDO 133), which agree within < 5% with the central values of i₀ used in the priors, namely, 37.4° and 36.9° for DDO 50 and DDO 133, respectively.

Fig. A.3.

Rotation curves in FDM for the inclination rogues. We represent the same information as in Fig. A.2.

A.3. Fits with a NFW profile

We present here the results of the fits to the rotation curves of the LTs catalog using the NFW profile in Eq. (10). We follow the same fit procedure as described in Sect. 3.2, where the theoretical circular velocity in Eq. (20) depends now on θ_NFW = (c, M_h) in the DM model instead of θ_FDM. The priors for the halo parameters and M_⋆ are equal to those described in Table 1. As the NFW profile only has three parameters (two for the DM contribution and one for the baryionic model), the MCMCs are faster than in FDM for the convergence tests. In Table A.1, we summarize the results for the core LTs galaxies plus, the two inclination rogues (DDO 50 and DDO 133) and NGC 6822. We confirm that the NFW profile fails in providing a good fit to most of the LTs galaxies, in contrast to the fits to a “cored” profile such as the FDM model.

Table A.1.

Median and 68% CL of the posterior distributions of the fitted parameters M_⋆, M_h, and c in the NFW + baryons model for the galaxies from the LTs catalog, including some of the rogues and NGC 6822.

Moreover, the values of the obtained halo parameters are inconsistent with the predictions of ΛCDM. In particular, in Fig. A.4 we test the population of DM halos and the M_⋆ − M_h relation via abundance matching (Behroozi et al. 2019) and the c − M_h relation (Dutton & Macciò 2014). We find that most of the galaxies in the LTs sample are beyond the 2σ bands in the case of abundance matching while the data cluster in regions below to the Dutton & Macciò relation for the c − M_h plane.

Fig. A.4.

Median values of posterior distributions for the NFW model compared with M_* − M_h (top) and c − M_h (bottom) relations. Uncertainties in the data points describe the 68% (solid) and 95% (dashed) CL ranges. The green bands describe ±1 and 2 σ scatter with σ = 0.3 dex for the top panel and σ = 0.16 dex for the bottom panel. Results for FDM are included for comparison.

A.4. Fitting linear scaling correlations

Here, we summarize the approach we use to quantify the linear correlations observed in our (variable mass) fits throughout Figs. 4 and 5, namely, in the (log₁₀m_a, log₁₀M_⋆), (log₁₀ρ_c, log₁₀r_c) and (log₁₀r_c, log₁₀M_c) planes. To do this, we follow the process described in Sects. 7 and 8 of Hogg et al. (2010).

Taking the case of the axion - stellar mass correlation as an example, we are interested in a fit of the form

$\begin{matrix} {log}_{10} m_{a} = β {log}_{10} M_{⋆} + α^{'}, \end{matrix}$ $\begin{aligned} \log _{10} m_{\rm a} = \beta \log _{10} {M_{\star }} + \alpha ^\prime , \end{aligned}$ (A.4)

to our data from the posterior distributions of such parameters (obtained from our previous fits using MCMCs), taking the median and mean 68% CL errors as the data to be fitted. The idea is, then, to formulate a likelihood function that we can implement with MCMCs in much the same way as before and derive the posterior distributions of α′ and β, with β being the key variable to assess the degree of correlation as quantified by the slope.

While the idea is simple, we need to account for the fact that we are dealing with errors in two dimensions. To this end, we compute residues as the orthogonal distance between the line and each of the observed central values and qw compute the corresponding projections of the errors, under the assumption that these are normal and independent. To outline the specific way these are computed, it is more efficient to parametrize by angle θ, where β ≡ tanθ, and the orthogonal distance between the line and the origin b_⊥ ≡ α′ cos θ, which are the ones we use internally.

The resulting orthogonal residual for a given point (x_i, y_i) of median observed values in the (log₁₀m_a, log₁₀M_⋆) plane is given by:

$\begin{matrix} Δ_{i} = y_{i} cos θ - x_{i} sin θ - b_{⊥}, \end{matrix}$ $\begin{aligned} \Delta _i = y_i \cos \theta - x_i \sin \theta - b_{\perp }, \end{aligned}$ (A.5)

while the resulting projected error Σ_i is given by

$\begin{matrix} Σ_{i}^{2} = δ x_{i}^{2} {sin}^{2} θ + δ y_{i}^{2} {cos}^{2} θ, \end{matrix}$ $\begin{aligned} \Sigma _i^2 = \delta x_i^2 \sin ^2 \theta + \delta y_i^2 \cos ^2 \theta , \end{aligned}$ (A.6)

where δx_i, δy_i denote the mean 68% CL errors from the median.

We also introduce an additional parameter v to account for potential systematic underestimation of errors. The final log-likelihood function is then given by:

$\begin{matrix} ln L (θ, b_{⊥}, v) = - \frac{1}{2} \sum_{i} \frac{Δ_{i}^{2}}{Σ_{i}^{2} + v} - \frac{1}{2} \sum_{i} ln (Σ_{i}^{2} + v) . \end{matrix}$ $\begin{aligned} \ln \mathcal{L} (\theta , b_{\perp }, v) = -\frac{1}{2} \sum _i \frac{\Delta _i^2 }{\Sigma _i^2 + v} -\frac{1}{2} \sum _i \ln \Big ( \Sigma _i^2 + v \Big ). \end{aligned}$ (A.7)

For the purposes of implementing this likelihood function into the MCMC fitting routine, we set the flat priors −π/2 < θ < π/2, −100 ≤ b_⊥ ≤ 100, and 0 ≤ v ≤ 100, which amply covers the set of expected values. Following the discussion in Section 3, we run the 48-walker chains for a suitably long period that allows for reliable sampling of the distribution, which in this case amounts to 100 000 steps with a 50 000-step burn-in period. This is significantly more than necessary for convergence purposes, but allowed us to probe more extremal values of the distribution that are useful to establish the statistical significance of the observed correlations.

During the course of the fits, the values obtained for v were quite small, indicating that the effect of such correction is not very significant. The results of these fits can be found in the Variable Mass/Relations/Scaling Correlations folder of the https://github.com/acastillodm/FuzzyDM GitHub repository.

Appendix B: Mass functions and abundance matching

Here we provide more details on the procedure we use for computing the abundance matching results and the associated mass function constraints from Fig. 6. Given a mass function $\frac{dn}{d ln M}$ $\frac{d n}{ d\ln M}$ we can define the total number of halos within a given mass range M₁ < M < M₂ and volume V as

$\begin{matrix} N = V \int_{ln M_{1}}^{ln M_{2}} d ln M \frac{dn}{d ln M} . \end{matrix}$ $\begin{aligned} N = V \int _{\ln M_1}^{ \ln M_2} d\ln M \frac{d n}{ d \ln M}. \end{aligned}$ (B.1)

For our purposes, we define the LGV consisting of a sphere of radius 1.8 Mpc. While we are not aware of FDM simulations being performed precisely in this environment and scale, we estimate the HMF of FDM by applying the suppression factor in Eq. (18) to the CDM mass function derived from Brook et al. (2014). While in this paper the mass function (as defined above) was not explicitly given, the cumulative mass function for the abundance of halos and subhalos above a given mass in the LGV was shown (Figure 1). We can accurately reproduce these results by introducing a normalization factor to the original Sheth-Tormen HMF for CDM (Sheth & Tormen 1999) to ensure that the total abundances are the same. This is shown in Fig. B.1, where we observe excellent agreement. We emphasize that a significant multiplicative factor is expected, since the LGV environment has a naturally higher abundance of galaxies than standard cosmological environments and we needs to account for the presence of subhalos.

Fig. B.1.

Comparison of accumulated abundances between the fitting function found in Brook et al. (2014) (black, dashed) from CDM simulations in the LGV and the one obtained from the Sheth-Tormen HMF for CDM (Sheth & Tormen 1999) normalized to match the same total abundance as the previous one (solid, red). We also include the original Sheth-Tormen abundance without the normalization.

With the CDM mass function, we can subsequently apply the suppression factor in Eq. (18) from Schive et al. (2016) to obtain the resulting FDM mass function that enables us to compute halo abundances via Eq. (B.1). We note that while this suppression factor was originally obtained from simulations in higher redshift environments, the factor itself does not have a redshift dependence.

Regarding the abundance matching procedure, in addition to the discussion in Section 4.3, we note that to obtain the effective cumulative SMF (as in Eq. (24), from a given M to ∞), we computed the cumulative HMF using the Sheth-Tormen HMF for each M_h and mapped it to the corresponding M_⋆ via the abundance matching relation of Behroozi et al. (2019). ¹¹ We then subsequently performed the abundance matching procedure with the HMF of FDM (without the normalization factor discussed earlier in the appendix). The Sheth-Tormen HMF is computed via the Colossus package in Python developed by Diemer (2018).

Appendix C: Comparison to the soliton-host halo relation

As discussed in Sect. 5, Bar et al. (2018) independently rederived a core-halo relation from Schive et al. (2014b) and applied it to the analysis of rotation curves of the SPARC database. An important outcome of these studies is the prediction of a “double-bump” in the rotation curves, the first of which should become prominent in the inner part of the galaxies for m_a ≳ 10⁻²² eV. In our analysis we implement the core-halo relations spanned by Eq. (15) from Chan et al. (2022), which also includes Eq. (14) (Schive et al. 2014b), and this structure is not generically present in our theoretical rotation curves (see Fig. 3). In order to understand this aspect, we first present a brief summary of the derivation of the core-halo relation in Bar et al. (2018). Schive et al. (2014a) found a relation connecting the core mass of their FDM distributions and the halo specific energy |E_h|/M_h (in natural units),

$\begin{matrix} M_{c} \approx (\frac{| E_{h} |}{M_{h}})^{\frac{1}{2}} \frac{m_{p}^{2}}{m_{a}} . \end{matrix}$ $\begin{aligned} M_{\rm c} \approx \Bigg (\frac{|E_{\rm h}|}{M_{\rm h}}\Bigg )^{\frac{1}{2}} \frac{m_p^2}{m_{\rm a}}. \end{aligned}$ (C.1)

This can be equivalently expressed as (Bar et al. 2018)

$\begin{matrix} \frac{E}{M} |_{soliton} \approx \frac{E}{M} |_{halo}, \end{matrix}$ $\begin{aligned} \frac{E}{M} \bigg |_{\text{soliton}} \approx \frac{E}{M} \bigg |_{\text{halo}}, \end{aligned}$ (C.2)

due to an analytic property of the soliton where one can express the core mass in terms of the soliton specific energy. Assuming that the halo is dominated by a virialized NFW profile, its specific energy can be computed as:

$\begin{matrix} \frac{E}{M} |_{halo} \approx \frac{π}{M_{h}} \int_{0}^{r_{vir}} d x x^{2} ρ_{NFW} (x) ϕ_{NFW} (x), \end{matrix}$ $\begin{aligned} \frac{E}{M} \bigg |_{\text{halo}} \approx \frac{\pi }{M_{\rm h}} \int _0^{r_{\text{vir}}} dx x^2 \rho _{\text{NFW}}(x)\phi _{\text{NFW}}(x), \end{aligned}$ (C.3)

where ϕ_NFW(x) is the NFW potential. In addition, Bar et al. (2018) used a convention to define the virial radius as the point, where the average density equals 200 times the cosmological critical density, with the corresponding convention for the halo mass denoted by M₂₀₀.

Replacing Eq. (C.3) in Eq. (C.1), one obtains the core-halo relation reported by Bar et al. (2018),

$\begin{matrix} M_{c} \approx 5.7 \times 10^{9} (\frac{m_{a}}{10^{- 23} eV})^{- 1} (\frac{H (z)}{H_{0}})^{\frac{1}{3}} (\frac{M_{200}}{10^{12} M_{⊙}})^{\frac{1}{3}} f (c) M_{⊙}, \end{matrix}$ $\begin{aligned} M_{\rm c} \approx 5.7 \times 10^9 \bigg (\frac{m_{\rm a}}{10^{-23} \text{ eV}} \bigg )^{-1} \bigg (\frac{H(z)}{H_0}\bigg )^{\frac{1}{3}} \bigg (\frac{M_{200}}{10^{12} \, {{M}_{\odot }}} \bigg )^{\frac{1}{3}} f(c)\,{{M}_{\odot }}, \end{aligned}$ (C.4)

where

$\begin{matrix} f (c) = 0.54 \sqrt{(\frac{c}{1 + c}) \frac{c - ln (1 + c)}{(ln (1 + c) - \frac{c}{1 + c})^{2}}}, \end{matrix}$ $\begin{aligned} f(c) = 0.54 \sqrt{\bigg (\frac{c}{1+c}\bigg ) \frac{c - \ln (1+c)}{\big (\ln (1+c) - \frac{c}{1+c}\big )^2}}, \end{aligned}$ (C.5)

is a function of c that varies between 0.9 ≲ f(c)≲1.1 for 5 ≤ c ≤ 30.

The core-halo relation in Eq. (C.4) is analogous to the one in Eq. (14) but it is defined using a different convention for the halo mass. It is straightforward to transform the NFW conventions from a critical spherical overdensity definition of 200 to one of 100, which is close to the one used by Schive et al. (2014b), Chan et al. (2022). One uses a transformation of (c₂₀₀, M₂₀₀)→(c₁₀₀, M₁₀₀) that parametrizes the same NFW density profile in Eq. (10). Transforming M₂₀₀ to M_h increases the given halo mass by 20 − 10% for c₂₀₀ in the range 5 − 30, yielding

$\begin{matrix} M_{c} \approx 5.5 \times 10^{9} (\frac{m_{a}}{10^{- 23} eV})^{- 1} (\frac{M_{h}}{10^{12} M_{⊙}})^{\frac{1}{3}} M_{⊙}, \end{matrix}$ $\begin{aligned} M_{\rm c} \approx 5.5 \times 10^9 \bigg (\frac{m_{\rm a}}{10^{-23} \text{ eV}} \bigg )^{-1} \bigg (\frac{M_{\rm h}}{10^{12} \, {{M}_{\odot }}} \bigg )^{\frac{1}{3}}\,{{M}_{\odot }}, \end{aligned}$ (C.6)

where we have taken z = 0 and f(c) = 1 for simplicity. Therefore, given a halo mass M_h, Eq. (C.6) predicts a core mass that is ≈80% heavier than the one predicted by Eq. (14).

The difference between the two relations is ultimately due to the way the specific energy of the halo is computed. In particular, as shown by Schive et al. (2014b), one arrives at the original relation in Eq. (14) by taking the halo specific energy to be that of a sphere of uniform density at the virial radius as

$\begin{matrix} \frac{E}{M} |_{halo} \approx - \frac{3 G M_{h}}{10 r_{vir}} \end{matrix}$ $\begin{aligned} \frac{E}{M} \bigg |_{\text{halo}} \approx - \frac{3 G M_{\rm h}}{10 r_{\text{vir}}} \end{aligned}$ (C.7)

and applying it to the empirical relation in Eq. (C.1). ¹² This implies an important physical difference, since the potential (and implicitly the virialized energy) of a uniform spherical distribution is not equivalent to that of a self-gravitating NFW profile. It appears that such energy distribution, in combination with Eq. (C.1), is the median of the results of the simulations, as seen in Figs. 2 and 4 of Schive et al. (2014b), although Eq. (C.6) is compatible with an upward scatter of the data.

This relation is also contained within the Eq. (15) from Chan et al. (2022) that is used in our analysis. However, this region of core-halo relations, leading to heavy solitons, is not favored by our fits for relatively heavy axion masses, m_a ≳ 10⁻²² eV. In particular, the rotation curve produced by a soliton has a maximum v_>c,max at a radius r_x,max scaling with the core mass as v_c,max ∝ M_c and r_c,max ∝ 1/M_c (Bar et al. 2018). A heavier M_c by a factor of ∼2 for a given halo mass M_h will lead to a factor of ∼2 higher rotation velocity peak at a distance ∼1/2 shorter from the center of the galaxy. This is illustrated in Fig. C.1 where we show our maximum posterior fit to the LTs galaxy DDO 154 for m_a = 10⁻²² eV compared to the predictions that one would obtain imposing Eq. (14) or (C.6) and using the same M_h. The larger core mass implied by the latter suffices to produce the type of bump in the rotation curve that was discussed in Bar et al. (2018).

Fig. C.1.

Comparison of the effect of imposing different core-halo relations in the rotation curve of DDO 154. Theoretical rotation curves are DM-only (i.e., excluding baryonic components) for visual clarity.

All Tables

Table 1.

Priors used to fit the DM and baryon contributions to the rotation curve data from the LTs catalog.

In the text

Table 2.

Median and 68% confidence level (CL) of the posterior distributions of the parameters m_a, M_c, c, M_h, and M_⋆, and the reduced chi-square of the maximum posterior fit to the LTs catalog, including some of the rogues.

In the text

Table 3.

Systematic correction in FDM mass central value as a result of introducing the modified soliton due to a background baryonic potential.

In the text

Table A.1.

Median and 68% CL of the posterior distributions of the fitted parameters M_⋆, M_h, and c in the NFW + baryons model for the galaxies from the LTs catalog, including some of the rogues and NGC 6822.

In the text

All Figures

	Fig. 1. Solutions to Eq. (4) for the ground state and first two excitations, where the scale invariance was used to set χ(0) = 1.
In the text

	Fig. 2. Rotation curves in FDM for two representative galaxies in the LTs catalog, WLM, and DDO 154. We plot the DM, gas, and stellar components obtained from the fits (maximum posterior values) and compared to the data. The red dashed lines are the median of the posterior distribution for the total velocity and the colored band is the corresponding uncertainty at 68% and 95% CL.
In the text

	Fig. 3. Rotation curves in FDM for two representative galaxies in the LTs catalog: WLM and DDO 154. We plot the total circular velocity and reduced chi-squares resulting from the maximum posterior fits to the data when fixing the axion mass at different values. We compare the results also to the fit to a NFW profile.
In the text

Fig. 4.

Axion masses obtained from the posterior distributions of the fits to the rotation curves of the LTs galaxies where we show medians and 68% and 95% CL errors. The 10 galaxies entering our fiducial analysis are represented by circles while the rogues are represented by triangles. The dot-dashed line and bands represent the mean, 68% and 95% CL regions of a weighted average to these determinations (see Eq. (23) and associated discussion). Left: Galaxies are ordered by the central value of the stellar mass reported in Read et al. (2017), which where found to be in good agreement with our posterior distributions. Stellar masses are increasing from left to right. Right: Empirical fit to a power law between m_a and M_⋆ showing evidence for a correlation between the two as extracted from the galaxies in our sample. The red solid line denotes the maximum posterior fit with the 68% and 95% CL regions.

In the text

Fig. 5.

Soliton scaling relations compared to the results of the fits to LTs, where each point is colored according to the value of M_⋆ of the corresponding galaxy. Dot-dashed lines and gray bands represent the mean and the 68% and 95% CL regions for the correlation expected from FDM for the average axion mass in Eq. (23). The solid line and surrounding bands represent the maximum posterior and 68% and 95% CL regions of the fit to a power law of the results extracted from the individual galaxies. Left: Results in the r_c − ρ_c plane compared to Eq. (8) and where we also show the results obtained independently for a different sample of galaxies by using a Burkert profile (Rodrigues et al. 2017; Deng et al. 2018). Right: Results in the r_c − M_c plane compared to Eq. (8).

In the text

Fig. 6.

Halo properties obtained from the posteriors of the FDM fits to the rotation curves of the LTs galaxies for a fixed axion mass m_a = 1.9 × 10⁻²³ eV. Gray-colored points are galaxies for which both median and maximum posterior values of r_t are within the observed range, while gold-colored are soliton-dominated by virtue of not satisfying that condition. Uncertainties in the data points describe the credible intervals at 68% (solid) and 95% (dashed) CL. Uncertainties in the theoretical bands refer to the 1 and 2σ scatter in the predictions. Top-left: M_c − M_h relation compared to the region allowed by Eq. (15) and the prediction in Eq. (14) from Schive et al. (2014b; dotted line). Top-right: c − M_h relation compared to Eq. (16) and to the prediction in CDM, Eq. (17) from Dutton & Macciò (2014) with σ = 0.16 dex (gray dashed line). Bottom-left: Stellar-to-halo mass relation from abundance matching compared to the FDM prediction using Eq. (18), and the CDM predictions from Behroozi et al. (2019; dashed gray line) and Brook et al. (2014; dot-dashed line) with a reference σ = 0.30 dex scatter in both cases. Bottom-right: Expected number of halos with M_h ≲ 10¹¹ M_⊙ in a Local Group volume of radius 1.8 Mpc predicted in FDM as a function of m_a, using Eq. (18) and calibrated with the CDM prediction of the HMF from Brook et al. (2014). This is compared to the HMF from Schive et al. (2016; dashed line). For the CDM prediction, we take 10⁹ ≲ M_h ≲ 10¹¹ M_⊙. The dot-dashed line and gray bands represent the mean and the 68% and 95% CL regions for the axion mass in Eq. (23).

In the text

Fig. 7.

Soliton density profiles for DDO 154 (left) and NGC 2366 (right) including the baryonic contributions to the gravitational potential. The red solid line denotes the original soliton matching the median values of the FDM mass and central density (with the latter being fixed in all cases). The blue dot-dashed line corresponds to the same FDM mass with the inclusion of baryons, while the green dashed line to the mass needed to replicate the same density profile as the original soliton. Lastly, the orange dotted line corresponds to a soliton without baryons fixed at the value of the global optimal mass from Eq. (23).

In the text

Fig. 8.

Inner slope vs stellar mass for simulations running with the optimal FDM mass of 1.9 × 10⁻²³ eV. In green appears the prediction of baryonic feedback simulations from Di Cintio et al. (2014), Tollet et al. (2016). We use a 1σ scatter of 0.3, matching closely with the results of Di Cintio et al. (2014), Tollet et al. (2016), while also adding an additional 2σ scatter band. The gray band represents a conservative scatter derived analytically for NFW profiles.

In the text

Fig. 9.

Bounds from cosmology and astrophysics on the axion mass. The constraint from HMF in the local group with the halos of nearby dwarf irregulars (dIrrs) from the LTs catalog is displayed with dark turquoise bars. On the dIrrs bar, we also show the optimal mass $m_{a} = 1 . 9_{- 0.4}^{+ 0.5} \times 10^{- 23}$ $m_{\mathrm{a}}=1.9^{+0.5}_{-0.4}\times10^{-23}$ eV (uncertainties at 2σ CL). Our constraints are compared to other cosmological and astrophysical limits (see main text for details).

In the text

Fig. A.1.

Corner and 1D posterior marginalized distributions for the fitting parameters θ_FDM and M_⋆ in the FDM plus baryons model for the two benchmark galaxies: WLM (top panel) and DDO 154 (bottom panel). Posteriors in the 1D plots show the 16%, 50%, and 84% CL ranges. The priors are shown in Table 1 and the blue lines and dots indicate the maximum posterior values of the parameters from the MCMC samples. Contours in the 2D plots display the ∼39.3%, 67.5%, and 86.4% CL regions, corresponding to the 1, 1.5, and 2σ regions of a 2D normal distribution. Cases where an additional inner contour is present denote an ∼11.8% CL or 0.5σ region.

In the text

Fig. A.2.

Rotation curves in FDM for the core galaxies in the LTs catalog including NGC 6822 and excluding WLM and DDO 154, shown in Fig. 2. We plot the DM, gas, and stellar components compared to data. The red dashed lines are the median of the posterior distribution for the total velocity and the colored band is the corresponding uncertainty at 68% and 95% CL. High-resolution figures can be found in the https://github.com/acastillodm/FuzzyDM GitHub repository.

In the text

	Fig. A.3. Rotation curves in FDM for the inclination rogues. We represent the same information as in Fig. A.2.
In the text

	Fig. A.4. Median values of posterior distributions for the NFW model compared with M_* − M_h (top) and c − M_h (bottom) relations. Uncertainties in the data points describe the 68% (solid) and 95% (dashed) CL ranges. The green bands describe ±1 and 2 σ scatter with σ = 0.3 dex for the top panel and σ = 0.16 dex for the bottom panel. Results for FDM are included for comparison.
In the text

	Fig. B.1. Comparison of accumulated abundances between the fitting function found in Brook et al. (2014) (black, dashed) from CDM simulations in the LGV and the one obtained from the Sheth-Tormen HMF for CDM (Sheth & Tormen 1999) normalized to match the same total abundance as the previous one (solid, red). We also include the original Sheth-Tormen abundance without the normalization.
In the text

	Fig. C.1. Comparison of the effect of imposing different core-halo relations in the rotation curve of DDO 154. Theoretical rotation curves are DM-only (i.e., excluding baryonic components) for visual clarity.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617 [Google Scholar]

[2] Armengaud, E., Palanque-Delabrouille, N., Yèche, C., Marsh, D. J. E., & Baur, J. 2017, MNRAS, 471, 4606 [NASA ADS] [CrossRef] [Google Scholar]

[3] Arvanitaki, A., & Dubovsky, S. 2011, Phys. Rev. D, 83, 044026 [NASA ADS] [CrossRef] [Google Scholar]

[4] Bar, N., Blas, D., Blum, K., & Sibiryakov, S. 2018, Phys. Rev. D, 98, 083027 [NASA ADS] [CrossRef] [Google Scholar]

[5] Bar, N., Blum, K., Eby, J., & Sato, R. 2019a, Phys. Rev. D, 99, 103020 [NASA ADS] [CrossRef] [Google Scholar]

[6] Bar, N., Blum, K., Lacroix, T., & Panci, P. 2019b, JCAP, 07, 045 [CrossRef] [Google Scholar]

[7] Bar, N., Blas, D., Blum, K., & Kim, H. 2021, Phys. Rev. D, 104, 043021 [NASA ADS] [CrossRef] [Google Scholar]

[8] Bar, N., Blum, K., & Sun, C. 2022, Phys. Rev. D, 105, 083015 [NASA ADS] [CrossRef] [Google Scholar]

[9] Bar-Or, B., Fouvry, J.-B., & Tremaine, S. 2019, ApJ, 871, 28 [Google Scholar]

[10] Baryakhtar, M., Lasenby, R., & Teo, M. 2017, Phys. Rev. D, 96, 035019 [NASA ADS] [CrossRef] [Google Scholar]

[11] Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]

[12] Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143 [NASA ADS] [CrossRef] [Google Scholar]

[13] Bernal, T., Fernández-Hernández, L. M., Matos, T., & Rodríguez-Meza, M. A. 2018, MNRAS, 475, 1447 [NASA ADS] [CrossRef] [Google Scholar]

[14] Bertone, G., & Hooper, D. 2018, Rev. Mod. Phys., 90, 064050 [NASA ADS] [CrossRef] [Google Scholar]

[15] Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]

[16] Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]

[17] Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93 [Google Scholar]

[18] Bosma, A. 1978, PhD thesis, University of Groningen, The Netherlands [Google Scholar]

[19] Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67 [Google Scholar]

[20] Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40 [NASA ADS] [CrossRef] [Google Scholar]

[21] Brandt, T. D. 2016, ApJ, 824, L31 [NASA ADS] [CrossRef] [Google Scholar]

[22] Brito, R., Ghosh, S., Barausse, E., et al. 2017, Phys. Rev. D, 96, 064050 [NASA ADS] [CrossRef] [Google Scholar]

[23] Brook, C. B., Di Cintio, A., Knebe, A., et al. 2014, ApJ, 784, L14 [NASA ADS] [CrossRef] [Google Scholar]

[24] Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]

[25] Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343 [Google Scholar]

[26] Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 539, 517 [CrossRef] [Google Scholar]

[27] Burkert, A. 2015, ApJ, 808, 158 [NASA ADS] [CrossRef] [Google Scholar]

[28] Burkert, A. 2020, ApJ, 904, 161 [NASA ADS] [CrossRef] [Google Scholar]

[29] Calabrese, E., & Spergel, D. N. 2016, MNRAS, 460, 4397 [NASA ADS] [CrossRef] [Google Scholar]

[30] Carballo-Bello, J. A., Sollima, A., Martinez-Delgado, D., et al. 2014, MNRAS, 445, 2971 [NASA ADS] [CrossRef] [Google Scholar]

[31] Chan, M. H., & Fai Yeung, C. 2021, ApJ, 913, 25 [CrossRef] [Google Scholar]

[32] Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981 [NASA ADS] [CrossRef] [Google Scholar]

[33] Chan, J. H. H., Schive, H.-Y., Woo, T.-P., & Chiueh, T. 2018, MNRAS, 478, 2686 [NASA ADS] [CrossRef] [Google Scholar]

[34] Chan, H. Y. J., Ferreira, E. G. M., May, S., Hayashi, K., & Chiba, M. 2022, MNRAS, 511, 943 [NASA ADS] [CrossRef] [Google Scholar]

[35] Chavanis, P.-H. 2011, Phys. Rev. D, 84, 043531 [NASA ADS] [CrossRef] [Google Scholar]

[36] Chavanis, P.-H. 2022, Phys. Rev. D, 106, 043538 [NASA ADS] [CrossRef] [Google Scholar]

[37] Chavanis, P. H., & Delfini, L. 2011, Phys. Rev. D, 84, 043532 [NASA ADS] [CrossRef] [Google Scholar]

[38] Chen, S.-R., Schive, H.-Y., & Chiueh, T. 2017, MNRAS, 468, 1338 [NASA ADS] [CrossRef] [Google Scholar]

[39] Chiang, B. T., Schive, H.-Y., & Chiueh, T. 2021, Phys. Rev. D, 103, 103019 [NASA ADS] [CrossRef] [Google Scholar]

[40] Chiang, B. T., Ostriker, J. P., & Schive, H.-Y. 2023, MNRAS, 518, 4045 [Google Scholar]

[41] Chou, A. S., Soares-Santos, M., Tait, T. M. P., et al. 2022, ArXiv e-prints [arXiv:2211.09978] [Google Scholar]

[42] Church, B. V., Mocz, P., & Ostriker, J. P. 2019, MNRAS, 485, 2861 [NASA ADS] [CrossRef] [Google Scholar]

[43] Correa, C. A., Schaller, M., Ploeckinger, S., et al. 2022, MNRAS, 517, 3045 [NASA ADS] [CrossRef] [Google Scholar]

[44] Dalal, N., & Kravtsov, A. 2022, Phys. Rev. D, 106, 063517 [NASA ADS] [CrossRef] [Google Scholar]

[45] Dave, B., & Goswami, G. 2023, JCAP, 2023, 015 [CrossRef] [Google Scholar]

[46] Davoudiasl, H., & Denton, P. B. 2019, Phys. Rev. Lett., 123, 021102 [Google Scholar]

[47] de Blok, W. J. G., & Bosma, A. 2002, A&A, 385, 816 [CrossRef] [EDP Sciences] [Google Scholar]

[48] Dehghani, R., Salucci, P., & Ghaffarnejad, H. 2020, A&A, 643, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[49] Delgado, V., & Muñoz Mateo, A. 2022, MNRAS, 518, 4064 [CrossRef] [Google Scholar]

[50] Della Monica, R., & de Martino, I. 2023, A&A, 670, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[51] Deng, H., Hertzberg, M. P., Namjoo, M. H., & Masoumi, A. 2018, Phys. Rev. D, 98, 023513 [NASA ADS] [CrossRef] [Google Scholar]

[52] Dentler, M., Marsh, D. J. E., Hložek, R., et al. 2022, MNRAS, 515, 5646 [NASA ADS] [CrossRef] [Google Scholar]

[53] Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415 [Google Scholar]

[54] Di Paolo, C., Salucci, P., & Erkurt, A. 2019, MNRAS, 490, 5451 [NASA ADS] [CrossRef] [Google Scholar]

[55] Diemer, B. 2018, ApJS, 239, 35 [NASA ADS] [CrossRef] [Google Scholar]

[56] Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108 [Google Scholar]

[57] Domcke, V., & Urbano, A. 2015, JCAP, 01, 002 [CrossRef] [Google Scholar]

[58] Donato, F., Gentile, G., Salucci, P., et al. 2009, MNRAS, 397, 1169 [Google Scholar]

[59] Du, X. 2018, Ph.D. Thesis, Georg-August-Universität Göttingen, Germany [Google Scholar]

[60] Du, X., Behrens, C., & Niemeyer, J. C. 2017, MNRAS, 465, 941 [NASA ADS] [CrossRef] [Google Scholar]

[61] Dubinski, J., & Carlberg, R. G. 1991, ApJ, 378, 496 [Google Scholar]

[62] Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]

[63] Efstathiou, G. 1992, MNRAS, 256, 43 [Google Scholar]

[64] Eggemeier, B., & Niemeyer, J. C. 2019, Phys. Rev. D, 100, 063528 [NASA ADS] [CrossRef] [Google Scholar]

[65] Ferreira, E. G. M. 2021, A&ARv, 29, 7 [NASA ADS] [CrossRef] [Google Scholar]

[66] Fitts, A., Michael, B.-K., Elbert, O. D., et al. 2017, MNRAS, 471, 3547 [NASA ADS] [CrossRef] [Google Scholar]

[67] Flores, R. A., & Primack, J. R. 1994, ApJ, 427, L1 [NASA ADS] [CrossRef] [Google Scholar]

[68] Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]

[69] Frenk, C. S., White, S. D. M., Davis, M., & Efstathiou, G. 1988, ApJ, 327, 507 [NASA ADS] [CrossRef] [Google Scholar]

[70] González-Morales, A. X., Marsh, D. J. E., Peñarrubia, J., & Ureña López, L. A. 2017, MNRAS, 472, 1346 [CrossRef] [Google Scholar]

[71] Goodman, J. 2000, New Astron., 5, 103 [NASA ADS] [CrossRef] [Google Scholar]

[72] Goodman, J., & Weare, J. 2013, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]

[73] Gorbunov, D., Khmelnitsky, A., & Rubakov, V. 2008, JCAP, 10, 041 [CrossRef] [Google Scholar]

[74] Guzman, F. S., & Urena-Lopez, L. A. 2004, Phys. Rev. D, 69, 124033 [NASA ADS] [CrossRef] [Google Scholar]

[75] Guzman, F. S., & Urena-Lopez, L. A. 2006, ApJ, 645, 814 [CrossRef] [Google Scholar]

[76] Hayashi, K., Ferreira, E. G. M., & Chan, H. Y. J. 2021, ApJ, 912, L3 [NASA ADS] [CrossRef] [Google Scholar]

[77] Hlozek, R., Grin, D., Marsh, D. J. E., & Ferreira, P. G. 2015, Phys. Rev. D, 91, 103512 [NASA ADS] [CrossRef] [Google Scholar]

[78] Hlozek, R., Marsh, D. J. E., & Grin, D. 2018, MNRAS, 476, 3063 [CrossRef] [Google Scholar]

[79] Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv e-prints [arXiv:1008.4686] [Google Scholar]

[80] Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [NASA ADS] [CrossRef] [Google Scholar]

[81] Hui, L. 2021, ARA&A, 59, 247 [NASA ADS] [CrossRef] [Google Scholar]

[82] Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541 [NASA ADS] [CrossRef] [Google Scholar]

[83] Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, AJ, 144, 134 [CrossRef] [Google Scholar]

[84] Iorio, G., Fraternali, F., Nipoti, C., et al. 2017, MNRAS, 466, 4159 [NASA ADS] [Google Scholar]

[85] Karukes, E. V., & Salucci, P. 2017, MNRAS, 465, 4703 [Google Scholar]

[86] Khelashvili, M., Rudakovskyi, A., & Hossenfelder, S. 2023, MNRAS, 523, 3393 [NASA ADS] [CrossRef] [Google Scholar]

[87] Khlopov, M. I., Malomed, B. A., & Zeldovich, I. B. 1985, MNRAS, 215, 575 [NASA ADS] [CrossRef] [Google Scholar]

[88] Kirby, E. N., Boylan-Kolchin, M., Cohen, J. G., et al. 2013, ApJ, 770, 16 [NASA ADS] [CrossRef] [Google Scholar]

[89] Kirby, E. N., Bullock, J. S., Boylan-Kolchin, M., Kaplinghat, M., & Cohen, J. G. 2014, MNRAS, 439, 1015 [NASA ADS] [CrossRef] [Google Scholar]

[90] Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82 [Google Scholar]

[91] Kobayashi, T., Murgia, R., De Simone, A., Iršič, V., & Viel, M. 2017, Phys. Rev. D, 96, 123514 [NASA ADS] [CrossRef] [Google Scholar]

[92] Kormendy, J., & Freeman, K. C. 2016, ApJ, 817, 84 [NASA ADS] [CrossRef] [Google Scholar]

[93] Krut, A., Argüelles, C. R., Chavanis, P. H., Rueda, J. A., & Ruffini, R. 2023, ApJ, 945, 1 [NASA ADS] [CrossRef] [Google Scholar]

[94] Lancaster, L., Giovanetti, C., Mocz, P., et al. 2020, JCAP, 01, 001 [NASA ADS] [Google Scholar]

[95] Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157 [Google Scholar]

[96] Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2020, ApJS, 247, 31 [Google Scholar]

[97] Li, X., Hui, L., & Yavetz, T. D. 2021, Phys. Rev. D, 103, 023508 [NASA ADS] [CrossRef] [Google Scholar]

[98] Lidz, A., & Hui, L. 2018, Phys. Rev. D, 98, 023011 [NASA ADS] [CrossRef] [Google Scholar]

[99] Macciò, A. V., Paduroiu, S., Anderhalden, D., Schneider, A., & Moore, B. 2012, MNRAS, 424, 1105 [CrossRef] [Google Scholar]

[100] Marsh, D. J. E. 2016, Phys. Rept., 643, 1 [NASA ADS] [CrossRef] [Google Scholar]

[101] Marsh, D. J. E., & Niemeyer, J. C. 2019, Phys. Rev. Lett., 123, 051103 [CrossRef] [Google Scholar]

[102] Matos, T., & Urena-Lopez, L. A. 2001, Phys. Rev. D, 63, 063506 [NASA ADS] [CrossRef] [Google Scholar]

[103] May, S., & Springel, V. 2021, MNRAS, 506, 2603 [NASA ADS] [CrossRef] [Google Scholar]

[104] May, S., & Springel, V. 2023, MNRAS, 524, 4256 [NASA ADS] [CrossRef] [Google Scholar]

[105] McConnachie, A. W. 2012, AJ, 144, 4 [Google Scholar]

[106] Menci, N., Merle, A., Totzauer, M., et al. 2017, ApJ, 836, 61 [NASA ADS] [CrossRef] [Google Scholar]

[107] Mina, M., Mota, D. F., & Winther, H. A. 2022, A&A, 662, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[108] Mocz, P., Lancaster, L., Fialkov, A., Becerra, F., & Chavanis, P.-H. 2018, Phys. Rev. D, 97, 083519 [NASA ADS] [CrossRef] [Google Scholar]

[109] Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2019, Phys. Rev. Lett., 123, 141301 [CrossRef] [Google Scholar]

[110] Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2020, MNRAS, 494, 2027 [NASA ADS] [CrossRef] [Google Scholar]

[111] Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2023, MNRAS, 521, 2608 [NASA ADS] [CrossRef] [Google Scholar]

[112] Moore, B. 1994, Nature, 370, 629 [NASA ADS] [CrossRef] [Google Scholar]

[113] Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19 [Google Scholar]

[114] Moroz, I. M., Penrose, R., & Tod, P. 1998, Class. Quant. Grav., 15, 2733 [NASA ADS] [CrossRef] [Google Scholar]

[115] Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903 [Google Scholar]

[116] Nadler, E. O., Gluscevic, V., Boddy, K. K., & Wechsler, R. H. 2019, ApJ, 878, L32, [Erratum: ApJ 897, L46 (2020), Erratum: ApJ 897, L46 (2020)] [NASA ADS] [CrossRef] [Google Scholar]

[117] Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996a, MNRAS, 283, L72 [NASA ADS] [Google Scholar]

[118] Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996b, ApJ, 462, 563 [Google Scholar]

[119] Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]

[120] Niemeyer, J. C. 2020, Prog. Part. Nucl. Phys., 113, 103787 [NASA ADS] [CrossRef] [Google Scholar]

[121] Nori, M., & Baldi, M. 2021, MNRAS, 501, 1539 [Google Scholar]

[122] Nori, M., Murgia, R., Iršič, V., Baldi, M., & Viel, M. 2019, MNRAS, 482, 3227 [NASA ADS] [CrossRef] [Google Scholar]

[123] Oh, S.-H., de Blok, W. J. G., Brinks, E., Walter, F., & Kennicutt, R. C., Jr 2011, AJ, 141, 193 [NASA ADS] [CrossRef] [Google Scholar]

[124] Oh, S.-H., Hunter, D. A., Brinks, E., et al. 2015, AJ, 149, 180 [CrossRef] [Google Scholar]

[125] Peebles, P. J. E. 2022, Ann. Phys., 447, 169159 [NASA ADS] [CrossRef] [Google Scholar]

[126] Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[127] Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[128] Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [NASA ADS] [CrossRef] [Google Scholar]

[129] Putman, M. E., Zheng, Y., Price-Whelan, A. M., et al. 2021, ApJ, 913, 53 [NASA ADS] [CrossRef] [Google Scholar]

[130] Read, J. I., Agertz, O., & Collins, M. L. M. 2016a, MNRAS, 459, 2573 [NASA ADS] [CrossRef] [Google Scholar]

[131] Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2016b, MNRAS, 462, 3628 [NASA ADS] [CrossRef] [Google Scholar]

[132] Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2017, MNRAS, 467, 2019 [NASA ADS] [Google Scholar]

[133] Ren, T., Kwa, A., Kaplinghat, M., & Yu, H.-B. 2019, Phys. Rev. X, 9, 031020 [NASA ADS] [Google Scholar]

[134] Robles, V. H., Bullock, J. S., & Boylan-Kolchin, M. 2019, MNRAS, 483, 289 [NASA ADS] [CrossRef] [Google Scholar]

[135] Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81 [NASA ADS] [CrossRef] [Google Scholar]

[136] Rodrigues, D. C., del Popolo, A., Marra, V., & de Oliveira, P. L. C. 2017, MNRAS, 470, 2410 [NASA ADS] [CrossRef] [Google Scholar]

[137] Rogers, K. K., & Peiris, H. V. 2021, Phys. Rev. Lett., 126, 071302 [NASA ADS] [CrossRef] [Google Scholar]

[138] Roper, F. A., Oman, K. A., Frenk, C. S., et al. 2023, MNRAS, 521, 1316 [NASA ADS] [CrossRef] [Google Scholar]

[139] Rubakov, V. A., & Gorbunov, D. S. 2017, Introduction to the Theory of the Early Universe: Hot big bang theory (Singapore: World Scientific) [Google Scholar]

[140] Rubin, V. C., Ford, Jr., W. K., & Thonnard, N. 1980, ApJ, 238, 471 [NASA ADS] [CrossRef] [Google Scholar]

[141] Ruffini, R., & Bonazzola, S. 1969, Phys. Rev., 187, 1767 [NASA ADS] [CrossRef] [Google Scholar]

[142] Safarzadeh, M., & Spergel, D. N. 2020, ApJ, 893, 21 [NASA ADS] [CrossRef] [Google Scholar]

[143] Sales, L. V., Wetzel, A., & Fattahi, A. 2022, Nat. Astron., 6, 897 [NASA ADS] [CrossRef] [Google Scholar]

[144] Sánchez Almeida, J., Trujillo, I., & Plastino, A. R. 2020, A&A, 642, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[145] Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931 [Google Scholar]

[146] Schive, H.-Y., Chiueh, T., & Broadhurst, T. 2014a, Nat. Phys., 10, 496 [NASA ADS] [CrossRef] [Google Scholar]

[147] Schive, H.-Y., Liao, M.-H., Woo, T.-P., et al. 2014b, Phys. Rev. Lett., 113, 261302 [NASA ADS] [CrossRef] [Google Scholar]

[148] Schive, H.-Y., Chiueh, T., Broadhurst, T., & Huang, K.-W. 2016, ApJ, 818, 89 [NASA ADS] [CrossRef] [Google Scholar]

[149] Schneider, A. 2018, Phys. Rev. D, 98, 063021 [NASA ADS] [CrossRef] [Google Scholar]

[150] Schutz, K. 2020, Phys. Rev. D, 101, 123026 [NASA ADS] [CrossRef] [Google Scholar]

[151] Schwabe, B., Niemeyer, J. C., & Engels, J. F. 2016, Phys. Rev. D, 94, 043513 [Google Scholar]

[152] Seidel, E., & Suen, W.-M. 1994, Phys. Rev. Lett., 72, 2516 [NASA ADS] [CrossRef] [Google Scholar]

[153] Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [Google Scholar]

[154] Silk, J., Moore, B., & Diemand, J. 2010, in Particle Dark Matter: Observations, Models and Searches, ed. G. Bertone (Cambridge: Cambridge Univ. Press) [Google Scholar]

[155] Simon, J. D., Geha, M., Minor, Q. E., et al. 2011, ApJ, 733, 46 [NASA ADS] [CrossRef] [Google Scholar]

[156] Sin, S.-J. 1994, Phys. Rev. D, 50, 3650 [NASA ADS] [CrossRef] [Google Scholar]

[157] Sommer-Larsen, J., & Dolgov, A. 2001, ApJ, 551, 608 [NASA ADS] [CrossRef] [Google Scholar]

[158] Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [NASA ADS] [CrossRef] [Google Scholar]

[159] Street, L., Gnedin, N. Y., & Wijewardhana, L. C. R. 2022, Phys. Rev. D, 106, 043007 [NASA ADS] [CrossRef] [Google Scholar]

[160] Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542 [CrossRef] [Google Scholar]

[161] Tremaine, S. D. 1976, ApJ, 203, 345 [NASA ADS] [CrossRef] [Google Scholar]

[162] Tremaine, S., & Gunn, J. E. 1979, Phys. Rev. Lett., 42, 407 [NASA ADS] [CrossRef] [Google Scholar]

[163] Trimble, V. 1987, ARA&A, 25, 425 [NASA ADS] [CrossRef] [Google Scholar]

[164] Tulin, S., Yu, H.-B., & Zurek, K. M. 2013, Phys. Rev. D, 87, 115007 [NASA ADS] [CrossRef] [Google Scholar]

[165] Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189 [Google Scholar]

[166] Veltmaat, J., Niemeyer, J. C., & Schwabe, B. 2018, Phys. Rev. D, 98, 043509 [NASA ADS] [CrossRef] [Google Scholar]

[167] Veltmaat, J., Schwabe, B., & Niemeyer, J. C. 2020, Phys. Rev. D, 101, 083518 [NASA ADS] [CrossRef] [Google Scholar]

[168] Vicens, J., Salvado, J., & Miralda-Escudé, J. 2018, ArXiv e-prints [arXiv:1802.10513] [Google Scholar]

[169] Vogelsberger, M., Zavala, J., & Loeb, A. 2012, MNRAS, 423, 3740 [NASA ADS] [CrossRef] [Google Scholar]

[170] Workman, R. L., Burkert, V. D., Crede, V., et al. 2022, PTEP, 2022, 083C01 [Google Scholar]

[171] Zhu, Q., Marinacci, F., Maji, M., et al. 2016, MNRAS, 458, 1559 [Google Scholar]

[172] Zoutendijk, S. L., Brinchmann, J., Bouché, N. F., et al. 2021, A&A, 651, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]