Galactic Rotation Curves vs. Ultra-Light Dark Matter: Implications of the Soliton -- Host Halo Relation

Bosonic ultra-light dark matter (ULDM) would form cored density distributions at the center of galaxies. These cores, seen in numerical simulations, admit analytic description as the lowest energy bound state solution ("soliton") of the Schroedinger-Poisson equations. Numerical simulations of ULDM galactic halos found empirical scaling relations between the mass of the large-scale host halo and the mass of the central soliton. We discuss how the simulation results of different groups can be understood in terms of the basic properties of the soliton. Importantly, simulations imply that the energy per unit mass in the soliton and in the virialised host halo should be approximately equal. This relation lends itself to observational tests, because it predicts that the peak circular velocity, measured for the host halo in the outskirts of the galaxy, should approximately repeat itself in the central region. Contrasting this prediction to the measured rotation curves of well-resolved near-by galaxies, we show that ULDM in the mass range $m\sim (10^{-22}\div 10^{-21})$ eV, which has been invoked as a possible solution to the small-scale puzzles of $\Lambda$CDM, is in tension with the data. We suggest that a dedicated analysis of the Milky Way inner gravitational potential could probe ULDM up to $m\lesssim 10^{-19}$ eV.

Axion-like particles with exponentially suppressed masses arise in the context of string theory [12,[14][15][16]. They are produced by misalignment and for mass in the range m ∼ (10 −22 ÷ 10 −18 ) eV their cosmic abundance could naturally match the observed DM density. ULDM with mass m ∼ (10 −22 ÷ 10 −21 ) eV is particularly motivated due to several puzzles, facing the standard WIMP paradigm on galactic scales (see, e.g. [12,17] for recent reviews). This range of m is in some tension with the matter power spectrum, inferred from Ly-α forest analyses [18][19][20][21][22], which yields a bound m 10 −21 eV. Nevertheless, since the strongest Ly-α bound constraints come from the smallest and most non-linear scales, where systematic effects are challenging, we believe it prudent to seek additional methods to constrain the model for m 10 −22 eV.
Numerical simulations [6,7,10,11,13] show that ULDM forms cored density profiles in the inner region of galactic halos, roughly within the de Broglie wavelength, The core, referred to as "soliton" in the literature [4,5,8,9], corresponds to a coherent quasi-stationary solution of the ULDM equations of motion and its density profile can be derived analytically 1 . This match between analytic and numerical simulation results offers a unique opportunity to understand and extend the simulations, making the ULDM model potentially more predictive than the WIMP paradigm. For near-by dwarf galaxies (v ∼ 10 km/s) or Milky Way (MW)-like galaxies (v ∼ 100 km/s), the soliton core could be resolved with current observational tools, offering a test of the model.
With some abuse of the word, we refer to it here as analytical to emphasize that the solution can be found by integrating a simple differential equation in a procedure that takes < 1 sec on a standard laptop.
complex structure, with the central soliton transiting into a large-scale host halo composed of an incoherent superposition of multiple scalar field wavepackets. Consequently, other works [8,[24][25][26] analysed galactic profiles using a soliton+host halo description.
A key question in comparing the soliton+host halo model to kinematical data, is how to model the transition between the central soliton and the host halo. Previous phenomenological analyses [8,[24][25][26] defined separate free parameters for the host halo and for the soliton, for each individual galaxy in the sample. On the numerical simulation side, several studies focused on finding a soliton-host halo relation [6,10,11,13]. The point in finding a soliton-host halo relation is, of course, that it could tie together the behaviour of ULDM in the largescale host halo to predict the central core with fewer free parameters.
In this work we consider the soliton-host halo relations, reported by different numerical simulation groups [6,7,10,11,13]. Our first observation is that properties of the analytic soliton solution can provide important insight on the numerical results. We show that: (i) the solitonhost halo relation, reported in [13], essentially attributes the total energy (kinetic+gravitational) of the halo to the dominant soliton. This energetic dominance of the soliton is unlikely to hold for realistic galaxies above a certain size; (ii) the soliton-host halo relation, reported in [7], essentially equates the energy per unit mass in the soliton to the energy per unit mass in the virialised halo. This relation can apply to real galaxies.
Assuming the soliton-host halo relation of [7], we show that it leads to a prediction: the peak circular velocity characterising the host halo on large scales (few kpc for typical (10 9 ÷ 10 10 ) M galaxies) should repeat itself in the core on small scales ( 1 kpc), insensitive to the details of the host halo density profile. This implies an observational constraint that can be tested without free parameters. Applying this test to high resolution rotation curves of late-type galaxies from [27,28], we find that ULDM in the mass range m ∼ (10 −22 ÷ 10 −21 ) eV is disfavored by the data. Baryonic physics is unlikely to cure the tension: the discrepancy between the predictions of ULDM (with the soliton-host halo relation) and the data is too large. In many galaxies in the sample we analyse, baryonic feedback would need to overcome a dark matter-to-baryon mass ratio of 10 : 1, to destroy the soliton. As a result, if the soliton-host halo relation of [6,7] holds for real systems, ULDM is disfavored below m ∼ 10 −21 eV and is unlikely to play a role in solving the small scale puzzles of ΛCDM.
We also discuss observational imprints of the soliton in big galaxies, such as the MW. In this case the shape of the soliton is modified due to the gravitational potential of baryonic matter and supermassive black hole (SMBH). Preliminary numerical simulations of ULDM halos with stars [29] indicate that in the presence of baryons, the soliton mass stays the same or increases compared to the soliton-host halo prediction of pure ULDM. If these results are confirmed, solitons formed by ULDM with masses m 10 −19 eV will produce order-one contribution into the mass budget of the inner MW. A dedicated analysis of inner MW kinematics -including simultaneous modeling of the baryonic mass and the soliton -could potentially test ULDM up to m 10 −19 eV.
The outline of this paper is as follows. In Sec. II we review some basic properties of the soliton. In Sec. III we discuss the soliton-host halo relations found in simulations, and show that these relations can be understood in terms of fundamental properties of the soliton and the halo. The simulations of [6,7] are considered in Sec. III A; they can be summarised by the statement (E/M )| soliton = (E/M )| halo , where E is the total energy (kinetic+gravitational, within the virial radius, for the halo) and M is the mass (again within the virial radius, for the halo). The simulations of [13] are considered in Sec. III B; their soliton-halo result can be summarised by the statement E| soliton = E| halo . We argue that this result is unlikely to represent realistic galaxies above a certain size.
In Sec. IV, adopting the soliton-host halo relation of [6,7], we work out its observational consequences. We show that in ULDM galaxies satisfying this relation, the rotation velocity in the inner core should be approximately as high as the peak rotation velocity in the outer part of the galaxy. In Sec. IV A we compare this analysis to numerical profiles taken directly from the published simulation results, finding good agreement. In Sec. IV B we compare this prediction to real galaxies from [27,28], finding tension for m ∼ (10 −22 ÷ 10 −21 ) eV. We show that the intrinsic scatter in the soliton-host halo relation does not resolve the discrepancy.
In Sec. V we calculate how baryonic effects could modify the soliton solution. In Sec. V A we consider a smooth fixed distribution of baryonic mass, deferring the case of a super-massive black hole to App. B. These results are applied to the MW in Sec. V B.
In Sec. VI we compare our results to previous literature and discuss some caveats to our analysis. Sec. VII contains a summary of our results and open questions.

II. SOLITON PROPERTIES: ANALYTIC CONSIDERATIONS
In this section we review the relevant properties of the soliton that help to understand results from numerical simulations.
We consider a real, massive, free 2 scalar field φ, satisfying the Klein-Gordon equation of motion and minimally coupled to gravity. In the non-relativistic regime it is convenient to decompose φ as with complex field ψ that varies slowly in space and time, such that |∇ψ| m|ψ| and |ψ| m|ψ|. The field ψ satisfies the Schroedinger-Poisson (SP) equations [32] i∂ t ψ = − 1 2m We look for a quasi-stationary phase-coherent solution, described by the ansatz 3 The ULDM mass density is The parameter γ is proportional to the ULDM energy per unit mass. Validity of the non-relativistic regime requires |γ| 1. Since we are looking for gravitationally bound configurations, γ < 0.
Assuming spherical symmetry and defining r = mx, the SP equations for χ and Φ are given by Finding the ground state solution amounts to solving Eqs. (7)(8) subject to χ(r → 0) = const, χ(r → ∞) = 0, with no nodes. Given the boundary value of χ at r → 0, the solution is found for a unique value of γ.
It is convenient to first solve Eqs. (7)(8) with the boundary condition χ(0) = 1. Let us call this auxiliary solution χ 1 (r), with γ 1 . A numerical calculation gives [4,5,8] and the solution is plotted in Fig. 1. The mass of the χ 1 soliton is Its core radius, defined as the radius where the mass density drops by a factor of 2 from its value at the origin, is FIG. 1. Profile of the "standard" χ1 soliton with λ = 1 (blue solid). We also show the corresponding gravitational potential (orange dashed) and circular velocity of a test particle (dotted green).
Other solutions of Eqs. (7)(8) can be obtained from χ 1 (r), Φ 1 (r) by a scale transformation. That is, the functions χ λ (r), Φ λ (r), together with the eigenvalue γ λ , given by also satisfy Eqs. (7)(8) with correct boundary conditions for any λ > 0. The soliton mass and core radius for χ λ are A mnemonic for the numerical value of λ is given by The product of the soliton mass and core radius is independent of λ, Formally, solutions exist for any positive value of λ and hence for any soliton mass. However, if we select λ 1 we reach |γ λ | > 1, outside of the regime of validity of the non-relativistic approximation. Thus, self-consistent solutions are limited to λ 1 and their eigenvalue |γ λ | = λ 2 |γ 1 | 1, consistent with the non-relativistic approximation.
The energy in an arbitrary non-relativistic ULDM configuration is
The first point to note is that soliton configurations, in a form close to the idealised form discussed in Sec. II, actually occur dynamically in the central region of the halo in the numerical simulations 4 . In Fig. 2 we collect representative density profiles from Ref. [6] (blue), Ref. [13] (orange), and Ref. [10] (green). We refer to those papers for more details on the specific set-ups in each simulation. To make Fig. 2, in each case, we find the λ parameter that takes the numerical result into the χ 1 soliton, rescale the numerical result accordingly and present it in comparison with the analytic χ 2 1 (r) profile.  While different groups agree that solitons form in the centers of halos, they do not appear to agree on the matching between the inner soliton profile and the host halo. Refs. [6,7] and Ref. [13] reported scaling relations between the central soliton and the host halo. As we show below, the scaling relations found by both groups are connected to properties of a single, isolated, self-gravitating soliton (part of these observations were made in [10,11]).
A. Soliton vs. host halo: the simulations of Ref. [6,7] At cosmological redshift z = 0, the numerical simulations of [6,7] yield approximately NFW-like halos which transit, in the central region, into a core with core radius and mass density   4 The first simulations of cosmological ULDM galaxies [33] did not have sufficient resolution to resolve the central core.
where M h is the virial mass of the host halo. As noted in [6][7][8], Eqs. (29)(30) are an excellent numerical fit for a soliton χ λ . The mass of this soliton is However, this is just Eq. (23), if we replace the halo energy E h by the energy of the soliton. Because the central density profile found in [13] was a χ λ soliton, to a good approximation, it must be the case that the total energy of the halo in the simulations of [13] was dominated by the central soliton contribution. This situation is unlikely to hold for realistic cosmological host halos with M h significantly above M h,min . How could this have happened? The initial conditions in the simulations of [13] were a collection of N solitons, which were then allowed to merge. It appears that these initial conditions were constructed such that one initial state soliton -the soliton of initially largest mass -grew to absorb the entire energy of the system. Differently from Ref. [7] that considered initial conditions of N identical initial solitons, the simulations of [13] initiated their N solitons with a random flat distribution in soliton radius. Such distribution would be skewed towards large soliton energy because E λ ∝ x −3 cλ . Considering the initial condition set-up as explained in [13], we find that the most massive initial state soliton typically needed to grow in mass by only a factor of 1.5 ÷ 2, to absorb the entire energy of the halo.
Note that energy dominance of the central soliton over the host halo, implied by Eq. (36), is not the same, of course, as equating the energy per unit mass of the soliton and the halo, implied by Eq. (35). Halos in [6,7] attained masses up to two orders of magnitude larger than the central soliton mass, meaning their halo energy was two orders of magnitude larger than the energy of the soliton.

C. Comments
As far as we can currently determine, Eq. (35) may indeed reflect a realistic soliton-host halo relation for large enough cosmological halos. In the following sections, we take a leap of faith and assume that the simulations of [6,7] produced the correct scaling relation. We stress that Eq. (35) is an empirical result, and was only tested in [6,7] for host halo masses ranging from ∼ 10 8 M to ∼ 10 11 M . Our key numerical analysis will concern systems in this range of mass.
We defer a theoretical study of the origin of Eq. (35) to future work. Here we give only a few comments. We stress that the discussion in the rest of this section does not affect any of our results.
For a soliton, E/M = γ/3. On the other hand, γm can be associated with the chemical potential of ULDM particles in the soliton (see e.g. [34] and references therein). This may appear to hint that Eq. (35) corresponds to thermodynamic equilibrium between the ULDM particles in the host halo and in the soliton. However, there is some evidence to the contrary from simulations.
Ref. [29] simulated ULDM, adding collisionless point particles ("stars"). The stars aggregated dynamically in a cuspy profile, resulting in a more massive soliton compared to the pure ULDM simulations [6,7] with a given host halo mass. Testing the reversibility of the system, Ref. [29] adiabatically "turned off" the stars after the initial system virialised. When eliminating the stars, the soliton+halo system did not relax back to Eq. (35). Instead, the excess ULDM mass that was contained in the soliton in the presence of stars remained captured in the soliton, and did not return to the host halo. The final state of the system was not described by Eq. (35): the soliton ended up containing larger (negative) E/M than the halo, and larger mass compared with Eq. (31).

IV. SOLITON-HOST HALO RELATION AND GALACTIC ROTATION CURVES
As we have seen, the soliton-host halo relation found in the simulations of [6,7] can be summarised by Eq. (35), equating the energy per unit mass of the virialised host halo to that in the soliton component. For a virialised system, the energy per unit mass maps to kinetic energy density: in particular, the characteristic circular velocity (or, up to an O(1) geometrical factor, the velocity dispersion) of test particles in the halo and in the soliton should match. The peak circular velocity of the soliton, given by Eqs. (27)(28), occurs deep in the inner part, x < 1 kpc, of the galaxy; while the peak circular velocity of an NFW-like halo occurs far out at x ∼ 2 R s , with R s the NFW characteristic radius, of order 10 kpc for a MW-like galaxy. Thus, if the scaling derived from the simulations of [6,7] is correct, ULDM predicts that the peak rotation velocity in the outskirts of a halo should approximately repeat itself in the deep inner region. We now discuss this result quantitatively.
Consider a halo with an NFW density profile where The profile has two parameters: the radius R s and the concentration parameter c = R 200 /R s , where R 200 is the radius where the average density of the halo equals 200 times the cosmological critical density, roughly indicating the virial radius of the halo. The gravitational potential of the halo is Near the origin, We can estimate the energy per unit mass of the virialised halo by This gives Typical values of the concentration parameter are in the range c ∼ 5 ÷ 30 [35]. In this range,c varies betweeñ c ∼ 0.55 ÷ 0.35, respectively. (For reference, fits of the MW outer rotation curve give c ∼ 10 ÷ 20 [36].) Plugging Eq. (42) into the soliton-host halo relation Eq. (35), the scaling parameter λ is fixed as Eq. (45) depends weakly on the NFW concentration parameter, via the factor f (c) that varies in the range 0.9 ÷ 1.1 for c = 5 ÷ 30. It agrees parametrically with the simulation result, Eq. (31) (including the redshift dependence, which we have suppressed in Eq. (31)). It also agrees quantitatively to about 20%; to see this, we need to account for the slightly different definition of the halo mass M h , used in [7], and our M 200 . We do this comparison in App. A. Consider the rotation velocity curve of an ULDM galaxy satisfying Eq. (35). The NFW rotation curve is given by This halo rotation curve peaks at x ≈ 2.16 R s with a peak value On the other hand, in the inner galaxy x R s , the circular velocity due to the soliton peaks to a local maximum of (48) where we used Eq. (44) to fix λ and Eq. (28) to relate it to maxV circ,λ .
As anticipated in the beginning of this section, Eq. (35) predicts approximately equal peak circular velocities for the inner soliton component and for the host halo,  While maxV circ,λ and the approximate equality Eq. (49) are m-independent, the soliton peak velocity occurs in an m-dependent location, In Fig. 4, to define the rotation velocity for the total system, we set the ULDM mass density for the total system to be ρ(x) = max {ρ λ (x), ρ N F W (x)}, calculate the resulting mass profile M (x), and use spherical symmetry to find V circ (x) = G M (x)/x. This prescription for matching between the soliton and NFW parts is ad-hoc and only roughly consistent with the simulations of [6,7]. The true transition region between the NFW part and the soliton part probably deviates from the pure NFW form. Ref. [37] considered this transition region and concluded that the density profile in this region should follow approximately ρ ∼ x − 5 3 , steeper than the usual inner NFW form ρ ∼ x −1 . This would affect the detailed shape of the rotation curve in the intermediate region between the two peaks, but not our general results.
Eq. (49) was derived for an NFW host halo, but it is the manifestation of Eq. (35) that is not tied to a particular parametrisation of the halo profile. Building on Eq. (35), we expect in general that for DM-dominated galaxies, the soliton peak circular velocity should roughly equal the peak circular velocity in the host halo. The NFW example demonstrates that details of the host halo profile affect this result at the 10% level or so.
In the rest of this paper, when we refer to Eq. (49), we set the RHS to unity. Approximating the RHS of Eq. (49) by unity, and replacing maxV circ,h instead of maxV circ,λ in Eq. (28), the peak circular velocity of a host halo allows to predict the scale parameter λ and thus the soliton relevant for that host halo.

A. Comparison to numerical simulations
In Fig. 5 we compare our results to two soliton+halo configurations from the simulations of [6] and [29] (for [6], we take the largest halo, and for [29] we take the initial state of Case C). To calculate the soliton, we read maxV circ,h from the large-scale peak (at x ∼ 20 kpc) of the numerically extracted halo rotation curves (solid lines). Following Eq. (49), we use maxV circ,h instead of maxV circ,λ in Eq. (28), and read off the value of λ. The predicted soliton bump is shown in dashed lines. It gives the correct soliton peak rotation velocity to ∼ 20% accuracy in both cases. In Fig. 6 we show the velocity profiles for 11 simulated halos, calculated for 6 halos from [7] (solid lines) and 5 halos from [6] (dashed lines). The rotation curves are scaled and normalised such that x peak,λ = 1 and maxV circ,λ = 1 in each case. All of the halos satisfy 0.65 < maxV circ,λ maxV circ,h < 1.4; for later reference, the shaded band highlights the range 0.5 < maxV circ,λ maxV circ,h < 1.5.

B. Comparison to real galaxies
We now consider some observational consequences of our analysis. We choose to do so by examining the rotation curves of nearby disc galaxies with halo masses in the range covered by the simulations of [6,7], and above the minimal mass of an ULDM halo, see Eq. (33).
In Figs. 7-10 we show the rotation curves of four galaxies taken from Ref. [27] (see Ref. [28] for a recent rendering of these and many other rotation curves), for which high-resolution kinematical data is available. The observed rotation curves are compared to the soliton contribution predicted by Eq. (49) We emphasize that in using Eq. (49) to predict the soliton, we set the RHS of that equation to unity, and thus we ignore any details of the shape of the host halo. We also neglect corrections due to baryons. On the one hand, as we have learned from the NFW analysis, this prescription for deriving the soliton profile would suffer O(10%) corrections from the detailed halo shape. The baryonic (stellar and gas) contribution to the gravitational potential affects the halo velocity at a similar level: in Ref. [28], the baryonic contribution to the halo peak rotation velocity of UGC 1281, UGC 4325, and NGC 100 was estimated from photometric data to be well bellow the DM contribution 6 , V  [7] (solid lines) and [6] (dashed lines), scaled and normalised such that x peak,λ = 1 and maxV circ,λ = 1. Shaded band highlights the range 0.5 < maxV circ,λ maxV circ,h < 1.5. 6 For this estimate we use 3.6µm mass-to-light ratio Υ * = 0.5 M /L . this simple procedure relieves us from the need to fit for the virial mass or other details of the host halo. All that is needed is the peak halo rotation velocity, a directly observable quantity 7 . Eq. (49) (represented by the dashed line in Figs. 7-10) corresponds to the central value of the soliton-host halo relation. Ref. [6,7] showed a scatter of about a factor of two around Eq. (34) between simulated halos. This translates to a factor of two scatter in the soliton λ parameter (we have illustrated this scatter for a sub-sample of simulated halos in Fig. 6). In Figs. 7-10 we represent this scatter by a shaded band, showing the results when the λ parameter inferred from Eq. (49) is changed by a factor of 2.
It is important to check if scatter between different galaxies could explain the discrepancy, with the four galaxies in Figs. 7-10 being accidental outliers. To address this question, we analyse the 175 rotation curves contained in the SPARC data base [28]. This sample includes, in particular, the galaxies UGC 1281, UGC 4325, and NGC 100, shown in Figs. 7-10.
Our SPARC analysis is as follows. For each galaxy, we make a crude estimate of the halo mass contained within the observed rotation curve profile by M gal ∼ RV 2 /G, where R is the radial distance of the last data point in the rotation curve, and V the corresponding velocity. We keep only galaxies with 5 × 10 11 M > M gal > 5 × 10 8 m/10 −22 eV −3/2 M . We do this in order to limit ourselves to galaxy masses that are comfortably above the minimal halo mass (33), and not above the range simulated in [6,7]. Our results are not sensitive to the details of this mass cut. Next, for each galaxy we determine the observed maximal halo rotation velocity maxV circ,h , and use it to compute the soliton prediction from Eq. (49). To avoid confusion between halo peak velocity and soliton peak velocity, we search for the halo peak velocity restricting to radial distance x > 3 m/10 −22 eV −1 kpc. Galaxies with no data above x = 3 m/10 −22 eV −1 kpc are discarded.
Our results are not sensitive to this criterion; defining the halo cut anywhere at 1 m/10 −22 eV −1 kpc guarantees that such confusion is avoided. SPARC galaxies come with photometric data, allowing to model the baryonic contribution to the gravitational potential [28]. We use this information to limit baryonic effects on our analysis, and to explore the sensitivity of our results to baryonic corrections. For each galaxy, we estimate the baryonic contribution to the observed rotation velocity using the mass models of [28] with 3.6µm mass-to-light ratio Υ * = 0.5 M /L . Setting V Our first pass on the data includes only galaxies for which the predicted soliton is resolved, namely, x peak,λ from Eq. (50), with maxV circ,λ = maxV (DM) circ,h , lies within the rotation curve data. For these galaxies, we compute from data the ratio Here, V circ, obs (x peak,λ ) is the measured velocity at the expected soliton peak position.
The results of this first pass on the data are shown in Fig. 11. Red, blue, and green histograms show the result when imposing f bar2DM < 1, 0.5, 0.33, respectively. For m = 10 −22 eV, we find 45, 26, and 5 galaxies that pass the resolved soliton cut for f bar2DM < 1, 0.5, 0.33. For m = 10 −21 eV, only 4 galaxies pass the resolved soliton cut for f bar2DM < 1, and none for f bar2DM < 0.5, 0.33.
Including only galaxies with a resolved soliton causes us to loose many rotation curves with discriminatory power. For example, UGC 4325 drops out of the analysis for m = 10 −21 eV, though it clearly constrains the model, as seen from the lower panel of Fig. 8. To overcome this without complicating the analysis, we perform a second pass on the data. Here, we allow galaxies with unresolved soliton, as long as the innermost data point is located not farther than 3 × x peak,λ . We need to correct for the fact that the soliton peak velocity is outside of the measurement resolution. To do this, we modify our observable as where x min,data is the radius of the first data point. This correction is conservative, because it takes the fall-off of the soliton gravitational porential at at x > x peak,λ to be the same as for a point mass. In reality, the potential decays slower and the soliton-induced velocity decreases slower. Keeping this caveat in mind, Fig. 12  with f bar2DM < 1, 0.5, we find 48 and 16 galaxies with unresolved soliton, that can be added to the sample of Fig. 11. No galaxy is added for f bar2DM < 0.33. For m = 10 −21 eV, 16 and 5 galaxies are added with f bar2DM < 1, 0.5, and none for f bar2DM < 0.33.
In Figs. 11-12, vertical dashed line indicates the soliton-host halo prediction. The shaded region shows the range of the prediction, modifying the RHS of Eq. (49) between 0.5 − 1.5, consistent with the scatter seen in the simulations.
We conclude that the four galaxies in Figs. 7-10 are not outliers: they are representative of a systematic discrepancy, that would be difficult to attribute to the scatter seen in the simulations. If the soliton-host halo relation of [6,7] is correct, then ULDM in the mass range  We have limited our attention to the range m = (10 −22 ÷ 10 −21 ) eV, for which we believe the results are clear. We leave a detailed study of the precise exclusion range to future work. We note that for lower particle mass, m 10 −23 eV, the soliton contribution extends over much of the velocity profile of many of the SPARC galaxies, leaving little room for a host halo. This limit, where the galaxies are essentially composed of a single giant soliton, was considered in other works. We do not pursue it further, one reason being that this range of small m is in significant tension with Ly-α data [18,19].
For higher particle mass, m 10 −21 eV, the soliton peak is pushed deep into the inner 100 pc of the rotation curve. Although high resolution data (e.g. NGC 1560, Fig. 9) is sensitive to and disfavours this situation, a more careful analysis would be needed to draw a definitive conclusion. Note that in this range of m, ULDM ceases to offer a solution to the small-scale puzzles of ΛCDM (see, e.g., review in [12]).

A. A fixed distribution of baryons
We now consider the soliton solution, and some aspects of the soliton-host halo relation, with a coexisting distribution of baryonic mass. In this section we consider a smooth mass distribution, deferring the analysis of the effect of a super-massive black hole (SMBH) to App. B.
The analysis of baryonic effects is qualitatively important, and quantitatively relevant to ULDM in, e.g., MWlike galaxies. We will see that in order to make a significant impact, baryons need to constitute an O(1) fraction of the total mass in the soliton region. For many of the galaxies that we analysed in Sec. IV B, the predicted soliton contribution to the rotation velocity exceeds the observed velocities by a factor of 3-5, implying that the soliton over-predicts the mass in the central region of the galaxy by about an order of magnitude. This leaves little room for baryons (and, we think, baryonic feedback) to significantly affect the dynamics of the inner ULDM halo.
Proceeding to the analysis, we denote the gravitational potential due to baryons by Φ b (r). We assume that the distribution of baryonic mass is spherically symmetric 8 and dies off at infinity sufficiently fast, so that It remains convenient to solve the problem using boundary conditions with χ(0) = 1. Let us denote this solution (satisfying χ(0) = 1) by χ 1 (r; Φ b ), accompanied by the soliton potential Φ 1 (r; Φ b ). The general solution χ λ (r; Φ b ) satisfying Eqs. (53-54) with boundary condition χ λ (0; Φ b ) = λ 2 is given by Defining the soliton mass and energy in the presence of the baryons by where the ULDM energy is (59) 8 The assumption of spherical symmetry is not realistic; some of the galaxies considered in Sec. IV B may in fact exhibit maximal discs. Nevertheless, we expect our simplified analysis to give us correct order of magnitude estimates, and defer the analysis of non-spherical baryonic+soliton+halo systems to future work.
Given the function Φ b , and plotting M λ (Φ b ) and E λ (Φ b ) vs. λ, we can find the value of λ and hence the profile for a soliton solution of any desired mass or energy. Solutions of Eqs. (55)(56) can be compared with the results of numerical simulations. Ref. [29] added a distribution of mass in the form of point particles ("stars") to ULDM simulations with m = 0.8 × 10 −22 eV. The first toy model they studied included an isolated soliton (without an ULDM host halo) of mass M = 3.3×10 8 M , to which a stellar distribution was added and evolved to a virialised state. From Fig. 1 of Ref. [29] we derive the baryonic contribution to the gravitational potential, and solve for the distorted soliton at the stated ULDM mass. In the top panel of Fig. 13 we show M λ vs. λ, derived from Eq. (57). Blue solid (green dashed) lines show M λ with (without) the stellar potential. Black horizontal line denotes the value of M that was fixed in the simulation: the intersection of the black and the blue lines gives the parameter λ describing the distorted soliton. In the bottom panel we show the density profiles. Red lines are taken from [29], while green lines show the analytic soliton solutions. Clearly, the numerical profile from [29] closely matches the analytically-derived distorted soliton of Eq. (55).
Ref. [29] also simulated the soliton in an ULDM halo, the output of a DM-only cosmological simulation. Then, a stellar mass distribution was introduced as before and the system allowed to evolve.
When embedded in a halo, the mass of the soliton is not constant anymore but can grow by absorbing mass from the halo. However, regardless of the large-scale halo, in the core region the perturbed soliton profile is still fixed by the SP equations up to the ambiguity of λ. We checked that in all of the virialised soliton+halo+stellar mass simulations, presented in Ref. [29], the soliton profile matches that of the analytic distorted soliton, once accounting for the stellar potential. In the top panel of Fig. 14 we illustrate this point for the t = t H snapshot of Case C, described in Fig. 5 of Ref. [29] With stars included, the total t = t H mass distribution, given by summing the red (DM) and green (stars) solid lines in the upper panel of Fig. 14, leads to the rotation curve shown by the red solid line in the lower panel of the same figure. This rotation curve is peaked at small radius due to the distorted soliton, that is more massive than the soliton-halo prediction of the DM-only simulations. The inner velocity exceeds the prediction of Eq. (49) by a factor of two. We also show, by green dashed line, the rotation velocity at t = t H due to ULDM alone. To compare, the blue dotted curve shows the rotation curve of the initial ULDM system, extracted from the cosmological simulations. This is the same curve we showed in Fig. 5; as we discussed, it satisfies Eq. (49) to 20%.
These numerical results suggest that the presence of baryonic matter tends to make the soliton somewhat more compact and massive than predicted by the pure ULDM soliton-host halo relation (31). This, in turn, in-creases the peak in the rotation curve due to the soliton, suggesting that our constraints on ULDM obtained in Sec. IV B are robust with respect to inclusion of baryons. We now consider an example, where baryonic effects are expected to be important: the Milky Way.

B. The Milky Way: Nuclear Bulge vs. Soliton
Refs. [6,7] commented that for a MW-like galaxy, Eq. (31)  as to what extent this effect can be important. Our goal in this section is to provide a preliminary study along these lines, using photometric baryonic mass estimates.
As an interesting outcome, we find that precision kinematical studies of the MW inner bulge could in principle be sensitive to ULDM with 10 −21 eV m 10 −19 eV, where the analysis along the lines of Sec. IV B may become challenging. Fig. 15 shows a spherically-averaged enclosed mass profile, derived for the MW via various dynamical tracers [38][39][40][41][42][43][44][45][46][47][48][49][50][51]. Flattening of the enclosed mass at small radii reflects the contribution of the SMBH. Soliton profiles as predicted by Eq.  Fig. 5 of Ref. [29]). Green and red solid lines show the stellar and ULDM densities, respectively, taken from the simulation. Blue dashed show a distorted soliton solution, computed including the stellar potential. Dotted black shows an un-distorted soliton with the same mass. Bottom: Rotation curve for the total ULDM+stellar system (red solid) and including only the ULDM contribution (green dashed), at t = tH . We also show the initial ULDM only distribution at t = 0 (blue dotted).
computed holding their mass fixed by Eq. (31), but including the effects of baryons as we discuss below. An NFW profile, fitted in Ref. [36] to r 10 kpc SDSS data, is shown in dashed black.  We stress that the purpose of Fig. 15 is to illustrate the possible signature of ULDM in the inner MW, and not for statistical analyses of the MW mass distribution. Modeling the inner MW is a complicated task. The measurement of inner kinematics of the galaxy, below a few kpc, is subject to large systematic uncertainties due, among other issues, to the effects of the Galactic bar and spiral arm structures [51], which impact tangent-point velocity measurements like those utilised in [52,53]. Our simplified derivation of the spherically-averaged mass profile in Fig. 15 combines many tracers with different systematics, and accounts for none of these subtleties.
Ref. [54] analysed the MW central gravitational potential using a large set of observational constraints. In addition to the classical bulge and disc, Ref. [54] found dynamical evidence for the presence of a mass component of ∼ 2 × 10 9 M extending to ∼ 250 pc. This mass component is visible as a mass bump in Fig. 15 (see, e.g. green data points extracted from [52]). Consistent with comments in [6,7], the bump is in tantalising agreement with the soliton prediction of Eq. (31) for m = 10 −22 eV (blue shaded band).
Unfortunately, there are about a billion stars in there, too: the bump in the mass profile at r ∼ 200 pc has been associated in the literature with the nuclear bulge (NB). Ref. [55] fitted the NB mass and light by a dense disc of stars, with mass density ρ * ∼ 200 M /pc 3 , scale hight 45 pc and scale radius ∼ 230 pc. In all, the NB is thought to contain (1.4 ± 0.6) × 10 9 M in stars, roughly enough to match the dynamically inferred mass. Subsequent kinematic detection supporting the stellar mass and disc-like morphology of this component was given in [56]. Microlensing analyses [57] lend further support to the results of [54][55][56] down to r 220 pc.
The photometrically-derived NB mass model of [55] is superimposed as purple line in Fig. 15. We stress that the photometric derivation is subject to large uncertainties due to the need to correct for very strong extinction and due to unknown stellar mass-to-light ratios. What we learn from this photometric mass model, therefore, is that stars could plausibly account for all of the kinematically inferred mass in this region.
Assuming that the NB is due to stars, we now use a toy model of this mass distribution to see its effect on an ULDM soliton. We replace the disc-like morphology of the NB in [55] by a spherical model with the same radially averaged mass. The nominal model, containing the NB and additional subleading components described in [55], contains ∼ 1.7 × 10 9 M in stars inside of r = 300 pc. Adding a SMBH of M BH = 4.3 × 10 6 M [38], we calculate soliton solutions in this baryonic potential. For M λ 3 M N B ∼ 5 × 10 9 M , the NB makes a negligible impact on the soliton. For larger ratio of the stellar to ULDM mass, M λ M N B , the NB becomes important, contracting the soliton profile. For the MW, this is the parametric region predicted by Eq. (31), implying that the solitons would receive significant distortion. In Fig. 15 we illustrated this effect by presenting, in shaded bands, the soliton mass profiles computed accounting for the nominal NB model. We observe that the solitons for m in the range (10 −22 ÷ 10 −19 ) eV are expected to affect the potential at an order unity level.
Thus, a dedicated analysis using our formalism to calculate steady state soliton profiles consistent with a given baryonic mass model could be sensitive to ULDM solitons all the way up to m ∼ 10 −19 eV. For larger m, the expected soliton becomes subdominant with respect to the SMBH.
Another limitation comes from absorption of the ULDM from the soliton by SMBH. As discusses in App. B 2, the accretion rate is negligible as long as m 5×10 −20 eV. However, for higher ULDM masses the time scale of accretion can become shorter than the age of the universe and our steady-state Newtonian analysis in this section may not apply.

A. Comparison to previous work
An earlier analysis of ULDM halos including the effect of baryons was done in [2] (see also [3]). This work attempted to fit the rotation curves of spiral galaxies from [58], assuming that the entire ULDM halo is contained in a giant soliton (see also [59] for a similar approach in modeling low surface brightness galaxies). This exercise led [2] to report m < 10 −23 eV, in tension with Ly-α data. This exercise is different from the current work. As suggested by numerical simulations, we expect heavier ULDM particles with m 10 −22 eV to produce galaxies with a large scale host halo following roughly the usual NFW profile, in which the soliton affects the rotation curve only in the inner part of the halo.
Refs. [6,8] performed a Jeans analysis, fitting a soliton+halo configuration to a small sample of dispersion-dominated MW-satellite dwarf spheroidal galaxies (dSph). They found that the modelling of stellar kinematics in the dSph is consistent with a cored inner region, as expected in ULDM, with good fits to the data for m 10 −22 eV. The fitted core found in [6] was consistent with expectations from the soliton-host halo relation. In [8], the transition between an assumed host halo NFW profile and the soliton was modelled using corresponding free parameters. To conclude, these analyses considered a small sample of dispersion-dominated galaxies, as opposed to our larger sample of rotation-dominated galaxies. Their best-fit ULDM mass m 10 −22 eV is disfavoured by our findings.
Ref. [26] fitted ULDM soliton+halo profiles to SPARC [28] galaxies. Without a soliton-host halo relation, Ref. [26] assigned separate free parameters to the NFW halo and to the soliton, for each galaxy. As a result of this freedom, it is difficult to understand from that analysis if ULDM is ruled out or not in the mass range m = (10 −22 ÷ 10 −21 ) eV. Nevertheless, Ref. [26] did state a best fit of m ≈ 0.5 × 10 −23 eV, and noted the tension with Ly-α data, in agreement with our findings. Our work here shows that the soliton-host halo relation reduces the modelling freedom, and leads to disagreement with rotation curve data.
Recently, Ref. [60] addressed the problem from a different angle. It builds on the results of Ref. [61] that fitted rotations curves of a sample of galaxies assuming that their halos are described by the Burkert profile. Ref. [60] points out a correlation between the parameters of the fitted profiles -the central dark matter density ρ c and the size of the core R c , with β ∼ 1. This is different from β = 4 predicted by ULDM if the cores are identified with the solitons, see Eq. (30). Thus, Ref. [60] concludes that ULDM cannot explain the origin of the halo cores, and cannot solve the core-cusp rpoblem of ΛCDM. Our results are consistent with these findings and strengthen them by showing that the soliton-host halo relation implies a tension between predictions of ULDM and the data for m ∼ (10 −22 ÷ 10 −21 ) eV, disfavoring this mass range altogether.
B. Caveats

Soliton formation time
Our analysis has been based on the assumption that ULDM solitons satisfying the soliton-host halo relation exist in the centres of all (or most) galactic halos. Formation of such solitons within the lifetime of the universe has been demonstrated in numerical simulations [6,62] for ULDM masses m ∼ 10 −22 eV and halos with mass (10 9 ÷ 10 11 ) M . Here we discuss whether such solitons have enough time to form for higher values of m and larger halos.
Recently, Ref. [63] considered the scenario, where solitons form as a result of gravitational relaxation of the inner part of the halo. In this case the relaxation time can be estimated as, where v and ρ are the mean velocity and density in the inner part of the halo, Λ = log(mvR) is the Coulomb logarithm depending on the size R of the relaxation region, and b is a numerical coefficient close to one. The formula has been confirmed in [63] by numerical simulations that found b to range between 0.7 ÷ 0.9 depending on the precise form of the initial conditions. Substituting into Eq.
We observe that the relaxation time has strong dependence on the velocity and the density in the central region of the halo, that have large uncertainties. Fixing the fiducial values of v = 100 km/s and ρ = 0.1M /pc 3 one would conclude that the relaxation takes longer than the age of the universe if m 3 × 10 −22 eV.
However, it can be incorrect to use Eq. (61) for the estimate of the soliton formation time inside ULDM halos in the cosmological context. Indeed, this equation has been derived and tested in the kinetic regime corresponding to large velocity dispersion of ULDM particles, such that their de Broglie wavelength is much shorter than the size of the system. Starting from initial conditions of a gas of wavepackets with large velocity dispersion in a box, Ref. [63] observed a fast formation of a virialised halo, followed by a long period of kinetic evolution populating the low-lying energy levels, until the system was able to form the soliton through Bose-Einstein condensation. The time scale (61) refers to the long second stage of the process. On the other hand, the cosmological initial conditions prior to virialisation are characterised by very low velocities, where the kinetic description does not apply. In this situation a long relaxation seems unnecessary and the formation of the soliton can proceed much faster, on the scale of the halo free-fall time, see discussion in [63]. This is supported by the simulations with cosmological initial conditions which appear to imply formation of the soliton in every halo at the moment of virialisation [6,62]. We believe that the issue of the soliton formation time requires further investigation, including possible effects of baryonic matter.

Non-gravitational interactions
In our analysis we have neglected any non-gravitational interactions of the ULDM particles. Let us discuss under which conditions this approximation is justified. As an example we consider quartic self-interaction of the scalar field φ with coupling constant κ, The self-interaction is negligible as long as the corresponding potential energy is much smaller than the gradient energy of the field, For a soliton with the scale parameter λ, this translates into a bound, Using the soliton-host halo relation in the form of Eqs. (29), (30) we obtain the condition, If this condition is satisfied, our analysis is applicable. Otherwise, the effects of the self-interaction on the soliton properties must be taken into account.
While the condition (66) appears stringent, it is naturally fulfilled in a broad class of theories containing axionlike particles. These typically have a periodic potential of the form, For the masses of interest m ∼ (10 −22 ÷ 10 −18 ) eV, the axion periodicity should be in the range f ∼ (10 16 ÷ 10 17 ) GeV to produce the right dark matter abundance via the misalignment mechanism (see e.g. [12]). Expanding the potential (67) in powers of φ we obtain the value of the coupling, We see that (66) is indeed satisfied for all galactic halos considered in this paper, independently of the ULDM mass.
Similar analysis can be performed in the case of nongravitational interactions between ULDM and baryonic matter, though in this case it will be strongly modeldependent.

VII. SUMMARY
Bosonic ultra-light dark matter (ULDM) is an interesting paradigm for dark matter. For particle mass m ∼ 10 −22 eV, ULDM would form solitonic cores in the center of galaxies, and it has been suggested that this could solve several puzzles of ΛCDM on small scales.
We analysed the results of numerical simulations of ULDM, which have found scaling relations between the mass of the central soliton and the mass or energy of the host halo. Simulations by different groups converge on a soliton profile in good agreement with the analytical solution of the Schroedinger-Poisson (SP) equation, admitting important analytic insight into the numerical results.
We showed that the simulations of Ref. [13] contain a central soliton that dominates the total energy of the entire soliton+halo system. This situation is unlikely to describe realistic galactic halos 9 with mass above ∼ 5 × 10 7 m/10 −22 eV −3/2 M .
We have demonstarted that the soliton-host halo relation found in the simulations of Refs. [6,7] can be summarised by the statement, that (E/M )| soliton = (E/M )| halo . The simulations of [7] show a small spread, less than a factor of two, around this relation. This E/M relation could apply to real galaxies.
The E/M soliton-host halo relation implies that the peak circular velocity of the large-scale halo should approximately reproduce itself in the inner soliton region. This can be tested against observations without free parameters. Contrasting this prediction with high resolution disc galaxy data, we showed that ULDM in the mass range m ∼ (10 −22 ÷ 10 −21 ) eV is disfavored by observations. Given that smaller m is in tension with cosmological measurements, this disfavores ULDM as a solution to the small scale puzzles of ΛCDM.
We also analysed the effect of a fixed background distribution of baryons and of a super-massive black hole (SMBH) on the soliton solution. We showed that this analysis would be important towards using highresolution kinematical data in the Milky Way or Milky Way-like galaxies. Mapping the Milky Way gravitational potential down to r ∼ 10 −2 pc may allow to probe ULDM with mass up to m ∼ 10 −19 eV.
It should be stressed that the E/M soliton-host halo relation (35) is an empirical result, deduced from numerical simulations without baryons and tested in a limited range of ULDM and halo masses. The importance of its phenomenological implications motivates further investigation to better understand the physics underlying this relation and map out its domain of validity.
Throughout the paper we have assumed that ULDM comprises the total amount of dark matter. We leave for future the study of scenarios where only a fraction of dark matter is in the form of ULDM.
Appendix B: Adding a super-massive black hole Most galaxies, if not all, host a SMBH, and we should check how it affects the ULDM soliton. This is an important point to verify and we do this analysis here in some detail. The upshot of our discussion is that a SMBH would not affect our results for small disc galaxies, whereas in the MW it becomes relevant at the upper end of considered ULDM masses, m 5 × 10 −20 eV.
We start with a Newtonian analysis, dealing with the ULDM configuration far from the BH Schwarzschild radius. Then we will consider the limitations of this analysis imposed by absorption of ULDM on SMBH that leads to the decay of the soliton+SMBH configuration.

Soliton shape in the presence of SMBH
Adding a SMBH, coincident with the soliton center of mass, changes Eqs. (7)(8) into where We chose the reference value for M BH to represent the case of the MW. It remains convenient to solve the problem using boundary conditions with χ(0) = 1. Let us denote this solution (satisfying χ(0) = 1) by χ 1 (r; A), accompanied by the potential Φ 1 (r; A). Having found χ 1 (r; A) and Φ 1 (r; A) for any value of A, the physical solution χ λ (r; A) satisfying Eqs. (B1-B2) with boundary condition χ λ (0; A) = λ 2 is given by Defining M λ (A), with obvious notation, the soliton mass is Fig. 17 shows the density profile χ 2 1 (r; A) for different values of A. For A 10 −2 the solution converges to the unperturbed χ 1 . For A 1 the solution approaches χ 1 (r; A → ∞) → e −Ar which is nothing, but the wavefunction of the ground state in the Coulomb potential. This means that for A/λ 1, solitons with different values of λ have the same profile up to over-all normalization: χ λ (r; A λ) → λ 2 e −Ar . In this limit the soliton core radius is r cλ (A) → ln(2)/(2A), and the mass is m . More generally, for any A the massradius relation is The mass-radius product is shown by the solid blue curve in Fig. 18 (top panel). The large-A analytic so- As discussed in Sec. III A, Eqs. (31) and (35) are equivalent when the ULDM soliton is self-gravitating, but this equivalence is broken by the SMBH gravitational potential. It is unlikely that Eq. (31) remains realistic when the SMBH becomes important. Instead, Eq. (35) may still be valid. In computing the ULDM energy we need to account for the SMBH contribution, which means that The large A limit (attained for a χ λ (r; 8A M 2 pl m . In this limit E λ /M λ tends to a λ-independent constant, E λ /M λ (A → ∞) → −A 2 /2, because the soliton is not held together by self gravity, but by the SMBH potential. The E/M relation for the χ 1 (A) soliton is shown in the bottom panel of Fig. 18.
In Fig. 19, solid lines show E/M for the soliton with different values of particle mass m, in the presence of a SMBH with M BH = 4 × 10 6 M . The dashed lines show the undistorted result, obtained from Eq. (25). The large A limit with constant E λ /M λ is attained when the soliton mass becomes smaller than the SMBH mass.
To further clarify this point we relate A/λ to the un- Next, we see from Eqs. (B4), (B5) that A/λ gives an effective A-parameter that must be substituted in the χ 1 -soliton to obtain the actual solution. Thus, whether the effects of SMBH are important or not can be read from the ratio of the SMBH mass to the mass of the undistorted soliton.
To analyse the rotation curve for the soliton+SMBH system, it is useful to view V circ for the χ 1 (A) soliton, obtained for different values of A. This is shown in For M λ 10 M BH , the soliton's circular velocity peak and the peak location are affected very little by the SMBH. However, at small x approaching the SMBH, the velocity curve turns back up reflecting the gravitational potential of the SMBH itself. This implies that the rotation curve of a galaxy hosting a SMBH may not decrease appreciably below x peak,λ , even when M BH is much smaller than M λ . For example, even for M BH = 0.05 M λ (A = 0.1 in Fig. 20), the circular velocity decreases by only ∼ 25% at x < x peak,λ before it goes back up due to the BH. For the same parameters, Fig. 18 shows that E λ /M λ is essentially unperturbed, meaning that the soliton-host halo relation of Eq. (35) and the rotation curve properties discussed in Sec. IV remain unaffected. Such a galaxy would have an approximately  flat rotation curve all the way in to the region of SMBH dominance.
The impact of a SMBH on the analysis of Sec. IV can be concluded as follows. For M λ 4M BH (corresponding, by using Eq. (B10), to effective A/λ ≈ 0.5), Fig. 18 shows that E λ /M λ is corrected by 25% compared to its unperturbed value. Eqs. (28) and (44) then imply that λ and maxV circ,λ are corrected by an insignificant ∼ 12%. For a large SMBH with M BH M λ , on the other hand, the SMBH gravitational potential itself must produce high rotation velocity at x peak,λ , comparable or higher than what the soliton itself would provide. As a result, besides from the fact that the rotation curve may not decrease towards lower x below the soliton peak, our analysis in and observational limits from Sec. IV should be robust against the effect of a SMBH 10 .

Absorption of soliton by SMBH
We now discuss accretion of ULDM from the soliton onto SMBH. To estimate the lifetime of the soli-ton+SMBH configuration we follow the approach of [12] (see also [65,66] for the study of scalar bound states in external Schwarzschild metric and [67] for a numerical analysis of the lifetime of soliton+BH system). Unruh [68] has derived the cross section for absorption of a scalar particle with mass m and momentum k by a Schwarzschild BH, whose size is much smaller than the Compton wavelength of the particle. In the nonrelativistic limit, k m, it has the form, where ζ = 2πGM BH m 2 /k .
As the cross section is s-wave dominated, it is appropriate for calculation of accretion from a spherically symmetric scalar field configuration. The growth of BH mass due to an infalling flux of particles with density ρ is, In the case of the soliton we can estimate the absorption rate by substituting into this expression the central density of the soliton and the characteristic momentum of particles in the soliton that can be read from the soliton profile (30),