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

,


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 et al. 2020;Alam & Ata 2017;Workman et al. 2022;Peebles 2022).However, observations probing the formation of structure at small galactic scales might pose a challenge to ΛCDM.
Stellar and gas kinematics indeed provide one of the most direct observational probes of the internal structure of DM halos (Bosma 1978;Rubin et al. 1980;Trimble 1987).While pure N-body CDM simulations generically predict density profiles peaking as ρ(r) ∝ 1/r at small radii (Dubinski & Carlberg 1991;Navarro et al. 1996Navarro et al. , 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. 2011Oh et al. , 2015;;Read et al. 2017;Li et al. 2020;Dehghani et al. 2020).This is known as the "core-cusp" puzzle that is among the small-scale problems related to ΛCDM (see Bullock & Boylan-Kolchin (2017) and Sales et al. (2022) for recent reviews).Other possible discrepancies between the predictions of the model and observations are referred to as the "missing satellites" problem that is related to the observation of many fewer low-mass galaxies in the Local Group than otherwise expected from CDM simulations (Klypin et al. 1999;Moore et al. 1999;McConnachie 2012).Then there is the "too-big-to-fail" problem, whereby the central masses are too low compared to the most massive (sub)halos predicted in ΛCDM (Boylan-Kolchin et al. 2011;Kirby et al. 2014).
These problems have prompted a flurry of activity over the past decade theorizing and analyzing phenomenologically models beyond CDM.One prominent candidate is fuzzy DM (FDM), which is composed of ultralight axions with masses in the range of ∼ 10 −23 −10 −21 eV (Sin 1994;Hu et al. 2000;Matos & Urena-Lopez 2001;Chavanis 2011;Schive et al. 2014a,b).This type of particle has been predicted by several theoretical models, preserving the features of CDM at large cosmological scales and which can be produced "naturally" in the early universe to match the DM relic abundance (Marsh 2016;Hui et al. 2017;Niemeyer 2020;Ferreira 2021;Hui 2021).On the other hand, on distances comparable to their de Broglie wavelength, which can be in the kiloparsec (kpc) scale, these particles populate the galactic halos with large occupation numbers and behave as self-gravitating DM waves.One of the consequences of this is the appearance of a pressure-like effect on macroscopic scales leading to the formation of a flat core, or "soliton," at the center of galaxies with a relatively marked transition to a less dense outer region that follows a CDM-like distribution (Schive et al. 2014a,b;Veltmaat et al. 2018;Mocz et al. 2018;Nori & Baldi 2021;Mina et al. 2020;Chan et al. 2022).This effect also produces a suppression of the linear power spectrum on small scales and to a decrease of the number of low-mass halos (Khlopov et al. 1985;Schive et al. 2016;Du et al. 2017;May & Springel 2021, 2022).The FDM provides, then, an elegant solution to the small-scale issues of ΛCDM, which has motivated several phenomenological studies and searches in astrophysical and cosmological probes of low-scale structure (see Ferreira (2021) for a review and Sect. 5 below for a discussion).
Other attempts to solve some of the small-scale issues of ΛCDM by modifying the fundamental properties of DM are warm DM (Sommer-Larsen & Dolgov 2001;Bode et al. 2001;Rubakov & Gorbunov 2017) including degenerate fermion DM (Tremaine & Gunn 1979;Gorbunov et al. 2008;Domcke & Urbano 2015;Bar et al. 2021;Chavanis 2022;Krut et al. 2023), self-interacting FDM Goodman (2000); Chavanis (2011); Chavanis & Delfini (2011); Delgado & Muñoz Mateo (2022) and, more generally, self-interacting DM (Spergel & Steinhardt 2000;Rocha et al. 2013;Vogelsberger et al. 2012;Tulin et al. 2013;Ren et al. 2019;Correa et al. 2022), as well as phenomenological approaches (see e.g., Sánchez Almeida et al. ( 2020)).However, these problems could also be the result of baryonic physics at play in galaxy formation and evolution within ΛCDM (Sales et al. 2022).For instance, the missing satellites could be due to a loss of efficiency in galaxy formation within low-mass halos as a result of the suppression of gas accretion during reionization (Efstathiou 1992;Bullock et al. 2000;Sawala et al. 2016).In addition, baryonic feedback from supernovae (SN) could provide a mechanism that is powerful enough to flatten CDM profiles at the core of dwarf galaxies (Navarro et al. 1996;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. 2016).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 ⊙ . 1he 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); Robles et al. (2019); Bar et al. (2022); Chan & Fai Yeung (2021); Street et al. (2022); Khelashvili et al. (2022); 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. (2016), which were obtained using state-of-the-art 3D reconstruction techniques to reanalyze the HI data cubes.Galaxies that could have a biased dynamical analysis, dubbed "rogues," were flagged and excluded from the core data set.This catalog, referred to as LITTLE THINGS in 3D (Iorio et al. 2017), constitutes a robust and homogeneous collection of high-resolution rotation curves of galaxies that lie on the edge of galaxy formation.It provides an excellent laboratory to investigate the small-scale issues of ΛCDM and their potential interconnections with the nature of DM.

FDM Generalities
We start by considering a real scalar field, ϕ(t, r), with a mass, m a , satisfying a Klein-Gordon equation coupled to the gravitational field, Φ(t, r).2In the non-relativistic limit, it is convenient to decompose it as ϕ(t, r) = 1/

