Chains of Boson Stars

We study axially symmetric multi-soliton solutions of a complex scalar field theory with a sextic potential, minimally coupled to Einstein's gravity. These solutions carry no angular momentum and can be classified by the number of nodes of the scalar field, $k_z$, along the symmetry axis; they are interpreted as chains with $k_z+1$ boson stars, bound by gravity, but kept apart by repulsive scalar interactions. Chains with an odd number of constituents show a spiraling behavior for their ADM mass (and Noether charge) in terms of their angular frequency, similarly to a single fundamental boson star, as long as the gravitational coupling is small; for larger coupling, however, the inner part of the spiral is replaced by a merging with the fundamental branch of radially excited spherical boson stars. Chains with an even number of constituents exhibit a truncated spiral pattern, with only two or three branches, ending at a limiting solution with finite values of ADM mass and Noether charge.


Introduction
Many non-linear physical systems support non-topological solitons, which represent spatially localized field configurations. One of the simplest examples in flat space is given by Q-balls, which are particle-like configurations in a model with a complex scalar field possessing a harmonic time dependence and a suitable self-interaction potential [1][2][3]. When Q-balls are coupled to gravity, the so-called Boson Stars (BSs) emerge, which represent solitonic solutions with a topologically trivial and globally regular geometry. The simplest such configurations are static and spherically symmetric, the scalar field possessing a mass term only, without self-interaction [4,5]. These solutions are usually dubbed mini -Boson Stars (mBS), being regarded as macroscopic quantum states, which are prevented from gravitationally collapsing by Heisenberg's uncertainty principle; also, they do not have a regular flat spacetime limit.
Both Q-balls and BSs carry a Noether charge associated with an unbroken continuous global U (1) symmetry. The charge Q is proportional to the angular frequency of the complex boson field and represents the boson particle number of the configurations [2,3].
In flat spacetime Q-balls exist only within a certain frequency range for the scalar field: between a maximal value ω max , which corresponds to the mass of the scalar excitations, and some lower non-zero critical value ω min , that depends on the form of the potential. On the one hand, as the frequency ω is approaching its extremal values, both the mass M and the Noether charge Q of the configurations diverge. On the other hand, M and Q attain a minimum at some critical value ω cr ∈ [ω min , ω max ]. Away from ω cr the mass and the charge of the configurations increase towards the divergent value at the boundary of the domain. Within [ω min , ω cr ] the configurations become unstable against decay.
The situation is different for BSs: gravity modifies the critical behavior pattern of the configurations. The fundamental branch of the BS solutions starts off from the perturbative excitations at ω ∼ ω max , at which M, Q trivialise (rather than diverge). Then, the BSs exhibit a spiraling (or oscillating) pattern of the frequency dependence of both charge and mass, where both tend to some finite limiting values at the centers of the spirals. Qualitatively, the appearance of the frequency-mass spiral may be related to oscillations in the force-balance between the repulsive scalar interaction and the gravitational attraction in equilibria. Indeed, radially excited rotating BSs do not exhibit a spiraling behavior; instead the second branch extends all the way back to the upper critical value of the frequency ω max , forming a loop [6].
The main purpose of this paper is to report on the existence of a new type of solutions, which correspond to chains of BSs. 1 These are static, axially symmetric equilibrium configurations interpreted as a number of BSs located symmetrically with respect to the origin, along the symmetry axis. We construct these solutions and investigate their physical properties for a choice of the scalar field potential with quartic and sextic selfinteraction terms, which was employed in most of the Q-balls literature. We note that similar configurations of chains of constituents are known to exist both for gravitating and flat space non-Abelian monopoles and dyons [8][9][10][11][12][13][14][15][16][17][18][19][20], Skyrmions [21][22][23], electroweak sphalerons [24][25][26][27][28], SU (2) non-self dual configurations [29,30] and Yang-Mills solitons in ADS 4 spacetime [31,32].
In these multi-component BSs configurations the constituents form a chain along the symmetry axis and, consequently, the scalar field is required to possess 'nodes' (zeros of the scalar field amplitude). Configurations with k z nodes on the symmetry axis possess k z + 1 constituents. Configurations with even and odd numbers of constituents can show a different pattern, when their domain of existence is mapped out. In particular, we find that the pattern exhibited by a mass-frequency diagram of the chains of BSs can differ both from the typical spiraling picture and from the closed loop scenario. For chains with an even number of constituents the pattern always terminates at critical solutions. For chains with an odd number of constituents, the pattern depends on the strength of the gravitational interaction. The configurations then either merge with the corresponding radial excitation of the fundamental solution, or the central constituent of the configurations exhibits oscillations while retaining smaller satellite constituents. This paper is organized as follows. In Section II the theoretical setting is specified. This includes the action, the equations of motion, the Ansatz and the boundary conditions for the BS chains. The numerical results for these new equilibrium configurations are shown in Section III . We give our conclusions in Section IV.

