Phase Diagrams of Charged Compact Boson Stars

TL;DR

This paper explores the complex phase structure of charged compact boson stars and shells in a gravity-coupled scalar field theory, revealing rich bifurcation phenomena through extended 2D and 3D phase diagrams.

Contribution

It extends previous phase diagram studies of charged boson stars to larger variable ranges, including 3D diagrams and detailed bifurcation analysis.

Findings

Identification of multiple bifurcation points in phase diagrams.

Discovery of an infinite series of bifurcation points.

Rich physics in the phase structure of charged boson stars.

Abstract

Compact boson stars, whose scalar field vanishes identically in the exterior region, arise in a theory involving a {\it massless} complex scalar field with a conical potential, when coupled to gravity. Their charged compact generalizations, obtained in the presence of a U(1) gauge field, exhibit further interesting features. On the one hand, charged compact boson shells can arise, whose scalar field vanishes also in the central region, while on the other hand, the domain of existence of charged compact boson stars exhibits bifurcation points. First 2D phase diagrams have been studied before. Here we extend these earlier studies to a larger range of the variables and study additional phase diagrams. We then extend these studies to obtain 3D phase diagrams and present these with a detailed discussion of their various regions with respect to the bifurcation points and argue, that there is an infinite series of such bifurcation points. Thus the theory is seen to contain rich physics in a particular domain of the phase diagrams. We also discuss the dependence of the fields on the dimensionless radial coordinate for some representative points of the phase trajectories in the phase diagrams of the theory.

Equations

$$S$$

$${\mathcal{L}}_M$$

$$D_{\mu }\Phi$$

$$F_{\mu \nu }$$

$$V(|\Phi |)$$

$$G_{\mu \nu }$$

$${\partial }_{\mu }\left(\sqrt{-g}{F}^{\mu \nu }\right)$$

$$D_{\mu }\left(\sqrt{-g}{D}^{\mu }\Phi \right)$$

$${\left[D_{\mu }\left(\sqrt{-g}{D}^{\mu }\Phi \right)\right]}^{\ast }$$

$$T_{\mu \nu }$$

$$ds^{2}=\biggl{[}-A^{2}Ndt^{2}+N^{-1}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta d\phi^{2})\biggr{]}\,,$$

$$G_t^t$$

$$G_r^r$$

$$G_{\theta }^{\theta }$$

$$\Phi (x^{\mu })=\varphi (r)e^{i\omega t},A_{\mu }(x^{\mu })d{x}^{\mu }=A_t(r)dt.$$

$$G_t^t=8\pi G{T}_t^t,$$

$$\displaystyle G_{\theta}^{\theta}=8\pi G\ T_{\theta}^{\theta}\ ,\ \$$

$$\frac{-1}{r^2}{\left[r\left(1-N\right)\right]}^{\prime }$$

$$\frac{2r{A}^{\prime }N-A{\left[r\left(1-N\right)\right]}^{\prime }}{A{r}^2}$$

$$\frac{2r{\left[r{A}^{\prime }N\right]}^{\prime }+{\left[A{r}^2{N}^{\prime }\right]}^{\prime }}{2A{r}^2}$$

$$\beta =\frac{\lambda e}{\sqrt{2}},{\alpha }^2(:=a)=\frac{4\pi G{\beta }^{2/3}}{e^2},$$

$$h(r)=\frac{(\sqrt{2}e\varphi (r))}{{\beta }^{1/3}},b(r)=\frac{(\omega +e{A}_t(r))}{{\beta }^{1/3}}.$$

$$h(\hat{r})=\frac{(\sqrt{2}e\varphi (\hat{r}))}{{\beta }^{1/3}},b(\hat{r})=\frac{(\omega +e{A}_t(\hat{r}))}{{\beta }^{1/3}}.$$

$${\left[AN{\hat{r}}^2{h}^{\prime }\right]}^{\prime }$$

$${\left[\frac{{\hat{r}}^2{b}^{\prime }}{A}\right]}^{\prime }$$

$$N^{\prime }$$

$$A^{\prime }$$

$$h^{\prime \prime }$$

$$b^{\prime \prime }$$

$$\hat{r}=x{\hat{r}}_o,0\le x\le 1,$$

$$\begin{gathered} A({\hat{r}}_o)=1 \\ , \end{gathered}$$

$$N(0)=1,b^{\prime }(0)=0,$$

$$h^{\prime }(0)=0,h({\hat{r}}_o)=0,h^{\prime }({\hat{r}}_o)=0.$$

$$j^{\mu }=-ie\left[\Phi (D^{\mu }\Phi {)}^{\ast }-{\Phi }^{\ast }(D^{\mu }\Phi )\right],j_{;\mu }^{\mu }=0,$$

$$Q=-\frac{1}{4\pi }{\int }_0^{{\hat{r}}_o}{j}^t\sqrt{-g}drd\theta d\varphi ,j^t=-\frac{h^2(\hat{r})b(\hat{r})}{A^2(\hat{r})N(\hat{r})}.$$

$$M=\biggl{(}1-N(\hat{r}_{o})+\frac{\alpha^{2}Q^{2}}{\hat{r}_{o}^{2}}\biggr{)}\frac{\hat{r}_{o}}{2}.$$

$$a_{c_n}=a_{c_{\infty }}+f{a}^n$$

$$P$$

Full text

††journal: Eur. Phys. J. C

1 Introduction

