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.


Introduction
Particle-like objects have a very long and broad history in science, arising long before Newton's corpuscles of light, and spanning the range from fundamental to astronomical. In the mid-1950s, John Wheeler sought to construct stable, particle-like solutions from only the smooth, classical fields of electromagnetism coupled to general relativity (Wheeler 1955;Power and Wheeler 1957). Such solutions would represent something of a "gravitational atom", but the solutions Wheeler found, which he called geons, were unstable (but see the discussion of geons in AdS in Sect. 6.3). However, in the following decade, Kaup replaced electromagnetism with a complex scalar field (Kaup 1968), and found Klein-Gordon geons that, in all their guises, have become well-known as today's boson stars (see Sect. II of Schunck and Mielke 2003 for a discussion of the naming history of boson stars).
As compact, stationary configurations of scalar field bound by gravity, boson stars are called upon to fill a number of different roles. Most obviously, could such solutions actually represent astrophysical objects, either observed directly or indirectly through its gravity? Instead, if constructed larger than a galaxy, could a boson star serve as the dark matter halo that explains the flat rotation curve observed for most galaxies?
The equations describing boson stars are relatively simple, and so even if they do not exist in nature, they still serve as a simple and important model for compact objects, ranging from particles to stars and galaxies. In all these cases, boson stars represent a balance between the dispersive nature of the scalar field and the attraction of gravity holding it together.
This review is organized as follows. The rest of this section describes some general features about boson stars. The system of equations describing the evolution of the scalar field and gravity (i.e., the Einstein-Klein-Gordon equations) are presented in Sect. 2. These equations are restricted to the spherical symmetric case (with a harmonic ansatz for the complex scalar field and a simple massive potential) to obtain a bosonstar family of solutions. To accommodate all their possible uses, a large variety of boson-star types have come into existence, many of which are described in more detail in Sect. 3. For example, one can vary the form of the scalar field potential to achieve a larger range of masses and compactnesses than with just a mass term in the potential. Certain types of potential admit soliton-like solutions even in the absence of gravity, leading to so-called Q-stars. One can adopt Newtonian gravity instead of general relativity, or construct solutions from a real scalar field instead of a complex one. It is also possible to find solutions coupled to an electromagnetic field or a perfect fluid, leading respectively to charged boson stars and fermion-boson stars. Rotating boson stars are found to have an angular momentum which is not arbitrary, but instead quantized, and can even coexist with a Kerr black hole. Multi-state boson stars with more than one complex scalar field are also considered. Recently, stars made of a massive vector field have been constructed which more closely match the original geon proposal because such a field has the same unit spin as Maxwell.
We discuss the dynamics of boson stars in Sect. 4. Arguably, the most important property of boson-star dynamics concerns their stability. Approaches to analyzing their stability include linear perturbation analysis, catastrophe theory, and fully nonlinear, numerical evolutions. The latter option allows for the study of the final state of perturbed stars. Possible endstates include dispersion to infinity of the scalar field, migration from unstable to stable configurations, and collapse to a black hole. There is also the question of formation of boson stars. Full numerical evolutions in 3D allow for the merger of binary boson stars, which display a large range of different behaviors as well producing distinct gravitational-wave signatures.
Finally, we review the impact of boson stars in astronomy in Sect. 5 (as astrophysical objects, black hole mimickers, gravitational-wave sources, and sources of dark matter) and in mathematics in Sect. 6 (appearing in critical behavior, the Hoop conjecture, other dimensions and anti-de Sitter spacetimes, and gravitational analogs). We conclude with some remarks and future directions.

