Dynamical $\ell$-boson stars: generic stability and evidence for non-spherical solutions

$\ell$-boson stars are static, spherical, multi-field self-gravitating solitons. They are asymptotically flat, finite energy solutions of Einstein's gravity minimally coupled to an odd number of massive, complex scalar fields. A previous study assessed the stability of $\ell$-boson stars under spherical perturbations, finding that there are both stable and unstable branches of solutions, as for single-field boson stars ($\ell=0$). In this work we probe the stability of $\ell$-boson stars against non-spherical perturbations by performing numerical evolutions of the Einstein-Klein-Gordon system, with a 3D code. For the timescales explored, the $\ell$-boson stars belonging to the spherical stable branch do not exhibit measurable growing modes. We find, however, evidence of zero modes; that is, non-spherical perturbations that neither grow nor decay. This suggests the branching off towards a larger family of equilibrium solutions: we conjecture that $\ell$-boson stars are the enhanced isometry point of a larger family of static (and possibly stationary), non-spherical multi-field self-gravitating solitons.


I. INTRODUCTION
Boson stars [1,2] (see [3,4] for reviews) are remarkable gravitational solitons. These self-gravitating, localised energy lumps of a complex, massive scalar field have appealing theoretical properties. A key one is their dynamical stability. For spherical boson stars there is a stable branch of solutions against perturbations. Indeed, a variety of studies including linear perturbation theory [5][6][7], catastrophe theory [8] and numerical simulations [9][10][11][12][13] agree that boson stars are perturbatively stable, as long as the amplitude of the scalar field is smaller than a critical value. When the latter is attained, boson stars acquire their maximum mass.
Being dynamically stable legitimates inquiring about the possible (astro)physical role of boson stars. Albeit exotic, lacking undisputed observational evidence, boson stars have found important applications in strong gravity and astrophysics. For instance, boson stars provide a common model for a black hole mimicker [14][15][16][17]. Being dynamically tractable, one can then compare dynamical spacetime properties, such as waveforms of binary boson star systems, with those of black holes [18][19][20]. This is particularly timely in view of the recently initiated gravitational-wave astronomy era [21,22], which provides data for both models to be compared with.
A second important application is in relation to a central mystery of contemporary science: the nature of dark matter. An increasing attention has been dedicated to models that consider dark matter as an ultra-light bosonic particle [23][24][25][26]. The bosonic nature allows this sort of dark matter to form co-herent macroscopic excitations. In this context, bosons stars can model, in particular, the core of dark matter galactic halos [27][28][29][30].
In their original guise, the Einstein-Klein-Gordon (EKG) model contains a massive, free scalar field, and the solitonic solutions are called mini-boson stars. A variety of generalisations ensued. Boson stars for scalar field theories with self-interactions have been reported, starting with the case of quartic self-interactions considered by Colpi et al. [31]. Spacetime angular momentum was introduced for mini-boson stars in [32,33], giving rise to stationary (but not static) selfgravitating solitons. A cousin model with a complex, massive vector (rather than scalar) field yields Proca stars [34]. These and other examples use single (complex) field models; however, multi-field boson stars have also been reported. One example is given by multi-state boson stars [28,29,35]. More recently, multi-field boson stars with an arbitrary odd number, N = 2 + 1, ∈ N 0 , of equal mass, uncoupled (except through gravity) complex scalar fields with harmonic time dependence were introduced [36]; they are dubbed -boson stars and they will be the focus of this paper.
-boson stars are described by spherically symmetric and static metrics. For = 0 they are simply the usual mini-boson stars. For 1, the 2 +1 scalar fields have an angular dependence given by the corresponding 2 + 1 spherical harmonics Y m . Then, if the radial dependence for all fields is the same, corresponding to choosing the amplitude of the spherical harmonics equal at all radial distances, static, spherical configurations are obtained, regardless of the energy-momentum tensor of each individual field being angular dependent. This is an example of symmetry non-inheritance: the (spherical) spacetime and the (non-spherical) individual matter fields do not share spherical symmetry. The usual boson stars already have a version of symmetry non-inheritance: the (time oscillating) scalar field and the (static) spacetime do not share time-translation symmetry. Consistency requires only that the spacetime geometry and the total energy-momentum tensor share the same symmetries, not the individual matter fields.
Generic -boson stars have been shown to exhibit similar properties to those of the standard = 0 stars. In particular,boson stars have a stable branch of solutions against spherical perturbations [37]. The main goal of this paper is to assess the stability of -boson stars (in this branch) against generic, non-spherical perturbations. As we shall see, our analysis will show that, in this respect, generic -boson star do not exactly mimic the = 0 case. Although no instabilities are observed, the analysis provides a glimpse of a larger landscape of solutions, of which -boson stars are just the enhanced symmetry point.
Departure from spherical symmetry is physically relevant. Firstly, spherical objects -such as -boson stars -need to be stable against non-spherical perturbations, in order to be dynamically viable. Secondly, astrophysical bodies are not, typically, perfectly spherical, in particular due to angular momentum. So one must assess if some perturbations actually deform -boson stars into acquiring new degrees of freedom. In this respect, it was recently proposed that multifield boson stars, in the non-relativistic regime, could have non-spherical stable configurations [38]. This provides an extra motivation to inquire about the behaviour of relativistic -boson stars under non-spherical perturbations Finally, assessing non-spherical configurations and perturbations often yields a richer phenomenology. As a fruitful example, it was recently found that spinning, single-field mini-boson stars are unstable against non-axisymmetric perturbations, either decaying into a non-rotating boson star or collapsing into a Kerr black hole [20,39]. By contrast, spinning Proca stars do not present instabilities under non-axisymmetric perturbations and furthermore, they can form dynamically [39]. This example shows how the study on non-spherical perturbations unveiled a new relevant dynamical property of boson stars.
We shall investigate the behaviour of -boson stars under non-spherical perturbations using fully non-linear numerical simulations of the corresponding Einstein-Klein-Gordon system. As initial data, we use configurations found in [36] which are then perturbed in two different ways. The first type of perturbation tests the stability against non-axially symmetric perturbations, targeting potential bar-mode instabilities. The second type of perturbation tests the stability against a relative change in the amplitude of the internal fields. In none of the two cases measurable growing modes were found, either by perturbing the total mass density or by perturbing each of the constituent fields, as long the -boson star belongs to the stable branch against spherical perturbations. By following the evolution of distortion parameters (defined below) we found, however, evidence for long-lived perturbations, which we interpret as zero modes. These modes, in turn, are interpreted as evidence for a larger family of equilibrium solutions.
Consider the Schwarzschild black hole of vacuum Gen-eral Relativity. It has been shown to be mode stable in the renowned works of Regge and Wheeler [40] and Zerilli [41].
No gravitational perturbations grow. However, a perturbation that carries angular momentum yield not decay. The Schwarzschild solution migrates to a small angular momentum Kerr solution and oscillates around this new ground state. Similarly, a perturbation which electric charge will not decay and the spacetime will oscillate around a small charge Reissner-Nordström solution. These special perturbations are zero modes. Such modes are often found when a spacetime is unstable against some sort of perturbations, at the threshold between stable and unstable modes. An example occurs for the superradiant instability of the Kerr spacetime due to a massive bosonic field. The zero modes indicate the bifurcation of the Kerr family towards a new family of black holes with bosonic hair [42,43]. But zero modes can also occur even if there is no instability, as in the Schwarzschild example, indicating, nonetheless, an enlarged family of solutions (Kerr or Reissner-Nordström), of which the initial spacetime (Schwarzschild) is a special case. Thus, one of the outcomes of our analysis is the conjecture that -boson stars are the enhanced isometry point of a larger family of static (and possibly stationary), non-spherical multi-field self-gravitating solitons.
In the rest of this work we will focus on configurations with = 1. Such -boson stars are described by N = 3 fields, with m = −1, 0, 1 respectively. In order to follow the dynamics of the perturbed system a numerical code that solves the Einstein-N -Klein-Gordon system is required. We have used the EINSTEIN TOOLKIT framework [44][45][46] with the CARPET package [47,48] for mesh-refinement capabilities to achieve our goal.
As a technical step we perform a Cauchy (3+1) decomposition on each scalar field that constitutes the star and solve the full Einstein-N -Klein-Gordon system. This is done implementing an arrangement in the EINSTEIN TOOLKIT, a thorn, to solve N scalar fields using finite differences [39]. This paper is organized as follows: Section II addresses the construction of initial data to set up perturbed -boson stars. Section III describes the diagnostic tools used to monitor the evolution and some aspects used to decide on whether instabilities are present. The numerical results are described in Section IV and in Section V our conclusions and final remarks are presented. In this work we use units where G = 1 = c.