The model 2.1 Action, field equations and global charges
We consider a self-interacting complex scalar field Φ, which is minimally coupled to Einstein's gravity in a (3 + 1)-dimensional space-time. The corresponding action of the system is where R is the Ricci scalar curvature, G is Newton's constant, the asterisk denotes complex conjugation, and U denotes the scalar field potential. Variation of the action (2.1) with respect to the metric leads to the Einstein equations is the stress-energy tensor of the scalar field.
The corresponding equation of motion of the scalar field is the non-linear Klein-Gordon equation where represents the covariant d'Alembert operator.
The solutions considered in this work have a static line-element (with a timelike Killing vector field ξ), being topologically trivial and globally regular, i.e. without an event horizon or conical singularities, while the scalar field is finite and smooth everywhere. Also, they approach asymptotically the Minkowski spacetime background. Their mass M can be obtained from the respective Komar expressions [33], Here Σ denotes a spacelike hypersurface (with the volume element dV ), while n µ is normal to Σ, with n µ n µ = −1 [33]. After replacing in (2.5) the Ricci tensor by the stress-energy tensor, via the Einstein equations (2.2), one finds the equivalent expression The action (2.1) is invariant with respect to the global U(1) transformations of the complex scalar field, Φ → Φe iα , where α is a constant. The following Noether 4-current is associated with this symmetry It follows that integrating the timelike component of this 4-current in a spacelike slice Σ yields a second conserved quantity -the Noether charge:

The ansatz and equations
Apart from being static, the configurations in this work are also axially symmetric 2 . Thus, in a system of adapted coordinates, the space-time possesses two commuting Killing vector fields ξ = ∂/∂t and η = ∂/∂ϕ, (2.9) with t and ϕ the time and azimuthal coordinates, respectively. Line elements with these symmetries are usually studied by employing a metric ansatz [34] where (ρ, z) correspond, asymptotically, to the usual cylindrical coordinates 3 . However, it the numerical treatment of the Einstein-Klein-Gordon equations, it is useful to employ 'quasi-isotropic' spherical coordinates (r, θ), defined by the coordinate transformation in (2.10) ρ = r sin θ, z = r cos θ , (2.11) with the usual range 0 r < ∞, 0 θ ≤ π. Then the metric can then be written in a Lewis-Papapetrou form, with The three metric functions f , l, and m are functions of the variables r and θ only, chosen such that the trivial angular and radial dependence of the (asymptotically flat) line element is already factorized. The relation between the metric functions in the above line-element and those in (2.10) is f = e −2U , lr 2 sin 2 θ = P 2 , m = e 2k . The symmetry axis of the spacetime is located at the set of points with vanishing norm of η, ||η|| = 0; it corresponds to the z-axis in (2.10) or the set with θ = 0, π in (2.12). The Minkowski spacetime background is approached for r → ∞, with f = l = m = 1.
For the scalar field we adopt an Ansatz with a harmonic time dependence, while the amplitude depends on (r, θ), where ω ≥ 0 is the angular frequency. The corresponding stress-energy tensor is static, with the nonvanishing components (2.14) After inserting the ansatz (2.12), (2.13) into the field equations (2.2), (2.4) we find a system of six coupled partial differential equations that needs to be solved. There are three equations for the metric functions f, l, m, found by taking the following suitable combinations of the Einstein equations, E r r + E θ θ = 0, E ϕ ϕ = 0 and E t t = 0, yielding and one equation for the scalar field amplitude Apart from these, there are two more Einstein equations, E r θ = 0, E r r − E θ θ = 0, which are not solved in practice, being treated as constraints and used to check the numerical accuracy of the solutions.
The mass M is computed from the Komar expression (2.5), where we insert the metric ansatz (2.12), with unit vector n = − √ f dt, the volume element dV = 1/ √ f √ −g dr dθ dϕ, and √ −g = r 2 sin θ √ lm f . In evaluating (2.5), we use the fact that R t t is a total derivative: .
Then it follows that M can be read off from the asymptotic expansion of the metric function f Alternatively, the mass M can be obtained by direct integration of (2.6) , In terms of the above Ansatz the Noether charge Q is given by Also, let us note that the solutions in this work have no horizon. Therefore they are zero entropy objects, without an intrinsic temperature. Still, in analogy to black holes, one may write a "first law of thermodynamics" [35], albeit without the entropy term, which reads dM = ωdQ . (2.20) This gives a relation between the mass and Noether charge of neighbouring BS solutions which can be used to check the numerical accuracy of the solutions.

