$\ell$-Boson stars

We present new, fully nonlinear numerical solutions to the static, spherically symmetric Einstein-Klein-Gordon system for a collection of an arbitrary odd number $N$ of complex scalar fields with an internal $U(N)$ symmetry and no self-interactions. These solutions, which we dub $\ell$-boson stars, are parametrized by an angular momentum number $\ell=(N-1)/2$, an excitation number $n$, and a continuous parameter representing the amplitude of the fields. They are regular at every point and possess a finite total mass. For $\ell = 0$ the standard spherically symmetric boson stars are recovered. We determine their generalizations for $\ell>0$, and show that they give rise to a large class of new static configurations which might have a much larger compactness ratio than $\ell=0$ stars.

It is within this line that we present a new class of static solutions to the Einstein-Klein-Gordon (EKG) system which generalize the known boson stars by considering a collection of an arbitrary odd number N of free complex scalar fields of the same mass. For N = 1 these configurations reduce to the standard, spherically symmetric boson star configurations. For N > 1, new static configurations are obtained which are still spherically symmetric if the collection of scalar fields are excited in an appropriate way.
Our approach is to consider that spacetime curvature is sourced by a collection of classical fields that compound a spherically symmetric configuration. The construction is based on a method introduced in Ref. [23], where it was shown that for the case in which the fields are real and massless, and as long as all the harmonics with given angular momentum number ℓ = (N − 1)/2 are excited with equal amplitude, the total stress energy-momentum tensor is spherically symmetric. In this way, it is possible to study the dynamics of scalar fields which, individually, have angular momentum, while keeping the spherical symmetry of the spacetime. As remarked in [23] this resembles a spherically symmetric kinetic gas where the individual particles may rotate but the configuration is spherical. The authors of Ref. [23] used this technique to study a possible effect of the angular momentum on the critical collapse of scalar field configurations.
In the present work we show that their method still works for a collection of complex, non-interacting massive scalar fields. The situation is similar to that in standard boson stars, where an internal U (1) symmetry "hides" the time dependency of the field avoiding Derrick's theorem [24,25]. If we extend the group of symmetry to higher values of N , then not only the time evolution, but also the angular dependency of the non-trivial harmonics with angular momentum number ℓ = (N − 1)/2 can be accommodated in the internal field space, leading to new static and spherically symmetric solutions to the EKG system. We refer to these objects as ℓ-boson stars in this paper. The main goal of this work is to determine their domain of existence and to compare their properties with the standard ℓ = 0 boson stars.
The remainder of this article is organized as follows. In Sec. II we derive the spherically symmetric field equations for the particular collection of scalar fields belonging to a fixed angular momentum number ℓ. Next, in Sec. III we demonstrate the local existence of solutions to the EKG system in the vicinity of the center. Based on a shooting algorithm starting from the local solution at the center, in Sec. IV we numerically solve the field equations to find the globally regular solutions describing the ℓ-boson stars, and discuss their behavior, including their density profile and mass curves for different values of ℓ. Conclusions and an outlook to future work are drawn in Sec. V, and relevant identities involving spherical harmonics that are needed in this article are derived in the Appendix at the end of the paper.
We use the signature (−, +, +, +) for the spacetime metric, and natural units such that = c = 1. The numerical solutions presented below are in Planck units, where we have also taken the gravitational constant equal to 1, G = 1.

II. FIELD EQUATIONS
In this section, we present the spherically symmetric EKG system that describes a collection of an arbitrary odd number of complex, non-interacting scalar fields Φ i , i = 1, . . . , N , of mass µ each, which are excited in an appropriate way consistent with the spacetime symmetries.
The stress energy-momentum tensor associated with such a collection of scalar fields is Here Φ * i denotes the complex conjugate of the field component Φ i , and ∇ µ is the covariant derivative with respect to the spacetime metric g µν . Notice that since the mass is the same for the N fields, and they do not interact between themselves, this expression is invariant under U (N ) transformations in the internal field space. We assume that the different fields in the configuration satisfy the Klein-Gordon equation individually, such that (∇ µ ∇ µ − µ 2 )Φ i = 0.
Following [23] we now consider solutions of the form where the angular momentum number ℓ = (N − 1)/2 is fixed in the configuration and m, which plays the role of the index i in Eq. (1), takes values m = −ℓ, −ℓ+1, . . . , ℓ. As usual Y ℓm denotes the standard spherical harmonics, and the amplitudes φ ℓ (t, r) are the same for all m. As shown in the Appendix this leads to a total stress energy-momentum tensor which is spherically symmetric. (The authors of Ref. [23] show a similar result valid for a collection of scalar fields which are, however, real and massless; for this reason we include here an independent proof.) Representing the spacetime metric in terms of the Misner-Sharp mass M (t, r) and the lapse function α(t, r), the line element is written as with r the areal radial coordinate and dΩ 2 the standard line element on the unit two-sphere. With this notation and the assumption in Eq. (2), the EKG system yieldṡ where κ ℓ := (2ℓ + 1) G is the (rescaled) gravitational coupling constant and where the source terms are given by 1 Here Π ℓ := α −1φ ℓ is the conjugate momentum associated with the scalar field, χ ℓ := γ −1 φ ′ ℓ , and a dot and a prime denote partial derivatives with respect to t and r, respectively.
In the following, we look for solutions with a harmonic time-dependency, with some real frequency ω (taken to be positive without loss of generality) and a real-valued function ψ ℓ (r). This provides a nonlinear eigenvalue problem for the quantities (M, αγ, ψ ℓ ) which can be written as For the particular case of N = 1, i.e. ℓ = 0, these equations reduce to the ones describing static boson stars with a zero angular momentum scalar field. These configurations have been studied extensively in the literature, and may be parametrized by the central value of the field and an excitation number [9][10][11][12].
In the following, we study the solutions for N > 1, i.e. ℓ > 0. In this case the presence of the centrifugal term ℓ(ℓ + 1)/r 2 in the EKG system requires ψ ℓ to vanish at the center of the configuration and to decay with a specific rate as r → 0. This is analyzed next.

