Dynamical boson stars

The idea of stable, localized bundles of energy has strong appeal as a model for particles. In the 1950s, John Wheeler envisioned such bundles as smooth configurations of electromagnetic energy that he called geons, but none were found. Instead, particle-like solutions were found in the late 1960s with the addition of a scalar field, and these were given the name boson stars. Since then, boson stars find use in a wide variety of models as sources of dark matter, as black hole mimickers, in simple models of binary systems, and as a tool in finding black holes in higher dimensions with only a single Killing vector. We discuss important varieties of boson stars, their dynamic properties, and some of their uses, concentrating on recent efforts.


The nature of a boson star
Boson stars are constructed with a complex scalar field coupled to gravity (as described in Section 2). A complex scalar field φ(t, r) can be decomposed into two real scalar fields φ R and φ I mapping every spacetime event to the complex plane φ(t, r) ≡ φ R (t, r) + iφ I (t, r). (1) Such a field possesses energy because of its spatial gradients and time derivatives and this energy gravitates holding the star together. Less clear is what supports the star against the force of gravity. Its constituent scalar field obeys a Klein-Gordon wave equation which tends to disperse fields. This is the same dispersion which underlies the Heisenberg Uncertainty Principle. Indeed, Kaup's original work [128] found energy eigenstates for a semi-classical, complex scalar field, discovering that gravitational collapse was not inevitable. Ruffini and Bonazzola [189] followed up on this work by quantizing a real scalar field representing some number of bosons and they found the same field equations.
None of this guarantees that such solutions balancing dispersion against gravitational attraction exist. In fact, a widely known theorem, Derrick's Theorem [69] (see also [187]), uses a clever scaling argument to show that no regular, static, nontopological localized scalar field solutions are stable in three (spatial) dimensional flat space. This constraint is avoided by adopting a harmonic ansatz for the complex scalar field φ(r, t) = φ 0 (r)e iωt (2) and by working with gravity. Although the field is no longer static, as shown in Section 2 the spacetime remains static. The star itself is a stationary, soliton-like solution as demonstrated in Figure 1.
There are, of course, many other soliton and soliton-like solutions in three dimensions finding a variety of ways to evade Derrick's Theorem. For example, the field-theory monopole of 't Hooft and Polyakov is a localized solution of a properly gauged triplet scalar field. Such a solution is a topological soliton because the monopole possesses false vacuum energy which is topologically trapped. The monopole is one among a number of different topological defects that requires an infinite amount of energy to "unwind" the potential energy trapped within (see [219] for a general introduction to defects and the introduction of [190] for a discussion of relevant classical field theory concepts).
In Section 2, we present the underlying equations and mathematical solutions, but here we are concerned with the physical nature of these boson stars. When searching for an actual boson star, we look not for a quantized wave function, or even a semiclassical one. Instead, we search for a fundamental scalar, say the long-sought Higgs boson. The Large Hadron Collider (LHC) hopes to determine the existence and nature of the Higgs, with evidence at the time of writing suggesting a Higgs boson with mass ≈ 125 GeV/c 2 [185]. If the Higgs does not ultimately appear, there are other candidates such as an axion particle. Boson stars are then either a collection of stable fundamental bosonic particles bound by gravity, or else a collection of unstable particles that, with the gravitational binding, have an inverse process efficient enough to reach an equilibrium. They can thus be considered a Bose-Einstein condensate (BEC), although boson stars can also exist in an excited state as well.
Indeed, applying the uncertainty principle to a boson star by assuming it to be a macroscopic quantum state results in an excellent estimate for the maximum mass of a BS. One begins with the Heisenberg uncertainty principle of quantum mechanics ∆p ∆x ≥ and assumes the BS is confined within some radius ∆x = R with a maximum momentum of ∆p = mc where m is the mass of the constituent particle mcR ≥ .
This inequality is consistent with the star being described by a Compton wavelength of λ C = h/ (mc). We look for the maximum possible mass M max for the boson star which will saturate the uncertainty bound and drive the radius of the star towards its Schwarzschild radius R S ≡ 2GM/c 2 . Substituting yields 2GmM max c = , which gives an expression for the maximum mass Recognizing the Planck mass M Planck ≡ c/G, we obtain the estimate of M max = 0.5M 2 Planck /m. This simple estimate indicates that the maximum mass of the BS is inversely related to the mass of the constituent scalar field. We will see below in Section 2 that this inverse relationship continues to hold with the explicit solution of the differential equations for a simple mass term in the potential, but can vary with the addition of self-interaction terms. Indeed depending on the strength of the coupling m and the other parameters of the self-interaction potential, the size and mass of the boson stars can vary from atomic to astrophysical scales.
Despite their connection to fundamental physics, one can also view boson stars in analogy with models of neutron stars. In particular, as we discuss in the following sections, both types of stars demonstrate somewhat similar mass versus radius curves for their solutions with a transition in stability at local maxima of the mass. There is also a correspondence between (massless) scalar fields and a stiff, perfect fluid (see Section 2.1 and Appendix A of Ref. [36]), but the correspondence does not mean that the two are equivalent [79]. More than just an analogy, boson stars can serve as a very useful model of a compact star, having certain advantages over a fluid neutron star model: (i) the equations governing its dynamics avoid developing discontinuities, in particular there is no sharp stellar surface, (ii) there is no concern about resolving turbulence, and (iii) one avoids uncertainties in the equation of state.

Other reviews
A number of other reviews of boson stars have appeared. Most recently, Schunck and Mielke [195] concentrate on the possibility of detecting BS, extending their previous reviews [163,164]. In 1992, a number of reviews appeared: Jetzer [124] concentrates on the astrophysical relevance of BS (in particular their relevance for explaining dark matter) while Liddle and Madsen [154] focus on their formation. Other reviews include [209,151]. Figure 1: Demonstration of the solitonic nature of the (mini-)boson star. Shown are snapshots of the magnitude squared of the complex scalar field for a head-on collision of two identical miniboson stars. The interacting stars display an interference pattern as they pass through each other, recovering their individual identities after the collision. However, note that the BSs have a larger amplitude after their interaction and so are not true solitons. The collision can therefore be considered inelastic. From [50]. See also [143] (e.g., Figure 5.12). 5

Solving for Boson Stars
In this section, we present the equations governing boson star solutions, namely the Einstein equations for the geometry description and the Klein-Gordon equation to represent the (complex) scalar field. We refer to this coupled system as the Einstein-Klein-Gordon (EKG) equations. The covariant equations describing boson stars are presented in Section 2.2, which is followed by choosing particular coordinates consistent with a 3+1 decomposition in Section 2.3. A form for the potential of the scalar field is then chosen and solutions are presented in Section 2.4.

Conventions
Throughout this review, Roman letters from the beginning of the alphabet a, b, c, . . . denote spacetime indices ranging from 0 to 3, while letters near the middle i, j, k, . . . range from 1 to 3, denoting spatial indices. Unless otherwise stated, we use units such that = c = 1 so that the Planck mass becomes M Planck = G −1/2 . We also use the signature convention (−, +, +, +) for the metric.

The Lagrangian, evolution equations and conserved quantities
The EKG evolution equations can be derived from the action [220] where R is the Ricci scalar of the spacetime represented by the metric g ab , and its determinant √ −g. The term L M describes the matter which here is that of a complex scalar field, φ whereφ is the complex conjugate of the field and V (|φ| 2 ) a potential depending only on the magnitude of the scalar field, consistent with the U (1) invariance of the field in the complex plane. Variation of the action in Eq. (7) with respect to the metric g ab leads to the well-known Einstein equations R ab − R 2 g ab = 8πGT ab (9) where R ab is the Ricci tensor and T ab is the real stress-energy tensor. Eqs. (9) form a system of 10 non-linear partial differential equations for the spacetime metric components g ab coupled to the scalar field via the stress-energy tensor given in Eq. (10). On the other hand, the variation of the action in Eq. (7) with respect to the scalar field φ, leads to the Klein-Gordon (KG) equation An equivalent equation is obtained when varying the action with respect to the complex conjugatē φ. The simplest potential leading to boson stars is the so-called free field case, where the potential takes the form V (|φ| 2 ) = m 2 |φ| 2 , with m a parameter that can be identified with the bare mass of the field theory.

6
According to Noether's theorem, the invariance of the Klein-Gordon Lagrangian in Eq. (8) under global U(1) transformations φ → φe iϕ implies the existence of a conserved current satisfying the conservation law The spatial integral of the time component of this current defines the conserved Noether charge, given by which can be associated with the total number of bosonic particles [189].

The 3+1 decomposition of the spacetime
Although the spacetime description of general relativity is very elegant, the covariant form of Einstein equations is not suitable to describe how an initial configuration evolves towards the future. It is therefore more intuitive to instead consider a succession of spacetime geometries, where the evolution of a given slice is given by the Einstein equations (for more detailed treatments see [5,25,35,98]). In order to convert the four-dimensional, covariant Einstein equations to a more intuitive "space+time" or 3+1 decomposition, the following steps are taken: • specify the choice of coordinates. The spacetime is foliated by a family of spacelike hypersurfaces, which are crossed by a congruence of time lines that will determine our observers (i.e., coordinates). This congruence is described by the vector field t a = αn a + β a , where α is the lapse function which measures the proper time of the observers, β a is the shift vector that measures the displacement of the observers between consecutive hypersurfaces and n a is the timelike unit vector normal to the spacelike hypersurfaces.
• decompose every 4D object into its 3+1 components. The choice of coordinates allows for the definition of a projection tensor γ a b ≡ δ a b + n a n b . Any four-dimensional tensor can be decomposed into 3+1 pieces using the spatial projector to obtain the spatial components, or contracting with n a for the time components. For instance, the line element can be written in a general form as The stress-energy tensor can then be decomposed into its various components as • write down the field equations in terms of the 3+1 components. Within the framework outlined here, the induced (or equivalently, the spatial 3D) metric γ ij and the scalar field φ are as yet still unknown (remember that the lapse and the shift just describe our choice of coordinates). In the original 3+1 decomposition (ADM formulation [10]) an additional geometrical tensor K ij ≡ − (1/2) L n γ ij = −1/ (2α) (∂ t − L β ) γ ij is introduced to describe the change of the induced metric along the congruence of observers. Loosely speaking, one can view the determination of γ ij and K ij as akin to the specification of a position and velocity for projectile motion. In terms of the extrinsic curvature and its trace, trK ≡ K i i , the Einstein equations can be written as In a similar fashion, one can introduce a quantity Q ≡ −L n φ for the Klein-Gordon equation which reduces it to an equation first order in time, second order in space • enforce any assumed symmetries. Although the boson star is found by a harmonic ansatz for the time dependence, here we choose to retain the full time-dependence. However, a considerably simplification is provided by assuming that the spacetime is spherically symmetric. Following [143], the most general metric in this case can be written in terms of spherical coordinates as where α(t, r) is the lapse function, β(t, r) is the radial component of the shift vector and a(t, r), b(t, r) represent components of the spatial metric, with dΩ 2 the metric of a unit twosphere. With this metric, the extrinsic curvature only has two independent components K i j = diag K r r , K θ θ , K θ θ . The constraint equations, Eqs. (18) and (19), can now be written as where we have defined the auxiliary scalar-field variables The evolution equations for the metric and extrinsic curvature components reduce to Similarly, the reduction of the Klein-Gordon equation to first order in time and space leads to the following set of evolution equations This set of equations, Eqs. (23) -(32), describes general, time-dependent, spherically symmetric solutions of a gravitationally-coupled complex scalar field. In the next section, we proceed to solve for the specific case of a boson star.

