Virialized Profiles and Oscillations of Self-interacting Fuzzy Dark Matter Solitons

We investigate the effect of self-interactions on the shape and oscillations of the solitonic core profile of condensed fuzzy dark matter systems without the backdrop of a halo, revealing universal features in terms of an appropriately scaled interaction strength characterizing the crossover between the weakly- and strongly-interacting regimes. Our semi-analytical results are further confirmed by spherically symmetric simulations of the Gross-Pitaevskii-Poisson equations. Inverting our obtained relations, we highlight a degeneracy that could significantly affect constraints on the boson mass in the presence of repulsive boson self-interactions and propose the simultaneous extraction of static and dynamical solitonic features as a way to uniquely constrain both the boson mass and self-interactions.


I. INTRODUCTION
Fuzzy dark matter (FDM) has lately emerged as a popular alternative to the conventional cold dark matter (CDM) model.Following the first proposition of an ultralight bosonic dark matter constituent [1], an ultralight particle within the mass range 10 −22 − 10 −20 eV, leads to a description based not on incoherent particles, but on a single coherent wavefunction, with a de Broglie wavelength of ∼ O(1) kpc.Numerical simulations [2][3][4] have shown that on large scales, FDM replicates the cosmic web structure of cold dark matter, whilst addressing some alleged challenges that CDM faces on galactic scales.Formation and growth of structures in such a dark matter background has been thoroughly investigated recently [5][6][7][8].The favourable look upon this type of dark matter in recent years has been motivated by the inherent properties which remedy ailments from which the CDM model allegedly suffers [9][10][11].Among these, most notably, is the cusp-core problem [12] which is resolved in FDM since the balance between the quantum pressure and gravitation in FDM naturally forms a core, embedded in a NFW-like halo [13].Furthermore, some tentative, more direct favourable evidence for FDM may have begun to emerge, see e.g.[14].Scrutiny of this model has become an active area of research -see e.g.[15][16][17] for recent reviews and references to the literature.
The cores in FDM halos, often termed FDM solitons, exhibit interesting dynamical behaviours.Given nonsymmetric initial conditions, the cores formed from gravitational collapse exhibit a random walk within the base of the gravitational potential [18,19].Furthermore, scalar field oscillations are observed in numerical simulations [2,3,[18][19][20][21][22][23], and significant advances have been made in understanding the parameter dependence of their frequency [24,25].These oscillations manifest in the form of homogeneous radial expansion and contraction of the soliton core, and have been found to send out density waves into the surrounding halo [20].Higher excited asymmetric states are also possible but they will not be discussed in this paper.In such non-symmetrically initialised simulations, the constructive interference of density waves and further stochastic density fluctuations [26,27] also generate granules with a typical scale of the order of the central soliton core [2,3,28].These granules have been stipulated to be the driving force behind the dynamical heating of stellar streams in the Milky Way [29], whilst density waves have been proposed as potential sources of disk heating [26].Although the relation and interactions of the solitons to their host halos are very interesting and still to be fully explored, in this paper we will examine the FDM soliton in isolation, without the backdrop of a complete halo.
Most FDM investigations have considered noninteracting bosons, for which the physical state of the system is described by a Schrödinger equation, coupled to the Poisson equation, in what is known as the Schrödinger-Poisson system of coupled equations (SPE).However, lacking any evidence to the contrary, the addition of particle-particle interactions renders the situation more interesting and dynamically complex -such a set-up, based on contact interactions adding a nonlinear contribution to the Schrödinger equation, has been the subject of many recent works [24,25,[30][31][32][33][34][35][36][37][38][39].Although some works focus on the limiting case when interactions dominate (and the kinetic energy generating the quantum pressure becomes negligible), e.g.[40][41][42], effort has also been made to obtain general behaviours of self interacting FDM systems [24,43].Furthermore, constraints upon the possible strength of self-interactions along with the boson mass have very recently been examined [31,32,35,38,44] and, as a result, typical values for a repulsive self-coupling of g ∼ (10 −32 − 10 −26 ) Jm 3 kg −1 are usually quoted within the context of interacting FDM literature, with the exact value depending on the chosen boson mass.Even in the lower ranges, such a small selfinteraction strength can still lead to clear and notable effects, which could be observationally relevant, with [44] presenting favourable observational evidence for the existence of a non-zero self-coupling.Moreover, although attractive FDM solutions are typically unstable beyond a small range of allowed values, potentially leading to collapse and/or the formation of a black hole [24,45], recent work [37] highlights the potentially favourable and desired role of attractive interactions in enhancing smallscale structure formation in large simulations of the cosmic web.Such examples of self interacting FDM would leave observational signatures like an increased number of solitons that form in areas of constructive interference in cosmic filaments, where instabilities are present under the attractive self interaction.
The aim of the present work is to characterize in detail -both analytically and numerically -the effect of repulsive interactions on the size and shape of the underlying virialized solitonic core, and to use such information to accurately predict the oscillations of the core -building on important earlier work on such observables [22,27,[46][47][48][49].Based on this we propose a scheme for uniquely identifying both the self-interaction strength and the boson mass by combining information of soliton shape details and oscillation frequency, which could potentially become available in future observational data analysis.Indeed, oscillations perturb the rotational velocity curves significantly within the soliton and in the region immediately adjacent to it [50].A statistical evaluation of the gravitational heating of galaxies could therefore give a preliminary indication of a relationship between frequency and other soliton parameters.
Firstly, building upon earlier work based on a Gaussian solitonic core of some interaction-dependent effective width [24], and making use of the virialization condition, we derive analytical expressions characterizing the change in shape parameters as a function of interaction strength in the context of the empirical model of solitonic cores [2,3,27,51].Further introducing a characteristic interaction strength marking the crossover between weakly-interacting and strongly-interacting regimes -defined by equating the interaction energy to the quantum kinetic energy (see also [24,34,42]) -we are able to obtain analytical results for appropriately scaled soliton parameters (radius, peak density), cast only in terms of the ratio of the interaction strength to such characteristic value; such results smoothly interpolate between the previously studied non-interacting and strongly-interacting (Thomas-Fermi) limits.
To test the accuracy of our extended analytical predictions, we also perform detailed spherically symmetric numerical simulations of isolated soliton cores generated through the very efficient imaginary time propagation method widely used in cold atom studies [52] (yet scarcely seen in cosmological works [53]).We compare and contrast our virialized solitonic shapes both in the absence and presence of repulsive self-interactions to the Gaussian [24], Thomas-Fermi [24] and empirical [2] profiles, with our findings providing a direct bridge between all such limiting cases.
Moreover, our scaled analytical formulae are numerically confirmed to be independent of the values of the bo-son mass (m) and total soliton mass (M ), at least within the probed ranges 10 −22 eV/c 2 ≤ m ≤ 10 −20 eV/c 2 and 10 6 M ⊙ ≤ M ≤ 10 8 M ⊙ .Such an analysis facilitates a critical test of the typically-used empirical profile, and provides a direct relation of the scaling of central core density and soliton width as a function of scaled interaction strength, smoothly interpolating between relevant previous findings.
Next we evaluate the oscillation frequency of the core, building upon the early analysis of Chavanis [24], but based here not on a Gaussian approximation to the soliton, but on the more accurate empirical profiles for which we explicitly derive a relationship between the peak density, ρ 0 , and the soliton oscillation frequency, f -thus supplementing previous related works based solely on numerical fitting [22,27,[46][47][48][49]. Implementing a perturbative approach on our accurately virialized numerically generated solitons, we obtain accurate numerical predictions for these oscillations which we contrast to our analytical predictions, and simpler previous ones.
Our numerical analysis of soliton shape (peak density, width) and oscillation frequency point to the existence of a degeneracy between boson mass and self interaction strength in determining the resulting soliton parameters.We discuss the implications of such findings for uniquely determining both interaction strength and boson mass from potential observational signatures.Through this analysis, we highlight that care is necessary when imposing constraints on the boson mass in FDM, as the addition of a free parameter in the self-interaction strength could significantly widen the allowed parameter space for the boson mass.We thus indicate that strong constraints such as those discussed in [21,27] should likely be revisited with this in mind.
The paper is structured as follows: After reviewing the governing equations of FDM and common soliton profile shapes discussed in the literature (Sec.II), we critically revisit the shapes of virialized soliton profiles (Sec.III).
Here we review the non-interacting results (Sec.III A), and focus on the modifications imposed by repulsive selfinteractions, which provide new extended relations in terms of a characteristic interaction strength approximately separating the weakly-and strongly-interacting limits: this is obtained by matching quantum kinetic and interaction energies (Sec.III B).Our numerical procedure and results are then shown (Sec.III C), with soliton oscillations discussed both analytically and numerically (Sec.IV).By inverting relations for soliton core density and oscillation frequency parameters previously obtained in terms of a fixed boson mass m, we obtain implicit relations for m(g) as a function of variable interactions for any radius/oscillation frequencies (Sec V), pointing to emerging degeneracies, whose inferred observational constrain implications are discussed (Sec.V B).Brief conclusions are then given in Sec.VI.

II. THE GROSS-PITAEVSKII-POISSON SYSTEM
In the framework of FDM, dark matter is composed of ultra-light bosonic particles and can be described by a classical complex field ψ(r, t).We will take its mass density to be represented by ρ = |ψ| 2 and the energy functional for the field can be written as where Θ, W and U are the kinetic energy, the gravitational energy and the interaction energy respectively.The Newtonian potential V(r,t) is a functional of ψ and obeys the Poisson equation, where the average mass density, ρ, has been subtracted as is necessary for the consistent description of continuous systems under Newtonian gravity [24,54].The self-interaction strength, g, is considered as a two-body contact (short-range) interaction which can be effectively associated with the s-wave scattering length a s through the relationship g = 4π 2 a s /m 2 .By varying Eq. ( 1), one arrives at a non-linear Schrödinger equation, a Gross-Pitaevskii type equation of the form + gρ(r, t) + mV (r, t) ψ(r, t) .
(3) which, along with equation (2) constitute the Gross-Pitaevskii-Poisson Equation system (GPPE) for the description of FDM.By applying the Madelung transformation, ψ = √ ρ exp(iS), and separating real and imaginary parts, we obtain the continuity equation, in which u = ∇S/m defines the velocity and the equation of motion takes the form of the generalised Bernoulli equation The last term of Eq. ( 5) is the quantum pressure, which originates from the uncertainty principle.In the noninteracting case (g = 0), equilibrium is achieved when the repulsion due to the quantum pressure balances the attraction due to gravity [55].
The virial theorem in this case states that It is customary to decompose the kinetic energy into classical and quantum parts, respectively defined by In equilibrium u = 0, and therefore the classical kinetic energy is zero, but the quantum kinetic energy, Θ Q , still contributes to the energy budget.The equilibrium system is static and obeys the virial theorem where only the quantum kinetic energy plays a part, thus

A. Common Soliton Profiles
A characteristic feature of the Gross-Pitaevskii-Poisson system is the coherent core-like soliton structure at the centre of the gravitational well of the halo [13,18,21,56].We discuss below characteristic expressions that are frequently used in the literature to describe the spatial profile of the soliton in various limits.
Firstly, in the absence of self-coupling (g = 0), the gravitational attraction is solely balanced by the quantum pressure term.In this limit, the soliton density is usually approximated by the empirical profile [2] where ρ 0 is the central (peak) core density and r c represents the core radius, defined as the radius at which ρ(r c ) = ρ 0 /2: this yields the exact value λ = 2 1/8 − 1 ≈ 0.091.We note here that while the above profile automatically satisfies the virial theorem 2Θ Q + W = 0, in order to ensure this it is in fact critical to use the exact (as opposed to the approximate) value of λ quoted above (see also subsequent related comments).We also stress here that ρ 0 and r c are not independent parameters, and a central density-radius relation has been reported [2,27,28,48,49] in the form Such an expression can also be obtained from variational energetic considerations [54].
Secondly, the opposite, strongly-interacting limit, corresponds to a hydrostatic equilibrium in which the gravitational attraction is balanced by repulsive interactions [24].In this limit, the solitonic core is instead approximated by a Thomas-Fermi profile [57][58][59][60] of characteristic spatial extent R T F .
For consistency with the notation used for the empirical profile, we introduce here the spatial extent, r T F c , corresponding to the point where the density drops to half its central value, i.e. ρ T F (r T F c ) = (1/2)ρ T F (r = 0), thus yielding More generally, but only accounting for global shape features without attention to detailed radial dependence, the soliton can also be approximated by a Gaussian [24]: while only approximate, and a poor description for both the g = 0 and g → ∞ limits, such description has two obvious benefits: (i) it facilitates an easy qualitative interpolation across these two limiting cases in the context of an effective interaction-dependent width, and (ii) it is easily amenable to analytical calculations.The form of such a Gaussian is The different arising density profiles and a comparison with our more accurate numerical results (discussed in Sec.III C) are shown in Fig. 1 across the entire range of interactions, with g * denoting a characteristic interaction strength (defined later) marking an approximate crossover from the non-interacting to the stronglyinteracting limits.

B. Generalised Soliton Profile Ansatz
In order to more accurately characterize the role of interactions on the soliton shape, we proceed by constructing a generalized ansatz.On physical grounds, we anticipate for such a general density profile to depend on the radial coordinate r, and the interaction strength g, only through the dimensionless ratios r/r c and g/g * , where g * is a characteristic interaction strength to be defined later.We thus introduce a generalised ansatz of the form where Following Chavanis [24] we start our analytical studies by reducing the expression for the mass and various system energies in terms of relevant physical variables Empirical [Eq.( 9)] Gaussian [Eq.( 13)] Thomas-Fermi [Eq.( 12)] Analytical: FIG. 1. Dependence of the shape of density profiles on interaction strength for characteristic examples, and comparison between numerical (generated using imaginary time propagation to find the ground state solution) and analytical predictions in such regimes, with interaction strength increasing from top to bottom curves.Compared to the noninteracting limit (top green line), we also show a profile with self-interactions of an adequate strength (g * ) to impact the properties of the soliton (orange line), and the corresponding profile in the Thomas-Fermi regime of strong self-interactions.The empirical profile is given in Eq. ( 9), the Thomas-Fermi and Gaussian and profiles are respectively defined by Eq. ( 12) and Eq. ( 13).In all cases, the profiles are normalised to M = 1 × 10 7 M⊙, and with m = 2 × 10 −22 eV/c 2 .The characteristic interaction strength g * corresponds to the value marking an approximate crossover from weakly-interacting to strongly-interacting, and is defined in Eq. ( 34) by balancing interaction and quantum kinetic energies -see subsequent discussion in the text.and appropriate dimensionless integrals, whose values can be calculated based on our generalised ansatz.Specifically we obtain the following expressions: Total Mass: Moment of Inertia: Quantum Kinetic Energy: Gravitational Energy: Interaction Energy: In the above expressions, η(Γ), α(Γ), σ(Γ), ν(Γ) and ζ(Γ) are parameters depending on the dimensionless interaction strength, obtained from a dimensionless integral directly related to the shape of the density profile of the soliton in the presence of self-interactions.We shall henceforth refer to such parameters as shape parameters.Such parameters were first introduced in the pioneering work of Chavanis [24] who used a Gaussian approximation to the soliton density profile, given by Eq. ( 13).We calculate the specific shape parameters for the three aforementioned profiles in Table I.
As evident, the values of these shape parameters are rather sensitive to the exact profile shape.Moreover, we emphasize that calculations of soliton static and dynamical parameters (e.g.radius, peak density, oscillation frequency) will be found to exhibit a high degree of sensitivity in changes to the shape of the profile -and thus the shape parameters -already in the non-interacting limit, an effect amplified when interactions are taken into consideration.

III. VIRIALIZED SOLITON PROFILES
A. Non-Interacting Virialized Profiles (g = 0) In the absence of interactions (g = 0), the Virial theorem between the quantum kinetic and the gravitational energies can be used, jointly with the above energy expressions to establish a direct relationship between peak density and core radius.
In terms of such notation, we thus obtain from Eqs. ( 21), ( 23) and ( 26) an expression for the soliton radius, r c in the non-interacting limit, in the form Substitution of this into the earlier reduced expression for the soliton mass Eq. ( 17), provides a relation between r c and ρ 0 , namely Interestingly, an r c ∝ ρ −1/4 0 relationship along with a coefficient of proportionality were numerically obtained from fitting [2,28] without any explicit reference to virialization.Our present analysis extends such numericallybased result by an analytical derivation of the coefficient of proportionality based on the empirical profile of Eq. ( 9), explicitly showcasing its dependence on a range of physical parameters.
It is now possible to rephrase the equation for soliton mass M in terms of a single parameter only, for example in terms of the peak density ρ 0 , in the form Using this relation, the empirical profile in Eq. ( 9) for the non-interacting limit can be reformulated as As the empirical profile provides the best approximation to the ground state of the GPPE in the non-interacting limit, Eqs. ( 27) and (30) reveal that, for any given pair of values for the soliton and boson masses (M, m), one can ab initio generate a virialised soliton profile for the noninteracting case, with shape parameters taken from the empirical profile.Such g = 0 profile is shown in Fig. 1, demonstrating the good agreement with the numerical solution of the GPPE discussed in the next section.The profile-specific shape parameters for the three commonly used profiles are given in Table .I, both in the absence and presence of interactions.

B. Modifications due to Interactions
Next, we consider the effect of repulsive selfinteractions (g > 0) on the soliton profile.These contribute to the system's energy via Eq.( 25).
In Sec.III A we have used the Virial theorem to obtain an expression for the core radius, r c , in the noninteracting limit, Eq. ( 27), expressed in terms of the non-interacting shape parameters σ 0 and ν 0 .Using now instead the generalized form of the virial theorem in the  I. Values for the various shape parameters arising from the relevant energy integrals depending on the exact profile shapes: Such values are obtained analytically for the cases of a Gaussian profile (an analysis already performed in Ref. [24] from which the values in the 3 rd column are quoted), and evaluated in this work also analytically for the empirical profile of Eq. ( 30) (4 th column).These are compared and contrasted to results from the numerical simulations in the non-interacting regime g = 0 (5 th column), and over a broad range of relevant interaction strengths smoothly interpolating between the non-interacting and strongly-interacting limits, showing the change in the shape parameters as a function of interactions (penultimate column).
To highlight the observation that the values of some shape parameters decrease (η, α), while others (σ, ν, ζ) increase with increasing interaction, we depict the limiting high interaction strength values by blue.For comparison, the final column shows the analytically predicted values in the strongly interacting Thomas-Fermi regime, based on the Thomas-Fermi profile of Eq. (12).A graphical representation of the dependence of such shape parameters on scaled interaction strength Γ is shown in Fig. 2.
presence of interactions, Eq. ( 8), we obtain the more general expression in the interacting case, namely (31) Note that as the shape of the soliton changes with increasing value of the self-coupling g, all shape parameters now become a priori unknown functions of g.However, we will show that such parameters become universal in terms of the dimensionless interaction strength Γ = g/g * , so in the above equation and henceforth we denote these as {η(Γ), α(Γ), σ(Γ), ν(Γ), ζ(Γ)}.
We note here that such an equation, but with constant values for the shape parameters based on a Gaussian ansatz, was already obtained in Ref. [24].Here, we have generalized such expression to reflect the changing soliton profile due to repulsive interactions through the Γdependent shape parameters.The above relation therefore describes an a priori unknown dependence of r c on g which is more complex than that implied solely by the explicit appearance of g on the rhs.
The corresponding expression for the peak density can be obtained from Eq. ( 17) as We can gain some immediate analytical insight for both core radius and peak density by defining corresponding limiting functions in the non-interacting and stronglyinteracting limits.To do this it is useful to first define an expression for the characteristic interaction strength, g * , separating these two limiting cases.

Characteristic Interaction Strength
A useful estimate of the (relative) strength of interactions can be obtained by defining a characteristic interaction strength scaling, g * , as the interaction strength value at which the interaction energy equals the quantum kinetic energy2 , i.e.
As the dependence of the shape parameters ζ, ν and σ on Γ = g/g * is a priori unknown, the above expression is in fact an implicit equation for g * .We can however simplify the procedure and obtain a well-defined explicit definition for the characteristic interaction strength g * by making the heuristic, but convenient, choice to evaluate the shape parameters using the empirical profile, Eq. ( 9), in the non-interacting (g = 0) limit.The important benefit of our choice is that it will allow us to present results in a universal manner with m and M scaled out.We would expect the crossover between the non-interacting and strongly-interacting regimes to take place around g ∼ O(g * ), i.e.Γ ∼ 1.
Based on such a choice, we define g * as We note that, although the actual value of g * would be slightly shifted following a different choice for the evaluation of the shape parameters, such choice would have no qualitative influence on the findings reported throughout this work.The effects of repulsive self interactions on the resulting soliton ground state are visible in Fig. 1.As one dials up the strength of the (repulsive) self interaction, the core begins to expand and the density drops, as already visible for g = g * (orange line).In this regime, neither the empirical profile nor the Thomas-Fermi profile adequately replicate the form of the density profile.However, once the interaction strength is large enough and the quantum kinetic energy is negligible in comparison to the interaction energy, the system is well modelled by the Thomas-Fermi profile, as shown there in the case of g = 50g * .

Analytical Expressions in Limiting Cases
The above expression for g * allows us to obtain useful analytical expressions for the leading order dependence on the scaled interaction strength in the weakly and strongly interacting limits: in each such case we simplify the expression of Eq. ( 31) by replacing the Γ-dependent shape parameters by the corresponding constant parameters obtained analytically in the limiting cases of zero and strong interactions (Thomas-Fermi).
In the weakly interacting case, we thus introduce where subscripts 0 are used to denote corresponding shape parameter values at g = 0.This can now be rewritten in terms of the dimensionless interaction strength, Γ, in the more compact form where r N I c describes the non-interacting value given previously by Eq. (27).
Correspondingly, in the strongly-interacting regime we introduce the notation where r T F c corresponds to the value of r c in the limit g → ∞, and g T F * is defined as in Eq. ( 34), but based on shape parameters evaluated analytically for the Thomas-Fermi profile, i.e.
The values of the shape parameters corresponding to the two limiting soliton profiles are given in Table I (columns 4 and 7 respectively for non-interacting and Thomas-Fermi limits).
Similar to the above discussion we define the limiting expressions for the interaction-strength dependence of the peak density, ρ 0 (Γ), in the non-interacting case; specifically we find whereas in the Thomas-Fermi limit In each above case, the shape parameters are evaluated at the corresponding limits.We also note here that

Crossover between Weakly-and Strongly-Interacting Regimes
The soliton size in the presence of interactions, r c (g), increases monotonically with g, always lying between these two limiting curves, i.e. r 0 c (g) ≤ r c (g) ≤ r T F c (g); this will be confirmed by numerical simulations below.Its exact behaviour can be directly extracted from Eq. (31) along with a numerical determination of the shape parameters σ(Γ), ζ(Γ) and ν(Γ) via the more compact relation where we have introduced a generalised dimensionless interaction strength parameter C(Γ) which also accounts for interaction-induced shape effects via The corresponding dependence of the central density on interactions is then given by where η(Γ) can again be determined from numerically obtained soliton profiles.
In the following subsection we perform a detailed numerical study which will confirm the validity of the above semi-analytical formulae.Indeed, we will see below that a numerical determination of the shape parameters, along with Eqs. ( 42) and (44), describes the size and central density of a soliton of mass M , made up of bosons of mass m, very well across all values of the self-coupling g from the non-interacting to the strongly interacting regimes.
C. Numerical Approach: Imaginary Time Propagation as a Tool for Ground State Profiles

Dimensionless GPPE
In order to numerically solve the GPPE we revert to dimensionless form by scaling all physical variables to appropriate values.This is done in terms of a time scale characteristic of an object crossing a uniform density configuration of some reference density ρ ref and the corresponding energy and length scales for some particle number N = M/m.Setting t = T t ′ and r = Lr ′ , the resulting dimensionless GPPE takes the form and where the interaction strength is written as g = g ′ √ Gρ ref /ρ sys , the gravitational potential as V = V ′ √ Gρ sys /m √ ρ ref and the wavefunction has been scaled to the mean density ρ sys within the computational box through the dimensionless form ψ ′ = ψ/ √ ρ sys .We choose the mean reference density of the system as ρ ref = 1.5 × 10 3 M ⊙ /kpc 3 which is commensurate with the current cosmic density.We will henceforth drop all primes for the sake of simplicity.
As the true ground state solution to the GPPE (which is what we are seeking numerically) is a spherically symmetric density distribution, we further simplify the problem to be solved numerically to a spherically symmetric form (thus disallowing any angular dependence/asymmetries to creep up also in subsequent dynamical analysis).A further significant benefit in doing so is that it makes the required numerical resolution of the soliton feasible, without unnecessarily strong demands on computational memory and time.
Using the further substitution ψ = φ/r, we obtain the radial (one-dimensional) GPE in the form which we solve using the implicit midpoint method.The Poisson equation to which the GPE is coupled to now takes the simpler form where we have made the substitution V = V/r and dropped the ρ term which is no longer necessary due to the fact that we are using a spherically symmetric regime with boundary conditions that go to zero.(Note that subtraction of the mean density would have been necessary if we were to include periodic boundary conditions.) The computational speed up acquired by translating to a spherically symmetric system and using this method is significant and allows for the study of higher mass solitons with larger resolution, where 3D approaches may suffer from a lack of computational power or simply from a lack of memory in being able to store large 3D arrays.

Imaginary Time Propagation & Numerical Ground State
To generate the ground state solution of the system, we employ a method widely used in the study of cold atomic condensates [52] by propagating our system in imaginary time.By setting dt → −idt, we reformulate the system in such a way that a 'forward' propagation in imaginary time results in all higher energy eigenstates decaying exponentially at a rate proportional to the energy of each eigenstate.As a result, the ground state decays comparatively more slowly than all other (excited) states.Thus, by additionally enforcing a renormalization of the particle number within the system at the end of each dt propagation, we ensure that the norm of φ remains constant through such a procedure.Such combination of imaginary time propagation and renormalization achieves the effect of gradually extracting energy from the system.As a direct consequence of this, the imaginary time propagation enforces a gradual transition to the ground state.
To obtain the ground (virialized) state to the desired accuracy, we need to introduce some numerical tolerances.In our calculations we assume such state is reached when the total energy of the system has converged in two consecutive timesteps to within one part in 10 5 .To ensure the size of our numerical grid has no influence on our findings, we also monitor the value of the density at r = ∆r, the inner-most point of the computational sphere, and also the ratio ρ(∆r)/ρ(L), where L is the edge of the computational sphere.By imposing a large value for this ratio (here ρ(∆r)/ρ(L) > 10 40 ), we ensure that the computational sphere is large enough in size to adequately contain the soliton without the boundaries interfering with its shape and dynamics.Moreover, in our numerical calculations, we need to balance the need for good accuracy of the soliton core (i.e.high spatial grid resolution) -which controls the extent to which the virial theorem is numerically satisfied -and the computational time required (which increases with increasing number of grid points).In practice we  found that the description can be sufficiently accurate upon including at least 5 grid points between r = 0 and r = r c , although a higher number of grid points will further enhance numerical accuracy.
Shape parameter values are extracted from simulations by performing the relevant numerical energy integrals on the generated ground state profile and then using the peak density and core radius of the numerical soliton to calculate them.Although all shape parameters should be uniquely identified in the case of g = 0, numerically we find a slight variance of up to few % in their values, as seen in the fifth column of Table I.Such variance accounts for a slight change due to a different spatial resolution and computational box size.However, in all simulated non-interacting cases, there is a clearly identified value for each shape parameter, with such values being in accordance with our analytical values calculated from the empirical formula (fourth column of Table I).The numerical average for each non-interacting shape parameter value is plotted on the left of each subplot in (subsequent) Fig. 2 as a filled green point at Γ = 0. Error bars are not visible in most cases, due to the small standard deviation in results.

Numerical Results for Variable Self-interaction Strength
Ground states were generated for a range of g values, for several different configurations of m and M .Shape parameter values were extracted from these numerical ground states.Such behaviour is shown for all 5 shape parameters (η, σ, ν, α and ζ) in Fig. 2. As evident, all such parameters follow a clear monotonic transition from the non-interacting limiting case (leftmost part, green points) to the Thomas-Fermi case (rightmost part), with some of them decreasing (η, α) and others increasing (σ, ν, ζ) with increasing (scaled) interaction strength Γ.
Interestingly, we find that the dependence of all such shape parameters on Γ can be excellently fit by a function of the form where u denotes any of these shape parameters, subscripts 0 and T F denote the corresponding analytical values in each regime, and a and b are fitting constants.
The specific chosen values of these fitting constants for each shape parameter are given in Table II, and vary slightly for each shape parameter.Nevertheless, in all cases, both a and b are of order unity.An important comment needs to be made regarding the shape 2. Dependence of the value of shape parameters on interaction strength, scaled to the characteristic value g * defined by Eq. (34).Such dependence is reproduced by a twoparameter fit based on Eq. ( 52) which has as limiting cases the non-interacting limit at Γ = 0 (left), and the Thomas-Fermi limit at Γ → ∞ (right), whose exact values for each parameter are given in Table I.The corresponding predictions of Ref. [24] based on a Gaussian approximation (reported in Column 3 of Table I) are shown here by the horizontal dashed line.Filled green points at the g = 0 line represent the averaged values of the corresponding shape parameters from numerical simulations (see Table I, fifth column).
parameter σ, which originates from the quantum pressure.Given that the Thomas-Fermi regime is defined by a negligible quantum pressure, the analytical value computed from the quantum kinetic energy integral for the Thomas-Fermi profile loses meaning.Moreover, the deviation of the numerical values for σ from the fit in Fig. 2 for large values of Γ > 10 3 , can be attributed to the much steeper decrease of the Thomas-Fermi profile towards a zero value (in stark contrast to the smoother decrease of the numerically generated profiles).
Next, we comment on the validity of the Gaussian approximation for different values of g.As is evident from Fig. 2, the Gaussian values taken from Ref. [24] always fall within our obtained numerical range: in fact, such points typically lie close to the midpoint between our presented analytical limiting values (except for η where the Gaussian value effectively corresponds to the Thomas-Fermi limit), and, as such, are good approximate values in the crossover regime.Note however that although such values are never off by more than a factor of 3, such analytically convenient choice of values cannot capture the true variable nature of the shape parameters.
Having fully quantified the universal dependence of the shape parameters on the dimensionless interaction strength Γ = g/g * , we can now discuss the corresponding, and observationally-relevant, dependence of the key soliton profile parameters of peak density and soliton radius on interactions.To present this in the most general manner, we consider the dependence of the rescaled central soliton density ρ 0 (Γ)/ρ 0 (0) and core radius r c (Γ)/r c (0) as a function of the rescaled self-interaction strength Γ = g/g * .The results of our numerical simulations are respectively shown in Fig. 3(a)-(b).
Analytical expressions for these curves can also be obtained as follows: The scaled soliton width dependence on Γ can be obtained through Eq. ( 42) as (53) The corresponding dependence of the peak density can be obtained from Eq. ( 44) via Inputting our numerically-evaluated universal shape parameter dependence on Γ based on Eq. ( 52) for σ(Γ), ν(Γ), ζ(Γ), and η(Γ) into the above equations yields the solid black lines in Fig. 3(a)-(b): such curves are clearly found to lie directly on top of all our numerical points for different (M, m) combinations.To clarify such observations, in each of the two cases we also plot the corresponding fully-analytical weaklyinteracting (r 0 c (Γ), ρ 0 0 (Γ)) [solid grey lines] and stronglyinteracting (Thomas-Fermi) limits (r T F c (Γ), ρ T F 0 (Γ)) [dotted orange lines], defined respectively by Eqs. ( 36), ( 39) and ( 37), (40).These curves delimit the region where all our data points lie and which we indicate by green shading.
In both cases we clearly see a transition from the noninteracting behaviour for small Γ to the Thomas-Fermi behaviour for large Γ, with such transition indeed occurring around Γ = 1, i.e. near the characteristic interaction strength g * of Eq. ( 34); such value has been highlighted by a vertical dashed red line for easier visualization.To make such transitions between the two limiting cases more transparent, we further highlight the weaklyand strongly-interacting limits in subplots (i) and (iii), respectively found to the left and right of the main plot.
We have thus fully characterized our soliton ansatz of Eq. ( 14) in terms of a universal dependence on the (scaled) interaction strength, with such behaviour obtained semi-analytically (by a combination of analytical expressions combined with the universal dependence of the shape parameters on Γ) and confirmed by full numerical simulations.

IV. SOLITON OSCILLATIONS
The simplest (lowest-energy) excitation of the soliton core comes in the form of a radial oscillation.The rate of such an excitation depends on the total mass of the soliton, the mass of the constituent boson and-in the interacting case-the interaction strength [24].Soliton core oscillations could in fact lead us to signals of FDM's existence in astronomical observations, as well as allow us to place constraints on the range of possible boson masses [21,27], and they may play a part in gravitational heating of galaxies.It is also possible that core oscillations may drive dynamical behaviour in the entire halo structure due to interference [50].
In this section we examine the dependence of the frequency of such small radial oscillations on the selfcoupling strength parameter in terms of the dimensionless ratio Γ.This will allow us to overlay data from a variety of solitons on a universal curve.

A. Isotropic oscillation
A proper study of the soliton's oscillations would involve a thorough analysis of its normal modes, an investigation that we leave for an upcoming publication.Here, as in [24], we will employ a perturbative method for studying the radial oscillatory dynamics of the soliton, making the further assumptions that: (i) small oscillations can be described by making r c , the parameter that sets the scale of the solitonic profile Eq. ( 14), timedependent, i.e. r c → r c (t), and (ii) the resulting velocity field takes the isotropic form As we now show, these two assumptions can satisfy the continuity equation, Eq. ( 4), for a specific form of the function f (t).34).All plots show numerically simulated data for 3 different (M , m) combinations and a continuous solid black line obtained from our analytical equations ( 53), ( 54), (76) using as input the Γ-dependent shape factors shown in Fig. (2).Panels (ii) [middle] reveal the entire crossover, with (i) left and (iii) right panels respectively highlighting the weakly-interacting and strongly-interacting limits.The two limiting boundaries of the green channels highlight the range of accessible predicted values, bounded from above and below by the respective characteristic universal curves: specifically, the non-interacting grey line is constructed from Eqs. ( 39), ( 36), (76)), with g * [Eq.(34)] computed using the non-interacting shape parameters; moreover, the dashed black Thomas-Fermi line arises from Eqs. ( 39), ( 36), (76) with g * [Eq.(34)] calculated instead using the Thomas-Fermi shape parameters.The vertical purple line highlights the characteristic interaction strength value g = g * (corresponding to Γ = 1).Note that a cross-over between the limiting non-interacting and Thomas-Fermi lines is in fact present in the case of the oscillation frequency [(c)(ii)], despite this being largely obscured by the very narrow accessible channel.We highlight that the excellent fit of our semi-analytical predictions (solid back lines) through all our numerical data demonstrates the importance of the correct incorporation of the universal variation of shape parameters on scaled interaction strength.Remarkably, such an approach accurately predicts the peak density, radius and frequency, even in the transition region between non-interacting and Thomas-Fermi limits.
Although the form of the soliton profile ρ(r) is not known for g = 0, we may make use of our ansatz ρ(r) = ρ 0 ϕ(r/r c , Γ), further assuming that the profile becomes time dependent only via r c → r c (t).We thus obtain and where ϕ ′ denotes a derivative w.r.t. the function's spatial argument.Clearly, the continuity equation is fulfilled if and therefore describing the oscillation by assuming that the whole profile evolves with time via r c (t) and with the velocity profile is consistent witn mass conservation.This allows us to calculate the classical kinetic energy, Θ C , as Therefore, solitons can exhibit oscillations which physically manifest as a uniform expansion and contraction of the soliton profile.The soliton mass M is constant in time and therefore the peak density must also experience oscillations as ρ 0 (t) ∝ r −3 c (t) -see Eq. ( 17).

B. Analytical Oscillation Frequencies
The total energy of the interacting system can be written as which can be thought of as an integral of the motion for the dynamical equation Hence, the fundamental, breathing oscillation mode of the soliton can be described by the one-dimensional motion of a non-relativistic, Newtonian particle moving in the potential V [24] defined above.A static solution corresponds to a g-dependent equilibrium radius, r * c , satisfying By performing a small perturbation about this equilibrium radius r c (t) = r * c + ε(t) in (66), we find where, to avoid notational clutter, we henceforth drop the * superscript.At this point we remind the reader that, as V (r c ) depends on both explicitly and implicitly via the shape parameters σ = σ(Γ), ν = ν(Γ) and ζ = ζ(Γ), such equilibrium value depends on g: in fact it corresponds to the g-dependent equilibrium value given by Eq. ( 31).We thus infer a frequency where we remind the reader that I = αM r 2 c is the moment of inertia.Accounting for the interacting virial condition [Eq.(8)] and using the expressions for the shape parameters we can rewrite this as (70) Note that r c is now a function of the self-coupling g.Such an expression -but with constant shape factorshas been previously analytically derived in Ref. [24] in the context of a Gaussian profile.Here we generalize this to include a broader range of profiles, with numerical differences hidden within the shape parameters.
In the g → 0 limit and by making use of Eq. ( 28), we can write the oscillation frequency in a form which is solely dependent on the peak density of the soliton, with the shape parameters ν 0 , η 0 and α 0 for the Γ = 0 empirical profile, Eq. ( 9), given in Table I.This form agrees with the f ∝ ρ 1/2 0 relationship expressed in previous literature [22,27,[46][47][48][49]. Specifically, we find while in e.g.[22] the coefficient is 10.94, a value obtained from analysing oscillations of the central soliton in a simulated core-halo system.Very similar values are quoted in the literature [27,[46][47][48][49].The difference is very small and could be due to the fact that values quoted in the literature are mostly extracted from solitons embedded in halos and not in isolation which may not be oscillating in their fundamental frequency only.
It is interesting to note that by looking at Eq. ( 27) and Eq. ( 72) we can compare our results to those from [24], finding that As a result, and perhaps somewhat unexpectedly, although the radius parameters between the Gaussian profile and the empirical profile ( 9) are clearly distinct, the balancing of these and the various shape parameter values in these limiting cases appearing in Table I result in almost perfect cancellation, thus leading to practically the same predictions for the frequency.Therefore, our extended present work confirms that the Gaussian ansatz approach can be considered a remarkably robust approach for analysing the soliton's oscillation frequency.

C. Numerical Results for Non-Interacting Solitons
A small-amplitude perturbation enabling the probing of such lowest-energy system oscillations can be very easily engineered in our numerical simulations, using the following trick: instead of running our imaginary time propagation until full system virialization has been achieved, we choose to terminate such process somewhat earlier to allow for the otherwise almost perfect ground state solution to be left with a 'natural' built-in perturbation.Subsequently propagating such perturbation in the (real) time domain we can extract the frequency of such an oscillatory mode by performing a Fourier analysis of the dynamics of the peak density value.A critical point to emphasise is the fact that these early terminations, or any perturbation of a soliton in fact, will result in the creation of an "atmosphere of bosons" about the soliton -a miniature version of a dark matter halo.Such atmospheres are difficult to minimise, and impossible to avoid as any addition of energy to the ground state results in excited bosons being ejected from the soliton.Therefore, it is crucial that one accounts for this when comparing analytical and idealised formulae to simulations.For that reason, we impose a check in our simulations by analysing the quantity of mass found within the ground state in the temporally averaged soliton across the entire dynamic simulation, and verify that minimal mass loss occurs.Specifically, we make use of the Penrose-Onsager mode to calculate the total mass in the ground state of a dynamically propagated simulation [54].Indeed, we find that in all our simulations, the losses to such atmospheres around the core are at a level of at most 10 −5 of the initial ground state solution and therefore data is not corrupted due to mass leakage from the soliton.
According to Eq. ( 27) the radius, and therefore the peak density, are both functions of boson mass.We can therefore rewrite Eq. (71) in terms of the boson mass and soliton mass, In order to compare results of the empirical and Gaussian profiles to our numerics, we evaluate the above frequency formula, Eq. (75), using the shape parameters obtained from the respective profiles.A detailed comparison between simulation data and the prediction from Eq. ( 75) is shown in Fig. 4 75) for the empirical (and Gaussian) analysis with our numerical data [Fig.4(a)].To better understand how these compare, and what the subtle differences between empirical and Gaussian results (concealed in Fig. 4(a)) may be, Fig. 4(b) plots their scaled residual differences in each case.
This reveals that, despite a very good overall agreement with our numerical simulation data, such results for both empirical (filled black circles) and Gaussian (hollow diamonds) consistently overshoot the numericallyobtained frequency (green line).Motivated by our preceeding detailed shape-parameter analysis, we thus proceed to calculate the shape parameters individually from each numerically generated ground state, and use theserather than the shape parameters from the empirical profile -in Eq. (75): as expected, this semi-analytical procedure (open purple circles) yields a much better agreement with the oscillation frequency extracted from simulations.It is therefore clear that the calculation of the oscillation frequency of a soliton is in fact also sensitive (on the 10% level) to the shape parameters and therefore the shape of the profile.

D. Results for Interacting Solitons
The analytical prediction (70) for the oscillation frequency in the case of g = 0 can be written in terms of the dimensionless parameter Γ as where we remind the reader that σ 2 (Γ) Γ -see Eq. (43).A comparison of the above formula, with  81) and Eq. ( 82) respectively.Here, m0 is defined as the boson mass at Γ = 0 for the respective soliton parameters (rc, f, ρ0).

V. PHYSICAL PARAMETER DEGENERACIES AND OBSERVATIONAL IMPLICATIONS
A. Physical Solitons and m − Γ Degeneracy Our preceding analysis for a soliton of fixed total mass M has been conducted under an implicit further assumption of a fixed boson mass m.As anticipated, this has revealed a dependence of static and dynamical properties on the boson interaction parameter g, which enters as a new parameter appearing both implicitly and explicitly in the equations for the soliton radius r c (g), [Eq.(31)], the peak soliton density ρ 0 (g) [Eq.(32)] and the soliton oscillation frequency f (g) [Eq.( 70)].
To explain this, let us focus below on the soliton radius r c in the absence, or presence of interactions, for which we remind the reader of the main functional dependence of the previous expressions (under the assumption of a fixed M ).These take the respective forms: For g = 0: where A 0 = 2(σ 0 /ν 0 )( 2 /GM ) depends on the noninteracting shape parameters (and the physical constants , G and constant soliton mass M ).For g > 0: where A(Γ) = 2(σ(Γ)/ν(Γ))( 2 /GM ) is defined in terms of the Γ-dependent shape parameters, and for conve-nience we have also introduced an expression C(Γ) ≥ 1, which has convenient limiting cases From Eq. (78), we see that the inter-dependence of parameters r c , m and Γ = g/g * implies that it takes a combination of the values of two of these parameters to fix the third one.Given that the main observationallyrelevant physical quantities are likely to be the soliton mass M (already assumed as constant in above discussion) and the soliton radius, we shall henceforth assume here a fixed value for the soliton radius r c (m, Γ) = R c , and consider the resulting inter-relation between the boson mass m and the scaled boson self-interaction Γ.Thus by inverting Eq. ( 79), we obtain for each value of R c the following dependence of m on Γ This expression shows clearly that fixing r c to the value R c and varying Γ results in a varying m(Γ) curve which becomes a locus of all values of the boson mass and dimensionless self-coupling corresponding to a soliton of fixed radius R c .The same logic can be independently applied for fixed values of each of the peak density and the frequency.In particular, the corresponding equations defining these degeneracy curves read: For a fixed peak density, ρ 0 , For fixed frequency, f , where we have introduced Equations ( 81) and (82) determine the loci of all values of m and Γ which result in solitons of the same central density and same oscillation frequency respectively.In all cases, such degeneracy curves reduce to their corresponding non-interacting limits when we take Γ = 0, being consistent with a single viable boson mass value.Such dependences of the boson mass m on Γ, resulting in (m, Γ) pairs that support solitons of fixed radius, peak density and oscillation frequency, are shown in Fig. 5.Note that the boson mass has been scaled by a reference mass m 0 = m(Γ = 0) to make the plot universal.In this plot, we see clearly that the curves corresponding to fixed soliton radius and fixed peak density overlap, whereas the scaling for fixed frequency reveals qualitative similarity but is distinctly identifiable.
In the above discussion, the mass has been obtained as a function of the dimensionless interaction strength Γ = g/g * for fixed soliton radius, peak density or oscillation frequency.It can be instructive to reconsider the problem in terms of the relation between the boson self-interaction strength g, and the boson mass m, for any fixed value of soliton radius, peak density or oscillation frequency, parameterised through a particular value of Γ.This can be done by writing where m(Γ) depends on r c , ρ 0 or f , as given in Eqs.(80), (81) and (82) respectively.Tracing over the parametric variable Γ, and for each fixed value of radius, peak density or frequency, we can thus obtain a parametric relation between g and m.Such dependence is plotted in the three panels of Fig. 6, in which the fixed value of each such parameter is represented by a different colour line.

B. Observational Relevance
For any single soliton with radius r c and peak density ρ 0 , we have explicitly shown that there exists an entire contour line of boson mass and interaction strength pairings that support its existence within FDM.This rules out basing any constraints on a high-resolution accurate observation of a singular galactic dark matter object as the degeneracy means that no meaningful limits can be established; instead one can only deduce which contour the object must lie on.However, multiple objects of different total soliton mass will have differing contour lines, Scheme for simultaneous identification of boson self-interaction and mass through the numerical overlap of key soliton g − m degeneracy lines for two different galaxies.Shown are heuristically chosen lines of constant radius (solid) and frequency (dashed) for soliton masses of M = 10 8 M⊙ (red) and M = 10 9 M⊙ (blue).Their overlap points to a unique combination of g and m values for such systems.Here we have chosen to show a 0.57 kpc radius and an oscillation frequency of 2.66 Gyr −1 for the 10 8 M⊙ soliton, whereas the 10 9 M⊙ soliton has a 0.23 kpc radius, and an oscillation frequency of 39.86 Gyr −1 .A single underlying FDM model would result in the degeneracy curves from a multitude of solitons with different masses, radii and oscillation frequencies all crossing at the same point (see text).
where g(m) gradients are dependent on their total mass.As a result, when one maps the contour lines of several galactic objects onto the same set of appropriaterly scaled axes, their intersection point will correspond to the true physical value of the boson mass and self-interaction strength which enables the formation of such solitons in our universe.Fig. 7 illustrates this point: First, we consider a specific pair of boson mass and self-interaction strength, m = 2.2 × 10 −22 eV/c 2 and g = 10 −29 Jm 3 /kg which, for a given soliton mass, will yield a specific set of soliton parameters -the core radius, peak density and oscillation frequency.We then note that any points along the corresponding m(g) degeneracy contour will yield an identical solitoni.e. the resulting solitons are observationally indistinguishable for points chosen anywhere along a single degeneracy contour line for a fixed total mass.
We can now design two solitons of different mass, chosen as 10 8 M ⊙ and 10 9 M ⊙ , and plot their degeneracy contour lines in the m − g parameter space, along which the radius and frequency remain constant (for fixed total mass).In this specific case, the 10 8 M ⊙ soliton has a 0.57 kpc radius and an oscillation frequency of 2.66 Gyr −1 , whereas the 10 9 M ⊙ soliton has a 0.23 kpc radius, and an oscillation frequency of 39.86 Gyr −1 .These values were chosen heuristically for the purpose of illustrating clearly how the information of the true boson mass and interaction strength may be extracted from the contour lines.The location of the intersect of these contour lines is the true boson mass and self-interaction strength which accommodates both of these solitons to exist in space.Once again, in our case the intersect values were arbitrarily chosen as an initialisation requirement.In a universe where the dark matter is fuzzy and described by a single scalar field, all degeneracy curves corresponding to candidate observed solitons must cross at the same point (within observational error bars) which would represent the single universal value for the mass m and self-coupling g that correspond to our universe.The existence of the aforementioned degeneracy could have an impact on current constraints on the boson mass which are likely to change if the parameter space is extended to include g. 3 Indeed, various limitations have been placed upon the allowed range of the the boson mass by probing both cosmological scales e.g.[33,61], as well as singular astronomical objects such as the star cluster orbiting the centre of Eridanus II [21,27] or rotational curves and stellar kinematics [62][63][64][65].All these works, apart form [33], do not include a self-interaction and in some cases clear differences between theoretical fits and observational data were identified [64], while in others it was further evident that a single boson mass could not adequately fit theoretical curves to observational data [62,65].On the other hand, [44] reports that a non-zero value of the self-coupling, along with a single boson mass, can fit the rotation curves of the dark matter dominated galaxies in the SPARC database, providing for the first time positive support for FDM solitons with a non-zero value of the self-interaction from rotation curves.The degeneracies discussed here would be relevant for all the above studies which make inferences about the boson mass. 4Clearly, all of the above are very preliminary observational investigations but provide strong motivation for precise, quantitative analysis of the model.

VI. CONCLUSIONS
In this paper we evaluated quantitatively the impact of a non-zero, repulsive self-interaction on the shape of virialized, fuzzy dark matter solitons and some of their key characteristics, namely their central density ρ 0 , radius r c and the frequency f of their small radial oscillations.We achieved this using a general ansatz for the density profile and quantified the change in its shape via the computation of 5 shape parameters, corresponding to dimensionless integrals that are associated to the different components of the soliton's energy, mass and moment of inertia.
By using a dimensionless measure of the selfinteraction, Γ, we mapped the transition of the shape parameters from the non-interacting to the strongly interacting regime with universal curves; all of the numerically generated solitons we explored (generated via imaginary time propagation with different masses M and composed of bosons of different mass m), can be placed on these curves.Knowledge of the shape parameters then informs the transition of ρ 0 , r c and f from the non-interacting to the strongly interacting regimes.Again this transition can be described by universal functions if appropriately scaled quantities are used.As an interesting sideobservation, we found that the accurate numerical determination of the shape parameters also allows for a better semi-analytical prediction of the oscillation frequency for g = 0 compared to the use of either the Gaussian or empirical density profiles.
Our results also led us to the notion of degeneracy curves, representing the loci of all those (m, g) pairs that result in solitons with identical r c , ρ 0 and f .The question then naturally arises of how such degeneracy might be broken via the observation of many possible solitonic objects.It would seem that such notions are highly topical given the increasing effort to place constraints on the parameters of FDM via confrontation with observations of astrophysical objects that could harbour FDM solitons.We leave further exploration of such confrontation to future work.
FIG.2.Dependence of the value of shape parameters on interaction strength, scaled to the characteristic value g * defined by Eq.(34).Such dependence is reproduced by a twoparameter fit based on Eq. (52) which has as limiting cases the non-interacting limit at Γ = 0 (left), and the Thomas-Fermi limit at Γ → ∞ (right), whose exact values for each parameter are given in TableI.The corresponding predictions of Ref.[24] based on a Gaussian approximation (reported in Column 3 of TableI) are shown here by the horizontal dashed line.Filled green points at the g = 0 line represent the averaged values of the corresponding shape parameters from numerical simulations (see TableI, fifth column).

FIG. 3 .
FIG.3.Universal dependence of key static and dynamical soliton properties on scaled dimensionless interaction strength.Shown are [from top to bottom] the cases for (a) peak soliton core density, (b) core radius, and (c) core oscillation frequency, with such quantities scaled to their corresponding non-interacting simulated values, and all self-interactions scaled to the characteristic value g * from Eq.(34).All plots show numerically simulated data for 3 different (M , m) combinations and a continuous solid black line obtained from our analytical equations (53), (54), (76) using as input the Γ-dependent shape factors shown in Fig.(2).Panels (ii) [middle] reveal the entire crossover, with (i) left and (iii) right panels respectively highlighting the weakly-interacting and strongly-interacting limits.The two limiting boundaries of the green channels highlight the range of accessible predicted values, bounded from above and below by the respective characteristic universal curves: specifically, the non-interacting grey line is constructed from Eqs. (39), (36), (76)), with g * [Eq.(34)] computed using the non-interacting shape parameters; moreover, the dashed black Thomas-Fermi line arises from Eqs. (39), (36), (76) with g * [Eq.(34)] calculated instead using the Thomas-Fermi shape parameters.The vertical purple line highlights the characteristic interaction strength value g = g * (corresponding to Γ = 1).Note that a cross-over between the limiting non-interacting and Thomas-Fermi lines is in fact present in the case of the oscillation frequency [(c)(ii)], despite this being largely obscured by the very narrow accessible channel.We highlight that the excellent fit of our semi-analytical predictions (solid back lines) through all our numerical data demonstrates the importance of the correct incorporation of the universal variation of shape parameters on scaled interaction strength.Remarkably, such an approach accurately predicts the peak density, radius and frequency, even in the transition region between non-interacting and Thomas-Fermi limits.

21 FIG. 6 .
FIG.6.Curves of constant radius, peak density and oscillation frequency (the value of which is shown by the colourbars), revealing the allowed parameter space of repulsive selfinteraction strength and boson mass based on such chosen parameters.These graphs can be extended towards the heavier boson mass regimes.The figure plotted for a soliton mass M = 10 8 M⊙.
FIG. 7.Scheme for simultaneous identification of boson self-interaction and mass through the numerical overlap of key soliton g − m degeneracy lines for two different galaxies.Shown are heuristically chosen lines of constant radius (solid) and frequency (dashed) for soliton masses of M = 10 8 M⊙ (red) and M = 10 9 M⊙ (blue).Their overlap points to a unique combination of g and m values for such systems.Here we have chosen to show a 0.57 kpc radius and an oscillation frequency of 2.66 Gyr −1 for the 10 8 M⊙ soliton, whereas the 10 9 M⊙ soliton has a 0.23 kpc radius, and an oscillation frequency of 39.86 Gyr −1 .A single underlying FDM model would result in the degeneracy curves from a multitude of solitons with different masses, radii and oscillation frequencies all crossing at the same point (see text).