Bayesian analysis for rotational curves with $\ell$-boson stars as a dark matter component

Using Low Brightness Surface Galaxies (LBSG) rotational curves we inferred the free parameters of $\ell$-boson stars as a dark matter component. The $\ell$-boson stars are numerical solutions to the non-relativistic limit of the Einstein-Klein-Gordon system, the Schr\"odinger-Poisson (SP) system. These solutions are parametrized by an angular momentum number $\ell = (N-1)/2$ and an excitation number $n$. We perform a bayesian analysis by modifying the SimpleMC code to perform the parameter inference, for the cases with $\ell = 0$, $\ell = 1$ and multi-states of $\ell$-boson stars. We used the Akaike information criterion (AIC), Bayesian information criterion and the Bayes factor to compare the excited state ($\ell$=1) and the multi-state case with the ground state ($\ell$=0) as the base model due to its simplicity. We found that the data in most galaxies in the sample favours the multi-states case and that the scalar field mass tends to be slightly bigger than the ground state case.


I. INTRODUCTION
Since the first high resolution observations of M31 by Vera Rubin and Kent Ford in 1970 [64], the galaxy rotation curves have become an important observable for testing dark matter models. Most recently the authors in [65], and references therein, have shown that Low Surface Brightness galaxies (LSBG) are suited candidates to test dark matter models due to their low visibility in the optical and HI photometry, therefore it is reliable to assume that rotational curve dynamics in those galaxies depends mainly on the dark matter component. Throughout the years there has been a plethora of proposals about the nature of this mysterious component, being the scalar field dark matter (SFDM) a candidate that have been proved to be viable and that display several advantages over the standard candidate Λ Cold Dark Matter (ΛCDM). For instance, one point in favor of the SFDM is that it does not suffer of the so-called "cusp/core" problem, -which states that the predicted density profiles, from simulations with ΛCDM, increase steeply and hence producing a "cuspy" dark matter distribution in small radii, while observations of dark matter density profiles for most of the dwarf galaxies indicate otherwise: they display flat central "cores" [66,67] -. Another possible discrepancy between observations and simulations with ΛCDM is the missing satellite problem (MSP). Simulations indicate the presence of more dwarf galaxies orbiting around a galaxy with the characteristics of the Milky Way, than those we are able to observe around us. In the past years some discussions have taken place about the MSP and how the high resolution observations and ΛCDM simulations could give evidence of a plausible solution, for more details consult [68]. One of the small scale discrepancies coming from the ΛCDM simulations is the "too big to fail" problem. This states that the galaxy satellites predicted by the model are too massive that it is impossible that they doesn't have visible stars, meaning that, the observed satellites of the Milky Way are not massive enough to be consistent with predictions from ΛCDM [69]. On the other hand, the SFDM model, -which in the simplest case, considers that dark matter can be modeled by a classical, spin-zero, massive, minimally coupled to gravity, complex scalar field -reproduces, at cosmological scales, the ΛCDM, background and linear density perturbations dynamics. An important difference is that the bosonic mass, which is the only free parameter of the model, fix a natural cut-off in the mass power spectrum at small scales [70], preventing the overproduction of small structures. The mass of the scalar field needs to be larger than 10 −23 eV in order to be compatible with power spectrum measurements [71,72]. At galactic scales, such a small mass implies cored galactic halos instead cuspy as ΛCDM produces. Structure formation simulations of SFDM [73] provide density profiles for dark matter halos with a soliton core and an outer Navarro-Frenk-White profile [74]. This density profile has been used to constrain dark matter features by using different systems such as the dwarf spheroidal galaxies (dSphs) [75,76]. The analysis has shown that a single state of the SFDM tightly constrained the scalar field mass, m, with Lyman-α observations of about (log 10 (m/eV) ∼ [- 23,-24]) [77,78]. Also, an analysis of different density profiles, by using non-parametric reconstruction of rotational curves, showed that 44% of the sample galaxies preferred the SFDM model [79]. As an effort to extend this scalar field dark matter approach, it has been suggested fields with self-interactions [80,81], axionic-like potential has been widely studied from a pure gravitational approach in [82,83], see also references therein. Similarly, in [84] they used the selfinteracting term in the lagrangian and LBSG to constrain the Bose-Einstein condensate (BEC) dark matter within the range (10 −6 , 10 −4 ) eV, while in [85] dark matter is composed of ordinary QCD axions,and the fact that QCD axions form a BEC is a consequence of their properties studied. Phenomenologically, dark matter described by a classical ultralight scalar field has shown to be a good candidate of dark matter, however the approach of this matter within the quantum field theory considering the quantum nature of the possible fundamental boson associated with it, has been less explored. A hint could come from Boson Stars (BS). BS where introduced by Kaup [86] as regular localized solutions of the Einstein Klein-Gordon equations for a massive minimally coupled classical complex scalar field. Soon after, Ruffini and Bonazzola presented them as static spherically symmetric self-gravitating configurations of quantum spin-zero particles in quantum field theory using the semiclassical gravity approximation [87]. For a general overview about boson stars we recommend [88,89]. Ruffini and Bonazzola's solutions coincide with the classical configurations if all particles populate a specific energy state however, the semiclassical approximation unveil more sophisticated stars. In particular, even within the spherical symmetry and staticity, configurations in which different energy states are populated simultaneously, naturally emerge. Such kind of configurations, in Kaup's approach, require introducing multiple independent classical fields, one for each occupied state, while in the semiclassical approximation different states are excitations of a single quantum field [90]. BS with all the particles in the ground state are viable as astrophysical objects, they are stable solutions [91], and furthermore, they have been shown to be atractor solutions under fairly general initial conditions [92,93]. In the SFDM scenario they would be the kind of structures formed in isolation boundary conditions [93,94] and have been used to model SFDM galactic halos to put constrains on the bosonic mass using different observables such as galactic rotation curves (RC) [95] and stellar cluster dynamics [96]. For halos modeled by BS in the ground state, it is known that RC of different galaxies require different masses of the boson to be fitted, see for instance results for LSBG [95]. It is possible to reproduce RC of different galaxies with the same bosonic mass using configurations in excited states (all the bosons not in the ground state), however those configurations are unstable [97]. Configurations with bosons populating different states provide a diversity of mass density profiles that could account for the dark matter distribution in galaxies with diverse characteristics. In such configurations each occupied state is labeled by quantum numbers n, ℓ and m, ground states mentioned in the above paragraph have n = 0, ℓ = m = 0. In [98], configurations with n ̸ = 0, ℓ = m = 0, for which n takes different values, the so called "multi-state" configurations, were proposed to obtain flat rotation curves within the SFDM. An important condition for the physical viability of those configurations is its stability, in [99,100] it was proven that there exist stable and virialized multi-state configurations. In this work, we explore rotation curves obtained from spherically symmetric self gravitating configurations with particles in states for which the angular momentum is non-zero. The simplest case of those configurations are the so called ℓ-boson stars, in them, all particles have fixed radial and total angular momentum numbers n and ℓ, with ℓ ̸ = 0, but are homogeneously distributed with respect to their magnetic number m. In [101] they were introduced in the classical regime as configurations composed by N = 2ℓ+1 classical independent complex scalar fields with the same mass. Each of these fields have a non trivial harmonic with an angular momentum number m = −ℓ, −ℓ + 1, . . . , ℓ, but all of them have the same temporal and radial dependence. It is interesting that this special combination gives a static configuration with zero total angular momentum even when independently the fields have angular momentum and are time dependent. At this point, in the classical approximation, the existence of different fields combined with such special characteristics could be seen as somehow artificial. This is different within the semiclassical approximation, where the same ℓ-boson star is an allowed state of a single, massive, real, free quantum scalar field. The state describes the excitation of N excitation modes of the quantum field. The corresponding Einstein-Klein-Gordon system of equations, considering the expectation value of the stress energy-momentum tensor operator takes the same form as its counterpart in a classical theory with N independent complex fields [90], (see also reference [102] in which the Newtonian limit of those configurations is studied). In the same sense the "multi-ℓ-boson stars", used to model dark matter halos in this work are allowed states of the quantum field. We consider configurations with n = ℓ + 1, m = −ℓ, −ℓ + 1, . . . , ℓ, for which ℓ takes different values and we take a phenomenological point of view, observations will telling us which configurations are preferred. To find a more fundamental explanation about which modes of the quantum field should be excited is an open, interesting and necessary research field. The stability of self-gravitating scalar field configurations have been widely studied as it is an important characteristic to determine their viability as astrophysical objects. In [103,104] and [105,106] it was shown that there exist stable ℓ-boson stars with n = ℓ+1, under radial perturba-tions, in the relativistic and Newtonian limit respectively. On the other hand, ℓ-boson stars are unstable under 3D perturbations, however it is possible to stabilize them by adding a sufficiently large fundamental, n = 1, ℓ = 0 boson star [107]. Multi-ℓ-boson stars are allowed states of the quantum field that contain information of the angular momentum preserving the spherical symmetry, its 3D stability is expected if n = ℓ + 1 and if the contribution of the ℓ = 0 state is large enough, furthermore their density profile differ significantly from their relatives with ℓ = 0, in particular for ℓ > 1 the central density of ℓ-boson stars is zero while for stars with ℓ = 0 the central density is a global maximum. These characteristics make multi-ℓboson stars suitable to be tested as part of the galactic halo, in this work we only considered configurations with n = ℓ + 1 and we found that best fits of RC are possible with multi-ℓ-boson stars in which the ℓ = 0 component is important, then those multi-ℓ-boson stars modeling dark matter galactic halos should be stable although a definite proof could be done in a future work. This paper is organized as follows, in section II we present the bases of the SFDM model in the non-relativistic limit, their characteristics and the three specific cases that we analyzed: the ground state or boson star (ℓ = 0), a boson star with an excited state and multi-states. In III we describe the galaxy data set used with the methods described in section IV. Our results, presented in section V, showed that most of the galaxies have a better fit to the data by using multi-states. Finally, in section VI we discuss our results and perspectives.