Mini-boson stars
The concept of a star entails a configuration of matter which remains localized. One therefore looks for a localized and time-independent matter configuration such that the gravitational field is stationary and regular everywhere. As shown in [91], such a configuration does not exist for a real scalar field. But since the stress-energy tensor depends only on the modulus of the scalar field and its gradients, one can relax the assumption of time-independence of the scalar while retaining a time-independent gravitational field. The key is to assume a harmonic ansatz for the scalar field where φ 0 is a real scalar which is the profile of the star and ω is a real constant denoting the angular frequency of the phase of the field in the complex plane. We consider spherically symmetric, equilibrium configurations corresponding to minimal energy solutions while requiring the space-time to be static. In Schwarzschild-like coordinates, the general, spherically symmetric, static metric can be written as in terms of two real metric functions, α and a. The coordinate r is an areal radius such that spheres of constant r have surface area 4πr 2 . For this reason, these coordinates are often called polar-areal coordinates. The equilibrium equations are obtained by substituting the metric of Eq. (34) and the harmonic ansatz of Eq. (33) into the spherically symmetric EKG system of Eqs. (27 -32) with β = 0, b = 1, resulting in three first order partial differential equations (PDEs) Notice that these equations hold for any stress-energy contributions and for a generic type of selfpotentials V (|φ| 2 ). In order to close the system of Eqs. (35 -37), we still have to prescribe this potential. The simplest case admitting localized solutions is the free field case of Eq. (12) for which the potential describes a field with mass m and for which the equations can be written as In order to obtain a physical solution of this system, we have to impose the following boundary conditions, which guarantee regularity at the origin and asymptotic flatness. For a given central value of the field {φ c }, we need only to adjust the eigenvalue {ω} to find a solution which matches the asymptotic behavior of Eqs. (44 -45). This system can be solved as a shooting problem by integrating from r = 0 towards the outer boundary r = r out . Eq. (39) is linear and homogeneous in α and one is therefore able to rescale the field consistent with Eq. (45). We can get rid of the constants in the equations by re-scaling the variables in the following manner Notice that the form of the metric in Eq. (34) resembles Schwarzschild allowing the association a 2 ≡ (1 − 2 M/r) −1 , where M is the ADM mass of the spacetime. This allow us to define a more general mass aspect function which measures the total mass contained in a coordinate sphere of radius r at time t.
In isotropic coordinates, the spherically symmetric metric can be written as where ψ is the conformal factor. A change of the radial coordinate R = R(r) can transform the solution obtained in Schwarzschild coordinates into isotropic ones, in particular where the first condition is the initial value to integrate the second equation backwards, obtained by imposing that far away from the boson star the spacetime resembles Schwarzschild solution.
As above, boson stars are spherically symmetric solutions of the Eqs. (38 -40) with asymptotic behaviour given by Eqs. (41 -45). For a given value of the central amplitude of the scalar field φ 0 (r = 0) = φ c , there exist configurations with some effective radius and a given mass satisfying the previous conditions for a different set of n discrete eigenvalues ω (n) . As n increases, one obtains solutions with an increasing number of nodes in φ 0 . The configuration without nodes is the ground state, while all those with any nodes are excited states. As the number of nodes increases, the distribution of the mass as a function of the radius becomes more homogeneous.
As the amplitude φ c increases, the stable configuration has a larger mass while its effective radius decreases. This trend indicates that the compactness of the boson star increases. However, at some point the mass instead decreases with increasing central amplitude. Similar to models of neutron stars (see Section 4 of [60]), this turnaround implies a maximum allowed mass for a boson star in the ground state, which numerically was found to be M max = 0.633M 2 Planck /m. The existence of a maximum mass for boson stars is a relativistic effect, which is not present in the Newtonian limit, while the maximum of baryonic stars is an intrinsic property.
Solutions for a few representative boson stars in the ground state are shown in Figure 2 in isotropic coordinates. The boson stars becomes more compact for higher values of φ c , implying narrower profiles for the scalar field, larger conformal factors, and smaller lapse functions, as the total mass increases.

Varieties of Boson Stars
Quite a number of different flavours of boson stars are present in the literature. They can have charge, a fermionic component, or rotation. They can be constructed with various potentials for the scalar field. The form of gravity which holds them together can even be modified to, say, Newtonian gravity or even no gravity at all (Q-balls). To a certain extent, such modifications are akin to varying the equation of state of a normal, fermionic star. Here we briefly review some of these variations, paying particular attention to recent work.

Self-interaction potentials
Originally, boson stars were constructed with a free-field potential without any kind of selfinteraction, obtaining a maximum mass with a dependence M ≈ M 2 Planck /m. This mass, for typical masses of bosonic particle candidates, is much smaller than the Chandrasekhar mass M Ch ≈ M 3 Planck /m 2 obtained for fermionic stars, and so they were known as mini-boson stars. In order to extend this limit and reach astrophysical masses comparable to the Chandrasekhar mass, the potential was generalized to include a self-interaction term that provided an extra pressure against gravitational collapse.
Although the first expansion to nonlinear potentials was considered in [162] including fourth and sixth power |φ|-terms, a deeper analysis was performed later considering a potential with only the quartic term [58] V |φ| with λ a dimensionless coupling constant. Written in terms of a general potential, the EKG equations remain the same. The families of gravitational equilibrium can be parametrized by the single dimensionless quantity Λ ≡ λ/ 4πGm 2 . The potential of Eq. (51) results in a maximum boson star mass that now scales as which is comparable to the Chandrasekhar mass for fermions with mass m/λ 1/4 [58]. This selfinteraction therefore allows much larger masses than the mini-boson stars as long as Λ ≫ 1, an inequality that may be satisfied even when λ ≪ 1 for reasonable scalar boson masses. The maximum mass as a function of the central value of the scalar field is shown in Figure 3 for different values of Λ. The compactness of the most massive stable stars was studied in [9], finding an upper bound M/R 0.16 for Λ ≫ 1. Figure 4 displays this compactness as a function of Λ along with the compactness of a Schwarzschild BH and non-spinning neutron star for comparison.
Many subsequent papers further analyze the EKG solutions with polynomial, or even more general non-polynomial, potentials. One work in particular [196] studied the properties of the galactic dark matter halos modeled with these boson stars. They found that a necessary condition to obtain stable, compact solutions with an exponential decrease of the scalar field, the series expansion of these potentials must contain the usual mass term m 2 |φ| 2 .
More exotic ideas similarly try to include a pressure to increase the mass of BSs. Ref.
[3] considers a form of repulsive self-interaction mediated by vector mesons within the mean-field approximation. However, the authors leave the solution of the fully nonlinear system of the Klein-Gordon and Proca equations to future work. Ref. [19] models stars made from the condensation of axions, using the semi-relativistic approach with two different potentials. Mathematically this approach involves averages such that the equations are equivalent to assuming the axion is constituted by a complex scalar field with harmonic time dependence.
Other generalizations of the potential allow for the presence of nontopological soliton solutions even in the absence of gravity, with characteristics quite different than those of the mini-boson  The compactness of a stable boson star (black solid line) as a function of the adimensional self-interaction parameter Λ ≡ λ/ 4πGm 2 . The compactness is shown for the most massive stable star (the most compact BS is unstable). This compactness asymptotes for Λ → ∞ to the value indicated by the red, dashed line. Also shown for comparison is the compactness of a Schwarzschild BH (green dot-dashed line), and the maximum compactness of a non-spinning neutron star (blue dotted line). From [9]. stars. In order to obtain these solutions the potential must satisfy two conditions. First, it must be a function of |φ| 2 to preserve the global U (1) invariance. Second, the potential should have an attractive term, bounded from below and positive for |φ| → ∞. These conditions imply a potential of at least sixth order, a condition that is satisfied by the typical degenerate vacuum form [149,91,92] V |φ| 2 = m 2 |φ| for which the potential has two degenerate minima at ±φ 0 . The case |φ| = 0 corresponds to the true vacuum state, while |φ| = φ 0 represents the degenerate vacuum state.
The resulting soliton solution can be split in three different regions. When gravity is negligible, the interior solution satisfies φ ≈ φ 0 , followed by a shell of width 1/m over which φ changes from φ 0 to zero, and an exterior that is essentially vacuum. This potential leads to a different scaling of the mass and the radius than that of the ground state [151] There is another type of non-topological soliton star, called Q-stars [157], which also admits soliton solutions in the absence of gravity (i.e., Q-balls [57,151]). The potential, besides being also a function of |φ| 2 , must satisfy the following conditions: it must behave like ≈ φ 2 near φ = 0, it has to be bounded < φ 2 in an intermediate region and must be larger > φ 2 for |φ| → ∞. The Q-stars also have three regions; an interior solution of radius R ≈ M Planck /φ 2 0 , (i.e., φ 0 ≈ m is the free particle inverse Compton wavelength) a very thin surface region of thickness 1/φ 0 , and finally the exterior solution without matter, which reduces to Schwarzschild in spherical symmetry. The mass of these Q-stars scales now as M 3 Planck /φ 2 0 . The stability of these Q-stars has been studied recently using catastrophe theory, such as [210,137]. Rotating, axisymmetric Q-balls were constructed in [135,136]. Related, rotating solutions in 2 + 1 with the signum-Gordon equation instead of the KG equation are found in [11].
Other interesting works have studied the formation of Q-balls by the Affleck-Dine mechanism [127], their dynamics in one, two and three spatial dimensions [23], and their viability as a self-interacting dark matter candidate [141].
Ref. [30] considers a chemical potential to construct BSs, arguing that the effect of the chemical potential is to reduce the parameter space of stable solutions. Related work modifies the kinetic term of the action instead of the potential. Ref. [2] studies the resulting BSs for a class of K field theories, finding solutions of two types: (i) compact balls possessing a naked singularity at their center and (ii) compact shells with a singular inner boundary which resemble black holes. Ref. [4] considers coherent states of a scalar field instead of a BS within k-essence in the context of explaining dark matter. Ref. [73] modifies the kinetic term with just a minus sign to convert the scalar field to a phantom field. Although, a regular real scalar field has no spherically symmetric, local static solutions, they find such solutions with a real phantom scalar field.

