{'fraction_non_alphanumeric': 0.03194551309452566, 'fraction_numerical': 0.011541420877147986, 'mean_word_length': 4.547265254505216, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-state Boson Stars

Motivated by the increasing interest in models which consider scalar fields as viable dark matter candidates, we have constructed a generalization of relativistic Boson Stars (BS) composed of two coexisting states of the scalar field, the ground state and the first excited state. We have studied the dynamical evolution of these Multi-state Boson Stars (MSBS) under radial perturbations, using numerical techniques. We show that stable MSBS can be constructed, when the number of particles in the first excited state, N2, is smaller than the number of particles in the ground state, N1. On the other hand, when N2>N1, the configurations are initially unstable. However, they evolve and settle down into stable configurations. In the stabilization process, the initially ground state is excited and ends in a first excited state, whereas the initially first excited state ends in a ground state. During this process, both states emit scalar field radiation, decreasing their number of particles. This behavior shows that even though BS in the first excited state are intrinsically unstable under finite perturbations, the configuration resulting from the combination of this state with the ground state produces stable objects. Finally we show in a qualitative way, that stable MSBS could be realistic models of dark matter galactic halos, as they produce rotation curves that are flatter at large radii than the rotation curves produced by BS with only one state.


I. INTRODUCTION
The existence of Dark Matter (DM) in the Universe is strongly supported by astronomical observations that range from galactic up to cosmological scales (see for example [1] and references therein). Observations indicate that stars rotate too fast around the center of the galaxy to be bound by Newtonian gravity if all matter is visible [2,3]. This issue, known as the Rotation Curves (RC) problem, implies within the context of Einstein's General Relativity, that a great amount of the matter in the galaxy is invisible. The nature of this dark matter, that has a negligible interaction with the visible matter and whose presence is only observed trough its gravitational effects, is still unknown. The most popular candidate are the so-called weakly interacting massive particles (WIMPs) [4,5], leading to the standard Cold Dark Matter model. This scenario is very successful at a cosmological level, as its predictions are in good agreement with the observational data [6,7]. However, it has difficulties in fitting the observations at a galactic level [8][9][10][11][12]. If DM is modeled by WIMPs, one obtains a cuspy density profile of the DM in the galaxy. But high resolution data of low surface brightness galaxies, which are composed mainly of DM, imply that their DM distribution has a flat core [10,11]. This model also fails in predicting the number of satellite galaxies around each galactic halo, exceeding far beyond what is observed around the Milky Way [12].
A different approach consists in describing the dark matter as a scalar field [13][14][15]. The Scalar Field Dark Matter (SFDM) model has been proved to be successful at cosmological scales [16]. This model can also avoid the problems that WIMPs present at a galactic level, producing a non-cuspy density profile [16][17][18] and explaining the dearth of satellite galaxies around each galactic halo [16]. Because of the viability presented by the SFDM model, it is stimulating to go further on testing it. For instance, the model has to reproduce the observed RC of galaxies. At this point Boson-Star-like objects could play an important role.
In the SFDM model, the dark matter particle is an ultra-light massive spinless boson (m ∼ 10 −23 eV [16]). These bosons could collapse forming gravitationally bounded structures. With such ultra-light mass, the boson's Compton wave length is of the order of kiloparsecs, and structures with comparable length scaleslike galactic halos -could be formed as condensates described by a coherent scalar field [19]. These condensates can be associated with Boson stars (BS), which are solutions of the Einstein-Klein-Gordon equations where the gravity attraction is balanced by the dispersive character of the scalar field [20,21]. BS were first studied by Kaup [22] and a year later by Ruffini and Bonazzola [23], who settled two different treatments. The first one, developed in [22], is a completely classical treatment with a massive complex scalar field minimally coupled to gravity. In the second one [23], a real quantized scalar field is introduced in order to describe a many boson system though maintaining the geometry as a classical entity (i.e., a semiclassical limit is adopted). The relevant quantity computed in this case is the mean value Q|T µν |Q of the energy momentum tensor operator, where |Q is the state of the system of many particles. When one considers a |Q for which all the particles are in the same state, it turns out that the mean value ofT µν generates the same energy-momentum tensor as the complex classical field and consequently, the same macroscopic results, i.e. both a quantized real scalar field and a classical complex scalar field yield to the same self-gravitating system.
Until very recently, only boson stars with all the particles in one state have been considered. However, Newtonian configurations with scalar fields coexisting in the ground and excited states were introduced in [24] in order to model dark matter halos. Previous studies in modeling dark matter halos using newtonian BS were done in [25][26][27][28][29]. However, these structures can not account for a realistic halo since the configurations in the ground state produce RC which are not flat enough at large radii. The case of a massless scalar field used in order to fit rotation curve data of several galaxies was considered in [30]. Nevertheless, it has been proved in [31,32] that no nonsingular self-gravitating solitonic objects can be formed with a massless scalar field. On the other hand, RC from excited BS are in better agreement with the astrophysical observations, but these structures are unstable [33].
In [34] it was shown that Newtonian mixed configurations could account for more realistic DM halos, as they are stable and could fit with better agreement the observed RCs even at large radii. In the present work, we are interested in the fully relativistic generalization of these Multi State Boson Stars (MSBS). The idea is to consider the possibility that the bosons are not all in the same state, but rather populating different coexisting states, as was already pointed out in [23]. It turns out (see [23] and Appendix A) that the resulting equations for a MSBS in the semiclassical approach are equivalent to the case where a collection of complex classical scalar fields is considered, one for each state, which are only coupled through gravity. Without any loss of generality, we can choose either the semiclassical formulation of such MSBS or its pure classical counterpart. We will follow the latter approach in order to investigate MSBS, analyzing in detail the properties and stability of configurations with two states, a ground state and a first excited state.
The previous stability studies of BS can be divided roughly in two categories, depending on the type of perturbations considered: 1. Studies where the perturbations preserve the number of particles (infinitesimal perturbations), which generally involve linear perturbation analysis [35,36] and catastrophe theory [37].
2. Those where the perturbations don't conserve the number of particles (finite perturbations), which have been addressed mainly by numerical studies. [33,38,39]. Furthermore, the late time evolution of unstable boson star under finite perturbation can only be followed by numerical simulations.
A consistent result coming from both type of studies is that BS in the ground state are stable against perturbations if the amplitude of the scalar field at the origin φ(0) is smaller than the critical value φ max (0) where the maximum mass M max is reached. In the case of excited BS there are some important differences. Although the linear stability analysis shows stability up to the critical value φ max (0) [21,36] when the number of particles is not conserved, excited BS are intrinsically unstable even for φ(0) ≤ φ max (0), since finite perturbations drive the star either to collapse to a black hole or to decay to the ground state [33,36]. From these results one could infer that the MSBS states would be unstable under perturbations when the number of particles is allowed to change, since they contain at least one excited state. Quite surprisingly, our numerical analysis shows that there is a region of the solution space with stable configurations. Roughly speaking, the ground state produces a deeper gravitational potential which can be enough to stabilize the excited state.
This paper is organized as follows. In Section II, we present the formalism used for the numerical evolution of the Einstein-Klein-Gordon system, describing a MSBS in the classical approximation. Section III describes how the initial data for a MSBS with two different states is constructed. In Section IV A, we present numerical results obtained from the evolution of the two-state boson stars. We study two features of the evolution, namely the stability and the late time behavior. In both cases, we add an small perturbation. In the first case, the perturbed MSBS is evolved for short time scales, in order to study the behavior of the perturbations. In the second case, we evolve unstable MSBS for longer times following the properties of the resulting configurations. We compute in section V RC from stable MSBS and discuss qualitatively why these RC are in better agreement with the observed RC of galaxies. We conclude in Section VI. The description of the semiclassical approach is presented in Appendix A, while Appendix B is devoted to a detailed description of the evolution equations used for numerical evolution.

II. THE EINSTEIN-KLEIN-GORDON SYSTEM
Let us consider a semi-classical real massive scalar field with P different excited states, which is equivalent to considering a collection of P classical complex scalar fields (one for each state) coupled only through gravity. In a curved spacetime, the dynamics of these MSBS can be described by the following Lagrangian density (adopting geometrical units, i.e. G = c = = 1), where R is the Ricci scalar, g ab is the spacetime metric, φ (n) are the scalar fields,φ (n) their complex con-jugate, and V (|φ (n) | 2 ) a potential depending only on |φ (n) | 2 . Throughout this paper, Roman letters from the beginning of the alphabet a, b, c, .. denote spacetime indices ranging from 0 to 3, while letters near the middle i, j, k, .. range from 1 to 3, denoting spatial indices. This Lagrangian gives rise to the equations determining the evolution of the metric (Einstein equations) and those governing the scalar fields behavior (Klein-Gordon equations). The variation of the action associated with the Lagrangian (1) with respect to the metric g ab , leads to the well-known Einstein equations where R ab is the Ricci tensor. T ab is the total stressenergy tensor, given by the addition of the single stressenergy tensors of each scalar field, namely The Einstein equations form a system of 10 non-linear partial differential equations for the spacetime metric components g ab . On the other hand, the variation of the Lagrangian (1) with respect to each scalar field φ (n) , leads to a set of Klein-Gordon (KG) equations which are only coupled through the gravity, where the box = g ab ∇ a ∇ b stands for the wave operator on a curved background. In the following, we will restrict ourselves to the free field case, where the potential takes the form with m a parameter that can be identified with the bare mass of the field theory. The matter Lagrangian is invariant under global U(1) transformations This symmetry implies that there is a set of Noether currents densities J satisfying for each n the conservation law ∇ a J (n) a = 0. The Noether charge contained in some radius is given by so that the total Noether charge of the system N is the sum of the total individual ones N (n) ≡ N (n) (∞), namely As discussed in [23], this quantity N can be associated with the total number of bosonic particles. Consequently, N (n) can be interpreted as the number of particles in the state labeled by n.

III. INITIAL DATA FOR MULTI-STATE BOSON STARS
The initial data for the MSBS configurations is computed in spherical symmetry with a one-dimensional code. We adopt the following harmonic ansatz for each scalar field, With this assumption, the source for the Einstein equations becomes time independent. Our goal is to find {φ n (r), ω n } and the metric coefficients, such that the spacetime generated by this matter configuration is static. We begin by considering the problem in polar-areal coordinates [22,40]. The line element in these coordinates takes the form Then the equilibrium equations, obtained by substituting the ansatz (11) and the metric (12) in the Einstein-Klein-Gordon system (2,5), are given by In order to obtain a solution of this system, we provide the following boundary conditions, motivated by the physical situation under study, which guarantee regularity at the origin and asymptotic flatness. For given central values of the fields {φ cn }, we only need to adjust the eigenvalues {ω n } and the value α(0) in order to generate a solution with the appropriate asymptotic behavior (20)(21). This is a shooting problem that we solve by integrating from r = 0 towards the outer boundary r = r out , with a second order shooting method. The boundary conditions for the scalar fields at r out are imposed considering that localized solutions decrease asymptotically as φ n ∼ exp − m 2 − ω 2 n r /r in a Schwarzschild-type asymptotic background. At the outer boundary, the conditions are The shooting procedure is performed for different values of r out . As r out is increasing, the shooting parameters converge, and we choose the solution as the one which satisfies the conditions (22) for some r out within a prescribed tolerance. From this point on, we match to the scalar fields and the metric coefficients their asymptotic behavior.
A qualitative characteristic of the radial functions φ n is their number of nodes (ie, how many times they do cross zero), which reflects the excited state of the boson star. If the radial function does not have any node, the boson star is in the ground state. When there is a node, the boson star is in the first excited state, and so on. In the next subsection we construct initial configurations with two scalar fields P = 2, one in the ground state and the other in the first excited state. Notice that this is the simplest non-trivial configuration, since the MSBS with two scalar fields in the ground state can be reduced to one scalar field solution by redefining the scalar fields. This is a consequence of the indistinguishably of the boson particles in the same state.
Once the solution is computed in this coordinate system, a change of coordinates is performed to maximal isotropic ones, which are more convenient for our numerical evolutions. Finally, a simple inspection on the system (13)(14)(15)(16) shows that the re-definitioñ leads to a set of equations that no longer have m on it. Hence, the selection of geometrical units (e.g. G = = c = 1) and the redefinition (24) give us dimensionless units forr,ω n , which is equivalent to choosing m = 1 in our equations. Throughout this paper, we use these dimensionaless coordinates. Let us consider the simplest non-trivial case with only two scalar fields P = 2, one with N (1) particles in the ground state, and the other with N (2) particles in the first excited state. A useful way to construct the initial data is specifying the fraction between the number of particles in each state of the configuration, In this case, we complete the system (13-16) with the differential expressions for the number of particles in each state with boundary conditions, N (n) (0) = 0. If η is specified, it is sufficient to prescribe as boundary conditions the central value of only one of the scalar fields, for instance, φ c1 . The new system of equations (22) and (26) becomes a shooting problem for the four parameters {ω 1 , ω 2 , α(0), φ 2 (0)}. For a specific fraction η, it is necessary to adjust the four parameters such that equation (22) and the condition N (2) (r out ) = ηN (1) (r out ) are satisfied. Fig. 1 shows an example of the radial profiles of the two scalar fields, lapse and conformal factor, for a MSBS with {η = 1, φ 1 (0) = 0.0197}.
Two important characteristics of MSBS are the total gravitational mass M and the radius R 99 . The first one is calculated as where m p the Planck's mass and m the mass of the boson associated to the scalar field. In Fig. 2 we have plotted the values of these two quantities for all the constructed MSBS initial configurations. On the top panel they are shown as functions of the central value of the scalar field in the ground state φ 1 (0) and the Noether fraction η, while in the bottom panel they are plotted as functions of the central value of the scalars fields.
These figures already show some differences between single BS and MSBS configurations: there are an infinite number of possible equilibrium configurations (i.e. solutions for the static EKG system (13)(14)(15)(16)) between the two extreme cases N (1) = 0 and N (2) = 0 that correspond precisely to the ground state BS and first excited state BS respectively. Fig. 2 shows a more complex behavior of the mass M and R 99 than for a single BS, which can be observed easily in Fig. 3 where some slices φ 1 (0) = constant and φ 2 (0) = constant of these surfaces are displayed. After constructing these configurations, the next problem is stability. It is clear that the known results of BS stability are not immediately applicable to MSBS.
The stability of the multi-states configurations under finite perturbations will be addressed in the next section.

IV. NUMERICAL SIMULATIONS
In this section, we present a numerical analysis of the dynamical properties of MSBS, focusing on long-term stability and the final state of the unstable configurations. To this purpose, we write the Einstein-Klein-Gordon system as a set of evolution equations for the scalar fields and the metric components. We consider a generic spherically symmetric spacetime with the line element where α is the lapse function and {g rr , g θθ } are the metric components. Notice the explicit dependence on the factor r 2 , such that the component g θθ is regular at the origin. This is a necessary condition in our implementation in order to deal with the coordinate singularity at r = 0. The evolution equations for the geometry are obtained by substituting the metric coefficients (29) into a particular formulation of the Einstein equations. In this study we have considered the Z3 formulation [41], which in-cludes the momentum constraint into the evolution system, by considering an additional vector Z i as an evolved field. Further details regarding the Z3 system in spherically symmetry can be found in [42], while the regularization of the coordinate singularity r = 0 is similar to the one described in [43].
In spherical symmetry, there are independent evolution equations only for the lapse α, the metric components {g rr , g θθ }, the extrinsic curvature {K r r , K θ θ } and the Z-vector components {Z r , Z θ }. All these evolution equations are prescribed by the Einstein equations, except the one corresponding to the lapse, which is related to the choice of coordinates and can be specified freely. A common choice, which give rise to a hyperbolic system of equations, is the harmonic slicing where trK = K r r + 2 K θ θ . On the other hand, the evolution equations for the scalar fields are obtained by substituting the spherically symmetric metric (29) in the Klein-Gordon equations (5).
A first order reduction in space can be performed by introducing as independent quantities the spatial derivatives of the metric and new fields related to the time and spatial derivatives of the scalar field, namely In this way, we obtain a fully first order system of evolution equations for the geometry and the scalar fields, with the following set of evolution variables {α, g rr , g θθ , K r r , K θ θ , A r , D rr r , D rθ θ , Z r , φ (n) , φ (n) r , φ (n) t }. This first order system can be written in balance law form which allows the use of advanced numerical methods based on Finite Volume algorithms. Details about the exact form of the evolution equations can be found in Appendix B.
We have implemented the equations using the Method of Lines, in order to separate the time and the spatial discretization. The time integration is performed with a third order Strong Stability Preserving Runge Kutta method [44]. The spatial discretization is based on a standard fourth order centered finite difference scheme, plus third order accurate dissipation [42].

A. Stability of the MSBS
The stability of MSBS configurations is a basic requirement for considering them suitable models of galaxy halos. For a single boson star, the stability has been previously studied both analytically and numerically, showing that it is stable if φ(0) ≤ φ max (0). In this section, we analyze the stability of MSBS in the range φ 1 (0) ≤ φ max (0), which are stable for η = 0, avoiding this way the toomassive unstable MSBS. The stability will be studied following only a numerical approach, by perturbing the MSBS and studying the evolution of this perturbation. We will restrict ourselves to study stability against spherically symmetric perturbations by using the equations described in the previous subsection. Notice that this is only a necessary condition for the most general case, since asymmetric perturbations may still be unstable.
In order to study numerically the stability of the MSBS configurations, we perform the following steps: • Construct different initial data sets for MSBS with a given φ 1 (0), by varying the Noether fraction η.
• Add a real scalar field far outside the radius R 99 of the MSBS, which will be coupled to the MSBS only through gravity. The energy density corresponding to the scalar field is only 0.01% of the total energy density of the MSBS, so it will act just as a small perturbation with negligible errors in the form of constraint violations.
• Perform evolutions of the Einstein-Klein-Gordon system and study the behavior of the MSBS. The scalar field perturbation will fall into the MSBS and later disperse to infinity. The gravitational interaction during that time is expected to excite the unstable modes, if any. In this way, the modes are excited sooner than only by numerical errors.
• Bracket the MSBS which lead to significant exponential growing modes. For the stable MSBS configurations, the perturbations will only oscillate without growing. We expect that the MSBS with low η will be stable, since the major contribution to the complete configuration comes from the stable ground state, while those with high η will correspond to unstable MSBS.
• Fit the growth rate of the unstable MSBS for each set of stars with the same φ 1 (0) by varying η. Extrapolate to find the maximum allowed Noether fraction η max which separates the stable and unstable states. This procedure allows us to obtain reliable estimations, without evolving every configuration in order to obtain η max . Moreover, it might be difficult to distinguish stable from unstable configurations when they are close to η max , since the exponential growth is very low in that region.
We are going to restrict the numerical stability analysis to only three different values of φ 1 (0) = {0.0143, 0.0197, 0.0423}. Fig. 4 displays the total mass and the radius R 99 for these configurations, as a function of η. In the simulations with η ≤ 1.2, we did not detect any unstable exponentially growing mode, or they were difficult to measure for some families of solutions. The results indicate as the upper bound η max < 1.2. In fig. 5, we show the typical behavior for one of these simulations, corresponding to stable and unstable MSBS. The perturbation has an exponentially growing behavior only for the unstable MSBS. The maximum of the scalar field φ 2 (r = 0) can be fitted with a function of the type A exp(σt) cos(ωt + ϕ), which allows us to compute the exponential growth rate σ. We performed fits for the unstable MSBS perturbations with η ≥ 1.2, marked with filled geometrical shapes in fig. 4. The results for σ are represented in fig. 7, with the extrapolation to the η max , which in principle could be a function of φ 1 (0). The three different families of solutions point to η max ≈ 1.
In order to show the robustness of these results we have repeated the simulations with more resolution and with a different amplitude of the perturbation. The results are almost identical, as shown in fig.6. The only significant difference is the unstable case with larger perturbation amplitude. The perturbation seem to interact non-linearly with the star, and the unstable exponentially growing mode is excited later. In spite of this delay, the growth rate is identical to the other cases with smaller perturbation and higher resolution.

B. Fate of the unstable states
Another question which arises from the previous stability analysis, refers to the final fate of the unstable MSBS with η > η max . We address this issue by performing long evolutions of unstable MSBS configurations, until they reach a stationary state. In order to accelerate the growth of the unstable modes, we perturbed the MSBS with a massless scalar field located far from R 99 . As explained in the previous section, this scalar field interacts gravitationally with the MSBS, perturbing it slightly and exciting the unstable modes. These modes grow exponentially, starting with a small amplitude, result that can be obtained also from a linear perturbation analysis. When the amplitude of these perturbations is larger, the nonlinear effects become important and the evolution can only be followed numerically in order to discern the final state of the MSBS. Fig. 8 displays the number of particles in different states and the total number of particles for two MSBS with the same φ 1 (0) = 0.0197 with a fraction given by η = 3 and η = 0.4. The total Noether charge remains almost constant in the stable case η = 0.4, showing the accuracy of the numerical code within a 0.004% error in this quantity. The unstable case η = 3 exhibits scalar field radiation during the evolution, producing a decrease of around 18% in the total Noether charge, as it can be seen in the convergence test presented in fig. 9. This radiation translates into a small change in the amplitude of the scalar fields. Taking a closer look at the maximum value of the scalar fields in the center and the frequen-  cies, displayed in figure 10, one can notice a change in the position of the node; the excited state has decayed to a ground one, while the ground one has jumped to the first excited state. With this "flip-flop" of the scalar fields, the final η is in the stable domain. A similar behavior is observed for all the unstable MSBS configurations included in the study, which indicates that this could be a The previous statements can be seen graphically in Fig.  11. The figure shows the collection of initial configurations, each one characterized by the number of particles in the ground state N 1 , and the number of particles in the excited state N 2 . Every pair (N 1 , N 2 ) corresponds to a configuration with a pair of eigen-values (ω 1 , ω 2 ), although only ω 1 is shown in the figure. On the top of these collection of initial configurations, we have plotted the time evolution of the previous simulations. The red dot corresponds to the configuration labeled with fraction η = 0.4 and the blue dots represent the time evolution for the configuration labeled with fraction η = 3.0. Applying the same perturbation for the two different states, we can discern two different behaviors: the fraction η = 0.4 is stable and remains static, while the configuration with η = 3.0 is unstable and evolves far from the original configuration (that lies on the sheet of static configurations) to a "forbidden" region. The jump is due to the flip of ground state → excited state and viceversa. Then, it starts to lose scalar field (i.e., evolves in (N 1 , N 2 )) and slowly approaches an equilibrium configuration on the sheet of equilibrium configurations. The final state of the MSBS with η = 3 can be inferred from the central value of the amplitude of the scalar fields (top panel Fig. 10). It is important to remind that now, due to the switch of the frequencies, the new ground state corresponds to φ 2 (0) while the new excited state corresponds to φ 1 (0). The amplitudes oscillate around a central value, namely Using these values as initial conditions for the new equilibrium configuration we can compute the eigenvalues and the number of particles shown in Table I. In the  same table are displayed the values obtained from the evolution of the MSBS configuration with η = 3. There is a remarkable agreement in the eigenvalues, but they still differ in the number of particles. This means that  the system will still loose particles at a slow rate during the evolution, as corroborated by Fig. 11. The slow loss of particles has been observed in finite perturbed system of excited states and even tough the final state has been inferred in the same way as we have done here [33,38]. In previous sections we have constructed MSBS with two states, a ground and a first excited state, and shown their stability. In this section, we come back to our initial motivation and we will illustrate how MSBS can lead to RCs which are in better agreement the RCs of galaxies in the context of SFDM. A detailed analysis of fitting the RC of galaxies including baryonic matter and experimental data is out of the scope of the present work. Instead, we will just present a comparison between the behavior of a test particle immersed in the gravitational potential produced by a MSBS and a single BS, showing that there is an improvement in the sense that the MSBS has a flatter profile far from the center.Neglecting the baryonic contribution is a reasonable assumption in some galaxies such as the Low Surface Brightness Galaxies.
For the static spherically symmetric metric considered here (29), the circular orbit geodesic obeys [45], As an illustrative example, the top panel of Fig. 12 presents a comparison between two rotational curves obtained for a single state boson star with φ 1 (0) = 0.0197 (ie, a MSBS with η = 0) and for a MSBS with the same amplitude φ 1 (0) = 0.0197 of the ground state scalar field, and the same number of particles in the first excited state (ie, η = 1). The region with a flat plateau is larger, suggesting that for higher excited states (or MSBS with several higher states) the region with constant rotational velocities could be extended to larger radii.
In order to understand this behavior, it will be helpful to see the mass density profile defined as and it is shown in the bottom panel of Fig. 12 for the same configurations mentioned above. We can see that for the single BS (η = 0) the density decays exponentially as r → ∞, making it difficult to fit the flat rotational curve profiles present in most galaxies. However, the MSBS configurations with large Noether fractions have a radius which is significantly larger than the one corresponding to the single ground state, with an exponential decay only in the tail of the excited state.
Another issue related with boson stars as dark matter models, was the lack of degrees of freedom to match the different sizes and masses of the observed galaxies. For a single boson star without self-interaction, the only free parameters are the mass of the boson particle m and the central value of the scalar field φ(r = 0), which determines the compactness of the object (ie, ratio of total mass over radius) in adimensional units. There have been several attempts to fit these parameters [24][25][26][27][28] with different levels of success. By allowing more general MSBS, there are extra free parameters coming from the different fractions between the ground and excited states. These parameters change not only the total mass, but also the compactness of the final object. The extra degrees of freedom may allow a better fit of the models to different galaxies.

VI. CONCLUSIONS
We have constructed generalized boson star configurations, where two coexisting states of the scalar field are present. Our initial data construction is based on two main quantities that describe them: the gravitational mass M and the radius R 99 . We have shown that these boson stars are stable under small radial perturbations, for a certain range of the fraction (η < 1) between the Noether charges. These results may sound counterintuitive, given the known fact that single BS in excited states are unstable under finite perturbations. Nevertheless, the addition of an extra scalar field allows for an infinite number of new equilibrium configurations. The known plot of M vs. φ 1 (0) for ground state or excited state single BSs is now extended, getting a different curve for each fixed value of φ 2 (0) = constant, as it is shown in the figure 13. The single BS in the ground state corresponds to the case φ 2 (0) = 0, which only has an extreme at the maximum allowed mass. The configurations on the right of that point (marked with a triangle) are unstable. For low values of φ 2 (0) these curves contain now two extremes. In addition to the maximum allowed mass, there is a minimum close to the fraction η ≈ 0.5. Although the presence of a new extreme in the curves could suggest a change in the stability around the fraction η ≈ 0.5, our numerical stability analysis presented in the section IV A indicates that the dynamical stability regime is extended beyond this point up to η max ≈ 1. Furthermore, we have found also stable and unstable configurations in regions with no extremes in the figure 13. This supports the idea that simple stability arguments from the single boson star case are not valid in this case and that a more detailed future study is necessary in order to confirm and understand completely this problem.
The unstable configurations evolve and settle down into stable configurations. MSBS allow to obtain a flat region in the velocity rotational curves, as shown by the examples in this paper. We considered cases with only two different states of the scalar field. As future work, we are planning to construct MSBS where several states are coexisting. These models allow more degrees of freedom, and could be used to fit accurately the rotational curves within the observational data.