Scalar charges and pulsar-timing observables in the presence of nonminimally coupled scalar fields
Raissa F. P. Mendes Tulio Ottoni





rfpmendes@id.uff.br
Instituto de Física, Universidade Federal Fluminense, Niterói, 
Rio de Janeiro, 24210-346, Brazil.

tulioottoni@id.uff.br
Instituto de Física, Universidade Federal Fluminense, Niterói, 
Rio de Janeiro, 24210-346, Brazil.





Pulsar-timing has become a celebrated tool for probing modifications to general relativity in the strong-field surroundings of neutron stars. Here we investigate whether scalar-tensor theories that incorporate a nonminimally coupled scalar degree of freedom may pass pulsar-timing tests, by computing the scalar charges entering such observables. In particular we show that for positive values of the nonminimal coupling ξ, pulsar-timing constraints may be evaded even in the presence of spontaneous scalarization.


04.50.Kd, 04.40.Dg, 04.80.Cc, 97.60.Gb 

Scalar charges and pulsar-timing observables in the presence of nonminimally coupled scalar fields
    Tulio Ottoni
    December 30, 2023
==================================================================================================






§ INTRODUCTION


Perhaps the simplest way in which general relativity (GR) can be modified is through the coupling to a new scalar degree of freedom Φ, as accomplished by the general scalar-tensor action

    S_g  = 1/16π G∫ d^4 x √(-g)[ F(Φ) R -  Z(Φ) ∇_μΦ∇^μΦ - V(Φ) ].

This class of theories naturally implements Dirac's idea of varying fundamental “constants” <cit.>, as G_eff≡ G/F(Φ) can be interpreted as a time- and space-dependent effective gravitational coupling. A well-motivated form for F(Φ) is the standard nonminimal coupling (SNMC),

    F(Φ) = 1 - 8πξΦ^2,

where ξ∈ℝ. Its motivation ranges from fundamental considerations—such as those arising from the quantization of the classical theory in a curved spacetime <cit.>—to its usefulness in cosmological model-building, specially in inflationary scenarios <cit.>. 

In the class of scalar-tensor theories (STTs) defined by Eq. (<ref>), post-Newtonian (PN) deviations from GR are proportional to (dF/dΦ)|_Φ_0, where Φ_0 is the scalar field value at the current cosmological epoch <cit.>. For the SNMC, (dF/dΦ)|_Φ_0∝Φ_0, and so the observed agreement <cit.> between solar system observations and GR's predictions implies that Φ_0 must be close to zero—but does not limit the viable range of ξ.
Interestingly, however, even if Φ_0 = 0, STTs may still differ considerably from GR in their predictions for neutron stars (NS), due to a nonperturbative, strong-field effect known as spontaneous scalarization <cit.>. This effect, which has long been known to happen for sufficiently negative values of the nonminimal coupling ξ, is characterized by the formation of a scalar cloud that modifies the star's equilibrium and perturbative properties <cit.>, and has dramatic implications, most notably for pulsar-timing observables <cit.>. 
Indeed, the inconsistency between pulsar-timing data and certain aspects of NS phenomenology in STTs, such as the existence of scalar dipole radiation, of scalar-field induced variations in the NS moment of inertia, and so on, can be used to rule out almost the entire range ξ≲ - 2.2 of field couplings allowing for spontaneous scalarization (in the case where V(Φ) = 0, which we will refer to as “massless” for simplicity).

More recently, it has been shown that positive values of the nonminimal coupling ξ can give rise to a similar spontaneous scalarization effect around sufficiently compact neutron stars, i.e., stars with G M/(R c^2) ≳ 0.26, where G is Newton's constant, c is the speed of light, and M and R denote the NS mass and radius <cit.>. Massless STTs with ξ>0 (which include the conformal coupling ξ = 1/6 as a particular case) are known to provide consistent cosmological scenarios <cit.>, but remain largely unconstrained by astrophysical observations. 

The purpose of the present work is to explore how the main post-Keplerian pulsar-timing observables—the Einstein time delay γ, the rate of periastron advance ω̇, and the rate of decay of the orbital period Ṗ_b—are modified around scalarized NSs in STTs, and investigate whether pulsar-timing data could also be used to constrain massless STTs with ξ >0. Interestingly, we find some crucial differences between the nature of spontaneous scalarization in the ξ <0 and ξ >0 cases, which reduces the effectiveness of pulsar timing observations in placing new constraints, even in the presence of spontaneous scalarization. In particular, the main scalar charge, α_A, entering these observables is typically suppressed in the ξ >0 case, and becomes progressively smaller as ξ increases. As a consequence, even though the presence of dipolar scalar radiation in STTs gives a contribution to Ṗ_b which is enhanced by a factor of (c/v)^2 (where v is the relative orbital velocity) with respect to GR, this term is suppressed by the smallness of the scalar charge, and Ṗ_b becomes dominated by the usual quadrupolar contribution. 
By exploring the dependence of α_A on the local scalar field environment, we also argue that the feedback mechanism responsible for the effect of dynamical scalarization found in some STTs <cit.> will likely be absent when ξ>0.

This work is organized as follows. In Sec. <ref> we discuss in more detail the framework we consider, including our choices of NS equations of state. In Sec. <ref> we discuss how the pulsar-timing observables {γ, ω̇, Ṗ_b} are modified in STTs and briefly review how to compute the scalar charges {α_A, β_A, k_A} that enter such observables. Our main results are presented in Sec. <ref> and final considerations are made in Sec. <ref>. A toy model presented in Appendix <ref> aims to elucidate in a simpler setting some of the features described in Sec. <ref>. We adopt natural units in which c = G = 1 unless specified.