Newtonian boson stars
The Newtonian limit of the Einstein-Klein-Gordon Eqs. (9)(10)(11) can be derived by assuming that the spacetime metric in the weak field approximation can be written as where V is the Newtonian gravitational potential. In this limit, the Einstein equations reduce to the Poisson equation Conversely, by assuming that in addition to the weak limit of Eq. (55), the Klein-Gordon equation reduces to which is just the Schrödinger equation with = 1. Therefore, the EKG system is reduced in the Newtonian limit to the Schrödinger-Poisson (SP) system [100].
The initial data is obtained by solving an eigenvalue problem very similar to the one for boson stars, with similar assumptions and boundary conditions. The solutions also share similar features and display a similar behaviour. A nice property of the Newtonian limit is that all the solutions can be obtained by rescaling from one known solution [100], where N ≡ m dx 3 φφ is the Newtonian number of particles. The possibility of including self-interaction terms in the potential was considered in [108], studying also the gravitational cooling (i.e., the relaxation and virialization through the emission of scalar field bursts) of spherical perturbations. Non-spherical perturbations were further studied in [28], showing that the final state is a spherically symmetric configuration. Single Newtonian boson stars were studied in [100], either when they are boosted with/without an external central potential. Rotating stars were first successfully constructed in [204] within the Newtonian approach. Numerical evolutions of binary boson stars in Newtonian gravity are discussed in Section 4.2.
Recent work by Chavanis with Newtonian gravity solves the Gross-Pitaevskii equation, a variant of Eq. (58) which involves a pseudo-potential for a Bose-Einstein condensate, to model either dark matter or compact alternatives to neutron stars [46,47,48].
Much recent work considers boson stars from a quantum perspective as a Bose-Einstein condensate involving some number, N , of scalar fields. Ref. [161] studies the collapse of boson stars mathematically in the mean field limit in which N → ∞. Ref. [132] argues for the existence of bosonic atoms instead of stars. Ref. [17] uses numerical methods to study the mean field dynamics of BSs.

Charged boson stars
Charged boson stars result from the coupling of the boson field to the electromagnetic field [125]. The coupling between gravity and a complex scalar field with a U (1) charge arises by considering the action of Eq. 7 with the following matter Lagrangian density where e is the gauge coupling constant. The Maxwell tensor F ab can be decomposed in terms of the vector potential A a The system of equations obtained by performing the variations on the action forms the Einstein-Maxwell-Klein-Gordon system, which contains the evolution equations for the complex scalar field φ, the vector potential A a and the spacetime metric g ab [180]. Because a charged BS may be relevant for a variety of scenarios, we detail the resulting equations. For example, cosmic strings are also constructed from a charged, complex scalar field and obeys these same equations. It is only when we choose the harmonic time dependence of the scalar field that we distinguish from the harmonic azimuth of the cosmic string [219]. The evolution equations for the scalar field and for the Maxwell tensor are Notice that the vector potential is not unique; we can still add any curl-free components without changing the Maxwell equations. The gauge freedom can be fixed by choosing, for instance, the Lorentz gauge ∇ a A a = 0. Within this choice, which sets the first time derivative of the time component A 0 , the Maxwell equations reduce to a set of wave equations in a curved background with a non-linear current. This gauge choice resembles the harmonic gauge condition, which casts the Einstein equations as a system of non-linear, wave equations [220]. Either from Noether's theorem or by taking an additional covariant derivative of Eq. (63), one obtains that the electric current J a follows a conservation law. The spatial integral of the time component of this current, which can be identified with the total charge Q, is conserved. This charge is proportional to the number of particles, Q = e N . The mass M and the total charge Q can be calculated by associating the asymptotic behavior of the metric with that of Reissner-Nordström metric, which is the unique solution at large distances for a scalar field with compact support. We look for a time independent metric by first assuming a harmonically varying scalar field as in Eq. (33). We work in spherical coordinates and assume spherical symmetry. With a proper gauge choice, the vector potential takes a particularly simple form with only a single, non-trivial component A a = (A 0 (r), 0, 0, 0). This choice implies an everywhere vanishing magnetic field so that the electromagnetic field is purely electric. The boundary conditions for the vector potential are obtained by requiring the electric field to vanish at the origin because of regularity, ∂ r A 0 (r = 0) = 0. Because the electromagnetic field depends only on derivatives of the potential, we can use this freedom to set A 0 (∞) = 0 [125].
With these conditions, it is possible to find numerical solutions in equilibrium as described in Ref. [125]. Solutions are found forẽ 2 ≡ e 2 M 2 Planck /(8 π m 2 ) < 1/2. Forẽ 2 > 1/2 the repulsive Coulomb force is bigger than the gravitational attraction and no solutions are found. This bound on the BS charge in terms of its mass ensures that one cannot construct an overcharged BS, in analogy to the overcharged monopoles of Ref. [156]. An overcharged monopole is one with more charge than mass and is therefore susceptible to gravitational collapse by accreting sufficient (neutral) mass. However, because its charge is higher than its mass, such collapse might lead to an extremal Reissner-Nordström BH, but BSs do not appear to allow for this possibility. Interestingly Ref. [191] finds that if one removes gravity, the obtained Q-balls may have no limit on their charge.
The mass and the number of particles are plotted as a function of φ c for different values of e in Figure 5. Trivially, forẽ = 0 the mini-boson stars of Section (2.4) are recovered. Excited solutions with nodes are qualitatively similar [125]. The stability of these objects has been studied in [121], showing that the equilibrium configurations with a mass larger than the critical mass are dynamically unstable, similar to uncharged BSs.
Recent work with charged BSs includes the publication of Maple [1] routines to study boson nebulae charge [64,170,169] Other work generalizes the Q-balls and Q-shells found with a certain potential which leads to the signum-Gordon equation for the scalar field [133,134]. In particular, shell solutions can be found with a black hole in its interior which has implications for black hole scalar hair (for a review of black hole uniqueness see [115]).
One can also consider Q-balls coupled to an electromagnetic field, a regime appropriate for particle physics. Within such a context, Ref. [77] studies the chiral magnetic effect arising from a Q-ball. Other work in Ref. [37] studies charged, spinning Q-balls.

Oscillatons
As mentioned earlier, it is not possible to find time-independent, spacetime solutions for a real scalar field. However, there are non-singular, time-dependent near-equilibrium configurations of self-gravitating real scalar fields, which are known as oscillatons [198]. These solutions are similar to boson stars, with the exception that the spacetime must also have a time dependence in order to avoid singularities.
In this case, the system is still described by the EKG Eqs. (27)(28)(29)(30)(31)(32), with the the additional simplification that the scalar field is strictly real, φ =φ. In order to find equilibrium configurations, one expands both metric components {A(r, t) ≡ a 2 , C(r, t) ≡ (a/α) 2 } and the scalar field φ(r, t) as a truncated Fourier series where ω is the fundamental frequency and j max is the mode at which the Fourier series are trun-cated. As noted in Ref. [217,6], the scalar field consists only of odd components while the metric terms consist only of even ones. Solutions are obtained by substituting the expansions of Eq. (65) into the spherically symmetric Eqs. (27)(28)(29)(30)(31)(32). By matching terms of the same frequency, the system of equations reduce to a set of coupled ODEs. The boundary conditions are determined by requiring regularity at the origin and that the fields become asymptotically flat at large radius. These form an eigenvalue problem for the coefficients {φ 2j−1 (r = 0), A 2j (r = 0), C 2j (r = 0)} corresponding to a given central value φ 1 (r = 0). As pointed out in [217], the frequency ω is determined by the coefficient C 0 (∞) and is therefore called an output value. Although the equations are non-linear, the Fourier series converges rapidly, and so a small value of j max usually suffices. A careful analysis of the high frequency components of this construction reveals difficulties in avoiding infinite total energy while maintaining the asymptotically flat boundary condition [173]. Therefore, the truncated solutions constructed above are not exactly time periodic. Indeed, very accurate numerical work has shown that the oscillatons radiate scalar field on extremely long time scales while their frequency increases [86,99]. This work finds a mass loss rate of just one part in 10 12 per oscillation period, much too small for most numerical simulations to observe. The solutions are therefore only near-equilibrium solutions and can be extremely long-lived. Planck /m) and fundamental frequency of an oscillaton as a function of the central value of the scalar field φ 1 (r = 0). The maximum mass is M max = 0.607M 2 Planck /m. (Bottom) Plot of the total mass versus the radius at which g rr achieves its maximum. From [6].
Although the geometry is oscillatory in nature, these oscillatons behave similar to BSs. In particular, they similarly transition from long-lived solutions to a dynamically unstable branch separated at the maximum mass M max = 0.607 M 2 Planck /m. Figure 6 displays the total mass curve which shows the mass as a function of central value. Compact solutions can be found in the Newtonian framework when the weak field limit is performed appropriately, reducing to the so-called Newtonian oscillations [217]. The dynamics produced by perturbations are also qualitatively similar, including gravitational cooling, migration to more dilute stars, and collapse to black holes [6]. More recently, these studies have been extended by considering the evolution in 3D of excited states [14] and by including a quartic self-interaction potential [218]. In [131], a variational approach is used to construct oscillatons in a reduced system similar to that of the sine-Gordon breather solution.
Closely related, are oscillons that exist in flatspace and that were first mentioned as "pulsons" in 1977 [34]. There is an extensive literature on such solutions, many of which appear in [82]. A series of papers establishes that oscillons similarly radiate on very long time scales [82,83,84,85]. An interesting numerical approach to evolving oscillons adopts coordinates that blueshift and damp outgoing radiation of the massive scalar field [117,119]. A detailed look at the long term dynamics of these solutions suggests the existence of a fractal boundary in parameter space between oscillatons that lead to expansion of a true-vacuum bubble and those that disperse [118].