II. INITIAL DATA
Following previous works on -boson stars [36,37], we consider a set of N = 2 + 1 complex scalar fields, with mass µ and no self-interaction within the Einstein theory of gravity, for which the energy-momentum tensor is given by: where the index i labels each field and the stress-energymomentum for each field is given by Complex conjugation is denoted by '*'. Following [36,49] we propose a set of scalar fields of the form where the angular momentum number is fixed, and m, which plays the role of index i in equation (2.3), takes the values m = − , − + 1, . . . , (hence the total number of fields needed for a fixed value of will be 2 +1), Y m are the spherical harmonics defined over the unitary 2D-sphere. Then we assume that the amplitudes ψ (r, t) are the same for all m. It was shown in [36] that if the N fields have all the same amplitude ψ , the stress-energy tensor (2.1) has spherical symmetry regardless if the fields have angular dependence. See also [49] and for a detailed discussion on the procedure, see [50]. Assuming the harmonic time dependence where φ (r) and the frequency ω are both real-valued, the stress-energy tensor becomes time independent. Under these assumptions it is possible to find self gravitating static, spherically symmetric equilibrium configurations by solving the EKG system of equations. Those configurations are parametrized by the angular momentum number , hence the name, -boson stars.
In order to obtain initial data suitable for numerical evolution, we construct equilibrium -boson stars, to be subsequently perturbed. Considering a spherically symmetric spacetime with a line element given by: where α and A, are functions of r, and the assumptions mentioned above for ψ , the EKG system yields (2.8)  Table I.
By studying the Klein-Gordon equation in the vicinity of r = 0 one finds that the scalar field behaves as φ ∼ φ 0 r in that region. For a fixed value of the angular momentum number , a given value of the parameter φ 0 , and the boundary condition at infinity requesting that φ decays exponentially, the system of equations (2.6-2.8) becomes a nonlinear eigenvalue problem for the frequency ω. We solve this set of equations in a finite size grid by means of a shooting method using the frequency ω as the shooting parameter. For numerical purposes we take the mass parameter µ = 1. Ref. [36] it was shown that -boson stars with > 0 have similar properties to those of single-field mini-boson stars, i.e. the = 0 case. For instance, given a value of , the mass M of the equilibrium configurations as a function of ω has a maximum, which gets larger as increases, yielding more compact stars. Furthermore, as in the case of 0-boson stars, the maximum value of the mass separates the space of solutions into two branches. These branches correspond to stable and unstable configurations against spherical perturbations, as shown in [37].
As mentioned above, the hypothesis that all the fields must have the same amplitude is essential to keep the spherical symmetry of the configuration. If one wants to consider different amplitudes of each constituent field, the assumption of spherical symmetry has to be relaxed. However, hitherto there has been no evidence that the resulting states may be equilibrium solutions of the Einstein-N -Klein-Gordon system. In this work we will show that deviations from spherical symmetry may indeed lead to new equilibrium solutions.
To proceed further with our non-spherical analysis we transform the solutions of the previous system of equations to Cartesian coordinates, x µ = (t, r, ϑ, ϕ) → x µ = (t, x, y, z). Then we perform a full non-linear numerical evolution of the perturbed stationary solutions.