Boson stars represent localized solutions of scalar fields coupled to gravity (see the reviews Lee:1991ax ; Jetzer:1991jr ; Mielke:2000mh ; Schunck:2003kk ; Liebling:2012fv ). In their original form, they are based on a complex scalar field with a mass term only Kaup:1968zz ; Feinblum:1968nwc ; Ruffini:1969qy , while later work has included self-interactions of the scalar field, allowing for much higher masses of the resulting boson stars Colpi:1986ye ; Friedberg:1986tq .

Since being first conceived half a century ago, boson stars have received tremendous interest. Much of this interest resides in their potential astrophysical applications. These range from being promising dark matter candidates all the way to representing highly compact supermassive objects in galactic centers Lee:1991ax ; Jetzer:1991jr ; Mielke:2000mh ; Schunck:2003kk ; Liebling:2012fv ). Among the many applications accretion disks around boson stars have been considered Meliani:2015zta ; Meliani:2017ktw or the gravitational-wave signatures of boson stars and binary systems of boson stars have been studied Macedo:2013jja ; Palenzuela:2017kcg . In addition, boson stars have found applications based on the AdS/CFT conjecture Hartmann:2012wa .

However, boson stars are also interesting from a more theoretical point of view. On the one hand, they allow a study of numerous features of more complicated systems, and thus represent a nice model to learn and analyze. On the other hand, they contain the freedom to allow for new fascinating properties, not encountered before – at least not under similar circumstances. One of the most impressive recent examples in this context represents the discovery of hairy black holes based on boson stars Herdeiro:2014goa ; Herdeiro:2015tia . Another recent surprise was the realization that rotating boson stars allow for static orbits Collodel:2017end .

When considering the properties of boson stars, first of all their mass, their particle number and their size, are of interest. The particle number of boson stars arises from the U(1) symmetry of the underlying theory, which gives rise to a conserved current and thus an associated charge. When the U(1) symmetry is made a local symmetry, charged boson stars arise Jetzer:1989av ; Jetzer:1989us ; Jetzer:1993nk ; Pugliese:2013gsa ; Delgado:2016jxq ; Collodel:2019ohy . While their particle number is proportional to their charge, their basic features seem to only straightforwardly generalize those of their uncharged brethren. This changes fundamentally, when a different type of scalar field potential is employed, namely a conical potential Arodz:2008jk ; Arodz:2008nm .

First of all, a conical potential allows for compact boson star solutions, where compact is not meant to describe the physical property of having a large mass residing in an area with a small radius, as typically employed in connection with neutron stars or black holes. Here compact is meant to describe a solution whose scalar field vanishes outside a compact region of space identically. Thus while ordinary boson stars do not feature sharp boundaries, since their scalar field exhibits only an exponential fall-off, in compact boson stars, the scalar field is completely confined to a finite region of space Kleihaus:2009kr ; Kleihaus:2010ep ; Hartmann:2012da ; Hartmann:2013kna ; Kumar:2014kna ; Kumar:2015sia ; Kumar:2016oop ; Kumar:2016sxx ; Kumar:2017zms .

In the presence of a gauge field, however, charged compact boson stars exhibit new exciting features. One of these corresponds to the fact, that charged compact boson shells arise Kleihaus:2009kr ; Kleihaus:2010ep . Here the boson field is finite only inside a shell and vanishes identically outside this shell. Interestingly, there does not even have to be a vacuum in the central region of the spacetime. Instead, a black hole may reside there, that will then be surrounded by a matter shell as discussed in detail before Kleihaus:2010ep . Another interesting feature of these charged compact boson stars represents the presence of bifurcations points in the domain of existence Kleihaus:2009kr ; Kumar:2017zms .

Studies of the charged compact boson stars in a theory involving a massive complex scalar field with a conical potential coupled to a U(1) gauge field and gravity were undertaken in Refs. Kumar:2014kna ; Kumar:2015sia and Kumar:2016oop ; Kumar:2016sxx . In this work, the 2D and 3D phase diagrams were investigated in rather great detail, and multiple bifurcation points showing the presence of rather rich physics were obtained for a particular range of the parameters of the theory. These studies were undertaken in the absence Kumar:2014kna ; Kumar:2015sia as well as in the presence Kumar:2016oop ; Kumar:2016sxx of a cosmological constant in the theory.

In this work, we will deepen our earlier studies Kleihaus:2009kr ; Kumar:2017zms of the domain of existence of the massless theory with a conical potential in terms of 2D phase diagrams. Firstly we consider a larger range of the variables involved as compared to our earlier studies, and secondly we study additional 2D phase diagrams of the theory. We then extend these studies to the 3D phase diagrams, and we present the corresponding 3D plots of the various phase diagrams with a detailed discussion of the various regions of the phase diagrams with respect to the bifurcation points. In particular, we argue that there is an infinite sequence of such bifurcation points, and present an approximate formula for this sequence. The theory is thus seen to contain rich physics in a particular domain of the phase diagrams.

Further, we also study the distribution of the field variables of the theory with respect to the dimensionless radial coordinate at some representative points of the phase trajectories in the phase diagrams of the theory. For this purpose, we choose some particular phase trajectories taking them as examples to illustrate the properties of the corresponding compact boson stars. However, the procedure could be used to study the distribution of fields at any given representative point of any given trajectory (to be fixed by fixing the free dimensionless parameter of the theory designated by $a$ at a later stage in the text).

