Galileons and strong gravity

In the context of a cubic Galileon model in which the Vainshtein mechanism suppresses the scalar field interactions with matter, we study low-density stars with slow rotation and static relativistic stars. We develop an expansion scheme to find approximated solutions inside the Vainshtein radius, and show that deviations from General Relativity (GR), while considering rotation, are also suppressed by the Vainshtein mechanism. In a quadratic coupling model, in which the scalarisation effect can significantly enhance deviations from GR in normal scalar tensor gravity, the Galileon term successfully suppresses the large deviations away from GR. Moreover, using a realistic equation of state, we construct solutions for a relativistic star, and show that deviations from GR are more suppressed for higher density objects. However, we found that the scalar field solution ceases to exist above a critical density, which roughly corresponds to the maximum mass of a neutron star. This indicates that, for a compact object described by a polytropic equation of state, the configuration that would collapse into a black hole cannot support a non-trivial scalar field.


Introduction
There is strong evidence that the present Universe is expanding in an accelerated way [1,2], and the standard model of cosmology presumes that a cosmological constant is responsible for this acceleration. The vacuum energy of the Universe would contribute to this cosmological constant, but a theoretical prediction for the vacuum energy is several orders of magnitude larger than the observed value to explain the acceleration. This theoretical complication has lead to alternative scenarios in order to explain the late time acceleration of the Universe, and also to solve the vacuum energy problem. As explained by Weinberg [3], one has to modify General Relativity (GR) in order to relax dynamically the value of the vacuum energy to smaller values. Following this line of thought, several infrared modified gravity theories have been proposed which, in principle, could achieve self-acceleration in the Universe due to gravitational degrees of freedom. However, there are stringent bounds on gravitational interactions which are different from GR in local (solar system) experiments. Therefore, for these alternatives to be physical, it is necessary to include a screening mechanism that "hides" the effects of such modifications to GR (see [4] for a review and references therein). One interesting mechanism, proposed by Vainshtein in 1972 [5], hides GR modifications due to non-linear derivative self-couplings. This mechanism can be embedded in the Galileon models [6,7].
The Galileon Lagrangian was initially introduced by Horndeski [8] few decades ago, as a subset of the most general scalar-tensor theory with second order equations of motion. This -1 -

JCAP10(2014)055
was rediscovered more recently by Nicolis et al. [6]. The Galileon Lagrangian has emerged as the effective theory of some theories such as the DGP model [9], the recently proposed ghostfree massive gravity [10,11], and vector theories that break in a specific way an Abelian symmetry [12][13][14]. Either within the embedding modified gravity theories or the simple Galileon Lagrangians, the Vainshtein mechanism has been intensively studied for static and spherically symmetric backgrounds in the weak field limit (see for example [15,[18][19][20][21][22][23]33]). However, little work has been done away from staticity and spherical symmetry, such as in two body systems [24][25][26] or cosmology [27][28][29][30][31][32][33]. In this work we try to give some insights into the implications of renouncing to staticity and to weak gravity approximations. For this purpose, we will consider the simplest Galileon model, which only includes the cubic Galileon term and find approximated solutions for slowly rotating low-density stars initially, and relativistic stars towards the end.
The paper is organised as follows. In section II, we present the model and the equations for spherically symmetric vacuum configurations. We first review the solution for a scalar field in the Minkowski background and develop a perturbation scheme that works inside the Vainshtein radius. We apply this method to find first order corrections to the vacuum solutions in GR due to the Galileon term. We also discuss the effect of adding rotation. In section III, we introduce a coupling to matter and discuss how the solution with a constant density source matches to these vacuum solutions. In section IV, we study the scalarisation phenomenon (strong amplification of the effective scalar-matter coupling) in a quadratic coupling model, and study whether this happens in the presence of the cubic Galileon term. We then study more realistic neutron stars and investigate how the Vainshtein mechanism suppresses the modification of gravity in the strong gravity regime. We also discuss the condition for existence of the scalar field solution.

Minimally coupled cubic Galileon action and equations of motion
Our starting point is the Einstein-Hilbert action, minimally coupled to a scalar field with the cubic Galileon term with no additional matter included where α is a dimensionless constant, M p is the reduced Planck mass and Λ is a constant with dimensions of mass. In models where the cubic Galileon term is associated with the late time acceleration of the Universe, Λ 3 α −2 is typically given by (1000 km) −3 . Varying this action with respect to the metric and the scalar field, we obtain the following equations of motion: One can understand the appearance of the Ricci tensor in (2.2b) from the commutation of the third order derivatives of the scalar field [∇ α , ∇ ν ]∇ β Φ = R ρ βαν ∇ ρ Φ.

JCAP10(2014)055
Let us consider a static and spherically symmetric spacetime, together with a scalar field that depends on the radial coordinate only, namely ds 2 = −e 2ν(r) dt 2 + e 2λ(r) dr 2 + r 2 (dθ 2 + sin 2 θdϕ 2 ), Φ = Φ(r). (2.3a) Later on we will depart from this assumption by considering a slowly rotating spacetime (with an off-diagonal dtdϕ term). Under the previous assumption, the equations of motion (2.2) then reduce to The equation for the scalar field (2.4d) can be integrated once, as can be expected from the fact that the Lagrangian does not depend on Φ itself but on its derivatives only. The integration of this equation results in e ν−λ r 2 Λ 3 Φ + 2α 2 re ν−3λ (2 + rν )Φ 2 = Λ 3 ζ, (2.5) where the dimensionless integration constant ζ must be determined by boundary conditions. This is an algebraic equation for Φ , with solutions Φ = −e ν−λ r 2 Λ 3 ± e 2(ν−λ) r 4 Λ 6 + 8α 2 Λ 3 ζre ν−3λ (2 + rν ) 4α 2 re ν−3λ (2 + rν ) . (2.6) From this expression, one notices that in the limit of large α (keeping all other quantities fixed) Φ scales as 1/α, plus subleading corrections. When substituting the scalar field solution into the Einstein equations, one finds that the Galileon contributions are suppressed by powers of 1/α, hence they can be neglected in the large α limit. This result suggests that we can analyse this Einstein-Galileon system, in the large α limit, by a perturbative expansion in powers of 1/α. This is the approach we will follow.

A test scalar field
Before studying the full system, it is illustrative to look at the behaviour of a test scalar field on a fixed background. If we consider the Minkowski spacetime as our background, the solution for Φ reduces to -3 -