III. DIAGNOSTICS
In order to test the stability of the static solutions we perform two different types of perturbations: (i) The first type consists in perturbing the energy density of the star given by ρ = n α n β T αβ , where n α is the four velocity of Eulerian observers in the 3+1 space time decomposition. The perturbed energy density is obtained by adding a non spherically symmetric small amplitude term to the homogeneous density in the following way [51]: where ρ 0 is the energy density of the equilibrium configuration, obtained from the solution of Eqs. (2.6)-(2.8), and R 99 is the radius enclosing 99% of the configuration's mass. In our simulations we choose κ = 0.1. This type of perturbations could trigger a potential bar-mode instability because it only affects the I xx and I yy components of the quadrupole moment defined as (3.2) Since the value of κ is small κ 1, this perturbation can be considered linear, initially; more importantly, it breaks the spherical symmetry of the original solution.
(ii) The second type of perturbations consist in varying separately the amplitude of each field. With these perturbations it is possible to study the stability of the stars against variations on each mode m and break the spherical symmetry. We choose the following form where φ is the unperturbed solution of the system of Eqs. (2.6)-(2.8). This perturbation introduces an additional constraint violation, besides the well known numerical error, but its magnitude is controlled by choosing a small , which, in general, depends on the , m-mode. Note that if is the same for all m, the perturbation is spherical. In order to assess the stability properties of the stars during the numerical simulation, we monitor the mass of the star, its angular momentum, and its density. We also follow the change in the quadrupole moment of the star, as shown below. Following the technique described in [51,52] to examine the stability of rotating neutron stars, we monitor the behaviour of the distortion parameter defined as which is a good measure of the magnitude of the bar-mode instability for perturbation (i). This parameter has been used to study the stability of rapidly, differentially rotating stars [51]. It has been observed that when the star is dynamically unstable, η z grows exponentially up to a maximum value; then the maximum value of η z remains constant on dynamical timescales. For stable stars, on the other hand, the maximum initial value of η z remains constant throughout the evolution. Thus the monitoring of η z provides a good tool to determine the properties of the star against bar-mode perturbations. In this work we also use as a measure of the deformation of the star. As further diagonostics, the maximum of the density and the lapse function are used to determine whether the configuration disperses or is undergoing a collapse. We have used the thorn AHFINDER [53] to follow the formation of an apparent horizon (AH) during the evolution. We have also computed the Hamiltonian constraint [54] to check the fourth order convergence of the implementation -see the Appendix for details on this procedure.