In our studies presented here, we construct the boson star solutions of this theory numerically. Our numerical procedure is based on the Newton-Raphson scheme with an adaptive stepsize Runge-Kutta method of order 4, and our numerical techniques have been calibrated by reproducing the previous work of Refs. Kleihaus:2009kr ; Kleihaus:2010ep ; Kumar:2017zms and Kumar:2014kna ; Kumar:2015sia ; Kumar:2016oop ; Kumar:2016sxx .

The paper is organized as follows. In the next section, we consider the mathematical formalism of the theory and derive the equations of motion. In Sec. 3, we rescale the constant parameters of the theory and the field variables to make them dimensionless. We further rescale the dimension-full radial coordinate $r$ to the dimensionless radial coordinate $\hat{r}$ and then express the equations of motion in terms of the dimensionless field variables as well as their derivatives expressed in terms of the dimensionless arguments. This then leads to a single parameter theory. In Sec. 4, we consider the boundary conditions. In Sec. 5, we lay down our scheme for obtaining the numerical solutions. In Sec. 6 and 7, we study and discuss in detail the 2D and 3D phase diagrams of boson stars, respectively. In Sec. 8, we study the mass $~{}M$ and charge $~{}Q$ of these gravitating objects. In Sec. 9, we study the distribution of the fields with respect to the dimensionless radial coordinate of the theory for some representative points in the phase diagrams. The summary and conclusions are presented in Sec. 10.

2 Mathematical Formalism and Field Equations

The theory that we study is defined by the following action:

[TABLE]

Here $R$ is the Ricci curvature scalar, and $G$ is Newton’s gravitational constant, $g=det(g_{\mu\nu})$, where $g_{\mu\nu}$ is the metric tensor, the asterisk denotes complex conjugation, and $\lambda$ is a constant parameter. The variational principle then leads to the equations of motion:

[TABLE]

[TABLE]

where the energy-momentum tensor $T_{\mu\nu}$ is given by

[TABLE]

We work with the static spherically symmetric metric in Schwarzschild-like coordinates

[TABLE]

where the arguments of the metric functions $A(r)$ and $N(r)$ have been suppressed. The components of Einstein tensor ($G_{\mu\nu}$) are then obtained as follows,

[TABLE]

where the prime denotes differentiation with respect to $~{}r$.

Assuming a vanishing magnetic field, we make the following Ansätze for the matter fields,

[TABLE]

With these Ansätze the Einstein equations

[TABLE]

(with the arguments of the field variables being suppressed) reduce to:

{strip}

[TABLE]

Here, the primes denote differentiation with respect to $~{}r$. Also, the equation $G_{\varphi}^{\varphi}=8\pi G\ T_{\varphi}^{\varphi}$ leads to an equation identical to Eq. (10).

3 Field Equations in Rescaled Variables

We now introduce the new constants,

[TABLE]

where $~{}a$ is dimensionless, and redefine the functions $\phi(r)$ and $A_{t}(r)$ through the following relations:

[TABLE]

For the sake of completeness, we mention that the metric functions $A(r)$ and $N(r)$ are already dimensionless field variables.

We also introduce a dimensionless coordinate $\hat{r}$ defined by $\hat{r}:=\beta^{1/3}\,{r}$ (which in turn implies $\frac{d}{d{r}}=\beta^{1/3}\frac{d}{d\hat{r}}$). Eq. (12) then gives:

[TABLE]

From now onwards, we will change the arguments of the field variables from the dimension-full radial coordinate $r$ to the dimensionless radial coordinate $\hat{r}$ and the primes will denote differentiation with respect to the dimensionless radial coordinate $\hat{r}$.

The matter field equations of motion involving the dimensionless field variables as well as their derivatives (with ${\rm sign}(h)$ denoting the usual signature function) then read:

[TABLE]

We thus obtain the set of equations of motion to be solved numerically as:

[TABLE]

In order to solve the above equations of motion numerically, we further introduce a new coordinate $x$ defined through the following relation:

[TABLE]

implying that the center of the boson star i.e. $\hat{r}=0$ corresponds to $x=0$ and the outer boundary of the boson star, i.e. $\hat{r}=\hat{r}_{o}$ corresponds to $x=1$.

4 Boundary Conditions

We now need to determine the boundary conditions for the functions describing the metric and the fields. For the metric function $A(\hat{r})$ we choose

[TABLE]

where $\hat{r}_{o}$ denotes the outer radius of the star. Thus $A$ assumes its asympototic value already at the star’s boundary.

In order to obtain globally regular ball-like boson star solutions, we impose

[TABLE]

Since the star has a sharp boundary, with the electro-vac equations holding in the exterior region $\hat{r}>\hat{r}_{o}$, we match the Reissner-Nordström solution at $\hat{r}=\hat{r}_{o}$.

The U(1) invariance of the theory leads to a conserved Noether current,

[TABLE]

whose time component corresponds to the charge density. The global charge $Q$ of the boson star is then given by

[TABLE]

We like to clarify here that the frequency of the scalar field $\omega$ shows up explicity in the Eqs. (8)-(10) however, it does not show up explicitley in Eqs.(14)-(19), merely because we introduced new dimensionaless field variables (with dimensionless arguments) through Eqs. (12)-(13). and the definition of $b(\hat{r})$ contains $\omega$ explicitly. The charge density given by Eq. (4) is seen to be proportional to $b(\hat{r})$ which in turn has its definition (in terms of $\omega$) given by Eqs. (12)-(13). One could thus visualize the dependence of the charge or of the charge density in terms of $\omega$ if one likes.