Rotating boson stars
Boson stars with rotation were not explored until the mid-1990s because of the lack of a strong astrophysical motivation and the technical problems with the regularization along the axis of symmetry. The first equilibrium solutions of rotating boson stars were obtained within the Newtonian gravity [204]. In order to generate axisymmetric time-independent solutions with angular momentum, one is naturally lead to the ansatz where φ 0 (r, θ) is a real scalar representing the profile of the star, ω is a real constant denoting the angular frequency of the field and k must be an integer so that the field φ is not multivalued in the azimuthal coordinate ϕ. This integer is commonly known as the rotational quantum number. General relativistic rotating boson stars were later found [194,223] with the same ansatz of Eq. (67). To obtain stationary axially symmetric solutions, two symmetries were imposed on the spacetime described by two commuting Killing vector fields ξ = ∂ t and η = ∂ ϕ in a system of adapted (cylindrical) coordinates {t, r, θ, ϕ}. In these coordinates, the metric is independent of t and ϕ and can be expressed in isotropic coordinates in the Lewis-Papapetrou form where f, l, g and Ω are metric functions depending only on r and θ. This means that we have to solve five coupled PDEs, four for the metric and one for the Klein-Gordon equation; these equations determine an elliptic quadratic eigenvalue problem in two spatial dimensions. Near the axis, the scalar field behaves as so that for k > 0 the field vanishes near the axis. Note that h k is some arbitrary function different for different values of k but no sum over k is implied in Eq. (69). This implies that the rotating star solutions have toroidal level surfaces instead of spheroidal ones as in the spherically symmetric case k = 0. In this case the metric coefficients are simplified, namely g = 1, Ω = 0 and f = f (r), l = l(r).
The entire family of solutions for k = 1 and part of k = 2 was computed using the selfconsistent field method [223], obtaining a maximum mass M max = 1.31M 2 Planck /m. Both families were completely computed in [143] using faster multigrid methods, although there were significant discrepancies in the maximum mass which indicates a problem with the regularity condition on the z-axis. The mass M and angular momentum J for stationary asymptotically flat spacetimes can be obtained from their respective Komar expressions. They can be read off from the asymptotic expansion of the metric functions f and Ω Alternatively, using the Tolman expressions for the angular momentum and the Noether charge relation in Eq. 15, one obtains an important quantization relation for the angular momentum [223] for integer values of k. This remarkable quantization condition for this classical solution also plays a role in the work of [70] discussed in Section 6.3. Figure 7 shows the scalar field for two different rotating BSs.
Recently, their stability properties were found to be similar to nonrotating stars [137]. Rotating boson stars have been shown to develop a strong ergoregion instability when rapidly spinning on short characteristic timescales (i.e., 0.1 seconds -1 week for objects with mass M = 1 − 10 6 M ⊙ and angular momentum J > 0.4GM 2 ), indicating that very compact objects with large rotation are probably black holes [45].
More discussion concerning the numerical methods and limitations of some of these approaches can be found in Lai's PhD thesis [143].

Fermionic-bosonic stars
The possibility of compact stellar objects made with a mixture of bosonic and fermionic matter was studied in Refs. [113,114]. In the simplest case, the boson component interacts with the fermionic component only via the gravitational field, although different couplings were suggested in [114] and have been further explored in [67,181]. Such a simple interaction is, at the very least, consistent with models of a bosonic dark matter coupling only gravitationally with visible matter, and the idea that such a bosonic component would become gravitationally bound within fermionic stars is arguably a natural expectation.
One can consider a perfect fluid as the fermionic component such that the stress-energy tensor takes the standard form where µ is the energy density, p is the pressure of the fluid and u a its four-velocity. Such a fluid requires an equation of state to close the system of equations (see Ref. [87] for more about fluids in relativity). In most work with fermionic-boson stars, the fluid is described by a degenerate, relativistic Fermi gas, so that the pressure is given by the parametric equation of state of Chandrasekhar where K = m 4 n /(32 π 2 ) and m n the mass of the fermion. The parameter t is given by where p o is the maximum value of the momentum in the Fermi distribution at radius r. The perfect fluid obeys relativistic versions of the Euler equations, which account for the conservation of the fluid energy and momentum, plus the conservation of the baryonic number (i.e., mass conservation). The complex scalar field representing the bosonic component is once again described by the Klein-Gordon equation. The spacetime is computed through the Einstein equations with a stress-energy tensor which is a combination of the complex scalar field and the perfect fluid After imposing the harmonic time dependence of Eq. (33) on the complex scalar field, assuming a static metric as in Eq. (34) and the static fluid u a = 0, one obtains the equations describing equilibrium fermion-boson configurations These equations can be written in adimensional form by rescaling the variables by introducing the following quantities x ≡ mr , σ(x) ≡ √ 4πGφ(0, r) , Ω ≡ ω/m 2 ,μ ≡ (4πG/m 2 )µ ,p ≡ (4πG/m 2 )p. (76) By varying the central value of the fermion density µ(r = 0) and the scalar field φ(r = 0), one finds stars dominated by either bosons or fermions, with a continuous spectrum in between. It was shown that the stability arguments made with boson stars can also be applied to these mixed objects [123]. The existence of slowly rotating fermion-boson stars was shown in [66], although no solutions were found in previous attempts [138]. Also see [74] for unstable solutions consisting of a real scalar field coupled to a perfect fluid with a polytropic equation of state.

Multi-state boson stars
It turns out that excited BSs, as dark matter halo candidates, provide for flatter, and hence more realistic, galactic rotation curves than ground state BSs. The problem is that they are generally unstable to decay to their ground state. Combining excited states with the ground states in what are called multi-state BSs is one way around this.
Although bosons in the same state are indistinguishable, it is possible to construct non-trivial configurations with bosons in different excited states. A system of bosons in P different states that only interact with each other gravitationally can be described by the following Lagrangian density where φ (n) is the particular complex scalar field representing the bosons in the n-state with n − 1 nodes. The equations of motion are very similar to the standard ones described in Section 2.2, with two peculiarities: (i) there are n independent KG equations (i.e., one for each state) and (ii) the stress-energy tensor is now the sum of contributions from each mode. Equilibrium configurations for this system were found in [26].
In the simplest case of a multi-state boson, one has the ground state and the first excited state. Such configurations are stable if the number of particles in the ground state is larger than the number of particles in the excited state [26,8] This result can be understood as the ground state deepens the gravitational potential of the excited state, and thereby stabilizing it. Unstable configurations migrate to a stable one via a flip-flop of the modes; the excited state decays, while the ground jumps to the first exited state, so that the condition (78) is satisfied. An example of this behaviour can be observed in Figure 8. Similar results were found in the Newtonian limit [216], however with a slightly higher stability limit N (1) ≥ 1.13 N (2) . This work stresses that combining several excited states makes it possible to obtain flatter rotation curves than only with ground state, producing better models for galactic dark matter halos (see also discussion of boson stars as an explanation of dark matter in Section 5.3).
Ref. [112] considers two scalar fields describing boson stars which are phase shifted in time with respect to each other, studying the dynamics numerically. In particular, one can consider multiple scalar fields with an explicit interaction (beyond just gravity) between them, say V |φ (1) | |φ (2) | . Refs. [38,39] construct such solutions, considering the individual particle-like configurations for each complex field as interacting with each other.

Alternative theories of gravity
Instead of modifying the scalar field potential, one can instead consider alternative theories of gravity. Constraints on such theories are already significant given the great success of general relativity [222]. However, the fast advance of electromagnetic observations and the anticipated gravitational wave observations promise much more in this area, in particular in the context of compact objects that probe strong-field gravity.
An ambitious effort is begun in Ref. [177] which studies a very general gravitational Lagrangian ("extended scalar-tensor theories") with both fluid stars and boson stars. The goal is for observations of compact stars to constrain such theories of gravity.
It has been found that scalar tensor theories allow for spontaneous scalarization in which the scalar component of the gravity theory transitions to a non-trivial configuration analogously to ferromagnetism with neutron stars [63]. Such scalarization is also found to occur in the context of boson star evolution [7].
Boson stars also occur within conformal gravity and with scalar-tensor extensions to it [41,42]. One motivation for alternative theories is to explain the apparent existence of dark matter without resorting to some unknown dark matter component. Perhaps the most well known of these is MOND (modified Newtonian dynamics) in which gravity is modified only at large distances [165,166] (for a review see [78]). Boson stars are studied within TeVeS (Tensor-Vector-Scalar), a relativistic generalization of MOND [59]. In particular, their evolutions of boson stars develop caustic singularities, and the authors propose modifications of the theory to avoid such problems.

