Boson stars in massless and massive scalar-tensor gravity

We study phenomenological features and stability of boson stars in massless and massive scalar-tensor theory of gravity with Damour-Esposito-Farese coupling. This coupling between the tensor and scalar sectors of the theory leads to a phenomenon called spontaneous scalarization, the onset of which we investigate by numerically computing families of boson-star models using shooting and relaxation algorithms. We systematically explore the effects of the theory's coupling, the mass of the gravitational scalar and the choice of the bosonic potential on the structure of weakly and strongly scalarized solutions. Scalarized boson-star models share many common features with neutron stars in the same scalar-tensor theory of gravity. In particular, scalarization can result in boson stars with significantly larger radii and masses, which tend to be energetically favored over their weakly or non-scalarized counterparts. Overall, we find that boson stars are not quite as susceptible to scalarization as neutron stars.

ST theories of gravity are one of the most natural and simplest extensions of GR that preserve the universality of a free fall [49].The key ingredient of ST gravity is the presence of a scalar field non-minimally coupled to the metric such that gravity is not only mediated by a spin-2 graviton but also by the spin-0 scalar field.The inclusion of the scalar field is a common by-product of the compactification of higher dimensional Kaluza-Klein type theories, like string theory, [50][51][52], making ST theories well-motivated gravity theories to study.
The theoretical implications of ST theories of gravity in various areas of fundamental physics are vast.For instance, on cosmological scales, ST theories have been proposed as models for dark energy, replacing the role of the cosmological constant [53,54].Furthermore, the additional scalar degree of freedom introduces a new channel of radiation emission in the form of scalar gravitational waves (GWs) also known as a "breathing mode".The perfect laboratories for tests and searches of these scalar waves are compact objects and binary systems composed thereof.In particular, it has been conjectured that in the process of spontaneous scalarization, the scalar radiation can be detected with future gravitational wave detectors [55].The phenomenon of spontaneous scalarization is caused by a linear tachyonic instability inside the star that triggers the growth of the scalar field around it.Once the instability is quenched, the star becomes endowed with "scalar hair"; see Ref. [56] for a comprehensive review.This remarkable feature was first identified in the pioneering work of Damour and Esposito-Farèse [57][58][59], who studied spontaneous scalarization of neutron stars in massless ST theory with a pair of parameters (α 0 , β 0 ), inspired by Post-Newtonian (PN) theory and controlling the onset of scalarization.
However, stringent constraints have been placed on massless ST theories by numerous observations, including the Cassini bound α 0 < 3.4 × 10 −3 [60] and pulsar-white dwarf binaries which rule out models with β 0 ≲ −5 [61,62]; for a summary of constraints on a broad range of modified theories of gravity obtained from gravitationalwave (GW) observations see also [63].As a consequence, massive ST theories have acquired a great deal of attention, mainly for two reasons.First, the screening of the massive field maintains their compatibility with current observations and leaves a significant portion of the parameter space unconstrained.Second, because the mass term endows the scalar GW signal with a highly conspicuous long-lived, inverse chirp character [64][65][66][67][68][69][70][71][72][73][74]; see also Refs.[75][76][77] for studies of massive ST theory with self-interacting fields.
So far NSs have been extensively explored in this class of ST theories as ideal compact matter objects, in whose environment scalarization can occur in this class of ST theories [78][79][80].An intriguing extension is to investigate how the presence of a non-minimally coupled scalar field affects the phenomenology of other compact matter objects, such as BSs.In this work, we tackle this question by exploring the scalarization of ground-state equilibrium BSs for both massless and massive gravitational scalar fields as well as the stars' stability for various BS potentials.
This paper is structured as follows.In Sec.II, we present an overview of ST gravity with the specific Damour-Esposito-Farèse coupling function employed in our work, and formulate the system of first-order differential equations with the bosonic field in the matter sector.In Sec.III, we present the results of our exploration of the parameter space and describe the main features of BS models in ST theory of gravity.We span the solution space for different parameters of the coupling function, mass values of the gravitational scalar field and BS potentials.In Sec.IV, we assess the stability of these models.Finally, we summarize our results together with an outlook on future research directions in Sec.V.In Appendix A, we study in detail the asymptotic behavior of BS solutions in ST theory of gravity and GR, highlighting the main inconsistencies that can arise in a flat-field treatment.Our numerical methods are described in Appendix B and in Appendix C, we illustrate the structure of BS families of solutions for a broad range of BS potentials.Unless specified otherwise, we employ units where the speed of light and Planck's constant satisfy c = 1 = ℏ, so that the Planck mass is given by M Pl = 1/ √ G.

II. THEORY AND FORMULATION A. Scalar-tensor action
We start by considering a general class of ST theories of gravity with a single, real-valued gravitational scalar field φ with potential W (φ), whose action in the Einstein frame is given by [81], Here, ḡµν is a conformal metric, with corresponding Ricci scalar R, related to the physical or Jordan-frame metric where F (φ) is the coupling function, which needs to be positive to ensure that the graviton carries positive energy [57].For the matter sources, represented in Eq. ( 1) by the term S M , we consider a relativistic boson condensate, modelled as a massive complex scalar field ψ = ψ R + iψ I with complex conjugate ψ * and potential V (ψ).Note that by Eq. ( 1) the matter moves according to the geometry of the physical metric g µν ; this is the price tag for achieving the minimal coupling between the gravitational scalar and the conformal metric ḡµν of the Einstein frame.Although the scalar field does not directly interact with the ordinary matter (which guarantees the weak equivalence principle), it affects the motion of test particles through its presence in Eq. ( 2).The appearance of two distinct scalar fields in our physical setup inevitably holds potential for confusion; whenever the meaning is not self-evident, we will therefore explicitly refer to the two fields as the gravitational scalar φ and the BS scalar ψ with amplitude A . .= |ψ|.
The covariant field equations are obtained from varying the action (1) with respect to ḡµν , φ and ψ.This leads to where an overbar denotes tensors and operators in the Einstein frame; in particular, the energy-momentum tensor is given by

B. Spherically symmetric solutions
In this work, we consider stationary, spherically symmetric BS models and employ a metric ansatz using ra-dial gauge and polar slicing.Specifically, we write where r is the areal radius in the Einstein frame and dΩ 2 . .= dϑ 2 + sin 2 ϑdϕ 2 the standard angular line element.In spherical symmetry, it is convenient to express the complex BS scalar field in terms of amplitude A and frequency ω, At this point, we have three mass scales in our problem: the two scalar masses, µ φ and µ for φ and ψ respectively, and the Planck mass M pl = 1/G that appears in the Einstein equation (3).All dimensional quantities have units of some power of µ or M Pl or a combination thereof.To see this explicitly, we still need to specify the potential functions V (ψ) and W (φ). The particular choices for these functions will be detailed in Sec.II D below, but we note already at this stage that they will take the form Pl .One can now eliminate all dimensional factors in the covariant equations (3) by introducing dimensionless variables according to These dimensionless variables can be converted back into SI units by using ℏc = 1.97327 × 10 −10 eV km , so that µ ℏc = µ 1.97327 × 10 −10 eV km .
Hence, if we set ℏ = 1 = c, a rescaled radius r = µr and frequency ω = ω/µ translate into SI units according to and likewise for the mass m and other dimensional quantities.The scalar field amplitudes, in turn, are measured in units of the Planck mass, A = ÂM Pl and φ = φM Pl .
In words, for a BS scalar mass µ = 1.97327 × 10 −10 eV, the numerical value of our radial coordinate gives the BS radius in km.For other scalar masses, the radius scales accordingly, becoming larger for smaller scalar masses and vice versa.From now on, we will work with dimensionless variables only and we will drop the caret in their notation for simplicity.