The mass can be read off the metric, as usual, making use of the fact, that the metric outside the compact star corresponds to a Reissner-Nordström metric. In the units employed, we then obtain for the mass $M$ of the boson star solutions

[TABLE]

An extremal Reissner-Nordström metric would satisfy a proportionality between the mass $~{}M$  and the charge $~{}Q~{}$ given by $~{}M=\alpha Q=\sqrt{a}Q~{}$, where $~{}\alpha(:=\sqrt{a})~{}$ enters here because of the units employed (whereas in the usual geometric units the extremal solution satisfies $~{}M=Q~{}$). We can now consider three different cases for the exterior solution ($r>r_{o}$) outside the bosonic matter distribution: (i) The case $~{}M/Q<\sqrt{a}~{}$ corresponds to a horizonless (naked) Reissner-Nordström solution. (We note, that naked solutions would also arise, when the mass and charge of known elementary particles would be inserted on the left hand side of this relation.) (ii) The case $~{}M/Q=\sqrt{a}~{}$ corresponds to an extremal Reissner-Nordström solution. (We note, that the sets of boson star solutions end, when such an extremal Reissner-Nordström solution is reached and $b(0)=0$.) (iii) The case $~{}M/Q>\sqrt{a}~{}$ corresponds to a solution where the radius of the boson star is greater than its putative event horizon radius $~{}{\hat{r}_{H}}~{}$ (implying that the boson star is well outside the range of becoming a black hole.

This discussion is reflected in Figs. 9, which show the variation of the various fields with respect to $~{}\hat{r}~{}$. Here it is evident that the metric function $~{}N(\hat{r})~{}$ of the (selected) solutions always has a non-zero positive value, (and also does not become zero at the outer radius of the boson star) implying thereby that the boson stars either have a radius $~{}\hat{r}>{\hat{r}_{H}}~{}$, or that they are smoothly matched to a naked Reissner-Nordström solution for $~{}M/Q<\sqrt{a}~{}$. The case $~{}M/Q=\sqrt{a}~{}$ is not exhibited in Figs. 9. It would correspond to a solution where a throat develops and $b(0)=0$.

In the present work, we study (i) the variation of $M$ with $\hat{r}_{0}$, (ii) the variation of $Q$ with $\hat{r}_{0}$, (iii) the variation of $M$ with $Q$, and (iv) the variation of $M/Q$ with $Q$, where the dimensionaless variable $~{}a~{}$ varies in the range $a=0.050$ to $a=0.250$ (to be seen in detail later in the discussion of the results).

5 Numerical Solutions

We now consider the numerical solutions of the coupled nonlinear differential equations, Eqs. (16)-(19), subject to the above set of boundary conditions: Eqs. (21)-(22) (cf. Refs. Kleihaus:2009kr ; Kleihaus:2010ep ; Hartmann:2012da ; Hartmann:2013kna ; Kumar:2014kna ; Kumar:2015sia ; Kumar:2016oop ; Kumar:2016sxx ; Kumar:2017zms ).

Going into the details, the solution of the equations is known outside the outer radius of the boson stars, since there we have the electo-vac equations with the known Reissner-Nordström solution. In the inner region the scalar field is non-trivial and vanishes only on the boundary. In this inner region the equations are then solved using the method of shooting, subject to the specified boundary conditions, compatible with the exterior Reissner-Nordström solution.

The shooting method is based on converting the two-point boundary value problems into an equivalent initial value problem and this method uses the methods used in solving initial value problems. The solution is done by assuming initial values at one boundary for the variables that would have been given if the ordinary differential equation were a initial value problem. The boundary values obtained at the other boundary is then compared with the actual boundary values. Using some intuitive approach, one tries to get as close to the boundary value as possible. In general, one finds discrepancies from the desired boundary values at the other boundary. Now we have a multidimensional root finding problem and we use a root finding method to find the adjustment of the assumed initial values at the starting point that makes the discrepancies at the other boundary vanish.

To solve the coupled nonlinear differential equations, Eqs. (16)-(19), subject to the above set of boundary conditions, Eqs. (21)-(22), we require initial values for the dependent variables $A(0)$, $b(0)$ and $h(0)$ at the inner boundary i.e. at $x=0$. We have to guess initial values in a way such that the boundary conditions at the outer boundary i.e. at $x=1$ are satisfied.

For $A(0)$ we can choose an initial guess between [math] and $1$ because $A^{\prime}(r)$ has a positive slope given by eq. (18) and $A(r)$ has an upper bound given by $A(r_{o})=1$. For $h(0)$ and $b(0)$ we can choose any finite value as an initial guess. For the parameter $a=0$ we can for example, start with an initial guess of $A(0)=1$, $h(0)=0$ and $b(0)>0$ and the obtained solution is used as an initial guess for the next solution.

In the shooting method we choose values for all the dependent variables, consistent with other boundary conditions, at the inner boundary, i.e. at $x=0$, and then integrate the initial value problem by the Runge-Kutta method of order 4 with an adaptive step size, where the step size is controlled by the estimate of errors at each step, until we arrive at the outer boundary, i.e. at $x=1$.

At each iteration step we use the Newton-Raphson method to find the zeros of the function representing the errors by which the solutions to the initial value problem fails to satisfy the boundary condition at $x=1$. The set of coupled equations is solved iteratively, until the prescribed accuracy and thus the stopping criterion is reached. We have been using this method here, requiring an accuracy of 6-7 digits. For further information on the shooting method for boundary value problems see Press:1992zz .

Employing this method, we determine the domain of existence of the solutions of the above equations by varying the parameter $~{}a$, the single remaining parameter of the theory after the redefinitions and rescalings of Sec. III.

In the following sections, we will study the 2D and 3D phase diagrams, the mass $~{}M$ and charge $~{}Q$ of these gravitating objects, and subsequently we will study the distribution of the fields with respect to the dimensionless radial coordinate at some representative points on phase trajectories in the phase diagrams of the theory.

6 2D Phase Diagrams

Studies of the 2D phase diagrams of the theory considered in the present work were first carried out by Kleihaus, Kunz, Lämmerzahl and List in Kleihaus:2009kr ; Kleihaus:2010ep , where a bifurcation point in the relevant phase diagrams of the theory was obtained.These studies were extended in our recent work Kumar:2017zms , where additional bifurcation points in the relevant 2D phase diagrams theory were found.

In this section we focus mainly on the additional bifurcation points in the 2D phase diagrams of the theory. In particular, we not only review our recent work Kumar:2017zms , which investigates the 2D phase diagrams, but we now study them for larger ranges of the variables, as depicted in the various results, as compared to our earlier studies Kumar:2017zms . Secondly, we also study some additional phase diagrams of the theory, which were not studied in our earlier work Kumar:2017zms , and thirdly, we discuss the next new bifurcation point and present a formula, indicating the presence of an infinite sequence of such bifurcation points.

For this, we first study and discuss the phase diagrams for boson stars as shown in Figs. 1, 1 and 1. In Fig. 1, we consider the functions $h(0)$ and $b(0)$, in Fig. 1 the functions $A(0)$ and $b(0)$, and in Fig. 1 the functions $A(0)$ and $h(0)$ (where $h(0)$, $b(0)$ and $A(0)$ denote the values of these functions at the centre of the star). Also, the figures shown in the insets in Fig. 1 to 1 exhibit particularly interesting areas of these figures with better precision. The asterisks shown in Fig. 1, corresponding to $h(0)=0$, represent the transition points from the boson stars to boson shells.

Let us now address the phase diagrams based on the values of the functions at the origin of the boson star, comprising the vector field function $b(0)$ and the scalar field function $h(0)$, that are obtained by studying the solutions for a sequence of values of the parameter $a$. Here we observe very interesting phenomena near specific values of $a$, where the system is seen to have bifurcation points, $B_{1}\,,\,B_{2}$ and $B_{3}$. These bifurcations arise at the following values of $a$, $a_{c_{1}}\simeq 0.198926,~{}a_{c_{2}}\simeq 0.169311$ and $a_{c_{3}}\simeq 0.168308$, respectively (and a possibility about the existance of further bifurcation points can not be ruled out in principle, if one could somehow increase the numerical accuracy beyond the one already achived by us in the present work (however, practically it is a rather very challenging problem )).

Thus, clearly the theory possesses rich physics in the domain $a=0.22$ to $a\simeq+0.16\,$.

For a better understanding, let us now divide the phase diagram in the vicinity of the first bifurcation point $B_{1}$ into four regions, which we denote by IA, IB, IIA and IIB (as depicted in Fig. 1). The asterisks in Fig. 1 reside on the axis $b(0)$, corresponding to $h(0)=0$, and they mark the transition points from the boson stars to boson shells Kumar:2017zms .

First we note that the regions IA, IIA and IIB do not possess any further bifurcation points, while the region IB does contain several more bifurcation points, which indicates that it contains rich physics. Consequently, we divide the region IB further into the regions IB1, IB2 and IB3 in the vicinity of the first additional bifurcation point $B_{2}$.

The region IB3 again contains a further bifucation point, which we denote $B_{3}$. This suggests a further subdivision of the phase diagram in the vicinity of $B_{3}$ into the regions IB3a, IB3b and IB3c, as seen in the inset of Fig. 1. Interestingly, the region IB3b contains closed loops, and the phase diagram in this region resembles the one of the region IB2.

From these figures we can conclude, that as the value of $a$ changes from $a=0.25$ to $a=0$, a lot of new rich physics is encountered. When varying $a$ from $a=0.25$ to the critical value $a=a_{c_{1}}$, the solutions exist in the two distinct domains, IIA and IIB (as depicted in Fig. 1). But when we decrease $a$ below this critical value $a=a_{c_{1}}$, the solutions reside in the regions IA and IB. Let us mention here, that for values of $a$ larger than $a=0.25$, the physics remains qualitatively the same as for the value $a=0.25\,$.

When the value of $a$ is decreased from the first critical value $a=a_{c_{1}}$ to the second critical value $a=a_{c_{2}}$, the curves in the region IA of the phase diagram we notice that the region IA in the phase diagram show a continuous deformation, leading to the region IB, which exhibits its own rich physics as explained above.

Decreasing the value of $a$ further below $a_{c_{2}}$, we see on the one hand in the region IA a continuous deformation of the curves towards $a=0$. However, on the other hand, a third bifurcation point is encountered in the region IB, dividing this region into IB1, IB2 and IB3. Again, there is a continous deformation of the curves in the region IB1, while the region IB2 again contains curves forming closed loops. Subdividing the region IB3 into the regions IB3a, IB3b and IB3c, there is again a continuous deformation of the curves in region IB3a, whereas the region IB3b again exhibits closed loops.

The regularity of the above described phenomenon makes it tempting to conjecture, that in fact a whole sequence of further bifurcation points may be waiting to be discovered, which might display a self-similar pattern. To support this conjecture, we have considered the set of bifurcation points $a_{c_{1}}$, $a_{c_{2}}$ and $a_{c_{3}}$, and contemplated how to extend this sequence to the next bifurcation point $a_{c_{4}}$. Since the difference between the bifurcation points $\Delta_{n}=a_{c_{n}}-a_{c_{n+1}}$ decreases strongly, an exponential approach toward a limiting value $a_{c_{\infty}}$ is expected. Thus an ansatz for the $n$-th bifurcation point $a_{c_{n}}$ of the type

[TABLE]

seems promising. Using the first three bifurcation points we have thus made a prediction for the fourth one and, indeed, found the fourth bifurcation point $a_{c_{4}}=0.16827861$, stretching our numerical calculations and accuracy to their current limits.

Since the numerical calculations become increasingly challenging, as one progresses from bifurcation point to bifurcation point, we have not been able to continue our quest for the next bifurcation points further, to support this fascinating expectation with even more data. But the fact, that our prediction was borne out, does lend considerable support to our conjecture. Indeed, a constantly increasing numerical accuracy is required to map out the domain of existence of the solutions. For instance, to analyze the bifurcation point B3, the value of $a$ needed already a specification of 6 decimal digits. The analysis of bifurcation point B4 needed two more digits. The expontial convergence of the values of the bifurcation points towards their limiting value, clearly requires a rapidly increasing number of accucarate digits of the calculations. Consequently, the global accuracy of the numerical scheme employed represents currently the barrier to continue these investigations of potential further bifurcation points.

In Fig. 1 we exhibit the radius $\hat{r}_{o}$ of the star versus the value of the vector field function at the center of the star $b(0)$. As in the previous figures, the region around the first bifurcation point $B_{1}$ can be divided into the four regions IA, IB and IIA, IIB, and the asterisks represent the transition points from the boson stars to boson shells. However, the oscillations (with respect to $b(0)$) seen in Figs. 1 and 1 in the regions IIA and IB, have transmuted in Fig. 1 into a spiralling behavior. The insets of Fig. 1 magnify parts of the physically interesting region IB.

7 3D Phase Diagrams

In this section, we present a study of the 3D phase diagrams of the charged compact boson stars. Figs. 3 to 3 show the 3D diagrams (with different viewing angles) of the scalar field function $h$ and the gauge field function $b$ versus the metric field function $A$ at the centre of the boson star for a large set of values of the constant $a$. The viewing angles are denoted by $(\theta,\phi)$, where we use the convention that $\theta$ denotes the angle of rotation about the $ox$-axis in the anticlockwise direction and takes values between $~{}0$ to $\pi$, and that $\phi$ is the angle of rotation about the $oz^{\prime}$-axis also in the anticlockwise direction and takes values between $~{}0$ to $2\pi$. Figs. (a) to (f) correspond, respectively, to the values $(\theta,\phi)\equiv(60,60),\ (60,175),\ (50,245),\ (135,35),\ (130,70)$ and $(135,340)$.

Figs. 4 to 4 show the 3D diagrams of the scalar field function $h(0)$ and the gauge field function $b(0)$ versus the radius of the star $\hat{r}_{o}$ for a large set of values of the constant $~{}a$. Figs. (a) to (f) correspond, respectively, to the values of the viewing angles $(\theta,\phi)\equiv(60,100),\ (115,15),\ (60,200),$ $\ (40,240),\ (140,60)$ and $(130,40)$.

8 Mass $M$ and Charge $Q$

Let us next consider the global properties of the solutions, their mass $M$ and their charge $Q$. We present the mass $M$ versus the radius $\hat{r}_{o}$ in Fig. 5, and the corresponding magnified region of the bifurcations is shown in Fig. 5. The charge $Q$ versus the radius $\hat{r}_{o}$ is depicted in Figs. 6 and 6, and it is clearly seen to possess a very similar dependence as the mass.

The stability of the boson stars can be investigated by considering the mass $M$ versus the charge $Q$, as shown in Fig. 7. Additionally it is instructive to consider the mass per unit charge $M/Q$ versus the charge, as shown in Fig. 8. In Fig. 7 we see that the curves of $M$ versus $Q$, that are located in region IA and concern the smaller values of $a$, all exhibit a monotonic increase from $M=Q=0$ towards their transition points, where the boson shells are formed, and which are marked by crosses. The boson stars forming these curves represent the fundamental boson star solutions (for the respective value of $a$). Therefore these boson stars should be stable and, in fact, the boson stars associated with all curves in region IA should be stable, since these solutions possess the lowest mass for a given charge (and parameter $~{}a$).

However, it is seen that above a certain value of $~{}a$, these curves no longer reach a boson shell, but instead their upper endpoint represents a solution, where a throat is formed. In that case the exterior space-time $~{}r>r_{0}$ corresponds to the exterior space-time of an extremal RN black hole. As was discussed in great detail in Kleihaus:2009kr ; Kleihaus:2010ep , this happens whenever the value $b(0)=0$ is encountered. Interestingly, the curves residing in region IIB encounter solutions with throats at both their endpoints, since both times the value $b(0)=0$ is reached. Because these solutions also possess the lowest mass for a given charge, they are also expected to be stable.

In contrast, in region IIA the solutions exhibit the typical oscillating/spiralling behavior, that is known to arise also for non-compact boson stars. Considered in a mass versus charge diagram, this behavior translates into a sequence of spikes, as seen in the insets of Figs. 7 and 8. In this case the solutions should only be stable along their fundamental branch, which ends when a first spike is encountered at a maximal value of the mass and the charge. Whenever a new spike is encountered, a new unstable mode should arise, analogously to the case of non-compact boson stars.

Addressing now the solutions in the bifurcation regions of the phase diagram, let us consider the region IB, and starting with the boson stars constituting the limiting curves. Adopting the value $a_{c_{1}}$, the two branches of solutions, which limit the region IA, possess lower masses than the two branches of solutions, that limit the region IB. Therefore the lower mass solutions are expected to be more stable.

Indeed, the two branches of solutions, which limit the region IB, might be classically stable, as well, at least until the first extrema of the mass and the charge are reached. From a quantum mechanical point of view, however, they would still be unstable, since the phenomenon of tunnelling might arise. Clearly, beyond these extrema, unstable modes should be present, and the solutions should also be classically unstable.

One can also extend these arguments to all the solutions in region IB. Quantum mechanically they should be unstable, since solutions with lower mass but the same value of the charge reside in the region IA (for all of them). But from a purely classical point of view the solutions with the lowest mass for a given $a$ within the region IB could still be stable. In contrast, the higher mass solutions should possess unstable modes and thus be classically unstable.

Figs. 7 and 7 depict the mass $M$ versus the charge $Q$ for the respective ranges of $a$ shown (where in Fig. 7 a part of Fig. 7 has been magnified). Similarly, Figs. 8 and 8 depict the mass to charge ratio $M/Q$ versus the charge $Q$ for the respective ranges of $a$ (where again in Fig. 8 a part of Fig. 8 has been magnified). These figures, in fact, illustrate that the solutions in the bifurcation region indeed correspond to higher mass solutions.

9 Distribution of Fields

Let us now consider the distribution of the fields of a representative set of solutions in order to understand how the physical configurations change, as we move along the phase trajectories in the phase digagrams. For that purpose we consider representative points $~{}P(\hat{r}_{o},h(0),b(0),A(0))$, located in the various regions of the phase diagrams. In particular, we choose the following set corresponding to two different sample values of $~{}a$, namely for (i) $~{}a=0.169$ and for (ii) $~{}a=0.25$.

For the case (i) we consider four particular representative points denoted by $P_{1},P_{2},P_{3},P_{4}$, defined as follows:

($P_{1}$): lying in region IA:

[TABLE]

($P_{2}$): lying in region IA (corresponding to the transition point to shells)

[TABLE]

($P_{3}$): lying in region IB2:

[TABLE]

($P_{4}$): lying in region IB3:

[TABLE]

The results for these four representative points are shown in Figs.9-9.

For the case (ii) we consider two representative points denoted by $P_{5}$ and $P_{6}$, defined as follows:

($P_{5}$): lying in region IIA:

[TABLE]

($P_{6}$): lying in region IIB:

[TABLE]

The results for these are shown in Figs. 9-9.

Figs. 9-9 describe the distribution of the four field functions $A(\hat{r})$, $h(\hat{r})$, $b(\hat{r})$ and $N(\hat{r})$) with respect to the dimensionless radial coordinate $~{}\hat{r}$, which are the solutions of the set of four coupled nonlinear differential Eqs. (16)-(19). Thus these figures show the values of these fields from the center of the star all the way not only up to the boundary of the compact boson star (where the scalar field $h(\hat{r})$ vanishes), but also beyond the boundary of the star into the exterior region. The electromagnetic field represented by $b(\hat{r})$ and the gravitational field represented by the metric functions $A(\hat{r})$ and $N(\hat{r})$ correspond to long range forces.