Gauged boson stars
In 1988, Bartnik and McKinnon published quite unexpected results showing the existence of particle-like solutions within SU (2) Yang-Mills coupled to gravity [21]. These solutions, although unstable, were unexpected because no particle-like solutions are found in either the Yang-Mills or gravity sectors in isolation. Recall also that no particle-like solutions were found with gravity coupled to electromagnetism in early efforts to find Wheeler's geon (however, see Section 6.3 for discussion of Ref. [71] which finds geons within AdS).
Bartnik and McKinnon generalize from the Abelian U (1) gauge group to the non-Abelian SU (2) group and thereby find these unexpected particle-like solutions. One can consider, as does Ref. [195] (see Section IIp), these globally regular solutions (and their generalizations to SU (n) for n > 2) as gauged boson stars even though these contain no scalar field. One can instead explicitly include a scalar field doublet coupled to the Yang-Mills gauge field [40] as perhaps a more direct generalization of the (U (1)) charged boson stars discussed in Section 3.3.
Ref. [75] studies BSs formed from a gauge condensate of an SU (3) gauge field, and Ref. [42] extends the Bartnik-McKinnon solutions to conformal gravity with a Higgs field [42]. 24 In this section, the formation, stability and dynamical evolution of boson stars are discussed. One approach to the question of stability considers small perturbations around an equilibrium configuration, so that the system remains in the linearized regime. Growing modes indicate instability. However, a solution can be linearly stable and yet have a nonlinear instability. One example is Minkowski space which, under small perturbations, relaxes back to flat, but, for sufficiently large perturbations, leads to black hole formation, decidedly not Minkowski. To study nonlinear stability, other methods are needed. In particular, full numerical evolutions of the Einstein-Klein-Gordon (EKG) equations are quite useful for understanding the dynamics of boson stars.

Gravitational stability
A linear stability analysis consists of studying the time evolution of infinitesimal perturbations about an equilibrium configuration, usually with the additional constraint that the total number of particles must be conserved. In the case of spherically symmetric, fermionic stars described by perfect fluid, it is possible to find an eigenvalue equation for the perturbations that determines the normal modes and frequencies of the radial oscillations (see, for example, Ref. [88]). Stability theorems also allow for a direct characterization of the stability branches of the equilibrium solutions [93,61]. Analogously, one can write a similar eigenvalue equation for boson stars and show the validity of similar stability theorems. In addition to these methods, the stability of boson stars has also been studied using two other, independent methods: by applying catastrophe theory and by solving numerically the time dependent Einstein-Klein-Gordon equations. All these methods agree with the results obtained in the linear stability analysis.

Linear stability analysis
Assume that a spherically symmetric boson star in an equilibrium configuration is perturbed only in the radial direction. The equations governing these small radial perturbations are obtained by linearizing the system of equations in the standard way; expand the metric and the scalar field functions to first order in the perturbation and neglect higher order terms in the equations [95,120]. Considering the collection of fields for the system f i , one expands them in terms of the background solution 0 f i and perturbation as which assumes harmonic time dependence for the perturbation. Substitution of this expansion into the system of equations then provides a linearized system which reduces to a set of coupled equations that determines the spectrum of modes 1 f i and eigenvalues σ 2 where L ij is a differential operator containing partial derivatives and M ij is a matrix depending on the background equilibrium fields 0 f i . Solving this system, known as the pulsation equation, produces the spectrum of eigenmodes and their eigenvalues σ. The stability of the star depends crucially on the sign of the smallest eigenvalue. Because of time reversal symmetry, only σ 2 enters the equations [150], and we label the smallest eigenvalue σ 2 0 . If it is negative, the eigenmode grows exponentially with time and the star is unstable. On the other hand, for positive eigenvalues the configuration has no unstable modes and is therefore stable. The critical point at which the stability transitions from stable to unstable therefore occurs when the smallest eigenvalue vanishes, σ 0 = 0. Equilibrium solutions can be parametrized with a single variable, such as the central value of the scalar field φ c , and so we can write M = M (φ c ) and N = N (φ c ). Stability theorems then indicate that transitions between stable and unstable configurations occur only at critical points of the parameterization (M ′ (φ c ) = N ′ (φ c ) = 0) [61,93,109,208]. Linear perturbation analysis provides a more detailed picture such as the growth rates and the eigenmodes of the perturbations.
Ref. [96] carries out such an analysis for perturbations that conserve mass and charge. They find the first three perturbative modes and their growth rates, and they identify at which precise values of φ c these modes become unstable. Starting from small values, they find that ground state BSs are stable up to the critical point of maximum mass. Further increases in the central value subsequently encounter additional unstable modes. This same type of analysis applied to excited state BSs showed that the same stability criterion applies for perturbations that conserve the total particle number [122]. For more general perturbations that do not conserve particle number, excited states are generally unstable to decaying to the ground state.
A more involved analysis by [150] uses a Hamiltonian formalism to study BS stability. Considering first order perturbations that conserve mass and charge (δN = 0), their results agree with those of [96,122]. However, they extend their approach to consider more general perturbations which do not conserve the total number of particles (i.e., δN = 0). To do so, they must work with the second order quantities. They found complex eigenvalues for the excited states that indicate that excited state boson stars are unstable. More detail and discussion on the different stability analysis can be found in Ref. [124].
Catastrophe theory is part of the study of dynamical systems that began in the 1960s and studies large changes in systems resulting from small changes to certain important parameters (for a physics-oriented review see [205]). Its use in the context of boson stars is to evaluate stability, and to do so one constructs a series of solutions in terms of a limited and appropriate set of parameters. Under certain conditions, such a series generates a curve smooth everywhere except for certain points. Within a given smooth expanse between such singular points, the solutions share the same stability properties. In other words, bifurcations occur at the singular points so that solutions after the singularity gain an additional, unstable mode. Much of the recent work in this area confirms the previous conclusions from linear perturbation analysis [210,211,212,213] and from earlier work with catastrophe theory [142].
Another recent work using catastrophe theory finds that rotating stars share a similar stability picture as nonrotating solutions [137]. However, only fast spinning stars are subject to an ergoregion instability [45].

Non-linear stability of single boson stars
The dynamical evolution of spherically symmetric perturbations of boson stars has also been studied by solving numerically the Einstein-Klein-Gordon equations (Section 2.3), or its Newtonian limit (Section 3.2), the Schrödinger-Poisson system. The first such work was Ref. [197] in which the stability of the ground state was studied by considering finite perturbations which may change the total mass and the particle number (i.e., δN = 0 and δM = 0). The results corroborated the linear stability analysis in the sense that they found a stable and an unstable branch with a transition between them at a critical value, φ crit , of the central scalar field corresponding to the maximal BS mass M max = 0.633M 2 Planck /m. The perturbed configurations of the stable branch may oscillate and emit scalar radiation maintaining a characteristic frequency ν, eventually settling into some other stable state with less mass than the original. This characteristic frequency can be approximated in the non-relativistic limit as [197] where R is the effective radius of the star and M its total mass. Scalar radiation is the only damping mechanism available because spherical symmetry does not allow for gravitational radiation and because the Klein-Gordon equation has no viscous or dissipative terms. This process was named gravitational cooling, and it is extremely important in the context of formation of compact bosonic objects [199] (see below). The behavior of perturbed solutions can be represented on a plot of frequency versus effective mass as in Figure 9. Perturbed stars will oscillate with a frequency below its corresponding solid line and they radiate scalar field to infinity. As they do so, they lose mass by oscillating at constant frequency, moving leftward on the plot until they settle on the stable branch of (unperturbed) solutions. Figure 9: Oscillation frequencies of various boson stars are plotted against their mass. Also shown are the oscillation frequencies of unstable BSs obtained from the fully nonlinear evolution of the dynamical system. Unstable BSs are observed maintaining a constant frequency as they approach a stable star configuration. From [197].
The perturbed unstable configurations will either collapse to a black hole or migrate to a stable configuration, depending on the nature of the initial perturbation. If the density of the star is increased, it will collapse to a black hole. On the other hand, if it is decreased, the star explodes, expanding quickly as it approaches the stable branch. Along with the expansion, energy in the form of scalar field is radiated away, leaving a very perturbed stable star, less massive than the original unstable one.
This analysis was extended to boson stars with self-interaction and to excited BSs in Ref. [16], showing that both branches of the excited states were intrinsically unstable under generic perturbations that do not preserve M and N . The low density excited stars, with masses close to the ground state configurations, will evolve to ground state boson stars when perturbed. The more massive configurations form a black hole if the binding energy E B = M − N m is negative, through a cascade of intermediate states. The kinetic energy of the stars increases as the configuration gets closer to E B = 0, so that for positive binding energies there is an excess of kinetic energy that tends to disperse the bosons to infinity. These results are summarized in Figure 10, which shows the time scale of the excited star to decay to one of these states.  Figure 10: The instability time scale of an excited boson star (the first excitation) to one of three end states: (i) decay to the ground state, (ii) collapse to a black hole, or (iii) dispersal. From [16].
More recently, the stability of the ground state was revisited with 3D simulations using a Cartesian grid [104]. The Einstein equations were written in terms of the BSSN formulation [201,24], which is one of the most commonly used formulations in numerical relativity. Intrinsic numerical error from discretization served to perturb the ground state for both stable and unstable stars. It was found that unstable stars with negative binding energy would collapse and form a black hole, while ones with positive binding energy would suffer an excess of kinetic energy and disperse to infinity.
That these unstable stars would disperse, instead of simply expanding into some less compact stable solution, disagrees with the previous results of Ref. [197], and was subsequently further analyzed in [106] in spherical symmetry with an explicit perturbation (i.e., a Gaussian shell of particles which increases the mass of the star around 0.1%). The spherically symmetric results corroborated the previous 3D calculations, suggesting that the slightly perturbed configurations of the unstable branch have three possible endstates: (i) collapse to BH, (ii) migration to a less dense stable solution, or (iii) dispersal to infinity, dependent on the sign of the binding energy.
Closely related is the work of Lai and Choptuik [144] studying BS critical behavior (discussed in Figure 11: Very long evolutions of a perturbed, slightly sub-critical, boson star with differing outer boundaries. The central magnitude of the scalar field is shown. At early times (t < 250 and the middle frame), the boson star demonstrates near-critical behavior with small-amplitude oscillations about an unstable solution. For late times (t > 250), the solution appears converged for the largest two outer boundaries and suggests that sub-critical boson stars are not dispersing. Instead, they execute large amplitude oscillations about low-density boson stars. From [144]. Section 6.1). They tune perturbations of boson stars so that dynamically the solution approaches some particular unstable solution for some finite time. They then study evolutions which ultimately do not collapse to BH, so-called sub-critical solutions, and find that they do not disperse to infinity, instead oscillating about some less compact, stable star. They show results with increasingly distant outer boundary that suggest that this behavior is not a finite-boundary-related effect (reproduced in Figure 11). They use a different form of perturbation than Ref. [106], and, being only slightly subcritical, may be working in a regime with non-positive binding energy. However, it is interesting to consider that if indeed there are three distinct end-states, then one might expect critical behavior in the transition among the different pairings. Non-spherical perturbations of boson stars have been studied numerically in [15] with a 3D code to analyze the emitted gravitational waves.
Much less is known about rotating BSs which are more difficult to construct and to evolve because they are, at most, axisymmetric, not spherically symmetric. However, as mentioned in Section 3.5, they appear to have both stable and unstable branches [137] and are subject to an ergoregion instability at high rotation rates [45]. To our knowledge, no one has evolved rotating BS initial data. However, as discussed in the next section, simulations of BS binaries [168,174] have found rotating boson stars as a result of merger.
The issue of formation of boson stars has been addressed in [199] by performing numerical evolutions of the EKG system with different initial Gaussian distributions describing unbound states (i.e., the kinetic energy is larger than the potential energy). Quite independent of the initial condition, the scalar field collapses and settles down to a bound state by ejecting some of the scalar energy during each bounce. The ejected scalar field carries away excess ever-decreasing amounts of kinetic energy, as the system becomes bounded. After a few free-fall times of the initial configuration, the scalar field has settled into a perturbed boson star on the stable branch. This process is the already mentioned gravitational cooling, and allows for the formation of compact soliton stars (boson stars for complex scalar fields and oscillatons for real scalar fields). Although these evolutions assumed spherical symmetry which does not include important processes such as fragmentation or the formation of pancakes, they demonstrate the feasibility of the formation mechanism; clouds of scalar field will collapse under their own self-gravity while shedding excess kinetic energy. The results also confirm the importance of the mass term in the potential. By removing the massive term in the simulations, the field collapses, rebounds and completely disperses to infinity, and no compact object forms. The evolution of the scalar field with and without the massive term is displayed in Figure 12. Figure 12: The evolution of r 2 ρ (where ρ is the energy density of the complex scalar field) with massive field (top) and massless (bottom). In the massive case, much of the scalar field collapses and a perturbed boson star is formed at the center, settling down by gravitational cooling. In the massless case, the scalar field bounces through the origin and then disperses without forming any compact object. From [199].