C. The field equations in spherical symmetry
Our physical system is completely described by φ, A, α and X, which are all functions of radius r.It is sometimes convenient to express the metric functions α and X in terms of the alternative variables which, when combined with Eqs.(3)-( 5), results in a set of ordinary differential equations (ODEs) that can be written as Here, we have introduced the auxiliary variables κ and θ in order to write the two scalar-field equations in firstorder form and used the notation .
The system ( 13) is complemented by the boundary conditions The vanishing of the scalar fields at infinity can only be achieved for specific values of ω and φ ctr which characterizes Eq. ( 13) as an eigenvalue problem.More specifically, the scalar fields have an asymptotic behaviour containing exponential modes where h and k depend on the scalar masses and frequency ω, and the unphysical, growing modes e +hr , e +kr only vanish for discrete values of φ ctr and ω.In the literature, the asymptotic behaviour of the BS scalar is sometimes given as but this only holds in the complete absence of gravity and is incorrect even for r → ∞ if the scalar field couples to gravity, either in GR or ST theory.As we show in Appendix A, instead of Eq. ( 16), the asymptotic behaviour with gravity is where h, k, δ and ϵ are given in Eqs.(A4) and (A5) for a massive and massless φ, respectively.The ground-state BS models, that we focus on in this study, are furthermore characterized by monotonic profiles A(r) with no zero crossings except for their vanishing at infinity.

D. Choice of model
The field equations derived in the previous section hold for arbitrary potentials W (φ), V (ψ) and conformal function F (φ).We now discuss the specific choices made in our work, starting with the conformal factor F (φ).We consider the class of ST theories given by the Damour-Esposito-Farèse function [58], with constants α 0 ≥ 0 and β 0 .The choice of Damour and Esposito-Farèse is in part motivated by the fact that the modifications of gravity at first post-Newtonian order are completely determined by the asymptotic values of the first and second derivatives of ln F in the far-field limit [57,59,82].In our notation, these derivatives are In particular, this implies an effective Newtonian gravitational constant G = G(1 + α 2 0 ).The parameters α 0 , β 0 furthermore completely determine the Eddington parametrized post-Newtonian parameters β PPN and γ PPN [49,83]; see also Eq. (3.7) and the surrounding discussion in Ref. [84].For β 0 = 0, we recover Brans-Dicke theory [85], whereas Damour and Esposito-Farèse's inclusion of the quadratic term leads to the spontaneousscalarization phenomenon for sufficiently negative β 0 .
The pioneering work on spontaneous scalarization focused on massless ST gravity, which, however, has been constrained significantly by the Cassini mission and binary pulsar observations, not least because of the dramatic effect of the spontaneous scalarization mechanism.This has resulted in an increasing amount of attention devoted to massive ST theories which, due to the screening effected by the mass term µ φ , largely remain compatible with the observations.Roughly speaking, the Yukawa fall-off of a massive scalar suppresses the impact of the scalar sector on the binary pulsar dynamics for sufficiently large separations of the binary constituents; this leaves ST theories with µ φ ≳ 10 −15 eV in agreement with present observations [64].We note, however, that µ φ may also be constrained through other types of observations.Specifically, a recent study on equilibrium neutron stars conjectures that the GW event 170817 is suggestive of a larger mass µ φ ≳ 10 −11 eV [86] while the observation of highly spinning black holes, combined with their superradiant interaction with surrounding massive fields, may exclude the mass range 10 −13 eV ≲ µ φ ≲ 10 −11 eV [87,88].Our analysis is largely motivated by theoretical considerations and the goal to obtain a comprehensive understanding of BS scalarization.While bearing in mind the possible constraints above, we therefore consider a wide range of mass values µ φ ≥ 0 as well as coupling parameters α 0 , β 0 with no further prejudice.In our theoretical formalism we include the gravitational scalar mass through a quadratic potential where µ φ ≥ 0 measures the gravitational scalar mass in units of the BS scalar mass parameter µ.As we will see later both the massless and the massive cases will share common phenomenological features.Finally, to describe the bosonic sector, we focus on repulsive and solitonic potentials,  (21).The non-interacting potential for mini BSs is recovered for σ0 → ∞ or λ4 = 0. Note the false vacuum state for solitonic potentials and the systematic increase of V (A) from small σ0 to the mini-BS limit and larger λ4.
where λ 4 and σ 0 are self-interaction terms of the two individual potentials.Note that both potential functions include the case of mini BSs in the respective limits λ 4 → 0 and σ 0 → ∞.We graphically illustrate the features of these potentials in Fig. 1 for two solitonic potentials with σ 0 = 0.2 and σ 0 = 0.4, the case of mini BSs and a repulsive potential with λ 4 = 50 over the range A ∈ [0, 0.2].We clearly see the false vacuum states for σ 0 = 0.2 and also note that in this range of A, the potential values systematically increase in the order σ 0 = 0.2, σ 0 = 0.4, the mini BS limit and λ 4 = 50.

E. Numerics and Diagnostics
Accurately capturing the exponential behaviour of the scalars and determining the eigenvalues φ ctr and ω is a non-trivial numerical challenge and we detail in Appendices A and B how we compute the solutions using shooting and relaxation techniques.The conceptual summary of our calculations is as follows.
3. Parametrize the sequence of BS solutions inside this framework by varying a control parameter, typically A ctr , but sometimes also ω or Φ ctr . .= Φ(0).

4.
For each value of this control parameter, we iteratively determine the eigenvalues φ ctr and ω leading to a physical (i.e.asymptotically flat) BS solution.
For some values of the control parameter, there exist multiple solutions with different degrees of gravitational scalarization.In order to find all possible solutions, we perform the iterative procedure using different initial guesses for ω and φ ctr , usually obtained from neighbouring models with similar degrees of scalarization.
The typical result of such a calculation for a given ST theory and potential function is a set of one-parameter branches of strongly and weakly scalarized BS solutions which can be displayed, for example, as curves in the mass-radius plane.Whereas the mass of a BS is well defined as the Arnowitt-Deser-Misner [89] or ADM mass, given in our formalism by m(r → ∞), the radius is not; BSs formally extend to infinity, albeit with exponentially decaying scalar amplitude A. In this work, we define the radius R J of a BS as the Jordan frame radial coordinate r/ √ F where the mass function m reaches 99 % of the ADM mass, i.e.where m(r) = 0.99 m(∞).Alternative measures for the radius, using a different percentage or other physical diagnostics, do not change our main results beyond minor differences in the radius values.

