Scalar Charges and Scaling Relations in Massless Scalar-Tensor Theories

TL;DR

This paper calculates scalar charges for neutron stars in massless scalar-tensor theories, providing tabulated data and scaling relations to facilitate tests of these theories with pulsar timing and gravitational wave observations.

Contribution

It introduces a comprehensive calculation and tabulation of scalar charges in scalar-tensor theories, along with analytic scaling relations for quick predictions across parameter space.

Findings

Scalar charges depend on neutron star equations of state.

Scaling relations enable rapid estimation of scalar charges.

Results are compatible with constraints from gravitational wave observations.

Abstract

The timing of binary pulsars allows us to place some of the tightest constraints on modified theories of gravity. Perhaps some of the most interesting and well-motivated extensions to General Relativity are scalar-tensor theories, in which gravity is mediated by the metric tensor and a scalar field. These theories predict large deviations from General Relativity in the presence of neutron stars through a phenomenon known as scalarization. Neutron stars in scalar-tensor theories develop scalar charges, which directly enter the timing model for binary pulsars. In this paper, we calculate and tabulate these scalar charges in two popular, massless scalar tensor theories for a collection of neutron star equations of state that are compatible with constraints placed by the recent, gravitational wave observations of a binary neutron star coalescence. We then study these scalar charges and explore analytic scaling relations that allow us to predict their value in a large region of parameter space. Our results allow for the quick evaluation of the scalar charge in a large region of scalar-tensor theory parameter space, which has applications for gravitational wave tests of scalar-tensor theories, as well as binary pulsar experiments.

Tables

Table 1. Table 1: The finite difference coefficients found in Eqs. ( LABEL:eq:fd_central )-( LABEL:eq:fd_backward ). For every method (central, forward, or backward) the columns denote the value of the coefficient for the corresponding subscripts found in Eqs. ( LABEL:eq:fd_central )-( LABEL:eq:fd_backward ) .

Method	-4	-3	-2	-1	0	1	2	3	4
central	—	—	1/12	-2/3	0	2/3	-1/12	—	—
forward	—	—	—	—	-25/4	4	-3	4/3	-1/4
backward	1/4	-4/3	3	-4	25/4	—	—	—	—

Table 2. Table 2: An example of the layout for the first 9 columns of the data file for MO theory and AP3 EOS. For a given value of β 0 subscript 𝛽 0 \beta_{0} we tabulated the data for every value of α 0 subscript 𝛼 0 \alpha_{0} , of which for every value of α 0 subscript 𝛼 0 \alpha_{0} we tabulate data for all values of m A subscript 𝑚 𝐴 m_{A} using the grid described in the text. This pattern repeats for all values of − 5 ≤ β 0 ≤ + 5 5 subscript 𝛽 0 5 -5\leq\beta_{0}\leq+5 that we sampled. The definitions of the quantities in the columns is described in the text. The actual data files we provide also have data for the central density, radius, central scalar field, and surface scalar field of the NSs.

${\alpha }_A$	${\beta }_A$	$k_A$	${\overline{m}}_A$	$m_A$	${\alpha }_0$	${\log }_{10}{\alpha }_0$	${\beta }_0$
0.444979	0.270288	0.931259	1.16903	1.000	1.	0.	-5.00
0.444787	0.270318	0.931345	1.17155	1.002	1.	0.	-5.00
0.444595	0.270350	0.931432	1.17408	1.004	1.	0.	-5.00
⋮	⋮	⋮	⋮	⋮	1.	0.	-5.00
0.331799	0.410289	1.101105	2.67262	2.100	1.	0.	-5.00
0.440000	0.314515	0.925646	1.15792	1.000	0.9	-0.045757	-5.00
0.439811	0.314451	0.925772	1.16041	1.002	0.9	-0.045757	-5.00
⋮	⋮	⋮	⋮	⋮	⋮	⋮	-5.00
0.120843	10.909670	4.323983	2.50472	2.100	0.000002	-5.69897	-5.00
0.444997	0.270390	0.931196	1.16931	1.000	1.	0.	-4.98
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
0.0000004	1.062219	-0.000001	2.50430	2.100	0.000002	-5.69897	5.00

Table 3. Table 3: Numbered points appearing in Fig.  1 and there associated relative errors as determined by Eq. ( 63 ). The columns with ρ ~ c subscript ~ 𝜌 𝑐 \tilde{\rho}_{c} in parentheses are the errors determined in Sec.  A.1 . Recall that Eq. ( 63 ) actually yields the percent error in the solutions, therefore 10 0 superscript 10 0 10^{0} actually represents a 1% error, which is the highest error we find from our results.

#	${\beta }_0$	${\alpha }_0$	${\alpha }_A\left({\tilde{\rho }}_c\right)$	${\beta }_A\left({\tilde{\rho }}_c\right)$	$k_A\left({\tilde{\rho }}_c\right)$
1	-6	$1.44\times {10}^{-2}$	—	—	—
2	-5	${10}^{-5}$	${10}^{-2}$	${10}^0$	${10}^0$
3	-5	${10}^{-4}$	${10}^{-2}$	${10}^0$	${10}^{-1}$
4	-5	${10}^{-3}$	${10}^{-3}$	${10}^{-2}$	${10}^{-2}$
5	-5	${10}^{-2}$	${10}^{-3}$	${10}^{-3}$	${10}^{-3}$
6	-5	${10}^{-1}$	${10}^{-3}$	${10}^{-1}$	${10}^{-1}$
7	-5	${10}^{-0}$	${10}^{-3}$	${10}^{-1}$	${10}^{-1}$
8	-3	${10}^{-3}$	${10}^{-6}$	${10}^{-5}$	${10}^{-2}$
9	-3	${10}^{-1}$	${10}^{-5}$	${10}^{-2}$	${10}^{-2}$
10	3	${10}^{-3}$	${10}^{-6}$	${10}^{-5}$	${10}^{-1}$
11	3	${10}^{-1}$	${10}^{-5}$	${10}^{-2}$	${10}^{-1}$
12	-4	${10}^{-4}$	${10}^{-5}$	${10}^{-4}$	${10}^0$
13	-4	${10}^{-3}$	${10}^{-5}$	${10}^{-4}$	${10}^{-2}$
14	-4	${10}^{-2}$	${10}^{-5}$	${10}^{-5}$	${10}^{-3}$

Equations