IV. TIME EVOLUTION AND NUMERICAL RESULTS
In this section we present the results from dynamical spacetime simulations from the perturbed -boson stars. We have compared the evolution of the equilibrium -boson stars with the perturbed stars.
We numerically integrate the EKG system using fourthorder spatial discretization within the EINSTEIN TOOLKIT framework. The EINSTEIN TOOLKIT solves the Einstein equations within the ADM 3+1 framework and evolves the spacetime using the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of the Einstein equations [55] through the MCLACHLAN thorn [56,57]. All the evolutions were made using the 1+log time slicing condition for the lapse α, and the Gamma-driver condition for the shift β i [58].
We use the Method of Lines thorn to solve the equations in time by using a fourth order Runge-Kutta scheme. The equations for the scalar fields are solved using a finite difference scheme of fourth order. We also employ the mesh refinement capabilities provided by the CARPET arrangements. The fixed mesh refinement grid hierarchy used consists of nested cubes with 3 levels of refinement. The finest is set in such a way that it covers the entire star.
We set the spatial resolution on the finest level to {dx, dy, dz} = 0.8 (and the coarsest to {dx, dy, dz} = 3.2) in order to fully capture the properties of the star. We follow the formation of an AH after the collapse of unstable stars.
More details on the resolution, as well as numerical convergence are given in the Appendix.
The three stationary configurations we chose to illustrate the general behaviour of the stars are represented with a square over the curve in Fig. 1 denoted as (M 1, M 2, M 3). Some of the properties of these stars are summarised in Table I. Both spherical and non-spherical perturbations to the stationary solutions are induced by increasing or decreasing the amplitude of the different constituent fields, see Eq.  perturbed. We use + or − to ascribe an increase ( > 0) or decrease ( < 0) of the amplitude, we use 0 to represent that no perturbation was introduced in that mode ( = 0). In this way, for instance, M 1 +0− means that M 1 has been perturbed in the following way: the first scalar field, φ 1,−1 , has been perturbed with > 0; the second field φ 1,0 has not been perturbed ( = 0), and the third field φ 1,1 has been perturbed with < 0.
In summary, we have perturbed M 1, M 2, M 3 in the following ways: perturbing all fields with the same amplitude as a test (spherical perturbation), introducing a non-axisymmetric bar-mode perturbation, and finally, we perturbed each constituent field using different amplitudes.