A. A prototypical example
We begin our discussion of the results with a representative family of BS models that illustrates the main features of their spontaneous scalarization.This set of models is obtained for ST parameters µ φ = 0.05, α 0 = 0 and β 0 = −11 as well as a solitonic potential V sol (21) with σ 0 = 0.2 for the BS matter.The main diagnostic quantities we display for this family of models are the ADM mass M ADM , the radius R J , the central and maximal values φ ctr , φ max of the gravitational scalar, the central magnitude of the complex BS scalar A ctr and the frequency ω.
The resulting BS models are graphically illustrated in the form of mass-radius diagrams in Fig. 2 with the other diagnostic quantities encoded in color.We note in this context that for α 0 = 0 the BS models are invariant under the change φ(r) → −φ(r) and the profile φ(r) has no zero crossings.Unless stated otherwise, we display in this section the solutions with positive φ = |φ|.The breaking of this degeneracy for non-zero α 0 will be explored in Sec.III D below.The main features visible in this figure are as follows.
1.The BS models of GR are recovered as a oneparameter family with vanishing gravitational scalar φ; this set of models is conspicuous as the dark blue branch in the top right panel of Fig. 2.
For the vanishing α 0 of our example, this branch exactly equals the BS models computed in GR as shown, for example, in Fig. 1 of Ref. [29] for the same σ 0 .
2. Scalarized models appear most prominently in the form of a second branch forming a loop starting from about a radius R J ≈ 5 towards larger radii and masses and eventually reconnecting with the GR branch at R J ≈ 3.As can be seen in the inset of the top right panel of Fig. 2, there exists a second segment of the scalarized branch extending to the right, a bit further down.As we will see in Sec.III B, the fine structure of the highcompactness (upper) end of the scalarized branch can change with β 0 .
3. Bearing in mind our sign convention for φ, the upper two panels in Fig. 2 demonstrate that for many of the scalarized models, the gravitational scalar φ reaches its maximum away from the origin: |φ| ctr is typically smaller than |φ| max .

4.
As also observed for neutron stars [58,90], the scalarized models can reach significantly larger masses and radii3 .
5. As the central BS scalar amplitude increases, the BS models become increasingly compact, as illustrated in the bottom left panel of Fig. 2.
6.As the BS compactness increases, the frequency ω generally decreases from lim Actr→0 ω = 1.This trend is only mildly reversed at the extreme end for highly compact models; cf. the inset in the bottom right panel.The smallest frequencies are obtained in the high-mass regime of the scalarized branch.
In Fig. 3, we plot the BS scalar amplitude A, the gravitational scalar φ, the mass aspect m as well as the trace T of the physical energy momentum tensor as functions of radius r for four selected models marked by circles in Fig. 2. Here, models 1 and 3 are located on the scalarized branch with compactness C . .= max r∈R m(r) r = 0.117 and C = 0.240, respectively.The BSs marked 2 and 4 represent GR stars with respectively equal compactness for comparison.The models exhibit quite similar main features despite their different locations in the mass-radius plane: (i) Whereas the BS scalar amplitude A is a monotonically decreasing function of r for all cases, the gravitational scalar of models 1 and 3 reaches its extremal value away from the origin.(ii) The mass function increases monotonically with radius, but levels off at relatively small radius due to the exponential decrease of A. (iii) The scalarization tends to push the BS scalar profile A(r) towards larger radii and results in larger total mass M ADM .(iv) The region of strongest scalarization coincides with the region of a negative trace of the energy momentum tensor as expected for the onset of a tachyonic instability [64].
In the following sections, we will explore in more detail the phenomenology of scalarized BS models in ST theory and will highlight how and when the main features illustrated by this prototype can alter for different potentials V (A) and ST parameters.