§ FRAMEWORK
 



 §.§ Field equations


Neutron stars can be studied in STTs by adding the contribution from the stellar fluid to the gravitational action (<ref>):

    S = S_g + S_m [Ξ_m; g_μν].

We assume that the energy-momentum tensor of the matter fields Ξ_m, T^μν≡ (2/√(-g)) δ S_m/δ g_μν, has the form of a perfect fluid:

    T^μν = (ϵ + p) u^μ u^ν + p g^μν,

where p and ϵ denote the pressure and energy density in the fluid's rest frame and u^μ is the four-velocity of fluid elements. 

For definiteness, in this work we will focus on massless scalar fields with no self-coupling, setting V(Φ) = 0 in Eq. (<ref>). A mass term, for instance, would have the effect of delaying the onset of spontaneous scalarization to larger values of |ξ|  <cit.>, while other choices of V(Φ) could lead to a different phenomenology, including screening mechanisms <cit.> (see also Refs. <cit.> in the context of pulsar-timing).
An additional simplification to Eq. (<ref>) results from exploiting its invariance under a scalar field redefinition, Φ→φ(Φ), to set Z(Φ) to a constant. Then, the only free parameter is the nonminimal coupling constant ξ, once F(Φ) is chosen as in Eq. (<ref>).

For numerical calculations, it is often convenient to define the conformally rescaled (Einstein-frame) metric

    g_μν^* ≡ F(Φ) g_μν

and redefine the scalar field, Φ→φ(Φ), so that 

    ( dφ/dΦ)^2 = 3/4 F(Φ)^2( d F(Φ)/dΦ)^2 + Z(Φ)/2 F(Φ).

These transformations turn Eq. (<ref>) into 

    S     = 1/16 π∫ d^4 x √(-g^*)( R^* - 2 g^μν_* ∂_μφ∂_νφ) 
       + S_m[Ξ_m; g^*_μν/F(Φ(φ))],

and therefore re-express the coupling between scalar field and geometry as a coupling to matter.
The field equations in the Einstein frame take a simpler form,

    G_μν^* - 2 ∂_μφ∂_νφ + g^*_μν g^σρ_* ∂_σφ∂_ρφ  = 8π T_μν F(φ)^-1 ,
       ∇_*^μ∇^*_μφ = - 4 π F(φ)^-2α(φ) T,
 
where T ≡ g_μν T^μν = 3 p - ϵ is the trace of the energy-momentum tensor (<ref>),

    α(φ) ≡ - 1/2d ln F(Φ(φ))/dφ,

and all quantities marked with an asterisk are computed from the metric (<ref>). We emphasize that although we restore to practical computations in the Einstein frame, our results can be easily translated to the Jordan frame and our physical conclusions are frame-independent.

For F(Φ) given by Eq. (<ref>), and adopting a canonical normalization, Z(Φ) = 8π, the field redefinition (<ref>) gives

    dφ/dΦ = 2 √(π)√(1-8πξ (1-6ξ) Φ^2)/1-8πξΦ^2,

from which we can see that for ξ >0 the finite domain Φ∈ (-Φ_cr,Φ_cr), with Φ_cr≡ 1/√(8πξ), is mapped into φ∈ (-∞, ∞). No such restriction arises for ξ <0. If one interprets G_eff = G/F(Φ) as an effective gravitational coupling, gravity becomes weaker in the presence of the scalar field when ξ <0, and stronger if ξ >0, becoming infinitely attractive as Φ→±Φ_cr.



Although Eq. (<ref>) can be integrated in terms of elementary functions, the inverse transformation Φ(φ) must be obtained numerically. Therefore, we will find it useful in this work to consider not only F(Φ) as given by Eq. (<ref>) but also the following hyperbolic coupling (HC) model, already expressed in terms of the Einstein-frame field:

    F(φ) = [ cosh(2 √(3)ξφ) ]^-1/(3ξ),

so that

    α(φ) = 1/√(3)tanh( 2√(3)ξφ).

This coupling function was suggested in Ref. <cit.> (with β = 2 ξ) as a useful analytical approximation to a standard nonminimally coupled scalar field (see Fig. <ref>), and has been considered in Refs. <cit.> in the context of pulsar-timing observations. Note that the SNMC and HC agree up to a cubic term in an expansion around φ = 0: In both cases F(φ) = 1 - 2 ξφ^2 + O(φ^4).

As mentioned in the introduction, solar system observations constrain the scalar field value at the current cosmological epoch, which we denote as Φ_0 (with φ_0 ≡φ(Φ_0) in the Einstein frame), to be very small. For instance, the parametrized post-Newtonian parameter γ_PPN <cit.>, which takes the form

    1- γ_PPN = (dF /dΦ)^2/ZF + 2 (dF/dΦ)^2

for the theories described by Eq. (<ref>), 
is subject to the Cassini bound  |1 - γ_PPN| ≲ 2.1 × 10^-5 <cit.>. The bound translates into |Φ_0| ≲ 4.6 × 10^-4 /|ξ|, which becomes more stringent as |ξ| increases. In what follows we will typically fix the asymptotic, cosmological value of the scalar field to be Φ_0 = 0 (or φ_0 = 0), but will also discuss some features of the Φ_0 ≠ 0 case.







 §.§ Equation of state


