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/00046361/202346686  
Published online  07 August 2023 
Confronting fuzzy dark matter with the rotation curves of nearby dwarf irregular galaxies
^{1}
Instituto de Astrofísica de Canarias, La Laguna, Tenerife 38200, Spain
email: abh@iac.es
^{2}
Departamento de Astrofísica, Universidad de La Laguna, La Laguna, Tenerife 38200, Spain
^{3}
Dipartimento di Fisica e Astronomia Galileo Galilei, Università di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy
^{4}
INFNPadova, Via Marzolo 8, 35131 Padova, Italy
^{5}
INAFPadova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
Received:
18
April
2023
Accepted:
18
May
2023
Aims. In this paper, we carry out a phenomenological investigation of the viability of fuzzy dark matter, which is composed of coherent waves of noninteracting ultralight axions with a mass of m_{a} ≈ 10^{−22} 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 highresolution rotation curves from the LITTLE THINGS in 3D catalog. This comprises a collection of isolated, dark matterdominated dwarfirregular 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 χ^{2} 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 corehalo 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^{−23} 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 radiusmass 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. GarciaFarieta 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 Nbody 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 DMdominated 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 “corecusp” puzzle that is among the smallscale problems related to ΛCDM (see Bullock & BoylanKolchin 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 lowmass 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 “toobigtofail” problem, whereby the central masses are too low compared to the most massive (sub)halos predicted in ΛCDM (BoylanKolchin 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^{−23} − 10^{−21} eV (Sin 1994; Hu et al. 2000; Matos & UrenaLopez 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 selfgravitating DM waves. One of the consequences of this is the appearance of a pressurelike 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 CDMlike 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 lowmass 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 smallscale issues of ΛCDM, which has motivated several phenomenological studies and searches in astrophysical and cosmological probes of lowscale structure (see Ferreira 2021 for a review and Sect. 5 below for a discussion).
Other attempts to solve some of the smallscale issues of ΛCDM by modifying the fundamental properties of DM are warm DM (SommerLarsen & 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), selfinteracting FDM Goodman 2000; Chavanis 2011; Chavanis & Delfini 2011; Delgado & Muñoz Mateo 2022 and, more generally, selfinteracting 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 lowmass 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^{8} M_{⊙} ≲ M_{⋆} ≲ 10^{9} M_{⊙} (Tollet et al. 2016) and it could still play a role for masses all the way down to M_{⋆} ∼ 10^{6} M_{⊙}^{1}.
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 lowmass 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 stateoftheart 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 highresolution rotation curves of galaxies that lie on the edge of galaxy formation. It provides an excellent laboratory to investigate the smallscale 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 KleinGordon equation coupled to the gravitational field, Φ(t, r)^{2}. In the nonrelativistic limit, it is convenient to decompose it as 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 . The field ψ verifies, then, the SchrödingerPoisson (SP) system of classical wave equations,
Space and time dependence of the field factorize in the quasistationary solutions of these equations,
We have chosen the normalization such that the physical density of the field is
and have introduced the Planck mass m_{P} = G^{−1/2} ≈ 10^{19} GeV for convenience. Assuming spherical symmetry, χ obeys a simpler version of the SP equations,
where we have rescaled the radial extension as x = m_{a}r (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,
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 groundstate 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:
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:
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):
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 & UrenaLopez 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 NavarroFrenkWhite (NFW) density profile characteristic of CDM (Navarro et al. 1997):
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:
where ρ_{m0} denotes the cosmological background matter density (ρ_{m0} = Ω_{m0}ρ_{c0}, with ), 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)^{3}. 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
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 “corehalo relation” was originally found by Schive et al. (2014a,b) and can be written as
which depends on the redshift, z. At the moment there is no comprehensive understanding of the physical mechanism underlying the corehalo 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 solitonmerger simulations. This resulted in a more general relation that at z = 0 can be parametrized as:
where and are the bestfit 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 NFWlike 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:
where γ_{1} = 15, γ_{2} = 0.3 and we have the characteristic mass scale M_{0} = 1.6 × 10^{10}(m_{a}/10^{−22} 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:
with a scatter of 0.16 dex^{4}.
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 NFWlike 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 NFWlike 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 smallscale 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^{10} M_{⊙}. This can be derived by a simple parametrization derived from numerical simulations (Schive et al. 2016; May & Springel 2021, 2023):
where, as previously for the concentration – mass relation, M_{0} = 1.6 × 10^{10}(m_{a}/10^{−22} eV)^{−4/3} M_{⊙}, α = 1.1 and β = 2.2. The characteristic mass M_{0} 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 highresolution 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 corecusp problem. We considered a representative and robust subset of these galaxies, with stellar masses in the range of 10^{5} M_{⊙} ≲ M_{⋆} ≲ 10^{8} 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^{5}. A fraction of these galaxies termed “rogues” were flagged and excluded from the core LTs data set because they were nearly faceon (five inclination rogues), their measured distance was not welldetermined (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^{8} 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^{6}. We also note that some galaxies in our sample have significant errors in the inclination values. However, this is selfconsistently 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:
We use the numerical groundstate 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 semianalytically 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 nontrivial 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,
In this equation , 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 fourvector 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),
with M_{⋆/gas} and R_{*/gas} being the mass and the exponential scale length of the disc, respectively, and with I(y) = y^{2}[I_{0}(y)K_{0}(y)−I_{1}(y)K_{1}(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 loglikelihood
(up to an additive constant), where V_{obs}(r_{i}) and δV_{obs}(r_{i}) are the velocity and corresponding uncertainty of the ith 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 affineinvariant ensemble sampler for Markov chain Monte Carlo (MCMC) (ForemanMackey et al. 2013), based on the algorithm proposed by Goodman & Weare (2013). To deal with possible multimodal posteriors and enhance sampling efficiency, we implemented a combination of MCMC moves, which significantly reduces autocorrelation times^{7}. 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 burnin 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 ForemanMackey 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 lognormal distribution. A similar type of prior was used for the angle i for inclination rogues, using a normal distribution instead of a lognormal 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 corehalo 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 solitondominated fits (for a given m_{a} and M_{c}). For the corehalo 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.
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 chisquare 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^{8}.
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. 
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 chisquare of the maximum posterior fit to the LTs catalog, including some of the rogues.
As shown by the reduced chisquare 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 chisquare 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 corehalo 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 solitonhalo 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 chisquares 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^{7} 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^{7} 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 fastexpanding 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^{−25} eV to 10^{−19} 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^{−23} eV − 10^{−22} 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 dotdashed 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_{10}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^{−23} eV, corresponding to an uncertainty of 0.019 dex, and a reduced chisquare . 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 (Workman et al. 2022). With this method we obtain a more conservative error:
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^{9}.
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, . We then performed a linear fit of and β_{⋆} to the 10 (log_{10}M_{⋆}, log_{10}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 , 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 , 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 socalled “column density,” Σ_{0} = ρ_{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 which undershoots slightly the empirical value pc^{−2} (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. Dotdashed 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 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 nontrivial 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^{−23} eV, disregarding the difficulties of the model to explain the data discussed in the previous section. First, we show on the topleft panel of Fig. 6 our results in the M_{h} − M_{c} plane compared to the corehalo 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 topright 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^{−23} eV. Graycolored points are galaxies for which both median and maximum posterior values of r_{t} are within the observed range, while goldcolored are solitondominated 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. Topleft: 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). Topright: 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). Bottomleft: Stellartohalo 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; dotdashed line) with a reference σ = 0.30 dex scatter in both cases. Bottomright: Expected number of halos with M_{h} ≲ 10^{11} 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^{9} ≲ M_{h} ≲ 10^{11} M_{⊙}. The dotdashed 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 lowmass 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 nonmonotonic behavior of the HMF in Eq. (18), we follow a method inspired by Vale & Ostriker (2004), whereby we find a M_{⋆} − M_{h} stellartohalo mass function by matching the cumulatives of the SMF and HMF,
In the bottomleft panel of Fig. 6 we show the results of abundance matching using the stellartohalo 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 stellartohalo mass function features a flattening starting at values as high as M_{h} ≈ 5 × 10^{11} 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 stellartohalo 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 bottomright panel of Fig. 6, we show the number of FDM halos with M_{h} ≲ 10^{11} 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 ShethTormen 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)^{10}. Using the former, we find a very conservative lower bound of m_{a} ≳ 4.3 × 10^{−23} 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^{−22} eV.
Finally, we can also study the population of dwarf galaxies in the MilkyWay’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 kpc (BlandHawthorn & Gerhard 2016)), we obtain the stronger lower bound m_{a} ≳ 7.3 × 10^{−22} 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 DMdominated 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 quasistationary (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
We can then apply the same quasistationary 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 DMdominated galaxies in Bar et al. (2019b), where axiallysymmetric baryonic distributions were implemented in the SP equations. The spherically averaged baryonic density is computed from the baryonic mass profile
using Eq. (21) so that
where we compute the derivative analytically.
Our initial approach to assess the selfconsistency 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 DMonly 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 dotdashed 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 DMonly soliton, for instance, the one approximated by Eq. (9). This does not imply necessarily that the two solitons (DMonly 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.
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 innerslope 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^{7} 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^{7} 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^{−23} 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 lowsurface 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 χ^{2} 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^{−23} eV. The soliton or solitonhalo 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 corehalo relation. In particular, the appearance of a “doublebump” structure in the rotation curves is predicted and then used to set an exclusion of FDM in the range of m_{a} ∼ 10^{−22} − 10^{−21} eV range. This structure is not a generic feature appearing in our analyses due to the broader range of corehalo 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 fullfledged fit to the data and uses exclusively the corehalo relation as a exclusion criterion for FDM in the range m_{a} ∼ 10^{−23} − 10^{−20} 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 corehalo 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 smallscale 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^{−24} eV.
The Lymanα forest is a powerful tracer of the linear matter power spectrum for subMpc scales, where the FDM suppression arises (Niemeyer 2020). With highresolution 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^{−21} 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^{−21} 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 (DESY1) 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^{−23} eV at 95% CL.
The abundances of lensed ultrafaint 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^{−22} at 3σ CL.
We now include results from astrophysical searches and bounds. For instance, in the regime of zero selfcoupling, the nondetection of blackhole 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^{−20} eV < m_{a} < 10^{−16} eV at the 95% CL, while stellar mass black holes yield an exclusion region 7 × 10^{−14} eV < m_{a} < 2 × 10^{−11} 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^{−21} eV < m_{a} < 4.6 × 10^{−21} eV (Davoudiasl & Denton 2019).
The FDM’s wave behavior would heat up stellar objects in galaxies (BarOr 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^{−23} 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^{−22} eV. By accounting for different sources of uncertainties, Chiang et al. (2023) finds that the range of m_{a} ∼ 0.5 − 0.7 × 10^{−22} eV is favored by the observed MilkyWay 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 socalled Fornax timing puzzle Tremaine (1976) as long as m_{a} > 10^{−21} eV (Lancaster et al. 2020).
Stellar streams are excellent tracers of the gravitational potential (CarballoBello et al. 2014) and can be used to test predictions of the subhalo mass function in different DM models. The fluctuations in the stellar streams lead to imposing the lower limit m_{a} > 2 × 10^{−21} 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^{−21} eV at a 95% CL.
The suppression of the cosmic growth in FDM reduces the number of lowmass halos, which impinges on the dwarfspheroidal 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^{−21} eV (Nadler et al. 2019).
By applying a Jeans analysis for Milky Waysatellite 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 eV (at 2σ). Comparing to simulated data on the lineofsight velocity dispersion (σ_{LOS}) of Fornax and Sculptor dSphs, GonzálezMorales et al. (2017) found a tighter upper bound m_{a} < 4 × 10^{−23} eV at 97.5% CL.
The existence of the subhalo to host Eridanus II lowsurface brightness dwarf galaxy and the survival of its star cluster can also constrain the FDM mass (BarOr et al. 2019; Brandt 2016). As shown in Marsh & Niemeyer (2019), the subhalo mass function in the Milky Way and the existence of Eridanus II implies m_{a} ≳ 8 × 10^{−22} 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^{−19} eV with the possible exception of the range 10^{−21} eV ≤ m_{a} ≤ 10^{−19} 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 ultrafaint 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^{−1} 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^{−19} 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 eV, which is consistent with the results obtained directly from the MUSEFaint survey (Zoutendijk et al. 2021). Finally, Calabrese & Spergel (2016) found that FDM particle with a mass of m_{a} ∼ 3.7 − 5.6 × 10^{−22} eV is in agreement with the data of halflight mass of Draco II and Triangulum II; although Safarzadeh & Spergel (2020) later found an upper limit m_{a} ≲ 10^{−21} 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 dwarfgalaxy 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 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 highresolution rotation curves from the LITTLE THINGS in 3D (or LTs for short in our paper) catalog. This comprises a collection of isolated, dark matterdominated dwarfirregular galaxies that provides an optimal benchmark for cosmological studies. The methodological basis of our study is a statistical framework based on the χ^{2} 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 corehalo 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^{−23} 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 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^{−23} 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^{−23} 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^{−22} 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 galaxydependent 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 corecusp problem or any of the other smallscale 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 socalled “catch22 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 corelike structures in dwarf galaxies by invoking a particle mass that would imply a suppression of smallscale 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, selfinteractions, as recently explored in Mocz et al. (2023) as a possible mechanism to enhance smallscale 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 firstprinciples 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 baryonicfeedback induced cores for dwarf galaxies while those from Read et al. (2016a) find the cores for all of them (see Bullock & BoylanKolchin 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 & UrenaLopez (2004), Moroz et al. (1998), Hui et al. (2017), Ruffini & Bonazzola (1969). In our formulae, we implicitly use natural units ℏ = c = 1.
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.
See https://emcee.readthedocs.io/en/stable/tutorials/moves/, where this is explained further using the same combination of moves.
We did this as such data were not explicitly available in Behroozi et al. (2019) due to their use of a nonparametric 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 smallscale 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íaFarieta, Théo Gayoux, SeHeon Oh, Marc HuertasCompany, Javier Olivares and HsiYu 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 PGC2018102016AI00 and J.M.C. is also funded by the “Ramón y Cajal” program RYC201620672. 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., PalanqueDelabrouille, 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]
 BarOr, 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ándezHernández, L. M., Matos, T., & RodríguezMeza, 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]
 BlandHawthorn, 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]
 BoylanKolchin, 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., & BoylanKolchin, 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]
 CarballoBello, J. A., Sollima, A., MartinezDelgado, 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., SoaresSantos, M., Tait, T. M. P., et al. 2022, ArXiv eprints [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, GeorgAugustUniversitä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]
 ForemanMackey, 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álezMorales, 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., & UrenaLopez, L. A. 2004, Phys. Rev. D, 69, 124033 [NASA ADS] [CrossRef] [Google Scholar]
 Guzman, F. S., & UrenaLopez, 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 eprints [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., FicutVicas, 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., BoylanKolchin, M., Cohen, J. G., et al. 2013, ApJ, 770, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Kirby, E. N., Bullock, J. S., BoylanKolchin, 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., & UrenaLopez, 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., PriceWhelan, 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., & BoylanKolchin, 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]
 SommerLarsen, 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., & MiraldaEscudé, J. 2018, ArXiv eprints [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 autocorrelation 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. Highresolution 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 lineofsight and the vector normal to the plane of the galaxy (Read et al. 2016b, 2017). Hence, a galaxy is edgeon when i = 90° and faceon 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:
where i_{0} is a “reference” angle. Therefore, for the regime of i < 40° (nearly faceon galaxies), small errors on i can have a strong impact on the reconstruction of the observed rotational velocity. The i_{0} 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},
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_{0} 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 loglikelihood is necessary since now the errors in eq. (A.1) are idependent,
where ℒ′(θ_{FDM},M_{⋆}) is the loglikelihood 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^{8} 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 (DDO 50) and (DDO 133), which agree within < 5% with the central values of i_{0} 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.
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_{10}m_{a}, log_{10}M_{⋆}), (log_{10}ρ_{c}, log_{10}r_{c}) and (log_{10}r_{c}, log_{10}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
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_{10}m_{a}, log_{10}M_{⋆}) plane is given by:
while the resulting projected error Σ_{i} is given by
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 loglikelihood function is then given by:
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 48walker 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 000step burnin 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 we can define the total number of halos within a given mass range M_{1} < M < M_{2} and volume V as
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 ShethTormen 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 ShethTormen HMF for CDM (Sheth & Tormen 1999) normalized to match the same total abundance as the previous one (solid, red). We also include the original ShethTormen 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 ShethTormen HMF for each M_{h} and mapped it to the corresponding M_{⋆} via the abundance matching relation of Behroozi et al. (2019). ^{11} We then subsequently performed the abundance matching procedure with the HMF of FDM (without the normalization factor discussed earlier in the appendix). The ShethTormen HMF is computed via the Colossus package in Python developed by Diemer (2018).
Appendix C: Comparison to the solitonhost halo relation
As discussed in Sect. 5, Bar et al. (2018) independently rederived a corehalo 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 “doublebump” in the rotation curves, the first of which should become prominent in the inner part of the galaxies for m_{a} ≳ 10^{−22} eV. In our analysis we implement the corehalo 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 corehalo 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),
This can be equivalently expressed as (Bar et al. 2018)
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:
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_{200}.
Replacing Eq. (C.3) in Eq. (C.1), one obtains the corehalo relation reported by Bar et al. (2018),
where
is a function of c that varies between 0.9 ≲ f(c)≲1.1 for 5 ≤ c ≤ 30.
The corehalo 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_{200}, M_{200})→(c_{100}, M_{100}) that parametrizes the same NFW density profile in Eq. (10). Transforming M_{200} to M_{h} increases the given halo mass by 20 − 10% for c_{200} in the range 5 − 30, yielding
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
and applying it to the empirical relation in Eq. (C.1). ^{12} 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 selfgravitating 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 corehalo relations, leading to heavy solitons, is not favored by our fits for relatively heavy axion masses, m_{a} ≳ 10^{−22} 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^{−22} 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 corehalo relations in the rotation curve of DDO 154. Theoretical rotation curves are DMonly (i.e., excluding baryonic components) for visual clarity. 
All Tables
Priors used to fit the DM and baryon contributions to the rotation curve data from the LTs catalog.
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 chisquare of the maximum posterior fit to the LTs catalog, including some of the rogues.
Systematic correction in FDM mass central value as a result of introducing the modified soliton due to a background baryonic potential.
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.
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 chisquares 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 dotdashed 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. Dotdashed 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^{−23} eV. Graycolored points are galaxies for which both median and maximum posterior values of r_{t} are within the observed range, while goldcolored are solitondominated 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. Topleft: 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). Topright: 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). Bottomleft: Stellartohalo 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; dotdashed line) with a reference σ = 0.30 dex scatter in both cases. Bottomright: Expected number of halos with M_{h} ≲ 10^{11} 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^{9} ≲ M_{h} ≲ 10^{11} M_{⊙}. The dotdashed 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 dotdashed 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^{−23} 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 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. Highresolution 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 ShethTormen HMF for CDM (Sheth & Tormen 1999) normalized to match the same total abundance as the previous one (solid, red). We also include the original ShethTormen abundance without the normalization. 

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

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