B. Scalar mass and the onset of scalarization
Having seen the main features of scalarized BSs in the previous section, we next study how the mass µ φ of the gravitational scalar affects the scalarization phenomenon and, in particular, at which values of β 0 scalarization sets in.For this purpose, we keep the BS potential function V (A) fixed at its solitonic shape with σ 0 = 0.2; the effect of alternative potentials will be explored in Sec.III C below.
The majority of spontaneously scalarized neutron-star The maximal value of the gravitational scalar, |φ|max, is shown as a function of the central BS scalar amplitude Actr for the same BS families.Note that for α0 = 0, the models for positive and negative φ are degenerate and we plot both extrema, ±|φ|max, to illustrate the loop structure.Again the loops shrink for larger µφ, but we also notice a split into two separate loops at µφ = 0.1.For µφ = 0.15, the second loop containing highly compact models with large Actr is no longer present.models have been found below a threshold for β 0 ≲ β 0,thr = −4.35where β 0,thr turns out to be remarkably robust against variations of other parameters such as the equation of state [64,[90][91][92][93].One striking exception to this rule are the strongly scalarized solutions found by Mendes and collaborators [94][95][96] for positive β 0 provided the trace of the energy momentum tensor is positive in the GR limit; see also Ref. [64] where this effect is explained in terms of a tachyonic instability.For the BS models considered in this study, we focus on negative β 0 , but we will see that the threshold value β 0,thr can differ from that for NSs and also depends more sensitively on the BS parameters.
As a general observation we find that increasing the mass µ φ weakens the effect of scalarization.This is illustrated in Fig. 4 where we display solitonic BS families for α 0 = 0, β 0 = −12 and σ 0 = 0.2, but varying µ φ from 0 to 0.15.The upper panel in this figure demonstrates that the scalarized branches significantly reduce in terms of the BS radius and mass, getting closer to the GR branch for larger µ φ .We note in this context that the rather large BS radii found for µ φ = 0 are a consequence of a slightly slower fall-off of A(r), likely caused by its coupling to the now long-ranging φ.We investigate this effect in more detail in Sec.III E below, but note that the outer regions even of BSs of radius R J ≈ 100 are tenuous with the small energy density merely amplified by the r 2 effect of the volume element.
The decrease in scalarization becomes even more evident in the lower panel of Fig. 4 where the same BS families are represented in the plane spanned by the central BS scalar amplitude A ctr and the maximal value of |φ(r)|.Recalling that the BS models are identical under φ → −φ, we display each BS in this figure in terms of its two extrema ±|φ| max in order to illustrate more clearly the loop like structure of the branches.Clearly, the maximal scalarization decreases for larger µ φ , but we also notice two qualitative changes occuring between 0.05 < µ φ < 0.1 and 0.1 < µ φ < 0.15, respectively.First, the single loop (long-dashed blue and dotted green) splits into two seprate loops (dashed orange) leaving a finite range 0.25 ≲ A ctr ≲ 0.31 where only GR models exist.For µ φ = 0.15, the second loop of highly compact (large A ctr ) models has disappeared; we cannot entirely rule out that it is shifted to very large A ctr , but all our searches for scalarized models in this regime for µ φ ≥ 0.15 have yielded no solutions.This behaviour is also reflected in the mass-radius diagram of the upper panel, where for µ φ = 0 or 0.05 the scalarized branch connects to the GR branch only at the low mass end, but remains separate at the high-mass end as shown in the inset.For µ φ = 0.1, instead, we have two separate scalarized branch segments, both connected to GR models: note the small dashed orange branch in the inset near the end of the black GR branch.For µ φ = 0.15 (solid magenta), this small extra branch is no longer present.
The lower panel of Fig. 4 resembles Fig. 3 of Ramanazoglu and Pretorius [97] which shows |φ| max obtained for scalarized neutron-stars.Their analysis of spontaneous scalarization in terms of a tachyonic instability leads to an approximate criterion λ φ > λ eff, star ∼ R J / C|β 0 | -their Eq. ( 6) -where C is an estimate of the star's compactness and λ φ is the Compton wavelength associated with µ φ .Translated into our variables, this relation becomes and suggests an approximately quadratic relation between the threshold for scalarization β 0,thr and the scalar mass µ φ .
We have numerically computed the onset of scalarization for the families of BS models with α 0 = 0, σ 0 = 0.2 by starting for each fixed µ φ with a sufficiently negative β 0 that admits scalarized models and then iteratively increasing β 0 towards 0 until scalarization disappears.In practice, we have stopped this search at 4 significant digits in β 0 and verified that the result remains unchanged under a doubling or halving of the number of grid points.The resulting threshold β 0,thr is plotted as a function of µ φ in Fig. 5 together with a quadratic fit.Typical values for radius and compactness of the BS models on the verge of scalarization are R J ∼ 5 and C ∼ 0.2, resulting in a coefficient R 2 J /(4π 2 C) ∼ 3 in Eq. (22).Bearing in mind its approximate character, the tachyonic analysis captures the quadratic contribution of our fit rather well.The constant and linear contributions to our fit, in turn, are likely a direct consequence of a finite threshold for scalarization in massless ST theory: using a series of self-gravity contributions, Damour and Esposito-Farèse conclude that the non-perturbative amplification effect of spontaneous scalarization can set in when β 0 ≲ −4.For the massless case considered in their analysis, this remains a remarkably good estimate for BSs.

C. Dependence on the bosonic potential
In this section, we consider the structure of BS solutions in ST theory of gravity for other potentials of the BS scalar field.We focus on models governed by the potentials given in Eq. (21), namely a solitonic one with σ 0 ≥ 0.2 and a repulsive potential with λ 4 ≥ 0. We find that the specific choices we consider here affect the structure of scalarized BS models, resulting in several differences to the prototypical example studied in Sec.III A. In particular, we find that the gravitational scalar changes its shape with stellar compactness and that strongly scalarized solutions may no longer have larger masses and radii than their GR counterparts.We analyse these features in more detail in the following subsections.

Massless gravitational scalar field
Gravitational scalar profiles: Previously we noted that a solitonic potential with σ 0 = 0.2 leads to the distinct we have ω1 = 0.8303 and ω2 = 0.7572.Note that the trace of the energy-momentum tensor has been scaled for presentation purposes.As supported by analysis in the main text, T becomes negative at small radii for models with a gravitational scalar peaking at the origin.Lower panels: Mass-radius plots with the colour bar indicating the central amplitude of the BS scalar field.The star marks the threshold at which the gravitational scalar starts to peak off centre; we ubiquitously observe that φ peaks at r = 0 for less compact stars (Actr is lower than for the threshold case, i.e. dark blue color) and φ peaks off centre for more compact stars (large Actr, i.e. cyan color).We exemplify this with models 1 (less compact with φ peaking at r = 0) and 2 (more compact with φ peaking off-center) for each of the potentials (left with λ4 = 0 and right with λ4 = 50).
feature of a gravitational scalar field φ(r) peaking offcentre; cf.Fig. 3.For repulsive (λ 4 > 0) and mini (λ 4 = 0) BS models, in contrast, the gravitational scalar more commonly reaches its global extremum exactly at the centre.We illustrate this in Fig. 6, where φ(r) gradually changes its behaviour from peaking at r = 0 to peaking off-centre, as one moves to models with higher central amplitudes A ctr of the BS scalar.The same behaviour has been observed for neutron stars of increasing compactness; cf.Fig. 7 in Ref. [90].In Fig. 6 we mark by a star the threshold, where this transition occurs.We additionally plot in each of the upper panels the profiles A(r) and φ(r) for two models marked 1 and 2 in the mass-radius diagram.These models are located on either side of the threshold: for model 1 (less compact, i.e. smaller A ctr than the threshold model) φ peaks at the centre r = 0 and for model 2 (more compact with larger A ctr ) it peaks off centre.
We can qualitatively understand this behaviour as follows: the shape of the gravitational scalar field near the origin is controlled by Eq. ( 13) for κ, which denotes a re-scaled version of the radial derivative for the gravitational scalar field ∂ r φ.The overall sign of ∂ r κ at small radii determines the concavity of φ near the origin, and, thus, whether it peaks off or at the centre.Without loss of generality, we assume that φ ctr > 0, then the sign of ∂ r κ will be determined by the two terms: 2αXF W ,φ 2 φ and ω 2 A 2 − 2α 2 V − F 2 θ 2 .The former is always positive or zero, whilst the latter can be negative, provided the magnitudes of the potential and the lapse (i.e.2α 2 V ) are large enough to counteract the positive contribution of ω 2 A 2 .The last term F 2 θ 2 is negligible near the origin since θ(0) = 0 by Eq. ( 14).The ingredients required to make ∂ r κ negative, and thus cause the gravitational scalar φ to peak at the origin, are as follows: 1. Large potential values V : for typical values of A ctr ∈ [0.1, 0.2] in the regime of scalarization, the non-interacting and the repulsive potentials take on larger values compared to a solitonic potential with small σ 0 ; cf.Fig. 1.
2. Small central amplitudes of the BS scalar: models with small A ctr result in less compact (or fluffy) stars.These fluffy configurations have shallower gravitational wells which manifest themselves in larger values of the lapse α near the origin 4 .Small A ctr and the ensuing large α therefore are more likely to make the overall contribution to ∂ r κ negative.Empirically, we find that mini and repulsive potentials, but rarely solitonic with σ 0 ≈ 0.2, allow for scalarized BSs with small A ctr , i.e. low compactness; cf.Fig. 7.
This behaviour of the gravitational scalar can also be explained using the trace of the energy-momentum tensor which is obtained from Eq. ( 4), Up to an overall factor, this expression equals the second source term of our above analysis that drives ∂ r κ via Eq.( 13).The sign of T can therefore be directly translated to the concavity of φ.Recalling that the tachyonic instability requires negative T , we expect strong scalarization to occur in the region where T < 0. Fig. 6 illustrates exactly this effect: for less compact models with T ≲ 0 at small r, φ peaks at the origin, and for compact models with a central T ≳ 0, φ peaks away from the origin; cf. also Fig. 3.The φ(r) profile of particularly compact scalarized models is therefore more shell -like.
The onset and degree of scalarization: The nature of the potential further determines the degree of scalarization of BS solutions and the global properties of such models.Figure 6 illustrates one example of this for mini and repulsive potentials in the mass-radius plane, where all of the scalarized branch falls underneath the GR one for λ 4 = 0, but only part of it for λ 4 = 50.Similarly, we see in Fig. 16 in Appendix C that many scalarized BSs for solitonic potentials with σ ≳ 0.3 have smaller masses than their GR counter parts with equal radius, very much in contrast to the solitonic σ 0 = 0.2 case of the prototypical example studied in Sec.III A.
Furthermore, the threshold for β 0 at which scalarization starts to occur varies with the potential; we summarize this dependence in Table I.Overall, we find that the repulsive potential requires less negative β 0 values for solutions to scalarize, which is then followed by the mini and solitonic potentials.For massless ST theory, the thresholds fall into a range −7.5 ≲ β 0,thr ≲ −4.96, a bit below the threshold of −4.35 typically found for neutron stars.
Even though solutions are less strongly scalarized for the solitonic potential with σ 0 ≳ 0.3, there are some subtle features in the dependence of scalarization on the selfinteraction term as illustrated in Fig. 7. First, the degree of scalarization varies non-monotonically with σ 0 .Typically smaller values of σ 0 ≈ 0.2 result in more strongly scalarized solutions than models with σ 0 ≈ 1.On the other hand, larger σ 0 lead to a wider range of A ctr for which scalarized solutions exist: in Fig. 7 the support of |φ| max (A ctr ) significantly increases for models with σ 0 ≳ 0.35.As σ 0 increases, the scalarized branches slowly converge to the scalarized solutions of mini BSs, which can be seen by a close overlap of the σ 0 = 2 and λ 4 = 0 curves.We present a more detailed illustration of this transition in Appendix C, where we compute families of solutions for a wider range of potentials.Finally, we find that for a fixed value of β 0 , the strongest scalarization occurs in the case of the repulsive potential with our largest λ 4 = 100.

Gravitational scalar profiles:
We next turn to analysing the features of scalarized BS solutions with other bosonic potentials in the massive case.A systematic investigation of the global maximum of |φ| across

I II
FIG. 8. Flow-chart showing the dependence of the gravitational scalar profile φ on the BS parameters, Actr and V (A), and the ST parameters, µφ and β0.As discussed in the main text, the sign of ∂rκ determines the behavior of φ at the origin.If ∂rκ > 0, φ peaks away from r = 0, whilst for ∂rκ < 0, φ peaks at r = 0.The two source terms, denoted by I and II in the diagram, play the leading role in determining the sign of ∂rκ at r = 0.By the boxes' colours we mark how the ST and BS parameters tend to impact the φ profiles -blue for peaking at the origin and pink for peaking off centre.
the parameter space confirms that for µ φ > 0, certain choices of the potential V (A) still result in the gravitational scalar profiles peaking at the origin.However, we find that an increasing mass µ φ tends to push the extremum of φ away from the origin.This effect can be explained using the same reasoning as in the massless case of Sec.III C 1 by assessing the sign of the righthand-side of Eq. ( 13) for ∂ r κ.Assuming that φ ctr > 0, we conclude that µ φ > 0 increases the overall positive contribution from the 2αXF W ,φ 2 φ term and empirically we find it also increases ω.Therefore, larger values of the gravitational mass push ∂ r κ towards positive values and the gravitational scalar φ tends to peak away from the origin.It is possible however to push the gravitational scalar field maximum towards r = 0 by making β 0 negative enough.The impact of a more negative β 0 parameter is two-fold: The combined effect of a more negative β 0 is therefore to make the final term of the ∂ r κ equation negative and large enough to counteract the non-zero positive contribution from 2αXF W ,φ 2 φ term.The sign of ∂ r κ becomes negative and φ peaks at r = 0.
We quantitatively illustrate the effect of µ φ and β 0 on the φ profile in Table II where we list for several combinations of (µ φ , β 0 ) the minimal σ 0 values at which we find models with φ peaking at r = 0. Below the respective σ 0,min , φ peaks off center for all scalarized BS models.As the table suggests, it is easier to find φ profiles peaking at the origin for smaller gravitational masses µ φ and more negative β 0 .In Fig. 8 we summarize the dependence of the gravitational scalar profile on the BS and ST parameters.21) with the indicated parameter values σ0 and λ4.The lines represent quadratic fits with the given coefficients.The case λ4 = 100, µφ = 1.5 exhibits one minor anomaly: here, the × symbol marks the disappearance of stable scalarized models, but extremely compact unstable BSs are still found up to the β0 value marked by the filled circle.Only the former (×) data point is used in the fit.
The onset and degree of scalarization: Finally, we make a few remarks on the effect of µ φ > 0 on the scalarization of BS solutions.As already seen in Fig. 5, the threshold β 0,thr for the onset of strong scalarization for solitonic BSs with σ 0 = 0.2 is well described by a quadratic function of the mass parameter µ φ for ST gravity with α 0 = 0.In Fig. 9 we plot the onset of scalarization for several other potential functions with parameters as listed in the legend.For all cases, the numerical data are in excellent agreement with a monotonically decreasing quadratic fit which confirms our general observation that an increasing µ φ weakens scalarization.In order to assess in more detail how different potential functions affect the onset of scalarization, we display in Fig. 10 the threshold values β 0,thr as a function of the parameters σ 0 and λ 4 for three fixed values of the mass µ φ = 0, 0.1 and 0.3.We recall that mini BSs are obtained in both limits, σ 0 → ∞ and λ 4 = 0, and can be regarded as the link between the solitonic and repulsive potential.The variation of β 0,thr in Figs. 9 and 10 shows that for µ φ ≲ 0.3, similar to the massless case, BSs with a repulsive potential are more susceptible to scalarization than solitonic and mini BSs, although not by a huge amount.For µ φ ≳ 0.3, however this trend is reversed and solitonic BSs scalarize most easily.