JCAP10(2014)055
where we have chosen the positive sign in (2.6) so that, in the limit r → ∞, we recover Φ flat = 0 modulus terms with O(r −2 ). In the large α limit, we can expand the previous solution to obtain This expansion is reliable only if the magnitude of every successive term is smaller than the previous one. Since the second term in the expansion arises from the kinetic term, the condition that the first term to dominate over the second term is equivalent to the condition that the contribution from the Galileon term is dominant with respect to the canonical kinetic term. The equality of the first two terms defines the 'Vainsthein' radius, For r < r V , both the validity of the series expansion and the dominance of the Galileon term over the kinetic term are guaranteed. Furthermore, the ratio of the solution (2.8) to that without the Galileon term (Φ can ∼ r −2 ) is given by This suggests that (2.8) can be written as Thus the α −1 expansion is equivalent to the expansion in terms of (r/r V ) 3/2 in this case. In the following, we will apply this expansion scheme to find a solution deep inside the Vainshtein radius r < r V . It was argued that in this regime, there exists a dual description of the theory that admits perturbative solutions [40,41]. If applied to a test scalar field in the Minkowski spacetime, the perturbative expansion in the dual formulation is controlled by (r/r V ) 3/2 , and thus agrees with our expansion. Although we will study the coupling to matter in the next section, in order to gain some insights into the integration constant ζ, it is worthwhile a brief digression to look at the scalar field solution in the presence of a test point source, coupled to the scalar field around a flat background. The point particle may be introduced in the original Lagrangian by adding the mass term The scalar field equation then becomes

JCAP10(2014)055
As a conclusion, the scalar field charge, which is proportional to ζ, plays essentially the same role as a mass. To close this digression, we also note that under this identification of ζ, the Vainshtein radius r V reads which agrees with the expressions in the literature, up to a redefinition of β, for the Vainshtein radius (see for example [44]).

Large α expansion
We now turn our attention to the study of the solutions in the full system using the large α expansion. At the leading order, we expect to recover GR solutions due to the Vainshtein mechanism. For a spherically symmetric system in GR, the solution is uniquely determined by the Schwarzschild metric. However, there is no Birkhoff's theorem in our system thus we cannot exclude a possibility to have solutions different from Schwarzschild. Here we are interested in a solution which is close to Schwarzschild within the Vainshtein radius to recover GR, hence we would like to see if it can can be consistently constructed using our large α expansion.
At the leading order, this is equivalent to consider a test Galileon field, which accretes into the Schwarzschild black hole, as it was done in [45]. Then, one could analyse the backreaction of such a field (see [46]). Here, we will calculate the metric corrections using our large α approximation. This approach should give some insights into the Vainshtein mechanism under strong gravitational fields. Unlike [45,46] where the existence of a black hole surrounded by a Galileon scalar field is discussed in the presence of non-trivial cosmological boundary conditions, here we are interested in finding the interior and exterior solutions of a massive body in an asymptotically flat and static spacetime. Therefore, we begin by constructing the vacuum exterior solution in this section. Later on, we will study the interior solution, and also discuss what happens when the mass of such a configuration approaches the limit where it might collapse and form a black hole.
In [44], the authors offered a demonstration of the Vainshtein mechanism around a thin-shell configuration of matter, showing that asymptotically flat solutions exist. Here, under the large α approximation, we construct explicitly solutions around a static matter configuration, and confirm that below the Vainshtein radius the behaviour of the scalar field is determined by the Galileon term. We also provide the explicit matching between the interior solution for a static source of constant density and the exterior Galileon solution. The transition from this solution to the asymptotically flat solution beyond r V is not done explicitly, but it is guaranteed by the sign choice made for Φ.
Our method is a generalization of what we did for a test scalar field in Minkowski spacetime. First we solve for the scalar field on a Schwarzschild background metric and then we look for the first order corrections to the metric in the α −1 expansion. We expand our functions in powers of α −1 as where r s is the Schwarzschild radius associated with a mass M . The first term of the scalar field expansion, Φ 1 , is fully determined by ν 0 and λ 0 in the scalar field equation. Altogether, ν 0 , λ 0 and Φ 1 determine ν 1 and λ 1 via the terms of order α −1 in the Einstein equations (2.4). After inserting these functions in the corresponding terms of the same order in the scalar field expansion (2.5), one can solve for Φ 2 . Finally, we can identify the contribution from the canonical kinetic term in Φ, and obtain an expression for the 'Vainshtein radius' as the radius at which the Galileon contribution becomes comparable to the contribution from the canonical kinetic term. Either by expanding (2.6), or by directly solving the scalar field equation (2.5) at the lowest order in α −1 , we find that the solution for Φ 1 is where ζ 0 has the same interpretation as before: the scalar charge developed by Φ on a given background metric. In the limit r s → 0, equation ( 16) and solve this equation at each order in the α −1 expansion. This shows that the total scalar charge is still a conserved quantity but it is artificially split into parts coming from different orders in the expansion (2.14). For example, ζ 0 in (2.15) has a contribution from the Schwarzschild part of the metric only. We now insert ν 0 , λ 0 and Φ 1 in the Einstein equations, obtaining two independent equations for λ 1 and ν 1 ; These equations can be solved exactly, giving us the dominant Galileon corrections to the Schwarzschild metric, which could then be used to explore strong gravity effects. The solutions have rather complicated expressions, hence we only display explicitly λ 1 but leave ν 1 unintegrated. The solutions are given by where λ 10 (and ν 10 that we will meet next) are integration constants associated with eqs. (2.17), 1 − k 2 sin 2 θdθ is an elliptic integral of the second kind. For future use, we present the weak gravity limit (r s /r 1) for λ 1 and ν 1 , Now we can find the solution for Φ 2 using (2.16). The resulting solution is exact, however the expression is too lengthy and it is sufficient to stress that the metric and scalar field solutions are all consistent with the Minkowski background solutions that we found in the previous section in the weak field limit r s /r 1. Moreover, in order to estimate the Vainshtein radius, we only need the contribution from the canonical kinetic term to Φ 2 , which is given by This expression should be compared to Φ 1 in order to obtain the radius at which the Galileon cubic term becomes of the same order as the canonical kinetic term. In other words, the following equality should be satisfied which reduces to Notice that the Vainshtein radius in the Schwarzschild spacetime is always larger than that in the flat background, (2.9), because the factor in front of r 3 is always less than one. This last claim however assumes that the integration constant ζ 0 has the same value both on a Schwarzschild and on a Minkowski background. In reality they may differ because they are given by matching different scalar field solutions for each background to the solutions inside the matter configuration under consideration. It is important to stress that in order to derive this results we did not use any weak gravity approximation, but only neglected contributions to Φ 2 coming from λ 1 and ν 1 . However, from the equations (2.17) it is possible to deduce that these metric corrections are, at most, of order α −1 , as in the Minkowski case.