II. MODEL
We take an spin 0 scalar field without self-interaction and use the non-relativistic limit of the Einstein-Klein Gordon system of equations as mentioned in [102]. For the non-relativistic ℓ-boson stars configurations the Schödinger-Poisson (SP) system is and Where G is the gravitational constant, m a is the mass of the scalar field, ℏ is the reduced Planck constant, γ nℓ0 is the frequency that will be determined by solving the eigenvalue problem, described in subsection II C, V 00 is the gravitational potential with spherical symmetry, the subscript indicates that the only term different to zero in the multipolar expansion is the monopolar one, which can be consulted in [102]; and ∇ 2 rℓ = ∂ 2 r + [2(ℓ + 1)/r]∂ r .
(2) By using the following expressions (3) we obtain a set of dimensionless equations, where the bar indicates the numerical solution, c is the speed of light and ϵ is a dimensionless quantity related with the speed of light and the rescaling amplitude for our rotational curve, being one of the free parameters of our model along with the mass of the scalar field. For simplicity, we omit the bar notation, giving us the following system of equations and To obtain the numerical solutions to the equation system, we add the following expression for the number of particles in each state [99] dN ℓ dr = ψ 2 nℓ0 r 2+2ℓ , and implemented the shooting method with a fourth order Runge Kutta. Because the system of equations only depends of the radial coordinate, we can use the expression for the circular speed (V c ) with spherical symmetry to obtain the rotational curve for the non-relativistic ℓ-boson stars where we used the change of variables for the dimensionless solutions (3) and performed the numerical integral from 0 to R to obtain the mass M (r) function In the above expression the density profile is given by where M ⊙ means solar masses. Therefore the circular speed (V c ) in terms of the new variables (3) is (9) Notice that the above equation only depends on ϵ and the radial distance has a dependency on ϵ and m a (see (3)). In particular, we choose three cases of states with zero nodes to assure stability, following the selection rule n − 1 − ℓ = 0 [102]: the ground state with ℓ = 0 (II A), an excited state with ℓ = 1 (II B) and the multi-state with ℓ = 0, 1, 2 (II C). In the next subsections we present the boundary conditions to solve the system of equations of each case and their characteristics (II D).
A. Ground state: ψ100 We solved the dimensionless equation system given by (4a) and (4b), with ℓ = 0 and n = 1, better known as a simple boson star, this solutions can be seen in the FIG 1; this case has been well studied in [97,108]. The boundary conditions to be determined by the shooting method are: γ 100 (r = 0) and V 00 (r = 0). It is important to mention that for the independent cases, i.e. the ground state II A and the excited state II B, ψ 100 (r = 0) is fixed to one, this is due to the rescaling properties discussed in subsection II D and the fact that the solutions form a family allowed us to use the parameter ϵ as the rescaling factor in the initial amplitude (ψ 100 (r = 0)). Therefore, the free parameters to be estimated for these cases are ϵ and m a .