√
2m a e −im a t ψ(t, r) + c.c. 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 | ψ| ≪ m|ψ|.The field ψ verifies, then, the Schrödinger-Poisson (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  4) for the ground state and first two excitations, where the scale invariance was used to set χ(0) = 1. 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.Eq. ( 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.
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 & Urena-Lopez 2006;Schwabe et al. 2016;Hui et al. 2017;Li et al. 2021), while the outer region is expected to form a virialized halo.This configuration of the FDM has been confirmed by simulations (Schive et al. 2014a,b;Veltmaat et al. 2018;Mocz et al. 2018;Nori & Baldi 2021;Chan et al. 2022).The outer halo follows the universal Navarro-Frenk-White (NFW) density profile characteristic of CDM (Navarro et al. 1997): 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 ρ c0 = 3H 2 0 /8G), 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 Ade et al. (2016), with h = 0.678, Ω m0 = 0.308 (all our calculations are at z = 0). 3An 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 "core-halo 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 core-halo relation, although some ideas have been put forward in Schive et al. (2014b); Hui et al. (2017); Bar et al. (2018); Eggemeier & Niemeyer (2019).In fact, subsequent simulations have found that there is a significant scatter that was studied in detail in Chan et al. (2022) by using a large sample of cosmological and soliton-merger simulations.This resulted in a more general relation that at z = 0 can be parametrized as: where β = 8.00 +0.52  −6.00 × 106 M ⊙ , log 10 (γ/ M ⊙ ) = −5.73+2.38  −8.38 and α = 0.515 +0.130  −0.189 are the best-fit parameters obtained in Chan et al. (2022) with their respective errors determining the scatter.
Another feature of FDM halos is that the characteristic flatness in their density profiles is not limited to the core region where the soliton dominates, but is also reflected in their 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 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. 4s a point of clarification, we note that the concentration parameter is defined in Eqs. ( 10) and ( 13) with respect to the NFW component of the density profile.This is irrespective of whether the point r = r s , which can be generally defined as the point where the logarithmic slope equals −2, is located in the NFW-like region or not, meaning that the two definitions are not necessarily equivalent.Nevertheless, as we discuss in Sect. 3 for the halos that we model in this work, we find that the median values of r s lie in the NFW-like regions and the two definitions are consistent with each other in our case.
A final property that is characteristic of FDM and is particularly useful in the context of small-scale structure issues of ΛCDM is the suppression of the primordial power spectrum that inhibits the formation of halos below a certain mass scale (Khlopov et al. 1985;Hu et al. 2000).Moreover, the subhalo population is expected to be further reduced by dynamical effects such as tidal disruption (Hui et al. 2017).This offers an alternative "exotic" solution to the missing satellites problem and, by extension, to the associated too big to fail problem.One observable statistic that is relevant to test this suppression of formation of structures at small scales is the halo mass function (HMF) for low halo masses of M h ≲ 10 10 M ⊙ .This can be derived by a simple parametrization derived from numerical simulations (Schive et al. 2016;May & Springel 2021, 2022): 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.

LITTLE THINGS in 3D
Local Irregulars That Trace Luminosity Extremes, The HI Nearby Galaxy Survey (LITTLE THINGS) is a Very Large Array HI survey which has produced high-resolution rotation curves for 26 nearby dwarf irregular galaxies within 11 Mpc of the local volume (Oh et al. 2015;Hunter et al. 2012).Most of these curves are able to resolve the inner part (within ∼ 1 kpc of the center) of the galaxies, making this survey an optimal laboratory to investigate the core-cusp problem.We considered a representative and robust subset of these galaxies, with stellar masses in the range of 105 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. 5A fraction of these galaxies termed "rogues" were flagged and excluded from the core LTs data set because they were nearly face-on (five inclination rogues), their measured distance was not well-determined (one distance rogue) or because they presented indications of being far from equilibrium (two disequilibrium rogues).Our fiducial cosmological analysis will include 11 of the remaining galaxies (first corresponding lines in Table 2 of Read et al. (2017)).For these galaxies, we will consider the measured rotation velocities in the whole measured range, which are quite regular even if some of them manifest the effect of HI holes, which could indicate recent (fast expanding) or old (quiescent) stellar activity.The disequilibrium rogues present irregular rotation curves while the distance rogue, DDO 101, could have a stellar mass quite above M ⋆ ∼ 10 8 M ⊙ , depending on the assumed distance, see Read et al. (2016).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. 6We also note that some galaxies in our sample have significant errors in the inclination values.However, this is self-consistently accounted for in the rotation curve errors (Iorio et al. 2017), so that it is unnecessary (and potentially inconsistent) to use the inclination angle as a fitting parameter for LTs.

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 ground-state solution to the SP Eqs.(4) for ρ sol (r) and the NFW density profile in Eq. ( 10) for ρ NFW (r).We Notes.The variables used directly in the MCMCs are log 10 m a , log 10 M c , log 10 c, log 10 M h , log 10 M ⋆ , and log 10 (sin i) in the case of inclination rogues, with all other variables being derived from these.M sol corresponds to the total mass of the soliton and it is related to M c via M c ≃ 0.24 M sol .r max denotes the maximum observable radius from the rotation curve data.In the case of stellar mass, M ⋆0 and δM ⋆± denote the central values with errors reported in Read et al. (2017) from population synthesis modeling, while i 0 , δi ± correspond to the same quantities for the case of the inclination that was obtained through their analysis of rotation curve data.
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 O(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. 2019Mocz et al. , 2020;;Veltmaat et al. 2020) with the general result that, for the same given m a , a soliton of a given mass is still formed but more compressed and denser than in its pure FDM configuration.This effect should be especially relevant for galaxies whose central masses are dominated by the baryons and it can be estimated semi-analytically as discussed in Bar et al. (2018Bar et al. ( , 2019a) ) and followed in Sect. 5.
Moreover, the structure of the boundary between the soliton and the halo at r t can be modeled differently.This can have a non-trivial effect (Bernal et al. 2018;Vicens et al. 2018) espe-cially 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 V 2 FDM (r; θ FDM ) = GM FDM (r; θ FDM )/r, where M FDM (r; θ FDM ) is the enclosed mass within radius r obtained from the density distribution in Eq. ( 19) and θ FDM = (m a , M c , c, M h ) is the four-vector of FDM fitted parameters.To describe the baryonic components in Eq. ( 20), we follow the reference Read et al. (2017) and model the distributions of gas and stars as exponentially thin discs (Binney & Tremaine 2008), 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 log-likelihood (up to an additive constant), where V obs (r i ) and δV obs (r i ) are the velocity and corresponding uncertainty of the i-th data point in a galaxy's rotation curve (the summation is over the radial data bins).We also assumed uncorrelated Gaussian errors for the data Table 2. Median and 68% confidence level (CL) of the posterior distributions of the parameters m a , M c , c, M h , and M ⋆ , and the reduced chi-square of the maximum posterior fit to the LTs catalog, including some of the rogues.Notes.We also include the derived core radius, transition radius and central density of the distribution.We flag the galaxy NGC 6822 with * because it was not initially selected as a rogue galaxy (see text for details).The transition radii with an upper uncertainty equal to 0 indicates that less than 84% of the posterior distribution lies below the maximum observable radius r max , while for those not showing an uncertainty band and marked by †, at least 84% of the posterior distribution lies above r max .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. 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).