In this work we will adopt three theoretical equations of state (EoS) for nuclear matter, the SLy <cit.>, ENG <cit.>, and MPA1 <cit.> models. In GR, the sequence of equilibrium configurations generated by these EoS is causal (in the sense that the speed of sound does not exceed the speed of light inside any stable star) and has a large enough maximum mass to accommodate the observation of a ∼ 2 M_⊙ NS <cit.>. Additionally, these EoS allow for stable NSs that are compact enough to trigger a spontaneous scalarization effect for both positive and negative values of the coupling constant ξ (see Sec. <ref> below). 

Instead of implementing these EoS through the interpolation of tabulated points, we shall approximate them by piecewise polytropes, adopting the parametrization developed in Ref. <cit.>. This parametrization was shown to accurately reproduce the main NS properties predicted by theoretical EoS <cit.>, as well as NS scalar charges in STTs <cit.>. 






 §.§ Spontaneous scalarization windows
 

One of the most interesting phenomenological aspects of STTs in astrophysical scenarios is the spontaneous scalarization effect <cit.>. 
In this section we briefly review the basic ideas behind this effect, and present the regions in parameter space where it takes place for the EoS employed in this work.

A careful inspection of the field equations (<ref>) and (<ref>) readily shows that a trivial scalar field profile, φ = 0, together with general-relativistic metric and matter configurations, form a solution of the field equations. However, for some relativistic stars, this trivial, GR-like solution may not be stable under scalar field perturbations <cit.>. This can be seen by expanding Eqs. (<ref>) and (<ref>) around φ = 0. To linear order in the perturbed quantities, the metric and fluid variables are not modified, while for the scalar field perturbation one obtains

    ^(0)δφ = - 8 πξ T^(0)δφ,

where the index (0) labels background quantities.
In the right-hand side of Eq. (<ref>), the combination m_eff^2 ≡ - 8πξ T^(0) can be loosely interpreted as an effective mass squared, and the fact that this can be negative signals the possible appearance of (tachyonic-like) instabilities (see Refs. <cit.> for a quantum analogue). 
The nonlinear development of this instability is the scalarization phenomenon: the spontaneous development of a cloud of scalar field around the star <cit.>. 



Since the trace of the energy-momentum tensor, T = 3p - ϵ, is typically negative (as energy density dominates over pressure), the scalarization effect would occur for ξ < 0 (so that m_eff^2 <0). Here, however, we will be mostly interested in studying pulsar-timing observables in the presence of positive values of the nonminimal coupling, since these are still unconstrained by astrophysical observations. In this case, in order to display the nontrivial phenomenology related to spontaneous scalarization, NSs must be sufficiently massive and compact that the trace of the energy-momentum tensor becomes positive in a region of their interior (so that m_eff^2 <0 in a sufficiently large region inside the star). In Ref. <cit.> the minimum compactness required for this property to hold was estimated to be C ≡ M/R = 0.262^+0.011_-0.017 (90% confidence interval), with the error bars accounting for our ignorance on the nuclear EoS. All the EoS considered in this work allow for stars more compact than this threshold value. 

In Fig. <ref> we show the regions of coupling parameters (ξ>0) and NS masses for which spontaneous scalarization takes place, for the three EoS considered in this work. To find these regions, we first solve the Tolman-Oppenheimer-Volkoff equations to construct a sequence of equilibrium configurations in GR and then determine whether unstable modes (δφ∝ e^Ω t, Ω>0) of Eq. (<ref>) can be found in that background. The lines delimiting the onset of instability coincide with the onset of spontaneous scalarization (as long as φ_0 =0) <cit.>. 
Notice that, although the realistic EoS considered in this work do not allow configurations that are compact enough to trigger this effect for the conformal coupling ξ = 1/6, more exotic structures could do so (see Fig. 3 of Ref. <cit.> for the case of a thin shell).




§ PULSAR-TIMING OBSERVABLES
 

Since the discovery of the first binary pulsar by Hulse and Taylor <cit.>, pulsar timing has become a major tool for testing GR (see <cit.> and references therein). The essence of pulsar timing lies in connecting the observed arrival time of the radio pulses to the proper time of emission. The resulting timing formula is obtained from a succession of steps, as follows <cit.>. First, by relating proper time in the pulsar's rest frame to a coordinate system attached to the binary center of mass (CM), one ends up with contributions coming from the pulsar motion around the CM (the transverse Doppler shift), as well as from a varying gravitational redshift that depends on the relative distance between the pulsar and its companion. These effects are combined in the “Einstein time delay”, usually expressed in terms of the measurable parameter γ. The next step consists in relating the coordinate time of emission to the coordinate time of arrival of a given pulse, by integrating along the null geodesics covered by the radiation. This picks up effects coming from dispersion in the interstellar medium, a geometrical contribution known as Roemer time delay, as well as from the Shapiro time delay due to the companion's gravitational well. Additional corrections are also accounted for, ranging from aberration effects due to the pulsar's rotation to corrections due to the motion of the Earth.  

At any instant, the orbit of each member of the binary system is tangent to a Keplerian ellipse (“osculating” orbit), characterized by six orbital parameters—say, the semilatus rectum p, eccentricity e, longitude of pericenter ω, time of pericenter passage T_0, inclination i, and angle of nodes Ω. In order to account for deviations from Newtonian dynamics, these parameters are allowed to undergo secular variations: x → x + ẋ t. In particular, the rate of periastron advance ω̇ and the rate of decay of the binary period Ṗ_b (derivable from the parameters above) are typically measurable.