B. Excited state: ψ210
Considering the case where the dark matter halo is only form by an excited state, which is related to the quadrupole symmetry, we solved the equation system given by (4a) and (4b), with ℓ = 1 and n = 2 different to zero. The boundary conditions are the same as the ones for the ground state II A, where γ 210 (r = 0) and V 00 (r = 0) will be determined by solving the eigenvalue problem.
C. Multi-states: ψ100, ψ210, ψ320 As it was mentioned before, we choose the states with zero nodes to assure stability and followed the selection rule n − 1 − ℓ = 0 [102]. For the multi-state case we took ℓ = 0, 1, 2. Therefore, the system of equations to be solved numerically is dV 00 dr + ψ 2 100 + 3r 2 ψ 2 210 + 5r 4 ψ 2 320 . (10d) Boundary conditions must guarantee that the solutions are regular and asymptotically flat. This implies that 10d becomes an eigenvalue problem for γ 100 , γ 210 and γ 320 . Regularity at the origin implies ψ 100 (r = 0) = C 1 , ψ 210 (r = 0) = C 2 , ψ 320 (r = 0) = C 3 and ψ ′ i00 (r = 0) = 0, for the potential V ′ 00 (r = 0) = 0. Asymptotically flatness implies ψ i00 (r → ∞) → 0 and we impose in this boundary V 00 /r + V ′ 00 = 0. To find the solutions we solved the eigenvalue problem to determine the frequencies γ 100 , γ 210 and γ 320 by using a shooting method. We took the Taylor expansions around zero to obtain the boundary conditions, finding that ψ ′′ and to avoid divergences we cancel out the terms 1/r for the integration in r = 0. The first integration starts from r = 0 to a close boundary and after finding a good frequency guess, the new frequency value is taken to the next step in the integration, until the frequencies converge to a desired precision. Therefore, ψ 100 (0), ψ 210 (0) and ψ 320 (0) became free parameters along with ϵ and m a . It is important to mention that adding the states central amplitudes to the free parameters implies the resolution of the eigenvalue problem for every step in the sampling algorithm, this makes the algorithm more computationally expensive. The expression for the number of particles in each state (5) is part of the system of equations, meaning that the solutions for N ℓ=0 , N ℓ=1 and N ℓ=2 are being found.