III. LOCAL SOLUTIONS NEAR THE CENTER
For ℓ = 0, a full proof for the existence of regular solutions of finite mass was given by Bizoń and Wasserman in Ref. [26]. Although it would be very interesting to generalize their results to arbitrary values of ℓ, a full existence proof clearly lies beyond the scope of this work. However, for the purpose of this paper we found it useful to generalize one partial result given in [26] (see Proposition 3.1 of that paper), regarding the local existence of solutions near the center. This provides a rigorous discussion for the existence and behavior of the solutions in the vicinity of r = 0, and allows us to set up the correct boundary conditions to start the numerical integration scheme described in the next section.
In order to discuss the solutions in the vicinity of r = 0, it is convenient to introduce the following dimensionless quantities: 2 1 Note that when multiplied by the constant factor (2ℓ + 1)/(8π), these quantities correspond to the energy density, ρ = n µ n ν Tµν , the radial stress S = P µ r P ν r Tµν , and the radial momentum flux j = −P µ r n ν Tµν , respectively, of matter as measured in a local orthonormal frame by the Eulerian observers moving along the direction normal to the spatial hypersurfaces, where n µ is the unit normal vector to these hypersurfaces, and P µ r the orthogonal projection operator onto them. in terms of which the system (7) can be rewritten as Here we have used the shorthand notation f x = df /dx, and it is understood that γ 2 should be substituted by For ℓ = 0, it was shown in Ref. [26] that for each positive values of Q 0 and a 0 , there exists a local solution of the form near the center x = 0. In the following, we generalize this local existence result to arbitrary values of ℓ. As mentioned above, the presence of the centrifugal term ℓ(ℓ + 1)/x 2 drastically modifies the behavior of the fields close to the center. In fact, a heuristic approximation assuming which has solutions of the form f ∼ x ℓ or f ∼ x −ℓ−1 . Since we require regularity a the center, we discard the second possibility. We now show the existence of a two-parameter family of solutions for which , in terms of the original variables. Notice that these expressions are only valid as long as ℓ > 0. In order to prove the existence of such solutions, we replace f and its first derivative with the new fields (G, H), such that Note that (up to a constant 2ℓ + 1) the two asymptotic solutions f (x) = x ℓ and f (x) = x −ℓ−1 correspond to the fields (G, H) = (1, 0) and (G, H) = x −2ℓ−1 (0, 1), respectively. In terms of the fields u := (B, Q, G, H), the system (9) can be rewritten in first-order form: with the source terms given by In these expressions, γ 2 should be replaced with [1 − xB] −1 and the terms x −ℓ f and x 1−ℓ f x should be rewritten in terms of G and H according to Eq. (13). Note that when rewritten in this way, these terms are regular at x = 0, which implies that the source terms F B , F Q and F G are regular as well at x = 0 for all ℓ ≥ 1.
Eqs. (14) can be reformulated as a fixed point problem of the form u = T u, with the nonlinear map T u = (T 1 u, T 2 u, T 3 u, T 4 u) given by Introducing the function space X of continuous, bounded functions u on some interval [0, x 0 ] with the sup-norm, and the closed subspace Y ⊂ X consisting of those fields u ∈ X whose distance to the point u 0 := (0, Q 0 , a 0 , 0) is no greater than 1, it is not difficult to show that T leaves Y invariant and defines a contraction on Y , provided x 0 is small enough. According to the contraction mapping principle, T possesses a unique fixed point in Y , and this fixed point provides the required solution. The solution may be constructed by means of the iteration u 0 , u 1 := T u 0 , u 2 := T u 1 , . . ., from which one easily finds the asymptotic behavior (12).