D. Dependence on α0
Previously we have considered ST theories with α 0 = 0, where the solutions are degenerate under the transformation φ → −φ.In the case, of α 0 ̸ = 0 we acquire additional structure through the splitting of the scalar- ized branches.More specifically, we find5 : 1.A larger branch (as viewed in the mass-radius plane), mostly with negative gravitational amplitudes φ ctr .At very high compactness (i.e.large A ctr ) these models start to have φ ctr > 0 and the φ(r) profiles have a zero-crossing.
2. A smaller closed loop, where we exclusively observe positive φ ctr .
These two branches contain strongly scalarized solutions, as well as weakly scalarized ones, in which case the BS models have similar mass and radius to their GR counterparts.This is illustrated in Fig. 11, where certain BS models of the α 0 > 0 sequences closely resemble the GR solutions obtained for α 0 = 0.For mild values of α 0 , the closed loop intersects itself and acquires a shape similar to the ∞ symbol as in the case of α 0 = 0.01.As we increase α 0 , first the ∞-shape and then the entire loop slowly shrink and disappear (see e.g.α 0 = 0.1).We note that certain values of α 0 we consider here are above the Cassini measurement constraint.This choice is made for visualization purposes: for smaller values of α 0 the two branches lie more closely to each other making their splitting less evident.Although here we discuss solitonic BSs, for other potentials we observe a qualitatively similar splitting of the solution branches.
The disappearance of one of the solution branches can be traced back to the conformal coupling function F (φ) = exp(−2α 0 φ − β 0 φ 2 ).When α 0 = 0, the value F is solely determined by the magnitude of β 0 : a more negative β 0 leads to stronger scalarization.A non-zero α 0 , on the other hand, leads to a non-zero −2α 0 φ contribution in the exponent of F .For φ > 0, this curbs the effect of β 0 and subsequent scalarization whilst for φ < 0 the α 0 term strengthens it: for α 0 > 0 we have F (φ > 0) < F (φ < 0).This effect is further illustrated by the relative size of the scalarized branches in Fig. 11: the solution branches for φ ctr > 0 are systematically smaller in size than their φ ctr < 0 counterparts.For sufficiently large values of α 0 with fixed β 0 , the solutions with φ ctr > 0 cease to exist and the branch disappears.Considering the case α 0 = 0.1 and β 0 = −7, for example, we can easily support this by a back-of-the-envelope calculation: assuming φ ∼ 0.1, the conformal coupling function is F (φ) ∼ exp(0.05),which is equivalent to setting α 0 = 0 and β 0 = −5.As discussed in Section III B (see also Fig. 5 and Table I), however, the onset of scalarization for µ φ = 0 and a solitonic potential with σ 0 = 0.2 occurs for β 0 ≳ −6.1.An effective β 0,eff ≈ −5 falls out of this range and is therefore compatible with the disappearance of positive solutions for α 0 = 0.1.
Besides the branch splitting, BS solutions with α 0 > 0 depend on β 0 and µ φ in a way very similar to the α 0 = 0 case discussed previously.In particular, a more negative β 0 results in stronger scalarization of solutions, with larger mass M ADM and radius R J (cf. the right panel of Fig. 11).Increasing the mass parameter µ φ , on the other hand, systematically reduces the scalarization effect in the same way we have observed in Sec.III B for α 0 = 0.

