Solar System constraints on massless scalar-tensor gravity with positive coupling constant upon cosmological evolution of the scalar field David Anderson and Nicolás Yunes eXtreme Gravity Institute, Department of Physics, Montana State University, Bozeman, Montana 59717, USA (Received 23 May 2017; published 22 September 2017) Scalar-tensor theories of gravity modify general relativity by introducing a scalar field that couples nonminimally to the metric tensor, while satisfying the weak-equivalence principle. These theories are interesting because they have the potential to simultaneously suppress modifications to Einstein’s theory on Solar System scales, while introducing large deviations in the strong field of neutron stars. Scalar-tensor theories can be classified through the choice of conformal factor, a scalar that regulates the coupling between matter and the metric in the Einstein frame. The class defined by a Gaussian conformal factor with a negative exponent has been studied the most because it leads to spontaneous scalarization (i.e. the sudden activation of the scalar field in neutron stars), which consequently leads to large deviations from general relativity in the strong field. This class, however, has recently been shown to be in conflict with Solar System observations when accounting for the cosmological evolution of the scalar field. We here study whether this remains the case when the exponent of the conformal factor is positive, as well as in another class of theories defined by a hyperbolic conformal factor. We find that in both of these scalar-tensor theories, Solar System tests are passed only in a very small subset of coupling parameter space, for a large set of initial conditions compatible with big bang nucleosynthesis. However, while we find that it is possible for neutron stars to scalarize, one must carefully select the coupling parameter to do so, and even then, the scalar charge is typically 2 orders of magnitude smaller than in the negative-exponent case. Our study suggests that future work on scalar-tensor gravity, for example in the context of tests of general relativity with gravitational waves from neutron star binaries, should be carried out within the positive coupling parameter class. DOI: 10.1103/PhysRevD.96.064037 I. INTRODUCTION General relativity (GR) has been shown to be consistent with all current observations, including those in the Solar System (SS) [1,2], with binary pulsars [3], and recently with gravitational waves through the merging of two black holes detected by advanced LIGO [4,5]. Scalar-tensor theories (STTs) of gravity are among the most natural extensions to GR [6–8] because they are both well posed [9] and well motivated from string theory [10,11], quantum field theory [12], and cosmology [13–16]. These theories have seen a recent revival given that they can lead to large deviations from GR in strong-field scenarios, making them a natural theory choice to test GR with gravitational-wave observations. STTs modify GR by introducing a scalar field that couples to the metric nonminimally, thus forcing matter to respond both to the metric tensor and to the scalar field. The easiest way to see this is to (conformally) transform the metric tensor into the Einstein frame, in which the action looks identical to the Einstein-Hilbert one, but with the matter sector responding to the product of the conformally transformed metric and the conformal factor. Therefore, STTs can be classified by the specific functional form of the conformal factor. One of the most studied classes is that originally proposed by Damour and Esposito-Farèse (DEF) [17,18], in which the conformal factor is a Gaussian in the Einstein-frame scalar field with an exponent proportional to a coupling constant β, i.e. AðφÞ ¼ expðβφ2=2Þ. Another class that has recently captured some attention is that in which the conformal factor is proportional to a certain power of a hyperbolic cosine with argument proportional to the product of the Einstein-frame scalar field and a constant β, i.e. AðφÞ3β ¼ coshð ffiffiffi 3 p βφÞ [19–22]. In this paper, we will study both classes, referring to the former as expo- nential and the latter as hyperbolic. Perhaps one of the most interesting consequences of STTs is that neutron stars (NSs) can undergo what is known as scalarization [17,18], in which the scalar field is amplified above its background value inside the NS. Three types of scalarization processes have been discov- ered, which can be classified by the physical process that causes the activation of the scalar field. For isolated NSs, spontaneous scalarization occurs once a critical density is reached inside the star [17,18,23,24]. In binary systems, NSs can undergo either dynamical or induced scalarization [25–28]. The former occurs when the binding energy of the orbit reaches a critical value, causing the violent activation of the scalar field, while the latter occurs when one member PHYSICAL REVIEW D 96, 064037 (2017) 2470-0010=2017=96(6)=064037(20) 064037-1 © 2017 American Physical Society https://doi.org/10.1103/PhysRevD.96.064037 https://doi.org/10.1103/PhysRevD.96.064037 https://doi.org/10.1103/PhysRevD.96.064037 https://doi.org/10.1103/PhysRevD.96.064037 of the binary is already scalarized and it induces a scalar charge in its companion. Until recently, scalarization was thought to only occur in the exponential class of STTs if β < 0, because then a certain scalar instability can occur for most equations of state (EoSs). But recently, a nonvanish- ing scalar charge has also been shown to occur inside NSs when β > 0 in both exponential and hyperbolic STTs with certain EoSs [19,20,29]. In these cases, however, the activation of the scalar charge [for example, as a function of Arnowitt-Deser-Misner (ADM) mass] is not nearly as sudden as in the β < 0 case and the maximum scalar charge is typically significantly smaller. STTs have also been popular because of the belief that they can be made to pass Solar System constraints by appropriately choosing the asymptotic value of the scalar field, but recently this belief was proven to be incorrect [30–33]. Solar System observations place a very strict constraint on exponential STTs because the asymptotic value of the scalar field in the Solar System cannot be chosen freely, but rather it must be chosen consistently with the field’s prior cosmological evolution. The cosmological evolution of the Universe forces the scalar field to grow rapidly with redshift when β < 0 (for all but very fine-tuned initial conditions), leading to clear violations of Solar System tests today, most notably in the consistency of the Shapiro time delay and the perihelion shift of Mercury with the GR predictions [34]. A previous study [30] attempted to remedy this problem by modifying the form of the exponential conformal factor by adding a term of higher order in the Einstein frame scalar field. While this was enough to ensure Solar System tests are passed today, it reduced the scalar charge in NSs by orders of magnitude. Given all of this, it is then natural to explore exponential and hyperbolic STTs with β > 0 in more detail to deter- mine if Solar System constraints can be placed, and if so, the degree to which NSs continue to scalarize within the allowed region of parameter space. We first study Solar System constraints by exploring the cosmological behavior of STTs from the time of big bang nucleosynthesis (BBN) until today. For any particular value of β and any particular set of initial conditions consistent with BBN constraints, the value of the scalar field today will either be consistent with Solar System observations or it will violate them, therefore determining if that value of β is part of the viable region of parameter space for those initial conditions. There are two distinct sets of initial conditions that we consider in detail. (1) BBN constraints on the scalar-field initial conditions are saturated such that deviations from GR are at a maximum. Figure 1 shows the β (blue) regions that allow hyperbolic STTs to pass Solar System tests when choosing initial conditions that saturate BBN constraints. For all viable initial conditions (including those that saturate BBN con- straints), hyperbolic STTs pass Solar System tests when 0 < β ≲ 17, while exponential STTs pass when 0 < β ≲ 24 (not shown in the figure). Moreover, there are special values of the initial conditions for which both theories pass Solar System tests above these critical β values, which are shown as thin blue regions in Fig. 1. This is because the cosmo- logical evolution of the scalar field has damped oscillations with redshift, thus allowing for the possibility that today the field happens to be in a state in which Solar System tests are passed, although this may not be the case in the future. (2) BBN constraints on the scalar-field initial conditions are relaxed such that deviations from GR are not at a maximum. The constraints reported above can be relaxed by fine-tuning the initial conditions at the time of BBN. Cosmological observations and the theory of nucleosynthesis require that the scalar field at the time of BBN be in a small region near FIG. 1. Maximum value of scalar charge αsc inside stable NSs as a function of the coupling parameter β in hyperbolic STTs with five different EoSs. In the background, blue (red) regions correspond to values of β that are consistent (inconsistent) with Solar System observations after cosmological evolution. Observe that the blue regions become thinner and more sparse as β increases. The bounds shown here are the most stringent ones one can place with initial conditions that deviate from GR as much as possible by saturating BBN constraints. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-2 the minimum of an effective potential. Saturating current constraints on the speed-up factor, which quantify deviations of the gravitational constantG from GR at the time of BBN, put the scalar as far away from the minimum as possible. In principle, however, one can relax this assumption and place the scalar closer to the minimum, making it easier for the theory to satisfy current Solar System constraints. A random distribution of viable initial conditions relaxes the bounds reported above to 0 < β ≲ 34 and 0 < β ≲ 25 in the expo- nential and hyperbolic cases respectively. Larger values of β would typically require fine-tuning of initial conditions for Solar System tests to be passed. It is worth noting that if these assumptions are completely relaxed such that there is no deviation fromGR at the time of BBN, then one cannot place any constraints on the value of β because the cosmological evolution of the Universe will be completely consistent with GR. With this at hand, we then explore whether spontaneous scalarization still occurs in exponential and hyperbolic STTs within the viable β regimes of parameter space. In the exponential STT case, scalarized NSs have already been shown to be unstable to gravitational collapse for such values of β [19]. In the hyperbolic case, however, scalarized NSs can be stable for certain EoSs, but the scalar charge is typically 2 orders of magnitude smaller than typical charges in the β < 0 exponential STT case. Figure 1 plots the maximum scalar charge for stable NSs in hyperbolic STTs for a large range of β using different EoSs. None of the EoSs that we consider here give stable scalarized NS solutions for β ≲ 17, but they do allow for stable scalarized stars when β ≳ 20. Moreover, there exist (small) periodic regions of viable parameter space that allow for scalarized NSs (blue regions that have nonzero scalar charge in them), although these separate and become thinner as β increases. These general conclusions are summarized in Table I. The results in Table I only hold for massless STTs and the addition of a potential term for the scalar field could allow the theories to completely avoid any constraints we place here. Our results are directly relevant to studies that constrain deviations from GR in the strong field, for example with gravitational waves. Most of the latter have so far focused on the massless exponential class of STTs with β < 0, but these have already been shown to be incompatible with Solar System observations [30,31]. We will here show that themassless exponential class of STTswith β > 0 does pass Solar System tests for a range of β, but these theories have already been shown to disallow stable scalarized NSs [19]. We will also show that the hyperbolic class of STTs with β > 0 also passes Solar System tests for a range of β, and these theories do allow for stable scalarized stars [19] making themmuchmore appropriate for gravitational-wave studies. Our results, however, also indicate that the amount of scalar charge in such theories is drastically smaller than what is obtained in the massless exponential case with β < 0. This raises the question of whether such small charges can be constrained in practice with second-gener- ation ground-based gravitational-wave detectors. A more detailed analysis is required to address this last question. The rest of this paperwill address the points above inmore detail. We first lay out our notation and the background of STTs in Sec. II. We next go into details of the cosmological evolution of the scalar field in both theories in Sec. III and place constraints on β from Solar System observations. Section IVaddresses NS solutions in detail, first for β < 0 to build a foundation and then for β > 0. Finally, we end with our conclusions in Sec. Vand discuss the implications of our results. We adopt units in which c ¼ 1 but restore cgs units in our discussion of NSs when appropriate. II. THE BASICS OF SCALAR-TENSOR THEORIES In this sectionwe present the details of the class of theories we investigate and establish notation, mostly following the presentation in Ref. [30]. The STTswe consider in this paper can be defined by the action SJ ¼ SJ;g þ SJ;mat, where the gravitational part is given by SJ;g ¼ Z d4x ffiffiffiffiffiffi−gp 2κ � ϕR − ωðϕÞ ϕ ∂μϕ∂μϕ � ; ð1Þ and where g and R are the determinant and Ricci scalar associated with the Jordan-frame metric gμν, ωðϕÞ is a coupling function for the scalar field ϕ, and κ ¼ 8πG where G is the bare gravitational constant. The matter action SJ;mat½χ; gμν� is a functional of the matter fields χ, which couple directly to the Jordan frame metric gμν, meaning that STTs of this form are metric theories of gravity. One can take the action in Eq. (1) and write it in a form resembling the Einstein-Hilbert action in GR through a conformal transformation gμν ¼ A2ðφÞg�μν where we refer to g�μν as the Einstein-frame metric. By choosing the conformal factor AðφÞ such that A2ðφÞ ¼ ϕ−1; ð2Þ and providing an explicit relationship between φ and ϕ via TABLE I. A summary of the various properties of the ex- ponential (exp.) and hyperbolic (hyp.) classes of STTs. The columns correspond to the following: Do they pass Solar System tests after cosmological evolution? What is the typical magnitude of the maximum scalar charge they can support inside NSs? Do they allow for stable and scalarized NS solutions? Theory SS tests? αsc;max Stable NS? β < 0 exp. ✗ Oð10−1Þ ✓ β > 0 exp. ✓ Oð10−3Þ ✗ β < 0 hyp. ✗ Oð10−1Þ ✓ β > 0 hyp. ✓ Oð10−3Þ ✓ SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-3 α2ðφÞ ¼ � d lnAðφÞ dφ � 2 ¼ 1 3þ 2ωðϕÞ ; ð3Þ the Einstein-frame action becomes SE ¼ Z d4x ffiffiffiffiffiffiffiffi−g� p 2κ ½R�−2gμν� ∂μφ∂νφ�þSE;mat½χ;A2ðφÞg�μν�; ð4Þ where g� and R� are the determinant and Ricci scalar associated with the Einstein-frame metric g�μν. Notice that the matter degrees of freedom χ now couple to the product of the conformal factor AðφÞ and the Einstein-frame metric. For later convenience, we define a conformal coupling potential via α ¼ ∂Vα=∂φ, such that VαðφÞ ¼ lnAðφÞ; ð5Þ and note that the curvature of this potential is then given by βðφÞ ¼ ∂2Vα ∂φ2 ¼ ∂αðφÞ ∂φ : ð6Þ The potential in Eq. (5) will later play the role of an effective potential that the scalar field evolves in cosmo- logically. For the rest of this paper we will strictly refer to AðφÞ as the conformal factor, αðφÞ as the (conformal) coupling, and VαðφÞ as the (conformal) coupling potential. The variation of the Einstein-frame action yields the field equations R� μν ¼ κ � T� μν − 1 2 g�μνT� � ; ð7Þ □�φ ¼ − κ 2 αðφÞTmat;�; ð8Þ where □� is the Einstein-frame covariant wave operator and Tmat;� is the trace of the Einstein-frame matter stress- energy tensor defined by Tmat;� μν ≡ 2ffiffiffiffiffiffiffiffi−g� p δSE;mat½χ; A2ðφÞg�μν� δgμν� : ð9Þ One can see that the Einstein-frame stress-energy tensor is related to the Jordan-frame one via Tμν mat;� ¼ A6Tμν mat by applying the conformal transformation to Eq. (9). The field equations also depend on the total stress-energy tensor, which is the sum of the matter and scalar field stress-energy tensors T� μν ¼ Tmat;� μν þ Tφ;� μν ; ð10Þ where κTφ;� μν ¼ 2ð∂μφÞð∂νφÞ − g�μνð∂σφÞð∂σφÞ: ð11Þ Our entire study, both cosmologically and for NSs, adopts a perfect fluid representation of the matter in the Jordan frame. However, this is also true in the Einstein frame since they are related through a conformal trans- formation, and thus we use Tmat;� μν ¼ ðϵ� þ p�Þu�μu�ν þ p�g�μν; ð12Þ where ϵ� is the energy density of the fluid, p� is the pressure, and u� is the 4-velocity, all in the Einstein frame. The normalization of the 4-velocity g�μνu μ �uν� ¼ −1 ¼ gμνuμuν leads to the relationship uμ� ¼ Auμ. Using this result in conjunction with Tμν � ¼ A6Tμν allows one to derive the direct relations ϵ� ¼ A4ϵ and p� ¼ A4p between the Einstein and Jordan frame density and pressure. The choice of AðφÞ, or consequently αðφÞ, defines a particular theory and therefore plays a critical role in testing that theory with observations. Specifically, these functions determine the local value of Newton’s gravita- tional constant [7] GN ¼ G½A2ð1þ α2Þ�φ0 ; ð13Þ where φ0 is the value of the scalar field today determined cosmologically. They also determine the local value of the parametrized post-Newtonian (PPN) parameters [35,36]. The γPPN parameter, which measures the spatial curvature due to a unit rest mass [1,37], is given in massless STTs by j1 − γPPNjφ0 ¼ 2α20 1þ α20 ; ð14Þ where α0 is the conformal coupling evaluated at the present- day value of the scalar field φ0. There exists another such expression for the βPPN parameter which is given by j1 − βPPNjφ0 ¼ 1 2 β0α 2 0 ð1þ α20Þ2 ; ð15Þ where β0 ¼ ∂α=∂φjφ0 . In GR γPPN ¼ βPPN ¼ 1, and the constraints j1 − γPPNj≲ 2.3 × 10−5 and j1 − βPPNj≲ 8 × 10−5 [1], from the verification of the Shapiro time delay by the Cassini spacecraft and the perihelion shift ofMercury respectively [37], provide the weak-field constraints on the theories we consider in this paper. III. COSMOLOGICAL EVOLUTION AND SOLAR SYSTEM CONSTRAINTS To determine if a particular theory is consistent with Solar System tests, we must first understand the cosmo- logical evolution of the scalar field. We here review the cosmological evolution equations originally presented in DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-4 Refs. [32,33] and continue to establish notation, again mostly following the presentation in Ref. [30]. We adopt a spatially flat Friedmann-Robertson-Walker metric in the Einstein frame ds2� ¼ −dt2� þ a2�ðdr2� þ r2dΩ2�Þ; ð16Þ where a� is the Einstein-frame scale factor. The Einstein- frame field equations then become 3H2� ¼ κϵ� þ _φ2; ð17Þ −3 ä� a� ¼ κ 2 ϵ�ð1þ 3λÞ þ 2 _φ2; ð18Þ φ̈þ 3H� _φ ¼ − κ 2 αϵ�ð1 − 3λÞ; ð19Þ where overdots stand for derivatives with respect to the Einstein-frame coordinate time t�, H� ¼ _a�=a� is the Einstein-frame Hubble parameter, ϵ� and p� are the energy density and pressure of all components of the Universe (matter, radiation, dark energy), and λ ¼ p�=ϵ� is the usual cosmological EoS parameter. The Einstein-frame variables can be transformed to Jordan-frame variables via the conformal transformation mentioned earlier, gμν ¼ A2ðφÞg�μν which yields dt ¼ AðφÞdt� and a ¼ AðφÞa�. Following Ref. [30], we assume that any pressure and energy density associated with φ will be negligible com- pared to other energy components of the Universe. Equations (17)–(19) do not lend themselves to a simple analytic solution for the evolution of the scalar field φ. However, by defining a new time coordinate dτ ¼ H�dt�, one can decouple the evolution equations. Under this time transformation one finds that the time derivatives become _φ ¼ H�φ0; ð20Þ φ̈ ¼ H2�φ00 þ _H�φ0; ð21Þ ä a ¼ _H� þH2�; ð22Þ where primes denote derivatives with respect to τ. Using these redefinitions, Eqs. (17)–(19) can be simplified into a single equation for the evolution of the scalar field 2 3 − φ02 φ 00ϵ� þ ϵ�ð1 − λÞφ0 ¼ −αðφÞϵ�ð1 − 3λÞ: ð23Þ It is important to note here that if one wishes to consider more than a single-component universe, ϵ� must not be divided out of the equation. In general, each term contain- ing ϵ� and λ will be a sum over all components of the Universe, and therefore, ϵ� cannot simply be removed from the equation. To have a complete description of Eq. (23) one must understand its evolution in τ. From conservation of energy in the Jordan frame one finds dðϵa3Þ ¼ −pðϵÞdða3Þ; ð24Þ which, using λ ¼ p=ϵ, can be written as ϵ ¼ a−3ð1þλÞϵ0; ð25Þ where ϵ0 is the value of the Jordan-frame energy density measured today (i.e. when a ¼ 1). To transform Eq. (25) to the Einstein frame we make use of the fact that dτ ¼ H�dt� to write a� ¼ eτ and then use a ¼ AðφÞa� to write a ¼ AðφÞeτ: ð26Þ We can then use ϵ� ¼ A4ϵ in combination with Eq. (26) to write the Einstein-frame energy density as ϵ� ¼ ϵ0e−3ð1þλÞτA1−3λ: ð27Þ If we divide Eq. (23) by the critical density ϵcrit and use the result of Eq. (27) the scalar field evolution equation becomes 2φ00 3 − φ02 C1ðφ; τÞ þ φ0C2ðφ; τÞ ¼ −αðφÞC3ðφ; τÞ; ð28Þ where we have defined C1ðφ; τÞ ¼ X i Ωi;0e−3ð1þλiÞτAð1−3λiÞ; ð29Þ C2ðφ; τÞ ¼ X i Ωi;0ð1 − λiÞe−3ð1þλiÞτAð1−3λiÞ; ð30Þ C3ðφ; τÞ ¼ X i Ωi;0ð1 − 3λiÞe−3ð1þλiÞτAð1−3λiÞ; ð31Þ where Ωi;0 ¼ ϵi;0=ϵcrit is the standard density parameter as seen today for the ith energy component. The evolution described in Eq. (28) is reminiscent of a damped oscillator. Assuming the density functions C1, C2, and C3 only contain a single energy component for simplicity, it becomes apparent that there is a velocity- dependent mass appearing in front of the second derivative—a damping term proportional to the velocity— and on the right-hand side of the equations there is a potential gradient influencing the overall evolution of the scalar field, as first pointed out in Refs. [32,33]. Even if the density functions are not constant in τ or φ, the damped oscillator picture still holds, as we will show numerically later in this section, and it will thus provide a useful analogy to extract a physical understanding from our numerical calculations. SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-5 A. Two scalar-tensor models Thus far we have remained model independent and have not made any assumptions about the form of αðφÞ. In this section we define the two STTs that we study in this paper. The first model we consider is the standard DEF theory, defined by AðφÞ ¼ e 1 2 βφ2 ; ð32Þ VαðφÞ ¼ 1 2 βφ2; ð33Þ αðφÞ ¼ βφ; ð34Þ βðφÞ ¼ β; ð35Þ which we will refer to as the exponential model. The other model we consider is that introduced by Mendes [19,20], defined by AðφÞ ¼ ½coshð ffiffiffi 3 p βφÞ� 13β; ð36Þ VαðφÞ ¼ 1 3β ln ½coshð ffiffiffi 3 p βφÞ�; ð37Þ αðφÞ ¼ tanhð ffiffiffi 3 p βφÞffiffiffi 3 p ; ð38Þ βðφÞ ¼ βsech2ð ffiffiffi 3 p βφÞ; ð39Þ which we will refer to as the hyperbolic model. Notice that β enters both models as a free parameter that will quantify the degree of departure from GR (with β ¼ 0 reducing the theory to GR). The conformal factor here is motivated as an analytic approximation to the conformal factor one would find if the action also included terms proportional to Rφ2 in the Einstein frame [19–22], motivated by quantum field theory considerations [12]. In the context of an oscillator, Vα is the potential that sources the scalar field evolution through its gradient αðφÞ. In the case of exponential STTs, this is simply a parabola, cf. the left panel of Fig. 2, and therefore the scalar field evolves analogously to a damped harmonic oscillator in τ. In the case of hyperbolic coupling, however, that potential is no longer a simple parabola and while we expect the evolution of the scalar field to still be oscillatory, it should be quantitatively different than in the exponential STT case. The key difference between the potentials in these theories lies in their behavior away from their minima (φ ¼ 0). In the exponential case, the gradient of the potential grows linearly as we increase φ (or β for that matter) where for hyperbolic coupling the gradient becomes constant; both behaviors can be seen in the right panel of Fig. 2. The behavior of these potentials away from their minima plays a key role in the behavior of the scalar field, both in cosmology and in NSs, as we will show in the following sections. B. BBN constraints and initial conditions To determine the evolution of the scalar field one must first determine a set of initial conditions that it must obey at some point in the past. The natural first choice is to determine initial conditions at the very beginning of the Universe, or at the very least at the end of inflation. However, there is substantial uncertainty in the true physics of the early Universe, and thus, any constraints would be subject to this uncertainty. Following the work in Ref. [30], we will use the observational constraints from BBN to determine the behavior of the scalar field at matter-radiation equality (Z ∼ 3600), which we will use as our “initial” point in time for our numerical integrations. FIG. 2. Schematic diagrams for both the conformal coupling potential Vα (left) and the conformal coupling αðφÞ (right). The exponential STT is represented by solid red curves and the hyperbolic theory by dashed blue curves. The horizontal gray lines at�1= ffiffiffi 3 p in the right panel correspond to the limiting values of αðφÞ in the hyperbolic theory. These figures were constructed with β ¼ 5, but the qualitative features remain the same for all positive values of β. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-6 Consider an era in which the Universe is radiation dominated, as was the case in the early Universe after inflation, when a ∼ 10−5. In this case we can neglect the energy density contributions from matter and dark energy and write Eqs. (28)–(31) as simply 2 3 − φ02 φ 00 þ 2 3 φ0 ¼ 0; ð40Þ where we have divided out the common energy density terms. Notice that the coupling function αðφÞ does not appear in this equation, meaning that the scalar field evolution during a purely radiation-dominated universe will always be the same, regardless of the coupling function of the theory. Equation (40) has a solution of the form (following the notation in Refs. [32,33]) φðτÞ ¼ φr − ffiffiffi 3 p ln ½Ke−τ þ ð1þ K2e−2τÞ1=2�; ð41Þ where K ¼ φi 0 ffiffiffi 3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − φi 02=3 p ; ð42Þ where φr is a constant and φi 0 is the scalar velocity upon exiting inflation. The only necessary assumption that we must make here is that the scalar velocity is not very close to its limiting value of ffiffiffi 3 p . This assumption is valid because if the scalar velocity ever does reach ffiffiffi 3 p then all dynamics of the evolution are removed, which can be seen from Eq. (28). When φ0 ¼ ffiffiffi 3 p the only solution to Eq. (28) is one in which the field grows linearly in time, which certainly violates Solar System constraints today since the right-hand side of Eq. (14) would asymptote to two. Thus, we can conclude that φ0 < ffiffiffi 3 p for all times in the past and from Eq. (41) we can conclude that φ0ðτÞ is exponentially damped to zero during the radiation-domination era. To constrain the initial value of the scalar field we use BBN constraints on the gravitational constant. The speed- up factor ξBBN ≔ H=HGR quantifies any differences between the actual Hubble parameter H and the Hubble parameter predicted by GR, HGR; these differences are caused by deviations from the standard gravitational con- stant. The observed Hubble parameter can be derived from Eq. (17) with _φ ¼ 0 (as we have just argued to be true immediately after the radiation era), and is given by H ¼ ð8π=3ÞGA2 RϵR where ϵR and AR are the Jordan-frame energy density and conformal factor, respectively, evalu- ated at the time of BBN. The expansion rate predicted in GR takes the form H2 GR ¼ ð8π=3ÞGNϵR, with GN given in Eq. (13) with α0 ¼ 0. The speed-up factor then becomes ξBBN ¼ � H HGR � ¼ � GA2 R GN � 1=2 ¼ 1ffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ α20 p AR A0 ; ð43Þ where the R and 0 subscripts indicate values at the end of the radiation era (or at the time of BBN more specifically) and today respectively. Observational bounds on the abundance of helium place constraints on the expansion rate of the Universe [38–40], which then implies that j1 − ξBBNj ≤ 1 8 ; ð44Þ which is a one-σ constraint. On the other hand, Solar System tests restrict α20 ∼ 10−5 which means that to a good approximation we can neglect its contributions in Eq. (43) and write ξBBN ∼ AR=A0. A constraint on AR can be placed through Eq. (44), giving ����1 − AR A0 ���� ≤ 1 8 : ð45Þ Saturating this constraint gives AR=A0¼7=8 or AR=A0¼9=8. From here the particle analogy can be used to understand what exactly this constraint tells us. Relating the constraint on the conformal factor A to the conformal coupling potential Vα via Eq. (5) one finds Vα;0 þ lnð7=8Þ ≤ Vα;R ≤ Vα;0 þ lnð9=8Þ: ð46Þ This tells us that the scalar field can, at most, be lnð9=8Þ higher in the potential than where it is today or lnð7=8Þ lower. However, both conformal potentials we consider here, Eqs. (33) and (37), have a global minimum which the scalar field must be close to today in order to satisfy the Cassini bound. Such a condition means that the scalar field cannot be lower in the potential than it is today and therefore only Vα;R ≤ lnð9=8Þ is a valid constraint (because Vα;0 ≈ 0Þ. For our calculations we will use Vα;R ¼ ζ lnð9=8Þ; ð47Þ where ζ is a parameter that ranges from 0 to 1 and will scale the BBN constraint to allow us to sample the full range of possible initial conditions. Now we must apply the BBN constraint in Eq. (47) to the conformal coupling potentials we consider in this paper. By equating Eq. (33) and Eq. (47) we find the initial conditions for the exponential STT to be φR;exp ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 β ζ ln � 9 8 �s ; ð48Þ φ0 R;exp ¼ 0: ð49Þ Likewise, in the hyperbolic case we have SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-7 φR;hyp ¼ 1ffiffiffi 3 p β cosh−1 �� 9 8 � 3βζ � ; ð50Þ φ0 R;hyp ¼ 0: ð51Þ The constraints on the conformal factor from BBN also provide information on how far in the past BBN occurred in τ time. Consider the definition dτ ¼ H�dt� and integrate it to get τ ¼ ln a�. A transformation back to the Jordan-frame scale factor via a ¼ AðφÞa� gives τ ¼ ln � 1 1þ Z � − lnðAÞ; ð52Þ where Z is the Jordan-frame redshift related to the Jordan- frame scale factor through a ¼ 1=ð1þ ZÞ. Because we set a ¼ 1 today, the value of τ today must be zero, i.e. τ0 ¼ 0. Finding the difference in τ time between now and when BBN occurred is simply τ0 − τR, or τR ¼ ln � 1 1þ ZR � − ζ ln � AR A0 � ; ð53Þ where we have made use of the fact that the redshift today is zero. Notice that the ratio AR=A0 is the same ratio that is constrained from BBN, and thus AR=A0 ≤ 9=8, which is valid for all values of β. Therefore, using ZR ¼ 3600 and ζ ¼ 1, one finds that BBN occurred at most τ ¼ τR ∼ −8.306. The value of τR determined in Eq. (53) with ZR ¼ 3600 and a given ζ will mark the initial time of our numerical calculations. C. Exponential theory We study the cosmological evolution of the scalar field in detail and place bounds on β. We start by considering the ζ ¼ 1 case, i.e. saturating the BBN constraints on the speed-up factor, to understand the nature of the evolution for 1 < β < 1000. From there, we relax this assumption on ζ and allow 0 < ζ < 1, and we use these results to place bounds on β by avoiding fine-tuning scenarios. In order to avoid confusion, we point out that, in general, Solar System constraints on γPPN and βPPN place relatively independent constraints on the value of β. For completeness, we present both results but in the end our conclusions will be based off of the PPN parameter that places the tightest constraints on β; in every case that we study, we find that βPPN con- sistently places the tightest constraints.1 1. Scalar-field evolution We solve Eq. (28) numerically from the end of the radiation-dominated era (Z ∼ 3600) to the present day. Describing the energy density of the Universe as a piece- wise function in τ time would only require one to consider a single energy component of the Universe at a time and thus only keep a single term in Eqs. (29)–(31) (this leads to a rather simple analytic solution for the scalar field [30–33]). Treating the energy densities in this manner, however, forces sharp transitions between different cosmological eras (i.e. from matter domination to dark energy domina- tion). To avoid any sharp transition and ensure a smooth evolution, we simply consider all energy components of the Universe for all time and evolve the full form of Eqs. (28)–(31) numerically. The left panel of Fig. 3 shows the evolution of 1 − γPPN in the exponential theory for the entire integration time and three different values of β > 0. A general feature of the evolution regardless of β is that there are oscillations inside FIG. 3. Left: 1 − γPPN evaluated during the evolution of the scalar field, with ζ ¼ 1 here, to show that an increase in β makes it harder to satisfy the Cassini bound today at τ0 ¼ 0. The dips here physically go to zero but due to finite numerical precision this is not reflected here. Right: 1 − γPPN (red) and 1 − βPPN (blue) evaluated today at τ0 ¼ 0 for a range of β. 1Technically, constraints on _GN=GN [as in Eq. (13)] can also constrain these theories because the oscillating scalar field induces an oscillating value of GN . However, we find that these constraints are insignificant when compared to γPPN and βPPN in the range of parameter space that we explore. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-8 of a decaying envelope, which should not come as a surprise. The scalar field starts at some height in the potential and oscillates about the minimum. Because there is a damping term in the evolution equation, the amplitude of the oscillations must decrease in time. The dips in this plot, which go to zero, are points when the scalar field passes through the minimum of the potential, i.e. when αðφÞ ¼ 0, and the peaks correspond to turning points in the evolution where φ0 instantaneously vanishes. There are a few features of note that are present for different values of β. As β increases, the period of oscillations about the minimum becomes smaller and therefore the motion is more rapid. This is primarily due to the steepness of the corresponding potential, i.e. when β is large the potential is steep and narrow, allowing the velocity to be large. This accounts for two features of the evolution: 1) the period of oscillation is shorter and 2) the damping is weaker so the decay is slower in τ time. The first point above is obvious: starting from a higher point in the potential means there is more “potential energy,” which converts to “kinetic energy,” leading to a faster velocity and therefore faster oscillations, but the second point may seem counterintuitive. Recalling Eq. (28), as φ0 becomes larger the term 3 − φ02 becomes smaller, meaning that the damp- ing term proportional to φ0 alone holds less significance during the evolution. This velocity-dependent mass term is what is responsible for the decreased damping at the beginning of the evolution as β increases. A final feature of the evolution worth noting is the different decay rates that occur in the left panel of Fig. 3 (most apparent for β ¼ 100). One would expect to see different decay rates for different intervals of τ during the evolution because damping is determined by the value of C2 in Eq. (30). As time passes, the dominant energy density of the Universe also changes, which causes C2 to be dominated by that energy and results in a distinct decay rate. Before τ ¼ −6 the effects of both radiation and matter are important, until about τ ¼ −1 matter seems to com- pletely dominate the damping, and for τ ≳ −0.5 the effect of dark energy becomes dominant causing the fastest decay. Although we do not show it explicitly, similar conclusions can be drawn from the behavior of βPPN throughout the evolution of the scalar field. The right panel of Fig. 3 shows the value of j1 − γPPNj (red) and j1 − βPPNj (blue) today for a subset of β that we considered. As one may expect from the previous results, as β becomes larger the damping becomes weaker and as a result the field may not settle to the minimum of the potential by today. There are, however, periodic values of β that do allow for Solar System tests to be passed even when β is large. These values of β lead to an evolution where the scalar field just happens to be passing by the minimum of the potential today and thus the Solar System bounds are satisfied. However, as time continues to progress the scalar field will leave the minimum of the potential in the future and may eventually fail to satisfy Solar System constraints. Such an outcome becomes unavoidable for β ≳ 30. The apparent linear β dependence of ð1 − γPPNÞφ0 can be explained by an approximate analytic solution to Eq. (28). For simplicity, consider an era that is dominated by matter and assume φ0 ≪ ffiffiffi 3 p such that the scalar-field evolution can be approximated by 2 3 φ00 þ φ0 ¼ −βφ; ð54Þ which is a classical damped harmonic oscillator. With the initial conditions in Eqs. (48)–(49), with ζ ¼ 1, the solution becomes φðτÞ ¼ φRe−TðτÞ � cosðBTðτÞÞ þ sinðBTðτÞÞ B � ð55Þ where TðτÞ ¼ 3ðτ þ τRÞ=4, B ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffið8=3Þβ − 1 p , and τR ¼ −8.306. Note that this solution is valid only for the underdamped case such that β > 3=8, which is precisely the region of parameter space we are interested in. For a more detailed discussion of the critically and overdamped cases see Ref. [33]. Because α2 ¼ β2φ2 is small today (it is of the same order as 1 − γPPN in most cases) we can approximate Eq. (14) as 1 − γPPN ¼ 2α2 ¼ 2β2φ2: ð56Þ Evaluating φ in Eq. (55) today and averaging over the small oscillations to extract the secular terms, one finds ð1 − γPPNÞφ0 ∝ β2φ2 R ∝ β; ð57Þ using Eq. (48) for φR. In other words, the overall growth of 1 − γPPN today is linear in β, as we see in the right panel of Fig. 3. Likewise, 1 − βPPN has an extra factor of β in the numerator, making it quadratic, and thus grows even faster than 1 − γPPN, leading to a stronger violation of Solar System constraints for large β. 2. Constraints on β We now turn our attention to placing quantitative constraints on β from the results of our numerical calcu- lations. Because of the oscillatory nature of the scalar field, there are different types of constraints that can be placed, which differ in their generality: (1) βcrit, where all β < βcrit pass Solar System tests at the present time but not necessarily at all future times. (2) βfail, where all β < βfail pass Solar System tests at the present time, and generically also at all future times (this was the method used in Ref. [30]). (3) βallfail, where all β > βallfail will generically fail Solar System tests either today or at some future time. In general, one finds that βallfail > βcrit > βfail and all of these possibilities lead to constraints that are consistent with each other, but for completeness we present all of them SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-9 in our analysis. Furthermore, we will show below that these constraints are relatively insensitive to the value of ζ, i.e. the initial conditions chosen at the time of BBN, provided that they are not finely tuned to be close to ζ ¼ 0 ðφ ¼ 0Þ. The discussion in the previous section specifically applied to the ζ ¼ 1 case only, so now we will focus on the consequences of allowing 0 < ζ < 1. For every set of initial conditions determined by ζ in Eqs. (48)–(51) there exists a critical value of β, denoted by βcrit, for which all β < βcrit lead to a value of φ0 that is consistent with Solar System observations. In general, there is a βcrit associated with each PPN parameter, corresponding to the intersec- tions of the constraints (red/blue horizontal lines) and corresponding curves (red/blue curves) in the right panel of Fig. 3. The smaller of these two will be the value of βcrit reported in this paper, since it places the tightest constraint. For ζ ¼ 1, one can see from the right panel of Fig. 3 that βcrit ¼ 24, corresponding to the black point where the blue horizontal line intersects the curve for j1 − βPPNj. The left panel of Fig. 4 shows the behavior2 of βcrit for all values of ζ, associated with both γPPN and βPPN. Small values of ζ lead to large values of βcrit, such that βcrit → ∞ in the limit that ζ → 0 for both PPN parameters. This should indeed be the case because when ζ → 0 the scalar field is sitting at the minimum of the potential at the time of BBN. If the field is at the minimum in the past, with no velocity, then it will certainly be at the minimum today, and thus it will satisfy Solar System tests with ease. Our results show that regardless of the value of ζ, current constraints on βPPN place the tightest constraints on the value of βcrit. The values of βfail and βallfail are determined in a similar fashion and have a very similar dependence on ζ as βcrit does. We list all of these values in Table II. Ultimately we want to place bounds on values of β and use those in our NS calculations, but ζ ∼ 0 does not restrict β at all according to the arguments presented above. We therefore restrict our calculations to ζ > 0.05 as a way of avoiding a large amount of fine-tuning; the choice of 0.05 as the lower bound is because in a randomly selected sample of ζ there is only a 5% chance that ζ < 0.05 would be selected. In another effort to avoid a finely tuned bound, we take a random distribution of ζ and select the median of the resulting distribution of βcrit to represent a measure of the “most probable” value of βcrit, denoted by βcrit;med, if ζ were chosen randomly. We find that in the exponential theory this median value is βcrit;med ∼ 34. The quantities βcrit and βcrit;med still do not necessarily tell the entire story of which values of β should be considered feasible because, as Fig. 3 shows, there is a (small) subset of β > βcrit (or β > βcrit;med in the case of more finely tuned initial conditions) that do indeed pass Solar System tests today. To quantify this, we define F ¼ #of β that pass Solar System tests #of β sampled ð58Þ and present these results as a function of ζ in the right panel of Fig. 4. The scalar F depends on the sampling of β and ζ, and thus, on their prior probability distribution. We here assume that all β and ζ are equally likely and choose flat distributions for their sampling, which is equivalent to discretizing the β and ζ spaces linearly; we are not aware of any theoretical argument to prefer any particular value of these couplings over others.3 As a reminder, ζ ¼ 1 FIG. 4. Left: Value of βcrit, as defined in the text, for every value of ζ which corresponds to the entire range of possible initial scalar positions at the time of BBN. Right: The total fraction out of all β values sampled (1 < β < 1000) that allow Solar System tests to be passed at τ ¼ 0, as a function of the initial conditions determined by ζ. The general trend of the top right panel tells us that this fraction would continue to decrease as we increase the upper bound of our sample size. 2The step-like nature of βcrit is related to the oscillatory behavior of the scalar field shown in Fig. 3, and it is not an artifact of numerical resolution. 3If a theoretical mechanism is devised in the future that suggests, for example, that ζ must be close to zero at the time of BBN, then the scalar F can be recalculated and the con- clusions found here revised. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-10 corresponds to the solutions displayed in Fig. 3 where βcrit ∼ 24, as determined by constraints on βPPN, and in this case F is only slightly over 8%. This means that there is roughly an 8% chance that a randomly selected value of β between 1 and 1000 will allow Solar System tests to be passed. Interestingly, F changes very slowly with ζ, such that even for very small values of ζ, i.e. for ζ > 0.05, F is still ≲25%. The above results may seem to indicate that there is a decent number of possible values of β above the critical or the median one for which Solar System tests are passed, provided that we sample larger values of β, but this conclusion would be premature. We have here greatly narrowed the scope of our study by only focusing on β < 1000, but in principle, β could take on any value between zero and infinity. Considering values of β > 1000, say up to β ¼ 106 for definiteness, would lead to two consequences: (i) the # of β that pass Solar System tests would increase slightly, and (ii) the # of β sampled would increase dramatically. Because of these effects, the relative fraction F would decrease significantly leaving a very small probability of passing Solar System tests without a very fine-tuned choice of β. Let us wrap up this subsection by summarizing the main results we have obtained, all of which can be inferred from Fig. 3 and Fig. 4. For the ζ ¼ 1 case, we see that for all β ≲ βcrit ∼ 24 the scalar field will evolve such that Solar System tests are passed today. For β ≳ βcrit ∼ 24 the solutions for which Solar System tests are passed is periodic in nature, but in general they do not pass Solar System tests. Moreover, the ones that do pass Solar System tests today do not generically do so in the future. In fact, all β ≳ 29 that happen to pass Solar System tests today will eventually fail them in the future. The monotonic growths seen in the right panel of Fig. 3 imply that this trend will continue as β increases and there will never be a turnover such that larger ranges of β will generally pass Solar System tests. Such conclusions turn out to be valid for all values of ζ > 0. A summary of these results can be found in Table II. D. Hyperbolic theory We now repeat the analysis described above but using the hyperbolic theory instead of the exponential one. The numerical results for the evolution of the scalar field are displayed in Fig. 5. There are several similar features in the scalar field evolution to what we found in the exponential case, as well as a few features exclusive to hyperbolic coupling that we discuss below. The upper panels of Fig. 5 show the evolution of the scalar field for ζ ¼ 1. As β increases, there is a distinct increase in the frequency of oscillation about the minimum of the potential. In the exponential theory, this was attributed to the fact that the potential became steeper with increasing β, generating steeper gradients and therefore larger accelerations. The potential gradient αðφÞ in the hyperbolic theory exponentially asymptotes to a value of 1= ffiffiffi 3 p for large values of βφ, cf. Fig. 2. This behavior limits the maximum value of 1 − γPPN via Eq. (14), and also restricts the coupling between the scalar field and matter, effectively placing an upper limit on the sourcing term of the oscillator. As a result, the scalar field does not gain as TABLE II. Bounds on β using different initial conditions (determined by ζ) and both PPN parameters (γPPN and βPPN) in both the exponential and hyperbolic theories. The parameters βcrit, βfail, βallfail, and F are those defined in the text. The first row shows data for the lowest value of ζ we consider, i.e. ζ ¼ 0.05. The second row shows the data for saturating BBN constraints (ζ ¼ 1). We present the median values of certain parameters in the third row, given by the definitions defined in the text. The last two rows illustrate the amount of fine-tuning of ζ, or alternatively φR, that would be necessary to force βcrit to either 50 or 100. Exponential Hyperbolic ð1 − γPPNÞφ0 ð1 − βPPNÞφ0 ð1 − γPPNÞφ0 ð1 − βPPNÞφ0 ζ ¼ 0.05 βcrit 722 104 395 91 βfail 722 103 394 91 βallfail 811 113 431 100 F 93% 25% 69% 21% ζ ¼ 1 βcrit 42 24 17 17 βfail 41 24 17 17 βallfail 60 29 25 19 F 25% 8% 10% 4% Median Values βcrit 76 34 38 25 βallfail 99 40 49 27 βcrit ¼ 50 ζ < 0.843 ζ < 0.229 ζ < 0.393 ζ < 0.162 φR < 0.063 φR < 0.032 φR < 0.088 φR < 0.041 βcrit ¼ 100 ζ < 0.376 ζ < 0.051 ζ < 0.196 ζ < 0.043 φR < 0.029 φR < 0.011 φR < 0.044 φR < 0.013 SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-11 much velocity as quickly here, and therefore, it spends more time away from the minimum of the potential initially, leading to minimal amounts of oscillatory motion and the saturation of 1 − γPPN occurring in Fig. 5 at early times. Again, while we do not show the time evolution of 1 − βPPN, its behavior is similar and ultimately leads to identical qualitative conclusions. Once enough time has passed and the scalar field has been damped enough to remain near the minimum of the potential, it finally begins to behave as in the exponential theory, but this has a significant effect on whether or not Solar System tests are passed today. Because αðφÞ remains at its limiting value for a substantial amount of time, which is evident for example when β > 10, the friction terms do not have a significant effect until later in the evolution, making it harder for the field to decay to a small enough value in order for either PPN parameter to be consistent with Solar System constraints. This is made clear in the top right panel of Fig. 5 where we see that βcrit ∼ 17, a value smaller than that found in the exponential case. We also see that it becomes even more difficult to pass Solar System tests for large β, i.e. the regions falling below the Solar System bounds become smaller as β becomes larger. We also find that all β ≳ 19 that allow Solar System tests to be passed today will ultimately fail to do so at some point in the future, if Solar System tests continue to show that the constraints on PPN parameters are time independent. Similar to the analysis in the previous section we can study what effect the initial conditions, parametrized via ζ, have on the scalar-field evolution, with the bottom panels of Fig. 5 showing βcrit and F as functions of ζ. The overall behavior of βcrit is similar to that in the exponential theory but it is always significantly lower. Just like we did previously, we can quantify a median value of βcrit corresponding to the median of the distribution appearing in the bottom left panel of Fig. 5 and we find that βcrit;med ∼ 25 for the hyperbolic theory. The fraction F is also smaller and is only around half of what it is in the exponential theory. As we argued before, this fraction will only decrease with an increased sample size of β, and to avoid fine-tuning scenarios one should try to avoid considering values of β larger than βcrit that satisfy Solar System tests today because they just happen to have αðφÞ ¼ 0 today. Table II contains a summary of all bounds on β from these calculations. E. Constraints on β: Summary Let us now briefly summarize our conclusions regarding the bounds we can place on β given Solar System FIG. 5. Same as Figs. 3 and 4 but for the hyperbolic class of STTs. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-12 observations and their implications. First, recall that we focused on the range 1< β< 1000 and 0 < ζ < 1, and that ζ close enough to zero, e.g. ζ < 0.05, corresponds to a finely tuned choice of initial conditions, i.e. choosing φ ∼ 0 at the time of BBN. Table II shows representative values of the bounds on β from the three choices discussed at the beginning of Sec. III C 2, for both PPN parameters in both the exponential and hyperbolic theories. Notice that in every case displayed, the bounds placed from βPPN are all tighter than the ones placed by γPPN. This is consistent with the structure of the PPN constraints because βPPN takes the curvature of the conformal potential into account, while γPPN does not. For any β ≳ 10 the curvature of the potential becomes significant near the minimum of the potential, in turn causing larger present day values of j1 − βPPNj which are inconsistent with Solar System observations. We there- fore focus explicitly on the bounds placed by βPPN. Let us first focus on the ζ ¼ 1 portion of Table II and recall the definitions of βcrit, βfail, and βallfail presented in Section III C 2. For both theories, βcrit, βfail, and βallfail all place nearly identical upper bounds on β. Furthermore,F is 8% and 4% for the exponential and hyperbolic theory respectively. Such small values of F suggest that it is extremely hard to randomly select a value of β > βcrit that will be consistent with Solar System constraints, and in fact this value only becomes smaller by considering values of β > 1000. This suggests that one would have to very precisely select a β > βcrit if one wants to remain consistent with Solar System observations, leading to a fine-tuning problem. Therefore, in order to naturally avoid this, one should only consider β < fβcrit; βfail; βallfailg, all of which are ∼20 in both theories. Let us now discuss the fine-tuning issue in a bit more detail, by focusing on the rest of Table II. First, notice that considering ζ ¼ 0.05 brings the upper bounds on β up to ∼100 and F up to ∼25% for both theories. While this increases the viable range of β, one should remember that forcing ζ to be small is equivalent to forcing φ → 0 at the time of BBN, resulting in a carefully selected set of initial conditions. To avoid this fine-tuning issue, we calculated the median values of the distributions of βcrit and βallfail, which represent the β’s that would be likely to pass Solar System tests if ζ were chosen randomly (these are located in the rows labeled “Median Values” in Table II). Ultimately, we see that this type of analysis only slightly relaxes the bounds on β relative to what is found when ζ ¼ 1, meaning that for a large subset of 0 < ζ < 1 the bounds on β are around ∼40 and ∼25 for the exponential and hyperbolic theory respectively. Finally, let us quantify the amount of fine-tuning that would be needed in order to raise βcrit to a region containing nonzero scalar charge in Fig. 1. This information is shown in the bottom two rows of Table II, but let us here only focus on the hyperbolic theory since we only find stable NSs with nonzero scalar charge there. In order to raise βcrit to 50 (100), ζ would have to be smaller than 0.162 (0.043). To achieve this one needs to demand that the scalar field be close to the minimum of the conformal potential and therefore be very similar to GR on cosmological scales. While allowing STTs to reduce to GR (ζ ∼ 0) on these scales would allow Solar System tests to be passed it is effectively a set of measure zero. For the reasons above, we conclude that the most natural way to evade Solar System constraints today is to consider bounds on β corresponding to ζ ¼ 1, and we note that only under a carefully selected set on initial conditions at BBN and β is it possible to pass Solar System tests otherwise. IV. NEUTRON STARS AND SCALARIZATION In this section we describe NSs in STTs and introduce the basics of scalarization. To reiterate the argument made in Ref. [30], we only study isolated NSs here because theories that do not allow for spontaneous scalarization are also not likely to allow for dynamical/induced scalarization. We first present the equations of structure and then describe the piecewise polytropic EoSs we use to model realistic tabulated EoSs [41]. Finally, once the framework for the numerics has been laid out we present results for NSs in the exponential theory and the hyperbolic theory, first for β < 0 and then for β > 0. For our study, we focus explicitly on isolated, static, and spherically symmetric spacetimes, which in general can be described by the Einstein-frame line element ds2� ¼ −eνðr�Þdt2� þ dr2� 1 − 2μðr�Þ=r� þ r2�dΩ2�; ð59Þ where ν and μ are only functions of the Einstein-frame coordinate radius r� and dΩ2� is the line element on the 2- sphere. To a good approximation, thematter inside coldNSs can be represented as a perfect fluid and thereforewe use the form of the matter stress-energy tensor given in Eq. (12). Applying conservation of stress-energy in the Jordan frame, ∇νTμν ¼ 0, along with the line element and perfect fluid assumptions above yield a set of first-order differential equations governing the structure of the NSs [17,18] μ0 ¼ 4πGr2�A4ðφÞϵþ 1 2 r�ðr� − 2μÞψ2; ð60Þ ν0 ¼ r�ψ2 þ 1 r�ðr� − 2μÞ ½2μþ 8πGr3�A4ðφÞp�; ð61Þ φ0 ¼ ψ ; ð62Þ ψ 0 ¼ 4πGrA4ðφÞ ðr� − 2μÞ ½αðφÞðϵ − 3pÞ þ r�ψðϵ − pÞ� ð63Þ − 2ψ ð1 − μ=r�Þ ðr� − 2μÞ ; p0 ¼ −ðϵþ pÞðν0=2þ αðφÞψÞ; ð64Þ SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-13 where primes denote derivatives with respect to r�, and all occurrences of p and ϵ are explicitly in the Jordan frame.4 To close the system of equations we need to prescribe an EoS. We use the piecewise polytropic EoSs studied in Ref. [41] in which the authors developed a four-parameter polytropic model that accurately approximates 34 tabulated EoSs. Let us briefly discuss this parametrization of the EoS. In each baryonic density region inside the star ρi−1 < ρ < ρi the pressure is given by the standard polytropic form p ¼ Kiρ Γi ; ð65Þ where Γi is the adiabatic index of the fluid in that region and Ki is the polytropic constant chosen to ensure con- tinuity at the boundary between density regions. The energy density ϵ appearing in the structure equations is then found by integrating the first law of thermodynamics d ϵ ρ ¼ −pd 1 ρ ; ð66Þ to find ϵ ¼ ρþ p Γ − 1 : ð67Þ We adopt a low-density crust region described by a degenerate relativistic gas with adiabatic constant Γ0 ¼ 4=3 and match it to a polytrope with adiabatic index Γ1. The polytropic constant K0 for this crust region is determined by forcing the crust to match the SLy EoS at a baryonic density of ρ0 ¼ 2.7 × 1014 g=cm3. A different choice of a fixed crust EoS only affects the physical quantities of the system by a few percent. More importantly, scalarization will only occur in the dense regions of the star, which are unaffected by the choice of crust used, so the simple choice we make here will suffice for our study. The boundary between the crust and the next polytropic region will therefore depend on the value of Γ1. At a fixed baryonic density of ρ1 ¼ 1014.7 g=cm3 and pressure p1 ¼ pðρ1Þ the EoS is matched to another polytrope with adiabatic index Γ2 and constant K2. Another boundary is set at ρ2 ¼ 1015 g=cm3 where the EoS is matched yet again to another polytrope with index Γ3 and constant K3. The set of parameters flogðp1Þ;Γ1;Γ2;Γ3g determines the best fit to tabulated EoSs and were determined in Ref. [41], which we refer the reader to for a more detailed description. Following Ref. [20] we restrict the scope of our inves- tigation to EoSs that allow for a maximum gravitational mass of more than two solar masses and do not violate causality. These constraints leave us with five EoSs, whose parameters are listed in Table III for convenience. We numerically solve the system of equations above parametrized by the central density ρc and coupling constant β. Following the methods employed in Ref. [30] we use MATHEMATICA’s default integrator which makes use of a LSODA approach (a variant of the Livermore Solver for Ordinary Differential Equations) [42]. Inside the star the pressure is governed by the piecewise polytropic EoS, while the pressure vanishes outside the star with the boundary where the pressure vanishes marking the surface of the star. We integrate the equations from the center of the star to an effective infinity that is far from the surface of the star, which allows us to determine the metric and the scalar field in the entire spacetime. Far from the star, the scalar field behaves as φ ¼ φ∞ þ ω r� þO � 1 r2� � ; ð68Þ where φ∞ is the asymptotic value of the scalar field determined from cosmology (which we fix to a small value and enforce as a boundary condition) andω is a constant that is proportional to the scalar charge. We define the scalar charge via αsc ¼ ω Gm� ; ð69Þ wherem� is the Einstein-frameADMmass of the star and the subscript “sc” is used to distinguish the scalar charge from the conformal couplingαðφÞ. Physically, the scalar charge is a measure of the effective coupling between the scalar field and the star, with a larger value meaning stronger coupling. The scalar charge is then determined by ωwhich we find by TABLE III. Parameters of the piecewise polytropes flogðp1Þ;Γ1;Γ2;Γ3g from Ref. [41] for the five EoSs we use for NS calculations and the values of ðβmin; βmax; αsc;maxÞ found here. For clarity, all ρ > ρcrit allow Tmat > 0 in the core, βmin is the lowest value of β allowing for stable scalarized solutions, βmax is where the maximal stable scalar charge occurs, and αsc;max is the maximal scalar charge, corresponding to the peak values of the curves in Fig. 1. EoS logðp1Þ Γ1 Γ2 Γ3 ρcrit=ρ0 βmin βmax αsc;max ENG 34.437 3.514 3.130 3.168 4.544 19.1 26.7 ∼1.4 × 10−3 MPA1 34.495 3.446 3.572 2.887 3.683 19.3 27.2 ∼1.4 × 10−3 MS1 34.858 3.224 3.033 1.325 3.067 35.6 46.3 ∼2.7 × 10−4 MS1b 34.855 3.456 3.011 1.425 3.092 35.8 48.3 ∼2.4 × 10−4 SLy 34.384 3.005 2.988 2.851 5.349 30.0 39.3 ∼4.2 × 10−4 4There is a missing factor of r� in the Eq. (63) equivalent appearing in Ref. [30]. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-14 extracting the 1=r piece of the external solution of the scalar field. A. β < 0: Standard spontaneous scalarization Choosing a sufficiently negative value of β (≲ − 4 for most EoSs) leads to an instability in the star that causes an activation of the scalar field above its background cosmo- logically determined value, a process known as sponta- neous scalarization [17,18,23]. In the context of cosmology, however, all massless STTs with β < 0 generi- cally lead to drastic disagreement with Solar System observations [30,31], at least for the subset of symmetric, monotonically increasing coupling potentials that we con- sider here. Thus, the NS solutions presented in this section should be strictly seen as a stepping stone towards giving a better understanding of what will be presented for the β > 0 cases. Figure 6 shows the results of our numerical calculations using both exponential and hyperbolic couplings for β ¼ −6 and β ¼ −5.5 as examples. For these calculations we employ the ENG piecewise polytropic EoS with parameters given in Table III, but it should be noted that the qualitative behavior seen here is general among all EoSs that we consider. At low densities, only the GR solution exists and the scalar field remains at its background value. There exists a critical density at which the scalar field begins to “turn on” and produces sudden deviations from GR as can be seen in the left panel of Fig. 6. The critical density at which this happens depends on the EoS and the value of β, but as we stated earlier, it is a generic feature among all EoSs considered here. Larger central densities eventually lead back to the GR solution and a trivial scalar field profile, just as in the low-density case. A more negative value of β, which in these theories explicitly means a stronger coupling tomatter, leads to larger deviations from GR and therefore larger values of scalar charge αsc (cf. the right panel of Fig. 6). This is consistent with the idea that scalar charge is the effective coupling between the star and the scalar field.Notice, however, that the exponential theory predicts larger deviations from GR than the hyperbolic theory does (i.e. the charge is roughly 3 times larger), regardless of what value of β we choose. This can be understood by looking at the first term in the square brackets in Eq. (63), or more explicitly the fact that ψ 0 ∝ αðφÞ, which provides a type of feedback loop. In the exponential theory, large values of φ produce large values of αðφÞ relative to the hyperbolic theory, which is obvious from the right panel of Fig. 2. Larger values of αðφÞ thus lead to a larger positive feedback into ψ 0 and so on, explaining why the exponential theory leads to a larger activation of the scalar field for given values of the central density. It is useful to provide an analytic explanation behind this spontaneous activation of the scalar field. In both theories, αðφÞ ¼ βφþOðφ2Þ for small values of βφ, which must be true due to the arguments presented in Sec. III. Following arguments and notation presented in Refs. [18,30], consider then the evolution of the scalar field in a weakly gravitating, spherically symmetric system such that □ → δij∇i∇j → ∇2 r , where ∇2 r is the radial part of the Laplacian in spherical coordinates. For simplicity we will assume T� mat is constant and while weakly gravitating systems would be expected to have T� mat < 0 we allow it to remain arbitrary in the interest of generality. These assumptions reduce the field equation for the scalar field to ∇2 rφ ¼ −K2φsignðβT� matÞ; ð70Þ FIG. 6. Left: Baryonic mass of the NS as a function of central baryonic density. Deviations from GR occur above some critical density which depends on β. A more negative β leads to larger deviations from GR, but notice that the exponential theory has more of an effect than the hyperbolic theory does. Right: Scalar charge as a function of the total baryonic mass. The “spontaneous” activation of the scalar field occurs at a particular mass, while the dependence on β appears to be independent of the theory. Again we see that the exponential theory causes a larger activation of the scalar field than the hyperbolic theory. The circles in both panels correspond to the solutions that have maximum scalar charge and the crosses indicate the maximum mass solutions, which also mark an instability in these branches of solutions [24,43]. SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-15 where K2 ¼ ðκ=2ÞjβT� matj for r� < R and must vanish outside the star. The solution to this must still obey the physical boundary conditions of the scalar field, i.e. φð0Þ ¼ φc and φ0ð0Þ ¼ 0 to ensure regularity at the center and the solution must be continuous and differentiable at the surface r� ¼ R. When the product βT� mat is positive the solution becomes φ ¼ φ∞ cosðKRÞ sinðKrÞ Kr : ð71Þ Even when φ∞ ≈ 0, which we emphasize is the case here, a nontrivial scalar-field solution can arise when KR ≈ π=2. This type of “resonance” between the scalar field and matter explains why spontaneous scalarization occurs in DEF-like theories and why there exist deviations from GR in Fig. 6. However, when βT� mat is negative the scalar field is exponentially suppressed to its background value, which can be seen by replacing the trig functions in Eq. (71) with their corresponding hyperbolic counterparts. It is then reasonable to expect that only certain ranges of density, which ultimately controls T� mat, will lead to an activation of the scalar field for any given β, and why for relatively small and large central densities we only see the GR solution in Fig. 6. We here focused explicitly on the sign of βT� mat rather than just β for a reason. Indeed, for the weakly gravitating system used in the example above T� mat < 0 and the sign in Eq. (71) would be completely determined by the sign of β. However, cosmological constraints force β > 0, so unless T� mat were also positive one would not expect scalarization to occur. This was studied recently in Refs [19,20,29] and will be the focus of the next section. B. β > 0: Go big or go home Studying NSs in STTs with β > 0 is important because this is the region of parameter space that is consistent with Solar System constraints upon cosmological evolution of the scalar field. Luckily, NSs are extremely dense and there can exist regions near the core where T� mat is indeed positive, allowing for possible scalarized branches of solutions. The authors of Ref. [19], however, demonstrated that NSs that scalarize in the exponential theory lead to gravitational collapse. We will start our discussion by summarizing the results presented in Refs. [19,29] and then present the results of our numerical study of NSs in the hyperbolic theory with β > 0. 1. Gravitational collapse in exponential theory The instability of the exponential theory can be under- stood through a linear stability analysis of Eq. (8). By allowing φ ¼ φ0 þ δφ and keeping terms up to first order one finds □�δφ ¼ − κ 2 β0T� matδφ; ð72Þ where β0 ¼ ð∂α=∂φÞφ0 is like the curvature of the coupling potential,5 and the box operator and the stress-energy tensor are those of GR. The authors of Ref. [19] performed a numerical search for solutions to Eq. (72) of the form δφ ¼ eΩtfðrÞwith a spherically symmetric background, a perfect fluid stress-energy tensor, and a polytropic EoS (in general, the location of the instability regions will depend on the EoS used). The results of this search are presented in Fig. 2 of Ref. [19] and we provide a rough schematic drawing of those results in Fig. 7.6 The results of the stability analysis are independent of the theory being studied and only depend on the value of β0 defined above. The red/blue regions in Fig. 7 correspond to regions that are unstable to scalar field perturbations, which could lead to a nontrivial scalar-field profile inside the star. White regions, on the other hand, are regions in the parameter space where only the GR solution is present and no excitations of the scalar field would be expected inside the star. We mark a line at β ¼ −6 to clearly show that as central density increases the instability arises but then as FIG. 7. A rough schematic representation of the instability regions from the linear stability analysis carried out in Ref. [19]. Red/blue regions are unstable to scalar-field perturbations and the black region corresponds to solutions that are unstable to gravitational collapse. The horizontal line marks the density above which T� mat > 0 at the center of the NS. We mark a vertical line at β ¼ −6 (corresponding to an example scenario we used in Fig. 6) to help demonstrate the behavior of the instability. 5Notice that Eq. (72) encompasses a large class of STTs, not just the exponential model. This analysis applies to STTs whose coupling functions αðφÞ can be approximated by βφþOðφ2Þ. 6Figure 7 is indeed a quantitative replication of the results formally presented in Ref. [19]. We adjusted the regions to match those of the ENG EoS studied here. We only present this figure here for completeness and to act as a visual aid to explain the behavior of the scalar instability. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-16 density continues to rise the instability eventually turns off. The turning on/off of the instability in the β < 0 regime is consistent with the results we previously presented in Fig. 6 for both theories. Recall that the instability is related to the sign of ðβT� matÞ as we argued in the previous section. For small ρc and β < 0 the sign of ðβT� matÞ is positive, causing the amplification of the scalar field inside the star. As ρc grows larger eventually T� mat becomes positive, causing ðβT� matÞ to become negative, ultimately suppressing the scalar field rather than amplifying it, which is why the instability eventually vanishes as ρc increases. In the β > 0 region of the parameter space the behavior of the instability is slightly different. The instability to scalar perturbations no longer vanishes for larger values of ρc, meaning that once the central density becomes large enough, i.e. in the blue regions of Fig. 7, the instabilitywill be present inside the star. The work in Refs. [19,29] showed that the final fate of scalarized NSs existing in the β > 0 region of Fig. 7 in the exponential theory was always gravitational collapse. Our understanding of these results is that NSs can exist in the exponential theorywith β > 0, provided that they live outside of the blue region of Fig. 7 and are identical to the trivial GR solutions. However, since these NSs would be identical toNSs inGR in these regions of parameter space no scalarized solutions would exist in nature. Because NSs in the exponential theory that could undergo scalarization would ultimately collapse to black holes we focus our attention on the study of NSs in the hyperbolic theory. 2. Scalarization in the hyperbolic theory Following the numerical procedure explained at the beginning of the section, we calculate NS solutions in all five EoSs given in Table III for 0 ≤ β ≲ 450 and β ¼ 1000. While we do not analyze the dynamical stability of the static solutions we find, we define stability based on a turning-point criterion in a sequence of equilibrium con- figurations described in detail in Refs. [24,43]. In short, every solution to the left of the maximum mass in a mass- density plot, like that in Fig. 8, would be considered stable and every solution to the right would be unstable. Figure 8 shows our numerical solutions using the ENG EoS for a small subset of β, with the top panel demonstrating stronger deviations from GR for larger values of β. There are three key features of note in these solutions: (1) There exist values of β that have no stable branch of solutions. (2) There appears to be an asymptotic branch of solutions for large values of β. (3) There will exist a maximum mass and stable scalar charge αsc;max for each β. Let us investigate the first point above. From our turning- point definition of stability, any solution existing on a branch to the right of the maximum will be unstable to gravitational collapse. Because small values of β do not lead to deviations from GR until after the maximummass is reached, it stands to reason that there exists a βmin such that only values of β > βmin allow for stable scalarized NS solutions in the hyperbolic theory. For example, in the top panel of Fig. 8 the β ¼ 10 curve clearly does not have stable solutions different from GR. In other words, suffi- ciently small values of β only lead to deviations from GR when the ρc is large and in the black region of Fig. 7 and would undergo gravitational collapse. When β ¼ 40 there exists a clear set of solutions different from GR that would indeed be stable up to the turning point, which we indicate with a circle in the figure for clarity. The values of βmin are listed in Table III for all EoSs studied. There also appears to be an asymptotic branch of solutions in the limit β → ∞, a qualitative feature occurring FIG. 8. Top: Baryonic mass as a function of central baryonic density for a wide range of β. There exists an asymptotic solution as β → ∞ (as β increases the branches only experience small changes from one to the next). Solid points indicate the maximum mass that is stable and the solid red line connects all maximum mass solutions for all β. Bottom: The scalar charge as a function of the central baryonic density for the same set of β’s as in the top panel. Solid points here correspond to the same solutions as the solid points in the top panel, with the solid red curve showing how these values change as a function of β. Even though the scalar charge grows monotonically there exists a maximum value that it can take for stable NS solutions. For large β, increasing β continues to decrease the maximum stable scalar charge, though it does not cause the charge to vanish. SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-17 in each EoS studied. Starting from β ¼ 0, increasing β initially induces dramatic changes in the corresponding branch of solutions; however, for β ≳ 40 for the ENG EoS, the relative change between solution branches changes only slightly even when β is increased by an order of magnitude. This idea can be explained by again referring back to the structure of αðφÞ in Fig. 2. As β → ∞ we find that αðφÞ → 1= ffiffiffi 3 p exponentially. Therefore, as long as φ does not identically vanish everywhere, the coupling to matter is saturated, and thus, there will exist a branch of solutions representing this maximal coupling case, which we see in Fig. 8 as β increases. This feature is absent in the exponential theory because αðφÞ is linear in β and therefore it has no asymptotic limit, meaning that the scalar field’s coupling to matter can never be saturated. Stability arguments [24,43] indicate that there exists a maximum mass, and therefore a maximum stable scalar charge, on each branch of solutions determined by β. However, the bottom panel of Fig. 8 indicates that the scalar charge grows monotonically as the central density increases. Using the turning point criterion allows us to extract the most massive stable solution and hence the corresponding scalar charge of that solution. The points marked in Fig. 8 correspond to identical solutions in each panel and for β ≳ 40 we see that the maximum stable scalar charge on each solution branch decreases with increasing β. There actually exists a βmax where αsc;maxðβÞ peaks, which is evident in Fig. 1 and whose values appear in Table III. Figure 1 shows the maximum stable scalar charge as a function of β for all EoSs in Table III. There is a general trend that EoSs that allow for larger compactness have large values of scalar charge. This complements the results presented in Ref. [20] where the author demonstrated a relation between the compactness and the onset of scala- rization for the same set of EoSs. We can also roughly quantify how rapid the scalar field “turns on”7 by meas- uring the slope of the charge in the bottom panel of Fig. 8 near the maximum stable solutions, i.e. the ones marked with points. Figure 9 shows the value of this slope, which we denote as α0sc;max ¼ ∂αsc;max=∂ρc, as a function of β. For small β there is no change in the charge because, as we previously stated, there are no stable scalarized solutions below βmin. As one may expect from Fig. 8 there is a clear monotonic decrease in α0sc;max as we increase β. A decreas- ing α0sc;max means that the scalar field grows much more gradually with central pressure as β becomes larger. To understand why the maximum scalar charge decreases for large β we must return to the field equations, particularly Eq. (63) which we rewrite for convenience as ψ 0 ¼ 4πGrA4ðφÞ ðr� − 2μÞ ½αðφÞðϵ − 3pÞ þ r�ψðϵ − pÞ� − 2ψ ð1 − μ=r�Þ ðr� − 2μÞ : This equation describes the curvature of the scalar field profile, with the first term in square brackets identified with ð−αðφÞT� matÞ. This term will explicitly be negative near the center of the NS because T� mat > 0. Increasing β increases the value of αðφÞ, albeit only slightly when β is large. This term then contributes to a more negative value of curvature, which causes the scalar field to decay faster inside the star as β increases. Numerically, we find that the NS solutions appearing in Fig. 1 have nearly identical central scalar field values for β > 100. Because the scalar field decays faster inside the star there will exist a smaller value of ω from Eq. (68) upon matching boundary conditions at the surface and therefore a smaller scalar charge. V. CONCLUSION Our calculations show that theories with an exponential conformal factor (identical to the ones first studied by Damour and Esposito-Farèse [17,18]) and theories with a hyperbolic conformal factor (motivated from quantum field theory [12,20,22]) always pass Solar System constraints provided that β ≲ βcrit ¼ 24 and β ≲ βcrit ¼ 17 respectively for all initial conditions consistent with BBN constraints. By considering a random distribution of initial conditions FIG. 9. Rate of change of the scalar charge evaluated at the maximum mass solutions (solutions identified by points in Fig. 8) as a function of β for the ENG EoS. There are no values for β ≲ 19 because there are no stable solutions with scalar charge, cf. the discussion in the text. Notice that the “rapidity” of the increase of scalar charge monotonically decreases with increasing β. The apparent kink near β ¼ 25 is a result of numerical error due to the fact that we only have limited precision and must interpolate between different solution branches. 7This activation, or “turn on,” is not dynamical. We are simply referring to how rapid the charge increases as ρc increases. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-18 below the BBN constraint, these bounds can be relaxed to βcrit ∼ 34 and βcrit ∼ 25 respectively. These results tell us that Solar System constraints place a rather tight bound on the β parameter if one wishes to avoid fine-tuning the initial conditions at the time of BBN. Recent work has suggested that scalarization is still possible in STTs with β > 0 if the star becomes dense enough such that the trace of the matter stress-energy tensor becomes positive inside the star, which is only possible for certain EoSs [19,20,29]. However, neutron stars in STTs with an exponential conformal factor have recently been shown to undergo gravitational collapse when β > 0 [19]; therefore, only NSs with Tmat ¼ −ρþ 3p < 0 would exist in this theory and these would not scalarize. After verifying the results in these recent studies, we constructed NSs in a theory with hyperbolic coupling and β > 0 for β regions consistent with Solar System observations. Our results show that there exists a minimum value of β, denoted as βmin in Table III, for which there exist stable scalarized NS solutions that are distinct from GR. However, in every EoS that we considered, all of these values of βmin were larger than the upper bound on β from Solar System observations when one saturates BBN initial conditions. The fact that βcrit < βmin for many EoSs does not mean that no scalarized NSs can exist. Indeed, one can choose more finely tuned initial conditions at the time of BBN, thus making βcrit slightly larger. The larger one wishes to make βcrit, the more one will need to fine-tune the initial conditions. Eventually, such a large fine-tuning of initial conditions would have to be justified through a theoretical physics mechanism that currently does not exist. Alternatively, one can pick a value of β > βcrit such that Solar System constraints are satisfied today but fail in the future. In this case, the likelihood of the existence of stable, scalarized NSs becomes smaller the larger β is above βcrit, as the bands of allowed β thin out and become more sparse as β increases. The results presented here and in Ref. [30] may beg the question of whether or not it is possible to engineer a conformal factor in massless STTs that generally passes Solar System tests and, furthermore, leads to scalarization in neutron stars. Solar System tests demand that, at least locally, the conformal coupling potential have positive curvature, i.e. βðφÞ > 0 in general, and therefore one must have an EoS that allows Tmat > 0 in order to have scalarized neutron stars. It seems that having βðφÞ > 0 and Tmat > 0 only allows for a scalar charge of Oð10−3Þ, whereas in the βðφÞ < 0 case the charge is of Oð10−1Þ. Therefore, these consistencies between Solar System observations and the cosmological evolution of the scalar field may limit the maximum size of the scalar charge that can occur in NSs for a large class of massless STTs. A promising avenue for future research would be the study of massive STTs [44,45] in this context or to perform an analysis like the one in this paper on tensor-multiscalar theories. If consistency between the bounds from Solar System tests and scalarized NSs were achieved, then NSs in these theories would have a different mass and radius than that predicted in GR and only deviate from GR after a critical density (or compactness) is achieved such that the trace of the stress-energy tensor is positive. A direct observation of such a system, perhaps with NASA’s NICER telescope [46,47], in conjunction with gravitational-wave observation of such systems, may allow further constraints to be placed on STTs. Both observations would be necessary in order to break degeneracies between the mass-radius relation, the tidal deformabilities and our ignorance of the equation of state through approximately universal relations [48–51]. Such studies would provide useful information on what the actual form of the conformal factor may be. We have seen that some models, like the exponential one, do not allow for scalarized NSs, while others, like the hyperbolic one, do. Therefore, any observation of neutron stars that could shed some light on the functional form of the conformal factor would be valuable. Similarly, it would also be interesting to investigate dynamical and induced scalarization in hyperbolic STTs with β > 0, in particular for gravitational-wave tests with ground-based instruments [52–54]. Given the results obtained in the spontaneous scalarization case, we expect the scalar charge to be dynamically generated in binaries, but its magnitude to be significantly smaller than in the exponential β < 0 case. The latter case, however, is already ruled out by Solar System observations, while we showed here that hyperbolic STTs with β > 0 can still be con- sistent. Studying scalarization in these theories would allow us to predict the signatures that would imprint on the gravitational waves emitted in neutron star binaries. Presumably, the activation of the scalar field will induce dipole radiation, which in turn will speed up the inspiral of the binary, leaving a signature on the gravitational waves emitted. If so, gravitational-wave observations with Advanced LIGO could place the first meaningful con- straints on such theories in the context of NSs. ACKNOWLEDGMENTS We thank Emanuele Berti, Enrico Barausse, Norbert Wex, Hector Okada da Silva, Raissa Mendes, and Nestor Ortiz for useful discussions on neutron stars in scalar-tensor theories, and for reading our manuscript carefully and providing useful suggestions. N. Y. acknowledges support from the NSF CAREER Grant No. PHY-1250636 and NASA Grant No. NNX16AB98G. SOLAR SYSTEM CONSTRAINTS ON MASSLESS SCALAR- … PHYSICAL REVIEW D 96, 064037 (2017) 064037-19 [1] C. M. Will, Living Rev. Relativ. 17, 4 (2014). [2] B. Bertotti, L. Iess, and P. Tortora, Nature (London) 425, 374 (2003). [3] P. C. C. Freire, N. Wex, G. Esposito-Farese, J. P. W. Verbiest, M.Bailes, B. A. Jacoby,M.Kramer, I. H. Stairs, J.Antoniadis, and G. H. Janssen, Mon. Not. R. Astron. Soc. 423, 3328 (2012). [4] B. P. Abbott et al. (Virgo and LIGO Scientific Collabora- tions), Phys. Rev. Lett. 116, 061102 (2016). [5] B. P. Abbott et al. (Virgo and LIGO Scientific Collabora- tions), Phys. Rev. Lett. 116, 221101 (2016). [6] R. V. Wagoner, Phys. Rev. D 1, 3209 (1970). [7] K. J. Nordtvedt, Astrophys. J. 161, 1059 (1970). [8] P. G. Bergmann, Int. J. Theor. Phys. 1, 25 (1968). [9] M. Salgado, Classical Quantum Gravity 23, 4719 (2006). [10] J. Polchinski, String Theory: Volume 1, An Introduction to the Bosonic String (Cambridge University Press, Cambridge, England, 1998). [11] M. Green, J. Schwarz, and E. Witten, Superstring Theory: Volume 1, Introduction (Cambridge University Press, Cambridge, England, 1988). [12] N. Birrell and P. Davies, Quantum Fields in Curved Space (Cambridge University Press, Cambridge, England, 1984). [13] S. Weinberg, Phys. Rev. D 77, 123541 (2008). [14] V. Faraoni, Cosmology in Scalar-Tensor Gravity (Springer, New York, 2004). [15] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. 513, 1 (2012). [16] J. Martin, C. Schimd, and J.-P. Uzan, Phys. Rev. Lett. 96, 061303 (2006). [17] T. Damour and G. Esposito-Farese, Classical Quantum Gravity 9, 2093 (1992). [18] T. Damour and G. Esposito-Farese, Phys. Rev. Lett. 70, 2220 (1993). [19] R. F. P. Mendes and N. Ortiz, Phys. Rev. D 93, 124035 (2016). [20] R. F. P. Mendes, Phys. Rev. D 91, 064024 (2015). [21] P. Pani, V. Cardoso, E. Berti, J. Read, and M. Salgado, Phys. Rev. D 83, 081501 (2011). [22] W. C. C. Lima, G. E. A. Matsas, and D. A. T. Vanzella, Phys. Rev. Lett. 105, 151102 (2010). [23] T. Damour and G. Esposito-Farese, Phys. Rev. D 54, 1474 (1996). [24] T. Harada, Phys. Rev. D 57, 4802 (1998). [25] C. Palenzuela, E. Barausse, M. Ponce, and L. Lehner, Phys. Rev. D 89, 044024 (2014). [26] E. Barausse, C. Palenzuela, M. Ponce, and L. Lehner, Phys. Rev. D 87, 081506 (2013). [27] M. Shibata, K. Taniguchi, H. Okawa, and A. Buonanno, Phys. Rev. D 89, 084005 (2014). [28] K. Taniguchi, M. Shibata, and A. Buonanno, Phys. Rev. D 91, 024033 (2015). [29] C. Palenzuela and S. L. Liebling, Phys. Rev. D 93, 044009 (2016). [30] D. Anderson, N. Yunes, and E. Barausse, Phys. Rev. D 94, 104064 (2016). [31] L. Sampson, N. Yunes, N. Cornish, M. Ponce, E. Barausse, A. Klein, C. Palenzuela, and L. Lehner, Phys. Rev. D 90, 124091 (2014). [32] T. Damour and K. Nordtvedt, Phys. Rev. Lett. 70, 2217 (1993). [33] T. Damour and K. Nordtvedt, Phys. Rev. D 48, 3436 (1993). [34] B. Bertotti, L. Iess, and P. Tortora, Nature (London) 425, 374 (2003). [35] C. M. Will and K. J. Nordtvedt, Astrophys. J. 177, 757 (1972). [36] K. J. Nordtvedt and C. M. Will, Astrophys. J. 177, 775 (1972). [37] C. Will, Theory and Experiment in Gravitational Physics (Cambridge University Press, Cambridge, England, 1993). [38] J.-P. Uzan, Living Rev. Relativ. 14, 2 (2011). [39] A. Coc, K. A. Olive, J.-P. Uzan, and E. Vangioni, Phys. Rev. D 73, 083525 (2006). [40] A. Coc, K. A. Olive, J.-P. Uzan, and E. Vangioni, Phys. Rev. D 79, 103512 (2009). [41] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009). [42] Wolfram Research Inc., MATHEMATICA 10.4.0 (2016). [43] T. Harada, Prog. Theor. Phys. 98, 359 (1997). [44] T. Arnoulx de Pirey Saint Alby and N. Yunes, arXiv:1703.06341. [45] F. M. Ramazanoğlu and F. Pretorius, Phys. Rev. D 93, 064005 (2016). [46] K. C. Gendreau, Z. Arzoumanian, and T. Okajima, Proc. SPIE Int. Soc. Opt. Eng. 8443, 844313 (2012). [47] Z. Arzoumanian et al., Proc. SPIE Int. Soc. Opt. Eng. 9144, 914420 (2014). [48] K. Yagi and N. Yunes, Science 341, 365 (2013). [49] K. Yagi and N. Yunes, Phys. Rev. D 88, 023009 (2013). [50] K. Yagi and N. Yunes, Classical Quantum Gravity 33, 13LT01 (2016). [51] K. Yagi and N. Yunes, Phys. Rep. 681, 1 (2017). [52] K. G. Arun and A. Pai, Int. J. Mod. Phys. D 22, 1341012 (2013). [53] X. Zhang, J. Yu, T. Liu, W. Zhao, and A. Wang, Phys. Rev. D 95, 124008 (2017). [54] L. Shao, N. Sennett, A. Buonanno, M. Kramer, and N. Wex, arXiv:1704.07561. DAVID ANDERSON and NICOLÁS YUNES PHYSICAL REVIEW D 96, 064037 (2017) 064037-20 https://doi.org/10.12942/lrr-2014-4 https://doi.org/10.1038/nature01997 https://doi.org/10.1038/nature01997 https://doi.org/10.1111/j.1365-2966.2012.21253.x https://doi.org/10.1111/j.1365-2966.2012.21253.x https://doi.org/10.1103/PhysRevLett.116.061102 https://doi.org/10.1103/PhysRevLett.116.221101 https://doi.org/10.1103/PhysRevD.1.3209 https://doi.org/10.1086/150607 https://doi.org/10.1007/BF00668828 https://doi.org/10.1088/0264-9381/23/14/010 https://doi.org/10.1103/PhysRevD.77.123541 https://doi.org/10.1016/j.physrep.2012.01.001 https://doi.org/10.1016/j.physrep.2012.01.001 https://doi.org/10.1103/PhysRevLett.96.061303 https://doi.org/10.1103/PhysRevLett.96.061303 https://doi.org/10.1088/0264-9381/9/9/015 https://doi.org/10.1088/0264-9381/9/9/015 https://doi.org/10.1103/PhysRevLett.70.2220 https://doi.org/10.1103/PhysRevLett.70.2220 https://doi.org/10.1103/PhysRevD.93.124035 https://doi.org/10.1103/PhysRevD.93.124035 https://doi.org/10.1103/PhysRevD.91.064024 https://doi.org/10.1103/PhysRevD.83.081501 https://doi.org/10.1103/PhysRevD.83.081501 https://doi.org/10.1103/PhysRevLett.105.151102 https://doi.org/10.1103/PhysRevD.54.1474 https://doi.org/10.1103/PhysRevD.54.1474 https://doi.org/10.1103/PhysRevD.57.4802 https://doi.org/10.1103/PhysRevD.89.044024 https://doi.org/10.1103/PhysRevD.89.044024 https://doi.org/10.1103/PhysRevD.87.081506 https://doi.org/10.1103/PhysRevD.87.081506 https://doi.org/10.1103/PhysRevD.89.084005 https://doi.org/10.1103/PhysRevD.91.024033 https://doi.org/10.1103/PhysRevD.91.024033 https://doi.org/10.1103/PhysRevD.93.044009 https://doi.org/10.1103/PhysRevD.93.044009 https://doi.org/10.1103/PhysRevD.94.104064 https://doi.org/10.1103/PhysRevD.94.104064 https://doi.org/10.1103/PhysRevD.90.124091 https://doi.org/10.1103/PhysRevD.90.124091 https://doi.org/10.1103/PhysRevLett.70.2217 https://doi.org/10.1103/PhysRevLett.70.2217 https://doi.org/10.1103/PhysRevD.48.3436 https://doi.org/10.1103/PhysRevD.48.3436 https://doi.org/10.1038/nature01997 https://doi.org/10.1038/nature01997 https://doi.org/10.1086/151754 https://doi.org/10.1086/151754 https://doi.org/10.1086/151755 https://doi.org/10.1086/151755 https://doi.org/10.12942/lrr-2011-2 https://doi.org/10.1103/PhysRevD.73.083525 https://doi.org/10.1103/PhysRevD.73.083525 https://doi.org/10.1103/PhysRevD.79.103512 https://doi.org/10.1103/PhysRevD.79.103512 https://doi.org/10.1103/PhysRevD.79.124032 https://doi.org/10.1143/PTP.98.359 http://arXiv.org/abs/1703.06341 https://doi.org/10.1103/PhysRevD.93.064005 https://doi.org/10.1103/PhysRevD.93.064005 https://doi.org/10.1117/12.926396 https://doi.org/10.1117/12.926396 https://doi.org/10.1117/12.2056811 https://doi.org/10.1117/12.2056811 https://doi.org/10.1126/science.1236462 https://doi.org/10.1103/PhysRevD.88.023009 https://doi.org/10.1103/PhysRevD.88.023009 https://doi.org/10.1088/0264-9381/33/13/13LT01 https://doi.org/10.1088/0264-9381/33/13/13LT01 https://doi.org/10.1016/j.physrep.2017.03.002 https://doi.org/10.1142/S0218271813410125 https://doi.org/10.1142/S0218271813410125 https://doi.org/10.1103/PhysRevD.95.124008 https://doi.org/10.1103/PhysRevD.95.124008 http://arXiv.org/abs/1704.07561