Structure of neutron stars in massive scalar-tensor gravity

We compute families of spherically symmetric neutron-star models in two-derivative scalar-tensor theories of gravity with a massive scalar field. The numerical approach we present allows us to compute the resulting spacetimes out to infinite radius using a relaxation algorithm on a compactified grid. We discuss the structure of the weakly and strongly scalarized branches of neutron-star models thus obtained and their dependence on the linear and quadratic coupling parameters $\alpha_0$, $\beta_0$ between the scalar and tensor sectors of the theory, as well as the scalar mass $\mu$. For highly negative values of $\beta_0$, we encounter configurations resembling a"gravitational atom", consisting of a highly compact baryon star surrounded by a scalar cloud. A stability analysis based on binding-energ calculations suggests that these configurations are unstable and we expect them to migrate to models with radially decreasing baryon density {\it and} scalar field strength.


Introduction
Ever since its formulation in 1915, Einstein's general relativity (GR) has been a tremendously successful theory of gravity, combining mathematical elegance with enormous predictive power. Phenomena ranging from Mercury's perihelion precession to the formation of black holes (BHs), the generation of gravitational waves and the big bang, find a mathematical description within this single theory. A wide range of lab-based experiments, solar-system tests and observations of astronomical phenomena have systematically scrutinized the accuracy of the theory's predictions and unanimously seen GR passing these tests with flying colors [1]. With the advent of gravitational-wave (GW) astronomy, marked by the detection of GW150914 by LIGO [2], new tests of GR have become possible in spacetime regions with strong and dynamical gravitational fields and sources moving at relativistic velocities. Once again, all GW observations so far are compatible with GR [3][4][5][6][7]. Notwithstanding GR's success, the search for possible alternative theories of gravity has for many decades been a highly active area of research [1,[8][9][10][11][12], motivated by important theoretical considerations, such as the incompatibility of GR with quantum theory at a fundamental level, as well as open questions domain, all the way out to infinity, while maintaining complete control over exponentially diverging solutions. Second, to present an in-depth analysis of the structure of the different solution branches and their dependence on the parameters of the ST theory. We begin this discussion in Sec. 2 with a review of the field equations governing the stars. In Sec. 3, we describe the numerical framework used for our computations. Our results on the structure of NS solutions in massive ST gravity are presented in Sec. 4 and we conclude in Sec. 5.

Formalism
The Bergmann-Wagoner class of ST theories, i.e. theories involving a single scalar field that are governed by two-derivative, covariant field equations and obey the weak equivalence principle, can be described in terms of the action [69] Here, ψ m collectively denotes the matter fields and S m represents their coupling to the spacetime geometry of the physical or Jordan metric g µν with determinant g and Ricci scalar R. The functions F(φ) and Z(φ) encapsulate the nonminimal coupling of the scalar field φ to the metric sector, and V(φ) is the potential function. As we shall see shortly, the function Z can be eliminated through an appropriate redefinition and is therefore often set to unity in the literature; see e.g. [70]. This class of theories is conveniently described in the so-called Einstein frame, obtained from the physical or Jordan frame through a conformal transformation of the metric and a redefinition of the scalar degree of freedom and its potential, where F ,φ = dF/dφ. In terms of these new functions, the action (1) becomes where an overbar distinguishes tensors in the Einstein frame from their Jordan counterparts. Henceforth, we use natural units where c = G = 1, unless stated otherwise. By transforming to the Einstein frame, we have eliminated the function Z. The equivalence (or lack thereof) of the Einstein and Jordan frame formulations has been the subject of a long standing debate (see e.g. [71][72][73] and references therein). Without entering this debate here, we merely note the extra freedom that the function Z introduces to the transformation (2) between the frames and henceforth follow the recommendation of Ref. [71] and work in the Einstein frame.
To complete the description of the gravitational theory we must specify the remaining free functions F and V. Following most previous studies in the literature, we write the conformal factor as 1 and take as our potential function the quadratic function which describes a non-self-interacting scalar field of mass µ. The field equations obtained by varying the Einstein action (3) with respect toḡ αβ , ϕ are given bȳ with the energy momentum tensor for which the Bianchi identity now implies the following conservation law From now on, we restrict our attention to spherically symmetric, time independent stellar models. More specifically, we employ polar slicing and radial gauge in the Einstein frame, so that the line element is of the form ds 2 =ḡ µν dx µ dx ν = −F α 2 dt 2 + F X 2 dr 2 + r 2 dΩ 2 , where α and X as well as the scalar field ϕ are functions of the radius r, and dΩ 2 denotes the standard line element on the unit 2-sphere. It is furthermore common practice to introduce the gravitational potential Φ and mass function m according to In this work we explore the behaviour of NSs in equilibrium and model their matter as a perfect fluid at zero temperature; the temperature of NS interiors in equilibrium, despite being of order 10 6 K, is well below the Fermi temperature O(10 11 ) K of matter at nuclear densities. The energy momentum tensor is then given in terms of the baryon density ρ(r), the specific enthalphy h(r) and the pressure P(r) bȳ and where is the specific internal energy. Inserting the Einstein frame metric (10) and the energy momentum tensor (12) into the field equations (6)-(9), we obtain the set of differential equations ∂ r P = −ρhFX 2 m r 2 + 4πr In order to close this set of equations, we need to relate the pressure and internal energy to the baryon density. In this work, we use cold polytropic EOSs with exponent Γ, The internal energy is then determined by the first law of thermodynamics dE =dQ − pdV for adiabatic processes withdQ = 0. For a total baryon number N and mass per baryon m b , the specific internal energy and baryon density are given by = E/(Nm b ) and ρ = m b N/V, respectively, and the first law results in The set of equations (13)- (17) can then be solved subject to the boundary conditions at r = 0 : The computation of solutions to this problem is complicated by three issues, which we list in increasing order of difficulty.
(i) The boundary conditions are specified at different locations of the domain, so that we have a two-point-boundary-value problem. (ii) For realistic values of the polytropic exponent Γ, the pressure will reach zero at a finite radius R S ; at this point, we need to match to an exterior solution with vanishing baryon density ρ. (iii) The asymptotic behaviour of the scalar field near infinity is determined by the scalar mass µ and is given by for constants A 1 , A 2 . We are only interested in bounded solutions with ϕ ∝ e −(µ/h)r /r. This exponential fall-off is responsible for the suppressed scalar contribution in the interaction of pulsar binaries in massive ST gravity and forms the key motivation for our study. From a purely numerical point of view, however, Eq. (21) creates a significant challenge. Numerical algorithms will pick up all possible modes of a solution -even if only through roundoff error.
We therefore seek an algorithm that provides us with explicit control over the asymptotic behaviour of our numerical solutions. In the next section, we will discuss how this can be achieved inside the more standard frameworks employed to address items (i) and (ii) of our above list.

Numerical framework
Numerical methods for solving two-point-boundary-value problems are well developed and fall into two main classes, shooting algorithms and relaxation schemes (including collocation methods) [74] To the best of our knowledge, all literature on static NS models in ST gravity has employed shooting algorithms; see e.g. [21,35]. This process integrates the differential equations from one end of the domain by supplementing the known boundary conditions at this point with appropriate trial values for the remaining variables. The resulting integration will typically not match the boundary conditions at the other end of the domain, but the degree of violation can be used, e.g. through a Newton-Raphson or a bisection method, to iteratively improve the trial values until all boundary conditions are satisfied within a user-specified threshold accuracy.
For the case of our system of differential equations (13)- (17) with boundary conditions (20), this would work as follows. We first note that the function Φ appears only in Eq. (13) and in the form of its spatial derivative. We can therefore set Φ(0) = 0 and add an arbitrary constant to match its boundary condition after having solved for all variables. Bearing in mind this freedom, we start the integration at the origin r = 0 by selecting values for the central baryon density ρ(0), the central scalar field amplitude ϕ(0) and the metric function Φ(0), additionally to the known η(0) = 0 and X(0) = 1/ F(ϕ(0)). The integration will reach zero pressure at a finite r which represents the NS radius. Beyond this radius, the integration continues setting ρ = P = 0 in Eqs. (13)- (17). In principle this surmounts the issue (ii) mentioned in the previous section. We note, though, that P = 0 is, in general, not realized on a grid point which adds a small discontinuity to the solution; the data on the outermost grid point inside the star and on the first point outside the star do not satisfy Eqs. (13)- (17). This discontinuity is typically not problematic, but we will see below how it can be simply eliminated in a relaxation approach. For a massless scalar field, it is even possible to analytically match the spacetime to an exterior vacuum 2 metric; cf. Eqs. (8), (9) in [18]. Integration beyond the neutron star radius is not required and the trial value ϕ(0) can be improved in accordance with the selected shooting algorithm.
Such an analytic matching is not known, however, for massive scalar fields. And now a more problematic issue arises as the integration is continued beyond the NS radius; no matter how accurate the central value ϕ(0) has been chosen, the numerical solution will contain an exponentially growing contribution from the asymptotic behaviour (21) and eventually blow up exponentially. Worse, this blowup prevents us from improving our trial value ϕ(0) through measuring the departure from the correct boundary condition at infinity; this departure is infinite and, hence, useless for numerical purposes. In shooting algorithms, this problem is circumvented by imposing the outer boundary conditions at a finite radius rather than infinity; cf. Sec. III A in [21].
While this method is still capable of generating accurate stellar models, a scheme covering the complete exterior and imposing the boundary conditions at infinity provides practical advantages besides the more rigorous boundary treatment. By extending all the way to infinity, our scheme can provide initial data for time evolutions on arbitrarily large computational domains (including compactified evolution schemes that incorporate spatial or null infinity) without resorting to adhoc procedures to extend results beyond the inevitably finite blow-up radius of shooting methods. We will also obtain a vacuum exterior that is matched to the NS interior on a grid point; in fact, the value of the NS surface, rather than the central density, will select the specific stellar model. Furthermore, the relaxation scheme provides an exceptionally elegant and simple way to implement the matching between the interior and exterior domain that we expect to be applicable to a wider range of problems, including extension to time evolutions.
For our method, we first introduce the NS radius R S as a free parameter. On the domain r ∈ [0, R S ], we use the differential equations (13)- (17) with boundary conditions (20) for η and X at r = 0; the condition 2 The term "vacuum" here refers to the baryonic matter; the scalar field is nonzero exterior to the star. FX 2 = 1 is formulated as an equation involving the unknown ϕ(0). In the exterior, we introduce a compactified radial coordinate set ρ = P = 0, and introduce rescaled variables By factoring the exponential dependence into our scalar field variables, we ensure that regular solutions σ and κ asymptote towards a ∝ y dependence at y = 0. We find this step crucial in achieving convergence of our relaxation scheme which struggles with the exponential fall-off of ϕ and η but copes smoothly with the benign linear behaviour of the rescaled σ and κ. We also notice a minor (but not crucial) improvement in the speed of convergence when switching from X to the mass function m of Eq. (11) and hence use the set of variables Φ, m, σ, κ in the exterior. The differential equations in the exterior domain y ∈ [0, y S ], and the matching conditions imposed at the surface of the NS are given by We formally also use the trivial equation ∂ y P = 0 in the exterior which allows us to use a constant number of five variables over the entire grid. This grid consists of N grid points in the interior and M points in the exterior. We discretize the differential equations using cell-centered second-order stencils which provides us with 5(N − 1) algebraic equations in the interior and 5(M − 1) equations in the exterior. The boundary conditions provide two further equations at r = 0 and three further equations at y = 0. The surface radius is represented twice on our grid, the outermost point r = R S of the interior and the innermost point y S of the exterior grid. The variables used on these points are related by the matching conditions (28) as well as the trivial P(R S ) = P(y S ) = 0. In total, we thus have 5(N + M) non-linear algebraic equations for the 5(N + M) unknown values of the variables on the grid points. Given an initial guess, we can linearize the equations around this trial solution which leads to a matrix equation with block-diagonal structure that is readily inverted to improve the guess iteratively; we use the algorithm of Ref. [74] and typically obtain convergence after about ten iterations. The initial guess is obtained by integrating Eqs. (13)- (17) up to R S , fixing σ and κ as linear functions ∝ y in the exterior and integrating Eqs. (24), (25) with these specified scalar sources. Note that in this calculation we set ρ = P = 0 in the exterior irrespective of whether or not they have reached zero at R S ; the discontinuity that may result at the matching point is removed in the ensuing relaxation process. Even for modest resolutions such as N = M = 401, this approach provides an accuracy of O(10 −4 ). All models discussed in the remainder of this paper have been computed with this code.

Overall phenomenology
We start this section by defining the terminology and diagnostic quantities as well as providing a qualitative review of the different branches of static NS models encountered in massive ST gravity. We then explore in the following subsections in more detail the impact of the ST parameters α 0 , β 0 and µ on the structure of these branches.
In the following, we will use the term "family" for the set of all NS models obtained for fixed EOS and ST parameters α 0 , β 0 and µ. We will use the term "branch" to denote a subset of solutions of a family that share some specific property, for example strong scalarization. A family thus consists of one or more branches. In some cases, we find a branch to have the shape of a closed loop disconnected from other branches, and we also refer to such a branch as a "loop". For reference, we note that the scalar mass µ introduces a Compton wavelength and characteristic frequency given by Unless stated otherwise, our numerical NS models in massive ST theory are computed with the polytropic EOS labelled "EOS1" in Ref. [56]. Translated into our notation, we therefore compute the pressure and specific internal energy from the baryon density ρ through Eqs. (18) and (19) with The families of solutions thus obtained are conveniently represented in a mass-radius diagram. For this purpose, we define the total baryon mass M b as the volume integral of the baryon number density multiplied by the mass per baryon m b . Translated into our baryon density ρ = m b n b , the expression becomes The motivation for using the baryon mass, rather than the gravitational mass of Eq. (11), arises from the conservation of the baryon number; if we consider the possibility that a NS might migrate dynamically from one branch to another, we expect it to do so at constant M b , whereas the binding energy and, hence, gravitational mass will, in general, change.
As in massless ST theories, all NS solutions can be classified as either weakly scalarized with scalar field profiles reaching a magnitude ϕ ∼ O(α 0 ) or strongly scalarized solutions where the scalar field reaches values ϕ ∼ O(1) [18]. In this work, we call these branches W (for weak) and S (for strong scalarization); see e. g. Fig. 1. The distinction between the two regimes naturally blurs for large values α 0 = O(1); in this work we consider only α 0 1 and thus retain a clear division between weakly and strongly scalarized NSs. As is well known, the NS solutions in GR 3 form a one-parameter family whose members can be characterized by its central density ρ c [75]. For the EOS (30), this leads to the black branch in the left panel 4 of Fig. 1. In ST gravity, the introduction of a scalar field leads to additional branches of NS solutions with a non-vanishing scalar-field profile ϕ(r). The shape of the extra branches depends on the values of the parameters α 0 , β 0 and µ. In agreement with the literature [48], we observe strongly scalarized solutions if β 0 −4.35. These new branches are displayed in the left panel of Fig. 1 in terms of a color code that denotes the central scalar field amplitude. In contrast, the weakly scalarized solutions we obtain for β 0 −4.35 have macroscopic properties that barely differ from those of the GR solutions and their branch would be indistinguishable from the GR family in the figure. For these cases, we have set α 0 = 10 −4 and µ = 4.8 × 10 −13 eV. We now discuss in more detail how the solutions and their branches vary when the parameter values are changed.

Dependence on µ
The variation of scalarized NS branches in massive ST theory has already been studied in Ref. [21] who generally observe that an increase in the scalar mass µ results in a weakening of the scalarization. By computing a sequence of models with equal gravitational or "Arnowitt-Deser-Misner" (ADM) [76] mass, they observe a monotonic decrease in the scalarization as the scalar mass is increased. Around 10 −12 eV, their scalar profile drops to negligible levels when β 0 = −4.5; cf. their Fig. 2. Our results exhibit a similar drop in scalarization. We illustrate this general behaviour by comparing the cases µ = 10 −15 eV and µ = 4.8 × 10 −13 eV in the right panel of Fig. 1; the color code of the branches represents the central scalar field amplitude and displays lower values for the larger µ. 3 We produce GR solutions with our code by simply setting α 0 = β 0 = 0. 4 The smooth curves in all our mass-radius plots are in fact made up of a large number of crosses which in most cases are not individually visible. We opt not to connect these with lines to avoid spurious cross-branch connections. In consequence, some curves appear to have breaks when the gradient becomes nearly vertical; these breaks are not physical. Whereas the S and W branches connect at two points when α 0 = 0, the S branch splits in two for nonzero α 0 with each part connecting to GR-like models in such a way that we obtain a "loop" of models separate from the main branch of solutions. We refer to the main branch as branch I and to the loop as branch I I.
We furthermore notice that the strongly scalarized NS branches reach out to smaller baryon mass and radii as we increase the scalar mass µ. As we shall discuss in more detail in Sec. 4.5, the stable NS model for a given baryon mass is that with the largest radius. For fixed M b , a larger scalar mass µ results in smaller and more compact stable NS models.

Dependence on α 0
When α 0 = 0, our system of equations (13)-(17) is invariant under the transformation 5 ϕ → −ϕ. In this case, each strongly scalarized model consists of two solutions that only differ by a minus sign in the scalar-field profile. Additionally to this degenerate scalarized branch, there exists a branch with zero scalar field, i.e. the set of models we also obtain in GR.
A nonzero α 0 breaks the degeneracy of the strongly scalarized branch which now splits into two branches with unequal macroscopic properties and whose scalar field magnitudes differ at a level O(α 0 ). This split is illustrated in Fig. 2 where we consider NS models in the M b − R S plane for µ = 4.8 × 10 −13 eV, β 0 = −4.5 and different value of α 0 . For α 0 = 0, we see that the scalarized branch directly connects to 5 Recall that ϕ → −ϕ implies η → −η and that F ,ϕ and V ,ϕ are linear in ϕ when α 0 = 0.   the GR branch. For α 0 = 0, we obtain a weakly scalarized branch W with ϕ = O(α 0 ) in place of the GR models. In terms of their mass M b and radius R S , however, these models are barely distinguishable from their GR counterparts, and we refer to them as "GR like" models. The strongly scalarized branch S, on the other hand, splits into two, each of them connecting to separate parts of the GR-like branch in such a way that we obtain one loop of models that is separated from the single, large branch; cf. the insets in the top right and bottom panels of Fig. 2. We find the gap between the loop and the main branch to be proportional to α 0 and independent of β 0 . For small but nonzero α 0 , we thus find two sets of solutions: A branch I that approximately follows the M b − R S curve of GR for small and for large central density ρ c , and a separate branch I I on a closed loop located in the region where strong scalarization occurs. The central baryon density strictly increases as we move along branch I, starting at (M b , R S ) = (0, 0). We note that branches I and I I both contain weakly as well as strongly scalarized models. Instead, their distinction arises from the separation of the loop from the main branch of models. As we increase α 0 , however, the loop of branch I I solutions shrinks and eventually disappears, leaving branch I as the only class of solutions. This remaining branch approximately overlaps with the GR family in the very high and low ρ c regime but shows a strong bulge of strongly scalarized solutions when ρ c has values comparable to nuclear density. In agreement with the literature, we observe that these scalarized models can reach significantly larger masses and radii than their GR counterparts. We illustrate these features in Fig. 3 for β 0 = −4.5 and −5, but note that this behaviour occurs universally in all cases we have studied.
We conclude the discussion of the α 0 dependence with a subtle observation we make throughout our computations: For all NS models of branch I I the central scalar field value has the same sign; ϕ c > 0 for our convention. For the vast majority (though not all) branch I solutions, ϕ c has the opposite sign; ϕ c < 0 in our case. We display this observation graphically in Fig. 4 for several combinations of β 0 and α 0 . We note, however, that branch I always contains a swap in sign(ϕ c ) at very large central baryon density: we always observe ϕ c > 0 as ρ c → ∞.

Dependence on β 0
Spontaneous scalarization is a non-linear phenomenon and driven by the quadratic coupling parameter when β 0 −4.35. It has already been remarked in [21], that this threshold value for strong scalarization is barely affected by the introduction of a non-zero scalar mass. Our results confirm this observation as is illustrated by the W and S branches of solutions shown in the left panel of Fig. 1 for β 0 = −4.5, −5 and −6 with fixed α 0 = 10 −4 and µ = 4.8 × 10 −13 eV. Note that the transition from weak to strong scalarization along a sequence of NS models, while strictly speaking continuous, is sufficiently abrupt to allow for a clear distinction between models belonging to branch W or S.
The three branches displayed in the left panel of Fig. 1 for β 0 = −4.5, −5 and −6 also demonstrate the increasing deviation in terms of mass and radius of branch S models from their GR-like counterparts. Increasingly negative values of β 0 allow for larger maximum mass and radius; cf. also Sec. IV A in [21]. This rather strong effect may provide opportunities for constraining β 0 through mass and radius measurement of NSs, although a reevaluation of the measurements in the framework of ST gravity (rather than assuming GR) will be required for this purpose. In Fig. 1, the strongly scalarized branch S models appear as an arc splitting off from the GR-like branch W. We always find branch S to have this qualitative shape and the size of the arc grows monotonically as β 0 takes on increasingly negative values. This is illustrated in Fig. 5, which displays branches W and S for several values of β 0 ≤ −4.5. This figure also demonstrates that branch S has a shape resembling an inverted 'S'; it initially splits off from branch W towards smaller radii (around R S ≈ 13.5 km and M b ≈ 1 M in the figure) before turning around and crossing branch W towards larger R S . Note also that for each choice of β 0 , branch S consists of two nearby but distinct curves; this splitting results from the relatively large value α 0 = 10 −2 as we have already seen in the last section; cf. the bottom right panel in Fig. 4.
For highly negative β 0 , we obtain NS models with yet larger radius and baryon mass as illustrated in Fig. 6, where we plot branches W and S for β 0 = −15, −17, −20, −25 and fixed α = 10 −4 , µ = 4.8 × 10 −13 eV. In this figure, we also notice a new effect: the upper (in the sense of larger ρ c ) end of branch S exhibits a more complex structure. Instead of connecting to branch W, branch S appears to remain separate and curl around; cf. the insets for β 0 = −15 and β 0 = −17. This behaviour becomes clearer for yet more negative β 0 : between β 0 = −20 and β 0 = −25, the intersection of branch S with branch W is lost and instead, branch S forms its own tail of NS models with very small central values of the scalarfield |ϕ c |; note the magenta color of this end of branch S. Contrary to what one might guess, the NS models on this tail of branch S are still strongly scalarized; the profile ϕ(r) merely reaches its maximum away from the center r = 0.
We explore this behaviour in more detail by comparing in Fig. 7 the sequence of models obtained for β 0 = −6 with that for β 0 = −17, keeping α 0 = 10 −4 and µ = 4.8 × 10 −13 eV fixed. The bottom panels show the respective families in the M b − R S diagram analogous to Fig. 6. Along the branches, we have marked several NSs by colored circles, and for these models we plot the baryon-density and scalar-field profiles ρ(r), ϕ(r) in the upper panels 6 . For β 0 = −6, we observe a simple pattern: At the lower branch point, we obtain a weakly scalarized model ("1") with comparatively small central baryon density. As we continue along branch S, the scalar field increases in strength, reaching a maximum at maximal radius (model "3"). Beyond that point, the central baryon density ρ c keeps increasing, but the scalar profile weakens. Eventually (model "6"), a weakly scalarized but highly compact NS marks the smooth reconnection with branch W; note that this weakly scalarized model is located on the unstable branch of the GR-like models. The analogous results for β 0 = −17 in the right column of Fig. 7 display a qualitatively similar behaviour near the lower branch point (model "1"); the central baryon-density and scalar-field values increase as we move along branch S (model "2"). Eventually, however, the scalar field profile changes its qualitative behaviour and peaks away from r = 0 while the central value ϕ c decreases. In consequence, the upper tail of branch S now consists of models with ϕ c ≈ 0 but strong scalarization at r > 0 and does not directly connect to branch W; compare model "8" for β 0 = −17 with model "6" for β 0 = −6. As branch S curls around, the scalar profile strengthens once again and we encounter models with a steep density cusp (model "9"). We cannot rigorously rule out that after further curling around, branch S might eventually connect with branch W, but our numerical results do not show any signs of this happening. We finally note the remarkable structure of these upper tail stars: A highly compact star of baryonic matter is surrounded by a shell of scalar-field (i.e. bosonic) matter. This structure is reminiscent of the atom like shape noticed for stars in massless ST gravity in [77] and scalarized black holes in modified gravity [78].

Stability of models
In the previous sections, we have seen many cases where for fixed ST parameters α 0 , β 0 , µ several equilibrium NS models with equal baryon mass exist; see e.g. M b = 2 M in the left panel of Fig. 1. We can analyze the stability of these models by comparing their binding energy. The model with the lowest ADM mass, i.e. strongest binding energy, represents the stable configuration and other models with equal . Several NS models are marked along the branches as colored circles. The top panels show the radial profiles of the baryon density ρ(r) and the scalar field ϕ(r) for these NSs using their respective color. The density profile always reaches a maximum at the origin; however, the scalar field profile in some cases reaches a peak at a non-zero radius.
baryon mass are expected to migrate to this configuration under perturbations. We note, however, that the physical relevance of this instability depends on the instability timescale as compared to other dynamical timescales under consideration.
Using this method, we classify in Fig. 8 stable and unstable NS models for several values of β 0 at fixed α = 10 −4 and µ = 4.8 × 10 −13 eV. The results confirm the theoretical prediction that the weakly scalarized branch W becomes unstable when strongly scalarized counterparts with equal baryon mass exist [48]. Note also that the scalarized branch S exhibits a stability structure analogous to the well known GR case: The maximum mass model separates stable from unstable stars and the stable models are those with larger radius. In Fig. 9, we analyze how the stable and unstable models spread among our "loop" branches I and I I of Sec. 4.3 for different values of α 0 . The stable NSs are the strongly scalarized models with the largest radius, whereas the NSs on branch I I (i.e. on the closed loop) are always unstable. As a general pattern in all our computations, we find the models with the strongest central scalar field value |ϕ c | to be the stable configurations. For our convention, these always turn out to be models with ϕ c < 0; for example compare Fig. 9 with Fig. 4.

Conclusions
In this study, we have numerically computed solutions of spherically symmetric NSs in massive ST theory using a numerical scheme that enables us to eliminate the exponentially growing modes from the scalar field. For this purpose, we split the domain into the NS interior and the exterior from the stellar surface to infinity and discretize the resulting equations with a second-order relaxation scheme. This method enables us to compute NS spacetimes extending all the way to infinity where we can prescribe the boundary conditions in simple Dirichlet form. This formalism also provides a trivially simple implementation of the matching conditions without the need to perform interpolation.
unstable stable α 0 = 3× 10 −2 , β 0 =−5.0 We have used the resulting code to compute solutions of static, spherically symmetric NSs in massive ST theory and explored in detail the structure of the resulting branches of solutions in the (baryon) mass-radius plane for combinations of the linear and quadratic coupling parameters α 0 , β 0 of the ST theory and the scalar mass µ. We summarize the main findings of our analysis as follows.
• In agreement with previous literature studies of NS equilibrium models in massive and massless ST gravity, we find larger values of α 0 and β 0 to result in larger deviations from the NS solutions in GR, whereas larger values of the scalar mass tend to reduce these deviations; cf. Figs. 1 and 3. • For α 0 = 0, the NS models of GR are also solutions of the field equations of massive ST gravity. For β −4.35, we find, additionally to the GR branch, the spontaneously scalarized class of NS solutions that Damour & Esposito-Farèse discovered in their original exploration of massless ST theory [18] and that were also identified in massive ST theory in [21]. These solutions are invariant under the scalar field transformation ϕ → −ϕ. • A non-zero α 0 breaks this degeneracy and results in a dissection of the branches around the branch points; instead of the two connected branches of scalarized and non-scalarized solutions for α 0 = 0, we now find a main branch I and a smaller loop of branch I I solutions; cf. Fig. 2. The solutions on branches I and I I are characterized by different signs of the central scalar-field value ϕ c ; cf. Fig. 4. • For sufficiently negative β 0 , roughly β 0 −15, we observe a qualitative change in the strongly scalarized branch S of solutions. Instead of smoothly approaching the weakly scalarized branch W as happens for milder β 0 , its upper (in the sense of increasing central baryon density) tail now either crosses or completely detaches from the W branch.
• For highly negative values of β 0 , we furthermore encounter a new type of strongly scalarized solutions at this upper end of the S branch: the maximum of the scalar field is located away from the stellar center; cf. Figs. 6, 7. In its most extreme form, these solutions are composed of highly compact NS models surrounded by a scalar shell; see e.g. [77,78] for similar "gravitational atom" like configurations in other theories of gravity. • Whenever multiple NS models with equal baryon mass exist, we find the scalarized model to be the stable configurations in the sense of minimal binding energy. Typically, though not always, this is the model with the largest radius; cf. Figs. 8, 9. We also observe that the stable configurations agree in the sign of the central scalar field value, ϕ c < 0 in our convention.
The behavior with respect to the scalar parameters seems to be universal as we have encountered the same M b − R S profile deviations with respect to GR for all other equations of state that we have studied. We have explored in a similar manner, though less exhaustively, the cold hybrid EOS1, EOS3 and EOSa [65], APR4 [79], 2H and HB [80] and observe qualitatively similar behaviour.