E. Thin-shell models
In Ref. [98], Collodel and Doneva computed an intriguing type of BS solutions in GR dubbed thin-shell models.The main characteristics of these BSs are (i) an approximately constant scalar profile A(r) extending out to rather large radius and (ii) a low frequency ω; see in particular their Fig. 6.As a consequence of these two features, the derivatives ∂ t A and ∂ r A become very small at small as well as very large radii, resulting in a nearly vanishing energy density except for a shell region where the scalar field eventually drops to zero.In Fig. 2 for our prototypical BS family, we can see that the scalarized stars tend to have lower frequencies than their GR counterparts.We have found this behaviour systematically in our investigations of the BS and ST parameter space (see Fig. 19), suggesting that scalarization can strengthen the thin-shell character of BSs.
We study this feature in more detail for the specific case of massless ST theory with α 0 = 0 and β 0 = −12 for the solitonic potential with σ 0 = 0.2.Due to different conventions, our σ 0 is related to the parameter σ used by Collodel and Doneva by σ = √ 2πσ 0 , so that our scenario corresponds to σ ≈ 0.5 in Ref. [98], well above the regime where they encounter thin-shell stars.We likewise find no indication of shell-like structure on the GR branch of our scenario which yields a minimal frequency of about ω = 0.428 in good agreement with Ref. [98]; this can be seen by comparing the bottom panel of our Fig. 12 with the left panel of their Fig. 2. Along the scalarized branch, however, we obtain BS models with much lower frequencies ω < 0.2.As indicated in the upper panel of Fig. 12, we also obtain significantly larger radii for these models, although we emphasize that this is largely due to the reduced fall-off of A(r) rather than its behaviour in inner regions of the star.The radial profiles displayed in Fig. 13 for the strongly scalarized BS with minimal ω (marked by circles in Fig. 12), however, exhibit an approximately constant A(r) out to r ≈ 7.In contrast, the GR model with minimal ω displays a distinctly "ordinary" profile where A(r) rapidly drops away from the origin.A similar investigation of neighbouring models along the GR and scalarized branches yields very similar profiles and, as evident in Fig. 12, frequency values, and we observe qualitatively the same behaviour for other values of σ 0 .Mathematically, we can explain this behaviour by considering Eq. ( 13) for ∂ r A and ∂ r θ.Thin-shell models can be regarded as a sequence of BSs that approach a uniformly constant scalar-field solution A = const extending to larger and larger radii.In Eq. ( 13), this limiting (Minkowski) solution is obtained for a zero right-hand side of ∂ r θ.Ignoring the geometric term −2θ/r, which merely incorporates the r 2 behaviour of a 3D volume, the key mathematical origin of thin-shell stars is a nearly vanishing For strongly scalarized BSs, on the other hand, F (φ) is dominated by the quadratic contribution in the exponential (2) and therefore F > 1 which reduces the term (24) relative to its GR value and hence drives BSs in the direction of thin-shell models.