Fig. 9 represents a typical compact boson star configuration, located in the phase diagram in region IA. The matter distribution function $h(\hat{r})$ is maximal at the center and decreases monotonically towards the outer boundary $\hat{r}_{o}$ of the star, where it vanishes. The metric function $N(\hat{r})$ first decreases from its central value $N(0)=1$, necessary for regularity, then it reaches a minimum and increases again towards its asymptotic value of one, necessary for asymptotic flatness. The metric function $A(\hat{r})$ starts from a finite value at the center and increases monotonically towards the boundary of the star, reaching its asymptotic value of one already at the boundary. The gauge field function $b(\hat{r})$, finally, also exhibits a monotonic increase from the center of the star towards a finite value at infinity. Recall, that the spacetime exterior to the star corresponds to a part of a Reissner-Nordström spacetime.

As we move along the phase trajectory, keeping $a$ fixed, and changing $h(0)$ smoothly, first increasing and the decreasing it, the configurations change smoothly, until $h(0)=0$ is reached. This configuration is shown in Fig. 9 and also located in the phase diagram in region IA. It corresponds to the transition point to boson shells. Therefore the matter has been depleted in the interior of the star, so that the matter function vanishes at the center, and then increases slowly towards the outer regions of the star, where it reaches a maximum. From this maximum the matter density then rapdidly drops to zero at the boundary of the star. At the same time, the metric function $N(\hat{r})$ has developed a deep minimum ($N(\hat{r}_{\rm min})=0.00811$) close to the outer boundary $\hat{r}_{o}$ of the star ($\hat{r}_{\rm min}=3.13054<\hat{r}_{o}=3.14508$), while the metric function $~{}A(\hat{r})$ shows a steep rise close to the boundary. The gauge field function $b(\hat{r})$ is almost constant inside the star.