Dynamics of binary boson stars
The dynamics of binary boson stars is sufficiently complicated that it generally requires numerical solutions. The necessary lack of symmetry and the resolution requirement dictated by the harmonic time dependence of the scalar field combine so that significant computational resources must be expended for such a study. However, boson stars serve as simple proxies for compact objects without the difficulties (shocks and surfaces) associated with perfect fluid stars, and, as such, binary BS systems have been studied in the two-body problem of general relativity. When sufficiently distant from each other, the precise structure of the star should be irrelevant as suggested by Damour's "effacement theorem" [62].
First attempts at binary boson star simulations assumed the Newtonian limit, since the SP system is simpler than the EKG one. Numerical evolutions of Newtonian binaries showed that in head-on collisions with small velocities, the stars merge forming a perturbed star [51]. With larger velocities, they demonstrate solitonic behaviour by passing through each other, producing an interference pattern during the interaction but roughly retaining their original shapes afterwards [52]. In [51] it was also compute preliminary simulation of coalescing binaries, although the lack of resolution in these 3D simulations prevents of strong statements on these results. The head-on case was re-visited in [27] with a 2D axisymmetric code. In particular, these evolutions show that final state will depend on the total energy of the system, that is, the addition of kinetic, gravitational and self-interaction energies. If the total energy is positive, the stars exhibit a solitonic behaviour both for identical (see Figure 13) as for non-identical stars. When the total energy is negative the gravitational force is the main driver of the dynamics of the system. This case produces a true collision, forming a single object with large perturbations which slowly decays by gravitational cooling, as displayed in Figure 14. Figure 13: Collision of identical boson stars with large kinetic energy in the Newtonian limit. The total energy (i.e., the addition of kinetic, gravitational and self-interaction) is positive and the collision displays solitonic behaviour. Contrast this with the gravity-dominated collision displayed in Figure 14. From [27].
The first simulations of boson stars with full general relativity were reported in [13], where the gravitational waves were computed for a head-on collision. The general behaviour is similar Figure 14: Collision of identical boson stars with small kinetic energy in the Newtonian limit. The total energy is dominated by the gravitational energy and is therefore negative. The collision leads to the formation of a single, gravitationally bound object, oscillating with large perturbations. This contrasts with the large kinetic energy case (and therefore positive total energy) displayed in Figure 13. From [27].
to the one displayed for the Newtonian limit; the stars attract each other through gravitational interaction and then merge to produce a largely perturbed boson star. However, in this case the merger of the binary was promptly followed by collapse to a black hole, an outcome not possible when working within Newtonian gravity instead of general relativity. Unfortunately, very little detail was given on the dynamics.
Much more elucidating was work in axisymmetry [143], in which head-on collisions of identical boson stars were studied in the context of critical collapse (discussed in Section 6.1) with general relativity. Stars with identical masses of M = 0.47 ≈ 0.75 M max were chosen, and so it is not surprising that for small initial momenta the stars merged together to form an unstable single star (i.e., its mass was larger than the maximum allowed mass, M max ). The unstable hypermassive star subsequently collapsed to a black hole. However, for large initial momentum the stars passed through each other, displaying a form of solitonic behaviour since the individual identities were recovered after the interaction. The stars showed a particular interference pattern during the overlap, much like that displayed in Figures 1 and 13.
Another study considered the very high speed, head-on collision of BSs [56]. Beginning with two identical boson stars boosted with Lorentz factors ranging as high as 4, the stars generally demonstrate solitonic behaviour upon collision, as shown in the insets of Figure 18. This work is further discussed in Section 6.2.
The interaction of non-identical boson stars was studied in [175] using a 3D Cartesian code to simulate head-on collisions of stars initially at rest. It was found that, for a given separation, the merger of two stars would produce an unstable star which collapses to a black hole if the initial individual mass were M ≥ 0.26 ≈ 0.4 M max . For smaller masses, the resulting star would avoid gravitational collapse and its features would strongly depend on the initial configuration. The parameterization of the initial data was written as a superposition of the single boson star solution φ 0 (r), located at different positions r 1 and r 2 Many different initial configurations are possible with this parameterization. The precise solution φ 0 is unaffected by changing the direction of rotation (within in the complex plane isospace) via ǫ = ±1 or by a phase shift δ. When ǫ = −1 the Noether charge changes sign and the compact object is then known as an anti-boson star. Three particular binary cases were studied in detail: (i) identical boson stars (ǫ = 1, δ = 0), (ii) the pair in phase opposition (ǫ = 1, δ = π), and (iii) a boson-anti-boson pair (ǫ = −1, δ = 0). The trajectories of the centers of the stars are displayed in Figure 15, together with a simple estimate of the expected trajectory assuming Newtonian gravity. The figure makes clear that the merger depends strongly on the kind of pair considered, that is, on the interaction between the scalar fields.
A simple energy argument is made in [175] to understand the differing behavior. In the weak gravity limit when the stars are well separated, one can consider the local energy density between the two stars. In addition to the contribution due to each star separately, a remaining term ∆ results from the interaction of the two stars and it is precisely this term that will depend on the parameters ǫ and δ. This term takes the simple form where ∆ 0 is a positive definite quantity. One then observes that the identical pair will have an increased energy density ∆ = +∆ 0 resulting in a deeper (and more attractive) gravitational well between the stars. In contrast, the pair with opposite phases has a decreased energy density ∆ = −∆ 0 between them, resulting in a gravitational well less attractive than the area surrounding it resulting in repulsion. The boson-anti-boson pair has an interaction that is harmonic in time ∆ = ∆ 0 cos (2ωt) and therefore sometimes positive and sometimes negative. However, if the time scale of interaction is not particularly fast, then the interaction averages to zero. Note that the boson-anti-boson pair trajectory is the closest to the simple Newtonian estimate. The qualitative behaviour agrees very well with the numerical results. The orbital case was later studied in [174]. This case is much more involved both from the computational point of view (i.e., there is less symmetry in the problem) and from the theoretical point of view, since for the final object to settle into a stationary, rotating boson star it must satisfy the additional quantization condition for the angular momentum of Eq. (71).
One simulation consisted of an identical pair each with individual mass M = 0.5, with small orbital angular momentum such that J ≤ N . In this case, the binary merges forming a rotating bar that oscillates for some time before ultimately splitting apart. This can be considered as a scattered interaction, which could not settle down to a stable boson star unless all the angular momentum was radiated.
In the case of boson-anti-boson pair, the total Noether charge is already trivial, and the final object resembles the structure of a rotating dipole. The pair in opposition of phase was not considered because of the repulsive effect from the interaction. The cases with very small angular momentum J ≪ N or with J ≤ N collapsed to a black hole soon after the merger. The trajectories for this latter case are displayed in Figure 16, indicating that the internal structure of the star is irrelevant (as per the effacement theorem [62]) until that the scalar fields overlap. A simple argument is made which qualitatively matches these numerical results, as discussed in Section 4.2. Also shown is the expected trajectory from a simple Newtonian two-body estimate. From [175].  Notice that the orbits are essentially identical at early times (and large separations), but that they start to deviate from each other on closer approach. This is consistent with the internal structure of each member of the binary being irrelevant at large separations. From [174].
Other simulations of orbiting, identical binaries have been performed within the conformally flat approximation instead of full GR, which neglects gravitational waves (GW) [168]. Three different qualitative behaviours were found. For high angular momentum, the stars orbit for comparatively long times around each other. For intermediate values, the stars merged and formed a pulsating and rotating boson star. For low angular momentum, the merger produces a black hole. No evidence was found of the stars splitting apart after the merger.

