CChains of Boson Stars
C. A. R. Herdeiro † , J. Kunz ‡ , I. Perapechka (cid:63) , E. Radu † and Y. Shnir (cid:5)† Department of Mathematics, University of Aveiro and CIDMA, Campus de Santiago, 3810-183 Aveiro, Portugal ‡ Institute of Physics, University of Oldenburg, Germany Oldenburg D-26111, Germany (cid:63)
Department of Theoretical Physics and Astrophysics, Belarusian State University, Minsk 220004, Belarus (cid:5)
BLTP, JINR, Dubna 141980, Moscow Region, Russia
January 2021
Abstract
We study axially symmetric multi-soliton solutions of a complex scalar field theory with a sexticpotential, minimally coupled to Einstein’s gravity. These solutions carry no angular momentum and canbe classified by the number of nodes of the scalar field, k z , along the symmetry axis; they are interpretedas 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 Noethercharge) in terms of their angular frequency, similarly to a single fundamental boson star, as long as thegravitational coupling is small; for larger coupling, however, the inner part of the spiral is replaced bya 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 alimiting solution with finite values of ADM mass and Noether charge. Many non-linear physical systems support non-topological solitons, which represent spatially localized fieldconfigurations. One of the simplest examples in flat space is given by Q -balls, which are particle-like con-figurations in a model with a complex scalar field possessing a harmonic time dependence and a suitableself-interaction potential [1–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 simplestsuch configurations are static and spherically symmetric, the scalar field possessing a mass term only, with-out 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 uncertaintyprinciple; 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) sym-metry. The charge Q is proportional to the angular frequency of the complex boson field and represents theboson particle number of the configurations [2, 3].In flat spacetime Q -balls exist only within a certain frequency range for the scalar field: between amaximal value ω max , which corresponds to the mass of the scalar excitations, and some lower non-zerocritical value ω min , that depends on the form of the potential. On the one hand, as the frequency ω isapproaching 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 thedomain. 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 , at1 a r X i v : . [ g r- q c ] J a n hich M, Q trivialise (rather than diverge). Then, the BSs exhibit a spiraling (or oscillating) pattern of thefrequency dependence of both charge and mass, where both tend to some finite limiting values at the centersof the spirals. Qualitatively, the appearance of the frequency-mass spiral may be related to oscillations in theforce-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 theway 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 correspondto chains of BSs. These are static, axially symmetric equilibrium configurations interpreted as a number ofBSs located symmetrically with respect to the origin, along the symmetry axis. We construct these solutionsand investigate their physical properties for a choice of the scalar field potential with quartic and sextic self-interaction terms, which was employed in most of the Q -balls literature. We note that similar configurationsof chains of constituents are known to exist both for gravitating and flat space non-Abelian monopoles anddyons [8–20], Skyrmions [21–23], electroweak sphalerons [24–28], SU (2) non-self dual configurations [29, 30]and Yang-Mills solitons in ADS 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). Configurationswith k z nodes on the symmetry axis possess k z + 1 constituents. Configurations with even and odd numbersof 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 thetypical spiraling picture and from the closed loop scenario. For chains with an even number of constituents thepattern always terminates at critical solutions. For chains with an odd number of constituents, the patterndepends on the strength of the gravitational interaction. The configurations then either merge with thecorresponding radial excitation of the fundamental solution, or the central constituent of the configurationsexhibits oscillations while retaining smaller satellite constituents.This paper is organized as follows. In Section II the theoretical setting is specified. This includes theaction, the equations of motion, the Ansatz and the boundary conditions for the BS chains. The numericalresults for these new equilibrium configurations are shown in Section III . We give our conclusions in SectionIV. 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 S = (cid:90) d x √− g (cid:20) R πG − g µν (cid:0) Φ ∗ , µ Φ , ν + Φ ∗ , ν Φ , µ (cid:1) − U ( | Φ | ) (cid:21) , (2.1)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 E µν ≡ R µν − g µν R − πG T µν = 0 , (2.2)where T µν ≡ Φ ∗ ,µ Φ ,ν + Φ ∗ ,ν Φ ,µ − g µν (cid:20) g στ (Φ ∗ ,σ Φ ,τ + Φ ∗ ,τ Φ ,σ ) + U ( | Φ | ) (cid:21) , (2.3)is the stress-energy tensor of the scalar field. Similar chains, but for a scalar field without self-interactions, were recently studied in [7], in the context of multipolar BSs.Here, we emphasise the interpretation of multiple BSs, rather than a single multipolar BS. (cid:18) (cid:3) − dUd | Φ | (cid:19) Φ = 0 , (2.4)where (cid:3) 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 Minkowskispacetime background. Their mass M can be obtained from the respective Komar expressions [33], M = 14 πG (cid:90) Σ R µν n µ ξ ν dV . (2.5)Here Σ denotes a spacelike hypersurface (with the volume element dV ), while n µ is normal to Σ, with n µ n µ = − M = 2 (cid:90) Σ (cid:18) T µν − g µν T γγ (cid:19) n µ ξ ν dV . (2.6)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 j µ = − i (Φ ∂ µ Φ ∗ − Φ ∗ ∂ µ Φ) . (2.7)It follows that integrating the timelike component of this 4-current in a spacelike slice Σ yields a secondconserved quantity – the Noether charge : Q = (cid:90) Σ j µ n µ dV. (2.8) Apart from being static, the configurations in this work are also axially symmetric . Thus, in a system ofadapted 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 areusually studied by employing a metric ansatz [34] ds = − e − U ( ρ,z ) dt + e U ( ρ,z ) (cid:104) e k ( ρ,z ) ( dρ + dz ) + P ( ρ, z ) dϕ (cid:105) , (2.10)where ( ρ, z ) correspond, asymptotically, to the usual cylindrical coordinates . However, it the numericaltreatment of the Einstein-Klein-Gordon equations, it is useful to employ ‘quasi-isotropic’ spherical coordi-nates ( r, θ ), defined by the coordinate transformation in (2.10) ρ = r sin θ, z = r cos θ , (2.11) The scalar field possesses a time dependence (with ∂ t Φ = − iω Φ), which disappears at the level of the energy-momentumtensor. However, the scalar field is axially symmetric, ∂ ϕ Φ = 0. In the Einstein-Maxwell theory, it is always possible to set P ≡ ρ , such that only two independent metric functions appearin the equations, with ( ρ, z ) the canonical Weyl coordinates [34]. For a (complex) scalar field matter content, however, thegeneric metric ansatz (2.10) with three independent functions is needed. (cid:54) r < ∞ , 0 (cid:54) θ ≤ π . Then the metric can then be written in a Lewis-Papapetrouform, with ds = − f dt + mf (cid:0) dr + r dθ (cid:1) + lf r sin θdϕ . (2.12)The three metric functions f , l , and m are functions of the variables r and θ only, chosen such that the trivialangular and radial dependence of the (asymptotically flat) line element is already factorized. The relationbetween the metric functions in the above line-element and those in (2.10) is f = e − U , lr sin θ = P , m = e k . 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 spacetimebackground 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 dependson ( r, θ ), Φ = φ ( r, θ ) e − iωt , (2.13)where ω ≥ T rr = fm (cid:18) φ ,r − φ ,θ r (cid:19) + ω φ f − U ( φ ) , T θθ = fm (cid:18) − φ ,r + φ ,θ r (cid:19) + ω φ f − U ( φ ) , T θr = 2 fr m φ ,r φ ,θ T ϕϕ = − fm (cid:18) φ ,r + φ ,θ r (cid:19) + ω φ f − U ( φ ) , T tt = − fm (cid:18) φ ,r + φ ,θ r (cid:19) − ω φ f − U ( φ ) . (2.14)After inserting the ansatz (2.12), (2.13) into the field equations (2.2), (2.4) we find a system of six coupledpartial 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 rr + E θθ = 0, E ϕϕ = 0and E tt = 0, yielding f ,rr + f ,θθ r + 2 f ,r r + cot θf ,θ r − f (cid:18) f ,r + f ,θ r (cid:19) + 12 l (cid:18) f ,r l ,r + f ,θ l ,θ r (cid:19) + 16 πG (cid:18) U ( φ ) − ω φ f (cid:19) m = 0 ,l ,rr + l ,θθ r + 3 l ,r r + 2 cot θl ,θ r − l (cid:18) l ,r + l ,θ r (cid:19) + 32 πG (cid:18) U ( φ ) − ω φ f (cid:19) lmf = 0 , (2.15) m ,rr + m ,θθ r + m ,r r + m f (cid:18) f ,r + f ,θ r (cid:19) − m (cid:18) m ,r + m ,θ r (cid:19) +16 πG (cid:34) fm (cid:18) φ ,r + φ ,θ r (cid:19) + U ( φ ) − ω φ f (cid:35) m f = 0 , and one equation for the scalar field amplitude φ ,rr + φ ,θθ r + (cid:18) r + l ,r l (cid:19) φ ,r + (cid:18) cot θ + l ,θ l (cid:19) φ ,θ r + ( ω f − ∂U∂φ ) mf φ = 0 . (2.16)Apart from these, there are two more Einstein equations, E rθ = 0 , E rr − E θθ = 0, which are not solved inpractice, 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 sin θ √ lmf . Inevaluating (2.5), we use the fact that R tt is a total derivative: √− gR tt = − ∂∂r (cid:32) r sin θ √ lf ,r f (cid:33) − ∂∂θ (cid:32) sin θ √ lf ,θ f (cid:33) . M can be read off from the asymptotic expansion of the metric function ff = 1 − M Gr + O (cid:18) r (cid:19) . (2.17)Alternatively, the mass M can be obtained by direct integration of (2.6) , M = (cid:90) (cid:0) T tt − T rr − T θθ − T ϕϕ (cid:1) √− g dr dθ dϕ = 4 π (cid:90) ∞ dr (cid:90) π dθ r sin θ √ lmf (cid:18) U ( φ ) − ω φ f (cid:19) . (2.18)In terms of the above Ansatz the Noether charge Q is given by Q = 4 πω ∞ (cid:90) dr π (cid:90) dθ √ lmf r sin θ φ . (2.19)Also, let us note that the solutions in this work have no horizon. Therefore they are zero entropyobjects, without an intrinsic temperature. Still, in analogy to black holes, one may write a “first law ofthermodynamics” [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 usedto check the numerical accuracy of the solutions. The solutions in this work were found for a potential originally proposed in [36, 37] , U ( | Φ | ) = ν | Φ | − λ | Φ | + µ | Φ | (2.21)which is chosen such that nontopological soliton solutions - Q -balls - exist in the absence of the Einsteinterm in the action (2.1), i.e. on a Minkowski spacetime background [38–41]. Also, at least in the sphericallysymmetric case, this choice of the potential follows for the existence of very massive and highly compactobjects 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 thefield | Φ | , besides the global minimum at | Φ | = 0.Two of the constants in (2.21) can actually be absorbed into a redefinition of the radial coordinatetogether with a rescaling of the scalar field, r → u µ r , φ → √ µν / √ u φ , (2.22)with u > ω → u µ ω. Then, the potential (2.21) becomes (up to an overall factor), U = φ − ¯ λφ + u φ , with ¯ λ = λu µ √ ν . (2.23)The choice employed in most of the Q -ball literature is u = 1 . , and ¯ λ = 2 , (2.24)which are the values used also in this work. 5or 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 thefollowing reduced action S red = (cid:90) drdθ (cid:0) L g − α L s (cid:1) , (2.25)with [for ease of notation, we denote ( ∇ S ) · ( ∇ T ) ≡ S ,r T ,r + r S ,θ T ,θ ]: L g = r sin θ √ l (cid:20) lm ( ∇ l ) · ( ∇ m ) − f ( ∇ f ) − rl (cid:18) l ,r + cot θr l ,θ (cid:19) + 1 rm (cid:18) m ,r + cot θr m ,θ (cid:19) (cid:21) , (2.26) L s = r sin θ m √ lf (cid:20) fm ( ∇ φ ) − ω φ f + φ − φ + 1 . φ (cid:21) . The (dimensionless) coupling constant reads α = 4 πGµ √ νu , (2.27)which determines the strength of the gravitational coupling of the solutions.As with the known spherically symmetric configurations, solutions exist for 0 (cid:54) α < ∞ . The limit α → i.e. Q -balls on a fixed Minkowski background. To understandthe large α limit, we define ˆ φ = φ/α . Then for large values of the effective gravitational coupling α , thenon-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 tothose found for a model with a non-self-interacting, massive complex scalar field ˆ φ . However, in this workwe shall restrict our study to the case of a finite, nonzero α . The basic properties of the mBS chains werestudied (in a more general context) in the recent work [7], while the issue of flat spacetime Q -ball chains willbe discussed elsewhere. The solutions studied in this work are globally regular and asymptotically flat, possessing finite mass andNoether 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 originrequires ∂ r f (cid:12)(cid:12) r =0 = ∂ r m (cid:12)(cid:12) r =0 = ∂ r l (cid:12)(cid:12) r =0 = 0 . (2.28)Demanding asymptotic flatness at spatial infinity yields f (cid:12)(cid:12) r →∞ = m (cid:12)(cid:12) r →∞ = l (cid:12)(cid:12) r →∞ = 1 . (2.29)Axial symmetry and regularity impose on the symmetry axis at θ = 0 , π the conditions ∂ θ φ (cid:12)(cid:12) θ =0 ,π = ∂ θ f (cid:12)(cid:12) θ =0 ,π = ∂ θ m (cid:12)(cid:12) θ =0 ,π = ∂ θ l (cid:12)(cid:12) θ =0 ,π = 0 . (2.30)Further, the condition of the absence of a conical singularity requires that the solutions should satisfy theconstraint m (cid:12)(cid:12) θ =0 ,π = l (cid:12)(cid:12) θ =0 ,π . In our numerical scheme we explicitly verified (within the numerical accuracy)this condition on the symmetry axis. In (2.26) we use the dimensionless radial variable r and the dimensionless frequency ω given in units set by µ . φ approachesasymptotically the global minimum, φ (cid:12)(cid:12) r →∞ = 0 , (2.31)while on the symmetry axis we impose ∂ θ φ (cid:12)(cid:12) θ =0 ,π = 0 . (2.32)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, f ( r, π − θ ) = f ( r, θ ) , l ( r, π − θ ) = l ( r, θ ) , m ( r, π − θ ) = m ( r, θ ) ; (2.33)for the scalar field, however, there are two possibilities: P = 1 (even parity) : φ ( r, π − θ ) = φ ( r, θ ) , (2.34) P = − φ ( r, π − θ ) = − φ ( r, θ ) . (2.35)We use this symmetry to reduce the domain of integration to [0 , ∞ ) × [0 , π/ ∂ θ f (cid:12)(cid:12) θ = π/ = ∂ θ m (cid:12)(cid:12) θ = π/ = ∂ θ l (cid:12)(cid:12) θ = π/ = 0 , (2.36)and P = 1 : ∂ θ φ (cid:12)(cid:12) θ = π/ = 0 , P = − φ (cid:12)(cid:12) θ = π/ = 0 , (2.37)together with P = 1 : ∂ r φ (cid:12)(cid:12) r =0 = 0 , P = − φ (cid:12)(cid:12) r =0 = 0 . (2.38)Finally, let us mention that, for both P = ±
1, one can formally construct an approximate expression ofthe solutions compatible with the above boundary conditions, e.g. by assuming the existence of a powerseries form close to r = 0. The resulting expressions are of little help in understanding the properties ofthe solutions, and a numerical approach is necessary . However, one important result of this study is thebound-state condition ω (cid:54) µ, (2.39)which emerges from the finite mass requirement. No similar result is found for the minimal value of frequency. The set of four coupled non-linear elliptic partial differential equations for the functions ( f, l, m ; φ ) has beensolved numerically subject to the boundary conditions defined above. In a first stage, a new compactifiedradial variable x is introduced, replacing r , with x ≡ rc + r , (2.40)with c an arbitrary parameter, typically of order one. With this choice, the semi-infinite region [0 , ∞ ) ismapped to the finite interval [0 , ×
89. The underlyinglinear system is solved with the Intel MKL PARDISO sparse direct solver [50] using the Newton-Raphson We recall that, more than fifty years after their discovery [4, 5], the static mBS are still not known in closed form. library. Inall cases, the typical errors are of order of 10 − .For the choice (2.24) of the potential’s parameters, the only input parameters are the gravitationalcoupling constant α together with the scalar field’s frequency ω . The parity of the solutions is imposed viathe boundary conditions (2.34), (2.38). The number of individual constituents results from the numericaloutput. Also, in all plots below, quantities are given in natural units set by µ, G . 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 thesymmetry axis by k z , the solutions can be classified by the number of nodes k z . The number of constituentsof the chains is given by k z + 1. The even-parity configurations ( P = 1) have an even k z , while the solutionswith an odd k z are found for P = −
1. For example, the spherically symmetric solutions have k z = 0 and onesingle constituent localized at the origin, r = 0. The simplest non-spherical configuration has k z = 1 andrepresents a pair of static BSs 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 dependson their relative phase [42, 43]. If the Q-balls are in phase, the interaction is attractive, if they are out ofphase, there is a repulsive force between them. Thus, an axially symmetric k z = 1 solution can be balancedby the gravitational attraction.Solutions with k z > k z = 5; however, theyare 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 repre-sentative 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 . ω, M )-diagram (as described below). These chains possess one tofive constituents (from top to bottom), as seen by the number of peaks of the charge density, shown in theleft panels. The middle panels represent the scalar field amplitude φ , and the right panels show the metricfunction 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 pairof BSs. The charge density has two symmetric peaks, the metric function has two symmetric troughs, whilethe scalar field function is anti-symmetric, featuring a peak and a trough. The triplet, quartet and quintetin the next few rows feature k z + 1 very similar (in size and shape) peaks for the charge density and troughsfor the metric function, while the scalar field shows alternating peaks and troughs, all located symmetricallyalong 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 isseen 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 − α = 0 . k z = 4, and ω/µ = { . , . , . , . , . } for k z = 0 −
4, respectively. As welook at the charge density of the configurations, we see a dominant peak at the center for odd chains and adominant inner pair for even chains, while the other peaks of the triplet, quartet and quintet have turnedinto slight elevations, hardly visible in the figures. In the metric functions the troughs at the center of the Complex Equations – Simple Domain partial differential equations SOLver is a C++ package being developed by one of us(I.P.). In principle, the k z = 1 solutions can be thought of as the static limit of negative parity spinning configurations consideredin [41]. Chains of BSs with one to five constituents (from top to bottom) on the first (a.k.a. fundamental ) branchfor α = 0 .
25 at ω/µ = 0 .
80: 3 d plots of the U (1) scalar charge distributions (left plots), the scalar field functions φ (middle plots) and the metric functions f (right plots) versus the coordinates ρ = r sin θ and z = r cos θ . odd chains dominate, and likewise the (almost merged) inner pair of troughs of the even chains. All othertroughs are weakened substantially. The scalar field itself, however, retains the outer peaks and troughs toa somewhat greater extent, still reflecting clearly the number of constituents of the chains. ω -dependence and the branch structure Recall the frequency dependence for the single BSs and for fixed coupling α [38]. The set of BSs emergesfrom 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 thefrequency ω spans the first or fundamental branch, which terminates at the first backbending of the curve,9igure 2: Chains of BSs with one to five constituents (from top to bottom) on the second branch for differing valuesof α and ω : 3 d plots of the U (1) scalar charge distributions (left plots), the scalar field functions φ (middle plots)and the metric functions f (right plots) versus the coordinates ρ = r sin θ and z = r cos θ . 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 andthe maximum of the scalar function φ shows damped oscillations. The set of solutions tends to a limitingsolution at the center of the spiral which has finite values of the mass and charge. However, the values of thescalar field function φ and the metric function f at the center of the star, which represent the maximal andminimal values of these functions, φ max and f min , respectively, do not seem to be finite, with φ max divergingand f min vanishing in the limit.Let us now consider the frequency dependence of the BS chains, when the coupling α is kept fixed. Likethe single BSs, all chains emerge similarly from the vacuum at the maximal frequency, given by the boson10ass ω 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 forthe single BS, for all chains this fundamental branch ends at a minimal value of the frequency, from wherea second branch arises. But then even and odd chains will in general exhibit different patterns, which willalso depend on the coupling strength α . Let us consider the BS pair and the higher even chains in more detail. We illustrate the ω -dependence foreven chains in Fig. 3, selecting the BS pair ( k z = 1) and the BS quartet ( k z = 3), and comparing withthe single BS. In the upper panels we show the k z = 1 pair (solid curves), the k z = 3 quartet (dash-dottedcurves), 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). Thedifferent colors refer to different values of the coupling α . The middle left panel shows in an analogousmanner the minimal value f min of the metric function f . Restricting to the k z = 1 pair, we then exhibit themaximal value φ max of the scalar field function φ on right middle panel, and the separation z d between thecomponents 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 quiteabruptly somewhere in the middle of a branch. The number of branches before this limiting configuration isencountered 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. 3for the BS pair.The mass and the charge exhibit only the onset of a spiraling behavior with two branches for the larger α values ( α = 0 . , α = 0 . , . , . f min the metric function (middle left) and the maximum φ max of the scalar function (middle right) exhibit onlytwo or three oscillations. The coordinate distance z d between the two components of the pair, as given bytwice the value of the z -coordinate of the maximum φ max , exhibits both types of behavior (lower). For small α ( α = 0 . z d shows three oscillations, while it decreases continuously. For larger α ( α = 0 .
25) it exhibitsthe 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 fieldfunction φ (right panels) on the z -axis, choosing a large (upper panels) and a small (lower panels) valueof the coupling α . Since for the large α , the coordinate distance z d does not decrease monotonically, butincreases again towards the limiting solution, we retain two well-separated components in the limit. Theminimum of the metric function f min tends to zero in the limit, while the maximum (minimum) of the scalarfield function φ max ( φ min ) becomes extremely sharp.In fact, the scalar field amplitude φ acquires the shape of two antisymmetric peak(on)s associated withthe locations of the minima of the metric function f , where the name peakons refers to field configurationswith extremely large absolute values of the second derivative at the maxima [52, 53]. Further numericalinvestigation 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) atthe origin alleviate the numerical problems considerably.For smaller α , the coordinate distance z d between the two components of the pair decreases monotonicallytowards the limiting solution. Fig. 4 shows that the closeness of the extrema of the scalar field functioncoincides 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 extremalbehavior, and thus a closer approach to the limiting solution. Still, the approach towards the limiting solutionis restricted by numerical accuracy.The k z = 3 quartet represents a bound state of four BSs, located symmetrically along the symmetryaxis. It may also be viewed as a bound state of two bound pairs of BSs, since this configuration is not fullysymmetric. The two inner BSs are slightly larger than the two outer BSs, though there is not too much11 .4 0.5 0.6 0.7 0.8 0.9 1.00510152025 =0.25 =0.5 =1 M Q f m i n m a x z d =0.25 =0.2 =0.15 Figure 3:
Comparison of the k z = 0 single BSs (dashed curves), the k z = 1 pair of BSs (solid curves), and the k z = 3quartet of BSs (dash-dotted curves): scaled ADM mass M (upper left panel), scaled charge Q (upper right panel),minimal value of the metric function f min (middle left panel), maximal value of the scalar field φ max of the k z = 1pair only (middle right panel), and separation z d between the two components of the k z = 1 pair only (lower panel) vs. frequency ω for several values of the coupling α . Note the quadratic scale for f min and z d . distinction between the stars as long as the configuration resides on the fundamental branch. As seen inFig. 3, these quartet configurations share most of the properties with k z = 1 pairs. For the chosen valuesof the coupling α we find two branches of solutions, the fundamental branch connected to the perturbative12igure 4: Profiles on the symmetry axis of the metric function f (left) and the scalar field function φ of the k z = 1(almost) critical solution on the second branch at ω = 0 .
85 and α = 1 (upper plots), and on the third branch at ω = 0 .
52 and α = 0 . z = r cos θ , with θ = 0 , π . excitations at ω → ω max , and the second branch leading to a limiting solution. While approaching thislimiting 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 pairtowards a limiting solution. We expect this scenario to represent a generic pattern seen for all chains withodd k z (although so far we have checked it for k z = 1 , , As noticed above, the odd chains always possess a BS centered at the origin, with the remaining BSs arelocated symmetrically with respect to the origin, on the symmetry axis. The presence of the central BSconstituent 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 inFig. 5 for four values of the coupling α (with different colors). Analogously to the even chains, we show inthe two upper panels the scaled ADM mass M (left) and the scaled charge Q (right). The middle panelsshow the minimal value f min of the metric function (left) and the maximal value φ max of the scalar fieldfunction (right), while the lower panel shows the separation z d between the neighboring components of the k z = 2 triplet. Again, all quantities are shown as a function of the frequency ω .Let us first consider small values of the coupling α . Then our expectation is borne out, and we observethe triplet forming spirals for the mass and charge, while damped oscillations for the extremal values f min .3 0.4 0.5 0.6 0.7 0.8 0.9 1.00510152025303540 =0.15 =0.25 =0.5 M =1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.001020304050607080 =1 =0.5 =0.25 =0.15 Q f m i n =1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.00.00.51.01.52.02.53.03.5 =0.15 =0.25=0.5 ( ) =10.3 0.4 0.5 0.6 0.7 0.8 0.9 1.00510152025303540 =0.15 =0.25 =0.5 z d =1 Figure 5: k z = 2 triplet of BSs: scaled ADM mass M (upper left panel), scaled charge Q (upper right panel),minimal value of the metric function f min (middle left panel), maximal value of the scalar field function φ max (middleright panel), and separation z d between neighboring components of the k z = 2 triplet (lower panel) vs. frequency ω for several values of the coupling α . Note the quadratic scale for Q , f min and z d . φ max . The extremal values always reside at the center of the configurations, i.e. , at the origin. Alsothe separation z d between the neighboring components of the k z = 2 triplet forms a spiral. Here, along thisspiral, the contribution of the central BS to the full configurations becomes dominant, while the outer BSscontributions diminish.When increasing the coupling constant beyond α (cid:38) .
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, beginningand ending at ω max . This may be interpreted as follows. Along the second branch the configurations againfeature a dominant central BS, but the ’satellite’ BSs dissolve into a (sort of) boson shell. Moreover, thesystem tends more and more towards spherical symmetry. Now we recall that a central BS surrounded bya boson shell is precisely what constitutes a radially excited spherically symmetric BS with a single radialnode, 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 theradially excited n r = 1 single BSs. Such a comparison is shown in Fig. 5, where the k z = 2 BS triplets aremarked by solid curves and the radially excited n r = 1 single BSs by dashed curves for two values of α . Theupper panels show the scaled ADM mass M (left) and the scaled charge Q (right), while the lower panelsshow the minimal value f min of the metric function f (left) and the maximal value φ max of the scalar fieldfunction φ (right).These figures are indeed very telling. The radially excited n r = 1 single BSs exhibit the typical curvepattern of spherically symmetric BSs. They emerge from the vacuum, form the fundamental branch and endin a spiraling/oscillating pattern. The k z = 2 triplet likewise emerges from the vacuum, forms a fundamentalbranch, and a second branch, but this second branch of the triplet nicely overlaps with the fundamentalbranch of the n r = 1 single BSs at some critical value of the frequency ω cr . The overlap happens when massand charge of the n r = 1 single BSs have already passed their maximal values and the radially excited starsare descending into the spiral.It is now clear how the domain of existence of odd chains with more constituents should be, for sufficientlylarge values of the coupling α . Let us consider a chain with k z nodes and thus k z + 1 constituents. Then thischain 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. Asthe system tends towards spherical symmetry, it will overlap with a n r = k z / k z = 4 BS quintet, which indeed overlaps withthe radially excited n r = 2 single BS. The main purpose of this paper was to report the existence of a new type of solitonic configurations for amodel 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 (cid:62) φ .The chains emerge from the vacuum φ = 0 at a maximal value of the boson field frequency ω , whichis given by the field’s mass. For any k z , a fundamental branch of solutions is found emerging from thisvacuum, 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 α . Thesubsequent solution curve depends on the number of constituents and the coupling α . A single spherical BShas been argued to form an infinite number of branches, leading to spirals or damped oscillations (dependingon the quantities considered) as its limiting configuration is approached. For even chains we do not see suchan endless spiraling/oscillating pattern. Instead we observe only two to three branches, depending on thecoupling α .As the limiting configuration is approached the even chains retain basically two of their k z +1 constituents,whose metric function g exhibits two sharp peaks, reaching a very small value, while the scalar field features15 .70 0.75 0.80 0.85 0.90 0.95 1.000246810121416 =0.5 M =1 0.70 0.75 0.80 0.85 0.90 0.95 1.000246810121416 =0.5 Q =10.70 0.75 0.80 0.85 0.90 0.95 1.0000.10.250.51 =0.5 f m i n =1 0.70 0.75 0.80 0.85 0.90 0.95 1.000.00.51.01.52.02.5 =0.5 ( ) =1 Figure 6:
Comparison of the radially excited n r = 1 single BSs (dashed curves), and k z = 2 triplet of BSs (solidcurves): scaled ADM mass M (upper left panel), scaled charge Q (upper right panel), minimal value of the metricfunction f min (lower left panel), maximal value of the scalar field function φ max (lower right panel) vs. frequency ω for two values of the coupling α . Note the quadratic scale for f min . two sharp opposite extrema located right at the location of these peaks. The resulting configurations thenfeature huge second derivatives of the functions, which impede further numerical analysis towards the limitingsolution.The odd chains show a similar pattern as the single BS, when the coupling α is small. This may beinterpreted 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 sphericalsingle BS with n r = k z / n r ‘boson shells’. These solutions can be generalized in various directions. The most obvious generalization is to includerotation. 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–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]. Wepredict 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 forbidtheir existence [59]. This result is, however, circumvented in the presence of spin, hairy generalizations ofthe Kerr black hole (BH) with a complex scalar field being reported in literature [44, 45, 60–65]. These hairy Boson shells without a central BS are also known [54, 55]. e.g. [46,47,66–76]. Synchronizedhairy BHs with an odd parity scalar field were obtained in [46, 47, 76]. While these solutions represent onlythe 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 atthe center of rotating configurations with more than two components, or to put a black hole at the centerof each of the components along the symmetry axis. Such configurations should correspond to hairy doubleKerr solutions or hairy multi-Kerr solutions. It would be interesting to see, whether the presence of thescalar hair can regularize such solutions, so that no conical singularity would be needed to hold them inequilibrium.Along similar lines, but replacing rotation by an electric charge, one could investigate chains of electro-magnetically charged BSs, generalizing the results for a single charged BSs in Einstein-Klein-Gordon-Maxwellmodel [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 constantsin the potential (2.21). Gravitating generalizations of these solutions should also exist, sharing some of theproperties 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 wouldbe again of particular interest to see whether they would support regular static multi-BH solutions. Acknowledgements
We gratefully acknowledge support by the DFG funded Research Training Group 1620 “Models of Gravity”and by the DAAD. We would also like to acknowledge networking support by the COST Actions CA16104and CA15117. Ya.S. gratefully acknowledges the support by the Ministry of Sceince and High Education ofRussian Federation, project FEWF-2020-0003. and by the BLTP JINR Heisenberg-Landau program 2020.This work is also supported by the Center for Research and Development in Mathematics and Applications(CIDMA) through the Portuguese Foundation for Science and Technology (FCT - Fundacao para a Ciˆenciae a Tecnologia), references UIDB/04106/2020 and UIDP/04106/2020 and by national funds (OE), throughFCT, I.P., in the scope of the framework contract foreseen in the numbers 4, 5 and 6 of the article 23, of theDecree-Law 57/2016, of August 29, changed by Law 57/2017, of July 19. We acknowledge support from theprojects PTDC/FIS-OUT/28407/2017, CERN/FIS-PAR/0027/2019 and PTDC/FIS-AST/3041/2020. Thiswork has further been supported by the European Union’s Horizon 2020 research and innovation (RISE)programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. The authors would like to acknowledgenetworking support by the COST Action CA16104.
References [1] G. Rosen, J. Math. Phys. (1968) 996, 999[2] R. Friedberg, T. D. Lee and A. Sirlin, Phys. Rev. D (1976) 2739[3] S. R. Coleman, Nucl. Phys. B (1985) 263 Erratum: [Nucl. Phys. B (1986) 744].[4] D. J. Kaup, Phys. Rev. (1968) 1331.[5] R. Ruffini and S. Bonazzola, Phys. Rev. (1969) 1767.[6] L. G. Collodel, B. Kleihaus and J. Kunz, Phys. Rev. D (2017) no.8, 084066[7] C. A. R. Herdeiro, J. Kunz, I. Perapechka, E. Radu and Y. Shnir, Phys. Lett. B (2021), 136027[arXiv:2008.10608 [gr-qc]]. 178] B. Kleihaus and J. Kunz, Phys. Rev. D (2000) 025003[9] B. Kleihaus and J. Kunz, Phys. Rev. Lett. (2000), 2430-2433[10] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Lett. B (2003) 237[11] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Rev. D (2003) 101701[12] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Rev. D (2004) 065010[13] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Rev. D (2005) 024013[14] R. Teh and K. Wong, J. Math. Phys. (2005), 082301[15] V. Paturyan, E. Radu and D. Tchrakian, Phys. Lett. B (2005), 360-366[16] B. Kleihaus, J. Kunz and U. Neemann, Phys. Lett. B (2005) 171[17] J. Kunz, U. Neemann and Y. Shnir, Phys. Lett. B (2006) 57[18] J. Kunz, U. Neemann and Y. Shnir, Phys. Rev. D (2007) 125008[19] K. Lim, R. Teh and K. Wong, J. Phys. G (2012), 025002[20] R. Teh, A. Soltanian and K. Wong, Phys. Rev. D (2014) no.4, 045018[21] S. Krusch and P. Sutcliffe, J. Phys. A (2004) 9037[22] Y. Shnir and D. H. Tchrakian, J. Phys. A (2010) 025401[23] Y. Shnir, Phys. Rev. D (2015) no.8, 085039[24] B. Kleihaus, J. Kunz and M. Leissner, Phys. Lett. B (2008), 438-444[25] R. Ibadov, B. Kleihaus, J. Kunz and M. Leissner, Phys. Lett. B (2008), 136-140[26] R. Ibadov, B. Kleihaus, J. Kunz and M. Leissner, Phys. Lett. B (2010), 298-306[27] R. Ibadov, B. Kleihaus, J. Kunz and M. Leissner, Phys. Rev. D (2010), 125037[28] R. Teh, B. Ng and K. Wong, Annals Phys. (2015), 170-195[29] E. Radu and D. H. Tchrakian, Phys. Lett. B (2006) 201[30] Y. M. Shnir, EPL (2007) no.2, 21001.[31] O. Kichakova, J. Kunz, E. Radu and Y. Shnir, Phys. Rev. D (2012) 104065.[32] O. Kichakova, J. Kunz, E. Radu and Y. Shnir, Phys. Rev. D (2014) no.12, 124012.[33] R. M. Wald, General Relativity , (University of Chicago Press, Chicago, 1984).[34] D. Kramer, H. Stephani, E. Herlt, and M. MacCallum,
Exact Solutions of Einstein’s Field Equations,
Cambridge University Press, Cambridge, (1980).[35] T. D. Lee and Y. Pang, Phys. Rept. (1992) 251.[36] W. Deppert and E. W. Mielke, Phys. Rev. D (1979) 1303.[37] E. W. Mielke and R. Scherzer, Phys. Rev. D , 2111 (1981).[38] R. Friedberg, T. Lee and Y. Pang, Phys. Rev. D (1987), 36581839] M.S. Volkov and E. Wohnert, Phys. Rev. D (2002) 085003.[40] B. Kleihaus, J. Kunz and M. List, Phys. Rev. D (2005) 064002 .[41] B. Kleihaus, J. Kunz, M. List and I. Schaffer, Phys. Rev. D (2008) 064025[42] R. Battye and P. Sutcliffe, Nucl. Phys. B , (2000) 329[43] P. Bowcock, D. Foster and P. Sutcliffe, J. Phys. A (2009), 085403[44] S. Hod, Phys. Rev. D (2012) 104026 Erratum: [Phys. Rev. D (2012) 129902][45] C. Herdeiro, E. Radu and H. Runarsson, Phys. Lett. B (2014) 302[46] J. Kunz, I. Perapechka and Y. Shnir, Phys. Rev. D (2019) no.6, 064032[47] Y. Q. Wang, Y. X. Liu and S. W. Wei, Phys. Rev. D (2019) no.6, 064036[48] Y. Brihaye and B. Hartmann, Phys. Rev. D (2009) 064013[49] A. Bernal, J. Barranco, D. Alic and C. Palenzuela, Phys. Rev. D (2010) 044031[50] N.I.M. Gould, J.A. Scott and Y. Hu, ACM Transactions on Mathematical Software (2007) 10;O. Schenk and K. G¨artner Future Generation Computer Systems (3) (2004) 475.[51] W. Sch¨onauer and R. Weiß, J. Comput. Appl. Math. 27, 279 (1989) 279;M. Schauder, R. Weiß and W. Sch¨onauer, Universit¨at Karlsruhe, Interner Bericht Nr. 46/92 (1992).[52] R. Camassa and D. D. Holm, Phys. Rev. Lett. (1993), 1661[53] J. Lis, Phys. Rev. D (2011), 085030[54] B. Kleihaus, J. Kunz, C. Lammerzahl and M. List, Phys. Lett. B (2009), 102-115[55] B. Kleihaus, J. Kunz, C. Lammerzahl and M. List, Phys. Rev. D (2010), 104050[56] F. E. Schunck and E. W. Mielke, Phys. Lett. A , 389 (1998).[57] S. Yoshida and Y. Eriguchi, Phys. Rev. D , 762 (1997).[58] F. D. Ryan, Phys. Rev. D , 6081 (1997).[59] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D (2015) no.09, 1542014[60] C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. (2014) 221101[61] C. Herdeiro and E. Radu, Phys. Rev. D (2014) no.12, 124018[62] C. Herdeiro and E. Radu, Class. Quant. Grav. (2015) no.14, 144001[63] S. Hod, Phys. Rev. D (2014) no.2, 024051[64] C. L. Benone, L. C. B. Crispino, C. Herdeiro and E. Radu, Phys. Rev. D , no. 10, 104024 (2014)[65] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D (2014) no.12, 1442014[66] B. Kleihaus, J. Kunz and S. Yazadjiev, Phys. Lett. B (2015) 406[67] C. A. R. Herdeiro, E. Radu and H. Runarsson, Phys. Rev. D (2015) no.8, 084059[68] C. Herdeiro, J. Kunz, E. Radu and B. Subagyo, Phys. Lett. B , 30 (2015)[69] C. Herdeiro, E. Radu and H. Runarsson, Class. Quant. Grav. (2016) no.15, 1540011970] Y. Brihaye, C. Herdeiro and E. Radu, Phys. Lett. B (2016) 279[71] S. Hod, Phys. Lett. B (2015) 177[72] C. Herdeiro, J. Kunz, E. Radu and B. Subagyo, Phys. Lett. B , 151 (2018)[73] C. Herdeiro, I. Perapechka, E. Radu and Y. Shnir, JHEP (2018) 119[74] C. Herdeiro, I. Perapechka, E. Radu and Y. Shnir, JHEP (2019), 111[75] J. F. M. Delgado, C. A. R. Herdeiro and E. Radu, Phys. Lett. B (2019), 436-444[76] J. Kunz, I. Perapechka and Y. Shnir, JHEP (2019), 109[77] B. Kleihaus and J. Kunz, Phys. Lett. B (2000), 130-134[78] P. Jetzer and J. van der Bij, Phys. Lett. B (1989), 341-346[79] C. A. R. Herdeiro and E. Radu, Eur. Phys. J. C (2020) no.5, 390 [arXiv:2004.00336 [gr-qc]].[80] J. P. Hong, M. Suzuki and M. Yamada, Phys. Rev. Lett. (2020) no.11, 111104 [arXiv:2004.03148[gr-qc]].[81] A. E. Mayo and J. D. Bekenstein, Phys. Rev. D54