D. The characteristics of the solutions
Solving numerically the SP equations we could observe interesting characteristics. Taking the independent cases (ψ 100 and ψ 210 ) we observe that they form a family, meaning that, using the solution for a given central amplitude ψ nℓ0 (0) we can use the expression to obtain the solution with a different central amplitude ψ nℓ0 (0) without solving the system of equations once again, see FIG 1. This is due to rescaling properties in the equations system that can be found using (r, ψ nℓ0,V00 ) → (rλ, ψ nℓ0 /λ ℓ+2 , V 00 /λ 2 ). Therefore, for the independent solutions the parameter λ is analogous to the parameter ϵ, this equivalence is broken for the multi-state solutions, this can be seen in the equation (8), where the parameter ϵ is multiplying the sum over n and ℓ, therefore, if the parameter ϵ changes it will change all the states over the sum equally.
In FIG 2 and FIG 3, we show the numeric density profiles for the system of equations described in sections II A, II B and II C, respectively. Comparing these figures we can observe that in the FIG 2 each numerical density has a bigger amplitude and radial extension than the numerical densities from FIG 3, therefore we can ascribe the differences to the gravitational interaction between the states and the coupled system of equations that equations (10a -10d) represent. Also, we can notice the differences between an independent contribution, i.e. solving independently each state and taking the superposition of them; and a coupled contribution of each state to the total rotational curve. As these solutions are related with the spherical harmonics, therefore we can see that for the multi-states solutions, each one of them has a multipole contribution. For further details about the ℓ-boson star characteristics and within other context consult [101,103,104].