The nature of a boson star
Boson stars (BS) are constructed with a complex scalar field coupled to gravity (as described in Sect. 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 (Kaup 1968) found energy eigenstates for a semi-classical, complex scalar field, discovering that gravitational collapse was not inevitable. Ruffini and Bonazzola (1969) 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 (Derrick 1964) (see also Rosen 1966 and its extension to the case of a general non-canonical scalar field Diez-Tejedor and Gonzalez-Morales 2013), 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 Sect. 2 the spacetime remains static. The star itself is a stationary, soliton-like solution as demonstrated in Fig. 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 Vilenkin and Shellard 1994 for a general introduction to defects and the introduction of Ryder 1996 for a discussion of relevant classical field theory concepts).
In Sect. 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 look to a fundamental scalar to provide the bosonic material of the star. Only recently has a scalar particle been experimentally found with the discovery by the Large Hadron Collider (LHC) of the standard model Higgs boson with a mass roughly 125 GeV/c 2 (Aad 2012;Chatrchyan 2012;Khachatryan et al. 2015). Of course, other proposed bosonic candidates remain, such as the 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 Fig. 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 mini-boson 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. Reprinted with permission from Choi et al. (2009). See also Lai (2004) (e.g., Figure 5.12) 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 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 ≥h.
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 ≡ 2 G M/c 2 . Substituting yields which gives an expression for the maximum mass Recognizing the Planck mass M Planck ≡ √h c/G, we obtain the estimate of M max = 0.5 M 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 Sect. 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 Sect. 2.1 and Appendix A of Brady et al. 2002), but the correspondence does not mean that the two are equivalent (Faraoni 2012). 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 (2003) concentrate on the possibility of detecting BS, extending their previous reviews (Mielke andSchunck 1999, 2002). In 1992, a number of reviews appeared: Jetzer (1992) concentrates on the astrophysical relevance of BS (in particular their relevance for explaining dark matter) while Liddle and Madsen (1992) focus on their formation. Other reviews include Straumann (1992); Lee and Pang (1992). Most recently, Mielke (2016) reviewed rotating boson stars, while Herdeiro and Radu (2015a) reviewed Kerr black holes with scalar hair.

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 Sect. 2.2, which is followed by choosing particular coordinates consistent with a 3 + 1 decomposition in Sect. 2.3. A form for the potential of the scalar field is then chosen and solutions are presented in Sect. 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 thath = 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 (Wald 1984) 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 ) is the 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 wellknown Einstein equations where R ab is the Ricci tensor and T ab is the real stress-energy tensor. Equation (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 conjugateφ. The simplest potential leading to boson stars is the so-called free field case, where the potential takes the form with m a parameter that can be identified with the bare mass of the field theory. According to Noether's theorem, the invariance of the Klein-Gordon Lagrangian in Eq. (8) under global U (1) transformations φ → φe iϕ (such that δφ = 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 (Ruffini and Bonazzola 1969). If one neglects the binding energy of the star, then the total mass can be expressed simply in terms of the bare mass as m N .

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 Alcubierre 2008; Baumgarte and Shapiro 2010;Bona et al. 2009;Gourgoulhon 2012). In order to convert the fourdimensional, 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 fourdimensional 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 τ ≡ T ab n a n b , 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 γ i j 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 Arnowitt et al. 1962) an additional geometrical tensor to describe the change of the induced metric along the congruence of observers. Loosely speaking, one can view the determination of γ i j and K i j 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 Lai (2004), 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 two-sphere. 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 Friedberg et al. (1987a), 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 spacetime to be static. In Schwarzschildlike 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 self-potentials 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 (see Dias et al. 2016 for a review on numerical methods to find stationary gravitational solutions). Equation (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. By comparing the angular metric coefficients, we also find that ψ = √ r/R. Further details can be found in Appendix D of Lai's thesis (Lai 2004). As above, boson stars are spherically symmetric solutions of the Eqs. (38-40) with asymptotic behavior given by Eqs. (41)(42)(43)(44)(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 Sect. 4 of Cook 2000), this turnaround implies a maximum allowed mass for a boson star in the ground state, which numerically was found to be M max = 0.633 M 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 Fig. 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 self-interaction, 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. To preserve the global U (1) invariance, and hence to retain a conserved particle number, such a potential should be a function of |φ|.
Although the first expansion to nonlinear potentials was considered in Mielke and Scherzer (1981) including fourth and sixth power |φ|-terms, a deeper analysis was performed later considering a potential with only the quartic term Colpi et al. (1986) 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 fermion ∼ m/λ 1/4 (Colpi et al. 1986). This self-interaction, 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 Fig. 3 for different values of Λ. The compactness of the most massive stable stars was studied in Amaro-Seoane et al. (2010), 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 black hole (BH) and non-spinning neutron star for comparison. The effect of repulsive (λ > 0) and attractive (λ < 0) quartic terms in the self-interaction potential have been studied in Eby et al. (2016).
Many subsequent papers further analyze the EKG solutions with polynomial, or even more general non-polynomial, potentials. One work in particular (Schunck and Torres 2000) 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. Agnihotri et al. (2009) 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. Barranco and Bernal (2011) 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 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 (Lee 1987;Friedberg et al. 1987b) 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 radius than that of the ground state (Lee and Pang 1992) There is another type of non-topological soliton star, called Q-stars (Lynn 1989), which also admits soliton solutions in the absence of gravity (i.e., Q-balls Coleman 1985; Lee and Pang 1992). The potential, besides being also a function of |φ|, 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 , and for some choices of the sixth order self-interaction potential the compactness of the boson star (defined with the expected value of R or R 2 ) can approach the black hole limit . The stability of these Q-stars has been studied recently using catastrophe theory, such as Tamaki and Sakai (2010); Kleihaus et al. (2012). Rotating, axisymmetric Q-balls were constructed in Kleihaus et al. (2005Kleihaus et al. ( , 2008. Related, rotating solutions in 2 + 1 with the signum-Gordon equation instead of the KG equation are found in Arodź et al. (2009). Other interesting works have studied the formation of Q-balls by the Affleck-Dine mechanism (Kasuya and Kawasaki 2000), their dynamics in one, two and three spatial dimensions (Battye and Sutcliffe 2000), and their viability as a self-interacting dark matter candidate (Kusenko and Steinhardt 2001).
It has been shown recently that very compact boson stars can also be found by using a V-shaped potential proportional to |φ| . The same V-shaped potential with an additional quadratic massive term has been considered in Kumar et al. (2015) Bhatt and Sreekanth (2009) considers a chemical potential to construct BSs, arguing that the effect of the chemical potential is to reduce the parameter space of stable solutions. Boson stars with a thermodynamically consistent equation of state, leading to an isotropic pressure, were considered in Chavanis and Harko (2012). The solutions, obtained by integrating the TOV equations, reached compactness smaller (but comparable) to neutron stars. The extension to boson stars with finite temperature was considered in Latifah et al. (2014).
Related work modifies the kinetic term of the action instead of the potential. Adam et al. (2010) 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. Akhoury and Gauthier (2008) considers coherent states of a scalar field instead of a BS within k-essence in the context of explaining dark matter. Dzhunushaliev et al. (2008) 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-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 withh = 1. Therefore, the EKG system is reduced in the Newtonian limit to the Schrödinger-Poisson (SP) system (Guenther 1995). 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 behavior. A nice property of the Newtonian limit is that all the solutions can be obtained by rescaling from one known solution (Guenther 1995), where N ≡ m dx 3 φφ is the Newtonian number of particles. The possibility of including self-interaction terms in the potential was considered in Guzmán and Ureña-López (2006), 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 Bernal and Guzmán (2006b), showing that the final state is a spherically symmetric configuration. Single Newtonian boson stars were studied in Guenther (1995), either when they are boosted with/without an external central potential. Rotating stars were first successfully constructed in Silveira and Sousa (1995) within the Newtonian approach, and, more recently, analytical approximate solutions for rotating boson stars in four and five dimensions has been achieved (Kan and Shiraishi 2016). Numerical evolutions of binary boson stars in Newtonian gravity are discussed in Sect. 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 (Chavanis 2011(Chavanis , 2012(Chavanis , 2015Chavanis and Harko 2012;Chavanis and Matos 2017). However, see a rebuttal to some of this work (Mukherjee et al. 2015).
Much recent work considers boson stars from a quantum perspective as a Bose-Einstein condensate involving some number, P, of scalar fields. Michelangeli and Schlein (2012) studies the collapse of boson stars mathematically in the mean field limit in which P → ∞. Kiessling (2009) argues for the existence of bosonic atoms instead of stars. Bao and Dong (2011) 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 (Jetzer and van der Bij 1989). 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 (Petryk 2006).
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 (Vilenkin and Shellard 1994). 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 (Wald 1984). 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 that the electric field vanishes 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 (Jetzer and van der Bij 1989).
With these conditions, it is possible to find numerical solutions in equilibrium as described in Jetzer and van der Bij (1989). It was shown that bound stable configurations exist only for values of the coupling constant less than or equal to a certain critical value, such that solutions are found forẽ 2 ≡ e 2 M 2 Planck /(8 π m 2 ) < 1/2. For e 2 > 1/2 the repulsive Coulomb force is bigger than the gravitational attraction and no solutions were found, although it has been shown recently that, due to the binding energy, solutions withẽ 2 = 1/2 and even slightly higher are also allowed (Pugliese et al. 2013). 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 Lue and Weinberg (2000). 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, Sakai and Tamaki (2012) 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ẽ in Fig. 5. Trivially, forẽ = 0 the mini-boson stars of Sect. 2.4 are recovered. Excited solutions with nodes are qualitatively similar (Jetzer and van der Bij 1989). The stability of these objects has been studied in Jetzer (1989b), 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 (Dariescu and Dariescu 2010;Murariu and Puscasu 2010;Murariu et al. 2008) and charged boson stars in the presence of a cosmological constant (Kumar et al. 2016). New regular solutions of charged scalar fields in a cavity are presented in Ponglertsakul et al. (2016), which are stable only when the radius of the mirror is sufficiently large.
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 (Kleihaus et al. 2009. 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 Chruściel et al. 2012).
One can also consider Q-balls coupled to an electromagnetic field, a regime appropriate for particle physics. Within such a context, Eto et al. (2011) studies the chiral magnetic effect arising from a Q-ball. Other work in  studies charged, spinning Q-balls.
Charged BSs in anti-de Sitter spacetimes have attracted some interest as noted at the end of Sect. 6.3.

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 nearequilibrium configurations of self-gravitating real scalar fields, which are known as oscillatons (Seidel and Suen 1991). 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 truncated. As noted in Ureña-López et al. (2002), Alcubierre et al. (2003), 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 -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 {φ 2 j−1 (r = 0), A 2 j (r = 0), C 2 j (r = 0)} corresponding to a given central value φ 1 (r = 0). As pointed out in Ureña-López et al. (2002), 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 (Page 2004). 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 (Fodor et al. 2010b;Grandclément et al. 2011). 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.
Although the geometry is oscillatory in nature, these oscillatons behave similarly 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 Planck /m. Bottom: Plot of the total mass versus the radius at which g rr achieves its maximum. Reprinted with permission from Alcubierre et al. (2003); copyright by IOP 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 (Ureña-López et al. 2002). The dynamics produced by perturbations are also qualitatively similar, including gravitational cooling, migration to more dilute stars, and collapse to black holes (Alcubierre et al. 2003). More recently, these studies have been extended by considering the evolution in 3D of excited states (Balakrishna et al. 2008) and by including a quartic self-interaction potential (Valdez-Alvarado et al. 2011). In Kichenassamy (2008), a variational approach is used to construct oscillatons in a reduced system similar to that of the sine-Gordon breather solution. Such localized solutions have also been constructed in AdS (see Sect. 6.3), and numerical evolutions suggest that they are stable below some critical density (Fodor et al. 2015).
Closely related, are oscillons that exist in flatspace and that were first mentioned as "pulsons" in Bogolyubskiȋ and Makhan'kov (1977). And so just as a Q-ball can be thought of as a BS without gravity, an oscillon is an oscillaton in the absence of gravity. Extensive literature studies such solutions, many of which appear in Fodor et al. (2008). A series of papers establishes that oscillons similarly radiate on very long time scales (Fodor et al. 2008(Fodor et al. , 2009a. An interesting numerical approach to evolving oscillons adopts coordinates that blueshift and damp outgoing radiation of the massive scalar field (Honda 2000;Honda and Choptuik 2002). 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 (Honda 2010). Dymnikova et al. (2000) examines the collision of two of these bubbles in the context of a first order phase transition. The reheating phase of inflationary cosmology generally feature oscillons which may produce observable gravitational waves (Antusch et al. 2017;Antusch and Orani 2016).

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 Newtonian gravity (Silveira and Sousa 1995). 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 (Schunck and Mielke 1996;Yoshida and Eriguchi 1997) 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 self-consistent field method (Yoshida and Eriguchi 1997), obtaining a maximum mass M max = 1.31 M 2 Planck /m. Both families were completely computed in Lai (2004) 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 (Yoshida and Eriguchi 1997) for integer values of k. This remarkable quantization condition for this classical solution also plays a role in the work of Dias et al. (2011) discussed in Sect. 6.3. Also, Smolić (2015) discusses the quantization condition of rotating BSs in the context of symmetry. Figure 7 shows the scalar field for two different rotating BSs. Spinning BS solutions with a quartic self-interacting potential have been found too, as well as their Kerr BH limit (Herdeiro et al. 2016b).

Fig. 7
The scalar field in cylindrical coordinates φ(ρ, z) for two rotating boson-star solutions: (left) k = 1 and (right) k = 2. The two solutions have roughly comparable amplitudes in scalar field. Note the toroidal shape. Reprinted with permission from Lai (2004) Recently, their stability properties were found to be similar to nonrotating stars . Rotating boson stars have been shown to develop a strong ergoregion instability when rapidly spinning on short characteristic timescales (i.e., 0.1 s-1 week for objects with mass M = 1 − 10 6 M and angular momentum J > 0.4 G M 2 ), indicating that very compact objects with large rotation are probably black holes (Cardoso et al. 2008). The presence of light rings (i.e., a region of space where photons are forced to travel in closed orbits) around rotating boson stars is studied in Grandclément (2016), while geodesics on the spacetime of these solutions are studied in Grandclément et al. (2014).
A recent review by Mielke (2016) focuses on rotating boson stars. Further discussion concerning the numerical methods and limitations of some of these approaches can also be found in the Ph.D. thesis by Lai (2004).

Fermion-boson stars
The possibility of compact stellar objects made with a mixture of bosonic and fermionic matter was studied in Henriques et al. (1989Henriques et al. ( , 1990a. In the simplest case, the bosonic component interacts with the fermionic component only via the gravitational field, although different couplings were suggested in Henriques et al. (1990a) and have been further explored in de Sousa et al. (1998), Pisano andTomazelli (1996). 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 stressenergy 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 Font 2008 for more about fluids in relativity). In most work with fermion-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 i = 0, one obtains the equations describing equilibrium fermion-boson configurations These equations can be written in adimensional form by rescaling the variables and introducing the following quantities By varying the central value of the fermion energy 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 be generalized to these mixed objects (Jetzer 1990). More recently, neutron stars with a bosonic component, sourced by dark matter accretion, have also been considered (Valdez-Alvarado et al. 2013;Brito et al. 2016). The fermionic matter for a cold star can be described easily by using simultaneously the polytropic and the ideal gas equation of state P = Kρ Γ = (Γ − 1)ρ , where ρ is the rest-mass density, its internal energy, K the polytropic constant, and Γ the adiabatic index (i.e., the energy density can be written then as μ = ρ(1 + )). For standard masses of the neutron star, stable configurations allow only about N ≈ 12%N F , where N F is the number of fermions.  (Kobayashi et al. 1994). Also see Dzhunushaliev et al. (2011) 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 nstate with n − 1 nodes. The equations of motion are very similar to the standard ones described in Sect. 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 Bernal et al. (2010).
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 Alic 2009) 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 behavior can be observed in Fig. 9. Similar results were found in the Newtonian limit (Ureña-López and Bernal 2010), 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 Sect. 5.4). Hawley and Choptuik (2003) considers two scalar fields describing boson stars that 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) | . , Brihaye and Hartmann (2009) construct such solutions, considering the individual particle-like configurations for each complex field as interacting with each other.

Proca stars
Boson stars can be understood as condensates of massive spin 0 bosonic particles modeled by a scalar field. Recently, analogous self-gravitating solutions, made of  Bernal et al. (2010); copyright by APS massive spin 1 particles, have been found in the novel work of Brito et al. (2016). These configurations, modeled by a massive complex vector field A a , are described by the Proca action for the matter sector where m is the mass of the Proca field and F ab the field strength satisfying The system of equations obtained by performing the variations on the action forms the Einstein-Proca system. The evolution equations for the Proca field are which implies that the Lorentz condition ∇ a A a = 0 is not a gauge choice like in Maxwell equations, but instead a dynamical requirement. The Einstein equations include now the stress-energy tensor Like in the scalar case, there is a global U (1) invariance of the action under transformations A a → e iθ A a , implying the existence of a conserved 4-current due to Noether's theorem In addition to carrying a conserved Noether charge, Proca stars share many other features with boson stars. Both have a harmonic time dependence but solutions exist only for a limited range of frequencies. Both achieve a maximum ADM mass, which for Proca stars is M max = 1.058M 2 Planck /m, larger, but of the same order, than those for (mini-)boson stars.

Kerr black holes with scalar hair and superradiance
Closely related to a BS, one can instead construct stable configurations of a complex scalar field around a rotating black hole (Hod 2012). Such solutions are akin to a BS with a black hole embedded at its center. As such, the scalar field serves as scalar hair (see a recent review about no-scalar-hair theorems by Herdeiro and Radu 2015b). To find such solutions, one proceeds in much the same fashion as the construction of rotating solutions (Sect. 3.5). In particular, because rotation is required to achieve a stable configuration, one works in axisymmetry and assumes a harmonic ansatz for both the internal and azimuthal rotations φ(r, t) = e −iωt e imφ ψ(r, θ).
Here ω is the (complex) angular frequency and m must be an integer (m = ± 1, ± 2, . . . ) for continuity in the azimuthal direction. Instead of solving the full system of equations, a first approximation can be obtained by solving the linearized scalar field equations on a fixed spacetime (Herdeiro and Radu 2015a). Within such a linear approximation, one finds that non-rotating (Schwarzschild) BHs do not allow for bound states with strictly real ω (Herdeiro and Radu 2014b). However, quasi-bound states can exist with (ω) < 0 in which the scalar field decays, infalling into the BH.
For a Kerr black hole with angular momentum J , mass M, and horizon radius in the equatorial plane r H , one can identify the angular velocity of the horizon as Ω H ≡ J/(2M 2 r H ). For such rotating BHs, there is a critical frequency ω c ≡ mΩ H separating disparate behavior. For ω = ω c , the frequency is strictly real allowing for regular bound states known as scalar clouds.
As ω increases above ω c , its imaginary part becomes negative, allowing again only for quasi-bound states with a time-decaying scalar field. In contrast, as ω decreases below ω c , its imaginary part becomes positive, indicating growth of the scalar field in Eq. (83). This growth of the massive field is called the superradiant instability (for a recent review of superradiance see Brito et al. 2015) and results in the extraction of energy, charge, and angular momentum from the black hole. For a rigorous treatment of this instability and a proof of boundedness see the work of Dafermos et al. (2014).
In Kühnel and Rampf (2014), an analog of a boson star (see Sect. 6.4 for physical analogs of BSs) is used to study superradiance. BSs have also been found as the zero radius limit of hairy black holes in AdS 4 (see Sect. 6.3 for BSs in AdS), and these hairy BHs are proposed as the end state of the superradiant instability (Dias and Masachs 2017).
These solutions persist when solving the fully nonlinear system in which the harmonic ansatz of Eq. (83) implies that the stress-energy tensor is independent of {t, φ}, and are generically known as Kerr BHs with scalar hair (Herdeiro and Radu 2014b). As reviewed by Herdeiro and Radu (2015a), solutions can be parametrized in such a way that connects pure Kerr BHs with pure BSs. In particular, defining q ≡ k N/J where N is the number of bosonic particles as in Eq. (15) and where k is the integer "quantum" number associated with the angular momentum as in Eq. (71), Kerr BHs are described by the vanishing of the scalar field, q = 0, and BSs are described by the vanishing of the horizon, q = 1. Fig. 11 shows the space of solutions interpolating between these two limits.
More recent work has extended these solutions. For example, a self-interacting potential with a quartic term was considered in Herdeiro et al. ( , 2016b, producing a larger amplitude scalar field but not a more massive black hole than with the non-self-interacting potential. Coupling the scalar field to the electromagnetic field allows for charged clouds (Delgado et al. 2016). Kerr black holes with Proca hair (see Sec. 3.8 for a description of Proca stars) were constructed in Herdeiro et al. (2016a). Superradiant instabilities are likely to be weaker for hairy black holes than for Kerr black holes with the same global charge (Herdeiro and Radu 2014a). A recent review on the physical properties of Kerr black holes with scalar hair can be found in Herdeiro and Radu (2015a), and prospects for testing whether BHs have hair is reviewed in Cardoso and Gualtieri (2016).
Chodosh and Shlapentokh-Rothman (2015) study the scalar cloud solutions analytically and demonstrate existence. They also consider certain uniqueness and stability properties of solutions close to Kerr and review past analytic work in this area.

Alternative theories of gravity
Instead of modifying the scalar field potential, one can consider alternative theories of gravity. Constraints on such theories are already significant given the great success of general relativity (Will 2014), and more strict bounds might be set with present and future astrophysical observations (Berti et al. 2015). However, the fast advance of electromagnetic observations and the rise of gravitational-wave astronomy promise much more in this area, in particular in the context of compact objects that probe strong-field gravity.

Fig. 11
Domain of existence for hairy black holes. The ADM mass of the solutions versus the oscillation frequency of the scalar field frequency. Solutions for a range of values of q interpolating between Kerr (q = 0) and BSs (q = 1) all with azimuthal quantum number m = 1. For 0 < q < 1, solutions describe rotating BHs surrounded by a scalar cloud, constituting scalar hair for the BH. Reprinted with permission from Herdeiro and Radu (2015a); copyright by IOP An ambitious effort is begun in Pani et al. (2011), 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. General theoretical bounds on the mass to radius ratio of stable compact objects (i.e., both neutron and boson stars) can be set for extended gravity theories, in particular for scalar tensor theories (Burikham et al. 2016).
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 (Damour and Esposito-Farèse 1996). Such scalarization is also found to occur in the context of boson-star evolution (Alcubierre et al. 2010) and scalarized hairy black holes (Kleihaus et al. 2015).
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 (Milgrom 1983(Milgrom , 2011) (for a review see Famaey and McGaugh 2012). A nonminimal coupling of the scalar field to the Ricci curvature scalar results in configurations that better resemble dark energy stars than ordinary boson stars (Horvat and Marunović 2013;Marunović 2015) . Boson stars are studied within TeVeS (Tensor-Vector-Scalar), a relativistic generalization of MOND (Contaldi et al. 2008). In particular, their evolutions of boson stars develop caustic singularities, and the authors propose modifications of the theory to avoid such problems. Bosons star solutions also exist in bi-scalar extensions of Horndeski gravity , in the framework of teleparallel gravity , and within conformal gravity and its scalar-tensor extensions Verbin 2009, 2010). Charged boson stars with torsion-coupled field have been considered in .
Recently there has been renewed interest in Einstein-Gauss-Bonet theory, which appears naturally in the low energy effective action of quantum gravity models. This theory only differs from General Relativity for dimensions D > 4, and so the easiest non-trivial case is to consider D = 5. Boson star have been found in (4 + 1)-dimensional Gauss-Bonnet gravity . Rotating configurations were constructed in Brihaye and Riedel (2014), and its classical instability and existence of ergoregions studied in Brihaye and Hartmann (2016). Rotating boson stars in odd-dimensional asymptotically anti-de Sitter spacetimes in Einstein-Gauss-Bonnet gravity are studied in Henderson et al. (2015). A non-minimal coupling between a complex scalar field and the Gauss-Bonnet term was studied in Baibhav and Maity (2017). Coupling Einstein gravity to a complex self-interacting boson field as well as a phantom field allows for new type of configurations, namely boson stars harboring a wormhole at their core (Dzhunushaliev et al. 2014).

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 (Bartnik and McKinnon 1988). 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 Sect. 6.3 for discussion of , 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 Schunck and Mielke (2003) (see Sect. 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 (Brihaye et al. 2005) as perhaps a more direct generalization of the (U (1)) charged boson stars discussed in Sect. 3.3. Dzhunushaliev et al. (2007) studies BSs formed from a gauge condensate of an SU (3) gauge field, and Brihaye and Verbin (2010) extends the Bartnik-McKinnon solutions to conformal gravity with a Higgs field.

Dynamics of boson stars
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 a 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, Font et al. 2002). Stability theorems also allow for a direct characterization of the stability branches of the equilibrium solutions (Friedman et al. 1988;Cook et al. 1994). 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 mainly two other, independent methods: by applying catastrophe theory and by solving numerically the time dependent Einstein-Klein-Gordon equations. Recently, a method utilizing information theory shows promise in analyzing the stability of equilibrium configurations. 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 (Gleiser 1988;Jetzer 1989a). 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 i j is a differential operator containing partial derivatives and M i j 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 σ . Recently, several powerful techniques have been introduced to compute the quasinormal modes of compact objects in complicated configurations, such as in the presence of interacting fields .
The stability of the star depends crucially on the sign of the smallest eigenvalue. Because of time reversal symmetry, only σ 2 enters the equations (Lee and Pang 1989), 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 of nonrotating BSs can be parametrized with a single variable, such as the central value of the scalar field φ c . We can therefore write the mass and particle number as M = M(φ c ) and N = N (φ c ), and stability theorems indicate that transitions between stable and unstable configurations occur only at the critical points in the parameter space such that These transitions in stability are completely analogous to those for neutron stars (Cook et al. 1994;Friedman et al. 1988;Harrison et al. 1965;Straumann 1984). One can generalize this result for fermion-boson stars which contain a number of fermions, N F , in addition to some number of bosons, N (see Sect. 3.6 for a discussion of fermion-boson stars). In particular, one looks for critical points in a higher dimensional parameter space by considering a vector of perturbations, n in a space spanned by the total mass at infinity, M, and the two particle numbers, N and N F . Following Henriques et al. (1990b), the critical points are such that the directional derivatives vanish where the subscript b means the value of the quantities at the critical point. The direction n at the stability boundary is tangential to the level curves of constant M and N ; formally speaking, the direction n is orthogonal to the gradient of the functions at the boundary, n ⊥ ∇(M, N , N F )| b . The condition expressed by Eq. 87 reduces to the stability condition of Eq. 86 when applied to single parameter solutions, but it allows for multi-parameter critical curves. Following the analysis of the fermion-boson star, the condition 87 implies that the equilibrium critical configurations manifest themselves at the extreme values of the number of particles when surveyed along a level curve of constant total mass (Valdez-Alvarado et al. 2013;Brito et al. 2016 where ρ c is the central density of the fermionic component. Linear perturbation analysis provides a more detailed picture such as the growth rates and the eigenmodes of the perturbations. For instance, Macedo et al. (2013a) studies the free oscillation spectra of different types of boson stars via perturbation theory. Gleiser and Watkins (1989) 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 (Jetzer 1989c). 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 Lee and Pang (1989) 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 Gleiser and Watkins (1989);Jetzer (1989c). 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 Jetzer (1992). 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 Stewart 1982). 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 Sakai 2010, 2011a, b, c) and from earlier work with catastrophe theory (Kusmartsev et al. 1991). Another recent work using catastrophe theory finds that rotating stars share a similar stability picture as nonrotating solutions ). However, only fast spinning stars are subject to an ergoregion instability (Cardoso et al. 2008).
A recent and promising alternative method to determine the stability bounds of self-gravitating astrophysical objects, and in particular of boson stars, makes use of a new measure of shape complexity known as configurational entropy (Gleiser and Jiang 2015). Their results for the critical stability region agree with those of traditional perturbation methods with an accuracy of a few percent or better.

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 (Sect. 2.3), or its Newtonian limit (Sect. 3.2), the Schrödinger-Poisson system. The first such work was Seidel and Suen (1990) 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.633 M 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 (Seidel and Suen 1990) 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 (Seidel and Suen 1994) (see below). The behavior of perturbed solutions can be represented on a plot of frequency versus effective mass as in Fig. 12. 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. 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 Fig. 12 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. Reprinted with permission from Seidel and Suen (1990); copyright by APS 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 Balakrishna et al. (1998), 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 Fig. 13, which shows the time scale of the excited star to decay to one of these states.
More recently, the stability of the ground state was revisited with 3D simulations using a Cartesian grid (Guzmán 2004). The Einstein equations were written in terms of the BSSN formulation (Shibata and Nakamura 1995; Baumgarte and Shapiro 1999),  Fig. 13 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. Reprinted with permission from Balakrishna et al. (1998); copyright by APS 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 Seidel and Suen (1990), and was subsequently further analyzed in Guzmán (2009) 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 (2007) studying BS critical behavior (discussed in Sect. 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 that 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 Fig. 14). They use a different form of perturbation than Guzmán (2009), 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 endstates, then one might expect critical behavior in the transition among the different pairings. Non-spherical perturbations of boson stars have been studied numerically in Balakrishna et al. (2006) with a 3D code to analyze the emitted gravitational waves.
The dynamics of non-standard boson stars have also been studied through numerical simulations in different scenarios. Boson stars in scalar-tensor theories of gravity were considered in Ruiz et al. (2012), focusing on the study of spontaneous and induced scalarization. Evolutions of fermion-boson stars have confirmed their stability properties and have found the normal modes of oscillations of neutron stars with a dark matter component (Valdez-Alvarado et al. 2013). Nonlinear evolutions of scalar clouds around black holes have been considered in Okawa (2015). Very recently, spherical Proca stars (see Sect. 3.8 for a discussion of such stars) have also been studied numerically (Sanchis-Gual et al. 2017), confirming that the evolutions of unstable solutions lead to outcomes analogous to those of boson stars (i.e., migration to the stable branch, total dispersion of the scalar field, or collapse to a black hole).
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 Sect. 3.5, they appear to have both stable and unstable branches  and are subject to an ergoregion instability at high rotation rates (Cardoso et al. 2008). To our knowledge, no one has evolved rotating BS initial data. However, as discussed in the next section, simulations of mini BS binaries 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. Reprinted with permission from Lai and Choptuik (2007) (Mundim 2010;Palenzuela et al. 2008) have found rotating boson stars as a result of merger.
The issue of formation of boson stars has been addressed in Seidel and Suen (1994) 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 and 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 Fig. 15.  Fig. 15 The evolution of r 2 ρ (where ρ is the energy density of the complex scalar field) with massive field (left) and massless (right). 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. Reprinted with permission from Seidel and Suen (1994); copyright by APS

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" (Damour 1987).
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 (Choi 1998). With larger velocities, they demonstrate solitonic behavior by passing through each other, producing an interference pattern during the interaction but roughly retaining their original shapes afterwards (Choi 2002). Choi (1998) simulated coalescing binaries, although the lack of resolution in these 3D simulations did not allow for strong conclusions.
The head-on case was revisited in Bernal and Guzmán (2006a) with a 2D axisymmetric code. In particular, these evolutions show that the final state will depend on the total energy of the system (e.g. the sum of kinetic, gravitational and self-interaction energies). If the total energy is positive, the stars exhibit solitonic behavior both for identical stars (see Fig. 16) and 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 Fig. 17.  Fig. 16 Collision of identical boson stars with large kinetic energy in the Newtonian limit. The total energy (i.e., the sum of kinetic, gravitational and self-interaction) is positive and the collision displays solitonic behavior. Contrast this with the gravity-dominated collision displayed in Fig. 17. Reprinted with permission from Bernal and Guzmán (2006a); copyright by APS The first simulations of boson stars with full general relativity were reported in Balakrishna (1999), where the gravitational waves were computed for a head-on collision. The general behavior is similar 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 (Lai 2004), in which head-on collisions of identical boson stars were studied in the context of critical collapse (discussed in Sect. 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 behavior since the individual identities Fig. 17 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 Fig. 16. Reprinted with permission from Bernal and Guzmán (2006a); copyright by APS were recovered after the interaction. The stars showed a particular interference pattern during the overlap, much like that displayed in Figs. 1 and 16. Another study considered the very high speed, head-on collision of BSs (Choptuik and Pretorius 2010). Beginning with two identical boson stars boosted with Lorentz factors ranging as high as 4, the stars generally demonstrate solitonic behavior upon collision, as shown in the insets of Fig. 25. This work is further discussed in Sect. 6.2.
The interaction of non-identical boson stars was studied in Palenzuela et al. (2007) 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 that 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 φ = (1) φ 0 (r 1 )e iωt + (2) φ 0 (r 2 )e i( ωt+θ) . (90)

Fig. 18 The position of the center of one BS in a head-on binary as a function of time for (i) [B-B] identical BSs, (ii) [B-poB] opposite phase pair, and (iii) [B-aB] a boson-anti-boson pair.
A simple argument is made which qualitatively matches these numerical results, as discussed in Sect. 4.2. Also shown is the expected trajectory from a simple Newtonian two-body estimate. Reprinted with permission from Palenzuela et al. (2007); copyright by APS Many different initial configurations are possible with this parameterization. The precise solution φ 0 is unaffected by changing the direction of rotation (within the complex plane) 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 Fig. 18, 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 Palenzuela et al. (2007) 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. This less attractive well results in an effective repulsion relative to the identical pairing. 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 behavior agrees very well with the numerical results. The orbital case was later studied in Palenzuela et al. (2008). 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 Fig. 19, indicating that the internal structure of the star is irrelevant (as per the effacement theorem Damour 1987) until the scalar fields overlap.
Other simulations of orbiting, identical binaries have been performed within the conformally flat approximation instead of full GR, which neglects gravitational waves (Mundim 2010). Three different qualitative behaviors 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.
Three dimensional simulations of solitonic core mergers colliding two or more boson stars in the Newtonian limit (Schrödinger-Poisson) are studied in the context of dark matter with different mass ratios, phases and orbital angular momentum (Schwabe et al. 2016). The final core mass does not depend strongly on the phase difference nor on the angular momentum. Cotner (2016) also studies collisions within the Schrödinger-Poisson system and discusses implications for dark matter. However, this work focuses on the head-on case and includes effects of different mass ratios, relative phases, selfcouplings, and separation distances. Interestingly, analytic estimates are compared to the numerical simulations (Cotner 2016). 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. Reprinted with permission from Palenzuela et al. (2008); copyright by APS The dynamics of particularly compact boson stars are interesting to contrast with the dynamics of black holes because, at least in part, we now have observations of the gravitational waves from binary BH mergers (discussed more in Sect. 5.3). To this end, the study of the head-on collision of solitonic boson stars (which can be quite compact  found the dynamics to be qualitatively similar to those observed previously with mini-boson stars (Palenzuela et al. 2007). However, the gravitational waves emitted displayed significant differences and, in some cases, closely resembled the signal from a binary black hole merger.
These studies have been extended to the orbital case in Bezares et al. (2017). Surprisingly, for stars not so massive as to collapse promptly, the merger does not lead to a rotating boson star but instead to a non-rotating perturbed BS (snapshots of some of these simulations are shown in Fig. 20). As apparent in Fig. 21, the system radiates most of its angular momentum via scalar radiation and gravitational waves soon after the merger.

Boson stars in astronomy
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. The mass is then a function of this parameter, and one finds the existence of a local maximum across which solutions transition from stable to unstable, just as is the case for neutron stars. Similarly, 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 Gundlach and Leveque (2011), Gundlach and Please (2009). 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 that 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 Sect. 3.6 for a discussion of fermion-boson 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 (Güver et al. 2014). 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 Kouvaris and Tinyakov 2010 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 (Kouvaris and Tinyakov 2010;de Lavallaz and Fairbairn 2010). If WIMPs do not self-annihilate, they will settle in the center of the star forming a fermion-boson star (as discussed in Sect. 3.6). The accretion of dark matter would then increase the star's compactness until the star collapses (de Lavallaz and Fairbairn 2010) (see discussion of BSs as a source of dark matter in Sect. 5.4). Núñez et al. (2011) 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.
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/or NSs. Differentiating BSs from other compact objects with gravitational-wave observations is discussed further in Sect. 5.3.
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 noninteracting 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 Psaltis (2008) 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 et al. (2004) consider the the viability of 10 M boson stars as BH candidates in X-ray binaries. 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). Guzmán and Rueda-Becerril (2009) considers a simplified model of accretion and searches for boson-star configurations that would mimic an accreting black hole. Although they find matches, they argue 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 (Lora-Clavijo et al. 2010). 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 Hod (2011). Differences between accretion structures surrounding boson stars and black holes are analyzed in Meliani et al. (2015), showing that the accretion tori around boson stars have different characteristics than in the vicinity of a black hole. Further studies on the subject include disk (Meliani et al. 2016) and supersonic winds (Gracia-Linares and Guzman 2016) accreting onto boson stars.
Some of the strongest evidence for the existence of BHs is found at the center of most galaxies. Observational evidence strongly suggests supermassive objects (of the order of millions of solar mass) occupying a small region (of order an astronomical unit), which is easily explained by a supermassive BH (Boehle et al. 2012). While definitive evidence for a BH horizon from conventional electromagnetic telescopes is perhaps just on the "horizon" (Johannsen et al. 2016;Broderick et al. 2011), there are those who argue for the viability of supermassive boson stars at galactic centers . 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 (Guzmán 2007). However, there are stringent constraints on BH alternatives to Sgr A* by the low luminosity in the near infrared (Broderick and Narayan 2006). 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. In particular, although a BS lacks a surface, any material it accretes would accumulate and that material would have a surface that would radiate thermally.
At least two possible ways to test the nature of astrophysical black hole candidates are apparent; either with X-ray observations or with millimeter very long baseline interferometry (VLBI).
The analysis of X-ray reflection spectroscopy with data provided by the current X-ray missions can only provide weak constraints on boson stars (Cao et al. 2016), Proca stars (Shen et al. 2017), and hairy Kerr BHs (Ni et al. 2016). The quasi-periodic oscillations (QPOs) observed in the X-ray flux emitted by accreting compact objects also provide a powerful tool both to constrain deviations from Kerr and to search for exotic compact objects. Therefore, a future eXTP mission or LOFT-like mission could set very stringent constraints on black holes with bosonic hair and on (scalar or Proca) boson stars (Franchini et al. 2017).
VLBI, on the other hand, may be able to resolve Sgr A*, the closest supermassive black hole located at the center of our galaxy. The Event Horizon Telescope (EHT) promises to resolve angular scales of the order of the horizon scale, and so soon there will be accurate images of the closest surroundings of the supermassive compact object at the center of the Galaxy. These images will allow the study of so-called BH shadows, Fig. 22 Computed images as might be expected from the EHT for: (left) a Kerr black hole and (right) a fast spinning boson star with accretion according to certain assumptions. The similarity in images indicates that ruling out a BS candidate in images of Sgr A* may prove difficult. Reprinted with permission from Vincent et al. (2016); copyright by IOP that is, the gravitational lensing and redshift effect due to the BH on the radiation from background sources.
Images of an accretion torus around Sgr A*, assuming this compact object is a boson star, are computed in Vincent et al. (2016). However, their results demonstrate that very relativistic rotating boson stars produce images extremely similar to Kerr black holes, making them difficult to distinguish from a black hole. Figure 22 displays images predicted from this work for both a BH and a BS which appear quite similar. The conclusion of Vincent et al. (2016) expresses a number of interesting caveats, and this study is also discussed as part of a more wide ranging paper about efforts to firmly establish Sgr A* as a BH (Eckart et al. 2017).
It has also been shown in Cunha et al. (2015) that hairy Kerr BHs can exhibit very distinct shadows from those of their vacuum counterparts when the light source is sufficiently far away from the BH. These differences remain, albeit less dramatically, when the BH is surrounded by an emitting torus of matter (Vincent et al. 2016).
Other studies have also studied the difference in appearance of a BS with that of the presumed BH in the center of our galaxy. Bin-Nun (2013) argues that, because BSs have an extended mass distribution that is transparent to electromagnetic radiation, the resulting strong gravitational lensing images of the S stars in the galactic center would yield much brighter images than a BH of similar mass.  studies BSs with a nonminimally coupled scalar field and makes a similar argument about bright images.
One can also consider differences between the motion of celestial bodies about BSs versus BHs. In particular, finding general geodesic motion of test particles in the space-time of boson stars generally requires numerical integration. Geodesics around non-compact boson star were studied in Diemer et al. (2013), finding additional bound orbits of massive test particles close to the center of the star that are not present in the Schwarzschild case and that could be used to make predictions about extreme-mass-ratio inspirals (EMRIs), such as the stars orbiting Sagittarius A*. One can also compute the mass parameters of compact objects from redshifts and blueshifts emitted by geodesic particles around them (Becerril et al. 2016). The motion of charged, massive test particles in the spacetime of charged boson stars was considered in .
There are other possible BH mimickers, and a popular recent one is the gravastar (Mazur and Mottola 2001). Common among all these alternatives is the lack of an event horizon. Both gravastars and BSs undergo an ergoregion instability for high spin J/(GM 2 ) > 0.4 (Cardoso et al. 2008). As mentioned above for BSs, gravitational waves may similarly be able to distinguish gravastars from BHs (Pani et al. 2009). In order to reach the high compactnesses needed to mimic a BH, one can adopt specialized potentials , but an alternative is to embed the BS within a global monopole as studied in Reid and Choptuik (2016) and Marunović and Murković (2014).

As source of gravitational waves
The era of gravitational-wave (GW) astronomy began in 2015, precisely 100 years after Einstein's development of GR. In particular, LIGO directly detected the gravitational waves from the inspiral, merger, and ringdown of a BH binary (Abbott 2016b). This observation has since been followed by a second BH binary and many more observations are expected from the facility (Abbott 2016a). The excitement about these first direct detections should also help ensure the completion of other gravitational wave observatories such as LISA (Armano 2017), KAGRA (Flaminio 2016), and next generation detectors (Abbott 2017). Now that we have actual GW observations in hand, it behooves us to extract as much science as possible from this new window on the Universe. Much work has already appeared examining the implications of these initial detections (Yunes et al. 2016;Yagi and Stein 2016;Abbott 2016c). In this paper, of course, we are concerned with the implications for BSs: (i) could these extent observations actually represent the signal from a pair of boson stars instead of BHs? (ii) might we observe a signal from boson stars, and, if so, what templates will we need? or (iii) can we place tight bounds excluding the existence of boson stars? A BS binary system is the most natural GW source. However, at early times, the precise structure of the stars is irrelevant and the signatures are largely the same whether the binary is composed of NSs, BHs, or BSs. However, during the late inspiral and merger, internal structure becomes important. In particular for boson stars, the relative phase determines the GW signature (Palenzuela et al. 2007(Palenzuela et al. , 2008. Gravitational waves may be an ideal messenger for revealing dark matter (discussed in Sect. 5.4). If new dark sector particles can form exotic compact objects (ECOs) of astronomical size, then the first evidence for such objects-and their underlying microphysical description-may arise in gravitational-wave observations. The relationship between the macroscopic properties of ECOs, such as their GW signatures, with their microscopic properties, and hence new particles, was studied in Giudice et al. (2016). The GW efficiency of compact binaries generally is examined in Hanna et al. (2017).  ; copyright by the authors Along the same lines, the tidal Love numbers for different ECOs, including different families of boson stars, are calculated in Cardoso et al. (2017). The tidal Love number, which encodes the deformability of a self-gravitating object within an external tidal field, depends significantly both on the object's internal structure and on the dynamics of the gravitational field. Present and future gravitational-wave detectors can potentially measure this quantity in a binary inspiral of compact objects and impose constraints on boson stars. Direct numerical simulations in head-on collision already have shown similarities in the gravitational waves emitted by black holes and boson stars in some cases . Fig. 23 compares the expected GW signal of a BH binary with various BS binaries.
One can also examine supermassive BHs and ask whether they could instead be some form of BS. In particular, the observation of gravitational waves from such objects may be able to distinguish BHs from BSs (Berti and Cardoso 2006). Such a test would occur in the bandwidth for a space-based observatory such as the LISA mission (Danzmann 2017). 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 (Kesden et al. 2005). In any case, observations of supermassive objects at the centers of galaxies can be used to constrain the scalar field parameters of possible mimickers (Barranco and Bernal 2011a). In Macedo et al. (2013a), the authors construct mini-boson, boson and solitonic boson stars and analyze the gravitational and scalar response of boson star spacetimes to an inspiralling stellar-mass object.

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 that holds the galaxy together despite its large rotation (see Feng 2010 for a review). 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 (Pena and Sudarsky 1997) (see Dias et al. 2011 as discussed in Sect. 6.3 for a discussion of rotating boson stars with embedded, rotating BH solutions). Instead, some form of boson star represents a possible candidate for providing the necessary dark mass.
Compact binaries are the primary target of LIGO, but instead of neutron stars or black holes, Soni and Zhang (2017) studies the expected signal from binaries consisting of SU (N ) glueball objects, one of the simplest models of dark matter. More discussion of the merger of two BSs and the production of GW can be found in Sect. 4.2. At the lower frequencies targeted by LISA, if galaxies generally possess some extended, supermassive configuration, then the inspiral of small compact body into this field will result in both dynamical friction and dark matter accretion, in addition to radiationreaction (Macedo et al. 2013b). These dynamical effects may potentially be encoded on observable gravitational waves from the inspiral.
Boson stars can be matched onto the observational constraints for galactic dark matter halos (Lee 2010a;Sharma et al. 2008). 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 (Ureña-López and . 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 (Rindler-Daller and Shapiro 2012).
Representing dark matter as BSs also offers certain computational benefits, avoiding some of the costs of modeling the particles themselves with an N -body scheme. For example, Davidson and Schwetz (2016) studies structure formation of an axion dark matter model with ground state solutions of the appropriate Schrödinger-Poisson system along with quantum pressure term (see Eq. 58). Even if dark matter consists of clumps of weakly interacting massive particles (WIMPs) instead of BSs, Mendes and Yang (Mendes and Yang 2016) map clumps of such particles to perturbed boson stars and study their tidal deformability, bypassing the large computational cost of studying the dynamics of these WIMPs with an N -body code. Tidal deformability of BSs was also studied recently in the context of testing strong-field general relativity (Cardoso et al. 2017).
Instead of galactic scale BSs, one could instead argue for the accumulation of bosonic field in neutron stars. Such solutions contain the "normal" fermionic matter as well as a bosonic component (discussed above in Sect. 3.6). However, the accumulation of additional mass in a neutron star, already the expected last stage before complete collapse to black hole, might conceivably lead to the star's collapse. If indeed collapse can be expected, then the existence of old neutron stars would place constraints on such a form of dark matter (Yz et al. 2012;Jamison 2013;Bramante et al. 2013). In the face of such arguments, Kouvaris and Tinyakov (2013), Bell et al. (2013) instead argue that a broad range of realistic models survive such constraints. Most recently, Brito et al. (2015) argue with perturbation and numerical methods that old stars are in fact stable to the accretion of light bosons by an efficient gravitational cooling mechanism (see also the Ph.D. thesis by Brito 2016). Another dark matter model arising from a scalar field is wave dark matter (Bray and Goetz 2014;Bray and Parry 2013;Goetz 2015a, b). In particular, they examine Tully-Fisher relationships predicted by this wave dark matter model (Bray and Goetz 2014;Goetz 2015b). High resolution simulations of a non-relativstic Bose-Einstein condensate within this model reproduce the large scale structure of standard cold dark matter while differing inside galaxies (Schive et al. 2014).
Other studies solve the Gross-Pitaevskii equation for a Bose-Einstein condensate as a model of dark matter stars and study its stability properties (Li et al. 2012;Madarassy and Toth 2015;Marsh and Pop 2015).
The solitonic nature of boson stars (see Fig. 1) lends itself naturally to the wonderful observation of dark matter in the Bullet Cluster (Lee et al. 2008). Lee and Lim (2010b) attempts to determine a minimum galactic mass from such a match.
Interestingly, Barranco et al. (2011) foregoes boson stars and instead looks for quasistationary 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. Amaro-Seoane et al. (2010) 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 (Doddato 2012).
Section 4.2 discusses the dynamics of boson stars including some references commenting on the implications of the dynamics for dark matter.

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 (Frank and Lenzmann 2009a), followed up by studying singularity formation (Lenzmann and Lewin 2011) and uniqueness (Frank and Lenzmann 2009b;Lenzmann 2009) for a certain boson star equation. In Cho et al. (2009), they study radial solutions of the semi-relativistic Hartree type equations in terms of global well-posedness. Bičák et al. (2010) demonstrates stationarity of time periodic scalar field solutions.
Already discussed in Sect. 3.9 has been the no hair conjecture in the context of BSs holding a central BH within. Beyond just existence, however, boson stars are often employed mathematically to study dynamics. Here, we concentrate on a few of these topics that 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 (Choptuik 1993), 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 that scales as a power law. See Gundlach et al. (2007) 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 (2000) 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 Fig. 24). For large perturbations, the real scalar field serves to compress the initial star and, after a period of oscillation 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 Fig. 14). The survival time τ obeys a power law in terms of the distance from criticality | p − p * |

Fig. 24
Evolution of a boson star (solid line) perturbed by a shell of scalar field (dashed line). 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. Reprinted with permission from Lai and Choptuik (2007) 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 (Rousseau 2003). Later, better resolution using adaptive mesh refinement within full general relativity was achieved by Lai (2004), Lai and Choptuik (2007), 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.
A very different type of critical behavior was also investigated by Lai (2004). 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 (Jin and Suen 2007) and subsequently confirmed by Kellermann et al. (2010). In the middle frames one sees the interference pattern characteristic of high kinetic energy BS collisions (as mentioned in Fig. 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. Reprinted with permission from Choptuik and Pretorius (2010); copyright by APS The gauged stars discussed in Sect. 3.11 also serve as critical solutions in spherical symmetry (Choptuik et al. 1996(Choptuik et al. , 1999Millward and Hirschmann 2003).

Hoop conjecture
An interesting use of boson stars was made by Choptuik and Pretorius (2010). 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, Landsberg 2006;Park 2012;Sirunyan 2017). Guidance on this question is provided by Thorne's Hoop Conjecture (Thorne 1972) 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.
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 Fig. 25). 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 (Choi et al. 2009;Choi 2002;Lai 2004;Mundim 2010). 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 (González and Guzmán 2011). Choptuik and Pretorius (2010) 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 and anti-de Sitter spacetime
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) (Maldacena 1998;McGreevy 2010;Polchinski 2010). Another source of motivation comes from the fact that higher dimensional black holes can have different properties than those in three spatial dimensions (Emparan and Reall 2008). 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 (Sakamoto and Shiraishi 1998a;Degura et al. 2001;Sakamoto and Shiraishi 1998). Higher dimensional scenarios were apparently first considered qualitatively a few years later in the context of brane world models (Stojkovic 2003). This discussion was followed with a detailed analysis of the 3, 4, and 5 dimensional AdS solutions (Astefanesei and Radu 2003).
More recently, Fodor et al. (2010b) considers oscillatons in higher dimensions and measures the scalar mass loss rate for dimensions 3, 4, and 5. They extend this work considering inflationary spacetimes (Fodor et al. 2010a).  and Herdeiro et al. (2015a) construct higher dimensional black hole solutions (Myers-Perry BHs) with scalar hair, and, in so doing, they find higher dimensional, rotating BS solutions.
The axisymmetric rotating BSs discussed in Sect. 3.5 satisfy a coupled set of nonlinear, elliptic PDEs in two dimensions. One might therefore suspect that adding other 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 (Kunz et al. 2006).
In Hartmann et al. (2010), 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 (93) 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 singlevaluedness 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 (Hartmann et al. 2010). Solutions have since been constructed in AdS 3 Stotyn et al. 2014a), in higher odd-dimensional AdS spacetimes (Stotyn et al. 2014b), and in Gauss-Bonnet gravity (Henderson et al. 2015) (see Sect. 3.10 for BS in alternative theories of gravity). The work of Dias et al. (2011) makes ingenious use of this 5D ansatz to construct rotating black holes with only a single Killing vector. They set the potential of Hartmann et al. (2010) to zero so that the scalar fields are massless and they add a (negative) cosmological constant to work in anti-de Sitter (AdS). Some of their solutions represent a black hole embedded inside a rotating BS. They find solutions for rotating black holes in 5D AdS that correspond to a bar mode for rotating neutron stars in 3D (see also Shibata and Yoshino 2010 for a numerical evolution of a black hole in higher dimensions which demonstrates such bar formation; see Emparan and Reall 2008 for a review of black holes in higher dimensions).
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. Very compact BSs constructed with this single Killing vector posses an ergoregion (Brihaye et al. 2015).
This construction can be extended to arbitrary odd-dimensional AdS spacetimes . 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 (Myers and Perry 1986) in different limits. Boson stars, along with neutron stars and black holes, in five dimensions are discussed in Brihaye and Delsate (2016), and see Emparan and Reall (2008) for a review of black holes in higher dimensions.
In AdS 4 this ansatz cannot be used, and the construction of spinning boson stars requires the solution of the appropriate multidimensional PDEs as is done in Radu and Subagyo (2012).
Interest in the dynamics of AdS spacetimes increased significantly with the work of Bizoń and Rostworowski (2011) who studied the collapse of a scalar field in spherically symmetric, global AdS 4 . They argued that a non-zero initial amplitude for the scalar field would generically result in gravitational collapse to black hole via turbulent instability. In particular, fully nonlinear numerical evolutions of small amplitude configurations of scalar field generically resulted in a continued sharpening of the initial pulse as it reflected off the AdS boundary. This instability in the bulk is considered the mechanism that achieves thermal equilibration in the conformal theory on the boundary.
Many studies followed trying to answer the many questions arising from this work. Did this instability extend to any initial amplitude? Was the instability tied to the precise structure of AdS or instead simply to the fact that the spacetime was bounded?
One question in particular concerned the implications of this instability for localized solutions which might naturally be expected to extend their stability in asymptotically flat spacetimes. To that end, Buchel et al. (2013) studied boson stars in AdS, and found that indeed they are stable. In the course of understanding how the boson stars were stable, this work found a whole class of initial data that appear immune to the instability. Later work added to this class, namely breather solutions in AdS (Fodor et al. 2015). Linear perturbation analysis of spherically symmetric Proca stars in AdS suggests that these too will be stable (Duarte and Brito 2016b).
The same authors of (Dias et al. 2011) also report on the existence of geons in 3 + 1 AdS "which can be viewed as gravitational analogs of boson stars" ) (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. The instability of these geons, black holes, and boson stars were studied in Dias et al. (2011) in the context of the turbulent instability reported in Bizoń and Rostworowski (2011), but later these authors argued for their nonlinear stability . Basu et al. (2010) 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. The stability of charged boson stars with a massive scalar field in five-dimensional AdS was studied in Brihaye et al. (2013). Also in AdS 5 , Buchel studies boson stars in a type IIB supergravity approximation to string theory in which the U (1) symmetry of the complex field is gauged instead of global (Buchel 2015;Buchel and Buchel 2015). A range of solutions, including Q-balls and shell solutions, for different values of the cosmological constant have similarly been constructed Riedel 2012, 2013). Basu et al. (2010) 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.
Earlier work with BSs in lower dimensional AdS was reported in Astefanesei and Radu (2003).
Boson stars in AdS with charge are constructed in Hu et al. (2012), and they are also used as the background for a study of entanglement entropy (Nogueira 2013) (for a review of holographic entanglement entropy see Rangamani and Takayanagi 2017). Charged boson stars with spin in AdS have also been studied (Kichakova et al. 2014). See Gentle et al. (2012) for a review of charged scalar solitons in AdS.