The phase trajectory for this value of $a$ consists of several disconnected parts, associated with several of the distinct regions of the phase diagram. Fig. 9 corresponds to a point in region IB2, i.e., in the region of the second bifucation point. Here the matter field $h(\hat{r})$ is concentrated much stronger at the center of the star, as compared to a typical star (as seen, e.g., in Fig. 9). The matter density then drops monotonically and rapidly towards the boundary of the star. The center remains regular, as seen from the behavior of the metric functions.

Fig. 9 illustrates the configurations in the region IB3 of the phase diagram, associated with the third bifurcation point. Now the matter is even more strongly localized at the center, dropping off even more steeply towards the boundary of the star. The strong localization of the matter at the center is reflected in the behavior of the metric functions, where the function $N(\hat{r})$ now also exhibits a strong decrease at the center (see inset to verify the regularity condition at the center), while the function $A(\hat{r})$ tends towards zero at the center (but remains finite with $A(0)=0.000677104$).

We conjecture, that for the further bifurcations points, the matter will be even more strongly localized at the center, and the metric function $N(\hat{r})$ will decrease even more strongly as well, while the metric function $A(\hat{r})$ will tend even closer towards zero at the center (but remain finite).

We now turn to another value of $a$, and thus another phase trajectory, which also has disconneted parts. Fig. 9 belongs to a boson star located in region IIA in the phase space diagram. This boson star looks very much like the one in Fig. 9. The maximum at the center has not yet increased much. However, it will increase strongly further along the phase trajectory, as soon as the oscillations will be encountered (which start at $h(0)=1.71$, $b(0)=1.9349$).