Rotation
In the previous section, we obtained static spherically symmetric solutions using our α −1 expansion without taking weak gravity limit. In this section, we will include rotation into the solution. It is difficult to find solutions with fast rotations analytically, so we focus on a slow rotation component only. Let us consider a slowly rotating spacetime characterized by a small function ω(r). To first order in this function, the metric is just the same as in the static case but with a new dϕ dt component in the metric, ds 2 = −e 2ν(r) dt 2 + e 2λ(r) dr 2 + r 2 (dθ 2 + sin 2 θdϕ 2 ) + 2r 2 ω(r) sin 2 θdϕdt.

JCAP10(2014)055
This additional metric component is determined by a new equation of motion, ξ ϕ t , which, because of ∂ ϕ Φ = 0, does not have explicit contributions from the scalar field and is therefore the same as in GR. In vacuum this equation can be integrated once and the resulting integration constant is just the angular momentum, (2.25) Expanding ν and λ as in (2.14), it is then natural to split ω in an O(α 0 ) part, which is just the usual Schwarzschild solution, and an O(α −1 ) part, which has to be determined from the metric corrections ν 1 and λ 1 given by (2.18). Also, the total angular momentum J is separated into contributions coming from different orders in α, J = n=0 α −n J n . Then where Ω 0 is the integration constant that appears when solving for ω in (2.25) to the lowest order in α −1 . Inserting the above expansions in (2.25) and using (2.18) it is straightforward to read off an expression for ω 1 ; resulting in which can be integrated once to obtain ω 1 explicitly. For the discussion, it is sufficient the weak gravity limit result, given by where Ω 10 and ν 10 are the integration constants for ω 1 and ν 1 , respectively. It might appear that for small r the corrections associated to ζ 0 become large, however we need to remember two things: i) this is valid only for r > r s , ii) ζ 0 is proportional to r s in the weak gravity limit (see (2.12)). In this limit, the term proportional to ζ 0 is suppressed compared to the J 0 /r 3 term from GR, by a factor (rM p ) −1/2 (α −2/3 Λr s ) 3/2 . Thus the contribution from the cubic Galileon term to the rotation is highly suppressed.

Coupling to matter
In the previous section we derived the vacuum solutions using the large α expansion. The integration constants in the solutions need to be determined by a coupling to matter. We include the matter coupling in the following way S M represents the action for some matter field Ψ m and A 2 (Φ) is the conformal factor relating the metric in the Jordan and Einstein frames, g J µν = A 2 (Φ)g µν . This conformal factor converts -8 -

JCAP10(2014)055
the non-minimal coupling between Φ and the curvature R into a Φ-dependent modification of the geodesic equations derived from (3.1); this is usually referred as a fifth-force. Actually for the simplest case of 2 ln A = −βΦ/M p , which corresponds to the Brans-Dicke theory when the Galileon term is absent, the Jordan frame version of (3.1) is given by where the pertinent redefinition of the scalar field and conformal transformation are In the next section, we will also consider the case of ln A ∼ Φ 2 /M 2 p , which, for a normal scalar field, provides an interesting example in which the effective coupling between the scalar field and matter can be dramatically amplified in the interior of high density stars, due to the so-called scalarisation phenomenon [35]. We will focus on these two choices of A(Φ) and refer them as the linearly and quadratically coupled models, respectively. Moreover, we work exclusively in the Einstein frame.
The equations of motion, derived from (3.1), with respect to Φ and g µν are We have defined the energy momentum tensor T µν (with trace T ) and the effective coupling strength between the scalar field and matter, γ(Φ), as The quantity γ traditionally plays a relevant role in constraining scalar-tensor theories within the Solar system using the parameterized post-Newtonian formalism (see for example [47]). However, the first post-Newtonian parameters depend on the asymptotic behaviour of the scalar field, which is uknown in the context of the Galileon model, because we do not know the exact exterior solution for the scalar field. The existence of an intermediate Vainshtein regime between the star and the asymptotically flat spacetime poses some difficulties in the use of asymptotic values to constrain the model (3.1). In the next section, we solve the equations (3.2) for a matter source described as a perfect fluid. First, we consider a simple model characterized by a constant and low density, and then a more realistic description of matter using a polytropic equation of state. For the first case, it is sufficient to consider the linearly coupled model since, as shown in [37], in weak gravity any deviation from GR is sensitive only to the cosmological value of γ, which is precisely the field-independent coupling strength of the linear model. However, for the second case we consider the quadratically coupled model to allow for the presence of strong gravity deviations from GR, analogous to the ones reported in [35,36]. We focus on the first case next.

Linearly coupled weak gravity model
Let us consider an incompressible source of matter described by a perfect fluid with pressure P (r) and constant density ρ = ρ 0 . Its energy momentum tensor is given by T µν = (ρ 0 + P )u µ u ν + P g µν , where the 4-velocity u µ of the fluid satisfies u µ u µ = −1. We consider the coupling function given by ln A 2 (Φ) = M −1 p βΦ. At the lowest order in the large α expansion, the equations (3.2) reduce exactly to those of GR, whose solution for the interior of a static and spherically symmetric distribution of matter with an exterior Schwarzschild spacetime is known as the Tolman-Oppenheimer-Volkoff (TOV) solution where the radius of the star, R, is defined as the radius where the pressure P T OV vanishes, and ν 0 is an integration constant used to match the metric with Minkowski at r = 0. Moreover, we have already used the matching at R between the radial component of the Schwarzschild exterior solution, 1 − rs r −1 , and the radial component of the solution above. This matching is given by An incompressible equation of state is one of the few cases in which the equations of motion admit an analytic solution, and despite of its simplicity, it illustrates some general features of more realistic solutions. In particular, one finds a lower bound for the massradius ratio of the star by requiring that the central pressure remains finite, as can be seen from (3.3b) by taking r = 0.
It is also possible to add a slow rotation to the TOV solution by adding the term 2r 2 ω(r) sin 2 θdϕdt to the line element (3.3a). We will do this in detail later, but we state in advance that the matching conditions determine the scalar charge and the angular momentum as follows: The first equation is obtained by integrating once the Euler-Lagrange equation for the scalar field inside the star and matching it with the scalar charge seen by the exterior solution.
In the second equation, J is the integration constant from the exterior solution as defined in (2.25). Note that these equations for J and ζ are valid at any order in α −1 .
Up to now, the only effect of the Galileon interaction is via the scalar field suppression inside the Vainshtein radius, which allows us to recover the GR solution for the metric at the -10 -