In this work we will focus on the three classical pulsar-timing observables ω̇, γ, and Ṗ_b, although many other post-Keplerian (PK) parameters can in principle be inferred from pulsar-timing observations <cit.>. A given theory of gravity will predict the value of these observables as a function of the Keplerian parameters and the masses of the binary components. If the masses (m_p and m_c) are unknown, as often is the case, one can use the measurement of two of these PK parameters to infer m_p and m_c, and perform a test of the gravitational theory with the third. This is typically portrayed by drawing, for each PK parameter, a level surface in the m_p-m_c diagram corresponding to the measured value of that parameter. The theory is consistent with observations if the resulting curves all intersect at the same point, the binary component masses <cit.>. 

In STTs, as in GR, the theoretical prediction for pulsar-timing observables is based on the PN description of the orbital motion. Explicitly, one has <cit.>:

    ω̇   = 3 n_b/1-e^2𝗏_b^2/c^2[
    1-α_p α_c/3/1+ α_p α_c - m_p β_c α_p^2 + m_c β_p α_c^2/6 M(1+α_p α_c)^2], 
    γ   = e/n_bm_c/M(1+α_c α_p)𝗏_b^2/c^2[ 1+ α_c k_p + (1+α_p α_c) m_c/M],
    Ṗ_b    = Ṗ_φ^monopole + Ṗ_φ^dipole + Ṗ_φ^quadrupole + Ṗ_g*^quadrupole,

where n_b≡ 2π/P_b, M ≡ m_p + m_c, and 𝗏_b ≡ (𝒢 M n_b)^1/3, with 𝒢≡ G (1+α_pα_c). The quantities α_A, β_A, and k_A, with the label A∈{p,c} denoting the pulsar or its companion, are functions of the (Einstein-frame) stellar mass m_A, and will be defined in Sec. <ref> below. Equation (<ref>) includes contributions related to monopole, dipole, and quadrupole scalar radiation, as well as the quadrupolar contribution from tensor waves familiar from GR. All terms are proportional to (𝗏_b/c)^5, and therefore of 2.5PN order, except from Ṗ_φ^dipole, that contributes already at 1.5PN. Since this is the dominant contribution to the energy loss and will be important for our discussion, we write this term explicitly:

    Ṗ_φ^dipole = -2 π m_p m_c/M^2(1+α_p α_c)𝗏_b^3/c^31 + e^2/2/(1 - e^2)^5/2 (α_p - α_c)^2.

For comparison,

    Ṗ_g*^quad = -192 π m_p m_c/5 M^2 (1+α_p α_c)(𝗏_b/c)^5/(1-e^2)^7/2( 1+ 73/24e^2 + 37/96e^4).

The expressions for the other terms can be found, e.g., in Refs. <cit.>. 



 §.§ Scalar charges


Crucially, in the pulsar timing observables (<ref>)-(<ref>) we find the appearance of the “gravitational form factors” (in the terminology of Ref. <cit.>) or “scalar charges” (in the terminology of Ref. <cit.>) α_A = α_A(m_A), β_A = β_A(m_A), and k_A = k_A(m_A) (with A ∈{p,c}). In this section we briefly review their definition, the rationale behind their appearance in pulsar-timing observables, and how to compute them in practice.

It is well known that, in GR, a body's orbital dynamics is governed by integral quantities such as its mass and spin, with corrections due to shape and internal structure appearing only at high PN orders.

For this reason, complex bodies can often be well approximated by point masses (insofar as their orbital motion is concerned). In STTs, the value of the scalar field around the star determines the strength of the effective gravitational coupling, and therefore influences its global properties. This feature can be effectively incorporated into a PN description based on point masses by allowing the mass of each body to be field dependent: For an N-body system, the action is taken to be <cit.>

    S_m = - ∑_A=1^N ∫ m_A (φ_A) dτ^*_A.

It is worthwhile to emphasize the dual role played by the function m_A(φ_A) above. If one assumes that the inter-body distance D between the binary components is much larger than their typical size R, and allow the ratio R/D to shrink to zero, then m_A(φ_A) expresses the stellar mass in terms of the scalar field value at the star's location. This is the “outer”, PN perspective. However, the  function m_A(φ_A) is not determined self-consistently within the PN scheme, but is assumed to be known from the matching to the “inner” problem, where the stellar structure is determined. From this “inner” perspective, the point particle limit corresponds to the matching sphere becoming infinitely large, and from this point of view the function m_A(φ_A) denotes the mass the star has when the asymptotic value of the scalar field is φ_A.

Once the field is expanded around its (cosmological) asymptotic value in the PN approximation, the effect of a field-dependent mass is encoded in the asymptotic value of its derivatives. In particular, one defines:

    α_A ≡. d log m_A/dφ_A|_φ_0,


    β_A ≡. dα_A/dφ_A|_φ_0,

and so on. 
The scalar charges above (closely related to the “sensitivities” used in the Jordan-frame description <cit.>) are the ones appearing at the Newtonian and post-Newtonian levels, relevant for the derivation of the pulsar timing observables (<ref>)-(<ref>).

The additional scalar charge k_A entering the Einstein time delay (<ref>) is defined as 

    k_A ≡ - .  ∂log I_A/∂φ_A|_φ_0,

where I_A is the moment of inertia of star A. It appears in the computation of the parameter γ when relating the intrinsic time of the pulsar clock to the proper time in a local inertial frame around the pulsar, due to the fact that, in STTs, the moment of inertia of the star depends on the local scalar field environment—which may fluctuate in a binary system, causing the angular velocity to fluctuate as well <cit.>. 