III. DATA
We use a set of Low Surface Brightness Galaxies (LSBG), these type of galaxies are disk galaxies characterized by their Tully-Fisher relationship, it is similar to the high brightness surface galaxies but with low content of stars [109], therefore it is assumed that most of their dynamics is due to the dark matter and make them good candidates to test the models.
We selected a set of 17 LSBG based on the good quality of the data classified by [110] and their mass model [111][112][113] for future comparisons. In the table I we showed the morphology of each galaxy in the chosen sample, for more details see [110]. Although their morphology could give us information about their formation history, we choose to present the results based on the radial extension and their linear behavior of the data.

IV. ANALYSIS
The Bayes theorem tell us that the probability of a model M with a parameter set θ, given an observed data set D is the posterior P where L is the likelihood, P (θ|M ) is the prior density of the parameter vector θ for a model M , containing the apriori information about the parameters of the model, and E is the evidence that will be explained in section IV C due to its importance to obtain the Bayes Factor and therefore to perform the model comparison.
As a first step in our analysis, we calculate the maximum likelihood for the independent solutions of the SP system of equations for each state (showed in FIG  2), first independently and then the superposition of each one (ψ 100 , ψ 210 and ψ 320 ), as an approximation, being the free parameters m a and ϵ for each state. This analysis gave us an insight of how each state contributes to the total rotational curve, concluding that most of the galaxies in the sample must include at least three states (the ground state and two excited). Because the results mentioned previously, we decided to center our work in the model described in section II C and used as priors in the Nested Sampler algorithm (NS) the results of the maximum likelihood for that case. As it is a coupled system of equations we cannot used the parameter ϵ of each state as a free parameter, in equation (9) we can see that ϵ became a global parameter for this case. Consequently we choose the central amplitude of each state (ψ nℓ (0)) as a free parameter, being part of the initial conditions needed to solved numerically the system of equations. This change made the parameter estimation more expensive computationally, since at each step in the NS the integration with the shooting method had to be made. The free parameters for the multi-state case are: m a , ϵ, ψ 100 (0), ψ 210 (0) and ψ 320 (0). We chose the following flat priors −26 ≤ log (m a [eV/c 2 ]) ≤ −20, −6 ≤ log (ϵ) ≤ −2, −5 ≤ log (ψ 100 (0)) ≤ 0, −6 ≤ log (ψ 210 (0)) ≤ 0 and −6 ≤ log (ψ 320 (0)) ≤ 0. For the independent cases, ground and excited states, we chose the same priors for the free parameters, m a and ϵ, respectively. For the number of live points necessary for the NS we followed the 50 × k rule, where k corresponds to the dimensionality of the free parameter vector, as a minimum. For the NS we modified the SimpleMC code that uses dynesty as an engine [114], this sampler allowed us to obtain the Bayesian evidence which was used to obtain the Bayes factor, for more details of how the NS works   [110].
consult [115]. We used the library fgivenx [116] with the output to obtain the 2σ contours in the rotational curves of FIG 9. In order to know which one of the three cases the data favours as a dark matter component, we computed the Akaike and Bayesian information criteria and the Bayes factor. The last one was compared with the ground state (ψ 100 , section II A) as the based model due to its simplicity and correspondence to the soliton profile obtained from the SFDM simulations [73], usually used in the literature with an outer NFW profile [75,76]. In addition, we take the bounds for the scalar field mass found in [113] for this case, 0.212 × 10 −23 < m a [eV/c 2 ] < 27.0 × 10 −23 . These criteria and the Bayes factor are defined in the following sections.