The potential and scaling properties
The solutions in this work were found for a potential originally proposed in [36,37] , which is chosen such that nontopological soliton solutions -Q-balls -exist in the absence of the Einstein term in the action (2.1), i.e. on a Minkowski spacetime background [38][39][40][41]. Also, at least in the spherically symmetric case, this choice of the potential follows for the existence of very massive and highly compact objects approaching the black hole limit [38]. The parameter µ in (2.21) corresponds to the mass of the scalar excitations around the |Φ| = 0 vacuum. Apart from that, the potential possesses two more parameters (ν, λ) > 0, determining the self-interactions, which are chosen in such a manner that it possesses a local minimum at some finite non-zero value of the field |Φ|, besides the global minimum at |Φ| = 0.
Two of the constants in (2.21) can actually be absorbed into a redefinition of the radial coordinate together with a rescaling of the scalar field, with u 0 > 0 an arbitrary constant. Note that the scalar field frequency changes accordingly, ω → u0 µ ω. Then, the potential (2.21) becomes (up to an overall factor), The choice employed in most of the Q-ball literature is which are the values used also in this work.
For completeness, let us mention that for the specific ansatz (2.12), (2.13) with the above scalings, the equations for gravitating Q-balls (which are effectively solved) can also be derived by extremizing the following reduced action 4 with [for ease of notation, we denote (∇S) · (∇T ) ≡ S ,r T ,r + 1 r 2 S ,θ T ,θ ]: The (dimensionless) coupling constant reads which determines the strength of the gravitational coupling of the solutions. As with the known spherically symmetric configurations, solutions exist for 0 α < ∞. The limit α → 0 corresponds to the non-backreacting case, i.e. Q-balls on a fixed Minkowski background. To understand the large α limit, we defineφ = φ/α. Then for large values of the effective gravitational coupling α, the non-linearity of the potential (2.21) becomes suppressed and the system approaches the usual Einstein-Klein-Gordon model [4,5], with its corresponding mBS solutions. That is, the resulting equations are identical to those found for a model with a non-self-interacting, massive complex scalar fieldφ. However, in this work we shall restrict our study to the case of a finite, nonzero α. The basic properties of the mBS chains were studied (in a more general context) in the recent work [7], while the issue of flat spacetime Q-ball chains will be discussed elsewhere.

Boundary conditions
The solutions studied in this work are globally regular and asymptotically flat, possessing finite mass and Noether charge. Appropriate boundary conditions guarantee that these conditions are satisfied.
Starting with the boundary conditions for the metric functions, regularity of the solutions at the origin requires Demanding asymptotic flatness at spatial infinity yields Axial symmetry and regularity impose on the symmetry axis at θ = 0, π the conditions Further, the condition of the absence of a conical singularity requires that the solutions should satisfy the constraint m θ=0,π = l θ=0,π . In our numerical scheme we explicitly verified (within the numerical accuracy) this condition on the symmetry axis.
Turning now to the boundary conditions for the scalar field amplitude, we mention first that φ approaches asymptotically the global minimum, while on the symmetry axis we impose The behaviour of the scalar field at the origin is more complicated, depending on the considered parity P.
As mentioned in the Introduction, the solutions split into two classes, distinguished by their behaviour w.r.t. a reflection along the equatorial plane θ = π/2. The geometry is left invariant under this transformation, for the scalar field, however, there are two possibilities: We use this symmetry to reduce the domain of integration to [0, ∞) and Finally, let us mention that, for both P = ±1, one can formally construct an approximate expression of the solutions compatible with the above boundary conditions, e.g. by assuming the existence of a power series form close to r = 0. The resulting expressions are of little help in understanding the properties of the solutions, and a numerical approach is necessary 5 . However, one important result of this study is the bound-state condition ω µ, (2.39) which emerges from the finite mass requirement. No similar result is found for the minimal value of frequency.