JCAP10(2014)055
lowest order in α −1 . In the next subsection, we show how the interior/exterior matching can be done consistently when the first α −1 corrections are included.

Corrections to the TOV solution
Inserting the TOV solution in the scalar field equation of motion and using the expansion in terms of α −1 , we can solve the scalar field equation to the next order. In order to include the first order Galileon corrections, we assume that the pressure has the following form and the metric is given by In this section, we will consider a low density star and take the weak gravity limit. We assume that the presence of the scalar field does not modify the incompressible nature of the star, so ρ = const.. Expanding the scalar field equation of motion (3.2) with Φ = α −1 Φ 1 (r) and the metric (3.7) we can solve for Φ 1 , Just as for the exterior scalar field there are two branches of solutions. We pick the positive root since we want to match the solution with the positive branch of the scalar field outside the star to ensure that the asymptotic flat solution exists. Performing the same expansion in the metric equations (3.2), we get a set of the equations for the metric and pressure corrections, which yield In (3.9) we have already fixed the integration constants for the interior solution by imposing regularity at r = 0 and vanishing P 1 at r = R, i.e., we impose the star has the same radius as it would be without the Galileon interactions. This means that the contributions from the scalar field to the pressure must change the density. It is convenient to write down explicitly the exterior solution with all the integration constants that we need to fix at the surface of the star. As we are considering a low density star we can take the limit r s /r 1. The metric -11 -

JCAP10(2014)055
and scalar field outside the star are given by We use the label ext to the integration constants of the exterior solutions (2.15) and (2.19). We note that once the Galileon interaction is turned on, there is no reason for the density of the star to be the same as in the GR solution. Therefore, it is convenient to split the energy density as where ρ 1 is the contribution to the density coming from the Galileon interaction and ρ 0 is equal to the density of the pure GR solution. Also note that ν ext 0 is just a normalization factor related to ν 0 in (3.3a), which is already fixed so that the interior metric becomes Minkowski in the limit r → 0, e 2Ψ r=0 = 1. Therefore, we are left with four free parameters r s , ζ 0 , λ ext 10 and ν ext 10 , which are determined by the matching conditions Φ 1 = Φ 1,ext , g int rr = g ext rr , g int tt = g ext tt and g int tt = g ext rtt , at r = R. The results are The correction to the mass due to the Galileon interaction is explicitly considered by the way in which we split the density ρ and this is reflected in the expression for r s above, where we see that the Schwarzschild radius is corrected (in pure GR, r s would be related only to ρ 0 ) even though λ ext 10 = 0. The matching for ζ 0 , which is performed directly from the scalar field solutions (3.8) and (3.11) here, confirms the result (3.5a) and naturally extends the interpretation of ζ 0 given before in terms of a point mass living in a Minkowski space-time. Indeed one can check that at the lowest order in r s /r, the result (2.13) for the Vainshtein radius in terms of the point mass is recovered. From (2.13) and the matching for r s given above we can estimate the condition for the TOV metric to be dominant over the corrections (3.9). In particular for the radial components we have (3.14) -12 -

JCAP10(2014)055
The left hand side can be at most of order 1/R while the right hand side is much larger than 1/r s as r V /r s 1, thus (3.14) is easily fulfilled. This closes our study of low density static solutions. The Vainshtein mechanism works well inside the matter configuration and suppresses the scalar field effects in order to recover the GR solution with a negligible backreaction from the Galileon term, and this solution can be matched smoothly to the exterior metric -the Schwarzschild solution with small Galileon corrections.

Corrections to rotation
Next we consider the solution for rotation in the weak field limit. Inside the star the equation of motion, 2e 2λ r(ρ + P )ω + M 2 has the solution where one of the integration constants has been fixed to zero in order to obtain a regular solution at r = 0 and Ω 0 represents a slowly rotating frame of reference with respect to the Minkowski spacetime at r = 0 obtained by the coordinate transformation t = t, r = r, θ = θ, ϕ = ϕ − Ω 0 t. The transition to the exterior slowly rotating solution, requires ω = ω ext and ω = ω ext at r = R, which in turn implies for an arbitrary Ω 0 = 0. Now we consider the correction to this solution due to the Galileon term. Fro this, consider the metric g µν dx µ dx ν = − e 2(Ψ+ 1 α ν 1) dt 2 + e 2(Ξ+ 1 α λ 1) dr 2 + r 2 (dθ 2 + sin 2 θdϕ 2 ) + 2r 2 ω(r) + 1 α ω 1 (r) sin 2 θdϕdt, (3.19) together with the splitting (3.6) and (3.12) for the density and pressure, respectively. The first perturbative term for the rotation, α −1 ω 1 , is completely determined by the corrections to the diagonal components of the metric and to the pressure (3.9), and by requiring that the solution is regular at the centre r → 0. As before, the scalar field contributes indirectly through λ 1 and ν 1 . Expanding (3.15) up to O(α −1 ) and solving for ω 1 , we obtain (3.20) On the other hand, the Galileon correction for the small rotation of the exterior solution is given by (2.28), and its integration constants are to be fixed by requiring g tϕ = g ext tϕ , -13 -

JCAP10(2014)055
(g tϕ ) = (g ext tϕ ) at the surface of the star, r = R. This is done straightforwardly and the result is The condition to guarantee α −1 ω 1 ω 0 is the same as (3.14). Thus for a slow rotation, the Vainshtein mechanism operates successfully to suppress deviations from GR.

Quadratically coupled model
In the previous section we studied a matter coupling characterized by ln A 2 (Φ) = M −1 p βΦ, or equivalently by a constant coupling strength 2γ = M −1 p β, which corresponds to a Brans-Dicke theory when the Galileon term is absent in the action. Now we go to the next case, that of a quadratic coupling, leading to an effective coupling strength γ ∼ βΦ/M 2 p with a dimensionless β. It has been shown in [35,36] that this is the simplest case where a configuration of matter, a neutron star in particular, can develop a scalar field large enough to produce significant deviations from GR in the strong gravity regime. This effect, called spontaneous scalarization, is present independently from the rotation. For this reason, we limit ourselves to the study of static solutions to the equations (3.2) in the presence of a matter source described by a perfect fluid with a polytropic equation of state to model high density stars.