IV. STABILITY
Depending on the ST parameters (µ φ , α 0 , β 0 ) and the potential V (A), there may exist up to seven (or even more for non-zero α 0 ) BS models with the same mass.This raises the question which of these models is energetically favored.We assess this by using binding energy arguments for a wide range of BS solutions.More specifically, we define the binding energy where Q BS denotes the conserved Noether charge, i.e. the number of bosonic particles making up the BS.It is defined as the spatial integral of the time component of the Noether current, so that Left: BS models in ST gravity with µφ = 0, β0 = −7 and varying α0 for a solitonic potential with σ0 = 0.2.The branches obtained for each value of α0 are presented in different shades of one color to differentiate between positive and negative central gravitational field amplitudes.For α0 = 0.01, for instance, dark green denotes models where φctr < 0 and light green those with φctr > 0. For reference, we also plot the solutions for α0 = 0 which are shown in grey.Right: Dependence of the BS solutions on β0 for α0 = 0.01, µφ = 0 and a solitonic potential with σ0 = 0.2.The inset highlights the self-intersecting φctr > 0 branch with ∞ shape for the sequence of models with β0 = −10.
For each set of BSs with equal Noether charge 6 , the model with the strongest binding energy E b , i.e. the smallest mass M ADM , is taken as the stable, energetically favored configuration.This does not necessarily imply that the other models are perturbatively unstable, but under sufficiently strong excitations, we expect them to either form a black hole, evaporate or migrate to a stable BS; we therefore refer to these other models as unstable in the following discussion.
We first apply this method to the prototypical example discussed in Sec.III A. This set of BS models is obtained for ST parameters µ φ = 0.05, α 0 = 0, β 0 = −11 and a solitonic potential with σ 0 = 0.2. Figure 14 shows the same mass-radius diagram as Fig. 2 but now color codes solutions in light (copper) for stable or dark (black) for unstable stars.We clearly see that the scalarized branch from the low-mass bifurcation point all the way up to its maximal mass at about (R J ≈ 15, M ADM ≈ 1.1) contains energetically favored models.This set of stable BSs is completed by a small piece of the GR branch around R J = 5 which contains low-mass but comparatively compact BSs for whom no scalarized counterparts exist.
Repeating the same analysis in Fig. 15 for other bosonic potentials and ST parameters, unravels a very similar pattern.For α 0 = 0 scalarized BSs are commonly energetically favored relative to their GR cousins with equal Noether charge.For non-zero α 0 , where we commonly encounter multiple scalarized branches, we observe that strongly scalarized BSs with larger radius tend to be stable.We identify some exceptions to these rules, however: scalarized BSs with relatively small mass or radius may be energetically disfavored relative to other scalarized solutions (see the left panel of Fig. 15) and, in some cases, even relative to GR models (see Fig. 20 in Appendix C, where we explore a wider range of potentials V ).In general, we observe that, similarly to BS models in GR, stable and unstable models in ST gravity are separated by the maximum-mass model.

V. CONCLUSIONS
In this work we have performed a systematic study of the structure of strongly and weakly scalarized boson-star solutions in scalar-tensor theory of gravity with Damour-Esposito-Farèse coupling.For this purpose we compute spherically symmetric ground-state BS models using shooting and relaxation algorithms.In order to control the challenging exponential behaviour of the scalar field solutions, we perform a comprehensive analysis of their asymptotic behaviour, identifying incompatibilities that arise from using the flat-field limit.
Our study suggests that there are many common features in the dependence of scalarized solutions on the ST parameters (µ φ , α 0 , β 0 ) as compared to the neutron star case.In particular, we observe the following main features: 1.The scalarized solutions form additional branches in the mass-radius plane.These models can have larger mass and radius compared to their GR counterparts, although exceptions to this rule exist, depending on the potential and ST parameters.4. Scalarized solutions are often energetically preferred (i.e.stable) when compared to their GR counterparts with the same Noether charge.We similarly find that a larger BS radius tends to favour stability.
Further to these main features, we summarize the specific traits inherent to scalarized BSs and their variation over the parameter space as follows.
Onset of scalarization: Compared with neutron stars, we require more negative values of β 0 to obtain scalarized BSs.Our analysis suggests that BSs with repulsive potentials with large λ 4 are most susceptible to scalarization in the massless case: the onset of scalarization happens for β 0 ≲ −4.96, compared to β 0 ≲ −4.35 for NSs.Mini and solitonic BSs are less liable to scalarize and typically require more negative β 0 .This trend also holds in the massive case up to µ φ ≈ 0.3, beyond which solitonic BSs become more receptive to scalarization.For all potentials, the threshold for scalarization β 0,thr takes on a quadratic dependence on the mass term µ φ ; cf.Fig. 9.
Degree of scalarization: The repulsive potential for large λ 4 results in the strongest scalarization, both in magnitude and range of BS models.For solitonic potentials, in contrast, we observe a non-monotonic dependence of the degree of scalarization on the self-interaction term σ 0 .In general small σ 0 ≈ 0.2 leads to stronger scalarization albeit over a narrower range of BS models than larger values σ 0 ≳ 1 and the mini BS limit σ 0 → ∞; cf.Fig. 7.
Gravitational scalar behaviour : The potential not only determines the onset and degree of scalarization, but also the resulting compactness of the star and the gravitational scalar behaviour near the origin; cf.Fig. 8. Generally, we find that less compact stars support a gravitational scalar whose radial profile peaks at the origin; this is most common for mini and repulsive potentials that result in less compact stars.As attractive self-interaction terms are included in the solitonic potential, the scalarized solutions become more compact and the gravitational scalar's extremum moves away from the origin.We also relate this behavior to the trace of the energymomentum tensor T , which needs to be negative in the regions of strong scalarization.We systematically observe that less compact models achieve negative values of T near the origin; cf.Fig. 6.For more compact models, in contrast, T only becomes negative further away from r = 0 resulting in shell-like gravitational scalar profiles.
Thin-shell like models: Strong self-interactions in the solitonic potential (small σ 0 ) paired with strong scalarization (small gravitational scalar mass µ φ and/or very negative β 0 ) result in BS models akin to the thin-shell BS models found in GR [98].These models correspond to highly massive BS solutions in their family and have a bosonic field amplitude A(r) approximating the shape of a Heaviside function.In GR, typically small values of σ 0 ≲ 0.15 are required to obtain such thin-shell like solutions.In ST theory, however, we easily find such models already for σ 0 = 0.2.
Our results demonstrate the rich structure of BSs in ST theory of gravity and open many opportunities for future work.A natural extension consists in further exploration of more extreme regions of the parameter space, as for example a thorough study of the onset of scalarization for thin-shell models with σ 0 < 0.2.Furthermore, dynamical evolutions of our models would provide a robust picture of their stability complementary to our binding energy estimates.Finally, it will be intriguing to see how the peculiarities identified in the structure of single BSs for different ST parameters and/or BS potentials affect the gravitational-wave emission of binary systems including the scalar or breathing polarization mode.
All frequencies we report in this work correspond to a unit lapse function at infinity, i.e.Φ 0 = 0; even if a non-zero Φ 0 is used in the numerical calculation, this is straightforwardly achieved by an a-posteriori rescaling of ω.The frequencies thus obtained ubiquitously fall into the range 0 < ω < 1.The functional behaviour k(ω) = (1 − ω 2 ) furthermore pushes k towards values rather close to unity in exactly the same way the Lorentz factor γ = 1/ √ 1 − v 2 reduces relativistic effects for all velocities but those close to the speed of light.This effect avoids extreme values of ϵ which we typically find to be in the range −2 ≲ ϵ ≲ 3.An unscalarized low-density mini BS with M = 0.345 and ω = 0.9785, for instance, has ϵ = −1.53whereas the strongly scalarized thin-shell BS with M = 2.39 and ω = 0.148 of the (µ φ = 0, α 0 = 0, β 0 = −12, σ 0 = 0.2) family discussed in Sec.III E has ϵ = 2.31.The exponential factor h for the asymptotic behaviour of the gravitational scalar, on the other hand, is directly given by the scalar mass µ φ and the exponent δ merely acquires an extra factor of the BS mass M ; δ vanishes in massless ST theory, but can reach values up to order unity for the BS sequences studied in this work.We similarly obtain values ϱ 1 = O(1) for strongly scalarized stars.
Finally, let us illustrate where and how the flat-field asymptotics go wrong, when used in the presence of gravity.Using the flat-field asymptotic behaviour leads to rescaled scalar variables