Numerical method
The set of four coupled non-linear elliptic partial differential equations for the functions (f, l, m; φ) has been solved numerically subject to the boundary conditions defined above. In a first stage, a new compactified radial variable x is introduced, replacing r, with with c an arbitrary parameter, typically of order one. With this choice, the semi-infinite region [0, ∞) is mapped to the finite interval [0, 1]. This avoids the use of a cutoff radius. The numerical calculations have been performed by using a using a sixth-order finite difference scheme. The system of equations is discretized on a grid with a typical number of points 129 × 89. The underlying linear system is solved with the Intel MKL PARDISO sparse direct solver [50] using the Newton-Raphson method. Calculations are performed with the packages FIDISOL/CADSOL [51] and CESDSOL 6 library. In all cases, the typical errors are of order of 10 −4 .
For the choice (2.24) of the potential's parameters, the only input parameters are the gravitational coupling constant α together with the scalar field's frequency ω. The parity of the solutions is imposed via the boundary conditions (2.34), (2.38). The number of individual constituents results from the numerical output. Also, in all plots below, quantities are given in natural units set by µ, G.

The solutions 3.1 Nodal structure and energy distribution
The choice of the parity P is related to the number of distinct constituents of the solutions (as resulting e.g. from the spatial distribution of the energy density). Moreover, denoting the number of nodes on the symmetry axis by k z , the solutions can be classified by the number of nodes k z . The number of constituents of the chains is given by k z + 1. The even-parity configurations (P = 1) have an even k z , while the solutions with an odd k z are found for P = −1. For example, the spherically symmetric solutions have k z = 0 and one single constituent localized at the origin, r = 0. The simplest non-spherical configuration has k z = 1 and represents a pair of static BSs 7 with opposite phases, the inversion of the sign of the scalar field function φ under reflections θ → π − θ corresponds to the shift of the phase ωt → ωt + π.
It was pointed out that the character of the interaction between Q-balls in Minkowski spacetime depends on their relative phase [42,43]. If the Q-balls are in phase, the interaction is attractive, if they are out of phase, there is a repulsive force between them. Thus, an axially symmetric k z = 1 solution can be balanced by the gravitational attraction.
Solutions with k z > 1 exist as well; the maximal value we have reached so far is k z = 5; however, they are likely to exist for an arbitrarily large k z .
To illustrate these aspects, we display in Fig. 1 several functions of interest for the five types of representative configurations. Both odd and even parity configurations are shown there, with the node numbers k z = 0 − 4; also the solutions have the same values of the input parameters, α = 0.25 and ω/µ = 0.8, being located on the 1st branch in the (ω, M )-diagram (as described below). These chains possess one to five constituents (from top to bottom), as seen by the number of peaks of the charge density, shown in the left panels. The middle panels represent the scalar field amplitude φ, and the right panels show the metric function f = −g tt . For the sake of clarity we have chosen to exhibit these figures in polar coordinates (ρ, z), as given by eq. (2.11).
The first row shows a single spherically symmetric BS for comparison. The second row exhibits the pair of BSs. The charge density has two symmetric peaks, the metric function has two symmetric troughs, while the scalar field function is anti-symmetric, featuring a peak and a trough. The triplet, quartet and quintet in the next few rows feature k z + 1 very similar (in size and shape) peaks for the charge density and troughs for the metric function, while the scalar field shows alternating peaks and troughs, all located symmetrically along the z-axis. Thus on the fundamental branch we basically encounter a chain consisting of k z + 1 BSs, all possessing similar size, shape and distance from their next neighbors.
This picture partially changes as we move along the domain of solutions, for a given coupling α. This is seen in Fig. 2, where these chains (k z = 0 − 4) are now shown for illustrative solutions sitting on the second branch of the M vs. ω domain of existence (see the discussion in the next subsection), with α = 0.25 for k z = 0 − 3, α = 0.5 for k z = 4, and ω/µ = {0.43, 0.47, 0.57, 0.7, 0.7} for k z = 0 − 4, respectively. As we look at the charge density of the configurations, we see a dominant peak at the center for odd chains and a dominant inner pair for even chains, while the other peaks of the triplet, quartet and quintet have turned into slight elevations, hardly visible in the figures. In the metric functions the troughs at the center of the odd chains dominate, and likewise the (almost merged) inner pair of troughs of the even chains. All other troughs are weakened substantially. The scalar field itself, however, retains the outer peaks and troughs to a somewhat greater extent, still reflecting clearly the number of constituents of the chains.