Galaxy
To effectively scan the parameter space of the FDM model, we used the emcee package, an affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) (Foreman-Mackey et al. 2013), based on the algorithm proposed by Goodman & Weare (2013).To deal with possible multi-modal posteriors and enhance sampling efficiency, we implemented a combination of MCMC moves, which significantly reduces autocorrelation times. 7We also initialized MCMCs under different local optima to ensure that the posterior distributions converged to the global one.We ran the MCMCs with 48 walkers for a total of 20 000 steps, discarding the first 10 000 as a conservative burn-in pe- Fig. 3. Rotation curves in FDM for two representative galaxies in the LTs catalog: WLM and DDO 154.We plot the total circular velocity and reduced chi-squares resulting from the maximum posterior fits to the data when fixing the axion mass at different values.We compare the results also to the fit to a NFW profile.
riod.We found these to produce consistent results under different runs, as well as meeting the criterion of exceeding 100 autocorrelation times, with acceptance fractions in the ∼ 0.2 − 0.4 range, in agreement with the criteria established in Foreman-Mackey et al. ( 2013).
The priors used for the initial analysis in this study (where the FDM mass is allowed to vary for each galaxy) are summarized in Table 1.The prior on M ⋆ is designed to yield results consistent with stellar population synthesis models, as reported in Read et al. (2016), using the same central values and defining the errors given as the 68 % CL bounds under a log-normal distribution.A similar type of prior was used for the angle i for inclination rogues, using a normal distribution instead of a 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 core-halo relation from Eq. ( 15) and demanding r t ≥ r c , (when r t < r max , where r max is the maximum radius where rotation curves are measured), which are both expected from simulations.Fits with r t ≥ r max would yield identical soliton-dominated fits (for a given m a and M c ).For the core-halo relation we set z = 0 using Eq. ( 15), where we have assumed that the uncertainty on the formation history of the halo, see Burkert (2020), is within the intrinsic scatter comprised by this formula.

