Ultra-light dark matter in disk galaxies

Analytic arguments and numerical simulations show that bosonic ultra-light dark matter (ULDM) would form cored density distributions (`solitons') at the center of galaxies. ULDM solitons offer a promising way to exclude or detect ULDM by looking for a distinctive feature in the central region of galactic rotation curves. Baryonic contributions to the gravitational potential pose an obstacle to such analyses, being (i) dynamically important in the inner galaxy and (ii) highly non-spherical in rotation-supported galaxies, resulting in non-spherical solitons. We present an algorithm for finding the ground state soliton solution in the presence of stationary non-spherical background baryonic mass distribution. We quantify the impact of baryons on the predicted ULDM soliton in the Milky Way and in low surface-brightness galaxies from the SPARC database.


I. INTRODUCTION
An ultra-light bosonic field oscillating around a minimum of its potential [1][2][3][4] can play the role of dark matter (DM). On cosmologically large scales ultra-light dark matter (ULDM) behaves similarly to cold weakly-interacting massive particle (WIMP) dark matter, reproducing its success with respect to the cosmic microwave background and largescale structure. On smaller scales comparable to the de Broglie wavelength, ULDM behaves differently to WIMPs. In particular, at the centre of galactic halos ULDM develops cored density profiles that lead to markedly different predictions than those found for ordinary WIMPs [1,. The cored ULDM distributions correspond to quasi-stationary minimum energy solutions of the equations of motion. We will follow common convention and refer to these solutions as "solitons".
Ref. [23] analysed the rotation curves of well-resolved low surface-brightness (LSB) disk galaxies from the SPARC database [31] and pointed out that these galaxies fail to show the soliton feature predicted by numerical simulations [9,10,32] 1 . This led to the bound m 10 −21 eV. A similar constraint 2 was found in [26] considering the dwarf spheroidal galaxy Eridanus-II. The matter power spectrum inferred from Ly-α forest analyses yields a comparable bound [34][35][36][37][38] 3 . These lower bounds on m are interesting because they probe DM using gravity alone, without requiring any direct interactions with SM fields; because they define how light DM could possibly be; and also because ULDM with m ∼ (10 −22 − 10 −21 ) eV was suggested as an explanation for puzzles facing the WIMP paradigm on small scales [16,47].
In the attempt to constrain (or detect) ULDM with galactic kinematics, an important issue is the modelling of the baryonic contribution to the gravitational potential which can distort the soliton 4 . Ref. [23] analysed the solution in the presence of a spherically-symmetric background potential, in order to estimate the size of the effect. That was found to be significant for the Milky Way (MW), but not significant for the relevant SPARC LSB galaxies. However, both in the MW and in rotation-supported SPARC galaxies, the baryonic mass distribution is non-spherical, following disk-like morphology. In a non-spherical system dynamics in the central region of the galaxy can be affected by the mass distribution at larger radii. It is therefore important to extend the soliton+baryon analysis to non-spherical configurations.
In this paper we present an algorithm to calculate the soliton solution in the presence of a non-spherical background gravitational potential. The algorithm is simple, fast and accurate and can replace the standard one-dimensional shooting method used for solving the spherically symmetric soliton.
Our goals in presenting this tool are twofold. First, it allows to perform a self-consistent analysis of the velocity profile in disk galaxies. Once the baryonic mass distribution is specified (presumably with input from photometry), the soliton contribution to the gravitational potential requires a single free parameter in the fit. This parameter can be chosen to be, e.g., the soliton mass. For example, stellar kinematics in the MW could provide a testing ground for ULDM up to m ∼ 10 −19 eV [23]. To this end, implementing the soliton in a self-consistent manner would be crucial and we expect that our tool would be useful.
Second, we revisit the analysis of Ref. [23] of baryonic effects in SPARC galaxies. As noted in [23], the solitonhalo relation predicted by DM-only numerical simulations strongly over-predicts the circular velocity in the centres of dozens of galaxies if m < 10 −21 eV. In many cases, the predicted soliton mass in the central ∼ 100 pc of the galaxy exceeds the observationally allowed total mass (bary-onic+DM) in that region by factors of order 10. This large mass mismatch led [23] to expect that baryonic effects are unlikely to change the constraints. Here, focusing on two sample galaxies, we determine the soliton solution while accounting for the non-spherical baryonic mass distribution. When noting that the gas and stellar distributions are highly non-spherical, and when naively counting the mass outside of the soliton region, the total baryonic mass in both galaxies is comparable to or larger than the soliton mass. Nevertheless, in both cases our analysis largely confirms the expectations of [23], showing that the baryonic mass external to the soliton region does not significantly affect the solution.
Recently Ref. [49] considered a non-spherical parametrisation of ULDM halos and solitons and applied it to a Jeans analysis of dwarf spheroidal galaxies. Differently to our work here, Ref. [49] did not base their parametrisation of the ULDM core on a solution of the equations of motion (EOM). To our view, one of the main points of beauty in the discussion of ULDM in galaxies is that numerical simulations -with and without stars -actually do consistently show soliton solutions that satisfy the EOM 5 . Using the tools we provide here, refining the analysis of [49] by solving for the self-consistent soliton would be straightforward.
The outline of this paper is as follows. In Sec. II we recall the ULDM equations of motion that define the soliton solution in the presence of an external (non-dynamical) background gravitational potential. The standard onedimensional shooting method, that can be used to solve the spherically-symmetric problem, becomes impractical (in general) once the background potential is not sphericallysymmetric because it requires a discrete infinity of shooting variables 6 . A simple numerical recipe to solve this problem, assuming an axisymmetric background potential, is detailed in App. A. While we do not pursue this here, extending the algorithm to full 3D is straightforward.
In Sec. II A we discuss the soliton-host halo mass relation found in DM-only numerical simulations [9,10,32]. Ref. [23] showed that this relation is equivalent to the statement, that the specific energy (total energy per unit mass) of the host halo is equal to the specific energy of a self-gravitating soliton. Here we point out that a more physical representation of the soliton-halo relation is obtained by equating the kinetic -rather than total -energy per unit mass of the soliton and the halo. This distinction is unimportant for DM-only simulations of massive halos, but becomes relevant once a background potential is introduced.
In Sec. III we take the MW as an illustrative example of a system where the potential due to baryons (mostly stars in this case) cannot be neglected in assessing the soliton properties. The analysis demonstrates the use of our numerical tool, but is not intended to provide constraints on ULDM: that would require a more comprehensive treatment that we postpone to future work.
In Sec. IV we consider two sample LSB galaxies from the SPARC database. For these galaxies, we reconstruct the circular velocity decomposition presented in the SPARC database using photometric data, reproducing the SPARC analysis. The baryonic mass models (stars+gas) derived in this way are used as input for the numerical non-spherical soliton solution, allowing us to revisit in detail the earlier rough estimates of Ref. [23]. We show that the total energy per unit mass, E/M , of the soliton is modified by the baryonic potential of these galaxies. However, the bulk of this effect is unphysical: it comes from a non-dynamical shift of the energy due to an external gravitational potential that is mostly constant throughout the relevant region of the galaxy. This is supported by the fact that -as we show -the specific kinetic energy, K/M , is essentially unaffected both for the soliton and the halo.
In Sec. V we summarise our results. A number of technical details are postponed to appendices. As mentioned above, App. A describes the nonspherical soliton-finding algorithm. App. A 1 specifies the steps required to implement a black hole in the code. In App. B we recall a convenient formula converting an axisymmetric mass distribution into the gravitational potential induced by it. In Sec. C we collect a useful auxiliary parametrisation for galactic discs, that we have found useful in modelling SPARC galaxies. App. D explains our reconstruction of the neutral gas distribution in UGC01281. Finally, our results in the main text are presented -for concretenessassuming ULDM particle mass of m = 10 −22 eV; in App. E we show relevant results for m = 10 −21 eV.

II. SOLITONS IN A NON-SPHERICAL BACKGROUND
We consider a real, massive, free 7 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 and satisfies the Schrödinger-Poisson equations (SPE) [53] i∂ t ψ = − 1 2m In Eq. (2) we include an external contribution to the gravitational potential, given by Φ b . We consider Φ b as the effect of a distribution of baryonic mass. Our working assumption is that Φ b should be constrained by external information such as photometry and microlensing measurements. We look for a quasi-stationary phase-coherent solution for the ULDM, described by the ansatz where M pl = 1/ √ G. The parameter γ is an eigenvalue of the SPE subject to the bound-state boundary conditions that we describe below.
We rescale the spatial coordinate, keeping this convention throughout the rest of the paper. Then in terms of the dimensionless χ and x the SPE are We assume cylindrical symmetry and parity symmetry (x 3 = z → −z), and define the radial coordinate in the plane R = x 2 1 + x 2 2 . At √ R 2 + z 2 → ∞ the potentials Φ and Φ b are assumed to decay ∝ 1/ √ R 2 + z 2 , implying that χ decays approximately exponentially ∝ e − √ 2|γ|(R 2 +z 2 ) . A given value of χ at the origin, specified by with λ a real positive number, fixes the minimal energy solution of Eqs. (6-7) consistent with the boundary conditions. In the case of vanishing Φ b , solutions of Eqs. (6-7) admit a scaling symmetry, the orbit of which can be parametrised by λ. This scaling symmetry is, in general, broken by Φ b = 0. It remains true, however, that varying the value of λ in Eq. (8) generates the continuous family of solutions of Eqs. (6-7). Thus, λ remains a useful tool to parameterise the mass, energy and any other property of the solution. For reference, the self-gravitating soliton (found for When baryons induce Φ b = 0 these relations are modified in a way that we will discuss below. We have developed a numerical relaxation method to find the ground state soliton solution for any axisymmetric background potential satisfying the boundary conditions described below Eq. (7). The algorithm is described in App. A, and is one of the main results of this paper.
In the next subsection we clarify some issues related to the soliton-host halo relation found in DM-only numerical simulations. Then, in the following sections we illustrate the use of the numerical tool of App. A by analysing the baryonic effects on the predicted ULDM soliton in the Milky Way and in two disk galaxies from the SPARC database.
Because Φ h is approximately constant over the region where the soliton is supported, the soliton shape is not distorted and its kinetic energy is not modified from its value for the self-gravitating solution. This means that for massive halos in DM-only simulations, Eq. (13) and Eq. (12) are indistinguishable.
Eq. (13) and Eq. (12) become distinguishable when we turn on Φ b = 0, with a nontrivial spatial profile such that Φ b is not constant throughout the large-scale halo.

III. APPLICATION: THE MILKY WAY
We now consider soliton solutions in the background of a gravitational potential Φ b , chosen to roughly mimic the inner region of the MW. Our goal is to illustrate the approximate size of the baryonic effects on the soliton, and not to characterise these effects in full; a detailed, accurate and precise modelling of the inner MW stellar and gas mass distributions is challenging and is postponed to future work. For concreteness, throughout this section we set m = 10 −22 eV.
The dominant contributions to the stellar mass profile of the MW inner few hundred pc were described in the photometric analysis of Launhardt et al [54] as a spherical nuclear stellar cluster (NSC) and a nuclear stellar disk (NSD), composing together the nuclear bulge (NB).
In addition to the stellar components, dynamics in the central ∼ 1 pc is dominated by a super-massive black hole (SMBH) with mass M BH ≈ 4 × 10 6 M ⊙ . Here we omit the SMBH contribution, which was studied in [23] and shown to have negligible impact on the soliton for m 10 −20 eV. We note that the numerical code in App. A is capable of handling the SMBH contribution via the procedure described in App. A 1. A gas torus at scale radius of ∼ 100 pc contributes ∼ 2 × 10 7 M ⊙ . For simplicity, the gas is also neglected here in comparison to the stellar components.
The NSC density profile was modelled as where r = √ R 2 + z 2 is stated in pc.ρ N SC = 3.3 × 10 6 M ⊙ /pc 3 for r < r 0 andρ N SC = 9.0 × 10 7 M ⊙ /pc 3 for r ≥ r 0 , with r 0 = 6 pc. The index n N SC = 2 for r < r 0 and n N SC = 3 for r ≥ r 0 (keeping the profile continuous at r 0 ). With these parameters we have 8 We parametrise the NSD stellar mass density as follows, 8 This NSC mass is larger than that quoted in [54] by a factor of ∼ 1.8. We are not sure of the reason for this mismatch, but it does not have an important effect on our results.
whereρ N SD = 330 M ⊙ /pc 3 and where z and R are stated in pc. This parametrisation approximately reproduces the NIR stellar volume emissivity model derived in [54] and yields an NSD mass M N SD ≃ 10 9 M ⊙ , consistent within the uncertainty with the value of (1.4 ± 0.6) × 10 9 M ⊙ quoted by [54]. A kinematic detection supporting the disk-like morphology of the NSD was given in [55], and the mass and approximate scale estimates are consistent with the dynamical modelling of [56] and with microlensing analyses [57] that probe the outer boundary of the NSD region.
In what follows we define Υ L ≡ Υ/Υ Launhardt as the mass-to-light ratio of the stellar distribution compared to the one used in Ref. [54]. We vary Υ L to explore the consequences of different total stellar mass in the NB region.
In Fig. 1 we plot the soliton mass vs. λ, which allows us to access different solutions. For λ 10 −3 we retrieve the self-gravitating soliton result, shown by the dashed line. For smaller λ we find M ∝ λ 4 [23] 9 . Fig. 1 can be compared to Fig. 16 in Ref. [23] which considered a spherically-averaged approximation to the same stellar mass model. It shows an O(1) difference in the M vs. λ relation in the phenomenologically interesting range λ ∼ 10 −4 − 10 −3 . In Fig. 2 we study the deformation in the soliton shape caused by the stellar mass distribution, at fixed soliton mass M ≈ 1.  It is instructive to consider the observable (in principle) soliton-induced effective circular velocity, In Fig. 3 we plot v eff on the plane of the disk (z = 0, dashed red) and transverse to the disk on the z-axis (R = 0, dot-dashed yellow). For comparison, we also plot v eff computed for a self-gravitating soliton with the same mass (solid blue). The main effect of the background stellar potential is to contract the soliton-induced peak velocity deeper into the inner halo, enhancing the peak velocity; this is an O(1) effect that cannot be ignored in realistic modelling of kinematic data. The deviation from radial symmetry is, however, small: a simplified treatment taking as input a radially-averaged baryonic mass distribution could suffice for practical purposes. For comparison, the result of such a procedure is plotted in the dotted line in Fig. 3.
In the top (bottom) panel of Fig. 4 we plot the total energy (kinetic energy) per unit mass as a function of soliton mass M . For (summarised by Eq. (13)) predicts K/M ≈ 5.5 × 10 −8 [23], shown by the black dot-dashed line. The shaded band denotes a factor of two spread around this prediction, motivated by the halo-to-halo spread seen in the simulations. Fig. 4 shows that because of the stellar-induced background potential, K/M for an actual soliton solution in this background is significantly deformed. This means that baryonic effects are likely to significantly modify the soliton properties, and the soliton-halo expectation from DM-only numerical simulations cannot be taken at face value. A consistent way to constrain (or possibly detect) an ULDM soliton in the MW, would be by a combined analysis of kinematical modelling and photometry, where the stellar potential constrained by photometry is used to self-consistently calculate the soliton shape and where the soliton mass is taken as a free parameter.
Ref. [58] argued for dynamical evidence in favour of an ULDM soliton in the MW, with m ≈ 10 −22 eV and M ≈ 10 9 M ⊙ in tantalising agreement with the expectations of DM-only numerical simulations. The dynamical evidence for a dense central mass component is consistent with earlier studies [23,[54][55][56]. Unfortunately, as we reviewed here and in [23] (see Sec. V.B there), there is room for and photometric evidence of about 10 9 M ⊙ in stars within the ∼ 200 pc would-be soliton region [54]. Thus, the central mass component could well be due to ordinary baryonic matter. Other systems, such as well-resolved LSB galaxies, offer much cleaner laboratories in which to look for ULDM solitons. We turn to such systems in the next section.

IV. APPLICATION: LOW SURFACE-BRIGHTNESS SPARC GALAXIES
Our second discussion of non-spherical solitons involves two low surface-brightness (LSB) disk galaxies from the SPARC database [31]: UGC01281 and F571-8. We choose these galaxies as representative examples of a larger sample including dozens of well-resolved LSB galaxies. For concreteness, throughout this section we set m = 10 −22 eV. Results for m = 10 −21 eV are collected in App. E.
The baryonic mass contributions in SPARC galaxies is divided into a spherical bulge component and axisymmetric disk and gas components. The stellar mass distribution is calibrated to match surface brightness data from Spitzer. The computation of the gravitational potential due to the disk is detailed in App. C. We focus here on galaxies that are consistent with negligible bulge.
The gas mass distribution for UGC01281 (not relevant for F571-8) is calibrated to approximately match the HI surface brightness data reported in [59], normalising to the total gas mass reported in [60]. We provide details on the gas fitting procedure in App. D.
In our computation we fix the total gas mass to match the total mass inferred from the photometry and vary the stellar mass-to-light ratio of the disk from Υ d = 0 up to larger values that saturate the observed kinematic velocity [61].
In Fig. 5 we plot the M − λ relation for a soliton in UGC01281. In the top (bottom) panel of Fig. 6 we plot the total energy (kinetic energy) per unit mass vs. M . The dashed black line denotes the soliton-halo prediction of DM-only numerical simulations. The shaded band shows a factor of two spread around this prediction. Inspecting Fig. 6 we see that in the neighbourhood of E/M values that conform to the DM-only simulation prediction, the actual E/M for a soliton in UGC01281 is significantly shifted compared to the self-gravitating solution. However, the effect on K/M is much less pronounced: the soliton shape is essentially unaffected.
We can also estimate the baryonic effect on the dynamics of the large-scale halo. To do this, we can compare the observed kinematic velocity at large distances (r ∼ 5 kpc in this example) with the contribution to the velocity that can be attributed to the baryons. The velocity decomposition is shown in the top panel of Fig. 9 (discussed in more detail at the end of this section). We find v 2 baryons /v 2 obs ∼ 0.26 (∼ 0.39), when adopting Υ d = 1.07 (Υ d = 2.14). This means that the baryonic potential distorts the ULDM largescale halo K/M by no more than 40%.
The next galaxy we consider is F571-8. Soliton properties for this galaxy are presented in Figs. 7 and 8. Here, for simplicity, we ignore the (negligible) gas contribution in computing the soliton. Again, K/M for a soliton in F571-8 is unaffected by baryons in the parameter region expected from DM-only simulations. The case of F571-8 is even clearer than UGC01281 because the baryonic effect on the dynamics of the large-scale halo, as seen by inspecting the rotation curve decomposition (bottom panel of Fig. 9), is not larger than ∼ 5%.
In the top (bottom) panel of Fig. 9 we show the rotation curve decomposition of UGC01281 and F571-8, as found in the SPARC database. The contribution due to soliton solutions with different values of λ (indicated in the plot) are overlaid in red, blue and black. The solitons are computed assuming different values of the disk stellar mass-to-light ratio Υ d , listed in the caption. The central value of λ (in blue) corresponds to the prediction of DM-only numerical simulations. For these predicted solitons the baryonic potential makes a negligible impact on the soliton shape, regardless  We conclude that if Eq. (13) correctly captures the soliton-halo relation of the simulations, then UGC01281 and F571-8 are clean systems in which to constrain the ULDM model, in the sense that the baryonic contribution to the gravitational potential is not important both for the largescale halo and for the central soliton. These conclusions stay unchanged when we consider more massive ULDM with m = 10 −21 eV (see App. E). Dozens of other comparably clean systems exist in the SPARC database. The constraints derived in Ref. [23] should therefore apply and ULDM with m < 10 −21 eV is in tension with the data.

V. SUMMARY
An ultra-light bosonic field oscillating around a minimum of its potential can play the role of dark matter (DM). On scales of order the effective de Broglie wavelength, wave mechanics dictates the dynamics of this ultra-light dark matter (ULDM) opening potential avenues to constrain (or detect) ULDM in various astrophysical and cosmological systems.
Stellar and gas kinematics of rotation-supported low surface-brightness (LSB) galaxies were used in Ref. [23] to derive the constraint m 10 −21 eV. This constraint relies on the validity of a soliton-host halo relation, found in DM-only numerical simulations. It is important to assess to what extent baryons could affect these results. For a nonspherical baryonic distribution, a new numerical tool was required in order to calculate the properties (shape, mass, energetics) of the non-spherical soliton obtained in the presence of the baryonic-induced background gravitational potential. In this paper we provided a simple algorithm (see Sec. II and App. A) that achieves this goal.
To illustrate the potential use of the non-spherical soliton solver, we estimated the impact of a Milky Way (MW) nuclear stellar disk (NSD) on an ULDM soliton. Adopting a plausible parameterisation of the stellar distribution, motivated by photometric measurements, we find that the NSD would distort the shape and energetics of an m = 10 −22 eV ULDM soliton at the O(1) level. Thus, an attempt to constrain ULDM in the MW should self-consistently account for the gravitational effect of stars. While we did not enter such an analysis, the numerical tool we provided is an important step in this direction. Having said that, we note that while the soliton can be compressed by an internal clump of stars it is not easily deformed into non-spherical shape. In the MW example, the highly non-spherical nuclear stellar disk (NSD) leads to a soliton that is significantly contracted but remains spherical to a good approximation. As a result, a spherical rearrangement of the stellar mass distribution (namely, replacing the disk-like baryonic distribution by a radially-averaged profile) would most likely be sufficient to calculate the soliton in a kinematical analysis.
Next, we revisited the SPARC galaxy analysis of [23]. Using two LSB galaxies as a concrete example, we modelled the baryonic potential consistent with photometric data and bracketed the possible impact on the shape and energetics of the predicted soliton. Our results reinforce the conclusions of [23], implying that baryons are not expected to change the constraints derived on ULDM based on rotationallysupported LSB SPARC galaxies.
Israel Science Foundation and by grant 1507/16 from the Israel Science Foundation. The work of JE was supported by the Zuckerman STEM Leadership Program.

Appendix A: Numerical algorithm for solitons in a non-spherical background potential
In what follows we describe a numerical method to find the ground state solution of Eqs. (6-7). This tool is one of the main results of our work: it is intended to be simple, fast, and robust enough to allow it to be used in detailed analyses of galactic kinematics with ULDM, in cases -such as the Milky Way galaxy -where the baryonic contribution to the gravitational potential in the soliton region cannot be neglected.
We assume that the baryonic-induced gravitational potential is a direct input to the code. Often, an input in terms of the stellar and gas mass density could be more natural. Converting an axisymmetric mass distribution into its corresponding gravitational potential is a straightforward exercise that we recall in App. B.
We use an N × N discretised lattice with physical size L × L in the R − z plane. The lattice spacing is δ = L/(N − 1). The physical coordinate of each point (recall that distance is measured in units of 1/m) is The Laplacian in cylindrical coordinates is We discretise it: Note that we do not need to define ∇ 2 Φ at i = N and/or j = N . We start by initialising Φ as zero everywhere, assigning an initial test profile of χ that is conveniently chosen as some numerical approximation of the known self-gravitating solution; see, e.g. [11]. Throughout the calculation we enforce The discretised Eq. (7) is then solved iteratively using the successive over-relaxation (SOR) method (see, e.g., ch.19.5 in [62]). In each iteration of the program, Φ i,j (i, j = 1, N ) is improved by the SOR method as where ω Φ is an auxiliary parameter 10 that we set O(1). For i = 1 and/or j = 1, the RHS of Eq. (A6) should be modified 10 To obtain our results in this paper, we have used ω Φ = 1.6 in all according to Eqs. (A3-A4). At i = N and/or j = N , Φ i,j is fixed by the following boundary conditions, Here, r i,N = 1 + (i − 1) 2 /(N − 1) 2 L and the dimensionless 11M is calculated as consistent with Gauss' Law (see Sec. A 1 below). Next, once Φ is fixed, the ground state solution of Eq. (6), χ 0 , can be found by considering the following imaginary In the large τ limit, the asymptotic behaviour of χ is Thus, in each iteration, χ i,j (i, j = 1, N ) is improved as where ω χ is an auxiliary parameter 12 that we set O(1). For i = 1 and/or j = 1, the RHS of Eq. (A11) should be modified according to the prescription in Eqs. (A3-A4). At i = N and/or j = N , χ i,j is fixed by Eq. (A5). We repeatedly update Φ and χ, using Eqs. (A6) and (A11), until convergence is attained. The eigenvalue γ is calculated as To calculate the total soliton energy, we use Eq. (10) (averaging over adjacent grid sites can be useful in order to reduce numerical error): [e i,j + e i,j+1 + e i+1,j + e i+1,j+1 ] 4 (A14) 12 To obtain our results in this paper, we have used ωχ = 0.8 in all computations. This value was chosen somewhat arbitrarily based on tests of the rate of convergence. We note that setting ωχ < ω Φ appears to be useful (see footnote 10).
with the integrand For i = 1 and/or j = 1, the RHS of Eq. (A15) should be modified according to the prescription in Eqs. (A3-A4). We do not need to include i = N and/or j = N , because there the integrand vanishes due to the boundary condition in Eq. (A5). This concludes the description of the numerical scheme.

Adding a black hole
Here we explain how a central black hole (BH) can be added to the discretised grid calculation of App. A. To this end we derive a discretised version of Gauss's Law.
Using Eqs. (A3-A4), for n < N , we obtain where κ i and λ j are defined as κ i = π/4 (i = 1) 2(i − 1)π (i = 1) , λ j = 1 (j = 1) 2 (j = 1) . (A18) From these equations we obtain where n i , n j < N . This becomes the usual Gauss's law Consider a black hole with physical mass M BH , translated in our conventions to M BH = (M 2 pl /4πm)M BH . It gives the potential Φ BH = −M BH /(4πx). The Poisson equation is ∇ 2 Φ = ρ; thus, the discretised ρ configuration that leads to Φ BH and that is consistent with Gauss's Law, Eq. (A19), is Finally, to interface to the code described in App. A, it is convenient to utilise the gravitational potential induced by the BH which is given on the axisymmetric grid as follows, The tricky point here is Φ BH,1,1 : this is determined by solving the discretised Poisson equation at the origin.
Appendix B: Gravitational potential of an axisymmetric mass distribution The solution of the Poisson equation in axisymmetry can be found directly using the method of Fourier-Bessel transform. Following [63], the gravitational potential is given by where the kernel K is given by with Q − 1 2 the Legendre function of the second kind of order − 1 2 . See also 6.612 (3) and 8.834 (1) in Ref. [64].
Our analysis is less sophisticated than that in [59], but captures the key features of the gas profile with sufficient accuracy. We model the gas density profile as a collection of K co-planar rings, with the mass density of each ring taking to be constant on the plane (z = 0) and decaying vertically with a Gaussian profile: where θ(x) is the Heaviside function. The gravitational potential due to this mass distribution is computed by the procedure given in App. B.
The surface brightness profile from this gas distribution is easily computed. Matching the model to the vertical profile reported in [59], we find a good fit for d z = 0.65 kpc. Considering the radial profile and matching (approximately, by eye) to the average profile shown in Fig. 2 of Ref. [59] (which averages the HI column density over a slab in the vertical direction), with find that a model of K = 50 rings of equal width ∆ k = 0.2 kpc, located with inner radii starting at R 1 = 0 kpc up to R 50 = 10 kpc, reproduces the brightness profile radial shape for the density assignment ρ k =ρ (0.5 + R k ) 1.2 exp − R k 1.5 1.4 , where R k are noted in kpc andρ is an over-all normalisation factor. We setρ = 3.9 × 10 6 M ⊙ /kpc 3 , so that the total gas mass (including a factor of 1.3 to account for He) is fixed to M gas = 3.2 × 10 8 M ⊙ , inferred in Ref. [60] from the total HI luminosity. The gas-induced rotation curve we find with this procedure is shown by the line in Fig. 11, compared to the velocity contribution attributed to the gas in the SPARC database (circles). The comparison is good enough for our purpose in the current work: as we show in the body of the work, the total baryonic effect (stars and gas combined) on the predicted soliton and on the large-scale halo of UGC01281 is small. We conclude this technical discussion with an amusing comment. The toroidal gas profile of UGC01281 prompted us to look for toroidal soliton solutions, that could co-exist in the background potential of such a baryonic mass distribution. Indeed, varying the gas mass and the soliton mass, we can find toroidal solitons; we show an example in Fig. 12. The parameters chosen to achieve this toroidal solution were: m = 10 −22 eV, with λ = 10 −5 and a gas mass 50 times larger than the observed one in UGC01281. These parameters do not represent an actual galaxy from SPARC: we merely bring it as an observation about deformed solitons and as demonstration of the versatility of the numerical code. Appendix E: Results with ULDM particle mass of m = 10 −21 eV Here we present a repetition of Figs. 5-8 from Sec. IV, done for ULDM particle mass m = 10 −21 eV.