IV. GLOBAL SOLUTIONS
In this section, we present numerical solutions to the EKG system described by Eqs. (7), that extend the local solutions of the previous section to provide global configurations that are regular everywhere and posses a finite mass.
From this point onwards, and in order to simplify the numerical analysis, we set G = 1, such that all quantities are dimensionless and measured in Planck units.
To proceed, we solve numerically the EKG system (7) expressed in the alternative form where u ℓ := ψ ℓ /r ℓ . The choice of appropriate boundary conditions must guarantee that the solutions are regular and asymptotically flat. Demanding regularity at the origin, i.e.
(see Sec. III for details), and a vanishing field amplitude at infinity, i.e. u ℓ (r → ∞) = 0, one obtains a nonlinear eigenvalue problem for the mode-frequency ω. Here u 0 ℓ is an arbitrary positive constant, and with no loss of generality we have fixed the value of the lapse function at the origin to one. Notice that since the system of equations is invariant under the transformation (α, ω) → λ(α, ω), with some positive arbitrary constant λ, one can always rescale the value of the lapse function at the end of the numerical integration in such a way that α(r → ∞) = 1, providing an asymptotic flat coordinate system at infinity. We always perform such a rescaling when reporting our results below.
A word about our numerical algorithm is on order here. The integration of the system is performed using a shooting algorithm to find the frequencies ω. Notice that, as mentioned above, we have rescaled the field as u ℓ := ψ ℓ /r ℓ , and we solve for u ℓ instead of ψ ℓ . The reason for this is that, for ℓ > 1, one finds that ψ ℓ and its first ℓ − 1 derivatives vanish at the origin, which results in a numerically ill behaved set of boundary conditions. On the other hand, u ℓ has a constant value at the origin and only its first derivative vanishes, which results in a numerically well behaved system of equations plus boundary conditions.
To find the solution one then integrates the system of Eqs. (17) outwards from the origin, with initial conditions given by Eqs. (18), and searches for the values of the frequency that match the required asymptotic behavior of the field (the solution should decay exponentially), until the shooting parameter converges to the desired accuracy. Notice that this process is somewhat delicate as for arbitrary ω there is always a solution that grows exponentially, and the shooting parameter must be precise enough to eliminate the growing solution. In practice we start with the outer boundary fairly close in, find a good initial guess for ω, and then slowly move the boundary outward, adjusting the value of ω to high numerical precision as we do so. We have in fact constructed two independent codes, with a fourth order and a fifth order Runge-Kutta numerical integration, and we have checked that the solutions of both codes converge with numerical resolution as expected, and that these solutions match.
For simplicity, in all our solutions we take the mass parameter µ = 1, although the solutions can be rescaled to an arbitrary value of µ using the invariance of the equation system under the transformation with the metric coefficients α and γ unchanged. We characterize the total mass of an ℓ-boson star in terms of the asymptotic value of the Misner-Sharp mass function, which is approximated by evaluating the metric coefficient γ(r) at the last grid point r max of the computational domain, such that Although ℓ-boson stars extend to infinity and do not posses a surface at a finite radius like usual fluid stars, one can associate to them an effective radius R(99%) defined as the areal radius of the object which contains 99% of the total mass. For a given angular momentum number ℓ, the equilibrium configurations are labeled by a continuous parameter representing the field amplitude, i.e. the constant u 0 ℓ in Eq. (18), and a discrete number n that labels the solutions  Fig. 1. Left panel: wave-function ψ1(r), and energy density ρ(r) = (2ℓ + 1)ρ (ℓ) (r)/(8π), as functions of the radial coordinate. Right panel: the metric coefficients γ(r) and α(r), also as a function of the radial coordinate. Notice that even though the wave-function vanishes at the origin, the energy density takes a value different from zero at that point. This is a characteristic of the ℓ = 1 configurations. See Fig. 3 for energy density profiles in the cases of ℓ = 0, 1 and 2. The regularity of these objects at the origin in evident in the profiles for the metric coefficients.   Fig. 1. Notice that the objects with a non-trivial angular momentum number are more compact than standard boson stars with ℓ = 0.
to the eigenvalue problem. In this paper we restrict our attention to ground state configurations only, n = 1, corresponding to the solutions with the lowest possible value of the frequency ω for a given angular momentum number and field amplitude. Each of these configurations have an associated frequency ω, a total mass M , and an effective radius R(99%). In Fig. 1 we show, for the equilibrium configurations with ℓ = 0, 1, 2, 3, 4, a plot for the total mass as a function of the effective radius (left panel) and frequency (right panel), respectively. Notice that the mass increases monotonically with the effective radius up to a maximum value, and then decreases, in analogy with what is known for standard ℓ = 0 boson stars. The configurations that correspond to this maximum value of the mass for the different values of ℓ have been labeled as A, B, C, D, and E in the figure. The main characteristics of those configurations are shown in Table I. Notice also that the compactness of these objects, defined as the ratio of the total mass to the effective radius, grows with the angular momentum number, at least for the first five possible values of ℓ considered in this paper. The relation of the mass with the frequency resembles some of the characteristics of rotating boson stars. Notice that since the centrifugal term decreases with increasing areal radius, the global characteristics of these objects approach those of standard bosons stars when they grow in size, the differences being only close to the center of the configuration.
A representative ℓ-boson star configuration appears in Fig. 2, where we show, for the configuration B of Fig. 1, the profile of the wave-function ψ 1 (r), introduced in Eq. (6), the energy density ρ(r) = (2ℓ + 1)ρ (ℓ) (r)/(8π), where ρ (ℓ) was given in Eq. (5a), and the metric coefficients γ(r) and α(r) of Eq. (3), as functions of the areal radius. Notice that, contrary to what happens for the standard ℓ = 0 boson stars and perfect fluid configurations, the energy density at the origin does not take a maximum value that decreases monotonically to zero when the angular momentum number is non-trivial. This is a consequence of the regularity conditions discussed in Sec. III. In particular, the energy density takes a local minimum at the origin for the configurations with ℓ = 1, and vanishes if ℓ > 1. See Fig. 3 for details.  Fig. 1, respectively. Notice that the vertical scale changes between the plots. It is interesting to emphasize that ρ(0) = 0 if ℓ = 0 or 1, but ρ(0) = 0 if ℓ > 1. Only if ℓ = 0 the energy density takes its maximum value at the origin. Right panel: Same profiles as in left panel, but multiplied by the square of the areal radius. Even if the amplitude of the energy density decreases as the angular momentum number increases, the configurations become more extended, and then the combination ρr 2 grows, leading to more massive objects; see Table I for details.