In order to explain how the scalar charges α_A, β_A, and k_A are computed, we now revert attention to the inner problem of an isolated, slowly rotating star, which is a suitable approximation for most of the observed pulsars. We therefore consider the spacetime of a slowly rotating body, with metric <cit.>

    ds_*^2    = - e^ν (r) dt^2 + e^λ(r) dr^2 + r^2 dθ^2 + r^2 sin^2 θ dϕ^2 
        - 2 ϖ(r) r^2 sin^2θ dt dϕ,

where ϖ(r) accounts for frame-dragging effects and is the only O(Ω) correction to the static case, with Ω denoting the star's angular velocity (see Refs. <cit.> for generalizations). 

The differential equations governing the metric functions ν(r), λ(r), ϖ(r), scalar field φ(r), and fluid variable p(r) can be found, e.g., in Refs. <cit.>, and are reproduced in Appendix <ref> for the sake of completeness. By imposing regularity conditions at r=0, asymptotic flatness at spatial infinity, as well as the condition φ(r) r →∞⟶φ_A, one can integrate these equations to obtain a one-parameter family of solutions (labeled, say, by the central value of the pressure). 

The scalar charges in Eqs. (<ref>)-(<ref>) measure the change of the stellar properties as φ_A changes, while keeping fixed the baryon mass inside the star (m̅_A = ∫ρ√(-g) u^0 d^3x, where ρ is the rest-mass density).

Global properties of the NS can be extracted from the asymptotic behavior of the metric functions. The total mass m_A is obtained from the 1/r behavior of e^ν(r)r→∞⟶ 1 - 2 m_A/r + O(1/r^2). The total angular momentum J_A is obtained from the 1/r^3 behavior of ϖ(r) r→∞⟶ 2J_A/r^3 + O(1/r^4), and enables the computation of the star's moment of inertia: I_A = J_A/Ω. Finally, from the leading 1/r contribution to the scalar field at spatial infinity, φr→∞⟶φ_A + ω_A/r + O(1/r^2), one can directly obtain the scalar charge α_A through α_A = -ω_A / m_A. That the quantity computed in this way is equivalent to the definition in Eq. (<ref>) is shown in the Appendix A of Ref. <cit.>. 



Contrary to the computation of α_A, which can be done straightforwardly, computing β_A and k_A is more involved, and here we proceed as follows. 
First, we construct three sequences of equilibrium solutions, for φ_A^(0) = φ_0, φ_A^(+) = φ_0 + Δφ and φ_A^(-) = φ_0 - Δφ, storing, for each value of the central pressure, the total mass m_A, baryonic mass m̅_A, and moment of inertia I_A of the resulting solution. This data is then used to define, by interpolation, functions m^(i)_A(m̅_A), α^(i)_A(m̅_A), and I^(i)_A(m̅_A), with i ∈{0,+,-}, from which we can compute the scalar charges by a finite difference approximation to the derivative operators in Eqs. (<ref>) and (<ref>). For our purposes, we find it enough to use a central stencil for the derivative operator, setting, for instance,

    k_A(m̅_A) ≈1/I^(0)_A(m̅_A)I^(+)(m̅_A) - I^(-)(m̅_A)/2Δφ.

More details on the numerical procedure and the accompanying errors can be found in Appendix <ref>. Note also that a bank of scalar charges, as well as a thorough account of the procedure to compute them, was recently provided in  Ref. <cit.>. The authors did not, however, explore the possibility of spontaneous scalarization for ξ>0, which is our main focus here.





§ RESULTS
 



 §.§ Scalar charges
 

In order to establish a point of comparison between the cases where ξ is negative and positive, in Fig. <ref> we show the behavior of the total mass m_A, moment of inertia I_A, and scalar charge α_A for a sequence of equilibrium solutions with ξ = -3 and -5 and ξ = 25, 50, and 100. The asymptotic value of the scalar field is taken to be φ_0 = 0.
As anticipated in Sec. <ref> (see Fig. <ref>), scalarized solutions for ξ >0 only exist for the most massive and compact stars, while scalarization happens for ξ <0 in a much wider range of masses. From the first row of Fig. <ref>, we see that as |ξ| increases the maximum mass of scalarized solutions increases when ξ <0 and decreases when ξ >0. This is consistent with the interpretation of G_eff = G/F(Φ) as an effective gravitational coupling. As ξ>0 increases, so does G_eff: gravity becomes stronger in the presence of the scalar field and less massive scalarized stars can be supported without undergoing gravitational collapse. The opposite happens for ξ <0.
For a given mass, the second row of Fig. <ref> shows that the moment of inertia is larger (smaller) than in GR when ξ < 0 (>0). This can also be understood as the effect of gravity becoming weaker (stronger) in each of these cases, and the stellar size becoming typically larger (smaller) than in GR when a nontrivial scalar field profile is present.

Since the field equations (<ref>) and (<ref>) are invariant under the transformation φ→ - φ for the coupling functions we consider, one always finds two twin scalarized solutions related by the aforementioned transformation, and with opposite scalar charges α_A. This is true as long as φ_0 = 0, as in Fig. <ref>, otherwise the boundary condition breaks the reflection symmetry of the solution.
In the third row of Fig. <ref> we show the scalar charge α_A for solutions with a positive scalar field profile. 

Two notable differences are found in the ξ>0 case (right panel), compared to when ξ<0 (left panel), namely that the magnitude of α_A (i) is now much smaller and (ii) typically decreases as |ξ| increases.
Both have to do with the fact that the scalar field tends to be amplified in the stellar region where ξ T > 0 and suppressed in the region where ξ T <0. For any realistic EoS, the trace of the energy momentum can be positive only in a small region in the stellar interior; therefore, a scalar field with ξ >0 is amplified in this inner region, but is necessarily suppressed in the outer layers of the star. As a consequence, although the central value of the scalar field may increase with increasing ξ, the scalar charge, measured asymptotically, ends up being smaller. In Appendix <ref> we present a simple, analytically solvable, toy model that illustrates these points. 