35
Scalar fields are often employed by astronomers and cosmologists in their efforts to model the Universe. Most models of inflation adopt a scalar field as the inflaton field, the vacuum energy of which drives the exponential inflation of the Universe. Dark energy also motivates many scalar field models, such as k-essence and phantom energy models. It is therefore not surprising that boson stars, as compact configurations of scalar field, are called upon to provide consequences similar to those observed.

As astrophysical stellar objects
We have already discussed a number of similarities between boson stars and models of neutron stars. Just as one can parameterize models of neutron stars by their central densities, one can consider a 1-parameter family of boson stars according to the central magnitude of the scalar field. Considering the mass as a function of this parameter, one finds the existence of a local maximum across which solutions transition from stable to unstable, just as is the case for neutron stars. Models of neutron stars can be constructed with different equations of state, whereas boson stars are constructed with differing scalar field potentials.
One difference of consequence concerns the stellar surface. Neutron stars of course have a surface at which the fluid density is discontinuous, as discussed for example in [101,103]. In contrast, the scalar field that constitutes the boson star is smooth everywhere and lacks a particular surface. In its place, one generally defines a radius which encompasses some percentage (e.g. 99%) of the stellar mass. Such a difference could have observational consequences when matter accretes onto either type of star.
It is still an open question whether some of the stars already observed and interpreted as neutron stars could instead be astrophysical boson stars. In a similar fashion, it is not known whether many, if not all, of the stars we observe already have a bosonic component that has settled into the gravitational well of the star (see Section 3.6 for a discussion of fermionic-bosonic stars). The bosonic contribution may arise from exotic matter which could appear at high densities inside the neutron star or from some sort of dark matter accretion. This possibility has gained popularity over the last years and there have been several attempts to constrain the properties of weakly interacting dark matter particles (WIMPs) by examining signatures related to their accretion and/or annihilation inside stars (for instance, see [139] and works cited in the introduction).
Recently, it was suggested that, due to the stronger gravitational field of neutron stars compared to other stars such as white dwarfs and main sequence stars, WIMPs will accrete more efficiently, leading to two different possibilities. If the dark matter is its own antiparticle, it will self-annihilate and heat the neutron star. This temperature increase could be observable in old stars, especially if they are close to the galactic center [139,65]. If WIMPs do not self-annihilate, they will settle in the center of the star forming a sort of fermionic-bosonic star. The accretion of dark matter would then increase the star's compactness until the star collapses [65].
Because of the similarities between boson stars and neutron stars, one finds that boson stars are often used in place of the other. This is especially so within numerical work because boson stars are easier to evolve than neutron star models. One can, for example compare the gravitational wave signature of a boson star merger with that of more conventional compact object binaries consisting of BHs and NSs (for a review of BH-NS binaries see [202]). At early times, the precise structure of the stars is irrelevant and the signatures are largely the same. However, for the late stages of merger, the relative phase of the boson stars determines the GW signature [175,174].
Ref. [172] follows such work by considering the result of a collision between a BH and a boson star. In particular, they consider the problem as a perturbation of a black hole via scalar accretion and analyze the resulting gravitational wave output. The hope is that observations of gravitational waves that are expected in a few years from ALIGO/VIRGO will be able to distinguish the type of matter accreting onto a BH.
With the continued advancement in observation, both in the electromagnetic and gravitational spectra, perhaps soon we will have evidence for these questions. At the same time, further study of boson stars can help identify possible distinguishing observational effects in these bands. One example where knowledge is lacking is the interaction between boson stars with a magnetic field. Whereas a neutron star can source its own magnetic field and a neutral star can obtain an induced charge when moving with respect to a magnetic field, we are aware of no studies of the interaction of boson stars with a magnetic field.

Compact alternatives to black holes
As a localized scalar field configuration, a boson star can be constructed as a non-interacting compact object, as long as one does not include any explicit coupling to any electromagnetic or other fields. In that respect, it resembles a BH, although it lacks a horizon. Can observations of purported BHs be fully explained by massive boson stars? See Ref. [184] for a review of such observations.
Neutron stars also lack horizons, but, in contrast to a boson star, have a hard surface. A hard surface is important because one would expect accretion onto such a surface to have observable consequences. Can a boson star avoid such consequences? Yuan, Narayan and Rees consider the the viability of 10M ⊙ boson stars as BH candidates in X-ray binaries [224]. They find that accreting gas collects not at the surface (which the star lacks), but instead at the center which ultimately should lead to Type I X-ray bursts. Because these bursts are not observed, the case against boson stars as black hole mimickers is weakened (at least for BH candidates in X-ray binaries).
Ref. [107] considers a simplified model of accretion and searches for boson star configurations which would mimic an accreting black hole. Although they find matches, they find that light deflection about a boson star will differ from the BH they mimic because of the lack of a photon sphere. Further work, studies the scalar field tails about boson stars and compares them to those of BHs [155]. If indeed a boson star collapses to a BH, then one could hope to observe the QNM of the massive scalar field, as described in [116].
Some of the strongest evidence for the existence of BHs is found at the center of most galaxies. The evidence suggests supermassive objects (of the order of millions of solar mass) occupying a small region (of order an astronomical unit) [33]. While definitive evidence for a BH horizon from conventional electromagnetic telescopes is perhaps just on the "horizon" [43], there are those who argue the viability of supermassive boson stars at galactic centers [215]. There could potentially be differences in the (electromagnetic) spectrum between a black hole and a boson star, but there is considerable freedom in adjusting the boson star potential to tweak the expected spectrum [105]. However, there are stringent constraints on BH alternatives to Sgr A* by the low luminosity in the near infrared [44]. In particular, the low luminosity implies a bound on the accretion rate assuming a hard surface radiating thermally, and therefore the observational evidence favors a black hole because it lacks such a surface.
However, the observation of gravitational waves from such objects may be able to distinguish BHs from BSs [29]. Such a test would occur in the bandwidth for a space-based observatory such as the beleaguered LISA mission. Because BSs allow for orbits within what would otherwise be a black hole event horizon, geodesics will exhibit extreme pericenter precession resulting in potentially distinguishable gravitational radiation [130]. In any case, observations of supermassive objects at the centers of galaxies can be used to constrain the scalar field parameters of possible mimickers [18].
There are other possible BH mimickers, and a popular recent one is the gravastar [159]. Common among all these alternatives, and most significantly, is the lack of an event horizon. Both gravastars and BSs undergo an ergoregion instability for high spin (J/(GM 2 ) > 0.4) [45]. As mentioned above for BSs, gravitational waves may similarly be able to distinguish gravastars from BHs [176].

As origin of dark matter
Studies of stellar orbits within various galaxies produce rotation curves which indicate galactic mass within the radius of the particular orbit. The discovery that these curves remain flat at large radius suggests the existence of a large halo of massive, yet dark, matter which holds the galaxy together despite its large rotation. However, what precise form of matter could fulfill the observational constraints is still very much unclear. Scalar fields are an often used tool in the cosmologist's toolkit, but one cannot have a regular, static configuration of scalar field to serve as the halo [179] (see [70] as discussed in Section 6.3 for a discussion of rotating boson stars with embedded, rotating BH solutions). Instead, galactic scale boson stars are one possible candidate.
Boson stars can be matched onto the observational constraints for galactic dark matter halos [146,200]. However, multi-state boson stars that superpose various boson star solutions (e.g. an unexcited solution with an excited solution) can perhaps find better fits to the constraints [216]. Boson stars at the galactic scale may not exhibit general relativistic effects and can be effectively considered as Bose-Einstein condensates (BEC) with angular momentum [186].
The solitonic nature of boson stars (see Figure 1) lends itself naturally to the wonderful observation of dark matter in the Bullet Cluster [147]. Ref. [148] attempts to determine a minimum galactic mass from such a match.
Interestingly Ref. [20] foregoes boson stars and instead looks for quasi-stationary scalar solutions about a Schwarzschild black hole that could conceivably survive for cosmological times. Another approach is to use scalar fields for both the dark matter halo and the supermassive, central object. Ref. [9] looks for such a match, but finds no suitable solutions. Quite a number of more exotic models viably fit within current constraints, including those using Q-balls [72].

Boson Stars in Mathematical Relativity
Although the experimental foundation for the existence of boson stars is completely lacking, on the theoretical and mathematical front, boson stars are well studied. Recent work includes a mathematical approach in terms of large and small data [89], followed up by studying singularity formation [153] and uniqueness [90,152] for a certain boson star equation. In Ref. [49], they study radial solutions of the semi-relativistic Hartree type equations in terms of global well posedness. Ref. [31] demonstrates stationarity of time periodic scalar field solutions.
Beyond just existence, however, boson stars are often employed mathematically to study dynamics. Here, we concentrate on a few of these topics which have attracted recent interest.