The ω-dependence and the branch structure
Recall the frequency dependence for the single BSs and for fixed coupling α [38]. The set of BSs emerges from the vacuum at the maximal frequency, given by the boson mass ω max = µ. Thus, unlike the case of Q-balls in flat space, where mass and charge diverge, these quantitites vanish in this limit. Decreasing the frequency ω spans the first or fundamental branch, which terminates at the first backbending of the curve, at which point it moves towards larger frequencies. The curve then follows a spiraling/oscillating pattern, with successive backbendings.
The mass and charge form a spiral, as ω is varied, while the minimum of the metric function f and the maximum of the scalar function φ shows damped oscillations. The set of solutions tends to a limiting solution at the center of the spiral which has finite values of the mass and charge. However, the values of the scalar field function φ and the metric function f at the center of the star, which represent the maximal and minimal values of these functions, φ max and f min , respectively, do not seem to be finite, with φ max diverging and f min vanishing in the limit.
Let us now consider the frequency dependence of the BS chains, when the coupling α is kept fixed. Like the single BSs, all chains emerge similarly from the vacuum at the maximal frequency, given by the boson mass ω max = µ, where their mass and charge vanish in the limit.
When ω is decreased, mass and charge rise and the chains follow along their fundamental branch. As for the single BS, for all chains this fundamental branch ends at a minimal value of the frequency, from where a second branch arises. But then even and odd chains will in general exhibit different patterns, which will also depend on the coupling strength α.

Even chains
Let us consider the BS pair and the higher even chains in more detail. We illustrate the ω-dependence for even chains in Fig. 3, selecting the BS pair (k z = 1) and the BS quartet (k z = 3), and comparing with the single BS. In the upper panels we show the k z = 1 pair (solid curves), the k z = 3 quartet (dash-dotted curves), as well as the k z = 0 single BS (dashed curves). For the latter we only show the first few branches. In the two upper panels we exhibit the scaled ADM mass M (left) and the scaled charge Q (right). The different colors refer to different values of the coupling α. The middle left panel shows in an analogous manner the minimal value f min of the metric function f . Restricting to the k z = 1 pair, we then exhibit the maximal value φ max of the scalar field function φ on right middle panel, and the separation z d between the components of the k z = 1 pair in the lower panel. All quantities are shown versus the frequency ω.
Whereas the single BSs constitute an infinite set of branches, that form a spiral or a damped oscillation, depending on the quantity of interest [38], the even chains seem to end in a limiting configuration quite abruptly somewhere in the middle of a branch. The number of branches before this limiting configuration is encountered depends on the strength of the gravitational coupling α. For small α there are more branches, for larger α, the limiting configuration is already encountered on the second branch, as illustrated in Fig. 3 for the BS pair.
The mass and the charge exhibit only the onset of a spiraling behavior with two branches for the larger α values (α = 0.5, 1), and three branches for the smaller ones (α = 0.15, 0.2, 0.25). While the minimum f min the metric function (middle left) and the maximum φ max of the scalar function (middle right) exhibit only two or three oscillations. The coordinate distance z d between the two components of the pair, as given by twice the value of the z-coordinate of the maximum φ max , exhibits both types of behavior (lower). For small α (α = 0.15) z d shows three oscillations, while it decreases continuously. For larger α (α = 0.25) it exhibits the onset of a spiral, with again larger values on the third branch.
Let us now address the limiting behavior of the pairs in more detail, and its dependence on the coupling α. To that end, we illustrate in Fig. 4 the profiles of the metric function f (left panels) and the scalar field function φ (right panels) on the z-axis, choosing a large (upper panels) and a small (lower panels) value of the coupling α. Since for the large α, the coordinate distance z d does not decrease monotonically, but increases again towards the limiting solution, we retain two well-separated components in the limit. The minimum of the metric function f min tends to zero in the limit, while the maximum (minimum) of the scalar field function φ max (φ min ) becomes extremely sharp.
In fact, the scalar field amplitude φ acquires the shape of two antisymmetric peak(on)s associated with the locations of the minima of the metric function f , where the name peakons refers to field configurations with extremely large absolute values of the second derivative at the maxima [52,53]. Further numerical investigation of such singular solutions is, unfortunately, plagued by severe numerical problems. We remark, that in the case of odd chains the choice of spherical coordinates and the presence of the major peak(on) at the origin alleviate the numerical problems considerably.
For smaller α, the coordinate distance z d between the two components of the pair decreases monotonically towards the limiting solution. Fig. 4 shows that the closeness of the extrema of the scalar field function coincides with a very steep rise of the scalar field function at the origin. The metric function, however, retains only a small peak at the origin. Here the numerical grid allows for better resolution of this extremal behavior, and thus a closer approach to the limiting solution. Still, the approach towards the limiting solution is restricted by numerical accuracy.
The k z = 3 quartet represents a bound state of four BSs, located symmetrically along the symmetry axis. It may also be viewed as a bound state of two bound pairs of BSs, since this configuration is not fully symmetric. The two inner BSs are slightly larger than the two outer BSs, though there is not too much  distinction between the stars as long as the configuration resides on the fundamental branch. As seen in Fig. 3, these quartet configurations share most of the properties with k z = 1 pairs. For the chosen values of the coupling α we find two branches of solutions, the fundamental branch connected to the perturbative excitations at ω → ω max , and the second branch leading to a limiting solution. While approaching this limiting solution the outer extrema of the metric and of the scalar field function become less pronouced, leaving basically the two inner extrema, which evolve completely analogously to the extrema of the pair towards a limiting solution. We expect this scenario to represent a generic pattern seen for all chains with odd k z (although so far we have checked it for k z = 1, 3, 5 only).

