Reduced-order surrogate models for scalar-tensor gravity in the strong field and applications to binary pulsars and GW170817

TL;DR

This paper develops fast, accurate surrogate models for scalar-tensor gravity in neutron stars, enabling efficient analysis of pulsar and gravitational wave data to test strong-field gravity theories.

Contribution

The authors create reduced-order surrogate models for scalar-tensor gravity predictions, significantly speeding up calculations while maintaining high accuracy, aiding in gravitational tests.

Findings

Models achieve ~1% accuracy

Speed up calculations by 100x

Constrain gravity parameters using pulsar and GW data

Abstract

We investigate the scalar-tensor gravity of Damour and Esposito-Far\`ese (DEF), which predicts non-trivial phenomena in the nonperturbative strong-field regime for neutron stars (NSs). Instead of solving the modified Tolman-Oppenheimer-Volkoff equations, we construct reduced-order surrogate models, coded in the pySTGROM package, to predict the relations of a NS radius, mass, and effective scalar coupling to its central density. Our models are accurate at $\sim1\%$ level and speed up large-scale calculations by two orders of magnitude. As an application, we use pySTGROM and Markov-chain Monte Carlo techniques to constrain parameters in the DEF theory, with five well-timed binary pulsars, the binary NS (BNS) inspiral GW170817, and a hypothetical BNS inspiral in the Advanced LIGO and next-generation GW detectors. In the future, as more binary pulsars and BNS mergers are detected, our surrogate models will be helpful in constraining strong-field gravity with essential speed and accuracy.

Tables

Table 1. Table 1: Binary parameters of the five NS-WD systems (PSRs J0348+0432  Antoniadis  et al. ( 2013 ) , J1012+5307  Lazaridis  et al. ( 2009 ); Antoniadis  et al. ( 2016 ); Desvignes  et al. ( 2016 ) , J1738+0333  Freire  et al. ( 2012 ) , J1909 − - 3744  Reardon  et al. ( 2016 ); Desvignes  et al. ( 2016 ); Arzoumanian  et al. ( 2018 ) and J2222 − - 0137  Cognard  et al. ( 2017 ) ) that we use to constrain the DEF theory. The observed time derivatives of the orbit period, P ˙ b obs superscript subscript ˙ 𝑃 𝑏 obs \dot{P}_{b}^{\rm obs} , will be corrected with the latest Galactic model in Ref.  McMillan ( 2017 ) . The intrinsic derivatives of the orbital period, P ˙ b int superscript subscript ˙ 𝑃 𝑏 int \dot{P}_{b}^{\rm int} , are obtained from P ˙ b obs superscript subscript ˙ 𝑃 𝑏 obs \dot{P}_{b}^{\rm obs} , by subtracting the other kinematic effects, such as the acceleration effect  Damour and Taylor ( 1991 ) and the “Shklovskii” effect  Shklovskii ( 1970 ) . For PSRs J0348+0432, J1012+5307, and J1738+0333, the pulsar masses, m p obs superscript subscript 𝑚 𝑝 obs m_{p}^{\rm obs} , are not derived from ratio observations. They are obtained from the companion masses, m c obs superscript subscript 𝑚 𝑐 obs m_{c}^{\rm obs} , and the mass ratios, q 𝑞 q , without assuming the validity of GR. Similar consideration was made for PSRs J1909 − - 3744 and J2222 − - 0137  Shao  et al. ( 2017 ) . The standard 1- σ 𝜎 \sigma errors in the least significant digit(s) are shown in parentheses.

Pulsar	J0348+0432 Antoniadis et al. (2013)	J1012+5307 Lazaridis et al. (2009); Antoniadis et al. (2016); Desvignes et al. (2016)	J1738+0333 Freire et al. (2012)	J1909$-$3744 Reardon et al. (2016); Desvignes et al. (2016); Arzoumanian et al. (2018)	J2222$-$0137 Cognard et al. (2017)
Orbital period, $P_b$ (d)	0.102424062722(7)	0.60467271355(3)	0.3547907398724(13)	1.533449474406(13)	2.44576454(18)
Eccentricity, $e$	$2.6\left(9\right)\times {10}^{-6}$	$1.2\left(3\right)\times {10}^{-6}$	$3.4\left(11\right)\times {10}^{-7}$	$1.14\left(10\right)\times {10}^{-7}$	0.00038096(4)
Observed ${\dot{P}}_b$, ${\dot{P}}_b^{\mathrm{obs}}$ ($\mathrm{fs}{s}^{-1}$)	$-273\left(45\right)$	$-50\left(14\right)$	$-17.0\left(31\right)$	$-503\left(6\right)$	$200\left(90\right)$
Intrinsic ${\dot{P}}_b$, ${\dot{P}}_b^{\mathrm{int}}$ ($\mathrm{fs}{s}^{-1}$)	$-274\left(45\right)$	$-5\left(9\right)$	$-27.72\left(64\right)$	$-6\left(15\right)$	$-60\left(90\right)$
Mass ratio, $q\equiv m_p/m_c$	11.70(13)	10.5(5)	8.1(2)	…	…
Pulsar mass, $m_p^{\mathrm{obs}}$ ($M_{\odot }$)	…	…	…	1.48(3)	1.76(6)
Companion mass, $m_c^{\mathrm{obs}}$ ($M_{\odot }$)	${0.1715}_{-0.0030}^{+0.0045}$	0.174(7)	${0.1817}_{-0.0054}^{+0.0073}$	0.208(2)	$1.293\left(25\right)$

Table 2. Table 2: Properties of GW170817 from the LIGO/Virgo Collaboration  Abbott  et al. ( 2017b , 2019c , 2018b , 2019b ) . The primary mass m 1 subscript 𝑚 1 m_{1} and secondary mass m 2 subscript 𝑚 2 m_{2} are determined with the low-spin prior assumption  Abbott  et al. ( 2019c ) . The radii of NSs, R 1 subscript 𝑅 1 R_{1} and R 2 subscript 𝑅 2 R_{2} , are derived with the EOS-insensitive relations  Yagi and Yunes ( 2017 ); Chatziioannou  et al. ( 2018 ) .

GW170817 Abbott et al. (2017b)	$m_1^{\mathrm{obs}}\left(M_{\odot }\right)$ Abbott et al. (2019c)	$m_2^{\mathrm{obs}}\left(M_{\odot }\right)$ Abbott et al. (2019c)	$R_1^{\mathrm{obs}}\left(\mathrm{km}\right)$ Abbott et al. (2018b)	$R_2^{\mathrm{obs}}\left(\mathrm{km}\right)$ Abbott et al. (2018b)
90% CL	$(1.36,1.60)$	$(1.16,1.36)$	$(9.1,12.8)$	$(9.2,12.8)$
68% CL	$(1.41,1.55)$	$(1.20,1.32)$	$(9.7,11.9)$	$(9.6,11.8)$

Table 3. Table 3: The different scenarios used in the paper to constrain the DEF theory. Different combinations of five binary pulsars (see  Table   1 ), GW170817 (see  Table   2 ), and a hypothetical BNS BNS@200Mpc , are investigated. Their corresponding log-likelihood functions, ln ⁡ ℒ PSR subscript ℒ PSR \ln\mathcal{L}_{\rm PSR} , ln ⁡ ℒ GW subscript ℒ GW \ln\mathcal{L}_{\rm GW} , and ln ⁡ ℒ BNS subscript ℒ BNS \ln\mathcal{L}_{\rm BNS} , are expressed in  Eqs.   20 , 19 and  18 . The scenarios (i) to (iii) correspond to real observations, while the scenarios (iv) to (vi) involve a hypothetical BNS@200Mpc to be observed by the Advanced LIGO, CE, and ET. Notice that the scenario (ii) was investigated in great detail in Ref.  Shao  et al. ( 2017 ) , while the other five scenarios include new extensions.

	Scenario	Log-likelihood function
(i)	GW170817	$\ln {ℒ}_{\mathrm{GW}}$	GW170817 only
(ii)	PSRs	$\ln {ℒ}_{\mathrm{PSR}}$	five pulsars
(iii)	PSRs+GW170817	$\ln {ℒ}_{\mathrm{PSR}}+\ln {ℒ}_{\mathrm{GW}}$	combining five pulsars and GW170817
(iv)	PSRs+aLIGO	$\ln {ℒ}_{\mathrm{PSR}}+\ln {ℒ}_{\mathrm{BNS},\mathrm{aLIGO}}$	combining five pulsars and a BNS@200Mpc observed by the Advanced LIGO
(v)	PSRs+CE	$\ln {ℒ}_{\mathrm{PSR}}+\ln {ℒ}_{\mathrm{BNS},\mathrm{CE}}$	combining five pulsars and a BNS@200Mpc observed by CE
(vi)	PSRs+ET	$\ln {ℒ}_{\mathrm{PSR}}+\ln {ℒ}_{\mathrm{BNS},\mathrm{ET}}$	combining five pulsars and a BNS@200Mpc observed by ET

Table 4. Table 4: The first column lists three ground-based GW detectors. The second column gives their frequencies at low end. The third and fourth columns are respectively the expected SNRs and the expected 1- σ 𝜎 \sigma uncertainties in the difference of the effective scalar couplings from the Fisher-matrix analysis  Shao  et al. ( 2017 ) , for a hypothetical event, BNS@200Mpc .

Detector	$f_{\mathrm{in}}$ (Hz)	$\rho$	$\sigma \left(\left|\Delta \alpha \right|\right)$
aLIGO	10	10.6	$8.8\times {10}^{-3}$
CE	5	450	$7.9\times {10}^{-4}$
ET	1	153	$4.0\times {10}^{-4}$

Equations

$$S$$

$$+S_m\left[{\psi }_m;A^2(\phi )g_{\mu \nu }^{\ast }\right].$$

$$R_{\mu \nu }^{\ast }$$

$${\square }_{g^{\ast }}\phi$$

$$\alpha (\phi )\equiv \frac{\partial \ln A(\phi )}{\partial \phi }.$$

$$\ln A(\phi )=\frac{1}{2}{\beta }_0{\phi }^2,$$

$${\beta }_0\equiv {\frac{{\partial }^2\ln A(\phi )}{\partial {\phi }^2}|}_{\phi ={\phi }_0}\lesssim -4.$$

$${\alpha }_A\equiv {\frac{\partial \ln m_A(\phi )}{\partial \phi }|}_{\phi ={\phi }_0},$$

$$V_{\mathrm{int}}=-\frac{G_{\ast }{m}_A{m}_B}{r_{AB}}\left(1+{\alpha }_A{\alpha }_B\right).$$

$${\dot{P}}_b^{\mathrm{dipole}}=-\frac{2\pi G_{\ast }}{c^3}g(e)\left(\frac{2\pi }{P_b}\right)\frac{m_p{m}_c}{m_p+m_c}{\left({\alpha }_p-{\alpha }_c\right)}^2,$$

$${\dot{P}}_b^{\mathrm{quad}}=-\frac{192\pi G_{\ast }^{5/3}}{5c^5}f(e){\left(\frac{2\pi }{P_b}\right)}^{5/3}\frac{m_p{m}_c}{{\left(m_p+m_c\right)}^{1/3}},$$

$$g(e)$$

$$f(e)$$

$${\sigma }_i\equiv \underset{h\in \mathbf{V}}{\max }{\Vert h(\cdot ;\mathbf{\lambda })-{\mathcal{P}}_{i}h(\cdot ;\mathbf{\lambda })\Vert }^2,$$

$$h(t;\mathbf{\lambda })\approx \sum _{i=1}^m{c}_i(\mathbf{\lambda })e_i(t)\approx \sum _{i=1}^m⟨h(\cdot ;\mathbf{\lambda }),e_i(\cdot )⟩e_i(t).$$

$$\mathbf{V}\equiv \mathcal{A}\times \mathcal{B}.$$

$$\mathcal{E}(X)=\left|\frac{X_{\mathrm{ROM}}-X_{\mathrm{mTOV}}}{X_{\mathrm{ROM}}+X_{\mathrm{mTOV}}}\right|,X\in \{m_A,R,{\alpha }_A\},$$

$$P\left({\alpha }_0,{\beta }_0|\mathcal{D},\mathcal{H},\mathcal{I}\right)=$$

$$\int \frac{P\left(\mathcal{D}|{\alpha }_0,{\beta }_0,\mathbf{\Xi },\mathcal{H},\mathcal{I}\right)P\left({\alpha }_0,{\beta }_0,\mathbf{\Xi }|\mathcal{H},\mathcal{I}\right)}{P\left(\mathcal{D}|\mathcal{H},\mathcal{I}\right)}\mathrm{d}\mathbf{\Xi },$$

$$\ln {\mathcal{L}}_{\mathrm{PSR}}=$$

$$-\frac{1}{2}\sum _{i=1}^{N_{\mathrm{PSR}}}{\left[{\left(\frac{{\dot{P}}_b^{\mathrm{th}}-{\dot{P}}_b^{\mathrm{int}}}{{\sigma }_{{\dot{P}}_b^{\mathrm{int}}}}\right)}^2+{\left(\frac{m_p-m_p^{\mathrm{obs}}}{{\sigma }_{m_p^{\mathrm{obs}}}}\right)}^2+{\left(\frac{m_c-m_c^{\mathrm{obs}}}{{\sigma }_{m_c^{\mathrm{obs}}}}\right)}^2\right]}_i,$$

$$\ln {\mathcal{L}}_{\mathrm{GW}}$$

$$+\sum _{i=1}^{N_{\mathrm{GW}}}{\mathtt{D}}_i\left(|{\alpha }_1-{\alpha }_2|,|\Delta \alpha {|}^{\mathrm{upper}}\right).$$

$$\ln {\mathcal{L}}_{\mathrm{BNS}}=-\frac{1}{2}\sum _{i=1}^{N_{\mathrm{GW}}^{\star }}{\left[{\left(\frac{{\alpha }_1-{\alpha }_2}{\sigma (|\Delta \alpha |)}\right)}^2\right]}_i.$$

$$\left({\tilde{h}}_1(f)\right|{\tilde{h}}_2(f))\equiv 2{\int }_{f_{\min }}^{f_{\max }}\frac{{\tilde{h}}_1^{\ast }(f){\tilde{h}}_2(f)+{\tilde{h}}_1(f){\tilde{h}}_2^{\ast }(f)}{S_n(f)}\mathrm{d}f.$$

$${\mathcal{I}}_{ij}=\left(\frac{\partial \tilde{h}}{\partial {\xi }_i}|\frac{\partial \tilde{h}}{\partial {\xi }_j}\right).$$

Full text

I Introduction

The theory of general relativity (GR), proposed by Albert Einstein in 1915, postulates that gravity is mediated only by a long-range, spin-2 tensor field, $g_{\mu\nu}$ Einstein (1915). For more than a century, tests of GR have never stopped. The bulk of accurate experimental tests from, (i) the Solar System Will (2014), (ii) the high-precision timing of binary pulsars Stairs (2003); Wex (2014); Shao and Wex (2016), and (iii) the gravitational-wave (GW) observations of coalescing binary black holes (BBHs) Abbott et al. (2016, 2017a, 2019a) and binary neutron stars (BNSs) Abbott et al. (2017b, c, 2019b), have all been proven to be in line with GR.

However, there are various motivations to look for theories beyond GR Will (2018). As one of the most natural alternatives, the scalar-tensor gravity, in addition to $g_{\mu\nu}$, adds a long-range, spin-0 scalar field, $\varphi$. This theory was conceived originally by Scherrer Goenner (2012) and Jordan Jordan (1949). Viewpoints similar to the one in the scalar-tensor gravity have also been sketched in the Kaluza–Klein theory, the string theory, and the brane theory Fujii and Maeda (2007). The extra scalar degree of freedom is potentially related to the dark energy, the inflation, and a possible unified theory of quantum gravity.

From 1960s to the present time, a healthy version of scalar-tensor gravity has played the most influential role (see Refs. Fujii and Maeda (2007); Will (2018) for reviews). We now call it the Jordan–Fierz–Brans–Dicke (JFBD) theory Jordan (1949, 1959); Fierz (1956); Brans and Dicke (1961). Inspired by the Mach’s principle, the gravitational constant $G$ is promoted to a time-varying dynamical field in the JFBD theory Brans and Dicke (1961). With the additional scalar field coupled nonminimally to the Einstein-Hilbert Lagrangian, the JFBD theory leads to a violation of the strong equivalence principle (SEP) Eardley (1975); Shao and Wex (2016). In this paper, we focus on a class of special mono-scalar-tensor gravity, formulated by Damour and Esposito-Farèse (DEF) Damour and Esposito-Farèse (1992, 1993, 1996). Relative to the JFBD theory, the DEF theory extends the conformal coupling function by including a quadratic term, which dictates the way matter couples to the scalar field. Within a certain parameter space, it significantly modifies the level where the SEP is violated for strongly self-gravitating NSs Damour and Esposito-Farèse (1993).

The theory of scalar-tensor gravity has been extensively investigated in the weak field, mainly from experiments in the Solar System Will (2014). The most stringent constraint comes from the Cassini probe Bertotti et al. (2003). In the parametrized post-Newtonian (PPN) framework Will (2014), it is verified to a high precision $\sim 10^{-5}$ that the DEF theory is very close to GR in the weak-field regime Damour (2007). In the nonperturbative strong-field regime, Damour and Esposito-Farèse (1993) noticed a sudden strong activation of the scalar field for NSs. This kind of intriguing feature, where the DEF theory can be reduced to GR in the weak field while being significantly different in the strong field, has caused enormous interests Freire et al. (2012); Antoniadis et al. (2013); Shao et al. (2017); Anderson et al. (2019). The process that causes significant differences from GR is referred as strong-field scalarization. It has been comprehensively investigated for decades in terms of spontaneous scalarization Damour and Esposito-Farèse (1993, 1996); Esposito-Farèse (2004); Anderson et al. (2016, 2019); Anderson and Yunes (2019), as well as induced and dynamical scalarizations Barausse et al. (2013); Shibata et al. (2014); Sennett and Buonanno (2016); Sennett et al. (2017); Shao et al. (2017).

Different kinds of scalarization phenomena are defined in the following ways.

Spontaneous scalarization occurs in an isolated and compact star whose compactness exceeds a critical value Damour and Esposito-Farèse (1993); Esposito-Farèse (2004). This phenomenon can be regarded as a phase transition Damour and Esposito-Farèse (1996); Sennett et al. (2017). 2. 2.

Induced and dynamical scalarizations are discovered in numerical-relativity (NR) simulations of BNSs in the DEF theory. They hasten the plunge and merger phases relative to GR Barausse et al. (2013). The induced scalarization occurs in those binary systems, where one component of the binary has already spontaneously scalarized while the other has not in the early inspiral Barausse et al. (2013). 3. 3.

Dynamical scalarization corresponds to the binary systems, where none of the binary components can spontaneously scalarize in isolation, but they get scalarized when the gravitational binding energy of the orbit exceeds a critical value Barausse et al. (2013).

In this paper, we pay particular attention to the strong-field region, where the spontaneous-scalarization phenomena are significant. Compared with an unscalarized system, the scalarized binary system brings the following manifest changes: (i) an additional gravitational binding energy for the orbit, (ii) an enhancement in the decay rate of binary’s orbital period and the energy flux by an extra dipolar radiation Damour (2007). As we know, dipolar radiation corresponds to the $-1$ post-Newtonian (PN) correction.111We refer $n\,$PN to $\mathcal{O}\left(v^{2n}/c^{2n}\right)$ corrections with respect to the Newtonian order, where $v$ is the characteristic relative speed in the binary. Here, we follow the convention in the GW phase, where the quadrupolar radiation reaction is denoted as $0\,$PN. Therefore, the dipolar radiation is at $-1\,$PN order. It enters at a lower order relative to the typical quadrupolar radiation in GR. This means that, the smaller the relative speed, the greater the relative radiating effect of the dipolar radiation.

With the dominant radiating component being the dipolar emission at early time, a binary system emits extra energy in addition to GR. In certain binary pulsar systems, it is a powerful means to probe the strength of dipolar contribution that can be caused by the spontaneous scalarization Damour and Esposito-Farèse (1996); Wex (2014). With a small characteristic speed $\sim 10^{-3}c$, the gravitational radiation of binary pulsar systems could be dominated by the dipolar emission.

A pulsar emits ratio signals like a lighthouse. In a binary pulsar system, the spin period of the recycled pulsar is usually extremely stable. For several years to decades, those periodic signals have been continuously monitored by large radio telescopes on the Earth. It makes pulsar a clock that can rival the best clocks for precision fundamental physics Lorimer and Kramer (2005); Guo et al. (2018). The accurate measurement technique, the so-called pulsar timing, models the times of arrival (TOAs) of pulses emitted from the pulsar and determines timing parameters to a high precision Stairs (2003); Lorimer and Kramer (2005).

In order to accommodate alternative gravity theories, the parametrized post-Keplerian (PPK) formalism was developed as a generic pulsar timing model Damour and Taylor (1992). A set of theory-independent Keplerian and post-Keplerian timing parameters are determined with high precision in a fit of the timing model to the TOAs Damour (2007). We can obtain extremely precise physical parameters to describe those systems. They can be used to place constraints on alternative gravity theories Stairs (2003); Wex (2014); Shao and Wex (2016). Binary pulsar is currently one of the best available strong-field testbeds for testing gravity Damour (2007).

Recently, GWs have started to compensate with binary pulsars in probing the strong-field gravity. The first GW event of coalescing BNSs, GW170817, was detected by the LIGO/Virgo Collaboration in August 2017 Abbott et al. (2017b). GW170817 provides a powerful laboratory in the highly dynamical strong field. The spacetime of BNSs is strongly curved and highly dynamical in the vicinity of NSs in the late inspiral. If the DEF theory correctly describes the gravity, GW phase evolution of BNSs is modified. For now, limited by the sensitivity of the LIGO/Virgo detectors below tens of Hz, the precision to constrain the dipolar radiation from the short duration of GW170817 is still less than binary pulsars Shao et al. (2017).

Observations of BNSs at lower frequency are beneficial in the dipolar-radiation test. To increase the detector sensitivity, LIGO/Virgo have once more upgraded their equipments, and recently started the observing run 3 (O3) on April 1, 2019. Meanwhile, as the first kilometer-scale underground GW detector, the Kamioka Gravitational Wave Detector (KAGRA) Kuroda (2010); Abbott et al. (2018a) is likely to join O3 before the end of 2019 as well Akutsu et al. (2019). The next-generation ground-based GW detectors, such as the Cosmic Explorer (CE) led by the United States, and the Einstein Telescope (ET) led by the Europe Hild et al. (2011), will further improve the sensitivity in the future. In particular, they extend the sensitivity bands to be below 10 Hz. At the time of the third-generation detectors, we expect to discover more BNS merger events with higher sensitivities and larger signal-to-noise ratios (SNRs). These events are to put more stringent limits on alternative gravity theories, with the DEF theory being an important example.

In deriving constraints on the scalar-tensor gravity, the structure of NSs needs to be solved Damour and Esposito-Farèse (1993, 1996). Thus, the equation of state (EOS) of NS matters plays a role. The EOS is used to infer a NS radius, mass, and the scalar charge, by integrating the modified Tolman-Oppenheimer-Volkoff (TOV) equations Damour and Esposito-Farèse (1993, 1996). There are still large uncertainties in the NS EOS. In this work we choose nine EOSs that are all consistent with the maximum mass of NSs larger that $2\,\mathrm{M_{\odot}}$. Since more observations are being made for pulsars at radio and X-ray wavelengths, and BNSs with an increasing statistics, the uncertainty in the nuclear EOS is to be reduced in the near future.

Given an EOS, with inputs from binary pulsars and the recent BNS observation, we carry out advanced algorithms to constrain the DEF theory in a statistically sound way. In particular, we employ Bayesian methods through Markov-chain Monte Carlo (MCMC) simulations. Those simulations update posterior probability distribution of parameters in the DEF theory by evaluating the likelihood function millions of times. Every step requires the corresponding NS properties (including mass, radius, and scalar charge), which are derived from the modified TOV equations iteratively. Being computationally intensive, such studies have been carried out already in Refs. Shao et al. (2017); Anderson et al. (2019).

Being statistically sound, MCMC simulations lead to a large number of iterative calculations however. Here we build a new model to reduce the computational burden. Instead of solving the modified TOV equations iteratively and repeatedly, we construct reduced-order surrogate models (ROMs) to predict NS properties. There are two parameters characterizing the DEF theory, $\alpha_{0}$ and $\beta_{0}$ (see the next section). We explore a sufficiently large parameter space for them in the regime of strong-field scalarization. With our surrogate models, the process of obtaining NS properties is no longer an iterative integration, but a linear algebraic operation. It costs a fixed amount of time, much shorter than that in the previous method. In practice, for a given DEF theory, we use the central matter density $\rho_{c}$ of a NS to predict its radius $R$, mass $m_{A}$, and the effective scalar coupling $\alpha_{A}$. Our models are numerically accurate at $\sim 1\%$ level for $\alpha_{A}$, and better than $0.01\%$ for $m_{A}$ and $R$. They accelerate the processes of parameter estimation significantly everywhere in the parameter space we explore. With the speedup of those models, one can perform MCMC simulations much more efficiently yet still accurately. According to our performance test, they speed up calculations at least one hundred times, compared with the previous method in Ref. Shao et al. (2017). Various practical examples with binary pulsars and BNS events are demonstrated in this paper.

The organization of the paper is as follows. In Sec. II, we briefly review the nonperturbative spontaneous-scalarization phenomena for isolated NSs. The additional dipolar radiation and the modification of mass-radius relations for different EOSs in the scalar-tensor gravity will be discussed. Section III analyzes the difficulties in solving the modified TOV equations with large-scale calculations. We develop a better numerical method, and code it streamlinedly in the pySTGROM package. We make it public for an easy use for the community.222https://github.com/BenjaminDbb/pySTGROM In Sec. IV, with the speedup from pySTGROM, we stringently constrain the DEF theory by combining the dipolar-radiation limits from observations of five NS-white dwarf (WD) systems, and that of GW170817, which also includes a modified mass-radius relation. Our results are in good agreement with that from Shao et al. (2017). We also forecast constraints involving a hypothetical BNS event, to be detected by the Advanced LIGO at its design sensitivity Aasi et al. (2015), and next-generation ground-based GW detectors. Finally, the main conclusions and discussions are given in Sec. V.

II Spontaneous scalarization in the DEF gravity

In our study, we concentrate on the DEF theory. It is defined by the following action in the Einstein frame Damour and Esposito-Farèse (1993, 1996),

[TABLE]

In Sec. II, $g_{*}\equiv{\rm det}\,g_{\mu\nu}^{*}$ denotes the determinant of the Einstein metric $g_{\mu\nu}^{*}$, $R_{*}$ is the Ricci curvature scalar of $g_{\mu\nu}^{*}$, $G_{*}$ is the bare gravitational coupling constant, $\varphi$ is the dynamical scalar field that is added to GR, $\psi_{m}$ describes any matter fields, and $A(\varphi)$ is the conformal coupling factor that determines how $\varphi$ couples to $\psi_{m}$ in the Einstein frame.

The potential of the scalar field $V(\varphi)$ can be neglected for a slowly varying $\varphi$ compared with the typical scale of the system. Therefore, in our study we set $V(\varphi)=0$. Alternatively, Refs. Alsing et al. (2012); Ramazanoğlu and Pretorius (2016) considered the effects of a massive scalar field, via $V(\varphi)\approx 2m_{\varphi}\varphi^{2}$. The field equations of the DEF theory can be derived by varying the action (II). They are Damour and Esposito-Farèse (1993, 1996),

[TABLE]

where $T_{*}^{\mu\nu}\equiv 2c\left(-g_{*}\right)^{-1/2}\delta S_{m}/\delta g^{*}_{\mu\nu}$ denotes the matter stress-energy tensor and,

[TABLE]

As Eq. 3 shows, $\alpha(\varphi)$ measures the field-dependent coupling strength between the scalar field $\varphi$ and the trace of the energy-momentum tensor of matter fields, $T_{*}\equiv g^{*}_{\mu\nu}T_{*}^{\mu\nu}$.

In the JFBD theory, $\ln A(\varphi)$ was chosen to be linear in $\varphi$, i.e., $\ln A(\varphi)=\alpha_{0}\varphi$. It is extended to the polynomial form up to the quadratic order in the DEF theory Damour and Esposito-Farèse (1993),

[TABLE]

with $\alpha(\varphi)\equiv\partial\ln A(\varphi)/\partial\varphi=\beta_{0}\varphi$. We denote $\alpha_{0}=\beta_{0}\varphi_{0}$ with $\varphi_{0}$ being the asymptotic (cosmological background) scalar field value of $\varphi$ at infinity. Therefore, there are only two extra parameters, $\alpha_{0}$ and $\beta_{0}$ (or equivalently, $\varphi_{0}$ and $\beta_{0}$), to describe a DEF theory uniquely. In GR, we have $\alpha_{0}=\beta_{0}=0$.

Within the PPN framework in the weak field, Solar-System experiments can be used to narrow down the parameter space in the DEF theory. However, generally only the $\alpha_{0}$ or the combination $\beta_{0}\alpha_{0}^{2}$ can be constrained (see Refs. Damour (2007); Will (2014)).

For NSs, the aforementioned nonperturbative scalarization phenomena happen when Damour and Esposito-Farèse (1993); Barausse et al. (2013),

[TABLE]

From Fig. 1 in Ref. Shao et al. (2017), it is clear that with certain condition a negative $\beta_{0}$ can trigger an instability in the scalar field Sennett et al. (2017). It describes the strength of nonperturbative phenomena. The more negative for $\beta_{0}$ from the critical value $-4.0$, the more manifest the spontaneous scalarization in the strong-field regime.

In the strong field, the effective scalar coupling for a NS $A$,

[TABLE]

measures the “sensitivity” of the coupling between the NS mass and variations in the background scalar field $\varphi_{0}$. It appears directly in the Keplerian binding energy between two stars, $A$ and $B$,

[TABLE]

Besides the gravitational attraction in GR, the effective scalar coupling, $\alpha_{A}$, brings an additional scalar interaction. It also affects the strength of the orbital period decay of binaries Damour and Esposito-Farèse (1996).

In the following, we investigate the dipolar contribution to the orbital decay from the scalar field, $\dot{P}_{b}^{\rm dipole}$, and the quadrupolar contribution from the tensor field, $\dot{P}_{b}^{\rm quad}$. They read Peters and Mathews (1963); Damour and Esposito-Farèse (1993)

[TABLE]

[TABLE]

where $P_{b}$ is the orbital period, “$p$” and “$c$” denote the pulsar and its companion respectively. Subdominant contributions to $\dot{P}_{b}^{\rm quad}$, that are totally negligible to our study, can be found in Eq. (6.52d) in Ref. Damour and Esposito-Farèse (1992). In above equations, $g(e)$ and $f(e)$ are functions of the orbital eccentricity $e$,

[TABLE]

In Eqs. 9 and 10, $G_{N}=G_{*}\left(1+\alpha_{0}^{2}\right)$ denotes the Newtonian gravitational coupling in the weak field Damour (2007).

The effective scalar coupling equals to zero for BHs due to the no-hair theorem Abbott et al. (2016); Berti et al. (2015), and approaches to $\alpha_{0}$ for weak-field WDs in the DEF theory. Therefore, the influences from dipolar radiation on the orbital period decay occur dominantly in NS-WD and NS-BH binaries, as well as in asymmetric NS-NS binaries.333In the BH-BH systems, scalar charges of binaries both equal to zero. They hardly contribute to the orbital decay. But in the most general scalar-tensor theory, the Horndeski gravity Horndeski (1974), such systems are able to have scalar hairs which, depending on the details of the theory, might be the case only for certain mass ranges. Since the dipolar contribution relates to $-1\,$PN contribution in the GW phase, instructively, a precise bound on the gravitational dipole emission can be derived with the early inspiral stage of relevant GW events.

For a comprehensive study on NS spontaneous scalarization in the strong field, we here turn our attention to the nuclear EOS and the other NS properties. Usually, the NS matter can be treated as a perfect fluid. In GR, given an EOS, we can solve the classical TOV equations Tolman (1939); Oppenheimer and Volkoff (1939) for NSs. The NS radius, $R$, and mass, $m_{A}$, can be derived when the central matter density, $\rho_{c}$, and the EOS are given. In the DEF theory, we involve the additional scalar field $\varphi$. Correspondingly, the modified TOV equations [Eq. (7) in Ref. Damour and Esposito-Farèse (1993) and Eqs. (3.6a) to (3.6f) in Ref. Damour and Esposito-Farèse (1996)] should be used. In the Jordan frame,444The Jordan frame, also known as the physical frame, equivalents to the Einstein frame by a conformal transformation with redefinitions of the metric and the scalar field Dicke (1962); Weyl (1993). physical quantities, $R$, $m_{A}$, and $\alpha_{A}$, can be obtained via integrating from $\rho_{c}$ and the value of the scalar field at the center of a NS, $\varphi_{c}$.

We show an example of mass-radius relation for NSs in the DEF theory with $\log_{10}|\alpha_{0}|=-5.0$ and $\beta_{0}=-4.5$ in Fig. 1 Shao (2019). Two massive pulsars, PSRs J0740$+$6620 Cromartie et al. (2019) and J0348$+$0432 Antoniadis et al. (2013), show that the maximum mass of NSs is most likely above $2\,\mathrm{M_{\odot}}$. Therefore, we choose the EOSs that satisfy this condition in our study. In total, we use nine EOSs. From Fig. 1 we can see that,

the mass-radius relation in the DEF theory is very close to that of GR in the majority of mass range; 2. 2.

there are “bumps” that correspond to the spontaneous scalarization phenomena and, when compared with GR, the bump for each EOS predicts a larger radius for a same NS mass; 3. 3.

in GR, the EOSs, SLy4 and H4, are getting more and more in tension with pulsar mass measurements, because their maximum masses hardly reach the 1-$\sigma$ lower limit of PSR J0740$+$6620 Cromartie et al. (2019).

It is worthy to note that, the EOSs H4 and PAL1 predict relatively larger radii for NSs. They are starting to be disfavored by the LIGO/Virgo’s observation of GW170817 Abbott et al. (2018b, 2019c).

III The reduced-order model

As discussed in the previous section, the initial values for the TOV integrator include the NS center matter density $\rho_{c}$, and the center value of the scalar field $\varphi_{c}$. The macroscopic quantities of NSs, namely $R$, $m_{A}$, and $\alpha_{A}$, can be obtained from integrating the modified TOV equations. From the viewpoint of the DEF gravity, instead of $\varphi_{c}$, it is more convenient to have the scalar field at infinity, $\varphi_{0}$, as an initial value. There is correspondence between $\varphi_{0}$ and $\varphi_{c}$. Thus, to obtain a desired $\varphi_{0}$, we need to find the value of $\varphi_{c}$. To solve such a boundary value problem, the shooting method turns out to be a practical way Shao et al. (2017). Nevertheless, sometimes it leads to a large number of iterations, especially around the region of spontaneous scalarization. Therefore, it is not efficient to perform large-scale calculations, such as the parameter estimation with the MCMC approach, based on the shooting method. It is helpful to find a faster and more efficient, but still accurate, surrogate algorithm for a better performance in solving the NS structure in the DEF theory.

Before discussing the surrogate algorithm for the DEF theory, we turn to GW data analysis for idea. A similar issue, that is equally computationally challenging and difficult, arises when considering GW data analysis Veitch et al. (2015). In the case of compact binary coalescences, the matched-filtering technique is used. It involves matching the data against a set of template waveforms (see e.g. Ref. Bohé et al. (2017)), to infer a potential astrophysical signal. The analysis could lead to a large computational cost. Therefore, effective approaches, the surrogate models, were built.

Distinct examples are developed by Pürrer (2014, 2016) and Field et al. (2014, 2011). This kind of model is constructed by a set of highly accurate sparse waveforms. It brings a new waveform model to provide fast and accurately compressed approximations. These accurate surrogate models are built without sacrificing the underlying accuracy for data analysis. Recently, several such methods have been proposed in the literature. Based on the singular value decomposition (SVD) method, time-domain inspiral surrogate waveforms are generated in Ref. Pürrer (2014). Another popular construction method is the greedy reduced basis method, which is usually combined with the empirical interpolation method (EIM) Chaturantabut and Sorensen (2010); Maday et al. (2009). It has been applied to GW waveforms Field et al. (2014) and BH ringdown Caudill et al. (2012).

Inspired by the above optimization for large-scale calculations in the GW data analysis, we adopt the greedy reduced basis method to the NS structure in the DEF theory. In our surrogate models, for each EOS, the NS properties are accurately and rapidly inferred with given NS central matter density $\rho_{c}$. In the following subsections, we briefly introduce the processes of general ROM, and describe in detail how to specialize this model to the DEF theory. Finally, we assess our ROMs’ accuracy, and find that the final results are consistent with our expectation. Our ROMs will be a helpful tool in generating NS properties with essential speed and accuracy for future studies.

III.1 A brief technical introduction to ROM

The theoretical aspects and the processes of building the ROM have been discussed extensively in Fig. 1 and Appendices A to D in Ref. Field et al. (2014). We here give a brief overview.

Let us use the GW data analysis as a prototype, and denote a curve $h(t;\bm{\lambda})$, representing a waveform with parameters $\bm{\lambda}$ (for GW waveforms, parameters in $\bm{\lambda}$ include the masses, the spins, and so on). In order to generate a ROM for $h(t;\bm{\lambda})$, we need the data of $h(t)$, produced with a set of given parameters $\left\{\bm{\lambda}_{i}\right\}$ on a grid. We denote the training space ${\mathbf{V}\equiv\left\{h(t;\bm{\lambda}_{i})\right\}}$. With the reduced basis (RB) method, we select a certain number of bases derived from the training space $\mathbf{V}$. Those bases are regarded as a set of representative bases for the remaining waveforms. Instead of choosing $m$ existing $\{\bm{\lambda}_{i}\}$ and their corresponding space $\{h(t;\bm{\lambda}_{i})\}$, we seek a more complete set of bases that are as independent as possible from each other.

The greedy algorithm is currently one of the methods to select $m$ orthonormal RBs (see Ref. Binev et al. (2011) and Appendix. A in Ref. Field et al. (2014) for details). The corresponding chosen space is $\mathbf{RV}=\left\{e_{i}\right\}_{i=1}^{m}$. Actually, the process of generating orthonormal RBs is mentioned as iterated Gram-Schmidt orthogonalization algorithm with greedy selection Ruhe (1983); Giraud et al. (2005); Hoffmann (1989); Antil et al. (2018). Firstly, a seed, which can be any waveform in $\mathbf{V}$, is chosen as the starting RB $(i=0)$. Then, we define the maximum projection error,

[TABLE]

and a user-specified tolerable error bound, $\Sigma$, to infer the next RB iteratively in the way that is explained here. In Eq. 13, $\mathcal{P}_{i}$ describes the projection of $h(t;\bm{\lambda})$ onto the span of the first $i$ RBs. The waveform corresponding to the maximum $\sigma_{i}$ is chosen as the next $(i+1)$-th RB, $e_{i+1}(t)$, after orthogonalized by the iterated Gram-Schmidt algorithm. In practice, $\Sigma$ is set as a threshold to terminate iterations of the greedy selection. When $\sigma_{m-1}\lesssim\Sigma$, the desired orthonormal $m$ bases are complete.

With RBs obtained above, every waveform in the training space is well approximated by an expansion of the form,

[TABLE]

In Eq. 14, $c_{i}(\bm{\lambda})$ is the corresponding coefficient of the $i$-th RB. It is calculated by a special “inner product” in the space $\mathbf{RV}$. We collect such coefficients derived from $\mathbf{V}$ as the greedy data for the ROM. After that, we perform the EIM algorithm (see Appendix B in Ref. Field et al. (2014) for details) to identify $m\,$ time-samples, which we call the empirical nodes. An interpolation can be built to accurately reconstruct any fiducial waveform by the RBs. Finally, at each empirical node, we perform a fit to the parameter space, $\{\bm{\lambda}_{i}\}$, and finish the construction of the ROM.555Particularly, a 2-dimensional 5th spline interpolating fit is used in our study. A different 2-dimensional fitting method works without practical difference. For convenience, the fitted data are saved in a model file.

Now, we assume that a waveform, parametrized by an unknown set of parameters $\bm{\lambda}_{\star}$, is required in the calculation. The parameters $\bm{\lambda}_{\star}$ are different from $\bm{\lambda}_{i}$, but within the boundary ranges of $\left\{\bm{\lambda}_{i}\right\}$. We can use the ROM and the fitted data on each empirical node to get a complete waveform as a function of $t$, $h(t;\bm{\lambda}_{\star})$. The value of $h(t_{\star};\bm{\lambda}_{\star})$ is derived for a given $t_{\star}$ in a fast and accurate manner.

III.2 Constructing ROMs for the DEF gravity

In this subsection, we specialize the greedy reduced basis method to build the ROMs for the DEF theory. The detailed approach is explained according to the steps outlined above. Because the idea of building ROMs is “borrowed” from the GW waveform studies, we use the word “waveform” to represent the desired functional forms. The parameters $\bm{\lambda}$ are specialized as $(\alpha_{0},\beta_{0})$.

Ideally we want to get the NS properties, $R$ and $\alpha_{A}$, as a function of the mass parameter, $m_{A}$. But in some parameter space pathological behaviors might happen, while simply integrating the modified TOV equations. We use the DEF theory with $\log_{10}|\alpha_{0}|=-5.3$ and $\beta_{0}=-4.8$ for the EOS AP4 as an example. We illustrate the different relations of NS effective scalar coupling $\alpha_{A}$, radius $R$, and center matter density $\rho_{c}$ in Fig. 2. In the left panel, the effective scalar coupling can obtain a value significantly larger than $0.1$ within the NS mass interval $m_{A}\in\left[1.4\,\mathrm{M_{\odot}},2.0\,\mathrm{M_{\odot}}\right]$. This strength of spontaneous scalarization in general decreases monotonically and rapidly when $m_{A}\gtrsim 2.0\,\mathrm{M_{\odot}}$. But, the multivalued relations between $\alpha_{A}$ and $m_{A}$ arise when $\beta_{0}\lesssim-4.6$, namely, multiple values for $\alpha_{A}$ are possible for a given $m_{A}$. The curve in the middle panel shows the mass-radius relation where the “bump” happens (see Fig. 1 for a similar behavior). There exists the multivalued relation as well. In the bump region, an $m_{A}$ can correspond to multiple values of $R$ within the interval $m_{A}\in\left[2.0\,\mathrm{M_{\odot}},2.05\,\mathrm{M_{\odot}}\right]$, as indicated in grey. Similarly, in the right panel, the mass of NS has a sudden drop when $\rho_{c}$ exceeds a certain value in the grey band.

In Fig. 2, in the interval $m_{A}\in\left[2.0\,\mathrm{M_{\odot}},2.05\,\mathrm{M_{\odot}}\right]$, a fixed NS mass $m_{A}$ can correspond to multiple values for $\alpha_{A}$ and $R$. The existences of those phenomena have also been confirmed in the other EOSs, when $\beta_{0}$ is below a critical value. For the EOS AP4, the critical value approximates to $-4.6$. Combining this observation and the relation of $m_{A}$ and $\rho_{c}$ in the right panel of Fig. 2, we have become fully aware where the pathologies come from. It is due to the excessive density at the center (about $\gtrsim 1.5\times 10^{15}\,{\rm g\,cm^{-3}}$ for the EOS AP4), which causes the NS to further collapse into a BH. In order to avoid dealing with those multivalued relations in the ROMs, we promote the implicit parameter $\rho_{c}$, to an independent variable which corresponds to the “time” $t$ in Sec. III.1. Therefore, we finally decide to construct three independent ROMs for the DEF theory, $m_{A}\left(\rho_{c}\right)$, $R\left(\rho_{c}\right)$, and $\alpha_{A}\left(\rho_{c}\right)$.666Actually, we use $\ln\left|\alpha_{A}\right|$ instead of $\alpha_{A}$ in the ROM for a better performance. Those relations are all single-valued. The NS properties, $(R,m_{A},\alpha_{A})$, are expected to be derived with a given parameter $\rho_{c}$ injectively in the ROMs.

For different EOSs, the ranges of $\rho_{c}$ are different. We choose the range of $\rho_{c}$ to have $m_{A}\in(1\,{\rm M}_{\odot},m_{A}^{\rm max})$ where the maximum NS mass $m_{A}^{\rm max}$ is EOS-dependent. The spacing $\Delta\rho_{c}$ is selected according to the steepness of the scalarization in different regions where we choose more points in rapidly changing parameter space. After building the ROMs, we can have results for all $\rho_{c}$ in its range directly from our ROMs.

In a short summary, we build three ROMs, $m_{A}(\rho_{c})$, $R(\rho_{c})$, and $\alpha_{A}(\rho_{c})$. Though with one more implicit parameter $\rho_{c}$. This choice has avoided the pathologies if we had built two ROMs, $R(m_{A})$ and $\alpha_{A}(m_{A})$.

For the sake of completeness, we should generate all training waveforms of interest at locations lying in the parameter space grid,

[TABLE]

In Eq. 15, $\mathcal{A}$ and $\mathcal{B}$ are 1-dimensional sets covering the desired parameter ranges. They are representing the parameters in the DEF theory, $\alpha_{0}$ and $\beta_{0}$, respectively. The spacing in the parameters of $\log_{10}|{\alpha_{0}}|$ and $-\beta_{0}$ is chosen in an uneven way, shown in Fig. 3. We concentrate more samples near the rapidly changing domain in the greedy data for each empirical node. Particularly, we choose $\mathcal{A}$ to cover the range of $\log_{10}|\alpha_{0}|\in[-5.3,\,-2.5]$ with 42 nodes. For $\beta_{0}$, 101 nodes are chosen to cover the interval $-\beta_{0}\in[4.0,4.8]$. Totally, we have involved $42\times 101=4242$ sparse waveforms to form the training space $\mathbf{V}$ for constructing the ROMs. For the boundary values in building the ROMs, the particular value $\alpha_{0}\approx 10^{-2.5}$ approximately equals to the upper limit given by the Cassini spacecraft Bertotti et al. (2003); Damour (2007), and the critical value $\beta_{0}\lesssim-4.0$ indicates places where the spontaneous scalarization happens in the DEF theory Damour and Esposito-Farèse (1993); Barausse et al. (2013). These choice will also serve as our priors in the parameter estimation that we will discuss in Sec. IV.

The greedy selection of RBs from the training space $\mathbf{V}$ is iteratively performed to generate the next RB until the desired accuracy is achieved, namely $\sigma_{i}\lesssim\Sigma$. In practice, $\Sigma$ is set to terminate the iteration of the greedy selection. In our ROMs, we choose $\Sigma=10^{-7}$ for $R$ and $m_{A}$, and $\Sigma=10^{-5}$ for $\alpha_{A}$.

In the iterative process, the relative projection errors, $\tilde{\sigma}_{i}\equiv\sigma_{i}/\sigma_{0}$, are recorded for convenience. We illustrate $\tilde{\sigma}_{i}$ as a function of the basis size in Fig. 4. In the figure, the most notable features are the decreasing speed for $\tilde{\sigma}_{i}$ with an increasing basis size. As the basis size increases, the $\tilde{\sigma}_{i}$’s of $R$ and $m_{A}$ quickly arrive at the level $\lesssim 10^{-8}$ within the basis size of 20. Particularly, the ROM of $m_{A}(\rho_{c})$ has the fastest decline. In contrast, $\tilde{\sigma}_{i}$ falls off smoothly to $\Sigma\sim 10^{-5}$ for about 150 steps in building the ROM of $\alpha_{A}$. This means that, more RBs are needed to ensure the accuracy of $\alpha_{A}$. Actually, considering the tolerable error involved by the shooting method when integrating the modified TOV equations, which is about $\sim\,1\%$, we find that the precision loss of the above processes in constructing ROMs almost negligible. This conclusion is verified when assessing the accuracy of the ROMs in Sec. III.3.

Finally, we follow the steps in Sec. III.1 to obtain RBs with the EIM algorithm, and build the entire ROMs. Those ROMs can generate NS properties, $R$, $m_{A}$, and $\alpha_{A}$, efficiently and rapidly as a function of $\rho_{c}$. We tested the ROMs with randomly generated parameters, and found that the time for one solve of the modified TOV equations improves from $\sim 300\,$milliseconds to $\sim 1\,$millisecond on our Intel Xeon E5 computers. We implement those three ROMs for the DEF theory in the pySTGROM package. As is shown with applications in Sec. IV, our ROMs can be performed at least two orders of magnitude faster than the shooting method.

III.3 Assessing the ROMs of the DEF gravity

After the ROMs of the DEF gravity are constructed, we need to assess their accuracy. We define,

[TABLE]

that measures the fractional accuracy of the ROMs. In Eq. 16, $\mathcal{E}(X)$ denotes the relative error with respect to the true value. The values we predict by the ROMs are denoted as $X_{\rm ROM}$; the values generated by the shooting algorithm (with a relative tolerable error $\sim 1\%$) are denoted as $X_{\rm mTOV}$.

We randomly generate $5\times 10^{5}$ sets of parameters for $\alpha_{0}$, $\beta_{0}$, and $\rho_{c}$ in the valid ranges of the ROMs. We then generate the corresponding values for $X_{\rm ROM}$ and $X_{\rm mTOV}$ from the ROMs and the shooting algorithm respectively. $\mathcal{E}(X)$ is calculated according to Eq. 16. The distributions of $\mathcal{E}(X)$ are illustrated in Fig. 5. As the figure shows, the relative errors of $R$ and $m_{A}$ are $\lesssim 10^{-5}$. In contrast, $\mathcal{E}(\alpha_{A})$ is approximately three orders of magnitude larger. The upper limit of $\mathcal{E}(\alpha_{A})$ is less than the tolerable error, $1\%$, of the shooting method. Therefore, these results are consistent with the relative maximum projection error discussed in Fig. 4. Notice that, the parameters used in making Fig. 5 are randomly generated, not necessarily on the grid in Fig. 3. Therefore, the results in Fig. 5 have included all errors for our ROMs, including the interpolating errors. Although the ROM of $\alpha_{A}$ leads to a larger error relative to those of $R$ and $m_{A}$, the error can be neglected when compared with the error from the shooting method. It has no noticeable influence in the parameter estimation that we are to discuss in Sec. IV.

IV Constraints from binary pulsars and gravitational waves

In this section, we apply our ROMs to various scenarios, and discuss the improvement in deriving NS properties. In these applications, we use the MCMC technique Brooks et al. (2011); Foreman-Mackey et al. (2013); Shao et al. (2017) to constrain the DEF theory.

IV.1 The setup

As mentioned before, two independent parameters, $\alpha_{0}$ and $\beta_{0}$, are needed to characterize the DEF theory. One of them, $\alpha_{0}$, has been well constrained by the Cassini mission Bertotti et al. (2003), that measures the Shapiro time delay in the weak-field regime. It gives a limit, $|\alpha_{0}|<3.4\times 10^{-3}$ at 68% confidence level (CL) Damour (2007). In the nonperturbative strong field, the DEF theory can be significantly different from GR. The constraints in the strong field were studied in Refs. Freire et al. (2012); Antoniadis et al. (2013); Wex (2014); Shao et al. (2017); Anderson et al. (2019).

NSs, whose gravitational binding energy can be as large as $\sim 20\%$ of their rest-mass energy, are powerful objects to test the strong-field gravity. We concentrate on the observations relevant to NSs. We consider (i) systems of binary pulsars in the quasi-stationary regime, and (ii) GWs from BNSs in the highly dynamical regime. With the construction of ROMs in Sec. III, we now have a numerically faster way to derive NS properties. The current and future constraints for the DEF theory are investigated in the following.

Some stringent limits via the MCMC approach to the spontaneous scalarization of the DEF theory have been obtained in Ref. Shao et al. (2017). Five asymmetric NS-WD binaries were used. They are PSRs J0348+0432 Antoniadis et al. (2013), J1012+5307 Lazaridis et al. (2009); Antoniadis et al. (2016); Desvignes et al. (2016), J1738+0333 Freire et al. (2012), J1909$-$3744 Reardon et al. (2016); Desvignes et al. (2016); Arzoumanian et al. (2018) and J2222$-$0137 Cognard et al. (2017). Relevant parameters for these binary pulsars are listed in Table 1. The WD companion corresponds to a weakly self-gravitating object. It has a tiny effective scalar coupling $\alpha_{c}\approx\alpha_{0}$. The small effective scalar coupling of the WD would lead to a large dipole term, $\propto(\alpha_{p}-\alpha_{c})^{2}$, in the parameter space we are investigating [see Eq. 9].

Pulsar parameters and orbit parameters are measured by the TOAs of pulses, including Keplerian and post-Keplerian parameters. Some parameters, such as the time derivative of the orbital period, $\dot{P}_{b}$, are functions of the masses of the pulsar and its companion Damour and Taylor (1992). The chosen systems are all well timed. The uncertainty in $\dot{P}_{b}$ can be very small for the pulsars we consider (see Table 1). Such precise measurements place strong bounds on various alternative theories of gravity Wex (2014).

On the other hand, a completely new era for testing highly dynamical strong field with NSs has began with GW170817 Abbott et al. (2017b). It offers an extraordinary opportunity, completely different from previously detected BBH merger events. From the phase of GW170817, one can derive the NS properties, for instance, the mass and radius of each NS. The LIGO/Virgo Collaboration used EOS-insensitive relations to derive those properties in Ref.  Abbott et al. (2019c). The results are shown in Table 2, and they are included to constrain the DEF theory in the MCMC approach.777Different from the LIGO/Virgo Collaboration, De et al. (2018) used another EOS-insensitive relation. They obtained similar results.

In binary pulsars, timing parameters are measured by radio techniques. Some of them have achieved extremely high precision from decades of observation. In contrast, the NS properties derived from GW170817, have a relatively poor precision. But, the radii of NSs were inferred from GW170817 Abbott et al. (2017b, 2019c). This measurement is unique to constrain the alternative theories of gravity.

In addition to the systems mentioned above, we expect that, more BNSs are to be detected with the future GW detectors. KAGRA, operated in Japan, has two 3 km orthonormal arms to form an underground GW interferometer. It has a similar designed sensitivity as that of the Advanced LIGO Akutsu et al. (2019). Moreover, the next-generation ground-based detectors are expected to detect more GW events with larger SNRs, due to their even higher sensitivities. CE Abbott et al. (2017d) is an L-shaped 40 km interferometer (compared to 4 km for LIGO, and 3 km for Virgo and KAGRA). It is to be built on the experience and success of current detectors. It is roughly ten times more sensitive than the Advanced LIGO. Another proposed next-generation GW detector, ET Hild et al. (2011), is a 10 km interferometer. It is supported by the European Commission. There are three arms forming an equilateral triangle, located underground to reduce seismic noises. It is also ten times more sensitive than the Advanced LIGO. Particularly, in the frequency band below $\sim 10$ Hz, ET can achieve a higher sensitivity than CE Abbott et al. (2017d). We want to find out how the DEF theory will be constrained with those future GW detectors.

Shibata et al. (2014) discovered the fact that, depending on the EOS one can still have strong scalarization in a mass range that is not yet constrained by pulsar experiments. After accounting for new pulsar tests, Shao et al. (2017) showed that there is a so-called “scalarization window” at $m_{p}\sim 1.7\,\mathrm{M_{\odot}}$. This is the place where there could still have a large deviation from GR, given all the current constraints. If the future GW detectors can observe asymmetric BNSs with masses around $1.6$–$1.7\,\mathrm{M_{\odot}}$, this window could be closed. In this paper, we assume that, a GW event from an inspiraling BNS with masses, $\left(1.65\,\mathrm{M_{\odot}},1.22\,\mathrm{M_{\odot}}\right)$, is detected at luminosity distance $D_{L}=200\,{\rm Mpc}$. This BNS system has similar NS masses as those of PSR J1913+1102 Lazarus et al. (2016). We denote this hypothetical BNS as BNS@200Mpc. Then, we perform MCMC simulations to investigate this hypothetical BNS. In the foreseeable future, BNS coalescence will be observed in large numbers. Correspondingly, a tight constraint on parameters of the DEF theory, $\alpha_{0}$ and $\beta_{0}$, can be obtained.

Some properties of NSs are governed by the EOS. NS EOS describes the relation between pressure and density of NS matters. It is involved in the modified TOV equation integration. However, because of our lack of knowledge about the inner structure of NSs, the EOS is still uncertain. Following our previous discussion, there exists a lower limit for the maximum mass of NSs, $M_{\rm max}\gtrsim 2\,\mathrm{M_{\odot}}$. Nine EOSs, AP3, AP4, ENG, H4, MPA1, PAL1, SLy4, WFF1, and WFF2, are adopted in this study (see Ref. Lattimer (2012) for a review). We have illustrated the mass-radius relations of NSs for these EOSs in Fig. 1. They are all consistent with the above maximum-mass limit. Moreover, for our studies, we believe that they are sufficient to investigate the EOS-dependent properties with spontaneous scalarization Shibata et al. (2014).

IV.2 The MCMC framework

Combining the observations of five pulsars, GW170817, and BNS@200Mpc, we perform MCMC simulations with each of these EOSs. MCMC techniques update posterior distributions of underlying parameters. After convergence, those distributions will be consistent with astrophysical observations by evaluating the likelihood function. We use the python implementation of an affine-invariant MCMC ensemble sampler, EMCEE888https://github.com/dfm/emcee. We use the ROMs, that is built in Sec. III and coded in pySTGROM package, to speed up the MCMC calculations for parameter estimation within the Bayesian framework Del Pozzo and Vecchio (2016); Shao et al. (2017).

In the Bayesian inference, given priors, the posterior distribution of $(\alpha_{0},\beta_{0})$ can be inferred with data, ${\cal D}$, and a hypothesis, ${\cal H}$. We use the formula of Bayes’ theorem,

[TABLE]

where $P\left(\alpha_{0},\beta_{0}|{\cal D},{\cal H},{\cal I}\right)$ is an updated (marginalized) posterior distribution of $(\alpha_{0},\beta_{0})$, $P\left({\cal D}|\alpha_{0},\beta_{0},\bm{\Xi},{\cal H},{\cal I}\right)\equiv\mathcal{L}$ is the likelihood function, $P\left(\alpha_{0},\beta_{0},\bm{\Xi}|{\cal H},{\cal I}\right)$ is the prior on parameters $(\alpha_{0},\,\beta_{0},\,\bm{\Xi})$, and $P\left({\cal D}|{\cal H},{\cal I}\right)$ is the model evidence. In Eq. 17, the hypothesis, ${\cal H}$, represents the DEF theory, ${\cal I}$ denotes all other relevant knowledge, $\bm{\Xi}$ collectively denotes all other unknown parameters in addition to $(\alpha_{0},\,\beta_{0})$.

For our studies of the DEF theory, the priors of $(\alpha_{0},\beta_{0})$ are chosen carefully to cover the interesting region where the spontaneous scalarization occurs. In order to speed up the MCMC simulations, the values of $(\log_{10}\left|\alpha_{0}\right|,\beta_{0})$ are restricted to the same rectangle region as in our ROMs that are built in Sec. III.2. We pick a uniform prior on $\log_{10}\left|\alpha_{0}\right|$ in the range $[-5.3,\,-2.5]$. The prior on $\beta_{0}$ is chosen to be uniform in the interval $\beta_{0}\in[-4.8,\,-4.0]$.

The parameters, $\alpha_{0}$ and $\beta_{0}$, can be constrained by evaluating the likelihood function. Due to their observational characteristics, the five pulsars in Table 1, GW170817, and BNS@200Mpc contribute differently to the likelihood function. Their corresponding log-likelihood functions are given separately as follows.

For binary pulsars, we have the log-likelihood function,

[TABLE]

for $N_{\rm PSR}$ binary pulsar systems (we have $N_{\rm PSR}=5$ in this study). We have assumed that observations with different binary pulsars are independent. For each pulsar, the intrinsic orbital decay, $\dot{P}_{b}^{\rm int}$, the pulsar mass, $m_{p}^{\rm obs}$, and the companion mass, $m_{c}^{\rm obs}$ are given in Table 1 when applicable. In some cases, we have the mass ratio $q\equiv m_{p}/m_{c}$ and $m_{c}^{\rm obs}$ instead; it is easy to obtain $m_{p}^{\rm obs}$. The 1-$\sigma$ uncertainty, $\sigma_{X}$ in Eq. 18, is the observational uncertainty for the quantity $X$, where $X\in\left\{\dot{P}_{b}^{\rm int},m_{p}^{\rm obs},m_{c}^{\rm obs}\right\}$. The predicted orbital decay, $\dot{P}_{b}^{\rm th}$, from the DEF theory equals to $\dot{P}_{b}^{\rm dipole}+\dot{P}_{b}^{\rm quad}$ [see Eqs. 9 and 10]. The quantities, $\dot{P}_{b}^{\rm th}$ and $m_{p}$, are implicitly dependent on the parameter set, $\left(\alpha_{0},\beta_{0},\bm{\Xi}\right)$. During the MCMC simulations, they are derived from those parameters. The companion mass, $m_{c}$, is picked randomly within its 1-$\sigma$ uncertainty.

For the full calculation, it is worth noting that, the orbital period, $P_{b}$, and the orbital eccentricity, $e$, should also be included in Eq. 18. Those quantities have been determined very well from the observations. We adopt their central values directly for simplifying the MCMC processes. It is verified that, the uncertainties of them have little effect on our final limits of $(\alpha_{0},\beta_{0})$.

In general, the log-likelihood function for $N_{\rm GW}$ GWs can be expressed as,

[TABLE]

Here, the properties of BNSs, such as the masses, $m_{j}$ $(j=1,2)$, the radii, $R_{j}$, are given in Table 2. The 1-$\sigma$ uncertainties of those properties, $\sigma_{m_{j}^{\rm obs}}$ and $\sigma_{R_{j}^{\rm obs}}$, are given at the 68% CL. The function for the dipole radiation, ${\tt D}_{i}\left(|\alpha_{1}-\alpha_{2}|,\,|\Delta\alpha|^{\rm upper}\right)$ , describes the contribution of the $i$-th GW’s dipolar radiation to $\ln\cal L_{\rm GW}$. It returns 0, when the effective scalar couplings of the BNS, $\alpha_{1}$ and $\alpha_{2}$, satisfy $\left|\alpha_{1}-\alpha_{2}\right|=\left|\Delta\alpha\right|\leq\left|\Delta\alpha\right|^{\rm upper}$, otherwise it returns $-\infty$. Here $\left|\Delta\alpha\right|^{\rm upper}$ is the upper limit of the absolute difference between $\alpha_{1}$ and $\alpha_{2}$ from observations. It is derived from the relation, $B=5\left|\Delta\alpha\right|^{2}/96$ Barausse et al. (2016), where $B$ is a theory-dependent parameter regulating the strength of the dipolar radiation in scalar-tensor theories. The presence of dipole radiation in GW170817 is constrained to be $B\leq 1.2\times 10^{-5}$ Abbott et al. (2019b). Therefore, the upper bound of $\left|\Delta\alpha\right|$, $\left|\Delta\alpha\right|^{\rm upper}\simeq 0.015$, is obtained for GW170817. Notice that we are using NS masses and radii derived from GR. Because here we are constraining the non-GR effects, we feel our approach sufficient. However, at the stage of discovering non-GR effects, we need a fully consistent approach that uses NS masses and radii derived from the scalar-tensor gravity, instead of GR.

In addition to GW170817, we introduce a hypothetical GW from an inspiraling BNS, BNS@200Mpc, for forecasting future constraints. It is meaningful to investigate the improvement with future GW detectors. The log-likelihood function of $N_{\rm GW}^{\star}$ hypothetical GWs is chosen as,

[TABLE]

In Eq. 20, $\sigma(|\Delta\alpha|)$ is the expected 1-$\sigma$ uncertainty of $\left|\Delta\alpha\right|$. It is obtained from the approach of the Fisher information matrix. We will introduce this approach briefly and perform the MCMC calculations with the log-likelihood function, $\ln\mathcal{L}_{\rm BNS}$, in Sec. IV.4. Compared with Eq. 19, we neglect the contributions from masses and radii, because, given the starting frequency of CE and ET, the dipolar radiation contribution is more constraining.

In a short summary, the likelihood function we use in Eq. 17 is the sum of Eqs. 18, 20 and 19 that have included all contributions from binary pulsars, GW170817, and future BNSs.

Now, we will explain how to employ the MCMC technique to get the posteriors from the priors on $(\alpha_{0},\beta_{0})$ and the log-likelihood functions. Combining the observations of binary pulsars and GWs, we have $N=N_{\rm PSR}+2\,\left(N_{\rm GW}+N_{\rm GW}^{\star}\right)$ NSs to constrain the $(\alpha_{0},\beta_{0})$ parameter space. To fully describe the contribution of the gravitational dipolar radiation of these systems, $N+2$ free parameters, denoted collectively as $\bm{\theta}$, are required. The parameters, $\bm{\theta}$, include $\alpha_{0}$, $\beta_{0}$, and $\rho_{c}^{(i)}(i=1,\ldots,N)$ Shao et al. (2017). The initial central matter densities, $\rho_{c}^{(i)}$, are needed in the Jordan frame. They are fed to our ROMs to derive the NS properties. Initially they will be sampled around their GR values, but they are allowed to explore a sufficiently large range in the MCMC process. In principle, if we were using the shooting method, the initial central values of the scalar field, $\varphi_{c}^{(i)}$, are needed as well. In our ROMs, they are no longer involved. Those NS properties can be obtained uniquely with $\alpha_{0}$, $\beta_{0}$, and $\rho_{c}^{(i)}$. Then, those derived properties are passed to the log-likelihood functions, Eqs. 18, 20 and 19, to evaluate the posteriors.

The constraints from the same five binary pulsars have been investigated in great detail in Ref. Shao et al. (2017). In this paper, in order to verify the validity of our newly built ROMs, we will reproduce their result. In addition to the five pulsars, we also use GW170817 and a hypothetical BNS, BNS@200Mpc. The BNS@200Mpc is assumed to be detected by the Advanced LIGO at its designed sensitivity, and next-generation detectors, CE and ET. Now, we can use those observations to constrain the DEF theory, and obtain the bounds in the $(\alpha_{0},\beta_{0})$ parameter space. Furthermore, we also quantify how much the tests will be improved with those next-generation detectors.

For each EOS, we perform six separate MCMC runs with different combinations of observations. These six scenarios are shown in Table 3. Here the investigations of the DEF theory are divided into two catalogs.

• •

Scenarios (i) to (iii) correspond to the present observations that contain five pulsars and GW170817. We will compare the constraining results of five pulsars with those of GW170817. The corresponding constraints from binary pulsars and GW170817 are given in Sec. IV.3.

• •

Scenarios (iv) to (vi) involve five pulsars and a hypothetical BNS@200Mpc to be observed by different GW detectors. We expect to obtain tighter constraints compared with present observations. Corresponding constraints from binary pulsars and a hypothetical BNS are investigated in Sec. IV.4.

In the following subsections, we use the EOS AP4 as an example. At the end of this section, a total of 54 MCMC runs (6 scenarios $\times$ 9 EOSs) are discussed.

For each run, we produce 2.6 millions of samples in total using multiple chains (26 chains $\times\,10^{5}$ samples for each chain). According to the guides in Refs. Brooks et al. (2011); Foreman-Mackey et al. (2013), the first half chain samples of these 54 runs are discarded as the burn-in phase. Then, we “thin” remaining samples, with a thinning factor of ten, to reduce the correlation of adjacent points. Finally, there are $1.3\times 10^{5}$ “thinned” samples remaining for each scenario. In our studies, it is verified that, those “thinned” samples have passed the Gelman–Rubin convergence diagnostic very well Gelman and Rubin (1992). It indicates that, all parameters in $\bm{\theta}$ have lost the memory of their initial values, such that they can be used to infer the parameters including $(\alpha_{0},\beta_{0})$.

IV.3 Constraints from binary pulsars and GW170817

In this subsection, we investigate three scenarios with real observations, GW170817, PSRs, and PSRs+GW170817. We use the EOS AP4 as an example. Following our previous discussion, we distribute initial values of $\log_{10}|\alpha_{0}|$ and $-\beta_{0}$ uniformly in the rectangle region of the parameter space. Then, after the MCMC simulations, we obtain the posterior distributions, from where we can get upper limits of theory parameters, $\alpha_{0}$ and $\beta_{0}$.

First, we investigate the scenario GW170817. The properties of GW170817 are given in Table 2. They are used to perform the MCMC simulations. The marginalized 2-dimensional distribution in the parameter space of $(\log_{10}|\alpha_{0}|,-\beta_{0})$ is illustrated in Fig. 6. It is evident that, GW170817 almost has no constraint on $\alpha_{0}$. Actually, it is within our expectation. Because of its short duration, compared with binary pulsar observations, the parameters are derived with lower precision. Current data are not accurate enough to give tight constraints. Nevertheless, the upper limit of $-\beta_{0}$ is constrained to be $\lesssim 4.5$ at 90% CL. It is consistent with the argument that, $\beta_{0}$ plays the major role in controlling the strength of scalarization. Therefore, $\beta_{0}$ is constrained.

We use five binary pulsars for the scenario PSRs. Different from the illustration in Fig. 6, the results of the scenario PSRs show very good constraints on the parameters of the DEF theory in Fig. 7. They are consistent with results in Ref. Shao et al. (2017) (in particular, we compared our results with Table II of Ref. Shao et al. (2017)). In Fig. 7, the posterior samples are gathered in the lower left corner in the marginalized 2-dimensional distribution. The corresponding parameters, $\log_{10}|\alpha_{0}|$ and $-\beta_{0}$, are constrained. Especially, the parameter $\alpha_{0}$ is constrained tightly in this scenario within our MCMC setting. The upper limit of it achieves the level of $\lesssim 10^{-3.7}$ at 90% CL. Compared with the constraint in Fig. 6, it dictates that, binary pulsar observations can lead to better constraints on $\alpha_{0}$ than GW170817. Moreover, the bound on the parameter $\beta_{0}$ becomes tighter as well. The right and upper parts of the parameter space have no support from MCMC samples.

It is worth noting that, the non-smoothness of the 1-dimensional marginalized distributions, e.g. in Fig. 7, is caused by the statistical fluctuations in the MCMC simulations. It has no statistical bearing, and will disappear if we have infinite samples. As we noted before, our samples have passed the Gelman–Rubin test for convergence, thus our limits on $\alpha_{0}$ and $\beta_{0}$ are reliable.

In Fig. 8, we illustrate the results for the scenario PSRs+GW170817. It involves the combination of the observations of five pulsars and GW170817. Basically, its result is consistent with the conclusion in Fig. 7. It is verified again, that the observation of GW170817 has little effect in constraining the DEF theory. In contrast to GW170817, the contemporary observation of radio pulsars gives us a quite powerful tool in probing the strong-field regime for gravity.

IV.4 Constraints from binary pulsars and a hypothetical BNS

In the previous subsection, we find that, compared with binary pulsars, GW170817 has little effect in constraining the DEF theory. It is interesting to investigate whether future ground-based GW detectors can surpass binary pulsar observations. We investigate the other three scenarios in Table 3, namely, PSRs+aLIGO, PSRs+CE, and PSRs+ET. We assume that, the hypothetical event, BNS@200Mpc, has component masses $(1.65\,\mathrm{M_{\odot}}$, $1.22\,\mathrm{M_{\odot}})$. One of the masses is chosen to be close to the “scalarization window” Shao et al. (2017) where, given the binary pulsar observations, there could still exist a large deviation from GR.

We investigate the power spectral density (PSD) $S_{n}(f)$ in GW detectors Abbott et al. (2017d). The SNR of a GW event is $\rho=\left(\tilde{h}(f)\right|\left.\tilde{h}(f)\right)^{1/2}$, where $\tilde{h}(f)$ is a Fourier-domain waveform, and the inner product is defined as Finn (1992),

[TABLE]

Here, the initial frequency $f_{\min}$ is chosen as the starting frequency, $f_{\rm in}$, for a GW detector, the final frequency is set to be twice of the frequency of the innermost stable circular orbit (see Ref. Shao et al. (2017)). Those starting frequencies and SNRs are listed in Table 4.

In Eq. 20, the expected 1-$\sigma$ uncertainties of $\left|\Delta\alpha\right|$, denoted as $\sigma\left(\left|\Delta\alpha\right|\right)$, are needed for different GW detectors. The technique of Fisher information matrix is used to calculate them. The Fisher information matrix is a measure of an experiment’s resolving power for the waveform parameters, collectively denoted as $\bm{\xi}$. It is defined as Finn (1992),

[TABLE]

See Appendix A in Ref. Shao et al. (2017) for more details.

For a parameter $\xi_{i}$, a lower bound on its variance expected from an experiment can be placed with the inequality of the Cramér-Rao bound Cramér (1946); Rao (1992), $\sigma\left(\xi_{i}\right)\geq\sqrt{\left(\mathcal{I}^{-1}\right)_{ii}}$. A lower bound on $\left|\Delta\alpha\right|$ can be got from the diagonal component of the inverse Fisher matrix. The corresponding 1-$\sigma$ uncertainties of $\left|\Delta\alpha\right|$ for the GW detectors are given in Table 4. The limits from the next-generation detectors will be better than that of the Advanced LIGO, due to better low-frequency sensitivities.

In the following, we give the results for the scenarios PSRs+aLIGO, PSRs+CE, and PSRs+ET.

For the scenario PSRs+aLIGO, in addition to the five binary pulsars, we use a hypothetical BNS@200Mpc, assumed to be observed by the Advanced LIGO. We illustrate the result of this scenario in Fig. 9. Compared with the scenario PSRs+GW170817 in Fig. 8, it is evident that the constraints on the DEF theory are not significantly improved.

The scenarios, PSRs+CE and PSRs+ET for the EOS AP4, are illustrated respectively in Fig. 10 and Fig. 11. The parameter space of the DEF theory is highly constrained with those next-generation detectors. We can obtain the tightest constraints on the parameters, $\alpha_{0}$ and $\beta_{0}$, in our studies. Here, the upper limit of $\left|\alpha_{0}\right|$ achieves the level of $\lesssim 10^{-4.0}$ at 90% CL, given our priors. The parameter, $-\beta_{0}$, can be constrained to $\lesssim 4.3$ at 90% CL. Due to its low-frequency sensitivity, ET can track a much longer time of a BNS inspiral than the Advanced LIGO and CE. It will constrain the DEF theory better.

In our study, we investigate all the scenarios mentioned above for nine EOSs. The upper limits of the parameters, $\left|\alpha_{0}\right|$, $-\beta_{0}$ at 90% CL are illustrated in Figs. 12 and 13 respectively. Figures 15 and 16 in Appendix A show the corresponding marginalized 1-dimensional KDE distributions.

We collect our marginalized limits in Figs. 12 and 13. In Fig. 12, for nine EOSs we adopted, the scenario GW170817 shows that $\left|\alpha_{0}\right|$ cannot be well constrained. With the five binary pulsars, the constraints are improved enormously down to the level of $\left|\alpha_{0}\right|\lesssim 2\times 10^{-4}$. It is evident that, the scenarios, PSRs, PSRs+GW170817, and PSRs+aLIGO, give similar upper bounds on $\left|\alpha_{0}\right|$. It dictates that, the constraints from GW170817 are weaker than those from binary pulsars. The constraints on $\left|\alpha_{0}\right|$ are expected to be greatly enhanced with the next-generation detectors, CE and ET. Their upper bounds on $\left|\alpha_{0}\right|$ can achieve the level of $\lesssim 1\times 10^{-4}$. Worthy to note that, the constraints on $\left|\alpha_{0}\right|$ given in Fig. 12 are heavily influenced by our priors (especially the one on $\beta_{0}$). A different prior will give a different limit, but the relative strength of these observational scenarios is fixed.

In Fig. 13, the bounds on $\beta_{0}$ from GW170817 are hardly constrained for most of the EOSs. Only for the EOS WFF1, $\beta_{0}$ can be limited to a meaningful level, $-\beta_{0}\lesssim 4.35$ at 90% CL. In contrast, for the scenarios, PSRs, PSRs+GW170817, and PSRs+aLIGO, the limits are improved up to the level of $\beta_{0}\approx-4.30$ for most of the EOSs (except for H4). When considering the scenarios PSRs+CE and PSRs+ET, we find that the bounds on $\beta_{0}$ have almost no improvement for all EOSs, except for AP4, SLy4, WFF1 and WFF2. The results of $\beta_{0}$ are different from those of $\alpha_{0}$. Particularly, the parameter $\beta_{0}$ plays the major role in controlling the strength of scalarization. We could not expect to obtain a tighter constraint on $\beta_{0}$ only by improving the precision of observations. The results in Fig. 13 are consistent with Fig. 7 in Ref. Shao et al. (2017). In particular, for the EOSs AP4, SLy4, WFF1 and WFF2, CE and ET have the potential to improve the current limit from binary pulsars.

We investigate the improvement from the scenario PSRs+ET, with respect to the scenario PSRs+GW170817. The upper bounds on the NS effective scalar coupling $\alpha_{A}$, as a function of the NS mass $m_{A}$, are illustrated in Fig. 14. Here, the upper bounds at 90% CL on theory parameters (see Figs. 12 and 13) are used for those curves. The area below the curve for each EOS corresponds to the unconstrained region for a NS. It is evident that, the allowed maximum effective scalar couplings are constrained strongly for all EOSs except for the EOS H4. The narrow range of the NS mass, that corresponds to the scalarization peak, is around $m_{A}\sim 1.9\,\mathrm{M_{\odot}}$ for the EOS H4. All NS masses from binary pulsars in Table 1, GW170817, and BNS@200Mpc are not in this region. It is also the reason why the constraints for the EOS H4 in Fig. 13 are hardly improved, with the next-generation detectors.

Combining the results in Figs. 14 and 13, we can understand some results qualitatively. In the previous discussion, the upper bounds on $\beta_{0}$ at 90% CL can be placed at the level of $\sim-4.35$ for the EOS WFF1. It is the tightest constraint in the scenario GW170817. Actually, it can be understood in Fig. 14. For the specific range of NS masses, the spontaneous scalarization is significant. For the EOS WFF1, the corresponding NS mass is around the value, $m_{A}\approx 1.4$–$1.7\,\mathrm{M_{\odot}}$. In contrast, the other EOSs allow NSs to scalarize strongly when $m_{A}\gtrsim 1.7\,\mathrm{M_{\odot}}$. One of the BNS masses derived from GW170817 is in the range $(1.36\,\mathrm{M_{\odot}},1.60\,\mathrm{M_{\odot}})$. It is within the scalarization mass range for the EOS WFF1. Therefore, solely with GW170817, we can constrain $\beta_{0}$ for this EOS.

In a short summary, the constraints on $\left|\alpha_{0}\right|$ improve with the precision of the observations. But, for $\beta_{0}$, not only the precision of the observations, but also the choice of the EOS can influence the limit. Different EOSs allow NSs to scalarize at different NS masses. We can use the observations, binary pulsars and BNSs, to constrain the DEF theory with different EOSs, if suitable systems are observed.

V Conclusion

In this paper, we studied the DEF theory in the nonperturbative strong field. In the parameter space that we investigate, NSs could develop a phenomenon called spontaneous scalarization. Instead of solving the modified TOV equations for NSs (as was done in earlier study Shao et al. (2017)), we have built efficient ROMs, implemented in the pySTGROM package, to predict NS properties. The code is made public for the community use. With the speedup of ROMs, MCMC calculations were carried out to constrain the parameter space, $(\alpha_{0},\beta_{0})$. Comparisons between various scenarios are discussed in detail. We summarize main points in the following.

To speed up large-scale calculations, pySTGROM is constructed for studying the DEF theory. Considering six scenarios that include currently available as well as future projected observations, we have tested our ROMs in practice. It turns out that, our ROMs are performed at least two orders of magnitude faster than the previous method. In the foreseeable future, we wish our ROMs to be used in relevant calculations for the DEF theory. 2. 2.

For a fixed EOS, a NS is allowed to scalarize strongly in a specific mass range of the NS. A binary pulsar and/or a BNS with a particular NS mass in this mass region may be observed. They can be used to constrain the DEF theory stringently for this EOS. 3. 3.

Using the uncertainty in $\left|\Delta\alpha\right|$ from the Fisher-matrix calculation Will (1994); Shao et al. (2017), the tightest constraints come from the scenario PSRs+ET in our studies. Given our specific priors, we can bound the parameters to be $\left|\alpha_{0}\right|\lesssim 10^{-4.0}$ and $-\beta_{0}\lesssim 4.3$ at 90% CL. Due to its low-frequency sensitivity, ET for sure will provide us with significant improvement over current constraints on the dipole radiation. 4. 4.

The spontaneous scalarization is mainly controlled by the $\beta_{0}$ parameter in the DEF theory. It is likely that, the constraints on the parameter $\left|\alpha_{0}\right|$ can be constrained more tightly with the improvement in the observational precision. Different from $\left|\alpha_{0}\right|$, the constraints on the scalarization parameter $\beta_{0}$ are related to the property of EOS, in addition to the precision of the observations.

The current and projected bounds were obtained in Sec. IV.3 and Sec. IV.4 respectively. It indicates that, the next-generation GW detectors, especially the ET, have the potential to further improve current limits, set by the binary pulsars and GW170817. Those current limits are to be improved over time, especially if suitable systems filling the “scalarization windows” are discovered Shao et al. (2017). Some binary pulsar systems, like PSRs J1012+5307 Lazaridis et al. (2009) and J1913+1102 Lazarus et al. (2016), still have large uncertainties in their masses. If the masses are constrained to be around the “scalarization window”, they may eliminate the possibility of a strong scalarization below $2\,\mathrm{M_{\odot}}$. The new large radio telescopes, such as the Five-hundred-meter Aperture Spherical radio Telescope (FAST) in China Nan et al. (2011); Li and Pan (2016); Li et al. (2019) and the Square Kilometre Array (SKA) in Australia and South Africa Kramer et al. (2004); Shao et al. (2015); Bull et al. (2018), can help to improve the timing precision. In addition, they may discover new systems that meet the requirements. On the other hand, O3 of LIGO/Virgo detectors has began in April 2019. It has been discovering a handful of GW candidates.999 GW candidates are collectively shown in the gravitational-wave candidate event database (GraceDB) by the LIGO/Virgo Collaboration in the link https://gracedb.ligo.org/latest. By the end of observing runs, some events may happen to be suitable systems to close the scalarization windows. Asymmetric BNSs (with $\alpha_{A}\neq\alpha_{B}$) as well as NS-BHs (with $\alpha_{\rm NS}\neq\alpha_{\rm BH}=0$) would give asymmetric systems and therefore could be particularly interesting. Furthermore, the next-generation ground-based GW detectors are expected to observe many more systems in the future, and they can be used to study alternative gravity theories in a more precise way. Our ROMs are built to meet the requirements of new observations to constrain the DEF theory in an efficient yet accurate way.

Acknowledgements.

We are grateful to Bin Hu and Michael Kramer for comments. We thank Norbert Wex for stimulating discussions and carefully reading the manuscript. This work was supported by the Young Elite Scientists Sponsorship Program by the China Association for Science and Technology (2018QNRC001), and was partially supported by the National Natural Science Foundation of China (11721303, 11475006, 11690023, 11622546), the Strategic Priority Research Program of the Chinese Academy of Sciences through the grant No. XDB23010200, the European Research Council (ERC) for the ERC Synergy Grant BlackHoleCam under Contract No. 610058, and the High-performance Computing Platform of Peking University. Z.C. was supported by the “Fundamental Research Funds for the Central Universities”.

Appendix A The MCMC results for nine EOSs

The marginalized 1-dimensional KDE distributions of the parameters, $\log_{10}\left|\alpha_{0}\right|$ and $-\beta_{0}$, for the EOSs in the set { AP3, AP4, ENG, H4, MPA1, PAL1, SLy4, WFF1, WFF2 }, are illustrated in Figs. 15 and 16. Their corresponding upper bounds at 90% CL are shown with dashed lines.