In Fig. 9, finally, a boson star belonging to region IIB is illustrated. This boson star is located on the rather short second part of the phase trajectory for this value of $a$. At both endpoints of this part of the trajectory a special solution is encountered, which is joining to an exterior extremal Reissner-Nordström solution, where $b(0)=0$ and $N(\hat{r}_{o})=0$. The boson star in this figure shows already clear signs of being close to such a special configuration. The metric function $N(\hat{r})$ has already a deep minimum ($N(\hat{r}_{\rm min})=0.00603$) close to the boundary of the star ($\hat{r}_{\rm min}=1.30994<\hat{r}_{o}=3.33932$), and the gauge field function $b(\hat{r})$ is already close to zero at the center of the star.

10 Summary and Conclusions

In this work, we have presented our studies of charged compact boson stars in a theory involving a massless complex scalar field with a conical potential coupled to a U(1) gauge field and gravity. Our presented work focuses mainly on the studies of the phase diagrams of this theory.

A study of the 2D phase diagrams of this theory was initially undertaken by Kleihaus, Kunz, Lämmerzahl and List in Refs. Kleihaus:2009kr ; Kleihaus:2010ep , where a bifurcation point in the relevant phase diagrams of the theory was obtained. These studies were extended recently in Ref. Kumar:2017zms , where some additional bifurcation points were obtained in the 2D phase diagrams Kumar:2017zms .