Fits to rotation curves
The results of the fits to the rotation curves of the 11 LTs galaxies and the 2 selected inclination rogues (DDO 50 and DDO 133) are shown in Table 2.For each galaxy we list the median and 68% CL region of the posterior distributions of the five fitted parameters, M ⋆ and θ FDM , the reduced chi-square and a subset of derived quantities characterizing the FDM distribution, including the core radius, r c , the transition radius, r t , and the central density, ρ c .In Fig. 2, we show the results of the fits for two representative galaxies in our sample, WLM and DDO 154.Besides the data, we show the median, the 68% and 95% CL of the posterior distribution of the total theoretical rotation curve and the individual contributions from the FDM and baryonic density distributions for the central values of the parameters.Further details, including the corner plots and mean autocorrelation times of the MCMC for the fits of these two galaxies are shown for reference in the Appendix A. The full set of rotation curves, corner plots, and mean autocorrelation times obtained from the fits in Table 2 are publicly available under the category Variable Axion Mass in the FuzzyDM GitHub repository.
As shown by the reduced chi-square values obtained in the fits, the performance of the FDM model describing the rotation curve data is, in general, excellent.In two cases, DDO 52 and UGC 8508, the value of the chi-square is too low suggesting there is overfitting or overestimated errors in the data.Our fits favor large solitons and with transition radii to the halo typically overlapping with r max .This reduces the power of the data to constrain the halo parameters c and M h , explaining their large uncertainties.The direct relation between the scale density of the NFW distribution and c in Eq. ( 13) imposes an upper bound on this parameter because, given a soliton describing the core of a galaxy, the density at r ≤ r t cannot be arbitrarily large.In case of the halo mass an increased constraining power is achieved by imposing the core-halo relation in the priors.On the other hand, the parameters of the soliton, m a and M c , are quite well determined, with relative uncertainties even reaching O(10%) in some cases.The sensitivity to the axion mass of these data is illustrated in Fig. 3 for the rotation curves of WLM and DDO 154.The fit is also compared to the to a pure NFW profile that is included for reference.The description of the data can be quite insufficent, and even worse than for NFW for some choices of m a .In these cases, the fit does not converge to the one of NFW (that would result from suppressing the soliton by effectively running r t to zero) because the soliton cannot be arbitrarily light due to the soliton-halo relation imposed in the priors.
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 Fig. 4. Axion masses obtained from the posterior distributions of the fits to the rotation curves of the LTs galaxies where we show medians and 68% and 95% CL errors.The 10 galaxies entering our fiducial analysis are represented by circles while the rogues are represented by triangles.The dot-dashed line and bands represent the mean, 68% and 95% CL regions of a weighted average to these determinations (see Eq. ( 23) and associated discussion).Left: Galaxies are ordered by the central value of the stellar mass reported in Read et al. (2017), which where found to be in good agreement with our posterior distributions.Stellar masses are increasing from left to right.Right: Empirical fit to a power law between m a and M ⋆ showing evidence for a correlation between the two as extracted from the galaxies in our sample.The red solid line denotes the maximum posterior fit with the 68% and 95% CL regions.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 fast-expanding HI bubble was identified in this galaxy for radii below 2.5 kpc (Read et al. 2017).For these reasons, we removed NGC 6822 from our core data sample, adding it to the list of rogue galaxies in Table 2 and for the purposes of the subsequent analyses.In the next subsections, we present different consistency and diagnostic tests on the behavior of the parameters obtained in the fitting procedure.

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.
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 chi-square χ 2 ν = 7.68.This relatively poor value stems from the small uncertainties derived for m a in some of the galaxies (see Fig. 2) and it could be due to unknown systematic uncertainties in the data or the model.Under this premise, there are several ways to account for them in the analysis (Workman et al. 2022;Hogg et al. 2010).We follow here the prescription followed by the particle data group for averages of experimental results in this case, which consists in increasing the final uncertainty by a scale factor S = (χ 2 ν ) 1/2 (Workman et al. 2022).With this method we obtain a more conservative error: Optimal axion mass: m a = 1.90 +0.24 −0.21 × 10 −23 eV, corresponding to an uncertainty of 0.052 dex.For completeness we have also obtained an alternative determination of this optimal mass by performing a global simultaneous fit to the parameters involving the rotation curves of all ten galaxies in the sample and using the same m a .This procedure is more robust as it accounts for possible correlations and features in the fit and posterior distributions.Nonetheless, we obtain the same central value although with a relatively smaller uncertainty.The full set of results obtained from the fits using Eq. ( 23) are publicly available under the category Benchmark analysis in the FuzzyDM GitHub repository.
On the other hand, the values of m a determined from the different galaxies do not seem to be randomly scattered around the optimal value but rather following a negative correlation with the stellar mass.Specifically, the lower the M ⋆ value, the higher these m a values tend to be.This correlation is visible in the left panel of Fig. 4, where we have ordered the values of m a according to increasing values of the stellar mass of the galaxies.Since m a is a fundamental constant in the model, this suggests that there are astrophysical processes influencing the observed cores beyond the structure predicted by a FDM soliton with the optimal mass in Eq. ( 23).We can try to quantify the statistical significance of the correlation by introducing an empirical power law between the two quantities, m a = α ⋆ M β ⋆ ⋆ .We then performed a linear fit of α ′ ⋆ = log 10 α ⋆ and β ⋆ to the 10 (log 10 M ⋆ , log 10 m a ) Fig. 5. Soliton scaling relations compared to the results of the fits to LTs, where each point is colored according to the value of M ⋆ of the corresponding galaxy.Dot-dashed lines and gray bands represent the mean and the 68% and 95% CL regions for the correlation expected from FDM for the average axion mass in Eq. ( 23).The solid line and surrounding bands represent the maximum posterior and 68% and 95% CL regions of the fit to a power law of the results extracted from the individual galaxies.Left: Results in the r c − ρ c plane compared to Eq. ( 8) and where we also show the results obtained independently for a different sample of galaxies by using a Burkert profile (Rodrigues et al. 2017;Deng et al. 2018).Right: Results in the r c − M c plane compared to Eq. ( 8).
data points, excluding the three rogue galaxies we have studied.For this, we use MCMCs and also include a possible source of unknown systematic uncertainties in the data points as described in Sect.8 of Hogg et al. (2010).The results of this fit are shown in the right panel of Fig. 4. We find that β ⋆ = −0.62 +0.21  −0.25 , corresponding to evidence at a ∼ 3.3σ level (with respect to β ⋆ = 0) of a correlation between the values of m a determined from the rotation curves of the galaxies and the relevant M ⋆ .We describe the implementation of these fits and the interpretation of the CL regions in Appendix A. As shown in Fig. 4 the two rogue galaxies NGC 6882 and DDO 133 seem to also follow the correlation, while DDO 50 is an outlier.We discuss the possible implications of this correlation in more detail in Sect. 5.
Another perspective on the consistency of the inferred axion masses can be obtained by testing the soliton scaling relations Eqs. ( 7) and (8).To quantify the agreement of the model with observations, we fit our results of ρ c and M c to power laws in r c with exponents β ρ and β M , respectively.These fits are done again with MCMCs, including possible systematic uncertainties from Hogg et al. (2010), and they are described in Appendix A. In the left panel of Fig. 5, we present our results in the r c − ρ c plane.The data show that the central densities fall with the size of the core, although at a pace much slower than predicted by FDM.Indeed, our fit to the data yields β ρ = −1.18+0.22  −0.30 , disagreeing with the prediction of FDM, β ρ = −4, with a significance larger than 3σ.It is remarkable that our result is compatible with β ρ ≈ −1.3, obtained by Deng et al. (2018) using a different sample of galaxies analyzed with a Burkert profile Rodrigues et al. (2017).Our results also agree within 1σ with β ρ ≈ −1, which is an empirical result derived from the fact that the so-called "column density," Σ 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), andDi Paolo et al. (2019)).In fact, we obtain Σ 0 = 51 +12 −9 M ⊙ pc −2 which undershoots slightly the empirical value Σ 0 = 75 +55 −45 M ⊙ pc −2 (Burkert 2020).
In the right panel of Fig. 5, we show our results in the r c − M c plane.In this case the data manifest a correlation which is almost orthogonal to the FDM prediction: The cores are heavier, the bigger (not the smaller) they are.Quantitatively we find β M = 0.45 +0.10 −0.08 which (at face value) disagrees with the FDM value β M = −1 with a significance larger than 5σ.In addition, and as shown by the color code used in the data points of Fig. 5, r c , and M c are also correlated with M ⋆ , as the heavier cores contain also more stars.This suggests again that there is a non-trivial interplay between the structure of the DM cores at the center of the dwarf galaxies in LTs and baryonic physics.

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 top-left panel of Fig. 6 our results in the M h − M c plane compared to the core-halo relation predicted from simulations in Eqs. ( 14) and ( 15).This relation is imposed as a prior in the fits and it provides the most important constraint on the halo masses for most of the galaxies.This is why the uncertainties of M h span almost the entire allowed range for a given M c .Next, we show on the top-right panel our results in the M h − c plane comparing them to the predictions of CDM and FDM.The uncertainties on the concentration parameter c obtained from the fits are large although, as discussed in Sect.4.1, this is constrained from above by the central densities reached by the soliton.This reduces systematically the values of c obtained from the data as compared to the predictions of CDM and agreeing with the expectations of FDM.
The most important test that can be performed on the properties of the FDM halos with our results is related to the abundance of dwarf galaxies.As discussed in Sect.2, FDM predicts a suppression of structure formation at small scales that translates into a low-mass cutoff in the HMF.One classical method to test observationally the HMF is via abundance matching methods with the observed stellar mass function (SMF) (Frenk et al. 1988;Vale & Ostriker 2004;Moster et al. 2010;Behroozi et al. 2013).Given the non-monotonic behavior of the HMF in Eq. ( 18), we follow a method inspired by Vale & Ostriker (2004), whereby we find a M ⋆ − M h stellar-to-halo mass function by matching the cumula-Article number, page 9 of 22 A&A proofs: manuscript no.main