A. Spherical perturbation test
First, we perform numerical evolutions of the models listed in Table I with spherical perturbations. We induce perturbations in each field of a -boson star with all the perturbations having the same amplitude. In this way we guarantee that the spherical symmetry is preserved. This type of perturbations is done in order to compare and validate our results with those found using a spherically symmetric 1D code reported in [37]. While perturbing the initial equilibrium configurations adding perturbations (with a positive or negative value for ) that preserve the spherical symmetry, we find that the configuration that was reported to be stable in Ref. [37] (model M 1) remains stable in the timescale we reach in the 3D simulations, run M 1 000 . The values of the amplitude of the perturbations for these perturbed configurations are reported in Table II as M 1 +++ and M 1 −−− , in which we perturb each field adding or subtracting | | = 0.01 to each mode.
Our results are also consistent with models of -boson stars that are unstable in spherical symmetry. According to the results in [37], the configuration M 3 000 is unstable in the 1D simulations. When we perturb the amplitudes of the fields adding (run M 3 +++ ) or subtracting (run M 3 −−− ) the same amount, the configuration collapses or migrates to the stable branch respectively, as described in Table II. We monitor the behaviour of the metric coefficients during the evolution, and, in particular, we use the lapse and the formation of an AH as an indicator of the collapse of the star and the formation of a black hole.
The model M 2 deserves special mention since it corresponds to the critical solution: the star with maximum mass. We found that perturbations increasing the amplitude of the field ( > 0, run M 2 +++ ) make the star collapse whereas perturbations that decrease the amplitude ( < 0, run M 2 −−− ) drive the configuration to a new stable state as described in Table II. These results are consistent with the results reported in [37] for perturbations that increase of decrease the mass of the star.

B. Non-spherical perturbation: perturbing the energy density
In order to determine whether -boson stars develop a bar mode instability, we took as initial data a stationary model and modified the energy density in accordance with eq. (3.1). We have performed this analysis for configurations with = 0 and = 1, both stable against spherical perturbations. In the case = 0 we have taken the equilibrium configuration corresponding to ω/µ = 0.937 and for = 1 the configuration with ω/µ = 0.882, M 1 in Table. I. By choosing κ = 0.01, the momenta of inertia I xx and I yy change by less than 0.5% with respect to the equilibrium solution, hence we consider that the induced initial perturbation is small. Then we evolve the perturbed system via the Einstein-N -Klein-Gordon equations and monitor the behaviour of η z .
In Fig. 2 we show η z as a function of time for perturbed and unperturbed configurations for = 0 (top panel) and = 1 (bottom panel). For = 0, the distortion η z oscillates around zero for the perturbed case, indicating the star maintains, essentially, the spherical symmetry.
On the other hand, for = 1 in the case where the perturbation was included, the initial perturbation induces a small deviation from spherical symmetry therefore η z acquires a nontrivial value during the evolution. This non-zero value of η z indicates that the shape of the star deviates from spherical and becomes oblate.
During the evolution time considered (t ∼ 3500) we did not find any signal of a bar-mode instability for the models considered: no exponential growth in η z was measured. Most importantly η z does not grow as it happens for unstable stars [39]. This change in the shape of the = 1 configuration, illustrated by η z = 0, is compared with the case where M 1 . Bottom panel: the unperturbed configuration is a multi-field boson star ( = 1). In neither perturbed case has a bar instability been observed. Instead, a long lived departure from spherical symmetry occurs for = 1, but not for = 0, as long as ηz = 0. In contrast, the unperturbed configurations preserve spherical symmetry as ηz oscillates around zero in both cases.
is not perturbed (black solid line). For the unperturbed case, η z simply oscillates around zero. The conclusion, therefore, is that the perturbed configuration lingers, neither collapsing nor dissipating, thus showing a non-spherical distribution that is either stable or long-lived, without signs of instability. It is worth emphasising the key difference with the = 0 case, for which the evolution oscillates around a spherical distribution, in agreement with the fact that such a distribution is the only equilibrium configuration.
We found that after some time, the stars acquire a small linear momentum due to the numerical error and thus the deviation parameters can not be obtained accurately. Once this becomes noticeable we stop the evolution.   Table (III) . Departure from spherical symmetry is shown for those configurations that have been perturbed differently in all the three fields φ1,−1, φ1,0 and φ1,1. On the contrary, the configuration M 1000 has the same perturbation for all the fields (spherical perturbation), and it shows ηy = 0 at all times.
C. Non-spherical perturbation: perturbing the amplitude of each mode In this section, we describe the evolutions we have performed implementing non-spherical perturbation by varying the amplitude of each field of the -boson star.