In this work, we have further extended our earlier studies of the 2D phase diagrams of the theory in two respects. Firstly we have considered a larger range of the variables involved as compared to our earlier studies, and secondly we have studied additional 2D phase diagrams.

In addition, we have further extended our above studies to present 3D phase diagrams of the charged compact boson stars. These 3D phase diagrams have been investigated, and a detailed discussion of the various regions of the phase diagrams with respect to the bifurcation points has been presented. The theory is seen to contain rich physics in a particular domain of the phase diagrams, namely in the domain $a=0.25$ to $a\simeq 0.16$, where we have identified four bifurcation points $B_{1},~{}B_{2},~{}B_{3}$ and $B_{4}$. In particular, we have predicted the fourth bifurcation point by analyzing the first three bifucation points and supposing and exponential law for the $n$-th bifurcation point, and then verifying its predicted value by direct computation. The existence of such a whole sequence of further bifurcation points following the exponential law seems therefore a well-supported conjecture, but it will remain numerically challenging to construct even further bifurcation points.

Physical properties of the solutions, including their mass, charge and radius have been investigated. The mass versus charge, and the mass per unit charge versus the charge, have been considered and the arguments concerning the stability of the solutions have been presented.

For any value of $a$, a fundamental branch of compact boson star solutions exists, which is expected to be stable, since the solutions on this branch represent the solutions with the lowest mass for a given charge. In that sense these solutions can be considered to represent the ground state. On the other hand, in the region of the bifurcations additional branches of solutions are present, which possess higher masses for a given charge. These solutions thus correspond to the excited states of the system. The lowest of these might be classically stable and unstable only quantum mechanically. To definitely answer this question, a mode stability analysis needs to be performed. However, such an analysis is beyond the scope of our present work and represents an interesting topic for future investigations.

In addition to this, we have also studied the distribution of field functions $A(\hat{r})$, $h(\hat{r})$, $b(\hat{r})$ and $N(\hat{r})$) versus the dimensionless radial coordinate $\hat{r}$ at some chosen representative points of some specific phase trajectories in the phase diagrams, in order to better understand the physical properties of the boson star configurations in the various parts of the phase diagram. This procedure can, in principle, be used to investigate the distribution of fields at any other representative points on any of the phase trajectories in the phase diagrams for any fixed value of $~{}a$.

Acknowledgements.

This work was supported in part by the DFG Research Training Group 1620 Models of Gravity as well as by and the COST Action CA16104 GWverse. SK would like to thank the CSIR, New Delhi, for the award of a Research Associateship during which a part of this work was carried out.