Figure <ref> shows the scalar charges α_A, β_A, and k_A for ξ = 25, 50, and 100, for both the SNMC and the HC models, and for three nuclear EoS describing the NS fluid. From the first row we see that the properties of α_A described in the preceding paragraph are not altered by the nuclear EoS. Indeed, the magnitude of all scalar charges is only mildly influenced by the EoS, the main effect of which is to change the range of masses where scalarization takes place. 

The second row of Fig. <ref> shows the scalar charge β_A = (dα_A/dφ_A)_φ_0 as a function of the stellar mass. This is a positive, order-of-unity quantity, which increases with increasing ξ. The fact that β_A is positive (and α_A negative) means that α_A receives a positive increment when the asymptotic value of the scalar field φ_A increases, and therefore typically diminishes in absolute value (as long as φ_A is close to 0). This is opposite to what happens when ξ <0, in which case β_A <0 and α_A increases in absolute value with increasing φ_A. This is shown more explicitly in Fig. <ref>, where the scalar charge α_A is represented for φ_A = 10^-4 and 10^-3, for ξ = 50 and ξ = -3. 



Interestingly, the features described above can have implications for scalarization in dynamical situations.
“Dynamical scalarization” was first observed in binary NS simulations in STTs <cit.>, wherein NSs that were not compact enough to scalarize in isolation suddenly developed a large scalar charge once their separation became sufficiently small, making the coalescence proceed in a faster timescale. It was then understood to be due to a kind of feedback mechanism <cit.>. 
Indeed, as discussed before, from the PN perspective the mass and other stellar properties are functions of the local scalar field value at the star's worldline, which is influenced by the presence of the companion star. In STTs with ξ <0, if the ambient scalar field value grows, it induces a growth in ω_A = -m_A α_A (see Fig. <ref>). This, in turn, increases the local value of the scalar field at the companion's location (given by φ_B ∼φ_0 + ω_A/r, where r is the separation distance), and the positive feedback proceeds until a fixed-point is reached. Although exploring this in detail is beyond the scope of the present work, we can anticipate that the opposite behavior is likely to occur in STTs with ξ >0: In this case β_A >0 and ω_A decreases with increasing φ_A; therefore the feedback mechanism present in the ξ <0 case will most likely be absent or reversed.

Finally, let us come back to the last row of Fig. <ref>, where k_A [cf. Eq. (<ref>)] is shown as a function of the stellar mass. The typical values of k_A are seen to decrease with increasing ξ, with a large spike close to the maximum allowed mass. This is a somewhat similar behavior to the ξ < 0 case (see, e.g., Fig. 4 in Ref. <cit.>), although the typical values in Fig. <ref> are quite small due to the relatively large values of ξ we consider. In all cases, the HC reproduces well the qualitative features of the SNMC.





 §.§ Implications for pulsar-timing observables
 

Having understood the properties of the scalar charges α_A(m_A), β_A(m_A), and k_A(m_A), we now turn to the pulsar timing observables ω̇, γ, and Ṗ_b, given in STTs as in Eqs. (<ref>)-(<ref>). 

Pulsar-timing observations can be used to put stringent bounds on ξ < 0 <cit.>. 
A crucial feature that makes this possible is that the scalar charge α_A, which enters in all pulsar-timing observables, increases in magnitude as ξ becomes more negative, making deviations from GR stronger. As a consequence, one is able to exclude the entire range of negative couplings up to a certain value. 

The situation for ξ>0 is much more permissive. On the one hand, since the scalar charge α_A typically decreases in magnitude as ξ>0 increases (cf. Fig. <ref>; see also Fig. 8 of Ref. <cit.>), deviations of pulsar-timing observables from GR are suppressed for large values of the nonminimal coupling. On the other hand, as we decrease the value of ξ in search of larger scalar charges, the range of NS masses allowing for spontaneous scalarization gets progressively smaller (cf. Fig. <ref>). Thus, STTs with ξ>0 tend to pass pulsar-timing tests, with the possible exception of a small range of couplings and NS masses. 

Having described the global picture, let us give a few more details about each of the observables ω̇, γ, and Ṗ_b. 
Since β_A = O(1) for ξ > 0 (cf. Fig. <ref>), we can see from Eq. (<ref>) that corrections to ω̇ are O(α_A^2), and therefore negligible due to the smallness of the scalar charge α_A in this case.
The same is true for the Einstein time delay γ, except for the term proportional to α_c k_p in Eq. (<ref>). Since k_A seems to diverge in limit of the maximum NS mass (cf. Fig. <ref>), this term still could give a non-negligible contribution as long as the companion's charge α_c is not vanishingly small. This would be restricted, however, to an exceptionally thin range of NS masses. As for the rate of decay of the orbital period, Ṗ_b, it is enough to compare the contribution in Eq. (<ref>) coming from dipolar scalar radiation to the usual quadrupolar contribution due to tensor waves, Eq. (<ref>). Although the first one dominates the second by a factor of (c/𝗏_b)^2, it is also suppressed by a factor of (α_p - α_c)^2. Typical orbital velocities in binary systems observed through pulsar-timing techniques are of the order of (𝗏_b/c) ∼ 10^-3. Interestingly, this is of the same order of magnitude as the largest values of |α_A| shown in Fig. <ref>. Therefore, it might still be possible to probe STTs with ξ≲ 25 with the radiation emitted by the most massive pulsars. However, as emphasized above, the range of masses allowing for spontaneous scalarization gets narrower as ξ decreases, and such effects, if present, would be restricted to very special systems. For the largest ranges of NS masses and nonminimal couplings, STTs with ξ>0 would still evade the sharp knife of pulsar-timing tests.