Non-spherical perturbation of M 1
In order to illustrate the procedure to perturb the star, let us consider first model M 1. Different perturbations have been applied to M 1 and all them are summarised in Table III Tables II and III respec-tively. The mass, as expected, is increased for those perturbations with > 0 and decreases when < 0.
scripts indicate which fields have been perturbed and if the amplitude is increased by > 0 (subscript +), decreased by < 0 (subscript −) or it has been left without perturbation = 0 (subscript 0), as described before. All our evolutions show that M 1 remains stable without collapsing (none AH was found), independently of the perturbation. Thus, from the results summarised in Table III we can conclude that the configurations in the stable branch (against spherical perturbations), that is, to the right of configuration M 2 in Fig. 1, are also stable under non-spherical perturbations, against collapse.
Let us now turn to another result that can be extracted by studying the distortion parameters η z and η y . For spherical configurations and for those configurations that are spherically perturbed, these are zero. On the other hand, for non-spherical perturbations a small deviation from spherical symmetry is induced. In other words, non-trivial values of η y are obtained throughout the evolution of M 1. This behaviour of η y as a function of time is shown in Fig. 3. Indeed, non-zero values of η y for the evolution of M 1 +00 , M 1 00,+ , M 1 −00 , M 1 00− , M 1 0+0 , and M 1 0−0 , are obtained. No instability is observed, but the deformation does not die off either. These long-lived deformed configurations, arising as dynamical solutions of the Einstein-N -Klein-Gordon with three fields with different m are not spherically symmetric. The corresponding perturbations appear to be zero modes, suggesting a larger family of solutions.
As expected, we observe that the equilibrium configuration M 1 000 has η y = 0 at all times of the evolution. Besides, we have found that the behaviour of η y and η z during the evolution is the same for perturbations in the modes m = 1 and m = −1 with the same values of . This suggests that the resulting configurations are axially symmetric.
Finally, we report in Fig. 4 the time evolution of the total mass of M 1 under different spherical and non-spherical perturbations. As expected, the mass of the perturbed config- urations decreases or increases when < 0 or > 0. Notice, however, the change in mass is the same whether we perturb the m = −1 or m = 1 modes, for the same the sign of . This fact supports the assertion that the resulting configurations are axially symmetric.
The total mass of the models decreases with time, showing a small drift, even for the unperturbed solution. We have checked that this drift is due to the numerical error, since it is reduced when the grid resolution is increased (see Appendix A). Fig. 5 displays a series of snapshots of projections in the planes xy and xz of the energy density for the run M 1 0(−)0 (With a large amplitude in the perturbation in the mode m = 0) . The initial perturbation is introduced in the mode m = 0, decreasing the mass of the star and inducing a small deformation. Notice that we have taken the largest perturbation presented in this section (model M 1 0(−)0 ), so that the deformation can be appreciated in these projections of the energy density. At later times the system evolves and settles down into a configuration without collapsing or exploding. Although the configuration looks almost spherical, the value of η y at late times is slightly different from zero (η y ∼ 0.05). We have evolved this configuration for t ∼ 10000 and remains in the same state not showing any signs of instability, or returning to a -boson star.
At this point it is important to mention that perturbations in the mode m = 0 do not modify the value of the total angu-   This particular result was also obtained for the perturbations of the models that will be presented below.