Black hole critical behavior
If one considers some initial distribution of energy and watches it evolve, generally one arrives at one of three states. If the energy is sufficiently weak in terms of its gravity, the energy might end up dispersing to larger and large distances. However, if the energy is instead quite large, then perhaps it will concentrate until a black hole is formed. Or, if the form of the energy supports it, some of the energy will condense into a stationary state.
In his seminal work [53], Choptuik considers a real, massless scalar field and numerically evolves various initial configurations, finding either dispersion or black hole formation. By parameterizing these initial configurations, say by the amplitude of an initial pulse p, and by tuning this parameter, he was able to study the threshold for black hole formation at which he found fascinating black hole critical behavior. In particular, his numerical work suggested that continued tuning could produce as small a black hole as one wished. This behavior is analogous to a phase transition in which the black hole mass serves as an order parameter. Similar to phase transitions, one can categorize two types of transition that distinguish between whether the black hole mass varies continuously (Type II) or discontinuously (Type I). For Choptuik's work with a massless field, the transition is therefore of Type II because the black hole mass varies from zero continuously to infinitesimal values.
Subsequent work has since established that this critical behavior can be considered as occurring in the neighborhood of a separatrix between the basins of attraction of the two end states. For p = p * , the system is precisely critical and remains on the (unstable) separatrix. Similarly other models find such threshold behavior occurring between a stationary state and black hole formation. Critical behavior about stationary solutions necessarily involve black hole formation "turning-on" at finite mass, and is therefore categorized as Type I critical behavior.
The critical surface therefore appears as a co-dimension 1 surface which evolutions increasingly approach as one tunes the parameter p. The distance from criticality |p − p * | serves as a measure of the extent to which a particular initial configuration has excited the unstable mode that drives solutions away from this surface. For Type II critical behavior, the mass of the resulting black hole mass scales as a power-law in this distance, whereas for Type I critical behavior, it is the survival time of the critical solution which scales as a power-law. See [102] for a recent review.
We have seen that boson stars represent stationary solutions of Einstein's equations and thus one would correctly guess that they may occur within Type I black hole critical behavior. To look for such behavior, Hawley and Choptuik [111] begin their evolutions with boson star solutions and then perturb them both dynamically and gravitationally. They therefore included in their evolutionary system a distinct, free, massless, real scalar field which couples to the boson star purely through its gravity.
The initial data therefore consisted of a boson star surrounded by a distant, surrounding shell of real scalar field parametrized by the amplitude of the shell. For small perturbations, the boson star oscillated about an unstable boson star before settling into a low mass, stable solution (see Figure 17). For large perturbations, the real scalar field serves to compress the initial star and, after a period of oscillating about an unstable boson star, the complex field collapses to a black hole. By tuning the initial perturbation, they find a longer and longer lived unstable boson star which serves as the critical solution (see Figure 11). The survival time τ obeys a power-law in terms of the distance from criticality |p − p * | where γ is a real constant that depends on the characteristic instability rate of the particular unstable boson star approached in the critical regime. One can also consider these BSs in axisymmetry in which non-spherically symmetric modes could potentially become important. A first step in this direction studied spherically symmetric BSs within conformally flat gravity (which does not allow for gravitational waves) in axisymmetry [188]. Later, better resolution using adaptive mesh refinement within full general relativity was achieved by [143,144] which upheld the results found within spherical symmetry. This work thus suggests that there are either no additional, unstable, axisymmetric modes or that such unstable modes are extremely slowly growing. Shown is the mass density ∂M/∂r for each contribution. By t ≈ 100 the real scalar field pulse has departed the central region and perturbed the boson star into an unstable, compact configuration. Contrast the t = 0 frame with that of t = 97.5 and note the increase in compaction. This unstable BS survives until t ≈ 500 only because the initial perturbation has been tuned to one part in 10 15 and indicates Type I critical behavior. From [144].
A very different type of critical behavior was also investigated by Lai [143]. By boosting identical boson stars toward each other and adjusting their initial momenta, he was able to tune to the threshold for black hole formation. At the threshold, he found that the time till black hole formation scaled consistent with Type I critical behavior and conjectured that the critical solution was itself an unstable boson star. This is one of the few fully nonlinear critical searches in less symmetry than spherical symmetry, and the first of Type I behavior in less symmetry. A related study colliding neutron stars instead of boson stars similarly finds Type I critical behavior [126] and subsequently confirmed by [129].
The gauged stars discussed in Section 3.9 also serve as critical solutions in spherical symmetry [54,55,167].

Hoop conjecture
An interesting use of boson stars was made by Choptuik and Pretorius [56]. They sought to answer classically whether the ultra-relativistic collision of two particles results in black hole formation. Such a question clearly has relevance to hopes of producing black holes at the LHC (see for example [145,178]). Guidance on this question is provided by Thorne's Hoop Conjecture [214] which suggests that, if one squeezes energy into some spherical space of dimension less than the Schwarzschild radius for that energy, then a black hole is formed. Figure 18: Evolutions of the head-on collisions of identical boson stars boosted toward each other with initial Lorentz factors γ as indicated. Time flows downward within each column and the top edge displays the axis of symmetry. The color-scale indicates the value of |φ|. In the middle frames one sees the interference pattern characteristic of high kinetic energy BS collisions (as mentioned in Figure 1). In the last column on the right, the collision produces a BH with apparent horizon indicated by the black oval in the third frame. From [56].
They therefore numerically collide boson stars head-on at relativistic energies to study black hole formation from just such dynamical "squeezing". Here, the nature of boson stars is largely irrelevant as they serve as simple bundles of energy that can be accelerated (see Figure 18). However, unlike using boosted black hole solutions, the choice of boson stars avoids any type of bias or predisposition to formation of a black hole. In addition, a number of previous studies of boson star head-on collisions showed interesting interference effects at energies below the threshold for black hole formation [50,52,143,168]. Indeed, it has been proposed that such an interference pattern could be evidence for the bosonic nature of dark matter because of evidence that an ideal fluid fails to produce such a pattern [97].
Choptuik and Pretorius [56] find that indeed black hole formation occurs at energies below that estimated by the Hoop Conjecture. This result is only a classical result consistent with the conjecture, but if it had not held, then there would have been no reason to expect a quantum theory to be consistent with it.

Other dimensions
Much work has been invested recently in considering physics in other dimensions. Motivation comes from various ideas including string theory (more dimensions) such as the AdS/CFT correspondence and holography (one fewer dimensions) [158,160,182]. Another source of motivation comes from the fact that higher dimensional black holes can have very different properties than those in three spatial dimensions [76]. Perhaps BSs will similarly display novel properties in other dimensions.
In lower dimensional AdS (2+1) spacetimes, early work in 1998 studied exact solutions of boson stars [192,68,193]. Higher dimensional scenarios were apparently first considered qualitatively a few years later in the context of brane world models [206]. This discussion was followed with a detailed analysis of the 3, 4, and 5 dimensional AdS solutions [12].
More recently, Ref. [81] considers oscillatons in higher dimensions and measure the scalar mass loss rate for dimensions 3, 4, and 5. They extend this work considering inflationary space-times [80].
The axisymmetric rotating BSs discussed in Section 3.5 satisfy a coupled set of nonlinear, elliptic PDEs in two dimensions. One might therefore suspect that adding another dimensions will only make things more difficult. As it turns out, however, moving to four spatial dimensions provides for another angular momentum, independent of the one along the z -direction (for example). Each of these angular momenta are associated with their own orthogonal plane of rotation. And so if one chooses solutions with equal magnitudes for each of these momenta, the solutions depend on only a single radial coordinate. This choice results in the remarkable simplification that one need only solve ODEs to find rotating solutions [140].
In [110], they extend this idea by assuming an ansatz for two complex scalar fields with equal magnitudes of angular momentum in the two independent directions. Letting the complex doublet be denoted by Φ, the ansatz takes the form Φ = φ(r)e iωt sin θe iϕ1 cos θe iϕ2 (85) in terms of the two angular coordinates ϕ 1 and ϕ 2 . One observes that the BS (i) retains a profile φ(r), (ii) possesses harmonic time dependence, and (iii) maintains single-valuedness in the two angles (the ansatz assumes a rotational quantum number of one). They find solutions that are both globally regular and asymptotically flat but these solutions appear only stable with weak gravitational coupling [110]. The work of [70] makes ingenious use of this 5D ansatz to construct rotating black holes with only a single killing vector. They set the potential of [110] to zero so that the scalar fields are massless and they add a (negative) cosmological constant to work in anti-de Sitter (AdS). They find solutions for rotating black holes in 5D AdS that correspond to a bar mode for rotating neutron stars in 3D (see also [203] for a numerical evolution of a black-hole in higher dimensions which demonstrates such bar formation). One might expect such a non-symmetric black hole to settle into a more symmetric state via the emission of gravitational waves. However, AdS provides for an essentially reflecting boundary in which the black hole can be in equilibrium. The distortion of the higher dimensional black hole also has a correspondence with the discrete values of the angular momentum of the corresponding boson star. For higher values of the rotational quantum number, the black hole develops multiple "lobes" about its center.
These solutions are extended to arbitrary odd-dimensional AdS spacetimes in [207]. Finding the solutions perturbatively, they explicitly show that these solutions approach (i) the boson star and (ii) the Myers-Perry black hole solutions in AdS [171] in different limits. See [76] for a review of black holes in higher dimensions.
The same authors of [70] also report on the existence of geons in 3+1 AdS "which can be viewed as gravitational analogs of boson stars" [71] (recall that boson stars themselves arose from Wheeler's desire to construct local electrovacuum solutions). These bundles of gravitational energy are stable to first order due to the confining boundary condition adopted with AdS. However, these geons and the black hole solutions above [70] are unstable at higher order because of the turbulent instability reported in [32].
Ref. [22] also studies black hole solutions in 5D AdS. They find solutions for black holes with scalar hair that resemble a boson star with a BH in its center. See Ref. [94] for a review of charged scalar solitons in AdS.

Final remarks
Boson stars have a long history as candidates for all manner of phenomena, from fundamental particle, to galactic dark matter. A huge variety of solutions have been found and their dynamics studied. Mathematically, BS are fascinating soliton-like solutions. Astrophysically, they represent possible explanations of black hole candidates and dark matter, with observations constraining BS properties.
Further constraints appear to be right around the corner in a few directions. The LHC is rapidly narrowing the possible mass range for a Higgs particle (the quantized version of a fundamental scalar field), and some hold hope that black holes might be produced that would indicate the existence of higher dimensions. COGENT, Pamela, and others are finding intriguing (and frustratingly contradictory) clues about the true nature of dark matter. And Planck and JWST, if it can overcome its funding and budgetary problems, promise to further refine our view of the cosmos and the role of scalar fields within. Advanced LIGO will be complete in a few years and the future of space-based GW telescopes such as LISA is currently being determined. GW observations will be ground-breaking on so many fronts, but could potentially help distinguish BSs from neutron stars or black holes.
Perhaps future work on boson stars will be experimental, if fundamental scalar fields are observed, or if evidence arises indicating the boson stars uniquely fit galactic dark matter. But regardless of any experimental results found by these remarkable experiments, there will always be regimes unexplored by experiments where boson stars will find a natural home.

Acknowledgments
It is our pleasure to thank Juan Barranco, Francisco Guzmán, and Luis Lehner for their helpful comments on the manuscript. We especially appreciate the careful and critical reading by Bruno Mundim. We also thank Gyula Fodor and Péter Forgács for their kind assistance with the section on oscillatons and oscillons. This research was supported by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. This work was also supported by NSF grants PHY-0969827 and PHY-0803624 to Long Island University.