50 FIG. 1 .
FIG.1.Four example cases of the BS scalar potential function V (A) defined in Eq.(21).The non-interacting potential for mini BSs is recovered for σ0 → ∞ or λ4 = 0. Note the false vacuum state for solitonic potentials and the systematic increase of V (A) from small σ0 to the mini-BS limit and larger λ4.

2 FIG. 2 .
FIG. 2. BS models in ST gravity with µφ = 0.05, α0 = 0 and β0 = −11 for a solitonic potential with σ0 = 0.2.Each point in these mass-radius diagrams represents a BS model.The location of the point represents the star's mass MADM and radius RJ, and the color coding in the respective panels displays the central gravitational scalar |φ|ctr, the maximum (over radius) gravitational scalar magnitude |φ|max, the central BS scalar amplitude Actr and the frequency ω.The four circles labelled 1 to 4 mark the models displayed in more detail in Fig. 3.
FIG.3.Radial profiles of the BS scalar amplitude A, the gravitational scalar φ, the mass function m and the trace T of the energy momentum tensor for four selected models marked 1 to 4 in Fig.2.The pair of models 1, 2 represents one scalarized and one GR star with equal compactness C = 0.117 and amplitudes Actr as indicated in the legend.Likewise, 3, 4 represent a scalarized and GR model with C = 0.240 each.As indicated by the dashed curves for GR and the solid lines for scalarized models, the scalarization tends to push the BS scalar A(r) towards larger radii and results in larger masses.In the region of strongest scalarization in φ, the energy momentum tensor has negative trace T .

2 FIG. 4 .
FIG.4.Upper panel: Mass radius diagram for four values of the scalar mass.The scalarized branches decrease in size for larger µφ.Lower panel: The maximal value of the gravitational scalar, |φ|max, is shown as a function of the central BS scalar amplitude Actr for the same BS families.Note that for α0 = 0, the models for positive and negative φ are degenerate and we plot both extrema, ±|φ|max, to illustrate the loop structure.Again the loops shrink for larger µφ, but we also notice a split into two separate loops at µφ = 0.1.For µφ = 0.15, the second loop containing highly compact models with large Actr is no longer present.

50 FIG. 6 .
FIG.6.Upper panels: Radial profiles of the BS scalar amplitude A, the gravitational scalar φ and the trace of the energymomentum tensor T for the BS Models 1 and 2 indicated by a circle and a triangle, respectively, in the panels below.The models with λ4 = 0 (left) have frequencies ω1 = 0.8518 and ω2 = 0.7729, whilst for models 1 and 2 with λ4 = 50 (right) we have ω1 = 0.8303 and ω2 = 0.7572.Note that the trace of the energy-momentum tensor has been scaled for presentation purposes.As supported by analysis in the main text, T becomes negative at small radii for models with a gravitational scalar peaking at the origin.Lower panels: Mass-radius plots with the colour bar indicating the central amplitude of the BS scalar field.The star marks the threshold at which the gravitational scalar starts to peak off centre; we ubiquitously observe that φ peaks at r = 0 for less compact stars (Actr is lower than for the threshold case, i.e. dark blue color) and φ peaks off centre for more compact stars (large Actr, i.e. cyan color).We exemplify this with models 1 (less compact with φ peaking at r = 0) and 2 (more compact with φ peaking off-center) for each of the potentials (left with λ4 = 0 and right with λ4 = 50).

10 FIG. 7 .
FIG.7.Maximum of the gravitational scalar field as a function of the central bosonic amplitude for various potentials considered in this work.We remark that repulsive potentials result in more strongly scalarized solution than any other potential considered here.

2 ϕ 2 ϕ α 0 = 0 FIG. 9 .
FIG.9.The threshold β 0,thr for the onset of scalarization as a function of µφ for different BS potentials V (A) as given in Eq. (21) with the indicated parameter values σ0 and λ4.The lines represent quadratic fits with the given coefficients.The case λ4 = 100, µφ = 1.5 exhibits one minor anomaly: here, the × symbol marks the disappearance of stable scalarized models, but extremely compact unstable BSs are still found up to the β0 value marked by the filled circle.Only the former (×) data point is used in the fit.

FIG. 10 .
FIG.10.The threshold β 0,thr for the onset of strong scalarization is shown for the different types of potentials and three values of the gravitational scalar mass µφ = 0, 0.1 and 0.3.The horizontal axis represents σ0 for the solitonic potential on the left and λ4 for the repulsive potential on the right.

− 7 ,
FIG.11.Left: BS models in ST gravity with µφ = 0, β0 = −7 and varying α0 for a solitonic potential with σ0 = 0.2.The branches obtained for each value of α0 are presented in different shades of one color to differentiate between positive and negative central gravitational field amplitudes.For α0 = 0.01, for instance, dark green denotes models where φctr < 0 and light green those with φctr > 0. For reference, we also plot the solutions for α0 = 0 which are shown in grey.Right: Dependence of the BS solutions on β0 for α0 = 0.01, µφ = 0 and a solitonic potential with σ0 = 0.2.The inset highlights the self-intersecting φctr > 0 branch with ∞ shape for the sequence of models with β0 = −10.

FIG. 12 .
FIG. 12. Top: Mass-radius diagram for solitonic BS models with µφ = 0, α0 = 0, β0 = −12 and σ0 = 0.2.Bottom: Frequency ω as a function of the central BS scalar amplitude Actr for the same family of models.The color represents the maximal gravitational scalar in both panels.

2 FIG. 14 .
FIG. 14. Mass-radius diagram for BS solutions for µφ = 0.05, α0 = 0, β0 = −11 and a solitonic potential with σ0 = 0.2.Stable models are determined by identifying among the set of all stars with equal Noether charge QBS the BS with the strongest binding energy E b = MADM − QBS.These stable models are displayed in light (copper) color and the remaining unstable stars in dark (black).

10 FIG. 16 .
FIG.16.Boson-star families in the mass-radius plane for fixed ST parameters α0 = 0, β0 = −10, µφ = 0 and BS potential functions V (A) as given in Eq.(21).The parameters σ0 or λ4 are specified in each panel for the corresponding BS family.The mini-BS case is recovered in either of the limits σ0 → ∞ or λ4 = 0.The color bar marks the degree of scalarization in terms of the maximum of |φ(r)|.

10 FIG. 20 .
FIG.20.Same as Fig.16but with the stability coded in color, light (copper) for stable BSs and black for unstable models.

TABLE I .
Summary of threshold values of β0 for the onset of scalarization for different potentials in the massless case. 2