V. DISCUSSION AND CONCLUSIONS
We have introduced new solutions to the static EKG system of equations that describe regular self-gravitating objects of finite mass. These configurations incorporate some of the effects of angular momentum while keeping the spherical symmetry of the spacetime metric, generalizing standard non-rotating boson stars. The technique requires an arbitrary odd number N of massive, complex, non-interacting scalar fields with an internal U (N ) symmetry, and are labeled by an angular momentum number, ℓ = (N − 1)/2, and the usual field amplitude and excitation number. This motivates the name ℓ-boson star. In this paper we have restricted our attention to ground state configurations only, although this can be easily generalized to the case of excited states.
We have shown that these objects possess similar properties to those of standard N = 1, ℓ = 0, boson stars. Specifically, for a given angular momentum number, the equilibrium configurations exhibit a maximum value of the mass, and this maximum grows as ℓ increases, leading to more compact objects. For the case of ℓ = 0, it is well known that the maximum mass configuration divides the stable and unstable branches. We have done some preliminary numerical evolutions for the case with ℓ = 1 and have found that the same appears to be true: for a given mass below the maximum, configurations with a large effective radius and higher frequencies (those to the right of the maximum in both panels of Fig. 1) remain stable when perturbed, while those with smaller effective radius and smaller frequencies (to the left of the maximum in Fig. 1) either collapse rapidly to a black hole, or oscillate with a very long period around a less compact state (depending on the precise way in which they are perturbed). We will report on the results of these simulations and those for ℓ > 1 elsewhere.
Even if in this paper we restricted our attention to the ground state n = 1 and only one single value of ℓ at a time in the configurations, the system of Eqs. (7) can easily be generalized to the case when several different modes (n, ℓ) are present simultaneously. This can be achieved by summing over all the excited states in the right hand side of Eqs. (7a) and (7b), similarly to the multi-state configurations constructed in [27,28] for ℓ = m = 0. Under the classical approach that we assume in this paper, this requires the inclusion of more fields, making it possible to use not only odd, but also even values of N , and to obtain more general density profiles than those reported here. We leave a more detailed analysis for a companion paper.
The main implication of this work is that the solution space of spherically symmetric, static boson stars is much larger and has a much richer structure when a collection of scalar fields is considered instead of a single complex scalar field that necessarily requires of the harmonic ℓ = m = 0. An alternative interpretation for the collection of scalar fields considered here, and potential applications for dark matter observations, will be explored in another paper.
Based on the identities (A1-A7), one easily finds the following expressions for the components of the stress energymomentum tensor given in Eq. (1): where here a, b = t, r, A, B = ϑ, ϕ andg ab refers to the components of the time-radial part −α 2 dt 2 + γ 2 dr 2 of the metric.