A. Akaike Information Criterion (AIC)
The AIC is based on information theory and it is a way to compare models for a given data [117] where k is the number of fitted parameters, n the number of data points and L is the maximum likelihood, calculated before for each model. The first term rewards the goodness of fit, while the second term penalizes the model by including extra free parameters, making it an increasing function. Therefore, the AIC discourages overfitting. By adding the last term it penalizes the fact of working with small data sets [118], which is our case.

B. Bayesian Information Criterion (BIC)
Similar to the AIC, the BIC is a model selection criterion [119] BIC = −2 ln L + k ln n, where the first term rewards the goodness of fit, and the second term penalizes the model by including the free parameters (k) and the number of data points used in the fit (n).

C. Bayes factor
Equivalent to the information criterion, the Bayes factor allows us to compare the fitness of two models, based on the Bayes theorem. The Bayes factor, B 12 , is the ratio between the posterior of a model (M 1 ) compared to another model (M 2 ), given certain data (D), in logarithmic (16) where E(D|M ) is the Bayesian evidence, defined by For a review about bayesian statistics and model selection see [120,121]. If log(B 12 ) is larger than the unity, the data slightly favours model M 1 , if the contrary occurs (log(B 12 ) is smaller than the unity), the data favours model M 2 . The table II contains the strength of evidence as the Jeffreys scale indicates [122].

V. RESULTS
Table III displays the parameter estimation obtained from the NS, the log (E) and the −2 ln L, for the cases II A and II B; the mean values are reported with 1σ confidence level. The −2 ln L and the log (E), in table III, show that the cases where the contribution of only the ground state (ψ 100 ) is present have a better fit than the cases with only the excited state (ψ 210 ) for the galaxy sample. Also, it is noticeable that for all galaxies, the mass of the scalar field (m a ) is bigger for the excited state than for the ground state. As mentioned in section III, we present our results based on the radial extension of the galaxies, they are divided in two sections: r < 10 kpc and r > 10 kpc; adding an additional restriction for the data with a linear behavior, meaning that the data points look similar to a straight line. The shaded regions correspond to the bounds for the scalar field mass found in [113], as mentioned previously. This classification is shown in FIG 4 and 5, for the independent cases, ψ 100 and ψ 210 , respectively. Here we notice that most of the galaxies with r < 10 kpc tend to have bigger masses while the galaxies with r > 10 kpc tend to prefer lighter masses. An important case raises in these contour plots, the galaxy UGC11616 has a radial extension of r = 9.6 kpc, due to the proximity to the r = 10 kpc value, it follows the r > 10 kpc behavior. On the other hand, the galaxy data with linear behavior and log  r < 10 kpc are highly correlated on the free parameters in the ground state case (ψ 100 ), this correlation seems to be broken on the excited state case (ψ 210 ). The table IV contains the parameter constraints, the log (E) and the −2 ln L for the multi-state case (II C). The mean values are reported along with with the 1σ confidence level. The contour plots related with these results are shown in FIG 6, where we have followed the classification mentioned before. The galaxies followed the same trend as the independent cases, those with r < 10 kpc tended to have bigger masses while the galaxies with r > 10 kpc tended to prefer lighter masses. For the galaxies with linear behavior the correlation seemed to be diminished between m a and ϵ. Although, for all galaxies, the correlation still remains between m a and the central amplitude of the first state (ψ 100 (0)), as we can see in FIG 7. To see the seventeen triangle plots for this case and the independent cases, see the repository 1 .
The plots in FIG 9 show the parameter estimation results reported in table IV at 1σ and 2σ as the gray color bar shows. We can observe that most galaxies have an interesting behavior regarding the contributions from each state to the total rotational curve (dark blue line), where the blue line corresponds to the ground state (ψ 100 ) shows a predominant contribution, ψ 210 (green line) contributes less than ψ 320 (orange line). Both of them contribute to the larger radius while the ground state (ψ 100 ) remains in the center, suggesting that the ψ 210 contribution is closer to zero. One can notice that galaxy UGC11583, one of the smallest with a r = 1.5 kpc radial extension has a different contribution from each state, where ψ 210 is the predominant one as ψ 320 has a smaller amplitude. Although we obtain an acceptable parameter inference for all the data sample according to the gray contours, we can observe that rotational curves for galaxies UGC11648 and UGC11748 the multi-state case does not fit the data too well. It is important to mention that unlike the NS results for the independent cases where the convergence is clear, in the multistate case, specifically for the central amplitudes of each state the convergence is not so clear and particularly the ground state central amplitude (ψ 100 (0)) seems to have a boundary that corresponds to the prior upper limit, log ψ 100 (0) = 0. The results for the NS used to obtain the plots and the plots themselves can be consulted in the repository 2 .
In table V, the AIC, BIC, log(B 12 ) and −2 ln L are reported for each case. The AIC and BIC values for the three cases are similar with a noticeable decrease for the ground state. As for the −2 ln L value, it is noticed that for all galaxies is bigger for the excited state (ψ 210 ) and the multi-state case smaller than the ground state, except for galaxy UGC11648 with a small increase. Particularly, the Bayes factor for each galaxy is represented in FIG 8, where the shaded regions correspond to the strength in Jeffrey's scale mentioned in table II. Galaxies with an asterisk ( * ) have smaller log (B 12 ) for the ψ 210 case and those (UGC11748) with a plus marker (+) have bigger log (B 12 ) for the multi-state case, that doesn't appear in the figure and can be consulted in table V. The position of the blue dots indicate that most galaxies prefer the ground state rather than the excited state ψ 210 , except for UGC11583 and ESO1200211 that seem to slightly prefer the excited state. On the other hand, the green stars position suggest that the data moderately favours the multi-state case even though the number of free parameters is bigger compare to the ground state. Some galaxies stand out, UGC11616 with a log (B 12 ) = −12.40 that supports the ground state and, galaxies UGC11648 and UGC11454 strongly supporting the multi-state case.