Analog gravity and physical systems
The study of the correspondence between gravitating systems and analogous physical systems goes by the name of analog gravity (Barceló et al. 2011). One example of such an analog is the acoustic or dumb hole, analogous to a black hole, that requires information to flow in a particular direction. For such a system the analog of Hawking radiation is expected, and, remarkably, such radiation may have already been measured (Unruh 2014).
Analogs exist for BS as well. Recent work of Roger et al. (2016) finds an interesting optical analog of Newtonian BSs. So far this analog appears to be mostly associated with corresponding equations of motion as opposed to some deep physical correspondence that might reveal critical insight.
A more concrete analog is the formation of a Bose-Einstein condensation such as studied in Kühnel and Rampf (2014) in the context of superradiance (see Sect. 3.9). However, note that as mentioned in Sect. 1.1, ground state BSs can be considered as condensed states of bosons without invoking any analogy (Chavanis 2015;Chavanis and Matos 2017).

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.
Remarkably, in just the last 5 years since the first version of this review, two incredibly significant experimental results have appeared, and a third may soon be on its way. The Higgs particle has been found by the LHC, the first scalar particle. While its instability makes it less than promising as the fundamental constituent of boson stars, perhaps its discovery heralds a new period of scalar discoveries. Far from the quantum particle regime of the LHC, LIGO has directly detected gravitational waves completely consistent with the merger of a binary black-hole system as predicted by general relativity. Not only does this put an end to the nagging questions about whether LIGO can really detect such extremely weak signals, but, as said often in the wake of these detections, it opens a new window into some of the most energetic events in the Universe. Although it is impossible to predict what new phenomena will be observed, one can hope that gravitational waves will further illuminate the nature of compact objects.
In the electromagnetic spectrum, the EHT just completed its 10 day observation, with an image expected by early next year (2018). The images from EHT are anxiously awaited because of their potential for demonstrating explicitly the presence of a horizon in Sgr A*. Of course, nature often surprises and so perhaps these images may instead lend credence to BSs.
With all of this experimental and observational data, physicists need to provide unambiguous tests and explicit predictions. Much work on that front is ongoing, trying to tease out observational differences from alternative models of gravity or alternatives to the standard compact objects (BHs and NSs) (Berti et al. 2015(Berti et al. , 2016Choptuik et al. 2015). Black holes were once exotic and disbelieved, but now BHs are the commonly accepted standard while BSs are proposed as just one of many exotic compact objects.
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. section on oscillatons and oscillons. SLL also thanks the Perimeter Institute for their hospitality where the first version of this work was completed. 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-1607291, PHY-1308621, PHY-0969827, and PHY-0803624 to Long Island University. CP acknowledges support from the Spanish Ministry of Economy and Competitiveness Grants FPA2013-41042-P and AYA2016-80289-P (AEI/FEDER, UE), as well as from the Spanish Ministry of Education and Science by the Ramon y Cajal Grant.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.