Toy model for scalarization in standard scalar-tensor gravity
Here we follow [35] to motivate the idea that a Φ−dependent effective coupling strength may lead to a peculiar behaviour of the scalar field inside a matter source. We only deal with the scalar field equation of motion, (3.2) with α = 0. Basically we are studying a test scalar field on a flat space-time in the presence of a matter source whose coupling is determined by β and the scalar field itself. Then we need to solve on a Minkowski background. Furthermore we assume where R and M are the radius and (ADM) mass of the star and the quantity s defined by the last equality as s = r s /R is the self-gravity of the star. In the weak gravity limit one can expand the effective coupling strength schematically as γ = γ 0 (1 + a 1 s + a 2 s 2 ), where a i are finite parameters. This would suggest that even when the self gravity of a star is large, the effective coupling γ remains small if the asymptotic γ 0 is small to pass the post-Newtonian constraints. However this result is obtained perturbatively in s and the toy model summarized here shows that it does not necessarily hold in the strong gravity regime.

JCAP10(2014)055
The scalar field profiles inside and outside a star are given by the solutions to The choice of sign for β determines, a posteriori, whether the scalarization phenomena takes place (β < 0) or not. Assuming β < 0, and after imposing the regularity at r = 0, the solution for r < R is given by where the central value of the scalar field, φ c , is determined by matching both Φ and its derivative with the exterior solution obtained from (4.3), Φ ext = −r −1 φ s + φ 0 , at the surface of the star. This gives a relation between φ c and the asymptotic value φ 0 , For non-zero β and s, |φ c | > |φ 0 |. In fact, if 3|β|s = π/2 the drastic amplification of the central scalar field happens and therefore the local coupling strength, |γ| = |β|φ c /M 2 p , is significantly large even if γ 0 in the perturbative expansion above is vanishingly small. Evidently this a simplistic model and the amplification may not be equally strong in a more realistic model, but the effect would still be present, enhancing deviations from GR in the strong gravity regime.

Scalarization -simplified Galileon model
Now, we explore the same toy model as before, but taking into account the Galileon term in the equation of motion for Φ, namely For the Minkowski metric the R µν term does not contribute. If we employ our α −1 expansion to solve this equation, we find with φ c = const.. At order α 0 , we are left only with the second and third terms while Φ is pushed to O(α −1 ), effectively suppressing the canonical kinetic part of Φ. The solutions for Φ 1 and Φ 2 will give rise to some additive integrations constants. However, it is important that φ c = 0, to keep a contribution from T at the lowest order in α −1 , since such a term is necessary to get a regular solution as r → 0. Moreover, the constant φ c gives the coupling strength γ = βφ c /M 2 p as before. Explicitly, the solutions obtained are as follows: • For the interior region, r < R, equation (4.6) gives whereβ = −sign(β)3|β|s is a dimensionless O(1) parameter. When solving for Φ 1 and Φ 2 , and for their respective integration constants by imposing the regularity of the scalar field at r = 0, we find The first order solution Φ 1 is similar to the first term of Φ 1 in the linearly coupled model (3.8). This behaviour suggests that the scalarization phenomena does not occur in this model. The second order solution Φ 2 is independent of the choice of sign for Φ 1 .
The total scalar field inside the star can be written as where the integration constants coming from integrating (4.9) have been set to zero. This is our interior solution.
• For r > R we are simply in vacuum and (4.6) is a quadratic equation that can be solved for Φ without the α −1 expansion, leading to We take only the upper sign since it corresponds to the asymptotically flat solution. Now we further split this solution into the two regimes separated by the Vainshtein radius.
i) r r V : here, and specially near the surface of the star, the solution is dominated by the Galileon contributions, which after using the α −1 expansion, reads This is the solution in the Vainshtein regime. As we approach to r V the scalar field dynamics coming from the canonical kinetic term becomes relevant and the solution must be connected to the asymptotically flat solution. However to see this within the α −1 series it is necessary to take several more terms in the expansion.
i) r > r V : far away from the source and beyond the Vainshtein radius, we have the asymptotic behaviour typical of a standard Scalar-Tensor theory, given by This is the exterior solution.
We want to relate φ c in the interior solution, to the value of Φ at R deep inside the Vainshtein radius and then to φ 0 in the exterior solution. The last step cannot be done easily in general as our large α approximation is valid only for r < r V and the asymptotic -16 -

JCAP10(2014)055
solution is valid for r > r V . Formally there is no overlaping region where the matching can be performed. For the particular case of Minkowski, φ 0 can be related to φ c using the exact solution. However once we consider the Schwarzschild exterior metric, this is no longer possible. Therefore, we will match the α −1 expanded solution (4.12) with the asymptotic solution (4.13) at r = r V , and asses the inaccuracies caused by this matching using the exact solution.
We impose the following matching conditions at r = r V and r = R on the solutions Φ I , Φ V and Φ E .
where the indices denote the region of validity of each solution: I for the interior, V for the Vainshtein regime and E for the exterior. From the first two equations we relate φ 0 and φ V at the Vainshtein radius, and from the third and fourth equations we relate φ V to φ c at the radius of the star, Putting these results together we have our final relationship between φ c and φ V , that basically depends on the difference Φ V | r V − Φ V | R plus the following explicit contribution from the selfgravity of the star Unlike the relation (4.5), the amplification of the central scalar field is not significant in this setup. Furthermore, the term proportional to r 2 V is always the dominant one, thus the difference between the scalar field at the centre of the star and its asymptotic value does not depend explicitly on the properties of the star. Since Λ 3 /α 2 ∼ ζ 0 r −3 V ∼ M p r s r −3 V , the dimensionless quantity M −1 p (φ c −φ 0 ) is of order r s /r V 1, as can be confirmed from figure 1. Therefore, the scalarisation does not happen in this case. As mentioned earlier, the matching at the Vainshtein radius is not well defined, but this is necessary since we cannot compute the solutions for a general metric without the large α 2 approximation. However for the Minkowski metric we do know the exact solution and therefore we can quantify the error due to this matching. If we perform the matching between the exact exterior solution given by (4.11) and the approximate interior solution (4.10), we find where Φ refers to the exact exterior solution (4.11).
To get an idea of the error's magnitude, we consider the values Λ/α 2/3 = (1000 Km) −1 , R = 10 8 m and r V = 10 16 m, which fix ζ 0 and also the Schwarzschild radius of the object through (2.12). The error in the relation between φ c and φ 0 is then given by This relative error tells us that when we use the matching at r V , we can trust the order of magnitudes of our results but not in the precise values. Figure 1 shows a comparison between the approximated solution and exact solution.