Non-spherical perturbation of M 2 and M 3
The results of the previous section indicate that those configurations (M 1) that are stable under spherical perturbations do not show non-spherical growing modes. Furthermore, perturbations to the fields φ 1,−1 and φ 1,1 applied to those configurations provide evidence for zero modes, producing new equilibrium configurations that are dynamically stable, and have small departures from the spherical configurations. Now we are interested in studying non-spherical perturbations on configurations that might undergo gravitational collapse to a black hole. Those configurations are M 2 and M 3.
In particular M 2, as we have mentioned, corresponds to the critical configuration with maximum ADM mass. Configuration M 2, as is shown in Section IV A and reported in the literature [37], divides stable from unstable configurations. The latter are those configurations that can collapse into a black hole. We now go one step further and study them under nonspherical perturbations.
The list of perturbation applied to each field of M 2 is summarised in Table IV and the result of the evolution is reported in the fifth column of the same Table. The results confirm that M 2 is the configuration that separates stable from unstable configurations. Indeed, as it can be observe in the results of Table IV, all those configurations which have perturbations that increased the total mass of the configuration undergo a collapse, while those configurations which have been perturbed and reduced the total mass of the configuration did not collapse to a black hole. In this respect, run M 2 −0+ is of special interest. The perturbation did not change the total mass of the configuration, and the result of their evolution is that it did not collapse. Finally, we have considered non-spherical perturbations of model M 3. The results of the evolution of the different models studied are summarised in Table V. This configuration is on the so called unstable branch. The results mimic, to some extent, those observed for M 2. Perturbations that reduced the total mass of the configuration led to a migration to the stable branch. On the other hand, those that increased the mass of the configuration let to a collapse into a black hole. But a key difference is seen for M 3 −0+ . This perturbation did not change the total mass of the configuration, and contrary to M 2 −0+ , it did collapse to a black hole. This result, combined with the spherical perturbations mentioned in IV A, further confirms the special status of the maximum mass configuration M 2: it marks the threshold of unstable configurations of -boson stars, for both spherical and non-spherical perturbations.

V. DISCUSSION AND OUTLOOK
In this paper we performed dynamical simulations in the fully non-linear EKG model to investigate the stability of -boson stars. Unlike previous works we have considered non-spherical perturbations. An expected result is that those configurations known to be unstable under spherical perturbations, are also unstable under more general perturbations. The most interesting question, however, was if the configurations known to be stable under spherical perturbations would remain stable under more general ones. Here, our conclusions are two-fold. Firstly, no growing modes have been measured in our simulations. In this sense -boson stars are stable against non-spherical perturbations. However, when deformed away from sphericity, -boson stars do not return to a spherical state. They appear to oscillate around a new (slightly) non spherical state. We take this as evidence that new, multi-field, equilibrium configurations of the Einstein-N -Klein Gordon system exist, which are non-spherical. This conjecture is our second conclusion.
If our conjecture is proven correct, the spherically symmetric -boson stars are only an enhanced isometry point of a larger family of solutions of the Einstein-N -Klein Gordon. As discussed in the introduction, this is analogous to the Schwarzschild BH being the isometry enhancement point of the Kerr family. It is well known that the Kerr solution brings about qualitatively novel features with respect to the Schwarzschild solution. So it will be quite interesting to understand the novelties brought by the enlarged family of solutions that this work is suggesting.
The conjecture on the existence of these new non-spherical, multi-field configurations can be tested by solving the Einstein-N -Klein Gordon system for static or stationary (i.e. spinning) configurations, without assuming spherical symmetry. Research in this direction is already ongoing.
For run M 1 000 we report the time evolution of the mass and violations of the Hamiltonian constraint, with different resolutions, {dx, dy, dz} = 3.2, {dx, dy, dz} = √ 2 1.6 and {dx, dy, dz} = 1.6, where dx, dy and dz are the sizes of the coarsest level of refinement The L2-norm of the Hamiltonian constrain, given by where N is the number of points in the grid, is shown in Fig. 6. Here we conclude that the constraint equations converge as H reduces when the resolution is increased, the black (solid) and the blue (dashed) line have been multiplied by the factors 4 and 16, showing fourth order convergence. The low resolution (solid line, coarsest grid We plot the mass for run M 1 000 . As the resolution is increased the mass converge to a constant value and the overall drift is reduced. In Fig. 6 and Fig. 7 we have plotted until time equal to 1300, however the ({dx, dy, dz} = 3.2) simulation extends up to t ∼ 3500, where the final total mass differs from the initial value by 0.6%.