Odd chains
As noticed above, the odd chains always possess a BS centered at the origin, with the remaining BSs are located symmetrically with respect to the origin, on the symmetry axis. The presence of the central BS constituent suggests that the (ω, M )-pattern of the odd chains could resemble that found for a single BS.
To scrutinize this conjecture, let us consider the branch structure for the k z = 2 triplet of BSs, exhibited in and φ max . The extremal values always reside at the center of the configurations, i.e., at the origin. Also the separation z d between the neighboring components of the k z = 2 triplet forms a spiral. Here, along this spiral, the contribution of the central BS to the full configurations becomes dominant, while the outer BSs contributions diminish.
When increasing the coupling constant beyond α 0.45, however, the scenario changes completely. While the configurations still follow a similar pattern along the fundamental branch, the further part of the (ω, M )-diagram does not involve a spiraling/oscillating behavior. Instead there is a single second branch, that leads all the way back to the vacuum configuration. Thus all physical quantities exhibit loops, beginning and ending at ω max . This may be interpreted as follows. Along the second branch the configurations again feature a dominant central BS, but the 'satellite' BSs dissolve into a (sort of) boson shell. Moreover, the system tends more and more towards spherical symmetry. Now we recall that a central BS surrounded by a boson shell is precisely what constitutes a radially excited spherically symmetric BS with a single radial node, n r = 1. This suggests that the k z = 2 triplet might merge with a n r = 1 single BS.
Let us therefore compare the M − ω diagram of the k z = 2 triplets at large coupling α with that of the radially excited n r = 1 single BSs. Such a comparison is shown in Fig. 5, where the k z = 2 BS triplets are marked by solid curves and the radially excited n r = 1 single BSs by dashed curves for two values of α. The upper panels show the scaled ADM mass M (left) and the scaled charge Q (right), while the lower panels show the minimal value f min of the metric function f (left) and the maximal value φ max of the scalar field function φ (right).
These figures are indeed very telling. The radially excited n r = 1 single BSs exhibit the typical curve pattern of spherically symmetric BSs. They emerge from the vacuum, form the fundamental branch and end in a spiraling/oscillating pattern. The k z = 2 triplet likewise emerges from the vacuum, forms a fundamental branch, and a second branch, but this second branch of the triplet nicely overlaps with the fundamental branch of the n r = 1 single BSs at some critical value of the frequency ω cr . The overlap happens when mass and charge of the n r = 1 single BSs have already passed their maximal values and the radially excited stars are descending into the spiral.
It is now clear how the domain of existence of odd chains with more constituents should be, for sufficiently large values of the coupling α. Let us consider a chain with k z nodes and thus k z + 1 constituents. Then this chain will feature on its fundamental branch k z + 1 more or less equal BSs, located on the symmetry axis. Subsequently the central BS will start to dominate while the satellite BSs will dissolve into boson shells. As the system tends towards spherical symmetry, it will overlap with a n r = k z /2 single BS at some critical value of the frequency. We have checked this behavior for the k z = 4 BS quintet, which indeed overlaps with the radially excited n r = 2 single BS.