§ CONCLUSION
 

It is well known that general relativity has passed with flying colors all strong-field tests imposed by pulsar-timing observations. Moreover, these observations have stripped many modified theories of gravity of a large portion of their parameter space, confining their predictions to the close vicinity of GR <cit.>. This is the case, in particular, of a class of scalar-tensor theories of gravity including the case of a massless scalar field with the standard, ξ R Φ^2, nonminimal coupling to gravity. In these theories a tachyonic-like instability develops for sufficiently compact stars and ξ≲ - 2.2, leading to the development of a scalar cloud around the star (see Sec. <ref>). If the scalar field is thus activated around a member of a binary system, it drives the emission of dipole scalar radiation, contributing to a steeper decrease of the orbital period in time. Since this effect—as well as the changes introduced in other observable quantities (see Sec. <ref>)—is incompatible with pulsar-timing data, almost the entire range of couplings allowing for spontaneous scalarization has now been excluded <cit.>.

Recently, it has been advocated that a similar scalarization effect might occur for positive values of the nonminimal coupling ξ, around the most massive neutron stars found in Nature, as long as they are also sufficiently compact <cit.> (which depends on the still unknown nuclear equation of state). This seems to open the possibility of using measurements of the most massive observed NS, such as the pulsar PSR J0348+0432 <cit.> or possibly PSR B1957+20 <cit.>, in order to also probe this range of couplings. 

In this work we carried out a study of the scalar charges α_A, β_A, and k_A, which determine the behavior of the main post-Keplerian pulsar timing observables, ω̇, γ, and Ṗ_b, in STTs with ξ >0 (see Sec. <ref>). We find that the scalar charges differ remarkably in the ξ >0 and ξ<0 cases. In particular, the main scalar charge α_A governing Newtonian and post-Newtonian deviations from GR is suppressed as ξ>0 becomes large, while the range of masses allowing for spontaneous scalarization decreases as ξ>0 becomes smaller. As discussed in detail in Sec. <ref>, this indicates that STTs with ξ >0, even in the presence of spontaneous scalarization, are able to pass pulsar-timing tests, with the possible exception of an exceedingly narrow range of couplings and NS masses. Moreover, the contrasting properties exhibited by the scalar charge β_A in the ξ>0 case suggest that the feedback mechanism responsible for the effect of dynamical scalarization might not be present for positive values of the nonminimal coupling. It remains an interesting research avenue to determine which NS properties are most sensitive to a scalar field background in this regime of small scalar charges (e.g. the NS oscillation spectrum as recently suggested by Ref. <cit.>).






We are glad to thank Gabriel Vidal for his involvement in the early phases of this project, and Néstor Ortiz for valuable comments on the manuscript.
R.M. acknowledges partial financial support by the National Council for Scientific and Technological Development–CNPq. T.O. acknowledges financial support from funding agencies CAPES and FAPERJ. 





§ A SIMPLE TOY MODEL
 

In order to gain some intuition about the behavior of the scalar charge α_A, we consider a simpler version of Eq. (<ref>) <cit.>, where we neglect the metric curvature and the details of the coupling function, keeping only its linear piece:

    Δφ = - ϵκ^2 φ.

Here Δ≡ d^2/dr^2 + (2/r) d/dr, κ≡√(8 π |ξ T|), and ϵ = sign (ξ T). Moreover we will assume that κ is constant throughout the star, but that ϵ can change sign, according to

    ϵ = 
    
    + 1,      0 ≤ r ≤ r^+ 
    
    - 1,      r^+ < r ≤ R

where R is the stellar radius and 0 ≤ r^+ ≤ R. The case most often considered in the literature, i.e., ξ <0 and T<0 throughout the star, would correspond to r^+ = R. On the other hand, we are most interested here in the case where ξ >0 and the trace T changes sign inside the star, as captured by the general form (<ref>).
Taking κ to be constant, as we do for simplicity, implies the additional assumption that the typical strength (say, the mean value) of T>0 in the region 0 ≤ r ≤ r^+ is comparable with the typical value of T<0 in the region r^+ < r ≤ R.

The general solution to this problem is

    φ =
    φ_c sin(κ r)/κ r,      0 ≤ r ≤ r^+ , 
    
    A sinh(κ r)/κ r + B cosh (κ r)/κ r,      r^+ < r ≤ R , 
    φ_0 + ω/r,      r > R ,

where φ_c, A, B, and ω can be explicitly computed in terms of φ_0, κ, r^+, and R by imposing continuity of the field and its derivative across r=r^+ and r=R. In particular, 

    φ_c = φ_0/sin(κ r^+) sinh(κ r^-) + cos(κ r^+) cosh(κ r^-)

and

    ω = - φ_0 R + φ_0/κtan(κ r^+) + tanh(κ r^-)/1 + tan(κ r^+) tanh(κ r^-),

where we defined r^- ≡ R - r^+. It is worth writing explicitly two particular cases: If r^+ = 0, we have

    φ_c = φ_0/cosh(κ R),   ω = - φ_0 R ( 1 - tanh(κ R)/κ R),

while if r^- = 0 we get

    φ_c = φ_0/cos(κ R),   ω = - φ_0 R ( 1 - tan(κ R)/κ R).