N
LGV FDM FDM (Schive + 16) LGV CDM (Brook + 14) 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.Gray-colored points are galaxies for which both median and maximum posterior values of r t are within the observed range, while gold-colored are soliton-dominated by virtue of not satisfying that condition.Uncertainties in the data points describe the credible intervals at 68% (solid) and 95% (dashed) CL.Uncertainties in the theoretical bands refer to the 1 and 2σ scatter in the predictions.Top-left: M c − M h relation compared to the region allowed by Eq. ( 15) and the prediction in Eq. ( 14 16) and to the prediction in CDM, Eq. ( 17) from Dutton & Macciò (2014) with σ = 0.16 dex (gray dashed line).Bottom-left: Stellar-to-halo mass relation from abundance matching compared to the FDM prediction using Eq. ( 18), and the CDM predictions from Behroozi et al. (2019) (dashed gray line) and Brook et al. (2014) (dot-dashed line) with a reference σ = 0.30 dex scatter in both cases.Bottom-right: Expected number of halos with M h ≲ 10 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 dot-dashed line and gray bands represent the mean and the 68% and 95% CL regions for the axion mass in Eq. ( 23).
tives of the SMF and HMF, In the bottom-left panel of Fig. 6 we show the results of abundance matching using the stellar-to-halo mass function reported in Behroozi et al. (2019) for CDM and corrected for FDM by Eq. ( 18) and the optimal axion mass in Eq. ( 23).We also add the prediction for CDM reported in Brook et al. (2014), which implements a HMF specific for the environment of the Local Group.The FDM prediction for the stellar-to-halo mass function features a flattening starting at values as high as M h ≈ 5 × 10 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 stellar-to-halo mass function predicted by CDM (Behroozi et al. 2019;Brook et al. 2014).
A more transparent way to illustrate this tension between FDM and observations is by directly comparing the expected number of halos that are predicted to exist in a specified local volume to the corresponding number of observed galaxies.Or, inversely, we can translate the abundance of observed galaxies to a lower bound on the FDM mass that is necessary to make the existence of such galaxies statistically plausible in the model.We take as a reference the simulations of galaxy formation in the Local Group within the ΛCDM that are reported in Brook et al. (2014)   The blue dot-dashed line corresponds to the same FDM mass with the inclusion of baryons, while the green dashed line to the mass needed to replicate the same density profile as the original soliton.Lastly, the orange dotted line corresponds to a soliton without baryons fixed at the value of the global optimal mass from Eq. ( 23).
(LGV), which is the one we use.In the bottom-right panel of Fig. 6, we show the number of FDM halos with M h ≲ 10 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 Sheth-Tormen prediction for the HMF of CDM (Sheth & Tormen 1999) and we normalize it to reproduce the results of Brook et al. (2014) for the LGV (see Appendix B).For reference, we also show the abundance that is obtained using the HMF originally reported in Schive et al. (2016).
The predicted abundances are to be compared with the three galaxies contained in the LGV that we have found in our analysis.More generally, we should also compare these predictions with the properties of 116 such galaxies that can be extracted from the most recent census Putman et al. (2021). 8sing 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 Milky-Way's halo using the subhalo mass function for FDM recently reported in Du (2018); Du et al. (2017).Using this, and focusing on the 59 galaxies from Putman et al. (2021) that lie within 400 kpc of the Milky Way (following Marsh & Niemeyer (2019) and which is roughly twice its virial radius r MW vir ≈ 200 kpc (Bland-Hawthorn & Gerhard 2016)), we obtain the stronger lower bound m a ≳ 7.3 × 10 −22 eV.

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. (2019Mocz et al. ( , 2020)); Veltmaat et al. (2020), finding that the inclusion of baryons in FDM halos does not impede their thermalization and formation of quasi-stationary (modified) solitons at their centers.This motivates the study of this solution even in the presence of baryonic effects.
A simple and minimal approach is to solve for the modified soliton solution under the assumption of a static background baryonic potential, which can be obtained in a very similar way to that of the original soliton.The resulting modified soliton solution has been studied in Bar et al. (2018) (and subsequently in Bar et al. (2019b)), which was found to accurately match simulations incorporating gravitational interactions with baryons in Chan et al. (2018).In order to proceed, we closely follow the approach of Bar et al. (2018).The methodology is very similar to that of the ordinary soliton in Eq. ( 1), except that in this case we introduce the baryonic density, ρ b , solving for the total (baryonic + DM) potential Φ tot and solve for a given central density ρ c .The resulting modified SP equation is given by We can then apply the same quasi-stationary ansatz as in Eq. ( 2), along with the boundary condition χ ′ (0) = 0, and proceed accordingly with the solution to the eigenvalue problem to obtain the modified soliton solution; that is, we solve for the unique, constantly decreasing solution with χ(r → ∞) = 0. We note that ρ b is taken to be static in time and spherically symmetric.Nonetheless, the validity of the spherical approximation was corroborated for DM-dominated galaxies in Bar et al. (2019b), where axially-symmetric baryonic distributions were Notes.m a is the median value of the original FDM mass found in Table 2 with corresponding errors, while ∆m a is the correction in m a needed to reproduce the original soliton without baryons, that is, the difference between the new FDM mass and the original one.
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 self-consistency of the model with baryons is to solve for the modified soliton for each galaxy fixing the axion mass and central density to the corresponding central values in Table 2.This leads to solitons more compact than the original ones, in agreement with Bar et al. (2018) and the results found in simulations (Mocz et al. 2019;Veltmaat et al. 2020).However, the differences between these solitons and their DM-only counterparts are less pronounced than in these simulations mainly because the relative proportion of baryons compared to DM is significantly lower in the LTs galaxies.This is illustrated in Fig. 7, where we show the effect of the gas and stars in the density distribution of the soliton for DDO 154 and NGC 2366.
We confirm the findings of Bar et al. (2018) that solitons formed including the baryonic effects closely follow the same universal distribution as the DM-only soliton, for instance, the one approximated by Eq. ( 9).This does not imply necessarily that the two solitons (DM-only and with baryons) are the same but rather that their difference is captured by a modification of the relation between ρ c and r c for a given m a (or of any other scaling relations such as r c and M c ); or, conversely, that ρ c and r c parametrize the same density profile regardless of whether baryons are included or not but where each case corresponds to a different value of m a .In particular, if baryons are included for a given empirical density profile, a lower m a should be obtained to compensate for the deepening of the gravitational potential.
Therefore, we can implement the inclusion of baryons in our analysis of the rotation curves of LTs by including a systematic negative shift in the FDM mass needed to replicate the original profile favored by the fits.In Table 3 we list the corrections to the value of m a (with respect to the results listed in Table 2) that is required to fit the rotation curves including the effects of the baryons.We find that the correction needed is typically within the original 68% CL errors, with the only exception of DDO 154 and NGC 2366 (both of which appear in Fig. 7) and NGC 6822 in case of the rogues.Based on these results, we conclude that the effect of including the background baryonic potential introduces a systematic drop of ≲ 15% in the axion masses which should be included as a correction in the determination of m a in Eq. ( 23).We note that this effect does not explain the empirical dependence of m a with M ⋆ that we found in our sample and does not help understanding the discrepancies of the data with the soliton scaling relations discussed in Sect.4.2.Furthermore, this effect of the baryons only exacerbates the tension found in the previous Sect.4.3 in relation to the excessive suppression of the HMF at small scales.
In connection to the baryonic effects, we also investigate the inner-slope parameter of the DM distribution of the LTs galaxies as a function of their mass.This quantity characterizes the flatness of the central DM distribution and it is a figure of merit for many studies of baryonic feedback within ΛCDM.We follow Di Cintio et al. ( 2014); Tollet et al. (2016) and define this parameter as the logarithmic slope of the density ρ(r) ∝ r α fitted to the data between 1 − 2% of the virial radius.In Fig. 8 we show the values of α for the LTs sample as a function of M ⋆ and compared to the prediction from simulations in ΛCDM (Di Cintio et al. 2014;Tollet et al. 2016).The M ⋆ range covered by the LTs galaxies overlaps with the region where baryonic feedback starts becoming relevant, M ⋆ ≳ 10 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.

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;Robles et al. 2019;Bar et al. 2022;Chan & Fai Yeung 2021;Street et al. 2022;Khelashvili et al. 2022;Dave & Goswami 2023), while Bernal et al. (2018) also studied a set of low-surface brightness (LSB) galaxies.In terms of methodology, our analysis presents some similarities with some of these works: In their work, Bernal et al. (2018) focuses on χ 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 soliton-halo scaling relations are not investigated nor the cosmological implications of the model for this mass.
In their work, Bar et al. (2018) focuses on the physical interpretation and phenomenological implications of the core-halo relation.In particular, the appearance of a "double-bump" structure in the rotation curves is predicted and then used to set an exclusion of FDM in the range of m a ∼ 10 −22 − 10 −21 eV range.This structure is not a generic feature appearing in our analyses due to the broader range of core-halo relations that we use, as described in detail in Appendix C. In their work, Bar et al. (2022) derived, for each m a , upper bounds on the mass of the soliton using a conservative method and the 175 galaxies in SPARC.However, this reference does not perform a full-fledged fit to the data and uses exclusively the core-halo relation as a exclusion criterion for FDM in the range m a ∼ 10 −23 − 10 −20 eV.In their work, Khelashvili et al. (2022) performs a Bayesian analysis of the 175 galaxies in SPARC and studies the statistical comparison between FDM and different CDM models.This yields that FDM is preferred by the data although no single value of m a gives a good fit to all galaxies.Moreover, the core-halo and soliton scaling relations are tested.
FDM has been studied in other different cosmological and astrophysical contexts and the derived limits on the axion mass have been recently compiled in Ferreira (2021).In the following paragraphs we include some of the cosmological bounds available in the literature.
The small-scale suppression given by FDM would impact the angular scale of the CMB acoustic peaks and anisotropies, as proposed in Hlozek et al. (2015Hlozek et al. ( , 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 sub-Mpc scales, where the FDM suppression arises (Niemeyer 2020).With high-resolution observations of the spectra on the Lyman-α forest, Nori et al. ( 2019 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 (DES-Y1) data and the Planck cosmic microwave background anisotropies information, Dentler et al. (2022) has reported a search of the effects caused by FDM in the cosmic shear.This sets m a > 10 −23 eV at 95% CL.
The abundances of lensed ultra-faint galaxies at high redshift provide 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.
The FDM's wave behavior would heat up stellar objects in galaxies (Bar-Or et al. 2019).In case of the Milky Way this effect can be search for through the velocity dispersion of stars in the solar system vicinity, leading to the lower bound m a > 6 × 10 −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 Milky-Way disc kinematics.
FDM affects dynamical friction and the effects can be searched for with the analysis of orbits of globular clusters in galaxies (Hui et al. 2017;Bar et al. 2021).In fact, FDM can help solving the so-called Fornax timing puzzle Tremaine (1976) as long as m a > 10 −21 eV (Lancaster et al. 2020).
Stellar streams are excellent tracers of the gravitational potential (Carballo-Bello et al. 2014) and can be used to test predictions of the sub-halo mass function in different DM models.The fluctuations in the stellar streams lead to imposing the lower limit m a > 2 × 10 −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 low-mass 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 Way-satellite dSphs, it is possible to study the density profile predicted by FDM.By considering a soliton core and kinematic data of the classical dSphs, Chen et al. (2017) obtained the value for m a = 1.79 +0.35  −0.33 × 10 −22 eV (at 2σ).Comparing to simulated data on the line-of-sight velocity dispersion (σ LOS ) of Fornax and Sculptor dSphs, González-Morales et al. (2017) found a tighter upper bound m a < 4 × 10 −23 eV at 97.5% CL.
The existence of the sub-halo to host Eridanus II low-surface brightness dwarf galaxy and the survival of its star cluster can also constrain the FDM mass (Bar-Or et al. 2019;Brandt 2016).As shown in Marsh & Niemeyer (2019), the sub-halo mass function in the Milky Way and the existence of Eridanus II implies m a ≳ 8 × 10 −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 ultra-faint dwarf (UFD) galaxies can be used to estimate the effects of wave interference produced by FDM.In particular, in order for the FDM's gravitational heating to be weak enough to explain the observed velocity dispersions of 2.5 − 3 km/s and sizes of 24 − 40 pc in Segue 1 and 2 (Simon et al. (2011); Kirby et al. (2013)), Dalal & Kravtsov (2022) derived a lower bound of m a > 3 × 10 −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 m a = 1.1 +8.3  −0.7 × 10 −19 eV, which is consistent with the results obtained directly from the MUSE-Faint survey (Zoutendijk et al. 2021).Finally, Calabrese & Spergel (2016) found that FDM particle with a mass of m a ∼ 3.7 − 5.6 × 10 −22 eV is in agreement with the data of half-light mass of Draco II and Triangulum II; although Safarzadeh & Spergel (2020) later found an upper limit m a ≲ 10 −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 dwarf-galaxy counts in the LGV.Nonetheless, before closing this section, we would like to point out that many of the analyses leading to these constraints depend on astrophysical modeling assumptions and are, thus, subject to significant uncertainties.We refer, for instance, to Ferreira (2021) and, in particular, to Sect. 4.3 of Chiang et al. (2022) for a recent critical reappraisal of some of these constraints.

Conclusions
In this work we have investigated and tested the predictions of FDM at galactic scales using a homogeneous and robust sample of high-resolution rotation curves from the LITTLE THINGS in 3D (or LTs for short in our paper) catalog.This comprises a collection of isolated, dark matter-dominated dwarf-irregular galaxies that provides an optimal benchmark for cosmological studies.The methodological basis of our study is a statistical framework based on the χ 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 core-halo relation The constraint from HMF in the local group with the halos of nearby dwarf irregulars (dIrrs) from the LTs catalog is displayed with dark turquoise bars.On the dIrrs bar, we also show the optimal mass m a = 1.9 +0.5 −0.4 × 10 −23 eV (uncertainties at 2σ CL).Our constraints are compared to other cosmological and astrophysical limits (see main text for details).
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 m a = 1.90 +0.24 −0.21 × 10 −23 eV.-However, we find that the individual determinations of m a are not scattered randomly around the average, but follow instead a mild correlation with the stellar mass of the galaxy with a significance of ∼ 3.3σ.Moreover, displaying the results of our fits in the r c − ρ c and r c − M c planes, we find scaling relations in the data that are at odds with the predictions of FDM.In case of the r c − M c relation, the data show a scaling that is almost orthogonal to the one predicted by FDM with a significance that (at face value) is ≳ 5σ.From this, we conclude that the cores we find in the LTs catalog do not have the properties of the standard FDM solitons.-In parallel, we investigated various cosmological implications related to the FDM halos for an axion mass m a ≈ 2 × 10 −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 galaxy-dependent negative correction to the axion mass.However, these corrections are typically too small to alleviate any of the problems of FDM with data that we found.We also briefly discuss the possible role of baryonic feedback and pointed out that our galaxy sample is in a region of stellar masses that overlaps with the one where these effects can be significant in CDM.
In summary, our analysis poses a serious challenge to FDM as a contending hypothesis for solving the core-cusp problem or any of the other small-scale issues related to ΛCDM.The tension that we have found between the favored FDM mass and the abundance of halos is yet another example of the problems of the model with cosmology, which is reminiscent of the so-called "catch-22 problem" found in warm DM (Macciò et al. 2012) and renewed recently in the context of FDM (Mocz et al. 2023).This arises when one wishes to replicate the core-like structures in dwarf galaxies by invoking a particle mass that would imply a suppression of small-scale structure that is much too strong to be consistent with observations.Baryonic physics within FDM could help relieving some of these tensions, although it seems unlikely in light of our findings and the results from pioneering FDM simulations with realistic baryonic feedback.An alternative approach would be extending the dark sector with, for instance, self-interactions, as recently explored in Mocz et al. (2023) as a possible mechanism to enhance small-scale structure formation, or by demanding that only a fraction of DM is FDM, as in, for instance, Kobayashi et al. (2017).Future work is expected to shed more light on these topics.On the observational side, a reanalysis of HI data from other isolated dwarf irregular galaxies, using the same methods as in the LTs catalog, would increase the statistical power of cosmological analyses, such as the one presented in this paper.of the asymmetric drift correction V a as: 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 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 The red dashed lines are the median of the posterior distribution for the total velocity and the colored band is the corresponding uncertainty at 68% and 95% CL.High-resolution figures can be found in the FuzzyDM GitHub repository.Notes.The last two columns include the reduced chi-square of the maximum posteriors for the NFW fit compared to the ones obtained for the FDM fits shown in Table 2.
algorithm, based on MCMC, uses the same attributes as those described in Sect.3.2.
An additional correction for the log-likelihood is necessary since now the errors in Eqs.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 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).10This implies an important physical difference, since the potential (and implicitly the virialized energy) of a uniform spherical distribution is not equivalent to that of a self-gravitating NFW profile.It appears that such energy distribution, in combination with Eq. (C.1), is the median of the results of the simulations, as seen in Figs. 2 and 4 of Schive et al. (2014b), although Eq. (C.6) is compatible with an upward scatter of the data.
This relation is also contained within the Eq. ( 15) from Chan et al. (2022) that is used in our analysis.However, this region of core-halo relations, leading to heavy solitons, is not favored by our fits for relatively heavy axion masses, m a ≳ 10 −22 eV.In particular, the rotation curve produced by a soliton has a maximum v c,max at a radius r c,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. 1 .
Fig. 1.Solutions to Eq. (4) for the ground state and first two excitations, where the scale invariance was used to set χ(0) = 1.

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

Fig. 8 .
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.

Fig. 9 .
Fig.9.Bounds from cosmology and astrophysics on the axion mass.The constraint from HMF in the local group with the halos of nearby dwarf irregulars (dIrrs) from the LTs catalog is displayed with dark turquoise bars.On the dIrrs bar, we also show the optimal mass m a = 1.9 +0.5 −0.4 × 10 −23 eV (uncertainties at 2σ CL).Our constraints are compared to other cosmological and astrophysical limits (see main text for details).
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 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 Table1and 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 The red dashed lines are the median of the posterior distribution for the total velocity and the colored band is the corresponding uncertainty at 68% and 95% CL.High-resolution figures can be found in the FuzzyDM GitHub repository.

Fig
Fig. A.3.Rotation curves in FDM for the inclination rogues.We represent the same information as in Fig. A.2.
(A.1) and (A.2) are i-dependent:ln L(θ FDM , M ⋆ , i) = ln L ′ (θ FDM , M ⋆ , i) − 1 2 j ln(δV ′2 c, j ), (A.3)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), ϕ 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 core-halo relation reported byBar et al. (2018), of c that varies between 0.9 ≲ f (c) ≲ 1.1 for 5 ≤ c ≤ 30.The core-halo relation in Eq. (C.4) is analogous to the one in Eq. (14) but it is defined using a different convention for the halo mass.It is straightforward to transform the NFW conventions from a critical spherical overdensity definition of 200 to one of 100, which is close to the one used bySchive 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 M c ≈ 5.5 × 10 9 Fig. C.1.Comparison of the effect of imposing different core-halo relations in the rotation curve of DDO 154.Theoretical rotation curves are DM-only (i.e., excluding baryonic components) for visual clarity.

Table 1 .
Priors used to fit the DM and baryon contributions to the rotation curve data from the LTs catalog.
using 1.8 Mpc for the radius of the Local Group Volume

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