Conclusions
The main purpose of this paper was to report the existence of a new type of solitonic configurations for a model with a gravitating self-interacting complex scalar field. These configurations represent chains of BSs, with k z + 1 constituents located symmetrically along the symmetry axis. The number k z 0 represents the number of nodes on the symmetry axis of the scalar field amplitude φ.
The chains emerge from the vacuum φ = 0 at a maximal value of the boson field frequency ω, which is given by the field's mass. For any k z , a fundamental branch of solutions is found emerging from this vacuum, in a (ω, M )-diagram (with M the ADM mass), ending at a minimal value of the frequency ω min , whose value is determined by the self-interaction potential and the gravitational coupling strength α. The subsequent solution curve depends on the number of constituents and the coupling α. A single spherical BS has been argued to form an infinite number of branches, leading to spirals or damped oscillations (depending on the quantities considered) as its limiting configuration is approached. For even chains we do not see such an endless spiraling/oscillating pattern. Instead we observe only two to three branches, depending on the coupling α.
As the limiting configuration is approached the even chains retain basically two of their k z +1 constituents, whose metric function g 00 exhibits two sharp peaks, reaching a very small value, while the scalar field features two sharp opposite extrema located right at the location of these peaks. The resulting configurations then feature huge second derivatives of the functions, which impede further numerical analysis towards the limiting solution.
The odd chains show a similar pattern as the single BS, when the coupling α is small. This may be interpreted as the central BS dominating the configurations on the higher branches. For larger α, however, the pattern changes totally, and the chains overlap on their second branch with a radially excited spherical single BS with n r = k z /2 (radial) nodes. In this case the central dominant BS will be surrounded by n r 'boson shells'. 8 These solutions can be generalized in various directions. The most obvious generalization is to include rotation. For the scalar field this means to include an explicit harmonic dependence on the azimuthal angle ϕ. The rotating single-BSs were obtained long ago [40,[56][57][58]. The rotating BSs with odd parity and k z = 1, representing the rotating generalizations of the pair of BSs, were also discussed in the literature [41]. We predict the existence of rotating generalizations for the triplet and the higher chains discussed in this work.
Single non-rotating BSs cannot be endowed with a black hole at the center; the no-hair theorems forbid their existence [59]. This result is, however, circumvented in the presence of spin, hairy generalizations of the Kerr black hole (BH) with a complex scalar field being reported in literature [44,45,[60][61][62][63][64][65]. These hairy BHs obey a synchronization condition relating the angular velocity of the event horizon and the boson field frequency. Most studies so far considered only an even parity scalar field, see e.g. [46,47,[66][67][68][69][70][71][72][73][74][75][76]. Synchronized hairy BHs with an odd parity scalar field were obtained in [46,47,76]. While these solutions represent only the simplest type of generalizations, containing a single black hole at the center of configurations with one (parity even) constituent or two (parity odd) constituents one can easily image to put a black hole either at the center of rotating configurations with more than two components, or to put a black hole at the center of each of the components along the symmetry axis. Such configurations should correspond to hairy double Kerr solutions or hairy multi-Kerr solutions. It would be interesting to see, whether the presence of the scalar hair can regularize such solutions, so that no conical singularity would be needed to hold them in equilibrium.
Along similar lines, but replacing rotation by an electric charge, one could investigate chains of electromagnetically charged BSs, generalizing the results for a single charged BSs in Einstein-Klein-Gordon-Maxwell model [54,78,79]. Some results in this direction were reported in the recent work [82], where (flat space) Q-chains were constructed in a model with a U (1) gauged scalar field, for a particular choice of the constants in the potential (2.21). Gravitating generalizations of these solutions should also exist, sharing some of the properties of the (ungauged) BS chains in this work. In this context, we mention the recent results in [79,80], showing that the no-hair theorem in [81] does not apply to a single charged static BS in a model with a Q-ball type potential (2.21), with the existence of BH generalizations. For chains of charged BSs, it would be again of particular interest to see whether they would support regular static multi-BH solutions.