Compact stars
In this section we solve the equations (3.2) under the α −1 expansion with a quadratic coupling for a stellar model consisting of a perfect fluid T µν = (ρ(r) + P (r))u µ u ν + P (r)g µν obeying a polytropic equation of state. The density and pressure are parametrised in terms of a dimensionless function χ(r), a polytropic exponent Γ, a polytropic constant K, and the formula From now on we stop using natural units, in order to give values of the densities in a way that is standard in the literature discussing compact stars in GR. For a polytropic exponent 2 Γ 3 the observational mass-radius curves of neutron stars are well reproduced within GR. For smaller values of Γ (around 5/3), the solution corresponds to lower density stars, commonly defined by c 2 n 0 m b χ < 5 × 10 14 gr/cm 3 ≡ρ 0 . For these low density stars, the value of K can be estimated from the non-relativistic limit of a Fermi gas model, which yields a polytropic equation of state with Γ low = 5/3 and a polytropic constant given by Compact stars can be modelled in a more sophisticated way by considering them as having both a low density and a high density region, i.e. by taking a two-component polytrope.
In this case the continuity of the pressure at the interface between these two components requires that for the high density region Although we consider a single component star of high density, the previous relation is useful to estimate the appropriate K for a given Γ. We use Γ = 2.34 because with this choice, the maximum mass of a neutron star, shown in the results below, is consistent with the maximum mass computed with a more realistic equation of state in agreement with observational data, e.g., EOS II in [39].

Equations of motion
In order to easily compare with the analysis of [36], it is convenient to write the metric explicitly in terms of a Schwarzschild-like radial component by redefining the O(α 0 ) part of the metric as e −2λ 0 = 1 − 2µ(r)/r. Since the exterior metric solution asymptotes a flat spacetime, µ(r → ∞) can be interpreted as the ADM mass. Up to O(α −1 ), the metric is written as We also redefine Φ to make it dimensionless, Φ → M p Φ, and write the pressure and density in terms of their Jordan frame counterparts, i.e., (ρ, p) → A(Φ) 4 (ρ, p) where A(Φ) = exp(βΦ 2 /2). As before, we expand the scalar field as Φ(r) = φ c + α −1 Φ 1 (r). Inserting the above expansions and redefinitions in the equations of motion (3.2), we got The last equation can be obtained either as a combination of the Einstein equations or directly from the conservation of the energy momentum tensor. All the dynamical effects of -19 -

JCAP10(2014)055
the scalar field are removed, and its only contribution is a constant rescaling of the density at the lowest order. As we saw in the toy model of the previous section, φ c is required to a obtain regular solution for the scalar field as r → 0, which introduces an effective Newton's constant given by e 2φ 2 c β G. The first order correction to the scalar field Φ 1 is determined again using (3.2), and results in Just like in the Minkowski case, the effective coupling strength between the scalar field and matter is constant, and given by φ c β. We approximate the initial conditions for µ and ν at r = r min ≈ 0 by taking the functions on the r.h.s. of (5.6) as constants over an infinitesimal interval r min + δr, and then we integrate with respect to r to obtain where all the subscripts c stand for quantities evaluated at r min . The same procedure for the scalar field requires that Φ 1 (r = 0) = 0. However, in order to avoid numerical singularities we have to set Φ 1 (r min ) ∼ 0, and in order to keep consistency with the large α approximation, we also need Φ 1 (r min ) φ c . We vary the initial condition for the density, χ c , in the range from 0.1 to 50, in units of c 2 /10m b n 0ρ0 . These choices are of order one below and above the typical density for a neutron star. Notice that the initial condition φ c can be different for each density, however and for a preliminary analysis, we keep it the same for all the densities, with an arbitrarily chosen value of φ c ∼ 0.05. Later, we will take into account the variation of φ c in order to match a given boundary condition (for example an asymptotic value of the scalar field at infinity). In addition, the possible variation of the initial condition for Φ 1 is neglected since it can be absorbed in φ c . Furthermore, to fix a value for the Galileon coupling constant we make use of the weak gravity estimate of Λ 3 /α 2 = M p H 2 0 , where H 0 is the current Hubble parameter. This estimate comes from the requirement that the Galileon modifications are relevant for the present day expansion of the Universe, and at the same time they are consistent with constraints on GR deviations within the Solar System [6]. This means Λ 3 /α 2 ∼ (1000 Km) −3 , giving a Vainshtein radius of the order of 10 18 m for the Sun, as can be computed from (2.11). For comparison, the aphelion of Neptune's orbit is of order 10 12 m, therefore the Solar System is well inside the region where the nonlinear Galileon dynamics takes place. From now on, we set α = 1 and by doing this we are treating Λ as a meaningful physical parameter. Note that, as shown after (2.10), our approximation is controlled by the ratio r/r V and not α or Λ themselves.
We organize our approach in three steps; first we look for static configurations in the presence of the constant term φ c in the scalar field solution. Then we study the dynamical part of the scalar field, Φ 1 , and finally we check the back-reaction of Φ 1 onto metric. To characterise the static solutions, we use two quantities that are standard in the study of A necessary condition for a star to be stable against radial perturbations is that its mass must increase as the central density increases, dM/dρ c > 0. This is known as the static stability criterion [42]. The stability changes at the turning point of the mass-radius curve, as can be seen by comparing figure 2 and figure 3; from low to high density, all the configurations before the turning point are stable (i.e. the fractional binding energy becomes larger as the density increases), while the ones after the turning point are unstable. The maximum mass corresponds to the turning point of the binding energy curve. From the same figures we also see that a large φ c would have a dramatic effect on the observable characteristic of a neutron star. This means that the leading order scalar field coupling βφ c must be small. In the next section we study the behaviour of the first order correction to the scalar field solution Φ 1 .

The dynamics of the scalar field inside the star
Let us focus on the configurations before the turning point of the binding energy curve (i.e. those with mass lower than the maximum mass) so that we have stable solutions. In figure 4 we show the solution for the scalar field for these stable configurations. As we mentioned before, φ c should be seen as an initial condition which can be different for every density, thus in principle, it is possible to engineer a set of solutions with small deviations from GR for low densities (i.e. a small φ c ) and noticeable deviations for high densities. However, as we show  below, under some physical requirements such a behaviour cannot be achieved. We assume the following: • There exists an asymptotically flat solution outside the Vainshtein radius. This is supported by the results of [44]. In such a region, typically on cosmological scales, a bound on the asymptotic value of the scalar field, φ 0 , can be imposed. Note that βφ 0 measure the coupling strength between the scalar field and matter in this region.
• The solution for the lowest density considered in our analysis saturates that bound.
We ask how large are the deviations from GR, i.e. what is the value of φ c , when we require to keep the same φ 0 as we move from lower to higher density solutions. This is basically the same set up that leads to the discovery of strong gravity effects in the Brans-Dicke gravity [35]. The difference is that in order to evaluate the asymptotic scalar field φ 0 we need to consider the intermediate Vainshtein regime between the star and the asymptotically flat region. This also implies that the constraint on φ 0 does not come from the Solar System test in our case unlike in Brans-Dicke gravity. In order to compute φ 0 we need explicit expressions for the matchings of the scalar field solution both at R and r V . Beyond the Vainshtein radius, the -22 -