VI. CONCLUSIONS AND DISCUSSION
This work analyzed the viability of ℓ-boson stars as a dark matter component in the rotational curves, using bayesian statistics tools such as the NS along with the Bayes factor, together with the information criteria to know which case the data favours, being the multi-state case. Seventeen LBSG were analyzed taking into account three ℓ-boson stars cases, ground state, which was taken as the base model to calculate the Bayes factor, a single excited state and multi-states. In summary all galaxies prefer to be made up of multi-states as FIG 8 suggests. The plot's contours in FIG 9 confirmed what the Bayes factor told us about the data, by obtaining a good parameter inference with the multi-state case. Although, there are some galaxies like UGC11648 and UGC11748, that are barely fitted but clearly they don't follow the model behavior. Additionally, it is noticed that for most galaxies the main contribution to the total rotational curve (dark blue line) is by the ground state (ψ 100 , blue line) and the second excited state (ψ 210 , orange line), this could be due to the spherical system assumption, meaning that more studies in this direction need to be done by solving the multi-state axial system of equations presented in [102]. Furthermore, it is important to mention that by adding more states and therefore, increasing the free parameters, the scalar field mass tends to become slightly bigger than the ground state, this could be seen in the first columns of the tables III and IV. Although it doesn't satisfy the constrains by Lyman-α [77,78]. Although, the NS convergence for the multi-state case is not clear compared to the independent results, specially for the ground state central amplitude (ψ 100 (0)) presenting an upper limit boundary at log (ψ 100 (0)) = 0, that could be interpreted as a necessity to extend the prior for this parameter. It can always be found a way to reparameterized it.
One of the main results that we can observe in FIG 4,  FIG 5 and FIG 6, is that the correlation between the parameters m a and ϵ seemed to be broken for the excited state (ψ 210 ) and the multi-state case. However, by looking at the posteriors in the triangle plots (FIG 7 and  those in the repository 3 ), we can see that this correlation between parameters remains, with ϵ and the ground state central amplitude, ψ 100 (0). Those correlations between parameters for dark matter density profiles for rotational curves have been studied by [123], where based on the mass discrepancy acceleration relation (MDAR) and that any DM halo will have a maximum acceleration they could break the correlation by a reparametrization, having just one free parameter. Another approach that could take place in the study to break the correlation between parameters is a reparametrization with the number of particles in each state by fixing the total number of particles and adding the ratios of the number of particles in different states with respect to the ground state [99] to the boundary conditions of the shooting method. However, by trying this approach we observed that the posterior distribution functions for the ratios of the number of particles where almost flat. Therefore, it is necessary to perform more studies in this regard, taking into account the approach in [123]. Perhaps, the two approaches mentioned above and adding the self interaction term in the potential could lead to a complete new analysis and to obtain a way that constrains the number of states. It has been mentioned in the literature by [124] and references therein, that adding more states could lead to a better parameter estimation but also, it is more expensive computationally. Because of this last reason we truncate the expansion to the third term, however we do a less expensive analysis, a χ 2 fit (see A), adding a fourth term in the expansion, for ESO3050090 and UGC11616 rotational curves. We find a better fitting and an important contribution to larger radius in the rotational curve extension by the latest state. We obtained χ 2 = 0.31 and χ 2 = 10.79 for ESO3050090 and UGC11616 respectively. Which compared to the maximum likelihood obtained in the multistate case (i.e. where the expansion is truncated to the third term) is smaller. Additionally, the value of the scalar field mass is slightly bigger (see Table VI). One can notice that in this case, the states that contribute the most to the rotational curve are the ground state (ψ 100 , light blue line) and the last excited state (ψ 430 , pink line), as we can observe in the Figure 10. While a physical explanation for this is still to be found, an statistical one could be obtain by comparing the AIC and BIC. By comparing the AIC and BIC, respectively, it is observed a small increment in their values when including the ψ 430 terms (see Table VII  FIG. 4: 2D marginalized posteriors distributions of the free parameters for the ground state case (ψ 100 ). The dashed vertical line represents m a = 1.11 × 10 −23 eV/c 2 needed to have a cut-off in the power spectrum [70]. The gray band represents the bounds for the mass found in [113]. Left plot contains all the galaxies with r < 10 kpc. Middle plot: all the galaxies with r > 10 kpc. Right plot: linear behavior, the three galaxies have r < 10 kpc and they are not included in the left plot.
for the galaxy's sub-sample the multistate case is favored. Although, a more detailed analysis is needed in that direction and will be addressed in an upcoming work altogether with an improvement in the code optimization.

DATA AVAILABILITY STATEMENT
The data underlying this article are available at [110].
The datasets were derived from sources in 23 FIG. 5: 2D marginalized posteriors distributions of the free parameters for the excited state case (ψ 210 ). The vertical dashed line represents the m a = 1.11 × 10 −23 eV/c 2 needed to have a cut-off in the power spectrum [70]. The gray band represents the bounds for the mass found in [113]. Left plot contains all the galaxies with r < 10 kpc.
Middle plot: all the galaxies with r > 10 kpc. Right plot: linear behavior, the three galaxies have r < 10 kpc and they are not included in the left plot.
Appendix A: χ 2 results with ψ100,ψ210,ψ320 and ψ430 In this section we show the results obtained with the χ 2 truncating the expansion in the fourth term in the system of equation SP (4a-4b), meaning that we are adding the terms and equations corresponding to ψ 430 , the free parameters are m a , ϵ, ψ 100 (0), ψ 210 (0), ψ 320 (0) and ψ 430 (0). Where we have chosen the galaxy ESO3050090 and UGC11616 due to the value obtained for the maximum likelihood in the multistate case, described previously and their radial extension. The Table VI displays the parameter estimation obtained from the χ 2 for the galaxies ESO3050090 and UGC11616, and the   (0)), log (ψ 210 (0)), log (ψ 320 (0)) and log (ψ 430 (0)).