$${\tilde{S}}_g=\int \frac{d^{4}x}{c}\frac{\sqrt{-\tilde{g}}}{4\kappa }\left[\varphi \tilde{R}-\frac{\omega (\varphi )}{\varphi }{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi \right],$$

$$S=\int \frac{d^{4}x}{c}\frac{\sqrt{-g}}{4\kappa }\left[R-2g^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi \right]+S_{\text{mbox}mat}\left[\chi ,A^2(\phi )g_{\mu \nu }\right],$$

$$A(\phi )={\varphi }^{-1},$$

$$\alpha (\phi {)}^2={\left(\frac{d\ln A(\phi )}{d\phi }\right)}^2=\frac{1}{3+2\omega (\varphi )}.$$

$$R_{\mu \nu }$$

$$\square \phi$$

$$T_{\mu \nu }^{\text{mbox}mat}\equiv \frac{2c}{\sqrt{-g}}\left(\frac{\delta S_m}{\delta g^{\mu \nu }}\right),$$

$$A(\phi )$$

$$V_{\alpha }(\phi )$$

$$\alpha (\phi )$$

$$\beta (\phi )$$

$$A(\phi )$$

$$V_{\alpha }(\phi )$$

$$\alpha (\phi )$$

$$\beta (\phi )$$

$$G_N=G\left[A_{\infty }^2\left(1+{\alpha }_{\infty }^2\right)\right].$$

$$\left|1-{\gamma }_{\text{mbox}PPN}\right|=\frac{2{\alpha }_{\infty }^2}{1+{\alpha }_{\infty }^2},$$

$$\left|1-{\beta }_{\text{mbox}PPN}\right|=\frac{{\beta }_{\infty }{\alpha }_{\infty }^2}{2{\left(1+{\alpha }_{\infty }^2\right)}^2},$$

$$\alpha (\phi )={\alpha }_0+{\beta }_0\phi ,$$

$${\alpha }_{\infty }^{\text{mbox}DEF}$$

$${\beta }_{\infty }^{\text{mbox}DEF}$$

$${\alpha }_{\infty }^{\text{mbox}MO}$$

$${\beta }_{\infty }^{\text{mbox}MO}$$

$${\omega }_A=-\frac{\partial m_A}{\partial {\phi }_{\infty }},$$

$${\alpha }_A={-\frac{{\omega }_A}{m_A}=\frac{\partial \ln m_A}{\partial {\phi }_{\infty }}|}_{{\bar{m}}_A},$$

$${\beta }_A={\frac{\partial {\alpha }_A}{\partial {\phi }_{\infty }}|}_{{\bar{m}}_A},$$

$$k_A={\frac{\partial \ln I_A}{\partial {\phi }_{\infty }}|}_{{\bar{m}}_A},$$

$$d{s}^2$$

$$e^{\lambda (\rho )}={\left(1-\frac{2\mu (\rho )}{\rho }\right)}^{-1}.$$

$$T_{\mu \nu }=(ϵ+p)u_{\mu }{u}_{\nu }+p{g}_{\mu \nu },$$

$${\mu }^{\prime }=\kappa {\rho }^2{A}^4(\phi )\widetilde{ϵ}+\frac{1}{2}\rho (\rho -2\mu ){\psi }^2,$$

$${\nu }^{\prime }=2\kappa \frac{{\rho }^2{A}^4(\phi )\tilde{p}}{\rho -2\mu }+\rho {\psi }^2+\frac{2\mu }{\rho (\rho -2\mu },$$

$${\phi }^{\prime }=\psi ,$$

$${\psi }^{\prime }$$

$${\tilde{p}}^{\prime }=-(\widetilde{ϵ}+\tilde{p})\left[\frac{{\nu }^{\prime }}{2}+\alpha (\phi )\psi \right],$$

$${\bar{m}}^{\prime }=4\pi G\tilde{\rho }{A}^3(\phi )\frac{{\rho }^2}{\sqrt{1-2\mu /\rho }},$$

$${\omega }^{\prime }=\varpi ,$$

$${\varpi }^{\prime }$$

Full text

I Introduction

Gravitational waves observations from aLIGO Abbott et al. (2016a, b, 2017a, 2017b, 2017c, 2017d) and the high precision timing of binary pulsars Freire et al. (2012); Kramer (2016); Weisberg et al. (2010); Wex (2014) has allow us to test Einstein’s theory of General Relativity (GR) in the most extreme environments Berti et al. (2015); Yunes and Siemens (2013). With more NS-NS merger events, we will be able to tightly constrain the equation of state (EOS) and other properties of neutron stars (NS) Abbott et al. (2018a, b). Furthermore, as radio astronomers continue to monitor binary pulsars systems, the errors in the timing model parameters will continue to decrease and help place tighter constraints on modified theories of gravity. Nonetheless, in order to investigate how these future observations will help us further test gravity, we first must understand the precise details of how observable predictions are modified.

A popular class of theories in the literature are scalar-tensor theories (STT) of gravity, in which gravity is not only mediated by the metric tensor but also by a long-range scalar field that is non-minimally coupled. Each theory in this class is defined by the choice of conformal coupling function, which mediates the degree of violation of the strong equivalence principle (SEP). Such theories were first studied by Jordan Jordan (1955, 1959), Fierz Fierz (1956), Brans Brans and Dicke (1961), and Dicke (JFBD) as the most natural alternatives to GR, and were later extended By Damour and Esposito-Farése (DEF) Damour and Esposito-Farese (1992, 1993) to include higher order effects. A more recent extension of these theories was introduced by Mendes and Ortiz Mendes and Ortiz (2016) (MO), which introduce a conformal coupling that replicates the behavior of including higher order scalar-field terms in the action Damour and Esposito-Farese (1996); Salgado et al. (1998); Lima et al. (2010); Pani et al. (2011).

Solar system observations, like that of the perihelion shift of Mercury and of the Shapiro time delay, are able to place tight bounds on the parameters of STTs Bertotti et al. (2003); Will (2014). However, these are only measurements in the weak field regime where the gravitational potential is small. STTs are able to satisfy weak field constraints and still produce strong field deviations from GR through a phenomenon known as scalarization Damour and Esposito-Farese (1992, 1993); Barausse et al. (2013); Palenzuela et al. (2014). Thus, one needs observations that probe the strong field, regions where non-linear effects like scalarization occur, in order to place tight constraints on STTs. While gravitational waves observations of binary NS coalescences with aLIGO and other detectors will be able to accomplish this regularly in the future, binary pulsar experiments are already able to probe and constrain STTs in the strong field to extremely high precision. By modeling the time of arrival (TOA) of pulses emitted from pulsar systems Damour and Deruelle (1985, 1986); Damour and Taylor (1992), one can take the observed data and place constraints on the underlying theory of gravity governing the motion of the binary. For this reason, binary pulsars are currently one of the best available testbeds for gravity.

To perform any test, however, one must first know precisely how observables are modified in STTs, and this depends on the so-called scalar charges, i.e. scalar quantities that determine how strongly a scalar field is sourced by an isolated neutron star. For example, the scalar charges enter directly into the parameterized-post-Keplerian (PPK) parameters used in testing STTs Damour and Taylor (1992); Damour and Esposito-Farese (1996); Horbatsch and Burgess (2012); Damour (2007); Esposito-Farese (2004) with binary pulsars. To find these scalar charges, one must numerically obtain NS solutions in STTs, subject to certain physical boundary conditions at the core of the star and at spatial infinity. This numerical process is, at first sight, simple, yet in practice it need not be, mainly for the following two reasons.

First, some of these charges can be numerically challenging to compute. The dominant scalar charge can indeed be read out easily from the leading 1/r piece of the scalar field at spatial infinity, given any numerical solution. Other scalar charges, however, depend on derivatives of the dominant charge with respect to the gravitational mass of the neutron star. During scalarization, the dominant scalar charge can change quite abruptly with respect to the asymptotic value of the scalar field, holding the baryonic mass of the star constant. This, in turn, leads to large spikes in the derivatives, which can be difficult to resolve if one is not careful.

Second, tests of STTs requires knowledge of the scalar charges everywhere in parameter space, and this can be computationally costly. Markov-Chain Monte-Carlo (MCMC) methods that explore the posterior probability distribution of parameters requires the evaluation of the likelihood function hundreds of thousands to millions of times. Every evaluation requires the scalar charges, and if these have to be numerically computed every time the likelihood is needed, the MCMC exploration becomes computationally prohibitive. Clearly then, the MCMC exploration of parameter space in STTs would greatly benefit from a “bank” of calculated scalar charges.

In this paper, we study, calculate, tabulate, and analytically explore the three main scalar charges ($\alpha_{A}$, $\beta_{A}$, and $k_{A}$ Damour and Esposito-Farese (1996)) that enter the PPK parameters of binary pulsar observations. We carry out this calculation both in the STT proposed by Damour-Esposito-Farése Damour and Esposito-Farese (1992, 1993) and the one studied by Mendes-Ortiz Mendes and Ortiz (2016) (hereafter referred to DEF theory and MO theory respectively). We explore a very large region of the parameter space spanned by the two coupling constants of these theories, using 11 different equations of state that are all consistent with neutron stars heavier than $2M_{\odot}$, including a few that are also consistent with the recent gravitational wave observation of a neutron star coalescence Abbott et al. (2017d).

Our main result is the construction of an accurate bank of scalar charges in these two theories that can now be used used in Bayesian model selection and parameter estimation studies of tests of STTs with binary pulsar observations. This bank is constructed both through direct numerical calculations, as well as through the exploration of certain analytic scaling relations. We determine the regime of parameter space in which the latter hold, and when they do, we use them to greatly accelerate the calculation of scalar charges in these regions of parameter space. The end result is a numerically accurate and dense bank of scalar charges that can be interpolated if necessary to provide charges everywhere in parameter space.

The remainder of this paper presents the details of the calculation summarized above. Section II.1 covers the basics of STTs along with observational constraints and includes a discussion of the scalar charges. Section III describes how NSs behave in STTs and begins to set up our numerical scheme to solve for the scalar charges. Section IV details the calculations behind the scalar charges and what numerical techniques are needed to accurately explore the parameter space. Section V provides a detailed description of our publicly available data files, along with instructions for how to use them and what their limitations are. Section VI concludes with a discussion of future work.

II Scalar-Tensor Theories

In this section, we introduce the basics of the class of STTs that we consider and establish the notation to be used in the rest of the paper. For completeness we present this class of theories from first principles using an action and provide a summary of the calculations needed to reach the field equations; we refer the reader to Will and M (1993) for further details. We then describe the two STTs we study in this paper and present the current constraints on these theories. We conclude this section with a discussion of scalarization and the definition of scalar charges in these theories.

II.1 Background and field equations

In general, the massless STTs that we consider can be described in the Jordan frame by an action of the form $\tilde{S}=\tilde{S}_{g}+\tilde{S}_{\mbox{\tiny mat}}$, with the gravitational part taking the form

[TABLE]

where $\tilde{g}$ and $\tilde{R}$ are the determinant and Ricci scalar of the metric $\tilde{g}_{\mu\nu}$ respectively, $\phi$ is a scalar field, $\omega(\phi)$ is a coupling function of the scalar field, and $\kappa=4\pi G/c^{4}$ with $G$ the bare gravitational constant. The matter part of the action, $\tilde{S}_{\mbox{\tiny mat}}[\chi,\,\tilde{g}^{\mu\nu}]$, is a functional of the matter fields $\chi$ that couple directly to the Jordan-frame metric. Therefore, the STTs we study in this paper are metric theories, and as such, laboratory clocks and rods measure time intervals and distances associated with $\tilde{g}_{\mu\nu}$.

While STTs can be completely described using the Jordan-frame action above, it is far more convenient to perform a conformal transformation that puts the action in a form reminiscent of the Einstein-Hilbert action. Let us then consider the transformation $\tilde{g}_{\mu\nu}=A(\varphi)g_{\mu\nu}$, where $g_{\mu\nu}$ is the Einstein-frame metric111From this point on, we use an overhead tilde to represent quantities that are specifically in the Jordan frame. Quantities without overhead tildes should be assumed to be in the Einstein frame., so that the action becomes

[TABLE]

where $g$ and $R$ are the now the determinant and Ricci scalar associated with the Einstein-frame metric $g_{\mu\nu}$. Notice that the matter fields now couple to $A^{2}(\varphi)g_{\mu\nu}$, and therefore, matter no longer falls along geodesics of the metric $g_{\mu\nu}$, but rather it is now also influenced by the scalar field $\varphi$.

The conformal transformation that takes Eq. (1) into Eq. (2) requires the conformal factor

[TABLE]

which then leads to a direct relation between $\varphi$ and $\phi$, given explicitly as

[TABLE]

One can think of $\alpha(\varphi)$ as the gradient of some “conformal potential” Anderson et al. (2016); Anderson and Yunes (2017); Damour and Nordtvedt (1993a, b), defined by $V_{\alpha}\equiv\ln A(\varphi)$, and one can denote the “curvature” of this potential as $\beta(\varphi)=d\alpha/d\varphi$. The conformal potential is a simple way to understand and visualize the coupling between matter and the scalar field, which will be directly quantified by its slope and curvature in the field equations. Clearly then, the choice of $A(\varphi)$, or any of the above quantities for that matter, defines a particular member of this general class of STTs, which is ultimately a choice of exactly how the scalar field affects matter.

The variation of the Einstein-frame action with respect to the dynamical fields, $g_{\mu\nu}$ and $\varphi$, yields the field equations

[TABLE]

where the Einstein frame stress-energy tensor is defined by

[TABLE]

and $T^{\mbox{\tiny mat}}\equiv g_{\mu\nu}T_{\mu\nu}^{\mbox{\tiny mat}}$ is its trace. The stress-energy tensor in the Einstein frame can be related to its Jordan-frame counterpart via the relation $T_{\mu\nu}^{\mbox{\tiny mat}}=A^{2}(\varphi)\tilde{T}_{\mu\nu}^{\mbox{\tiny mat}}$ by applying the conformal transformation to Eq. (7). Assuming a perfect fluid description of the stress-energy tensor [see e.g. Eq. (30) below] allows one to derive the relations $\epsilon=A^{4}\tilde{\epsilon}$ and $p=A^{4}\tilde{p}$ between the energy density and pressure of the fluid in the different frames.

II.2 Scalar-tensor models

Scalar-tensor theories of the form described in the previous subsection allow for deviations from GR because of the new coupling that exists between matter and the scalar field. The particular choice of the conformal factor $A(\varphi)$, and likewise the formal coupling $\alpha(\varphi)$, defines the theory and plays a crucial role in understanding the observable modifications a theory predicts. The first model we consider is DEF theory and it is defined by

[TABLE]

where $\beta_{0}$ is a free coupling parameter. Aside from the JFBD theory in which $A(\varphi)=e^{\alpha_{0}\varphi}$ , this is the simplest massless STT one can consider. One notices that the conformal potential is exactly a parabola whose curvature is precisely determined by the free parameter $\beta_{0}$.

The other model we consider, which has gained attention in the past few years, is MO theory and it is defined by

[TABLE]

where again $\beta_{0}$ is a free coupling parameter. The MO theory was introduced as an analytic approximation to a more fundamental theory that includes quadratic terms of the scalar field coupled to curvature in the action Damour and Esposito-Farese (1996); Salgado et al. (1998); Birrell and Davies (1984); Pani et al. (2011); Mendes and Ortiz (2016). This theory is functionally equivalent to DEF theory in the limit that $\varphi\rightarrow 0$, but it has strictly different behavior when the combination $\beta_{0}\varphi\neq 0$. Therefore, these theories have distinctly different properties, and therefore, modify observables in different ways Mendes and Ortiz (2016); Anderson and Yunes (2017).

II.3 Solar System Constraints

In principle, observations we make, whether they be in the solar system Will (2014); Bertotti et al. (2003) or of binary pulsar systems Kramer (2016); Wex (2014); Damour (2007), constrain the free parameters of the theory. Let us then consider how observables are modified in STTs. As an example, let us first consider the local value of Newton’s gravitational constant. This quantity is given by

[TABLE]

where $G$ is the bare gravitational constant appearing in the action, and an $\infty$ subscript denotes quantities evaluated at $\varphi=\varphi_{\infty}$, e.g. $A_{\infty}\,=\,A(\varphi_{\infty})$. The correction to the gravitational constant causes bodies to accelerate differently depending on the magnitude of scalar field, the parameters of the theory, and the bodies’ composition through violations of the strong-equivalence principle.

The choice of coupling parameter also determines the local value of the parameterized post-Newtonian (PPN) parameters Nordtvedt (1970); Nordtvedt and Will (1972); Will and Nordtvedt (1972). Scalar-tensor theories are a class of fully conservative theories and therefore only pose modifications to the $\gamma_{\mbox{\tiny PPN}}$ and $\beta_{\mbox{\tiny PPN}}$ parameters Will and M (1993); Will (2014). The former is given by

[TABLE]

while the latter is given by

[TABLE]

where as before $\alpha_{\infty}=\alpha(\varphi_{\infty})$ and $\beta_{\infty}=\beta(\varphi_{\infty})$. The $\gamma_{{\mbox{\tiny PPN}}}$ parameter is a measure of the spatial curvature induced by a unit rest mass, while the $\beta_{{\mbox{\tiny PPN}}}$ parameter is a measure of the amount of non-linearity in the superposition law for gravity. The $\gamma_{{\mbox{\tiny PPN}}}$ parameter has been measured from the Shapiro time delay observed by the Cassini spacecraft Bertotti et al. (2003); Will (2014), and it is constrained to $\left|1-\gamma_{\mbox{\tiny PPN}}\right|\lesssim 2.3\times 10^{-5}$. The $\beta_{{\mbox{\tiny PPN}}}$ parameter is measured from observations of the perihelion shift of Mercury Will (2014), and it is constrained to $\left|1-\beta_{\mbox{\tiny PPN}}\right|\lesssim 8\times 10^{-5}$.

The notation we have used above is slightly different from what is found in the literature so let us clarify this here. Typically, instead of using $\alpha_{\infty}$ and $\beta_{\infty}$, some papers that studied DEF theory used a different set of parameters $\{\alpha_{0},\,\beta_{0}\}$. This is because if one modifies Eq. (10) in DEF theory to

[TABLE]

and sets $\varphi_{\infty}=0$, then $\alpha(\varphi_{\infty})=\alpha_{0}$ and $\beta(\varphi_{\infty})=\beta_{0}$, and all observables can be entirely parameterized by the set $\{\alpha_{0},\beta_{0}\}$. This parameterization is identical to our description of DEF theory appearing in Eq. (8), provided that one enforces $\varphi_{\infty}=\alpha_{0}/\beta_{0}$ Damour (2007), which is the choice we make in this paper. When considering MO theory, however, $\alpha(\varphi_{\infty})\neq\alpha_{0}=\varphi_{\infty}\beta_{0}$ and $\beta(\varphi_{\infty})\neq\beta_{0}$, as one can easily see from Eqs. (12)-(15).

In this paper, we want both theories to share the same free parameters $\{\alpha_{0},\,\beta_{0}\}$ and, therefore, the quantities that enter the PPN parameters, $\{\alpha_{\infty},\,\beta_{\infty}\}$, are different functions of $\{\alpha_{0},\,\beta_{0}\}$ in the two theories. These functions are

[TABLE]

in DEF theory, and

[TABLE]

in MO theory. These choices have the advantage of reducing $(\alpha_{\infty},\beta_{\infty})$ to the known relations of DEF theory, while properly generalizing them to MO theory.

II.4 Scalarization and binary pulsar constraints

Solar system observations have the ability to place tight constraints on STTs through weak field observations Will (2014). STTs, however, are able to satisfy these constraints and still deviate substantially from GR inside and near NSs Damour and Esposito-Farese (1992, 1993); Barausse et al. (2013); Palenzuela et al. (2014). The strong field deviations are caused by a phenomenon known as scalarization, in which the scalar field can grow rapidly towards order unity inside NSs even when the asymptotic value, that which is constrained by Solar System observations, approaches zero.

The key behind this rapid growth is the existence of an instability in the field equations when a star reaches a sufficiently large compactness. The onset of this instability is analogous to spontaneous magnetization in ferromagnets Damour and Esposito-Farese (1996). To understand this, consider the external scalar field far from a neutron star, labeled $A$, $\varphi\,=\,\varphi_{\infty}+G\,\omega_{A}/r+\mathcal{O}(1/r^{2})$, where $\omega_{A}$ is a type of “charge” that is energetically conjugate to the external scalar field,

[TABLE]

with $m_{A}$ the total gravitational mass of the NS. When one considers a sequence of neutron stars of masses $m_{A}$, $\omega_{A}$ can become suddenly non-zero at a critical value of the mass or compactness of the star. This sudden activation of the scalar field is what is referred to as spontaneous scalarization. When a NS is scalarized, $\omega_{A}\neq 0$ and the scalar field is excited above its background value, leading to local gravitational effects that will generically be different than those in GR.

For binary pulsar tests, it is convenient to introduce certain quantities that enter the PPK parameters of the binary pulsar timing model. We call these parameters scalar charges in this paper, the first of which is defined by

[TABLE]

which plays the role of an effective coupling between the scalar field and the $A$th NS in the binary. This quantity is the strong field counterpart of the $\alpha_{\infty}$ parameter introduced in the previous section. Similarly, there is a strong field counterpart to the $\beta_{\infty}$ parameter of the previous section, namely

[TABLE]

which encodes higher order effects associated with the exchange of multiple scalar particles between the binary components. Lastly, there is one more charge that enters the PPK parameters, namely

[TABLE]

where $I_{A}$ is the moment of inertia of the NS. Similarly to how $\alpha_{A}$ was an effective coupling between the mass of the NS and the scalar field, this quantity acts as a coupling between the field and the star’s spin angular momentum, and it describes how the NS’s inertia reacts to the presence of an external scalar field. All of these scalar charges must be calculated while holding the baryonic mass of the star constant, as they measure the “sensitivity” of a star to the external scalar field. The main goal of this work is to calculate these scalar charges and make them publicly available. As such, we will provide the details of these calculations later in §IV after we introduce the relevant framework needed to understand NSs in STTs.

III Compact Stars in Scalar-Tensor Gravity

In this section we discuss how compact stars behave in STTs and introduce some of the basics of scalarization from an analytic perspective. We focus our attention on isolated, slowly-rotating stars because the binary pulsar systems that we wish to provide scalar charges for are widely separated. We begin with a discussion of the interior spacetime of slowly-rotating compact stars and then discuss the exterior spacetime and how one connects it to the interior. Following the discussion of the field equations, we discuss the different types of equations of state one can use and conclude with an analytic discussion of scalarization in the different regions of parameter space.

Before we begin, however, let us discuss how we will describe compact stars in STTs. We focus on a stationary, axisymmetric spacetime, which allows for a description of a star that is slowly rotating. Following closely the work of Damour and Esposito-Farese (1996), we make the metric ansatz proposed by Hartle Hartle (1967)

[TABLE]

in which one only keeps terms to first order in the star’s angular velocity $\Omega=u_{\phi}/u_{t}$, where $u^{\mu}$ is the fluid’s four velocity. In this ansatz, the metric functions $(\nu,\lambda)$ are zeroth-order in rotation and functions of the radial coordinate $\rho$ only, while the metric function $\omega$ is first-order in rotation and depends both on $\rho$ and the polar angle $\theta$. For practical and physical convenience, we also redefine the $g_{rr}$ component of the metric through the interior mass function $\mu(\rho)$ defined via

[TABLE]

Moreover, we model the NS matter with a perfect fluid stress-energy tensor given by

[TABLE]

where $\epsilon$ is the fluid’s energy density and $p$ is the pressure, both in the Jordan frame.

III.1 Interior Spacetime

The scalar-tensor field equations with this metric ansatz are similar to those in GR. At zeroth-order in rotation, the field equations for the diagonal components of the metric are identical to those in the spherical case, which have already been studied in detail in the literature Damour and Esposito-Farese (1996). At first order in rotation, one finds a second-order equation for the $\omega$ metric function, which can be converted into an ordinary, second-order differential equation in $\rho$ through a Legendre decomposition in $\theta$ Damour and Esposito-Farese (1996). Such a decomposition reveals that only the $\ell=1$ mode in the Legendre decomposition has support. In particular, the field equations can be written in the first-order form Damour and Esposito-Farese (1996)

[TABLE]

where primes denote derivatives with respect to the radial coordinate $\rho$. Notice that we have also included an equation for the enclosed baryonic mass of the star $\bar{m}(\rho)$, which gives the total baryonic mass through $\bar{m}_{A}=\bar{m}(R)=\int\tilde{\rho}\,u^{t}\sqrt{-\tilde{g}}\,d^{3}x$, where $R$ marks the surface of the star, where by definition the pressure vanishes. In Eq. (31) we have explicitly used the Jordan-frame fluid variables $\tilde{\epsilon}$, $\tilde{\rho}$, and $\tilde{p}$ since these are the physical quantities that are measured by observations and appear in the EOSs that are discussed below.

From a numerical standpoint these equations pose a problem near the center of the star at $\rho=0$ since the equations diverge there. The proper way to deal with this is to expand the equations about $\rho=0$ and start the numerical integration at some arbitrary small distance away from the center, say $\rho_{\min}$. Expanding Eqs. (31) and evaluating at $\rho_{\min}$ gives the boundary conditions

[TABLE]

where $\tilde{\epsilon}_{c}=\tilde{\epsilon}(\tilde{\rho}_{c})$ and $\tilde{p}_{c}=\tilde{p}(\tilde{\rho}_{c})$ are the central values of the Jordan-frame energy density and pressure respectively and are defined by the EOS. The values of $\varphi_{c}$ and $\tilde{\rho}_{c}$ are chosen independently and define a particular solution to the field equations. Equations (31) and (32h) allow one to integrate the field equations from $\rho=\rho_{\min}$ to any arbitrary radius, even outside the star as the field equations describe the entire spacetime. A more computationally efficient method, however, is to use an analytic solution in vacuum that is valid in the exterior of the star, and then, match it to the interior solution at $\rho=R$, as we describe in the next section.

III.2 Exterior spacetime

The exterior solution to the field equations [Eqs.(31)] were found by Just in the late 1950s Horbatsch and Burgess (1959). In the coordinates introduced by Just, the exterior metric takes the form

[TABLE]

where one has

[TABLE]

while the scalar field takes the form

[TABLE]

where $a$, $b$, and $d$ are integration constants and are constrained by the relation $a^{2}-b^{2}=4d^{2}$. The integration constants can all be expressed in terms of the gravitational mass of the star $m_{A}$ and some effective coupling constant $\alpha_{A}$, which in this context plays the role of the scalar charge in Eq. (25). These relations take the form

[TABLE]

The coordinates used in the Just spacetime above are not the same as those used in the Hartle ansatz of the interior metric. Comparing the two line elements, one finds that

[TABLE]

Given these relations, we can now read off the total gravitational mass of the star $m_{A}$ from the $1/\rho$ behavior of $g_{tt}$, or $g_{\rho\rho}$, and the star’s $z$-component of the total angular momentum $J_{A}$ from the $1/\rho^{2}$ portion of $g_{t\phi}$. We can recast the later in terms of the $1/\rho^{3}$ behavior of $\omega$ as

[TABLE]

in which case the moment of inertia follows as

[TABLE]

By directly integrating the equation for $\omega(\rho)$ and inserting Eqs. (37b) into Eqs. (31), we arrive at a set of relations that are valid at the stellar surface Damour and Esposito-Farese (1996) (with a subscript $s$ denoting values at the surface of the star), namely

[TABLE]

These set of relations allow us to extract important observables at $\rho=\infty$ by simply knowing their appropriate values at the surface of the star.

III.3 Equations of State

The description of the problem is not complete without first knowing how the fluid properties $\{\tilde{\rho},\,\tilde{p},\,\tilde{\epsilon}\}$ depend on one another. This is accomplished by an EOS and it allows us to close the system of equations. In this subection we describe the three types of EOSs considered in this paper in detail.

III.3.1 Polytropes

The most simple EOS to consider for NSs is a polytropic equation of state in which the the pressure and baryonic density are related through a power law, i.e.

[TABLE]

where $K$ is the polytropic constant and $\Gamma$ is the adiabatic index of the fluid. A particular choice for $K$ and $\Gamma$ define the EOS as well as some macroscopic properties of the NS, such as the maximum mass and compactness. In Damour and Esposito-Farese (1996), a polytropic EOS was used in the calculation of the “gravitational form factors” (what we call the scalar charges in this paper) and therefore we will use the same polytropic EOS here to validate our computational algorithm. In particular, we will choose

[TABLE]

with a fiducial baryonic mass density $\tilde{\rho}_{0}=1.66\times 10^{14}$ g cm*-3*, to make comparisons to the results in that paper, which we present later in Fig. 3.

The baryonic mass density only appears in the integral for calculating the baryonic mass of the star, and therefore, we need another relation relating pressure and baryonic density to the total energy density. The first law of thermodynamics provides such a relation:

[TABLE]

Equations (42) and (43) can now be inserted directly into Eqs. (31) to provide a complete description of the interior of the NS.

While polytopes are very simple analytic EOSs that facilitate quick numerical calculations, they are an oversimplification of the true microphysics occurring inside the NS. Polytropes assume the same functional dependence between pressure and density throughout the entire star and do not account for any real differences that exist in different density regimes. Due to this reason, we only use polytopes as a proof of concept for the existence of scalarization and to make valid comparison between our results and those presented in Ref. Damour and Esposito-Farese (1996) to ensure numerical consistency.

III.3.2 Tabulated Equations of state

A complete description of the microphysics occurring inside NSs requires the full modeling of $N$-body quantum systems at extremely high pressures and densities. The calculations required to solve for these relations is very expensive and not practical on the fly for every density and pressure inside a NS. Moreover, because NS type densities cannot be observed in laboratories on Earth there is uncertainty on what physics is actually taking place at these densities. There have been numerous models proposed to describe matter at supra-nuclear densities and they have been tabulated such that one can interpolate them as needed.

For the purpose of this paper, we consider a wide range of tabulated EOSs that produce NSs with masses that are consistent with observations, most notably that of a near $2M_{\odot}$ pulsar in J0348+0432. We consider 11 different tabulated EOS Lattimer and Prakash (2001) that satisfy this constraint: AP3-4 Akmal et al. (1998), ENG Engvik et al. (1996), H4 Lackey et al. (2006), MPA1 Muther et al. (1987), MS0 Mueller and Serot (1996), MS2 Mueller and Serot (1996), PAL1 Prakash et al. (1988), SLy4 Douchin and Haensel (2001) 222The SLy4 EOS we use here is commonly denoted as simply SLy in the literature. , and WFF1-2 Wiringa et al. (1988). All but one (H4) of these EOSs contain plain nuclear matter and do not contain any form of strange matter. Because these EOS arise from numerical calculations that include true microphysics (or at least various justified approximations) we consider these to be the most physically relevant of the EOSs that we consider and are the ones we use for our main results. More importantly, however, many of these EOSs are consistent with aLIGO’s constraints placed from the observation of coalescing NSs Abbott et al. (2018a).

III.3.3 Piecewise Polytropes

A useful compromise between the simple polytropic and tabulated EOSs is a piecewise polytropic model that stitches together multiple polytropes in different density regions inside the NS. In particular, we consider the piecewise polytropic EOSs studied in Ref. Read et al. (2009) in which the authors developed a parameterized model that can accurately capture the feature of tabulated EOS. Of the 34 EOSs that the authors fit their model to, 8 of them overlap with the set of tabulated EOSs that we consider in this paper333One might notice that, aside from the PAL1 EOS that is explicitly not contained in Ref. Read et al. (2009), there are 10 EOSs listed in the previous section that correspond to the EOSs that were fit in this paper. The data we received from Norbert Wex came from the original work by Lattimer and Prakash Lattimer and Prakash (2001) and in fact the EOSs labeled MS0 and MS2 are different from those of the same name appearing in Ref. Read et al. (2009). There appears to be a mismatch in nomenclature that we feel is worth pointing out.. While these approximations have been used extensively in the literature for their convenience, we later investigate how these approximations affect the scalar charges that are developed in NSs. Hence, let us briefly discuss these EOSs, referring to Ref. Read et al. (2009) for a more detailed and complete description.

Similarly to the standard polytropic EOS, the various regions inside of the NS are described by a polytrope of the form

[TABLE]

where now $\Gamma_{i}$ is the adiabatic index for the $i$th region of the NS and $K_{i}$ is the polytopic constant chosen to ensure continuity at the boundaries between regions. Similarly to the single polytrope case, the first law of thermodynamics in Eq. (43) is used to find the energy density in each region

[TABLE]

where

[TABLE]

is an integration constant that must be present in order to ensure that all fluid variable are continuous across the boundaries between regions. For the single polytrope case, the requirement that $\tilde{\epsilon}/\tilde{\rho}=1$ in the limit that $\tilde{\rho}\rightarrow 0$ forces $a=0$ and thus reduces Eq. (45) to Eq. (43).

We adopt a low-density EOS for the crust of the NS that is identical to the one presented in Table II of Ref. Read et al. (2009) where the SLy (SLy4 as appearing in this paper) is independently fit for $\tilde{\rho}\lesssim 10^{12}$ g/cm3. The boundary between the crust and high-density EOS is then determined by the intersection of the respective polytropes and is ultimately determined by the value of $\Gamma_{1}$. Then, at a fixed baryonic density of $\tilde{\rho}=10^{14.7}$ g/cm3 and best fit pressure $\tilde{p}_{1}=\tilde{p}(\tilde{\rho}_{i})$ the first region is matched to a second region with adiabatic index of $\Gamma_{2}$ and polytropic constant $K_{2}$. Another match is then performed at another boundary $\tilde{\rho}=10^{15}$ g/cm3 to an even higher density region with constants $\Gamma_{3}$ and $K_{3}$. The complete set of parameters $\{\log(\tilde{p}),\,\Gamma_{1},\,\Gamma_{2},\,\Gamma_{3}\}$ represents the best fit values found in Table III of Ref. Read et al. (2009). For our purposes, we focus exclusively on the approximations of AP3 and SLy4 to compare scalar charges between tabulated and piecewise polytropic EOSs.

IV Calculating the Scalar Charges

Now that we have a full description of the problem at hand, we solve the equations numerically to obtain the relevant scalar charges appearing in Eqs (25)-(27). We begin this section with an analytic description of the scalar charges, and then proceed with a description of our numerical methods for solving the field equations and extracting the various scalar charges. We then recap, but in more detail, the different regions of parameter space that we are concerned with and discuss their importance. We first present a comparison between some of our results and those originally found in Ref. Damour and Esposito-Farese (1996), and then discuss our full results in each region of parameter space.

IV.1 Analytic investigation of scalarization

The phenomenon of scalarization has been well studied in the literature over the past decades, particularly in the context of spontaneous scalarization occurring when $\beta_{0}\lesssim-4.3$. However, it is useful to review what happens outside of this regime as well since there are still non-linear effects coming into play, particularly when $\beta_{0}\gtrsim-3.5$. In this subsection we review the analytics that help guide our calculations in the different regions of parameter space.

Scalarization can be understood from an analytic standpoint when one makes a few simple approximations. In both theories we consider here, the conformal coupling takes the form $\alpha(\varphi)=\beta_{0}\varphi+\mathcal{O}(\varphi^{2})$ in the limit that $\beta_{0}\varphi$ is relatively small compared to unity. Let us now consider the field equation for the scalar field in Eq. (6), but instead of considering the full nonlinear equation, we make a weak field approximation such that $\Box\rightarrow\delta_{ij}\nabla_{i}\nabla_{j}\rightarrow\nabla_{r}^{2}$ with the last term being the radial portion of the flat space Laplacian in spherical coordinates. We assume that $T_{\mbox{\tiny mat}}$ is constant since we are considering weakly gravitating systems. As mentioned in Damour and Esposito-Farese (1993), we do not expect the trace of the stress-energy tensor to be negative for weakly gravitating systems, but it is fruitful to leave its sign general and consider the full breadth of parameter space with the same analytics.

The assumptions made thus far allow one to write the equation of motion for the scalar field as

[TABLE]

where we have introduced the constant $C^{2}=\kappa|\beta_{0}\,T_{\mbox{\tiny mat}}|$ for $r<R$, which vanishes when $r>R$. The field equation is still subject to the same boundary conditions as before, and thus, $\varphi(r=0)=\varphi_{c}$ and $\varphi^{\prime}(r=0)=0$ to ensure regularity at the center and to ensure the scalar field is continuous and differentiable at the surface.

Now we must solve Eq. (47) inside and outside the star, subject to the boundary conditions above. The exterior solution takes the form

[TABLE]

where $\varphi_{\infty}$ and $\omega_{A}$ are integration constants. There are two interesting scenarios that can occur within the star, namely when the product $\beta_{0}\,T_{\mbox{\tiny mat}}$ is positive and when it is negative. In the positive case, the interior solution for the scalar field takes the form

[TABLE]

The asymptotic value of the scalar field at infinity $\varphi_{\infty}$ can made chosen to be arbitrarily small, but even in the case of $\varphi_{\infty}\rightarrow 0$, the central value of the scalar field can still be non-zero if $\cos(CR)\to 0$ at the same rate. While this is an over-simplified description of the problem, one in which we essentially are assuming a constant density inside the star and ignoring non-linear effects, it does demonstrate that there can be non-trivial scalar field solutions even when one forces the asymptotic value of the scalar field to vanish.

In the situation where the product $\beta_{0}\,T_{\mbox{\tiny mat}}$ is negative, there exists an opposite effect inside the star: any deviations from GR are exponentially suppressed. The negative case leads to an interior solutions of the form

[TABLE]

in which case any non-vanishing value of $\varphi_{\infty}$ is suppressed even further by $\cosh\left(CR\right)$. The suppression mechanism drives the STT solution to GR inside the star when $\beta_{0}\,T_{\mbox{\tiny mat}}<0$.

In this linear $\varphi$ regime, the quantity $G\omega_{A}$ appearing in Eq. (48) can be expressed as

[TABLE]

in the case of $\beta_{0}\,T_{\mbox{\tiny mat}}>0$ and with the tangent exchanged for a hyperbolic tangent when $\beta_{0}\,T_{\mbox{\tiny mat}}<0$. Recalling from Sec. II that $\alpha_{A}\,=\,-\omega_{A}/m_{A}$ and that $\varphi_{\infty}=\alpha_{0}/\beta_{0}$ one finds that $\alpha_{A}\propto f(\beta_{0},m_{A})\alpha_{0}$444Technically, this relation should read $\alpha_{A}\propto f(\beta_{0},m_{A})\alpha_{\infty}$ as the scalar charge should always reduce to its weak field counter part $\alpha_{\infty}$ in the absence of strongly self-gravitating matter. However, these relations are already derived under weak field assumptions and in both theories $\alpha_{\infty}\eqsim\alpha_{0}$ in this regime., as long as one can neglect non-linear interactions of the scalar field. While many of the most interesting effects of STTs, like spontaneous scalarization, occur in the most non-linear regions of parameter space, this simple relations provides valuable insight into the types of solutions one would expect in other regions of parameter space.

IV.2 A classification of parameter space

Let us now use the analytic insight described above to guide us in our numerical exploration of the $\{\alpha_{0},\,\beta_{0}\}$ parameter space shown in Fig. 1. We will break our investigation into six distinct regions in parameter space, each of them having distinct features that need to be handled differently when solving for the scalar charges numerically. In all of these regions we must set up some numerical grid in $\alpha_{0}$, $\beta_{0}$, and $m_{A}$, and these grids are precisely how each of these regions differ from one another. We discuss the various regions in Fig, 1 in detail below, relying on the analytic insight above and our numerical results, and in the next subsections we will present the numerical techniques used in each region and our results.

There exist three regions, I–III (brown, red, and orange respectively) in Fig. 1 in which spontaneous scalarization occurs, i.e. when $\beta_{0}\leq-4.3$. In all of these regions there exists a phase transition in the scalar field, with the sharpness of the transition being determined by the value of $\alpha_{0}$ (lower values leading to more steep transitions). Due to the varying level of steepness in the phase transitions and the numerical limitations of taking derivatives, we have broken this region of $\beta_{0}$ into 3 distinct subregions in which we use different numerical techniques to most efficiently explore the parameter space. In short, region I in Fig. 1 contains the sharpest transitions, and therefore requires a finer numerical grid in $m_{A}$ to resolve the relevant features of interest. Region II contains less sharp transitions and it can be accurately explored with less grid points. Regions I and II both contain the same grid spacing in $\log_{10}(\alpha_{0})$ and $\beta_{0}$, but a different spacing in $m_{A}$. Lastly, region III contains the same grid spacing in $m_{A}$ and $\beta_{0}$ as region II, but it contains more grid points in $\alpha_{0}$ to allow us to accurately calculate the numerical derivatives necessary for the scalar charges.

Region V (cyan) in Fig. 1 is where the non-linearity of the scalar field can effectively be neglected and hence the scalar charge turns out to scale as in Eq. (51). When we solve the full set of field equations we make no approximations, but the resulting scalar charges do indeed follow the scaling relation $\alpha_{A}\propto f(\beta_{0},m_{A})\alpha_{0}$. Thus, in region V we solve for a single set of solutions at a single value of $\alpha_{0}$, lying on the black dashed lines in Fig. 1 at $\alpha_{0}=10^{-3}$ for $\beta_{0}>-3.5$ and $\alpha_{0}=10^{-4}$ for $\beta_{0}<-3.5$, and we use the scaling relation to populate the entire region. We numerically verify in Sec. A.2 that these scaling relations are indeed accurate when compared to the full numerical exploration of this region of parameter space. This scalable region does not cover the entire parameter space where spontaneous scalarization does not occur. We have found numerically that the scaling does not hold when $\beta_{0}<-3.5$ and when $\alpha_{0}\in(10^{-4},10^{-2})$, which is why we have isolated this part of parameter space to region IV. In this region, we investigate the solutions as if we expected spontaneous scalarization.

For region VI in Fig. 1 we find that the scaling relations previously discussed no longer exist because of how large $\alpha_{0}$ can become. The lack of quasi-analytic solutions here should not come as a surprise because this region has such large values of $\alpha_{0}$ that STT modifications can be easily constrained with solar system observations. For completeness, however, we still investigate this region extensively as some parts of this parameter space are actually useful when placing binary pulsar constraints555In Ref. Freire et al. (2012), for example, there exists a region near $\beta_{0}\sim-2$ in which constraints on scalar dipole radiation fail to constrain STTs better than Cassini and other weak field tests. This “horn” appearing in the binary pulsar constraints occurs for NS-WD binaries in which the quantity $(\alpha_{NS}-\alpha_{WD})^{2}\sim(\alpha_{NS}-\alpha_{0})^{2}$ vanishes, which tends to happen near $\beta_{0}\sim-2$.. We have also removed a portion, region VII in gray, from the parameter space because we are not able to calculate NS solutions here. This has been noted in the literature before when considering DEF theory Mendes and Ortiz (2016); Anderson and Yunes (2017), but in those papers, the authors only investigate small values of $\alpha_{0}$ and focused on values of $\beta_{0}$ considerably larger than what we consider here. Nonetheless, we find results very similar in this gray region of parameter space for both DEF and MO theory and it becomes impossible to extract any useful information from our numerical calculations. Therefore, since binary pulsar typically do not probe this region and Solar System tests have already ruled it out, we are justified in neglecting it.

The final point to discuss regarding the parameter space is the vertical dashed line and numbered points appearing in Fig. 1. The black vertical dashed line at $\beta_{0}=-6$ marks a set of special solutions we have calculated to compare our results to those original found in Ref. Damour and Esposito-Farese (1996). While we do not explore this region of parameter space in depth it provides a useful comparison to validate our code and make comparisons between known results in DEF theory and new results in MO. Details about the numbered points in Fig. 1 can be found in Table 3 and they represent a set of points we use to investigate the error in our numerical results. The details of this error analysis can be found in Sec. A.

IV.3 Numerical methodology

We parameterize our NS solutions by a choice of $\{\tilde{\rho}_{c},\,\alpha_{0},\,\beta_{0}\}$ which ultimately determines the star’s gravitational mass $m_{A}$, baryonic mass $\bar{m}_{A}$, and the asymptotic value of the scalar field $\varphi_{\infty}=\alpha_{0}/\beta_{0}$. We take the approach of solving Eqs. (31) starting from the center of the NS, using Eqs. (32h) to start our numerical integration away from the singularity at $\rho=\rho_{\min}$. We follow the methods employed in Refs. Anderson et al. (2016); Anderson and Yunes (2017) and use Mathematica’s default ODE solver666The default method used is LSODA which is a variant of the original LSODE (Livermore Solver for Ordinary Differential Equations) approach to solving a wide class of differential equations. to integrate the equations to the surface of the NS where we then use Eqs. (40) to extract values at spatial infinity.

In order to begin the integrations, however, we must make a choice of $\tilde{\rho}_{c}$ and $\varphi_{c}$ as these are the two independent parameters appearing in the boundary conditions. It is near impossible to guess the correct value of $\varphi_{c}$ that correspond to a given $\varphi_{\infty}=\alpha_{0}/\beta_{0}$ to within a small tolerance. Therefore, we use a Newton-Raphson shooting method to converge onto the correct value of $\varphi_{c}$. In particular, we solve the equations again with a new central value of the scalar field $\varphi_{c,n+1}=\varphi_{c,n}+\Delta\varphi_{c}$, which gives a slightly different value of $\varphi_{\infty}$. The difference of the extracted values of $\varphi_{\infty}$ allow us to construct a simple difference equation

[TABLE]

where $n$ is the iteration number. The equation above allows us to predict a new value of $\varphi_{c}$ that gives a value of $\varphi_{\infty}$ that is closer to the desired result, and we iterate this process until the resulting value of $\varphi_{\infty}$ is equivalent to $\alpha_{0}/\beta_{0}$ to within some numerical tolerance. At the subsequent point in $\tilde{\rho}_{c}$ we use the previous value of $\varphi_{c}$ as the starting point for the shooting, which typically allows convergence to within numerical tolerance in about 2-3 iterations.

The method described above provides NS solutions corresponding to a single combination of $\{\tilde{\rho}_{c},\,\alpha_{0},\,\beta_{0}\}$ for any choice of theory and EOS. It is convenient that the scalar charge $\alpha_{A}$ can be extracted directly from the boundary conditions at the surface. The quantities $\beta_{A}$ and $k_{A}$, however, must be calculated by taking derivatives across multiple solutions while keeping the baryonic mass constant. The most accurate way to take these derivatives would involve parameterizing the NS solutions by $\{\bar{m}_{A},\,\alpha_{0},\,\beta_{0}\}$ instead and shooting in both $\tilde{\rho}_{c}$ and $\varphi_{c}$. Such an approach would allow one to construct multiple NS solutions with identical values of $\bar{m}_{A}$ and varying values of $\varphi_{\infty}$, corresponding to different values of $\alpha_{0}$, which is required for the derivatives needed to calculate $\beta_{A}$ and $k_{A}$. While this approach works extremely well, it is very computationally expensive since one must now shoot in two dimensions multiple times just to calculate the scalar charges for a single combination of $\{\bar{m}_{A},\,\alpha_{0},\,\beta_{0}\}$. To completely populate the parameter space of interest one must compute the charges for roughly $10^{5}$ combinations of $\{\bar{m}_{A},\,\alpha_{0},\,\beta_{0}\}$ just for a single EOS and theory choice. This quickly becomes cumbersome from a computational standpoint so we decided to take a different approach.

Our computational method is as follows. We continue to parameterize the NS solutions with $\{\tilde{\rho}_{c},\,\alpha_{0},\,\beta_{0}\}$, but rather than focusing on a single value of $\bar{m}_{A}$, we calculate an entire mass-radius (MR) curve of solutions, corresponding to a set of $\{\tilde{\rho}_{c,i}\}$, for a large discrete set of $\{\alpha_{0},\,\beta_{0}\}$ values. This approach allows us to then interpolate the scalar charge $\alpha_{A}$ as a function of the baryonic mass, thus generating curves like those in Fig. 2. Because we have a finely discretized grid in $\bar{m}_{A}$, we can interpolate between points and extract $\alpha_{A}(\bar{m}_{A})$ for any value of baryonic mass, and we can do the same for every curve we calculate. Therefore, we can compute numerical derivatives of the scalar charge (or any other quantity) at any value of baryonic mass in a computationally efficient way. While this method is prone to more numerical error than the previous one, it allows us to sample the entire parameter space very finely and it can be carried out orders of magnitude faster. A discussion of the errors associated with our methods is presented in Sec. A.2

Let us now discuss the way we take numerical derivatives. We choose to use a fourth-order accurate finite difference scheme to calculate the derivatives in Eqs. (26)-(27). For reasons discussed below, we have to use central, forward, and backward finite difference schemes in order to most effectively utilize our numerical grid, and they take the forms

[TABLE]

where $F_{\pm n}=F(\varphi_{0}\pm n\Delta\varphi_{0})$, and $c_{i}$, $f_{i}$, and $b_{i}$ are the corresponding finite difference coefficients for central, forward, and backward derivatives respectively found in Table 1.

Let us now briefly discuss the numerical grid in parameter space. Each region of parameter space in Fig. 1 uses a different numerical grid in $\{m_{A},\,\alpha_{0},\,\beta_{0}\}$. The spacing in $\alpha_{0}$ is determined by the level of accuracy we want when using the various finite difference schemes for the derivatives. Since $\Delta\varphi\propto\Delta\alpha_{0}$ we need to choose our spacing such that $\Delta\alpha_{0}$ between consecutive solutions branches is not too large. Finally, the $\beta_{0}$ grid is determined strictly by the presence of spontaneous scalarization. Therefore, if $\beta<-4.3$, the grid spacing is $\Delta\beta_{0}=0.02$ and otherwise it is $\Delta\beta_{0}=0.1$. The details of the grids in each subspace are presented in the next subsection.

The spacing in $m_{A}$ in each region is determined by the presence of any sharp features that may appear in the solutions, such as spontaneous scalarization. The green regions appearing in the left panels of Fig. 2 represent the lowest resolution portions of our grid, in which we have a grid point every $0.1\rho_{0}$; we call this value $\Delta\rho_{c}^{GR}$ since it is dense enough to accurately reproduce a MR curve in GR. Spontaneous scalarization turns “on” and “off” in the red regions in Fig. 2, and in order to capture the phase transition effectively, we increase our resolution to $\Delta\rho_{c}^{PT}=\Delta\rho_{c}^{GR}/800$, where $PT$ stands for phase transition777We only increase this resolution by a factor of 800 when we consider $\alpha_{0}\leq 10^{-4}$ and $\beta_{0}<-4.3$. When $\alpha_{0}>10^{-4}$, we find that we do not need to sample as many points to fully capture the features of the phase transition, hence we allow $\Delta\rho_{c}^{PT}=\Delta\rho_{c}^{GR}/200$. In regions where spontaneous scalarization does not occur, i.e. $\beta>-4.3$, there is no phase transition and we simply use $\Delta\rho_{c}=\Delta\rho_{c}^{GR}/20$ everywhere.. The blue region in Fig. 2 between the phase transitions is where there are no sharp features, but where we still want increased resolution since non-linear effects do come into play; in these regions, we use a grid spacing $\Delta\rho_{c}^{scal}=\Delta\rho_{c}^{GR}/20$ where $scal$ stands for scalarization. Figure 2 illustrates the density of grid points in these different regions by the number of orange circles appearing along each curve. We limit our solutions to values of $\tilde{\rho}_{c}$ that give a $0.5M_{\odot}$ NS in GR and the $\tilde{\rho}_{c}$ that predicts the maximum mass NS in GR, for each EOS we consider. As in the $(\alpha_{0},\beta_{0})$ grid spacing case, the details of the $\tilde{\rho}_{c}$ grid are presented in the next subsection.

IV.4 Numerical results

IV.4.1 Code Verification: the $\beta_{0}=-6.0$ case

With a basic idea of our numerical grid and numerical methods for solving NS solutions in hand, let us present a comparison between our results and the ones found in Damour and Esposito-Farese (1996). In that study, the authors used the polytropic equation of state described in Sec.  III.3.1 and they show results for the $\beta_{0}=-6$ and $|\varphi_{\infty}|=2.4\times 10^{-3}$ (or $\alpha_{0}=1.44\times 10^{-2}$ in our framework) case, which corresponds to point 1 in Fig. 1. In addition to using the polytropic EOS of Damour and Esposito-Farese (1996), we will also show here results for the tabulated and piecewise polytropic versions of AP3 for both the DEF and MO STTs.

Let us first take a look at the top row of Fig. 3. Comparing our result to those in Damour and Esposito-Farese (1996), we see that the black curves, those for DEF theory, are in great agreement888Noticed that we have used gravitational mass instead of baryonic mass on the horizontal axis.. The red curves in Fig. 3 are the results for the MO theory and we see that in every case the magnitude of the charge is less than those of DEF theory and that they also do not reach very high masses. This feature is general for MO theory, i.e. this feature is not a result of a special choice of $\alpha_{0}$ and $\beta_{0}$. The reasoning behind this is that the curvature of the conformal potential appearing in Eq. (13) is not as large as that of DEF theory, and therefore, the excitation of the scalar field is never as strong once the instability occurs. The curves for $\beta_{A}$ and $k_{A}$ formally diverge at higher masses, a phenomenon due to the fact that these are derivatives of quantities that “turn over” on themselves, e.g. in the right panels Fig. 2 one can see that the slope in $\alpha_{A}$ and $\log_{10}I_{A}$ become infinite at some point as the mass increases.

The bottom row in Fig. 3 shows the same scalar charges but for realistic equations of state, namely AP3 here, and its piecewise polytropic approximate. The maximum values of $\alpha_{A}$ are very weakly affected by the EOS in both theories. There is, however, a strong dependence on the EOS when it comes to the location of the critical mass, $m_{\mbox{\tiny crit}}$, at which spontaneous scalarization occurs, and the maximum mass above which stable NSs do not exist. The other scalar charges, $\beta_{A}$ and $k_{A}$, are quite different between the different EOSs and theories. The phase transition is now less sharp, and therefore, the magnitude of $\beta_{A}$ is smaller than in the polytropic case. Moreover, it happens to be the case here that in MO theory there is no formal divergence in $\beta_{A}$ at large masses. Finally, the “negative spike” in $k_{A}$ that usually occurs is greatly suppressed for realistic EOSs.

The piecewise polytropic EOS leads to scalar charges that are very similar to those found with a tabulated EOS. Aside from a very slight shift in the masses, the structure and magnitude of the scalar charges are nearly identical. The small differences are likely due to the fact that polytropes are just too simple to accurately capture the different physics that occurs at different densities inside NSs, which can over/under exaggerate features we find in the scalar charges. Using the piecewise polytropes, however, can speed up numerical calculations immensely since they are analytic. A more detailed investigation of the differences between tabulated and piecewise polytrope results can be found in Sec. A.3

IV.4.2 Spontaneous scalarization: $\beta_{0}\leq-4.3$

As is well-known in the literature, spontaneous scalarization occurs in STTs when $\beta<-4.3$ regardless of the EOS. As demonstrated in the previous section, there exists a $m_{\mbox{\tiny crit}}$ at which a phase transition of the scalar field occurs, but the precise value of this mass, however, depends on the theory, the EOS, and the value of $\beta_{0}$. These dependencies require that we determine what this critical mass is, and to place a finer grid in $\tilde{\rho}_{c}$ centered near this location to ensure we capture the details of the phase transition with enough accuracy and precision to calculate the scalar charges at these points, cf. Fig. 2. Spontaneous scalarization occurs in regions I–III (brown, red, and orange) of Fig. 1, each of which has a slightly different grid that we discuss in detail next.

Region I (Brown):

This region has the finest grid in $\tilde{\rho}_{c}$ because the phase transition is most sharp here, and near the transition we decrease our spacing in $\tilde{\rho}_{c}$ by a factor of 800 relative to $\Delta\rho_{c}^{GR}$, as mentioned earlier. This level of resolution requires around 600 NS solutions in total for each value of $\alpha_{0}$ and $\beta_{0}$. In this region, we use $\alpha\in\{1.2,\,1.1,\,\cdots,\,0.2\}\times 10^{\text{mag}}$ where $\text{mag}\in\{-4,\,-5\}$, and a spacing in $\beta_{0}$ of $\Delta\beta_{0}=0.02$. Spacing $\alpha_{0}$ in this manner gives a spacing in $\varphi$ of $\Delta\varphi=10^{\text{mag}-1}/\beta_{0}\sim 2\times 10^{\text{mag}}$ for this range of $\beta_{0}$. Since our finite difference schemes are all fourth-order accurate, our spacing in $\alpha_{0}$ is small enough to confidently calculate the needed derivatives999Going to smaller values of $\alpha_{0}$ is difficult because we can never let it change signs when calculating the derivatives. This means that $\Delta\alpha_{0}$ must always be smaller than $\alpha_{0}$, and it is computationally expensive to calculate solutions with enough accuracy and precision for $\Delta\alpha_{0}\lesssim 10^{-6}$.. This particular spacing allows us to use the central finite difference scheme for mantissa of $\alpha_{0}$ from 1.0 to 0.4, the forward finite difference scheme for 0.2 and 0.3, and the backward scheme for 1.2 and 1.1. Since we need to calculate 1.2, 1.1, 0.3, and 0.2 for using the central finite difference scheme, we automatically get the derivatives at these extra points for free just by changing the finite differencing.

Region II (Red):

In this region, we are able to use only a factor of 200 more points near the transition regions, i.e. $\Delta\rho_{c}^{GR}/200$. We use the same grid in $\beta_{0}$ here as in region I above, i.e. $\Delta\beta_{0}=0.02$. Our grid in $\alpha_{0}$ is set up in a similar way as well, i.e. $\alpha\in\{1.2,\,1.1,\,\cdots,\,0.2\}\times 10^{\text{mag}}$ where $\text{mag}\in\{-2,\,-3\}$, and are both set up to make efficient use of the finite differences introduced above. As one can see in Fig. 3 for example, the phase transitions are not as sharp in this region as they are in region I (see Fig. 2), and this allows us to confidently under-sample $m_{A}$ relative to these smaller values of $\alpha_{0}$.

Region III (Orange):

In this region, the spacing in $\alpha_{0}$ is finer than what we used before. Overall, we adopt a similar scheme as before, but we now require a finer grid centered around the main values of $\alpha_{0}$ we are interested in. Here we use $\alpha_{0}\in\{1.0,\,0.9,\,\cdots,\,0.2\}\times 10^{\text{mag}}$ with $\text{mag}\in\{0,\,-1\}$ as the main grid points, and around each of these grid points we choose the set $\alpha_{0}\in\{+2,+1,-1,-2\}\,\times 10^{\text{mag}-3}$ to give us the other necessary points needed for the finite differences. These choices roughly enforce the same level of accuracy in our results when compared to the grids used for smaller values of $\alpha_{0}$.

Region IV (Blue):

This region of Fig. 1 extends to larger values of $\beta_{0}$ than where spontaneous scalarization occurs, but there are still non-linear effects here that come into play and prevent us from using the quasi-analytic relations. In this region, we only sample in intervals of $\Delta\beta_{0}=0.1$, but we must sample $\alpha_{0}$ on a grid like that used in region I. Thus, region IV is essentially a transition region between spontaneous scalarization and the rest of parameter space, and thus it requires special consideration.

Scalar Charges:

Some representative results for the various scalar charges in MO theory can be found in Fig. 4, for multiple values of $\alpha_{0}$ and $\beta_{0}$ described in the caption. The first three rows in Fig. 4 show the behavior of the scalar charges for $-5\leq\beta_{0}\leq-4.0$, ranging from most red to most blue respectively, and three orders of magnitude in $\alpha_{0}$. Notice that as $\alpha_{0}$ decreases, the growth of $\alpha_{A}$ becomes more rapid as the phase transition of the scalar field becomes more sudden. As a result of this, the peaks in $\beta_{A}$ and $k_{A}$ increase in magnitude and decrease in width because higher order effects become more localized in $m_{A}$. One also notices that the location of $m_{\mbox{\tiny crit}}$, the mass at which spontaneous scalarization turns “on”, moves towards larger masses as $\beta_{0}$ become less negative, as one would expect Damour and Esposito-Farese (1996).

The maximum value of $\beta_{A}$ at $m_{\mbox{\tiny crit}}$ decreases as we decrease $|\beta_{0}|$, while the max values of the $k_{A}$ tends to increase. A possible explanation for this is related to the location of $m_{\mbox{\tiny crit}}$. For the situations where $m_{\mbox{\tiny crit}}$ is larger, the NSs at this mass have a larger moment of inertia, and therefore, it is possible the NS’s inertia is more sensitive to the external scalar field, and hence the increase in $k_{A}$. This reasoning also explains why the right peaks of $k_{A}$ are larger than the left peaks101010The appearance of this second peaks is somewhat unique to MO theory for this range of $\beta_{0}$. As one can see in Fig. 3, $\beta_{A}$ and $k_{A}$ diverge for high masses in DEF theory but not in MO theory. For more negative values of $\beta_{0}$, however, MO theory also exhibits such divergences for large masses..

We have determined numerically that the max values $\beta_{A,\max}$ and $k_{A,\max}$ follow simple linear relations in $\log_{10}$ space for $\alpha_{0}$ and a fixed $\beta_{0}$. The coefficient of these relations remain EOS and theory dependent, but they take the general form

[TABLE]

where $B_{0},\,B_{1},\,K_{0},$ and $K_{1}$ are coefficients that depend on the value of $\beta_{0}$. The quantities $B_{1}$ and $K_{1}$ are always negative for all values of $\beta_{0}$ that we consider, and there is no reason to think that this would be any different for more negative values of $\beta_{0}$. This means that $\beta_{A,\max}$ and $k_{A,\max}$ always increase (decrease) with decreasing (increasing) $\alpha_{0}$ as one might expect. Moreover, for our data we find that $B_{1}$ is always less than 2, and it is typically less than unity. This is important because the combination $\alpha_{0}^{2}\beta_{A}$ appears directly in the PPK parameter $\dot{\omega}$ for binary pulsars with a companion white dwarf. While it is hard to numerically investigate regions of parameter space with $\alpha_{0}<10^{-5}$, these linear relations tell us that the combination $\alpha_{0}^{2}\beta_{A}$ will always decrease with decreasing $\alpha_{0}$ and have negligible effects on the PPK parameters. These relations only hold for $\alpha_{0}<10^{-2}$ but they provide a convenient way to estimate the maximum value of the scalar charges without numerically solving the field equations.

The bottom row of Fig. 4 shows the behavior of the charges for a constant $\beta_{0}=-5$, but multiple orders of magnitude in $\alpha_{0}$, red corresponding to the smallest values and blue to the largest. As one may expect, the overall magnitude of $\alpha_{A}$ is determined by the values of $\beta_{0}$, while its growth rate near the critical mass is determined by $\alpha_{0}$ and is correlated directly with the magnitude of $\beta_{A}$. A similar statement can be made for $k_{A}$ in that $\alpha_{0}$ determines how sudden the growth of the scalar field is, and therefore, it leads to larger values of $k_{A}$ near the critical mass at which these transitions occur.

IV.4.3 No Spontaneous Scalarization: $\beta_{0}>-4.3$

Region V (Cyan):

This region of parameter space is perhaps the simplest and easiest to explore numerically. As mentioned in Sec. IV.1, as long as $\alpha_{0}$ is not too large, there exist scaling relations that we can use to calculate the scalar charges. Reiterating what we explained in that section, in this region of parameter space we can assume the scalar charge $\alpha_{A}$ takes the form

[TABLE]

which then allows us to derive simple relations for the other scalar charges. Taking the derivative of this scalar charge, and applying the chain rule to the definition of $\beta_{A}$, we find

[TABLE]

Equation (59) tells us that, for any value of $\beta_{0}$, $\beta_{A}(m_{A})$ is the same for all values of $\alpha_{0}\lesssim 10^{-2}$. Moreover, this equation also tells us that $\beta_{A}$ is always directly proportional to $\alpha_{A}$, meaning that we technically do not even need to take any derivatives to determine it.

We can find a similar relation to that in Eq. (58) for the inertial charge $k_{A}$, but the derivation is slightly more complicated. Starting with Eqs. (40i)-(40j) and assuming weak fields everywhere such that $e^{\hat{\nu}}\sim 1$ one finds that the moment of inertia becomes

[TABLE]

where we have also neglected terms of order $(Gm/Rc^{2})^{2}$. Using this in the definition for $k_{A}$ in Eq. (27) we find

[TABLE]

where we have made use of the definitions of the other scalar charges $\alpha_{A}$ and $\beta_{A}$. Then, making use of Eqs. (58)-(59) we find that

[TABLE]

where $g$ is a function independent of $\alpha_{0}$.

The scaling relations we have introduced make the exploration of this region in $\beta_{0}$ almost trivial. We simply calculate the three scalar charges for points in parameter space lying on the horizontal dashed lines of Fig. 1 (at $\alpha_{0}=10^{-4}$ and $\alpha_{0}=10^{-3}$), and then, we rescale the results to find the charges for any other point in region V, holding $\beta_{0}$ constant. For the solutions we do calculate directly, we use a grid in $\beta_{0}$ given by $\Delta\beta_{0}=0.01$ and a grid in $\tilde{\rho}_{c}$ given by $\Delta\rho_{c}^{GR}/20$. Because there are no sharp features in the scalar charges in this region of parameter space, we sample considerably less central densities when compared to when spontaneous scalarization occurs.

Region VI (Yellow):

This region presents the same numerical difficulties as region III, expect that there is a lack of spontaneous scalarization in region VI. We use the same grids in $\beta_{0}$ and $\alpha_{0}$ as those used in region III, and the same grid in $\tilde{\rho}_{c}$ as that used in region V. Note, however, the the gray region in Fig. 1 cuts out a significant portion of region VI. As we mentioned earlier, this is because we cannot find stable NS solutions in the gray region, and we must thus omit them from our analysis since they are extremely unlikely to affect binary pulsars constraints.

Scalar Charges:

Some representative results of the scalar charges in region V are shown in Fig. 5. In this figure, we hold $m_{A}$ constant and plot the scalar charges as a function of $\beta_{0}$ spanning the entire region for $\beta_{0}\geq-3.5$. Each curve in Fig. 5 represents a different value of $\alpha_{0}$ and it becomes clear that $\log_{10}\alpha_{A}$ and $\log_{10}k_{A}$ scale directly with $\log_{10}\alpha_{0}$, showing that the relations in Eqs. (58) and (62) are indeed accurate. Moreover, the value of $\beta_{A}$ shown in the middle panel of Fig. 6 shows no dependence on $\alpha_{0}$, verifying that Eq. 59 holds. Figure 6 shows similar scalar charges but as a function of the gravitational mass of the NSs for a representative set of $\beta_{0}$ from region V. One will notice that NSs in this region of parameter space begin to “de-scalarize” as the mass of the NS increases, i.e. $\alpha_{A}$ becomes smaller than $\alpha_{0}$ as $m_{A}$ increases. As one might expect from Fig. 5, all three scalar charges monotonically decrease in magnitude, for all masses, as one increases $\beta_{0}$.

V Using the data file

Now that we have discussed how we calculate the charges and presented some of the results, let us discuss how one can use the end product of this analysis: the data generated for all the scalar charges. This section explains how the master data file is generated, what its properties and limitations are and

V.1 The generation of the master data file

We first set up our numerical grid according to Sec. IV and subsections therein. From each NS solution we extract the boundary conditions in Eqs. (40) and save them to file for post-processing. The previous step requires the bulk of the computational time, as we need to calculate on the order of $10^{5}$ different NS solutions (1 for each combination of $\{\tilde{\rho}_{c},\,\alpha_{0},\,\beta_{0}\}$) for each combination of theory (DEF or MO) and EOS in Sec. III.3.2.

With the full set of data in hand for a theory-EOS combination, we now process it to extract the scalar charges. As we mentioned in Sec. IV, we interpolate the raw data in order to extract information from more masses than we actually sampled. To do this interpolation, however, we need to proceed with caution when dealing with NS that spontaneous scalarize, like those in Fig. 3. One notices that $\alpha_{A}$ “turns over” on itself for large masses. While this feature is not present for every set of NSs that undergoes spontaneous scalarization, it does present a problem for interpolation since the function is multi-valued. To avoid this issue, we remove the data points that lie on the unstable branch of solutions, i.e. the ones that coincidentally make $\alpha_{A}$ double valued (cf. the dashed points in the top left panel of Fig. 3). This is possible because, at least when spontaneous scalarization occurs, there is always NS solutions that reach maximum masses that are at least as large as the maximum mass in GR111111This can be seen from a mass-density curve like that in Fig 2, in which case the scalarized branch of solutions “departs” from the GR curve and eventually “return” for larger values of $\tilde{\rho}_{c}$. Therefore, even if the scalarized branch does not produce a NS with mass greater than the maximum mass in GR, the GR branch will.. With the unstable points of the solution removed we simply continue with the interpolation as described in the previous paragraph.

In total we must interpolate 4 separate functions to give us the data we need for constructing the data files, which are $\{m_{A}(\bar{m}_{A}),\,\bar{m}_{A}(m_{A}),\,\alpha_{A}(\bar{m}_{A}),\,\log I_{A}(\bar{m}_{A})\}$. We need $\alpha_{A}$ and $\log I_{A}$ as function of the baryonic mass for each $\varphi_{\infty}$ in order to take the relevant derivatives in Eqs. (25)-(27) and we need $\{m_{A}(\bar{m}_{A}),\,\bar{m}_{A}(m_{A})\}$ in order to freely switch back and forth between baryonic and gravitational mass121212The gravitational mass is the one appearing in parameterized-post-Keplerian parameters that get constrained from binary pulsar experiments. Therefore, we need the baryonic mass to take the appropriate derivatives for the scalar charges, and the gravitational mass to link our results to binary pulsar experiments.. For the two masses, we interpolate them with a simple linear method to remove any possible artifacts that arise from the interpolation itself. For $\alpha_{A}$ we implement different interpolation schemes depending on if spontaneous scalarization occurs. If spontaneous scalarization is absent, we simply use a cubic spline on $\alpha_{A}$ and this does great since the curves are smooth and generally free of any numerical anomalies. If spontaneous scalarization is present, however, then we use a cubic spline on $\log_{10}\alpha_{A}$ as this helps us better handle the rapid growth of the scalar field, especially when $\alpha_{0}<10^{-4}$. We find that the errors, discussed in Sec. A.1, are significantly smaller when we interpolate $\log_{10}\alpha_{A}$ instead of just $\alpha_{A}$. Lastly, we interpolate $\log I_{A}$ with a cubic spline as well, and while this may introduce some error for low masses, it better suites the data for larger masses, c.f Sec. A.1 for a discussion.

What follows after the interpolation of the raw data is the calculation of the scalar charges $\beta_{A}$ and $k_{A}$ according to Eq. (26) and Eq. (27) respectively, making use of the various finite difference schemes in Eqs. (LABEL:eq:fd_central)-(LABEL:eq:fd_backward). At this point we are able to produce the data we have shown in our figures thus far, and we have the ability to sample our results as finely as we need to in $m_{A}$ in order to produce the most accurate results. However, when producing the master data files for each theory-EOS combination we do not have the luxury of over sampling in $m_{A}$ otherwise each individual data file would be far too large in size. Therefore, we decided to sample in the region $1M_{\odot}<m_{A}<M_{\text{max, GR}}$, with a spacing $\Delta m_{A}=0.002M_{\odot}$, since this is a generous mass range in which we expect to observe pulsars.

V.2 Properties of the master data file

The data files we have generated contain nine columns of data and they have the structure that appears in Table 2. The first three columns of the data correspond to the scalar charges $\alpha_{A}$, $\beta_{A}$, and $k_{A}$ respectively. The fourth and fifth columns contain the values for baryonic mass and gravitational mass respectively, the latter of which falls on the grid described in the previous paragraph. Columns six and seven contain $\alpha_{0}$ and $\log_{10}\alpha_{0}$, of which the former lay on the grid described in section Sec. IV.4.2. Column eight contains the value of $\beta_{0}$, which ranges from $-5\leq\beta_{0}\leq+5$ and lies on a grid with spacing $\Delta\beta_{0}=0.02$ for $\beta_{0}<-4.3$ and $\Delta\beta_{0}=0.1$ for $\beta_{0}>-4.3$.

While not shown in Table 2, we have included 4 more columns of additional information in the master data files that some readers might find useful: the ninth column contains the central density, the tenth contains the radius, the eleventh contains the central value of the scalar field, and the twelfth contains the surface value of the scalar field. However, we point out that in region V not all of these quantities are calculated explicitly because we make use of the aforementioned scaling relations and therefore some assumptions have been made. Because the scalar field is relatively (compared to unity) small in region V, the radius and central density can be assumed to obey the same functional relationship to the gravitational mass as those solutions found on the horizontal black dashed line appearing in Fig. 1, i.e. the set of solutions we apply the scaling relations to. One can see from the boundary condition in Eq. (40g), that the surface value of the scalar field, $\varphi_{s}$, scales directly with $\varphi_{\infty}$, at least to first order in $\varphi_{\infty}$, and thus we have made use of this in the construction of the data files. However, the central value of the scalar field cannot be assumed to obey the same relations because is is one of the free parameters the we must numerically determine through a shooting method. While it would be reasonable that this would also obey the same relations in the limit of weak scalar fields, we could not verify this from our numerical data and have therefore given them a value of “0” in the data files. We emphasize that the true value of the central scalar field is not actual zero, but for the sake of providing a complete data file that can be easily interpolated we have included these null values as a place holder.

The files we have generated can be used in a variety of fields where scalar charges appear, not just binary pulsar experiments, although this is the primary target of this work. The data can be accessed though a git repository https://github.com/eXtremeGravityInstitute/scalar˙charges and is structured on a uniform grid such that one can either use the data as is, or interpolate it if desired. For interpolation purposes, we have also included three separate files containing the grid points in $\{m_{A},\,\alpha_{0},\,\beta_{0}\}$ that we used and are labeled appropriately in the git repository. We have included both a Python and Mathematica script that is capable of reading in the data, setting up the numerical grid, and interpolating the data using SciPy’s RegularGridInterpolator function for Python and Mathematica’s native Interpolation function, both of which make use of a linear order method for three-dimensional data.

Recall that for the gray region of parameter space in Fig. 1 we were not able to find NS solutions that were of any use. As we have mentioned earlier, excluding this data from our investigation is not necessarily a shortcoming since these regions of parameter space are so heavily constrained by solar system observations and they are generally excluded from any analysis. However, for the sake of enforcing that all data lie on a uniform grid and can therefore be interpolated with ease, we found it beneficial to include this region of parameter space in our data. To do this, we artificially give values to the three scalar charges as a sort of place holder in the data. Since these regions of parameter space are generally excluded by observations we have given the scalar charges all an unphysical value of $10^{5}$. We chose this value more out of convenience in an effort to force, say, an MCMC to avoid these points in parameter space because NSs do not exist and therefore the charges technically do not exist either.

V.3 Comparison to full data

One questions we need to concern ourselves with is whether or not the data in these files accurately reproduces the raw data from which they came since we have had to, in some cases, undersample the data in $m_{A}$. Here we discuss how the data recovers the results presented in Sec. IV, which were produced with full numerical data. The points in parameter space we consider here lie on our numerical grid and therefore should show great agreement with the full data, provided we sampled the charges finely enough in $m_{A}$. A more detailed discussion of the limitations of our data can be found in Appendix B where we point out some known issues.

Figure 7 shows a comparison between some of the full data found in Fig. 4, i.e. points labeled 2-5 in Table 3, and the data file for MO theory and AP3 EOS. As one can see, there is remarkable consistency between the two data sets, even in the most non-linear regions of parameter space, e.g. $\beta_{0}=-5$ and $\alpha_{0}=10^{-5}$. In fact, the only real error that is noticeable is in $\beta_{A}$ for these parameters, and only near the sharpest part of the peak does the data files deviate from the full data

Even deviations as large as 10% should not have a significant effect on binary pulsar constraints, and the reason is threefold. First, the peaks in $\beta_{A}$ are extremely isolated in $m_{A}$ and there may not be a pulsar with that particular mass. Second, even if this value of $m_{A}$ were important, the value of $\beta_{A}$ is still so large that it would be immediately ruled out. Third, this error only occurs in the regions of parameter space that are already tightly constrained Damour and Esposito-Farese (1996); Freire et al. (2012); Anderson et al. (2016); Anderson and Yunes (2017). The astute reader may point out that the location of the peaks in $\beta_{A}$ depend on the EOS and $\beta_{0}$, which would make it quite possible for an observed pulsar to lie on one of these peaks at some points in the $\{\alpha_{0},\,\beta_{0}\}$ plane. While this is true, our second point above still holds and, in fact, the differences between the tabulated and full data tends to decrease as $\beta_{0}$ becomes less negative and $\alpha_{0}$ becomes larger.

VI Conclusion

We have extended the original work in Ref. Damour and Esposito-Farese (1996) and calculated the scalar charges $\{\alpha_{A},\,\beta_{A},\,k_{A}\}$ for a large region of the $\{\alpha_{0},\,\beta_{0}\}$ parameter space and 11 physical equations of state, in two distinct scalar tensor theories (that of Damour-Esposito-Farése and Mendes-Ortiz). We have presented the numerical schemes we implemented to complete these calculations and presented our results. Our goal was to calculate the scalar charges and tabulate them so that they could be of use in the future, particularly in the application of binary pulsar tests of gravity. We have investigated both the error of our numerical solutions, as well as the ability of certain scaling relations to reproduce full numerical results. Through this paper, the data is made fully available to the community.

Future work that utilizes this data include tests of STTs with binary pulsars and gravitational waves. In particular, this data makes it possible to perform a Bayesian analysis on the PPK parameters of binary pulsar system. Such analysis would require the use of an MCMC in which case one would need knowledge of the scalar charges in order to compute the likelihood. Instead of calculating the scalar charges on the fly, which we have demonstrated to extremely computationally expensive, one can make use of the data file provided in this paper to significantly speed up these likelihood evaluations. We intend to perform such an investigation in future work.

LIGO and VIRGO will inevitably detect more NS mergers in the future, some of which are likely to be NS-BH systems. Such systems are ideal for testing STTs because BHs do not develop scalar charges in these theories and therefore the emission of dipolar GWs would be maximized in these scenarios. The data provided in this paper would again be necessary for one to perform a Bayesian analysis of the data through MCMC simulations. Furthermore, if one wishes to study the constraints the future GW detectors can place on STTs, high resolution calculations of the scalar charges would be a crucial for such an analysis.

Acknowledgements.

We would like to thank Norbert Wex, Hector O’ Silva, and Travis Robson for useful insight and discussions. DA would like to thank Paulo Freire and the Max-Planck-Institut für Radioastronomie for their hospitality during part of this work. We would also like to acknowledge the support of the Research Group at Montana State University through their High Performance Computer Cluster Hyalite. NY and DA acknowledge support from NSF grant PHY-1759615 and NASA grants NNX16AB98G and 80NSSC17M0041.

Appendix A Error Analysis of Results

Here we attempt to quantify the numerical error in our results. We start with an investigation of the error associated with our grid spacing in $\tilde{\rho}_{c}$ and show that we can trust our results to within 1% in the worst of cases for the regions of parameter space that we have explored. We then look into how reliable the scaling relations are and how much error we introduce by using these rather than populating the entire parameter space fully numerically. We then finish with a discussion of the effects of using piecewise polytrope versus their tabulated counterparts. For our analysis, we focus on the numbered points appearing in Fig. 1, which we have detailed below in Table 3. For each of the analyses we perform, we investigate the error at these points in parameter space in detail in an attempt to quantify our errors across the entire parameter space.

A.1 Grid in central density

In order to obtain the most accurate and precise results from our numerical calculations one would have to use an infinitely dense grid in central density $\tilde{\rho}_{c}$ such that interpolating between points in Fig. 2 leaves no room for error. However, this is not computationally feasible, and thus, by reducing the number of points in this grid we introduce numerical error that is not associated with our methods of solving the field equations. The grid we have used in the end, for each individual region of parameter space, was chosen to achieve sufficient confidence in our results for the amount of computation time needed for the calculations. While we are confident in the grids we have chosen, there is still numerical error associated with our choices and we quantify those here.

In order to assess our errors we decided to double the number of points in our central density grid and calculate the relative errors between these solutions and the ones we have calculated using our original grids, which we have done for every point in parameter space detailed in Table 3. While we analyzed each of these points in detail, we only discuss the results for the points with the worst errors and give explanation for why the errors arise. For concreteness, to calculate relative percent error for any of the scalar charges, denoted by “$\chi$” here, we use

[TABLE]

where $\chi_{2}$ is the solution with twice as many grid points as $\chi_{1}$. We evaluate this for all values of $m_{A}$ for the specific combination of $\alpha_{0}$ and $\beta_{0}$ and quote the largest values of relative error in the columns of Table 3 with $\tilde{\rho}_{c}$ in parentheses.

Figure 8 shows a representative sample of our errors, including points labeled 2, 8, 10, and 13 as they show special features worth discussing. The first row in Fig. 8, showing the errors in the charges for point 2 in Table 3, represents what we consider to be one of the most error prone regions of parameter space. Spontaneous scalarization occurs here, and since $\alpha_{0}$ is so small the phase transition in the scalar field is extremely sharp. However, our results for $\alpha_{A}$ are good to withing 0.01%, with some of the worst error appearing exactly when the phase transition occurs ($\sim 1.3M_{\odot}$). The errors in $\beta_{A}$ are slightly worse as one might expect since a numerical derivative is involved, but surprisingly enough the largest error is not associated with the phase transition. The spike in the error at $\sim 1.4M_{\odot}$ is related to a mishap in the interpolation of the finer grid data, and while this does not happen often, our solutions are sometimes prone to this type of error in this part of parameter space.

The error in $k_{A}$ may seem alarming at first for low masses, but this only occurs because $k_{A}\sim 0$ here and even tiny numerical noise in the results can generate a large relative error in the solution. This error is not a consequence of our numerical grid, but rather is an error associated with a the fact that numerical integration has finite precision and how we interpolate our results. For comparison we have also included the absolute error for $k_{A}$ plotted by a dashed red line in Fig. 8. An obvious downfall here might be that we should have used a finer grid for the low masses, which would have most certainly increases the accuracy of our interpolations somewhat. The other issue, however, could be the method of interpolation we used. As described in Sec. V, when calculating $k_{A}$ we must interpolate $\log I_{A}(\bar{m}_{A})$ for multiple values of $\alpha_{0}$ in order to calculate derivatives, and we use a cubic spline to do so. Using a cubic spline on potentially noisy data like this is a good way to introduce extra error, which is what we are seeing here. However, switching to a linear order interpolation method significantly increases our errors for larger masses, precisely in the more interesting regions of parameter space where pulsar masses tend to lie. For this reason, we sacrifice precision on the lower end of masses in order to increase it elsewhere for more relevant masses.

The second, third, and fourth rows of Fig. 8 tell a different story than the first. For all these points in parameter space we find great agreement between solutions and therefore very small relative error. As expected, error in $\alpha_{A}$ are extremely small and this is due to the fact that there is no spontaneous scalarization and it is easy to extract $\alpha_{A}$ from the NS solutions. The errors in $\beta_{A}$ are slightly worse than those of $\alpha_{A}$ but still show exceptional agreement. As we saw with the with point 1, the errors in $k_{A}$ tend to be much worse than those of the other scalar charges.

A.2 Analytic Scaling

One of the greatest properties about the cyan region of parameter space in Fig. 1 is that we can make use of the scaling relations presented in Sec. IV.4.3 to significantly reduce the number of NS solutions needed to explore this region. There is, however, some error associated with these scaling relations since we are, effectively, ignoring some of the non-linearities that appear in the field equations. While we expect these relations to hold in the small $\alpha_{0}$ regime, we need to show numerically that this is indeed the case. To investigate this particular kind of error we have decided to explore the cyan region of parameter space with the same numerical grid for the blue region, described in Sec. IV.4.2, but only for AP3 and SLy4 EOS to get a sense of the error.

To make use of the scaling relations, we simply calculate all the scalar charges on the dashed horizontal lines in Fig. 1 and use Eq. (58) to find the function $f(\beta_{0},m_{A})$. Once we have found $f(\beta_{0},m_{A})$, then $\alpha_{A}$ can be solved for all other $\alpha_{0}<10^{-2}$ by substituting $f$ and the new $\alpha_{0}$ back into Eq. (58). Since we have already calculated the actual, unscaled scalar charges for the AP3 and SLy4 EOSs can can compare the scaled versions to the full numerically solved ones.

Figure 9 shows some of the results from our error analysis for $\beta_{0}=\pm 3$, where the scaling relations are expected to hold. We see that the relations have the most error for solutions with $\alpha_{0}=10^{-2}$ (blue curves) but are never more than 1%. For all values of $\alpha_{0}<10^{-2}$ the error continues to decrease and would presumably go to zero in the limit that $\alpha_{0}\rightarrow 0$ if we had infinite precision in our numerics. In general, we find errors of the same order of magnitude as the ones presented in Fig. 9 across the entire cyan region of parameter in Fig. 1, with the errors being slightly larger in MO theory than in DEF theory. One notices that the errors for $\alpha_{A}$ and $\beta_{A}$ look very similar, and they should be this way according to Eq. (59) since $\beta_{A}$ is a derivative of $\alpha_{A}$. The inertial charge $k_{A}$ also has a similar structure to the other charges but is polluted with more numerical noise, which is extremely evident for $\beta_{0}=+3$. The distinct difference in the function form of the errors between $\beta_{0}=-3$ and $\beta_{0}=+3$ can be attributed to the fact that the effects of the non-linearities in the scalar field are more prominent for negative values of $\beta_{0}$.

A.3 Piecewise Polytropes

The piecewise polytropic approximation to tabulated EOSs in Ref. Read et al. (2009) allows one to considerable speed up numerical calculations involving NSs because of their analytic nature. The data files we provide were all generated using the full tabulated data, but for others who wish to perform future calculation, it would be nice to have an idea how these approximations affect the final results. In Fig. 3 we briefly compared how charges calculated using the piecewise polytropes compare to those calculated with the full tabulated EOS, and one can see that aside from slight apparent shifts in the masses, the curves seem nearly identical. We have verified, using the points in Table 3 as a representative sample, that this is indeed consistent across the entire parameter space. The similarity between the charges should not come as a surprise considering how well the piecewise polytropes match the actual data, c.f Table III in Ref. Read et al. (2009). For all of the equations of state we consider here, the residuals in Ref. Read et al. (2009) for both the mass and moment of inertia of the NSs are less than 2%, which are the two most important quantities when calculating the scalar charges. Having such small deviations between tabulated and approximated EOSs leads to very small deviations when calculating $\alpha_{A}$, $\beta_{A}$, and $k_{A}$ from Eqs. (25)-(27).

Appendix B Limitations of Data Files

We have discovered a few limitations of the data we have provide and these are discussed in this appendix. We should point out, however, that these complications arise when one wishes to interpolate the data files we have generated, in order to determine that charges for points that do not lie on our numerical grid. As far as we are able to tell, the data files are able to reproduce the full numerical data used to make the files to great level of accuracy (see Sec. V) for all points that lie on our grid.

Figure 10 shows how the interpolation of the data behaves for points that lie on (solid) and off (dashed) numerical grid we have established for $\alpha_{0}$ and $\beta_{0}$. One notices that the solid curves appear as expected, according to the results presented in Figs. 3 and 4; $\alpha_{A}$ experiences a smooth rapid growth while $\beta_{A}$ and $k_{A}$ have peaks near $m_{\mbox{\tiny crit}}$ for which spontaneous scalarization turns “on” and “off”. However, the dashed curves, which lie off the numerical grid, have somewhat significantly different features. The curves for $\alpha_{A}$ appear to be good, but one notices that the dashed curves develop slight instantaneous “discontinuities” in the slope during the growth of the scalar field. As a result of this apparent discontinuities, $\beta_{A}$ develops a double peak near the critical mass at which the phase transition occurs. Likewise, a very similar feature develops in $k_{A}$ in which the single peak that is present in the solid curves turns into two peaks that are not as large in magnitude.

The issue we are seeing in Fig. 10 is a direct consequence of the nature of the scalar charges, particularly the presence of the phase transition, and trying to interpolate them. Consider $\beta_{A}$ for example, in which case we expect the magnitude of the peak at the critical mass to decrease as $\beta_{0}$ becomes less negative, cf. the plots of $\beta_{A}$ in Fig. 4. We also expect the critical mass to shift to higher masses as we allow $\beta_{0}$ to become less negative. Therefore, there are multiple features changing in the solutions with the variation of just a single parameter, in this case $\beta_{0}$. Because of this dependence, the linear interpolator has issues when it tries to interpolate between these peaks because the algorithm considers changes in $\beta_{0}$ and $m_{A}$ at the same time, which physically are in effect co-dependent on each other in a non-linear way. As it might be expected, it is hard to interpolate any function with extremely sharp features like we have here and taking the log of the data does not seem to improve anything in this situation.

The features described thus far are most prominent in the $\alpha_{0}\leq 10^{-5}$ case since the peaks are the most narrow here. However, this artifact of the interpolation does arise for other values of $\alpha_{0}$ when spontaneous scalarization, but since the effects are less localized, i.e. the phase transition is more smooth for larger values of $\alpha_{0}$, the discrepancy is less severe. Once out of the regime of spontaneous scalarization these artifacts no longer appear to be present and the interpolation of the data does an excellent job of giving us information between our grid points.