In the first case, Eq. (<ref>), the field is suppressed with respect to its asymptotic value and |ω| is bounded by |φ_0| R, going to zero as φ_0 → 0. This would be the picture if, say, ξ >0 and T<0 throughout the star. In the second case, Eq. (<ref>), the field is amplified with respect to its asymptotic value and |ω| is enhanced with respect to |φ_0| R. Indeed, these quantities may have a nonzero limit even when φ_0 → 0, as long as κ R = π/2. 
This gives a heuristic picture of the spontaneous scalarization effect when ξ < 0 <cit.>.

Take now the full expressions (<ref>) and (<ref>), which are the relevant ones when ξ >0. We see that there is a competition between enhancement and suppression effects. In particular, the onset of scalarization, which was determined by the condition (κ R) = 0 for Eq. (<ref>), is delayed to higher values of κ R, with the relevant condition becoming (κ r^+) = - tanh (κ r^-). 
Moreover, a number of features displayed by the scalar charge α_A for scalarized solutions (cf. Sec. <ref>) are already exhibited in the expressions above (outside the scalarization regime). For instance, the 1/κ dependency in Eq. (<ref>) indicates that |ω| decreases with increasing κ, except near the resonances mentioned above. This is in contrast with the case (<ref>), where |ω| increases with κ.





§ NUMERICAL SETUP
 

With the metric ansatz (<ref>), the field equations (<ref>) and (<ref>) yield

    m'    = 4 π r^2 F^-2ϵ + 1/2 r(r-2m) (φ')^2, 
    ν'     = 8π r^2 F^-2 p/r-2m + r(φ')^2 + 2m/r(r-2m), 
    φ”   = 4π r F^-2/r-2m[ α (ϵ - 3 p) + r φ' (ϵ - p) ] - 2(r-m)/r(r-2m)φ',  
    
    p'    = - (ϵ + p) [ 4π r^2 F^-2 p/r-2m + r(φ')^2/2 + m/r(r-2m) + αφ' ],  
    ϖ”   = ϖ' [ -4/r + r (φ')^2 + 4π r^2 (ϵ + p)/F^2 (r-2m)] 
        + 16π r (ϵ + p)/F^2 (r-2m) (Ω - ϖ),

where primes denote radial derivatives, Ω≡ u^ϕ/u^t is the fluid angular velocity, and we have set e^λ(r)≡ 1/(1-2m(r)/r). When supplemented by a choice of EoS, the set of differential equations above can be integrated numerically by standard methods. The relevant boundary conditions to the problem are regularity at r=0, which requires that m(0) = 0, ϖ'(0) = 0, and φ'(0) = 0, and asymptotic flatness, which requires that lim_r→∞ν(r) = 0, lim_r→∞ϖ(r) = 0 and lim_r→∞φ(r) = φ_A. The condition p(R) = 0 determines the stellar radius R. As described in Sec <ref> the stellar mass m_A and angular momentum J_A are determined by the asymptotic behavior of the metric functions. Alternatively, Eqs. (<ref>)-(<ref>) can be integrated up to the stellar radius and matched to the known form of the exterior solution <cit.>, from which global quantities like m_A and J_A can be inferred.

As described in Sec. <ref>, the computation of the scalar charges, 

    α_A = d log m_A/dφ_A,   β_A = d α_A/dφ_0,   
    k_A = - d log I_A/dφ_0,

is somewhat involved, since the derivatives must be evaluated for a fixed value of the baryonic mass. For this purpose, given an EoS and a value for ξ, we construct three sequences of equilibrium solutions, with the asymptotic value of the scalar field given by φ_A^0 = φ_0, φ_A^+ =φ_0 + Δφ, and φ_A^- =φ_0 - Δφ. For most of the results presented in this work we employ φ_0 = 0 and Δφ = 0.0005, except for Fig. <ref> where different values of φ_0 are used. From the data thus generated, we construct by interpolation the functions m_A^(i)(m̅_A), α_A^(i) (m̅_A), and I_A^(i)(m̅_A), with i ∈{0,+,-} and estimate the scalar charges (<ref>) through a simple central finite difference approximation to the derivative operator:

    α_A^approx (m̅_A)    = 1/m_A^(0)(m̅_A)m_A^(+)(m̅_A)-m_A^(-)(m̅_A)/2Δφ,
    
    β_A^approx  (m̅_A)    = α_A^(+)(m̅_A)-α_A^(-)(m̅_A)/2Δφ, 
    
    k_A^approx(m̅_A)     = 1/I^(0)_A(m̅_A)I^(+)(m̅_A) - I^(-)(m̅_A)/2Δφ.


Recall that the scalar charge α_A can either be computed from the expression above, or directly from the asymptotic behavior of the scalar field, through α_A = -ω_A/m_A. The latter procedure yields a much more reliable estimate for α_A, and so the comparison between α_A^(0) and α_A^approx enables us to estimate the error incurred in the finite difference approximations above. This is shown in Fig. <ref>. For Δφ = 5 × 10^-4, which is typically used in this work, we see that the relative error in the scalar charge is of the order of 0.01%, but may increase near the boundaries of the mass interval where scalarized solutions exist. Errors of the same order are obtained for other values of ξ as well.



Notice that, when presenting our results, we restrict consideration to the range of baryonic masses for which both (+) and (-) quantities exist. This may exclude from our consideration a narrow range of masses at the borders of the mass interval. These boundary values could be taken care of by a one-sided approximation to the derivative operator, which would lead, in particular, to a better resolution of the spikes in Fig. <ref>. These fine details are not, however, too relevant for our analysis.





