Charged boson stars revisited

We consider again stationary solutions to the spherically symmetric Einstein–Maxwell–Klein–Gordon system, commonly known as “charged boson stars”, originally studied by Jetzer and Van Der Bij (Phys Lett B 227:341, 1989). We construct families of charged boson stars in the ground state, for different values of the charge parameter q, and different values of the central scalar field. Following Jetzer and Van Der Bij, one can define a critical value for the charge q=qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=q_c$$\end{document} that corresponds to the value for which the Coulomb repulsion of the bosonic particles exactly cancels their newtonian gravitational attraction. We confirm the claim made by Pugliese et al. (Phys Rev D, 2013. https://doi.org/10.1103/physrevd.88.024053) that super-critical solutions exist for a limited range of charges above the critical value q>qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q>q_c$$\end{document} (though we find an even smaller range of q for which this is possible). Our analysis indicates, however, that all such super-critical solutions are gravitationally unbound, and are therefore expected to be unstable. One of the main results of our analysis is the fact that, even though we do find a family of slightly super-critical solutions in the sense that q>qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q>q_c$$\end{document}, there are no super-critical solutions in the sense that the total charge Q is larger than the total mass M of the system.


I. INTRODUCTION
Boson stars are self-gravitating soliton-type configurations for a massive complex scalar field in general relativity (see [3][4][5][6][7] for reviews).These compact objects have been extensively studied since the pioneering works of Kaup in 1968 [8], where they are referred to as a "Klein-Gordon Geon", and of Ruffini and Bonazzola in 1969 [9].In 1987, the same compact object was studied by Friedberg, Lee, and Pang [10], where it was called a "mini-soliton star".Boson Stars can be considered descendants of the self-gravitating photonic configurations called geons (gravitational electromagnetic units) proposed by Wheeler in 1955 [11].Standard boson stars (also known as mini-boson stars) are described by a massive complex scalar field that is localized in a compact spatial region and is supported by its self-gravity.Complex scalar fields allow compatibility with a static geometry of spacetime, such that the scalar field contains two degrees of freedom, both oscillating harmonically in time, but out of phase.In this way, it is possible to evade the no-soliton Derrick Theorem [12].The case of real scalar fields has also been studied (see e.g.[13]), though in this case one can only form quasi-stationary configurations known as "oscillatons".
Although boson stars remain as purely theoretical, interest in these self-gravitating compact objects has recently increased due to developments in both particle physics and cosmology, such as the confirmation of the Higgs boson [14], suggesting that in the early stages of the universe boson stars may have formed out of fundamental scalar fields and could play a role in understanding the origin of dark matter.Furthermore, it has been shown that boson stars can be considered as candidates for black hole mimickers [15,16].Different types of bosonic structures have also been studied, such as the so-called "Proca stars", which are analogous to boson stars but for the case of a massive complex abelian vector field minimally coupled to gravity.In a similar way to the scalar field, the massive complex abelian vector field can also form self-gravitating solutions to the Einstein-Proca system [17][18][19].Another class of static solutions to the Einstein-Klein-Gordon equations recently studied are the -boson stars, which incorporate the effects of angular momentum while maintaining the spherical symmetry of the spacetime [20,21].
In this work we focus on the case of charged boson stars, which are spherically symmetric self-gravitating solutions for a massive complex scalar field coupled to the gauge group U (1): the Einstein-Maxwell-Klein-Gordon (EMKG) system.Charged boson stars were first studied by Jetzer and Van Der Bij in [1], and later by Pugliese et al. in [2].More recently, the EMKG system has also been studied in different contexts, such as the case of fermion-charged-boson stars in [22], or the gravitational collapse of charged matter in [23].Charged boson stars can be characterized by the charge parameter q associated to the bosonic conserved current.In [1] Jetzer and Van Der Bij argue that there is a critical value for the charge q c such that charged boson stars can only exist for q < q c .In the original units of [1] (see discussion at the end of Section II below), the critical value is given by q c = 1/ √ 2(m/M P L ) 0.707(m/M P L ), with m the mass parameter of the complex scalar field and M P L = c/G the Planck mass.In newtonian terms, this critical value for the charge corresponds to the case when the Coulomb repulsion exactly cancels the gravitational attraction of the fundamental bosonic particles, so that for super-critical values of the charge one would not expect to find stationary solutions.However, it was claimed in [2] that solutions with super-critical charge up to q ∼ 0.8 (m/M P L ) > q c can in fact be found.
In this work we study again in detail the solutions for charged boson stars, and show that indeed there exits a family of solutions corresponding to charged boson stars with a super-critical charge q > q c .We find, however, that some of the super-critical solutions of Pugliese et al. are not correct, apparently due to the fact that the spatial range they considered did not extend far enough, so that when the spatial domain is extended the scalar field does not decay exponentially and the resulting space-time is not asymptotically flat.In our case we can only find super-critical solutions up to q ∼ 0.739 (m/M P L ).However, one of our main results is the fact that, even though we do find a family of slightly super-critical solutions in the sense that q > q c , there are no super-critical solutions in the sense that the total charge Q is larger than the total mass M of the system (see Section VI below).
In order to solve the EMKG system we use a 3+1 decomposition of the electromagnetic field as described in [23,24], and construct families of stationary solutions for different values of the charge parameter q.For a given value of q, the system of equations can then be cast as an eigenvalue problem for the oscillation frequency of the scalar field ω, with the value of the scalar field at the origin as a free parameter.We have also performed short time evolutions of some of our solutions in order to verify that the frequencies obtained in our analysis do correspond to the frequencies observed during a dynamical evolution.For these evolutions we have used a fully non-linear time evolution code for numerical relativity in spherical symmetry, the OllinSphere code previously described in [23,25,26], that uses the Baumgarte-Shapiro-Shibata-Nakamura formulation adapted to spherical symmetry [27][28][29][30].However, we will leave a full analysis of these time evolutions for a future work.
This paper is organized as follows.In Section II we present the Einstein-Maxwell-Klein-Gordon system, and also discuss the units conventions taken for the electromagnetic field, as well as the normalization that we use for the charge parameter q.In Section III we derive the spherically symmetric field equations for stationary configurations corresponding to the charged boson stars.In Section IV we define the total mass, total charge, and binding energy of the system.Section V describes the boundary conditions and the numerical methods that we have used to ensure that the solutions obtained are asymptotically flat.Section VI presents our results for the different families of solutions.We conclude in Section VII.

II. THE EINSTEIN-MAXWELL-KLEIN-GORDON SYSTEM
A self-gravitating charged massive complex scalar field φ is described by the following action (we use a metric signature (−, +, +, +) and Plank units such that G = c = = 1): where g is the determinant of the space-time metric, R its associated Ricci scalar, m is the scalar field mass parameter, F µν is the electromagnetic Faraday field tensor (2-form) given in terms of the potential 1-form A µ as: and D µ is the gauge invariant covariant derivative: with ∇ µ the usual space-time covariant derivative and q the scalar field charge parameter.The action (1) is invariant under local U (1) gauge transformations of the form: with θ(x α ) a local gauge function which depends on the space-time coordinates.This implies the existence of a conserved Noether 4-current, which acts as a source for the electromagnetic field and is given by: The energy-momentum tensor, which is the source of the gravitational field, includes contributions from both the complex scalar field and the electromagnetic fields, T µν = (T φ ) µν + (T F ) µν , where the scalar and electromagnetic contributions are respectively: Variation of the action (1) with respect to the different fields results in the Einstein field equations: with G µν the Einstein curvature tensor of the space-time, plus the Mawxell equations: where F * µν := − µναβ F αβ /2 is the dual electromagnetic tensor (using the convention that 0123 = −1/ √ −g and 0123 = + √ −g), and the Klein-Gordon equation for the scalar field: Notice that this is the usual Klein-Gordon equation except for the fact that it involves the gauge covariant derivatives defined above in Eq. (3).
A word about our conventions for the electromagnetic field is in order here.While quite common in the general relativity community, our conventions are by no means universal and many references, including the work of Jetzer and Van Der Bij in [1], define the electromagnetic Lagrangian with a factor 1/4 instead of the 1/16π used in the action (1) above.This implies that the definitions for the potential 1-form and the Faraday tensor, as well as the electric and magnetic fields that we will introduce below, have an extra factor of 1/(4π) 1/2 with respect to our convention, Âµ = A µ /(4π) 1/2 .In order to keep expressions like the gauge covariant derivative unchanged this requires one to also rescale the charge as q = (4π) 1/2 q.This implies, for example, that while with our convention the Newtonian force between two charged particles is simply F = q 1 q 2 /r 2 , with the alternative convention it would be F = q1 q2 /(4πr 2 ).Similarly, while with our convention a maximally charged Reissner-Nordström black hole would have a total charge equal to its total mass, Q = M , with the alternative convention this would correspond to Q/(4π) 1/2 = M .In [1] Jetzer and Van Der Bij also do not use the 1/2 factor in the Lagrangian for the scalar field, so that their scalar field is rescaled with respect to ours as φ = φ/ √ 2 (they also use an opposite signature for the metric).Finally, in order to simplify further their expressions, Jetzer and Van Der Bij later rescale the charge as q = (q/m)/[(8π) 1/2 ] = (q/m)/ √ 2 (equation (14) in reference [1]), which implies that while the critical charge for boson stars is expected to be q c = m with our convention, in their convention this corresponds to qc = 1/ √ 2 (remember that we use Planck units so that M P L = 1).This convention has later been followed by other authors [2,22], so we will use it below when we present our numerical results en section VI.

III. CHARGED BOSON STARS
In order to study stationary configurations corresponding to charged boson stars we assume spherical symmetry and the usual harmonic ansatz for the time dependence of the scalar field: with ω a real constant corresponding to the frequency of oscillation of the scalar field, and φ 0 (r) a real-valued radial function which corresponds to the profile of the charged boson star [1,2,5].Both contributions to the energymomentum tensor given by Eqs. ( 6) and ( 7) are independent of time under this harmonic ansatz.As a consequence, the space-time metric turns out to be static.This implies that one can write the metric in spherical coordinates (r, θ, ϕ) in the polar-areal gauge as: with α and A functions only of the radial coordinate r, and where dΩ 2 = dθ 2 + sin 2 (θ)dϕ 2 is the standard solid angle element.
In the following we will consider the 3+1 formalism of general relativity, where the space-time (M, g µν ) is assumed to be globally hyperbolic, so that it can be foliated by a family of spacelike hypersurfaces Σ t that are parametrized by a global time function t, that is M ∼ = R × Σ t [31,32].In this formalism the metric given by ( 12) corresponds to taking a null shift vector β i = 0 and a spatial metric of the form γ ij = diag(A, r 2 , r 2 sin(θ)), with α(r) the lapse function [23].With this form of the metric each spacelike hypersurface Σ t has a normal timelike unit vector field n µ such that: The timelike unit vector n µ can be identified with the 4-velocity of observers moving along the normal direction to the spacelike hypersurfaces, the so-called Eulerian observers.

A. Gauss constraint
Within the 3+1 formalism the covariant equations for the electromagnetic field can be formulated in terms of the electric and magnetic fields as measured by the Eulerian observers (see e.g.[24]).In terms of the Faraday tensor F µν and its dual, the electric E µ and magnetic B µ fields are defined as: Since the electric and magnetic fields are purely spatial one can consider only their spatial components E i and magnetic B i .Moreover, in spherical symmetry the electric field only has a non-zero component in the radial direction E i = (E, 0, 0), while the magnetic field vanishes identically B i = 0.This implies that the potential 1-form can be taken to be of the form A µ = (A 0 , 0, 0, 0) [1,22].
Projecting the covariant Maxwell equations ( 9) onto the normal vector n µ we obtain two constraint equations to solve for the initial data [24].However, the constraint equation for the magnetic field is trivial since this field vanishes.The only Maxwell equation to solve is therefore the Gauss constraint, which now takes the form [23]: where D i is the covariant derivative compatible with the spatial metric γ ij , and e := −n µ j µ is the electric charge density measured by the Eulerian observers: Using now our harmonic ansatz for the scalar field, Eq. ( 11), the Gauss constraint can be written in terms of the metric functions and the electric field explicitly as: where we have defined the electric scalar potential as F := αΦ, with Φ := −n µ A µ .Since the potential 1-form has only a non-zero time component we immediately find F = −A 0 .From this and Eq. ( 2) one also finds: This is just the generalization of the usual expression for E as minus the gradient of the electric potential to the case of our curved space-time.

B. Klein-Gordon equation
The Klein-Gordon equation ( 10) can be reduced to first order form by defining the variables: with P µ ν = δ µ ν + n µ n ν the projection operator onto the spatial hypersurfaces.Using the metric (12) and the harmonic ansatz (11) one now finds: where now: From this, the Klein-Gordon equation can be rewritten as a first order differential equation for χ:

C. Hamiltonian constraint and slicing condition
Since the space-time is static the extrinsic curvature vanishes K ij = 0.This implies that the momentum constraints are trivial so we only need to solve the Hamiltonian constraint, which in this case reduces to: where (3) R is the Ricci scalar associated with the spatial metric γ ij , and ρ := n µ n ν T µν is the local energy density measured by an Eulerian observers.Using the energy-momentum tensor given by Eqs. ( 7) and ( 6), the Hamiltonian constraint takes the form of a first order differential equation for the radial metric function A(r): We still need to find an equation for the lapse function α(r).Since we are working in the polar-areal gauge [33,34], the natural slicing condition is precisely the so-called polar slicing condition ∂ t K θθ = K θθ = 0, which in this case reduces to:

IV. TOTAL MASS AND CHARGE
When working in terms of the areal radius it is well known that the total mass M of the spacetime can simply be found as integral over a flat volume element of the form: where ρ is the energy density of matter that we introduced before, ρ = n µ n ν T µν .The above result can be shown directly from the Hamiltonian constraint.In our case, the energy density ρ has contributions both from the scalar field and the electromagnetic field, and takes the form: The mass integral (26) defined above, though correct, nevertheless has a serious drawback due to the fact that for a charged boson star (or indeed for any charged particle) the electric potential F decays as 1/r, so that the integral converges to the total mass M rather slowly with r.We will come back to this problem below.
The total electric charge Q is provided by the conserved Noether charge which is defined by the local U (1) symmetry.Integrating the time component of the conserved current j µ we obtain: with γ the determinant of the spatial metric.Substituting the value of j 0 using Eq. ( 5), and the harmonic ansatz, we obtain the following explicit expression for the charge integral: In contrast to the total mass, the charge integral above converges very rapidly since for a boson star the scalar field φ 0 decays exponentially (see following Section).It is also interesting to notice that while the mass integral (26) involves a flat volume element, the charge integral (29) involves the full curved space volume element, hence the factor A 1/2 that appears in (29) but not in (26).
Since the total charge is related to the total number of particles N as Q = qN , the same integration allows us to find N [1,7].This is quite useful as one can use the total mass M and the total number of particles N to define a binding energy for the star as: In order to understand this definition notice first that the total mass M includes all possible contributions to the energy of the star, that is it includes the rest-mass plus the kinetic and potential energies.So, if we subtract from M the total rest-mass given by mN , we are left with just the kinetic (positive) and potential (negative) contributions, which is precisely the definition of the binding energy.If the binding energy is negative the system is gravitationally bound, while if it is positive the system is not bound, and even very small perturbations can cause it to dissipate to infinity.
One can also find the total mass of the system in an alternative way by assuming that far away the metric reduces to the Reissner-Nordström metric, so that: Solving for M we then find: Having first found the total charge Q using (29), we can use the above expression to find the total mass M .It turns out that this expression in fact converges very rapidly with r, since once we are in a region where the scalar field is negligible the space-time reduces to the pure electro-vac Reissner-Nordström solution.Because of this we prefer to use (32) to obtain the total mass instead of the integral ( 26), but we have indeed checked that for our solutions both expressions agree for very large r.
As already mentioned, for a boson star the scalar field decays exponentially so that there is no real surface.However, we can use the charge integral above to define an effective radius R 99 as that which contains 99% of the total charge.We could in principle do the same with the mass integral and define a (somewhat different) effective radius R 99 that contains 99% of the total mass, indeed this is what is usually done for boson stars with no electric charge.However, in light of the discussion above regarding the convergence rates of the mass and charge integrals, for charged boson stars it is in fact much better to define R 99 in terms of the charge integral, and this is what we will do when we report our results below.

V. BOUNDARY CONDITIONS, RESCALING, AND NUMERICAL METHODS
The system of equations to be solved in order to construct charged boson star configurations consists on the Gauss constraint (17), the Klein-Gordon equation (22), the Hamiltonian constraint (24), and the polar slicing condition (25), together with the equation for the scalar potential (18), and the definition of χ (21), for the six functions {E, F, φ 0 , χ, A, α}.
To solve this system we must also choose appropriate boundary conditions to ensure that the solutions are regular at the origin and that the space-time is asymptotically flat.For the boundary conditions at the origin we take: with k > 0 a positive real constant.The vanishing of the radial derivatives at the origin is due to the spherical symmetry.Asking for A(0) = 1 is required in order to guarantee that the space-time is locally flat there, while the constant value of φ 0 (0) = k is our free parameter.On the other hand, choosing α(0) = 1 and F (0) = 0 is done just for simplicity, since we don't know the correct values of those variables there (though one could argue that there are no "correct" values there since these are just gauge functions).But we will have more to comment on these choices below.Notice in particular that with these conditions we also have For solutions that represent an isolated star the scalar field must also vanish at infinity, that is φ 0 (r) → 0 for r → ∞.For each choice of φ 0 (0), our system of equations has solutions that decay exponentially at infinity only for certain frequencies ω.This means that, given a central value of the scalar field, we must solve an eigenvalue problem in order to find the frequency ω.Notice that the Klein-Gordon equation ( 22) implies, in particular, that for large r the following condition must be satisfied in order to have exponentially decaying solutions for the scalar field: with α ∞ and F ∞ the asymptotic values of α and F .If this condition is not satisfied one would have instead sinusoidal solutions for the scalar field for large r, that are simply not compatible with an asymptotically flat space-time.Given a value of φ 0 (0) = k as a free parameter, we choose a trial value of the frequency ω and integrate our system of equations outwards from the origin using standard fourth order Runge-Kutta.We then use a shooting algorithm to find the correct value of ω that corresponds to exponential decay of the scalar field far away (the value of ω found in this way is typically such that ω > m, but this changes once we apply the rescaling described below).We also look for solutions with no nodes in the scalar field, corresponding to the ground state of our charged boson stars.One can also solve for excited states with one or mode nodes, but we will discuss such solutions elsewhere.
Going back to our boundary condition for α at the origin, we now notice that in the final solution we in fact do not want α(0) = 1, but rather α(r) → 1 at infinity, corresponding to Minkowski space-time.But this is no problem as the polar slicing condition (25) is linear in α, so we can always just rescale the lapse.However, in order not to affect the solution, we must also rescale the frequency ω and the scalar electric potential F by the same factor, as the whole system of equations is easily shown to be invariant under the change: with C 1 an arbitrary constant.We choose the constant C 1 by extrapolating the value of the lapse α at infinity assuming an asymptotic behavior of the form α ∼ α ∞ + cte/r, so that C 1 = α ∞ .This ensures that after rescaling the lapse will now go to 1 at infinity.For a non-charged boson star this rescaling gives us the final "physical" frequency, but for the charged case we are not yet done.Our solution was also found by asking for the electric potential to satisfy F (0) = 0, and this remains true even after the rescaling above.But it would seem much more natural to ask instead for F (r) → 0 at infinity.We can fix this by making a gauge transformation as in Eq. ( 4), with an appropriately chosen gauge function of the form θ = C 2 t, with C 2 constant.The gauge transformation in this case simplifies to: As before, we now choose the constant C 2 by extrapolating the value of F at infinity assuming an asymptotic behavior of the form F ∼ F ∞ + cte/r, so that C 2 = −F ∞ .The gauge transformation above clearly implies that the frequency must also be transformed as: After these two transformations our final solution is now such that α(r) → 1 and F (r) → 0 at infinity, as desired.
For a given value of φ 0 (0) = k we then have three different values for the frequency: an initial value ω 1 obtained by the shooting algorithm using our original boundary conditions α(0) = 1 and F (0) = 0; a second value ω 2 obtained after rescaling the lapse; and a final "physical" value ω 3 obtained after performing the gauge transformation.In all our figures and tables below we always report this final value for the frequency.Notice also that since in our final solution the lapse function α goes to 1 at infinity and the electric potential F goes to 0, condition (34) reduces simply to ω ≤ m.We find that indeed this is always the case for all our solutions.M q = 0.0 q = 0.2 q = 0.4 q = 0.6 q = 0.7 q = 1/ √ 2 q = 0.72 q = 0.73 q = 0.735 q = 0.739 amplitudes.As we increase the charge the allowed range of values for φ 0 (0) becomes smaller and smaller, until for q 0.739 it disappears completely and no more solutions are found.
We should emphasize here that when finding solutions one must really make sure that they decay exponentially far away by moving the numerical integration boundary to large radius.Is is quite easy to "find" what seems to be a nice decaying solution with the boundary relatively close by, only to discover that when we move the boundary further out the apparent exponential decay turns into a sinusoidal oscillation with small amplitude and large wavelength.Reference [2], for example, reports finding supercritical solutions with charge as high as q = 0.8, which according to our results would require ω > m, and which we simply have been unable to reproduce.We believe that such apparent solutions are no real solutions at all, and would have an asymptotic sinusoidal behavior if they were to be extended to large radius.The solutions in [2] are also shown only up to a radius r ∼ 10, which is still quite small for many of the configurations studied here.This is confounded by the fact that they seem to have rescaled their frequencies using the behavior of the lapse and scalar potential at a finite radius (presumably the boundary of their computational domain) instead of their asymptotic behavior at infinity, which makes it very difficult to reproduce the frequencies they report.
Figure 2 shows the total mass M (in units of M 2 P L /m) for the different families of solutions as a function of the central value of the scalar field φ 0 (0).For boson stars with no charge, q = 0, it is well known since the work of Kaup in 1968 [8], and Ruffini and Bonazzolla in 1969 [9], that the mass has a maximum value M Kaup ≈ 0.633 M 2 PL /m.For central values of the field to the right of this maximum the total mass oscillates slightly and then converges to a value M ∼ 0.37 M 2 PL /m for very large φ 0 (0).From the Figure it is to easy see that this behavior of the mass is in fact very similar for charged boson stars as long as q < 1/ √ 2. The main change is that the maximum mass increases with q, while its position moves to lower values of φ 0 (0) (this was also found by Jetzer and Van Der Bij in [1]).Perhaps more interesting is the fact that the position of the maximum moves to φ 0 (0) → 0 as the charge approaches the critical charge q → 1/ √ 2, while the value of the maximum mass diverges M max → ∞.This might seem surprising at first, but we must remember that as the central value of the field approaches zero the effective radius of the boson stars increases without bound, so even if the energy density decreases the total mass can still grow.As mentioned before, we also find solutions for slightly super-critical charges in a narrow range of φ 0 (0).For such super-critical solutions the mass has a local minimum within this allowed interval, and no local maximum is found.
For boson stars with no electric charge, q = 0, the central value of φ 0 (0) for which the mass reaches its maximum is known to separate stable configurations (to the left) from unstable ones (to the right).One would expect a similar thing to happen for charged boson stars.This suggests that all our super-critical solutions, with no maximum and only a local minimum, should be unstable.Of course, in order to be sure of this one would need to do either a linear stability analysis or a full non-linear dynamical evolution (we will consider the dynamical simulation of charged boson stars in a future work).
In Figure 3 we show the relation between the total mass M and the effective radius R 99 .We have divided this Figure into two separate plots in order to better appreciate the changes in behavior close to and above the critical charge.As before, from the Figure it is clear that although the mass and effective radius increase with q, the behavior is similar for charges such that 0 ≤ q < 1/ √ 2 (lower panel).On the other hand, for super-critical charges 1/ √ 2 ≤ q < 0.739 PL /m q = 0.7 q = 1/ √ 2 q = 0.72 q = 0.73 q = 0.735 q = 0.739 M q = 0.0 q = 0.2 q = 0.4 q = 0.6 Figure 3.Total mass (in units of M 2 P L /m) versus effective radius R99 for the different families of boson stars (notice that the vertical axis is on a logarithmic scale).Lower panel: sub-critical charge.Upper panel: super-critical charge.We have included the sub-critical case q = 0.7 in the upper panel in order to show that for q 0.7 the total mass M for all configurations is always larger than the maximum for a q = 0 boson star MKaup.
the behavior changes and there in no longer a local maximum for the mass (upper panel).The system with q = 0.7 was included in the upper panel in order to show that for q 0.7 the total mass M for all configurations is always larger than the maximum mass for a q = 0 boson star M Kaup .
Figure 4 shows the compactness of the boson stars defined as C 99 := M/R 99 , as a function of the central value of the scalar field φ 0 (0).Again, for sub-critical charges the behavior is very similar to that for standard q = 0 boson stars.Notice, however, that for the family with critical charge q = 1/ √ 2 the compactness reaches a maximum of C 99 → 0.25 as φ 0 (0) → 0 (compare this with the compactness for a Schwarzschild black hole C = 1/2, and for a maximally charged Reissner-Nordström black hole C = 1).This might seem counter intuitive as in the case φ 0 (0) = 0 there is no boson star, but remember that as we approach that limit the mass and effective radius both diverge.This maximum compactness turns out to be an upper limit for all configurations.Notice also that for super-critical solutions with q > 1/ √ 2 the maximum compactness falls again.In Figure 5 we show the binding energy E B := M − mN for the different families of solutions, with M the total mass and N the total number of bosons.Solutions with negative binding energy are gravitationally bound, while those with positive binding energy are unbound.Notice that gravitationally bound solutions can be either stable or unstable depending on whether the value of φ 0 (0) is to the left or the right of that which corresponds to the maximum mass for that family.On the other hand, unbound solutions with E B > 0 are all expected to be unstable.From the Figure we can see that all sub-critical families have a region with negative binding energy, though this region becomes smaller and smaller as we approach the critical charge.However, for all solutions with q ≥ 1/ √ 2 the binding energy is always positive indicating that all such solutions are gravitationally unbound and therefore almost certainly unstable.From dynamical simulations of standard q = 0 boson stars we know that unstable but bound configurations can either collapse to a black hole or migrate to a stable solution when perturbed, whereas unstable and unbound configurations either collapse to a black hole or disperse (explode) away to infinity.We expect a similar behavior for C 99 q = 0.0 q = 0.2 q = 0.4 q = 0.6 q = 0.7 q = 1/ √ 2 q = 0.72 q = 0.73 q = 0.735 q = 0.739 q = 0.0 q = 0.2 q = 0.4 q = 0.6 q = 0.7 q = 1/ √ 2 q = 0.72 q = 0.73 q = 0.735 q = 0.739 Figure 5. Binding energy (in units of M 2 P L /m) for the charged boson stars configurations as a function of the central value of the scalar field φ0(0).We can see that all sub-critical families have a region with negative binding energy, while for configurations with q ≥ 1/ √ 2 this region disappears.
charged boson stars, but this will be studied in detail in a future work.There is another interesting relation one can find.Notice that the total number of bosons N can be written in terms of the total charge Q as N = Q/q, so that the binding energy becomes . This means that for gravitationally bound solutions we must have Q/M > q/m, while for the unbound solutions we have Q/M < q/m.For consistency, we have also renormalized the total charge with a factor of √ 2M , so that we now define Q := (Q/M )/ √ 2. Notice that q has already been defined as q = (q/m)/ √ 2, so that comparing Q/M with q/m is now equivalent to comparing Q directly with q.
Figure 6 shows the renormalized charge Q as a function of the central value of the scalar field φ 0 (0).We have divided Figure 6 into two separate plots.In the right panel corresponds to families of boson stars with a sub-critical charge q < 1/ √ 2. The behavior is similar for all cases: there is region with Q > q to the left of the plot, and another with Q < q to the right.These regions correspond to the gravitationally bound and unbounded solutions respectively, as mentioned above.On the other hand, the left panel of the Figure shows families of boson stars with a super-critical charge q > 1/ √ 2, including the sub-critical case q = 0.7.As we can expect, the behavior for the case with q = 0.7 is still the same as for the other sub-critical cases on the right panel, but the region with Q > q corresponding to gravitationally bound solutions is now almost negligible.For the critical case q = 1/ √ 2 the region with Q > q has now vanished completely and all solutions are unbound.Furthermore, we find that the limit determined by the critical 2 as a function of the central field amplitude φ0(0).Right panel: Subcritical cases with q < 1/ √ 2. Left panel: Super-critical cases with q ≥ 1/ √ 2 .We have included the sub-critical case q = 0.7 in the left panel to show that the critical charge determines a upper bound of Q/M for all solutions.
charge qc = 1/ √ 2 (black line) acts as an upper bound for all solutions, that is we find that Q ≤ 1/ √ 2 for all the families of charged boson stars.From the definition of Q we can in fact see that this upper bound corresponds to having precisely Q = M .This means that even though we have found super-critical solutions in the sense that q > m (q > 1/ √ 2), there are no super-critical solutions in the sense that , which is precisely what one would have expected on physical grounds.That is, in general relativity it is not so much the ratio between the charge and mass parameters q and m the one that determines the existence of solutions, but rather the ratio between the total charge Q and the total mass M .

B. Some particular configurations
We will now show examples of some particular configurations for charged boson stars in order to better understand the behavior of the different functions in terms of the radial coordinate r. Figure 7 shows solutions for four boson stars with charge q = 0.4, and different central values of the scalar field φ 0 (0) = 0.025, 0.05, 0.1, 0.15.The top panel shows the scalar field φ 0 (r) as a function of the radial coordinate r.Notice that the scalar field is plotted in a logarithmic scale in order to show that the field does decay exponentially far away, which in a logarithmic scale corresponds to a straight line.The shape of the field does not change significantly for the different central values φ 0 (0), but the solution becomes more compact and cuspy for larger values of φ 0 (0).Notice also that for φ 0 (0) = 0.025 the exponential decay only becomes apparent for very large values of r, as already mentioned above.The middle panel of the Figure shows the lapse function α(r), and the lower panel shows the radial metric function A(r).From the Figure we can clearly see that the lapse and radial metric satisfy the boundary conditions at infinity, α(r) → 1 and A(r) → 1, and also A(r = 0) = 1.Notice that as φ 0 (0) increases the lapse becomes smaller at the origin, while the radial function has a larger and larger maximum that also moves closer to the origin.
Figure 8 shows solutions for configurations with the critical charge q = 1/ √ 2, and different central values of the scalar field φ 0 (0) = 0.025, 0.05, 0.1, 0.15.Although the overall shape is very similar to that of the previous Figure, we can notice that the scalar field distribution is in general wider than in the solutions with q = 0.4, i.e. the effective radius is larger (notice the change of scale on the horizontal axis).This implies that for smaller values of the central amplitude we in fact need to integrate to much larger values of r in order to see the correct exponential decay.The increase in the width of the solution is most probably due to the electromagnetic repulsion.
As our last example, in Figure 9 we show solutions for configurations with a super-critical charge q = 0.735.The behavior of the different functions is again similar to the previous cases, though we can see that the scalar field distribution is now even wider.From our previous discussion we know that super-critical solutions are only allowed for a small range of values of the central amplitude, which in this case corresponds to 0.22 φ 0 (0) 0.328.The plot then shows the cases with φ 0 (0) = 0.225, 0.2625, 0.3, 0.327.We can see that in all these cases the scalar field does indeed decay exponentially for large r.

C. Time evolutions
The numerical solutions for charged boson stars can be taken as initial data for dynamical evolutions.We have performed short evolutions of the different configurations described in the previous section using the OllinSphere code, which is a fully non-linear time evolution code for numerical relativity in spherical symmetry previously described in [23,25,26].We evolve the unperturbed data in order to verify that the different frequencies obtained in our analysis do correspond to the frequencies observed during dynamical evolution (we will leave a detailed study of the evolution of perturbed solutions for a future work).The numerical evolution code integrates the Einstein equations in time, using the Baumgarte-Shapiro-Shibata-Nakamura formulation adapted to spherical symmetry [27][28][29][30], coupled with the Klein-Gordon equation for the complex scalar field, and the Maxwell equations for the electromagnetic field as described in [23].
Figure 10 shows the time evolution of the (unperturbed) real part of the scalar field evaluated at the origin, Re(φ(t, r = 0)), for two of the boson star families described in the previous section.In particular, we consider the families with charges q = 0.4 and q = 1/ √ 2, corresponding to Figures 7 and 8.The evolutions presented here were performed with a resolutions of ∆r = 0.005 and a Courant parameter ∆t/∆r = 0.5.
Similarly, Figure 11 shows the time evolution of the (unperturbed) real part of the scalar field evaluated at the origin for the boson star models corresponding to Figure 9, which are solutions with a super-critical charge q = 0.735.As discussed above, these super-critical solutions are all gravitationally unbound and are therefore expected to be unstable.Indeed, in our numerical simulations we have found that even a small discretization error is sufficient to trigger either the collapse of these solutions to a black hole, or else their dispersion to infinity, in a relatively short time (we will discuss this in detail in a future publication).Because of this, in order to reduce the discretization error, for these evolutions we have increased the resolution so that we now take ∆r = 0.001 (with the same Courant parameter as before).Also, in the Figure we only plot the time evolutions up to a time t ∼ 32 (whereas the time evolutions of Figure 10 are plotted to much larger times).  .Time evolution of the real part the scalar field evaluated at the origin, Re(φ(t, r = 0)), for the different boson star models shown in Figure 9 with supercritical charge q = 0.735.The inset shows the frequencies obtained by performing a Fourier transform of these time evolution data.

D. Summary of properties of our particular configurations
Finally, in Table I we show the frequencies ω obtained for each of these specific configurations discussed above (after the rescaling and gauge transformation described in Section V), as well as the total charge Q, total mass M , total particle number N , effective radius R 99 , and binding energy E B .Do notice that the table shows the total charge Q = qN = √ 2 qN , and not the rescaled charge Q, this is in order to make it easier to compare Q directly with the  I. Charged boson stars models for different values of the charge q and different central amplitudes φ0(0).We show the frequency obtained from the solution of the eigenvalue problem ω (after the rescaling and gauge transformation described in the text), as well as the total charge Q, total mass M , total particle number N , effective radius R99, and binding energy EB.We also show the frequencies obtained from a Fourier transform of the dynamical evolution of the scalar field ωF .total mass M .From the table one can see that in all cases we have Q < M , with the particular configuration with critical charge q = 1/ √ 2 and small amplitude φ 0 (0) = 0.025 having Q very close to M , though still smaller.The table also includes the frequencies of oscillation ω F obtained from a Fourier transform of the time evolution data from Figures 10 and 11 (shown in the insets of those Figures).One can see that the frequencies obtained from the solution of the eigenvalue problem do indeed correspond, to several decimal places, with the frequencies obtained from the Fourier transform of the time evolution.
As we have argued above, the disappearance of the local maximum of the total mass M for super-critical configurations would seem to indicate that all such solutions are unstable.This conclusion is strengthened by considering the binding energy E B for the different configurations.Indeed, we find that for boson stars with a sub-critical charge there always exists a region where the solutions are gravitationally bound such that their binding energy is negative, E B < 0. On the other hand, all solutions with super-critical charges have a positive binding energy E B > 0, so they are gravitationally unbound.In relation to this, one can show that the binding energy is directly related to the total charge to mass ratio Q/M as E B = M [1 − (m/q)(Q/M )] = M (1 − Q/q), where Q is now defined as Q = (Q/M )/ √ 2. From this we find that all super-critical solutions are such that Q < q.
But from our results we in fact find the much stronger conclusion that all charged boson star configurations, both sub-critical and super-critical, are such that Q ≤ M ( Q ≤ 1/ √ 2), with the equality only achieved for the specific case of a critical charge q = 1/ √ 2 when φ 0 (0) → 0. This implies that, even though we do find a family of super-critical solutions in the sense that q > m (q > 1/ √ 2), there are in fact no super-critical solutions in the sense that Q > M .In other words, in general relativity it is not the ratio between the charge and mass parameters q and m the one that determines the existence of solutions, but rather the ratio between the total charge Q and the total mass M , as one would have expected on physical grounds.
Finally, we have also performed some preliminary time evolutions for unperturbed sub-critical and super-critical configurations in order to verify that the different frequencies obtained from our solution to the eigenvalue problem in fact correspond to the frequencies observed during a dynamical evolution, and find that the frequencies do coincide to several decimal places, indicating that our solution to the eigenvalue problem is correct.Our numerical evolutions also indicate that super-critical configurations are indeed unstable, as even a (small) numerical discretization error is enough to trigger either collapse to a black hole, or dispersion to infinity (though we will report on this in detail elsewhere).

Figure 2 .
Figure 2. Total mass (in units of M 2P L /m) as a function of the central value of the field φ0(0) for families of boson star solutions with different values of the charge q (notice that the vertical axis is on a logarithmic scale).

Figure 4 .
Figure 4. Compactness C99 := M/R99 as a function of the central value of the scalar field φ0(0).The compactness for boson stars with critical charge behaves as an upper bound for all other configurations.

Figure 7 .Figure 8 .
Figure 7. Charged boson stars solutions corresponding to configurations with charge q = 0.4 and different central amplitudes φ0(0) = 0.025, 0.05, 0.1, 0.15.The top, middle and lower panels show respectively the scalar field φ0(r) (in a logarithmic scale), the lapse function α, and the radial metric A(r).The insets show a closeup or the behavior of the different functions close to the origin.

Figure 9 .
Figure 9. Charged boson stars solutions corresponding to configurations with a super-crticial charge q = 0.735, and different central amplitudes φ0(0) = 0.225, 0.2625, 0.3, 0.327.The top, middle and lower panels show respectively the scalar field φ0(r) (in a logarithmic scale), the lapse function α, and the radial metric A(r).

Figure 10 .
Figure 10.Time evolution of the real part of the scalar field evaluated at the origin, Re(φ(t, r = 0)), for the different boson star models described in the previous Section.Top panel: Models with charge q = 0.4 corresponding to Figure7.Bottom panel: Models with charge q = 1/ √ 2 corresponding to Figure8.The inset shows the frequencies obtained by performing a Fourier transform of these time evolution data.

Figure 11
Figure 11.Time evolution of the real part the scalar field evaluated at the origin, Re(φ(t, r = 0)), for the different boson star models shown in Figure9with supercritical charge q = 0.735.The inset shows the frequencies obtained by performing a Fourier transform of these time evolution data.