JCAP10(2014)055
usual scalar tensor theory solution, without Galileon, can be used, hence for r > r V , we assume the general static, spherically symmetric vacuum solution is ds 2 = −e νext dt 2 + e −νext dr 2 + (r 2 − ar)(dθ 2 + sin 2 θdϕ 2 ) , (5.10) where a, b and d are integration constants subjected to the condition a 2 −b 2 = 4d 2 . Note that this solution is in a gauge grr = g −1 tt , whereas the solutions for the star and the Vainshtein region are in the gauge g θθ = r 2 . To match quantities computed in those different gauges we need the relation between r andr, given by For the scalar field in the Vainshtein region, we use equation (2.15), with the positive sign to ensure that it is compatible with an asymptotically flat solution. The procedure is summarised as follows: we solve numerically for the scalar field inside the star, and then we match the solution at r = R with the scalar field solution in the Vainshtein regime. We compute the Vainshtein radius using (2.23) and by mathcing the solution at r = r V to the exterior solution, we obtain φ 0 . After doing this for the lowest density solution, we require that all the solutions for higher densities have the same φ 0 . This is achieved by changing the initial condition φ c at every different density. With this procedure, we obtain where Φ I 1 | R is evaluated using the numerical solution for the interior of the star, and all the other Φ refer to the scalar field solution in the Vainshtein region.
The results are shown in figure 5, where large values of φ c and Λ were chosen so that the small effects for higher density solutions are visible. A detailed understanding of these results can be achieved with figure 6. Since Φ in the Vainshtein region grows as √ r in weak gravity limit, then the dominant term in φ 0 is Φ| r In figure 6, we see that at low densities r V grows monotonically with both, the density and φ c , so in a rough approximation φ 0 ∝ M −4 p ρφ c . For this reason, if we want to keep φ 0 constant as we increase the density, φ c has to decrease, pushing the solutions towards the GR limit (remember φ c measures the GR deviations inside the star). At higher densities, the dependence of r V on ρ changes drastically due to strong gravity effects, and the Vainshtein radius begins to decrease as the density grows. However this is not enough to trigger large deviations from GR in our results.
At this point we have learned that in the cubic Galileon model a neutron star can carry a scalar field and still be indistinguishable from a neutron star in GR. The corrections to the metric of order α −1 are expected to be small too. In the next section we will confirm that this is the case. how large can the backreaction from the scalar field to the metric be. For this, we assume that the product φ c β is small, such that at the lowest order in α −1 , the solutions are close to GR, and we ask whether the first order corrections can cause large deviations or not. We assume that the large α expansion holds for the scalar field, i.e. we must have α −1 φ 1 < φ c everywhere inside the star. The question we ask is: how large are the corrections to the metric when the large α approximation for the scalar field approaches its limit of validity? There are two ways to approach such a limit; one is to make Λ larger so that the profiles for Φ 1 are amplified, and the other is to take an initial condition Φ 1 (r min ) close to φ c . We explore both possibilities to answer the previous question.

Corrections to the metric
In figure 7, we show Φ 1 and the respective corrections to g rr for Λ 3 /α 2 ∼ (10 −11 Km) −3 . The suppression of the large α 2 is not valid any more, as can be seen from the fact that  Figure 7. The left panel shows some scalar field profiles in the extreme case when Λ is increased up to the limit when Φ 1 (R) reaches values comparable to φ c , which is fixed to be φ c = 0.0005. Even in this case, the corrections to g rr are still suppressed, as can be seen in the right panel. Φ 1 becomes comparable to φ c near the radius of the star. Despite this, the effects on the metric, characterized by λ 1 in (5.5) are minimal. We do not need to go further to conclude that this corrections do not alter the results displayed above for the binding energy and the mass-radius curves.
As we mentioned before, the other possibility to approach the limit where the largeα approximation stops being reliable, is to increase the value of the initial conditions for Φ 1 , taking it closer to φ c . The results are shown in figure 8 for Φ 1 (r min ) = 0.0001, while φ c = 0.0005. The deviations from GR are again negligible, although they are larger than the deviations caused by a large Λ. This can be understood as a result of Φ 1 being of order φ c all the way across the star, while for the Λ-induced corrections, Φ 1 becomes large only near the surface of the star.
Until now, everything seems to indicate that a gravitational configuration with a Galileon scalar field is hardly distinguishable from a pure GR solution. The O(α 0 ) term of the scalar field can be arbitrarily small for low densities and, under a reasonable physical requirement, we showed that it tends to decrease even more for higher densities. The Vainshtein radius derived from the first order corrections to the scalar field is in a good agreement with the result in the weak gravity limit, and the effects of the scalar field Φ 1 on the metric are highly suppressed. However, there is still an issue with the solutions for the scalar field beyond the critical density. We comment on this in the next subsection.

Critical density for the scalar field
For general static solutions in vacuum, one of the beauties of the equation of motion for the Galileon scalar field is that it can be integrated once and it can be written as a simple quadratic algebraic equation for Φ , (2.5). In the presence of matter, this cannot be done as the density and pressure are unknown functions of r. However, it turns out that the quadratic nature of the scalar field equation eventually shows up. The value of the density for which the scalar field solution ceases to exist is not the density corresponding to the maximum mass for neutron stars, but the density at which the trace of the energy momentum tensor changes sign, causing the scalar field to become complex. A reason for this behaviour can be inferred from the the scalar field solution (4.10), whereβ plays the same role as the trace of the energy-momentum tensor, T . The parameterβ appears inside the square root, and if it changes sign at some radius, the solution ceases to exist. This is what happens inside neutron stars when T changes sign. Figure 9 shows the dependence of T on the density at r min . The existence of a critical density, ρ crit , where T changes sign is generic for relativistic matter, as we will explain further below.
A naive way around would be to change the sign of the scalar field coupling constant in the hope of finding a branch of solutions for ρ > ρ crit . This certainly works in the flat space toy model, as can be seen from equation (4.10). However, as shown in figure 10, for the configurations beyond ρ crit , the trace of the energy momentum tensor shows a further change of sign inside the star (this is independent of the sign of β), thus a continuous solution for the scalar field does not seem to exist.
The exact value of ρ crit depends upon the particular polytropic exponent and polytropic constant for the stellar model, and it can be larger or smaller than the density corresponding to the maximum mass. The maximum mass of a neutron star is still uncertain both from a theoretical and an observational point of view (see e.g., [34] for a review). However, it is known for sure that it cannot be less than 1.97 ± 0.04M , which corresponds to the mass of the pulsar J1614 − 2230 [38]. The particular polytropic parameters that we are using are adjusted in such a way that the predicted maximum mass is in agreement with this last value. Coincidentally, the onset of radial instabilities, signalled by the maximum mass and the change of sign for T , happen roughly at the same density, as can be seen in the left panel of figure 11. One may ask what happens if the maximum mass of neutron stars turns out to be larger than the current observational value? Comparing the left and right panel of figure 11, we learn that as the maximum mass increases by changing the polytropic parameters, the critical density for the existence of the scalar field decreases. This allows us to state that for a neutron star described by a polytropic equation of state, which is in agreement with the (current) observational maximum mass, the configurations that would collapse and form a black hole cannot support a non-trivial scalar field.
The (non-)existence of high density stars in modified gravity models has been discussed in the literature extensively. In the context of massive gravity, it was first claimed that there is no non-singular solution that connects the asymptotic solutions to the solutions inside a source [48]. Later ref. [49,50] found that this was a numerical artefact and they constructed numerically and analytically a solution that connects smoothly the solution inside a star to the asymptotic flat spacetime. However, it was found that, when the density of the object increases, the numerics becomes unstable and singularities are found to appear. It is still not clear whether those singularities are physical or just numerical artefacts. The existence of neutron stars was also debated in f (R) gravity. It was claimed that there is an upper bound for the mass of neutron stars as the solution could not be found for high density stars [51,52].  However, it was shown later that this was also due to numerical artefacts [53][54][55]. At the same time, it was found that static solutions can be found only when P < ρ/3 [53,54]. This is closely related to our finding although the reason for the non-existence of the static solutions is quite different. In f (R) gravity case, there is a potential for the scalar field. The coupling to the trace of energy-momentum tensor leads to the density (and pressure) dependent effective potential. When P < ρ/3, the effective potential does not have a minimum and stable static solutions cannot exist. In our case, there is no potential. Instead, the non-liner derivative coupling is responsible for suppressing the scalar field. We argued that the sign change of −T = ρ − 3P makes the solution complex and the solution ceases to exist. We note that a similar problem was found even in weak gravity when the density becomes negative in voids [33].

Conclusions
In this paper, we developed an expansion scheme to solve the non-linear equations inside the Vainshtein radius in a system with the cubic Galileon term. If applied to a test scalar field in the Minkowski spacetime, this expansion is controlled by (r/r V ) 3/2 , which is the same expansion obtained in a dual theory introduced in [40,41]. We applied this expansion method and found the first order corrections to the metric and the scalar field due to the Galileon term for spherically symmetric configurations. We also included a slow rotation in the system and derived again the first order correction. The integration constants in the vacuum solutions can be determined by introducing a star. In order to couple matter, we choose two different couplings, one linear and other queadratic in the scalar field. We start from the linearly coupled model and derived corrections to a low-density TOV solution. This solution could be connected smoothly to the exterior Schwarzschild solution with small corrections due to the Galileon term. We confirm that all deviations from the GR solutions are strongly suppressed within the Vanshtein radius, including corrections to the rotation. In the quadratically coupled models, in which "scalarisation" can strongly enhance deviations from GR in standard scalar tensor gravity, we showed that the Vainshtein mechanism suppressed the difference between the scalar field at the centre of the star and its asymptotic value, hence the scalarisation effect does not take place. Finally, we found, numerically, the corrections to relativistic star solutions due to the Galileon term. As expected, the higher the density, deviations from GR in the star are more suppressed for a given asymptotic value of the scalar field. Although in the strong gravity regime, the Vainshtein radius does not grow as the density increases, unlike in the weak gravity regime, this is not enough to have detectable deviations from GR. One caveat of our approach is that our expansion scheme is valid only inside the Vainshtein radius while the asymptotic solution is valid outside the Vainshtein radius. Thus strictly speaking, our matching at the Vainshtein radius is not well justified. Using the exact solution for a scalar field in the Minkowski spacetime, we checked that the inaccuracies caused by this procedure does not change orders of magnitude. Still we cannot guarantee that the matching always exists for the full solution. Full numerical studies of the solution are left for future investigations. Our approximated solutions will be useful to find exact solutions as it is a challenging task to construct a relativistic star numerically due to the hierarchy of scales in the system (i.e. the size of the star and the Vainshtein radius).
It is well know that these Vainshtein solutions do present, generically, sub-and superluminal propagation [6]. Our solution for the scalar field becomes the same as the flat spacetime solution at large radii r r s . Thus we expect that we cannot avoid superluminality. The situation may differ near the Schwarzschild radius due to the effects of the background curvature. This requires us to study time dependent fluctuations around our static solutions. This is an interesting problem and we leave this for future work. Finally, we found that the existence of the scalar field solution is not always guaranteed in the strong gravity regime. When the trace of the energy momentum tensor −T = ρ − 3P changes sign, the scalar field solution ceases to exist. We argued that this was due to the quadratic nature of the scalar field equation and the solution became complex at this point. For the equation of state that we used in this paper, this happens close to the maximum mass that corresponds to the onset of the gravitational instability. This equation of state is adjusted to explain the observed maximum mass of neutron stars. This indicates that unstable neutron stars cannot carry a scalar charge although we cannot exclude a possibility of dynamical solutions or non-spherically symmetric solutions with a non-trivial scalar configuration. By changing the equation state, the maximum mass does not coincide with the sign change of the trace of the energy-momentum tensor. However, we found that if the maximum mass increases, the critical density, above which the scalar field solution ceases to exist, decreases. This implies that, for a neutron star described by a polytropic equation of state in agreement with the observational maximum mass, a cubic Galileon scalar field cannot exist for the configurations that would collapse and form a black hole. Along the same line, there is a no-hair theorem for black holes in the presence of the Galileon scalar field [43]. This theorem proves that static, spherically symmetric black hole solutions cannot sustain non-trivial scalar profiles. It is still an open question what is the fate of the scalar field when a star's density increases dynamically and eventually collapses to become a black hole.