Constraining nonperturbative strong-field effects in scalar-tensor gravity by combining pulsar timing and laser-interferometer gravitational-wave detectors

Pulsar timing and gravitational-wave (GW) detectors are superb laboratories to study gravity theories in the strong-field regime. Here we combine those tools to test the mono-scalar-tensor theory of Damour and Esposito-Far{\`e}se (DEF), which predicts nonperturbative scalarization phenomena for neutron stars (NSs). First, applying Markov-chain Monte Carlo techniques, we use the absence of dipolar radiation in the pulsar-timing observations of five binary systems composed of a NS and a white dwarf, and eleven equations of state (EOSs) for NSs, to derive the most stringent constraints on the two free parameters of the DEF scalar-tensor theory. Since the binary-pulsar bounds depend on the NS mass and the EOS, we find that current pulsar-timing observations leave scalarization windows, i.e., regions of parameter space where scalarization can still be prominent. Then, we investigate if these scalarization windows could be closed and if pulsar-timing constraints could be improved by laser-interferometer GW detectors, when spontaneous (or dynamical) scalarization sets in during the early (or late) stages of a binary NS (BNS) evolution. For the early inspiral of a BNS carrying constant scalar charge, we employ a Fisher matrix analysis to show that Advanced LIGO can improve pulsar-timing constraints for some EOSs, and next-generation detectors, such as the Cosmic Explorer and Einstein Telescope, will be able to improve those bounds for all eleven EOSs. Using the late inspiral of a BNS, we estimate that for some of the EOSs under consideration the onset of dynamical scalarization can happen early enough to improve the constraints on the DEF parameters obtained by combining the five binary pulsars. Thus, in the near future the complementarity of pulsar timing and direct observations of GWs on the ground will be extremely valuable in probing gravity theories in the strong-field regime.


I. INTRODUCTION
In general relativity (GR), gravity is mediated solely by a rank-2 tensor, namely the spacetime metric g µν . Scalar-tensor theories of gravity, which add a scalar component to the gravitational interaction, are popular alternatives to GR. Though first proposed in 1921 [1], contemporary interest in these theories has been spurred by their potential connection to inflation and dark energy, as well as possible unified theories of quantum gravity [2]. A modern framework for the class of scalar-tensor theories we consider was developed in Refs. [3][4][5][6][7][8] (see also more generic Horndeski scalar-tensor theories in Ref. [9]).
Ultimately, the existence (or absence) of scalar degrees of freedom in gravity will be decided by experiments. Most scalar-tensor theories are designed to be metric theories of gravity, that is, they respect the Einstein equivalence principle [8,10]. Therefore, precision tests of the weak-equivalence principle, the local Lorentz invariance, and the local position invariance in flat spacetime are unable to constrain them [8,[10][11][12]. However, such theories generally violate the strong-equivalence principle. Tests of the strong-equivalence principle with self-gravitating bodies provide an ideal window to experimentally search for (or rule out) the scalar sector of gravity [8,10,13]. * lijing.shao@aei.mpg.de Particularly prominent violations of the strong-equivalence principle are known to arise in the class of massless monoscalar-tensor theories, studied by Damour and Esposito-Farèse in the form of nonperturbative strong-field effects in neutron stars (NSs) [14][15][16]. In this paper, we investigate the extent to which pulsar timing and ground-based gravitationalwave (GW) observations can constrain these phenomena (space-based GW experiments [17][18][19] are beyond the scope of this paper). Our results demonstrate that, depending on the parameters of binary systems and NS equations of state (EOSs), these two types of experiments can provide complementary bounds on scalar-tensor theories [10-12, 20, 21]. These results are especially timely as new instruments come online in the upcoming years in both fields [22,23].
The paper is organized as follows. In the next section, we briefly review two nonperturbative phenomena, notably spontaneous scalarization [14,15] and dynamical scalarization [20,[24][25][26][27], that arise in certain scalar-tensor theories of gravity. Then, in Sec. III, we derive stringent constraints on these theories by combining state-of-the-art pulsar observations of five NS-white dwarf (WD) systems. In Sec. IV, we employ these constraints and investigate the potential detectability of nonperturbative effects in binary NS (BNS) systems using the Advanced Laser Interferometer Gravitationalwave Observatory (LIGO) [28] and next-generation groundbased detectors. Finally, in Sec. V, we discuss the main results and implications of our finding, and give perspectives for future observations.

II. NONPERTURBATIVE STRONG-FIELD PHENOMENA IN SCALAR-TENSOR GRAVITY
In this paper we focus on the class of mono-scalar-tensor theories that are defined by the following action in the Einstein-frame [4,5,14,15], where G * is the bare gravitational coupling constant, g * µν is the Einstein metric with its determinant g * , R * ≡ g µν * R * µν is the Ricci scalar, ψ m collectively denotes the matter content, and A(ϕ) is the (conformal) coupling function that depends on the scalar field, ϕ. Henceforth, for simplicity, we assume that the potential, V(ϕ), is a slowly varying function that changes on scales much larger than typical length scales of the system that we consider, thus, we set V(ϕ) = 0 in our calculation.
Using a perfect-fluid description of the energy-momentum tensor for NSs in the Jordan frame, in 1993 Damour and Esposito-Farèse derived the Tolman-Oppenheimer-Volkoff (TOV) equations [14] for a NS in their scalar-tensor gravity theory. Interestingly, they discovered a phase-transition phenomenon when β 0 −4, largely irrespective of the α 0 value (a nonzero α 0 tends to smooth the phase transition [15]). The phenomenon was named spontaneous scalarization. With a suitable (α 0 , β 0 ), the "effective scalar coupling" that a NS develops, α A ≡ ∂ ln m A /∂ϕ 0 (the baryonic mass of NS is fixed while taking the derivative), could be O(1) when the NS mass, m A , is within a certain EOS-dependent range 1 . For masses 1 For black holes, the effective scalar coupling equals to zero. Therefore, the  Table I. We can clearly see a scalarization window at m A ∼ 1.7 M . below and above this range, the effective scalar coupling is much smaller 2 . In Fig. 1 we show an example of spontaneous scalarization for a NS with the realistic EOS SLy4, and compare it to existing individual binary-pulsar constraints.
In general, if two compact bodies in a binary have effective scalar couplings, α A and α B , they produce gravitational dipolar radiation ∝ (∆α) 2 , with ∆α ≡ α A − α B , which is at a lower post-Newtonian (PN) order than the canonical quadrupolar radiation in GR [15] 3 . In Ref. [16], Damour and Esposito-Farèse for the first time compared limits on the DEF gravity arising from Solar system and binary pulsar experiments with expected limits from ground-based GW detectors like LIGO and Virgo. The analysis in Ref. [16] is based on soft (by now excluded [31,32]), medium and stiff EOSs, and for the LIGO/Virgo experiment it assumes a BNS merger with PSR B1913+16 like masses (1.44 M and 1.39 M ), as well as a 1.4 M -10 M NS-BH merger. Damour and Esposito-Farèse come to the conclusion that binary-pulsar experiments would generally be expected to put more stringent constraints tests performed with binary black holes [30] do not directly apply to the DEF theory. 2 For sufficiently negative β 0 ( −4.6), NSs do not de-scalarize before reaching their maximum mass, i.e. spontaneous scalarization is found for all NSs above a certain critical mass, which depends on the actual value of β 0 and the EOS [14,15]. 3 In this paper, generally we denote with nPN the O(v 2n /c 2n ) corrections to the leading Newtonian dynamics (equations of motion). Therefore, the gravitational dipolar radiation reaction is at 1.5 PN, and the quadrupolar radiation is at 2.5 PN. In the GW phasing, when there is no potential confusion we sometimes refer to the quadrupolar (dipolar) radiation as 0 PN (−1 PN), as typically done in the literature. on the parameters (α 0 , β 0 ) than ground-based detectors, such as LIGO and Virgo. Since then, several analyses have followed [32][33][34][35][36], but typically those studies did not probe a large set of NS masses and EOSs. Considering advances in our knowledge of NSs and more sensitive current and future ground-based detectors, we revisit this study here. Quite interestingly, as pointed out in a first study in Ref. [20], the constraints on the parameters (α 0 , β 0 ) from binary pulsars depend quite crucially on the EOSs and the masses of the NSs, in particular in the parameter space that allows for spontaneous scalarization. By taking into account this dependence when setting bounds from pulsar timing, we shall find that current and future GW detectors on the Earth might still be able to exclude certain specific regions of the parameter space (α 0 , β 0 ) that are not probed by binary pulsars yet.
Twenty years after the discovery of spontaneous scalarization, Barausse et al. [24] found another interesting nonperturbative phenomenon in a certain parameter space of the DEF theory. This time the scalarization does not take place for a NS in isolation, but for NSs in a merging binary. Indeed, modeling the BNS evolution in numerical relativity, Barausse et al. found that the two NSs can scalarize even if initially, at large separation, they are not scalarized. This phenomenon is called dynamical scalarization, and its onset is determined by the binary compactness instead of the NS compactness.
Reference [24] also demonstrated that a spontaneously scalarized NS can generate scalar hair on its initially unscalarized NS companion in a binary system through a process known as induced scalarization. Dynamical and induced scalarization cause BNSs to merge earlier [20,24,26] than in GR, resulting in a significant modification to the GW phasing that is potentially detectable by ground-based GW interferometers [24,25,27,37].
Finally, it is important to note that cosmological solutions in the strictly massless limit of the DEF theory are known to evolve away from GR when β 0 is negative [37][38][39][40]; to be consistent with current Solar-system observations, such cosmologies require significant fine tuning of initial conditions. 4 Various modifications to the theory have been considered to cure this fine-tuning problem, e.g. adding higher order polynomial terms to ln A(ϕ) [40] or including a mass term V(ϕ) = 2m ϕ ϕ 2 in the action [15,[42][43][44]. To date, none of these proposals have produced scalar-tensor theories that: (i) satisfy cosmological and weak-field gravity constraints, (ii) generate an asymptotic field ϕ 0 that is stable over timescales relevant to binary pulsars and GW sources, and (iii) give rise to the nonperturbative phenomena present in DEF theory. As is commonly done in the literature [7,14,15,20,[24][25][26][27]37], we will ignore these cosmological concerns in this paper and focus only on (massless) DEF theory. The mass-radius relation of NSs (in GR) for the 11 EOSs that are adopted in the study. The mass constraint (with 1-σ uncertainty) from PSR J0348+0432 [32] is depicted in grey. The color coding for different EOSs is kept consistent for all figures in this paper.

III. CONSTRAINTS FROM BINARY PULSARS
Until now, binary pulsars have provided the most stringent limits to the DEF theory [15,32,35,36,47]. These limits were usually obtained with individual pulsar systems and with representative EOSs [36] 5 . Here, by contrast, we combine observational results from multiple pulsar systems employing Markov-chain Monte Carlo (MCMC) simulations [53]. In particular, we pick the five NS-WD binaries that are the most constraining systems in testing spontaneous scalarization: PSRs J0348+0432 [32], J1012+5307 [45], J1738+0333 [35], J1909−3744 [46], and J2222−0137 [47]. We choose these five binaries basing on the binary nature (namely, NS-WD binaries), the timing precision that has been achieved, and the NS masses. These aspects are important to the study here, and see Refs. [11,12,52] for more discussion. For convenience, we list the parameters of these binaries in Table I, and notice that it is the combination of theirṖ int b and the NS mass that makes them particularly suitable for the test of spontaneous scalarization. We obtain the limits using 11 different EOSs that have the maximum NS mass above 2 M [54]. The names of these EOSs are AP3, AP4, ENG, H4, MPA1, MS0, MS2, PAL1, SLy4, WFF1, and WFF2 (see Refs. [54,55] for reviews). Figure 2 shows the mass-radius relation of NSs in GR for these EOSs. As evidenced by the spread of the curves in the figure, we believe that these EOSs are sufficient to cover the different EOS-dependent properties of spontaneous scalarization, and at the same time satisfy the two-solar-mass limit from pulsar-timing observations [31,32].
Markov-chain Monte Carlo techniques allow us to preform  [32,35,[45][46][47]. The observed time derivatives of the orbit period P b are corrected using the latest Galactic potential of Ref. [48]. For PSRs J0348+0432, J1012+5307 and J1738+0333, the mass ratios were obtained combining radio timing and optical high-resolution spectroscopy, while the companion masses are determined from the Balmer lines of the WD spectra based on WD models [32,[49][50][51]. For PSRs J1909−3744 and J2222−0137, the masses were calculated from the Shapiro delay, where the range of the Shapiro delay gives directly the companion mass, and the pulsar mass is then being derived from the mass function, using the shape of the observed Shapiro delay to determine the orbital inclination [46,47]. The masses below are based on GR as the underlying gravity theory. However, since the companion WD is a weakly self-gravitating body, they are practically the same in the DEF theory (with a difference 10 −4 ). We give in parentheses the standard 1-σ errors in units of the least significant digit(s).  (25) parameter estimation within the Bayesian framework. These methods provide the posterior distributions of the underlying parameters that are consistent with observations. In Bayesian analysis, given data D and a hypothesis H (here, the DEF theory), the marginalized posterior distribution of (α 0 , β 0 ) is given by [53], where I is all other relevant prior background knowledge and Ξ collectively denotes all other unknown parameters besides (α 0 , β 0 ), which are marginalized over to obtain the marginalized posterior distributions for just (α 0 , β 0 ) [see below for more details]. In the above equation, given H and I, P (α 0 , β 0 , Ξ|H, I) is the prior on (α 0 , β 0 , Ξ), P (D|α 0 , β 0 , Ξ, H, I) ≡ L is the likelihood, and P (D|H, I) is the model evidence. As said, we use MCMC techniques to explore the posterior in Eq. (4). We discuss below our choices for the priors and the likelihood function [see Eq. (9)]. We assume that observations with different binary pulsars are independent.
We now explain how we employ the MCMC technique to get the posterior by combining binary pulsar systems. Let us assume that N pulsars (N = 1, 2, 5, see below) are used to constrain the (α 0 , β 0 ) parameter space. To obtain a complete description of the gravitational dipolar radiation of these systems in the DEF theory, we need N + 2 free parameters in the MCMC runs, which are θ = α 0 , β 0 ,ρ (i) c , whereρ (i) c (i = 1, · · · , N) is the Jordan-frame central matter density of pulsar i 6 . As an initial value to the TOV solver, we also need 6 Actually, in the full calculation we need some other quantities, as well, for example, the orbital period, P b , and the orbital eccentricity, e. Those quantities are observationally very well determined (see Table I), thus we use their central values and find that our constraints on (α 0 , β 0 ) do not change on a relevant scale when we take into account the errors on those quantities. the value of the scalar field in the center of a NS, ϕ (i) c , but the latter is fixed iteratively by requiring that all pulsars have a common asymptotic value of ϕ, ϕ 0 ≡ α 0 /β 0 . Givenρ (i) c and ϕ (i) c for pulsar i, we integrate the modified TOV equations (see Eq. (7) in Ref. [14] or Eq. (3.6) in Ref. [15]) with initial conditions given by Eq. (3.14) in Ref. [15]. During the integration, we use tabulated data of EOSs, and linearly interpolate them in the logarithmic space of the matter density,ρ, the pressure, p, and the number density,ñ [54]. Note that only one quantity among {ρ,p,ñ} is free, while the others are determined by the EOS. The end products of the integration provide us, for each pulsar, the gravitational mass, m (i) A , the baryonic mass,m (i) A , the NS radius, R (i) , and the effective scalar coupling, α (i) A [15]. For the MCMC runs we use a uniform prior on log 10 |α 0 | for |α 0 | ∈ [5 × 10 −6 , 3.4 × 10 −3 ], where 3.4 × 10 −3 is the limit obtained from the Cassini spacecraft [29,56]. We pick the parameter β 0 uniformly in the range [−5, −4], which corresponds to a sufficiently large parameter space where the scalarization phenomena can take place [14,24]. During the exploration of the parameter space, we restrict the values of (α 0 , β 0 ) to this rectangle region, as well, in order to avoid overusing computational time in uninteresting regions. The initial central matter densities, ρ (i) c , are picked around their GR values, but they are allowed to explore a very large range in the simulations. As we discuss below, we perform convergence tests and verify that when evolving the chains all parameters in θ quickly lose memory of their initial values.
During the MCMC runs, we evolve the N + 2 free parameters according to an affine-invariant ensemble sampler, which was implemented in the emcee package [57,58] 7 . At every step, we solve the N sets of modified TOV equations on the fly, using for the companion masses of the binary pulsars the values listed in Table I 8 . Then, for the decay of the binary's orbital period, which enters the likelihood function [see Eq. (9)], we use the dipolar contribution from the scalar field and the quadrupolar contribution from the tensor field as given by the following, well known, formulae [7,59], with f (e) ≡ 1 + 73 24 We find that higher order terms, as well as the subdominant scalar quadrupolar radiation, give negligible contributions to this study. Notice that in Eq. (5), we have replaced the effective scalar coupling of the WD companion with the linear matter-scalar coupling constant, since for a weakly selfgravitating WD α A α 0 in the β 0 range of interest. 9 Furthermore, we can approximate the bare gravitational constant G * in the above equations with the Newtonian gravitational constant G N = G * (1 + α 2 0 ), since |α 0 | 1 (e.g., from the Cassini spacecraft [29,56]).
We construct the logarithmic likelihood for the MCMC runs as, where for PSRs J1909−3744 and J2222−0137 we replace the second term in the squared brackets with m p − m obs p /σ obs m p 2 . 7 http://dan.iel.fm/emcee 8 The masses in PSRs J0348+0432, J1012+5307, and J1738+0333 are based on a combination of radio timing of the pulsars and optical spectroscopic observation of the WDs. The derivation of the masses only depends on the well-understood WD atmosphere model in combination with gravity at Newtonian order, and the mass ratio q, which is free of any explicit strongfield effects [29]. Therefore, even within the DEF theory, these masses are valid [32,35]. For PSRs J1909−3744 and J2222−0137, the masses are derived from the range and shape of the Shapiro delay [15,36]. Since for the weakly self-gravitating WD companion |α B | ≈ |α 0 | 1, these masses are practically identical to the GR masses in Table I. 9 In this context, see footnote "d" in Ref. [33], concerning WDs and very large (positive) β 0 . 4.8 The marginalized 2-d distribution of (log 10 |α 0 |, −β 0 ) from MCMC runs on the five pulsars listed in Table I,  In Eq. (9), the predicted orbital decay from the theory isṖ th , and σ obs X is the observational uncertainty for X ∈ Ṗ int b , q, m p , as given in Table I. Note thatṖ th b and m p implicitly depend on (α 0 , β 0 , Ξ), through direct integration of TOV equations in the DEF theory.
We pick J0348 and J1738 due to their mass difference (2.01 M and 1.46 M respectively), and their high timing precision (see Table I), which leads to interesting differences in the constraints on the DEF parameters, especially on β 0 . For each run, we accumulate sufficient MCMC samples to guarantee the convergence of MCMC runs. By using the Gelman-Rubin statistic [60], we find that, for each EOS, 200,000 samples for cases J0348 and J1738, and 400,000 samples for cases 2PSRs and 5PSRs, are enough, respectively. We discard the first half chain points of these 44 runs (4 cases × 11 EOSs) as the burn-in phase [58,61], while we use the remaining samples to do inference on the parameters (α 0 , β 0 ).
As an example, we show in Fig. 3 the marginalized 2-d distribution in the parameter space of (log 10 |α 0 |, −β 0 ) for the case TABLE II. Limits on the parameters of the massless mono-scalar-tensor DEF theory for different EOSs when applying the MCMC analysis to the five pulsars J0348+0432, J1012+5307, J1738+0333, J1909−3744, and J2222−0137. Results at 68% and 90% CLs are listed. |α A | max is the maximum effective scalar coupling that a NS could still possess without violating the limits, and m max A is the corresponding (gravitational) mass at this maximum effective scalar coupling (see Figure 5).  5PSRs and the EOS SLy4. As mentioned above, we distribute the initial values of log 10 |α 0 | and −β 0 uniformly in the rectangle region of Fig. 3. We see that after MCMC simulations, the region with large |α 0 | or large (negative) β 0 is no longer populated, and only a small corner is consistent with the observational results of the five NS-WD binary pulsars.
Furthermore, we extract the upper limits of log 10 |α 0 | and −β 0 from their marginalized 1-d distributions at 68% and 90% CLs. We summarize the upper limits at 90% CL from all 44 runs in Fig. 4. It is interesting to observe the following facts. First, for all EOSs, J1738 gives a more constraining limit on α 0 than J0348. This result might be due to the fact that the σ obṡ of J1738 is about two orders of magnitude smaller than that of J0348, thus giving a tighter limit on α 2 0 by roughly the same order of magnitude. Second, the constraints on β 0 from J0348 and J1738 are extremely EOS-dependent. This should be a consequence of the masses of the NSs, which are (in GR) 1.46 M for J1738, and 2.01 M for J0348. For EOSs that favour spontaneous scalarization at around 1.4-1.5 M , J1738 gives a better limit, while for EOSs that favour spontaneous scalarization at around 2 M , J0348 gives a better limit. This trend is also consistent with Fig. 5 (to be introduced below). Third, by combining two pulsars (2PSRs), NSs are limited to scalarize at neither 1.4-1.5 M nor ∼ 2 M . Therefore, almost for all EOSs, β 0 is well constrained. This result demonstrates the power of properly using multiple pulsars with different NS masses to constrain the DEF parameter space for any EOS. Fourth, we obtain the most stringent constraints with five pulsars (5PSRs). This is especially true for β 0 , which is constrained at the level of ∼ −4.2 (68% CL) and ∼ −4.3 (90% CL) for all EOSs. Finally, we list in Table II the marginalized 1-d limits for 5PSRs. We shall use them in the next section when combining binary pulsars with laser-interferometer GW observations.
Considering the results that we have obtained when combining the five pulsars (5PSRs), one could wonder whether isolated NSs can still be strongly scalarized. To address this question, we use the limits on (α 0 , β 0 ) and calculate the effective scalar coupling that a NS can still develop as a function of the NS mass, for the 11 EOSs used in this paper. The results at 90% CL are summarized in Fig. 5, while in Table II we list the maximally allowed effective scalar couplings at 68% and 90% CLs, and their corresponding (gravitational) NS masses (marked as dots in Fig. 5).  Table II). Furthermore, quite remarkably Fig. 5 shows that there are scalarization windows (this feature could also be seen in Fig. 1 Table II). The point of the maximum |α A | is marked with a dot, and the values (and the corresponding masses) are listed in Table II. For these specific masses, using the 11 EOSs that can give rise to spontaneous scalarization, we have constrained stringently the DEF parameters. However, some EOSs can still allow NS to scalarize strongly (i.e., acquire large effective scalar couplings) for other values of the masses. As Fig. 5  Those scalarization windows could be closed in the future if binary pulsars with these masses are discovered and their gravitational dipolar radiation is constrained by pulsar timing. As we shall discuss in Sec. IV, the presence of scalarization windows also opens the interesting possibility to close these gaps with future GW observations from BNSs, if the NS's masses lie in the scalarization window.

IV. PROJECTED SENSITIVITIES FOR LASER-INTERFEROMETER GRAVITATIONAL-WAVE DETECTORS
Having determined constraints on the DEF's parameter space from binary pulsars (see Table II) and found scalarization windows, we now address the question of whether present and future laser-interferometer GW observations on the ground can still improve these limits and close the gaps. Two scenarios are considered: (i) asymmetric BNS systems, equipped with separation-independent effective scalar couplings, whose gravitational dipolar radiation during the inspiral can modify the GW phasing [16,62], and (ii) BNS systems that dynamically develop scalarization during the late stage of the inspiral, leading to significant, nonperturbative changes in the GW signal [20,24,26,27]. A complete description of BNSs in the DEF theory should include both effects. However, complete waveform models from the theory are still not available, so here we investigate the two scenarios separately to obtain some conservative understanding of the whole picture.

A. Dipole radiation for binary neutron-star inspirals
The presence of a scalar field can significantly modify the inspiral of an asymmetric BNS system, due to the additional energy radiated off by the scalar degree of freedom. The most prominent effect is a modification of the phase evolution in GW signals. For two NSs with effective scalar couplings α A and α B respectively, one finds for the evolution of the orbital frequency, Ω, up to 2.5 PN order [62][63][64], where ∆α ≡ α A − α B , η ≡ m A m B /M 2 , and the (dimensionless) "characteristic" velocity, The quantity κ is given in Refs. [7,32]. In GR one has α A = α B = 0 and κ = 1. Note that there is also a subdominant contribution from the scalar quadrupolar waves at 2.5 PN order, which however can be absorbed by a 1% change in the mass parameters. Here we assume that α A and α B are constant during the inspiral, and their values are obtained from isolated NSs. This assumption is valid as long as the induced or dynamical scalarization mechanisms are not triggered.
For an asymmetric compact binary where α A α B , the most prominent deviation from the GR phase evolution is determined by the dipole term in Eq. (10), i.e., the contribution ∝ V 3 . To leading order, the offset in the number of GW cycles in band until merger due to the dipole term is given by where V in corresponds to V in Eq. (11) when the merging system enters the band of the GW detector, i.e., when Ω = π f in (see Refs. [16,62] for details). In the above equation we have used the approximation κ 1 and the fact that V in is much smaller than V just before merger. Within the approximation of Eq. (12) one can use V in (G N Mπ f in ) 1/3 /c, i.e., replacing the effective gravitational constant G * (1 + α A α B ) in Eq. (11) by the Newtonian gravitational constant G N ≡ G * 1 + α 2 0 [7,15]. Again, we stress that Eq. (12) is based on the assumption that the effective scalar couplings of the two NSs,  ) for frequencies f > f in , and its change due to the dipole radiation in the DEF's theory, ∆N dipole , assuming |∆α| ≡ |α A − α B | 0.0199, which comes from the maximally allowed effective scalar couplings for the EOS SLy4 at 90% CL (see Figure 5). The limits on the contributions from leading-order spin-orbit and spin-spin terms, ∆N β and |∆N σ |, are listed where the (dimensionless) spins of the double pulsar (when it merges in 86 Myr) are used. For ∆N β and |∆N σ |, we also give in parentheses when both NSs are spinning at the maximal spin that we have ever observed (in an eclipsing binary pulsar J1748−2446ad).
To obtain a rough understanding of the effects of dipolar radiation, let us calculate the dephasing from GR by an asymmetric BNS inspiral with m A = 1.25 M and m B = 1.7 M . According to Fig. 5, at present, binary-pulsar experiments cannot exclude |α A | as large as 10 −2 -10 −1 for NSs of a certain mass range, which depends on the EOS. For the EOS SLy4 we find from the corresponding (dark green) curve in Fig. 5 |α A | 0.0007 and |α B | 0.0206, hence |∆α| ≡ |α A − α B | 0.0199.
In our study we consider the Advanced LIGO (aLIGO) detectors at design sensitivity [28], and future ground-based detectors, such as the Cosmic Explorer (CE), and the Einstein Telescope (ET). We use the starting frequencies, f in = 10 Hz for aLIGO, f in = 5 Hz for CE, and f in = 1 Hz for ET [23,65,66]. In Table III we list the number of GW cycles as predicted by GR, N GR , and the change in the number of cycles caused by the existence of a dipole radiation for a BNS signal, ∆N dipole . From Table III we can already see that, given the BNS parameters, current bounds by pulsars still leave room for significant time-domain phasing modifications in BNS mergers, in particular if one of the NSs falls into the scalarization window of ∼ 1.7 M (for EOSs AP4, SLy4, and WFF2) to ∼ 1.9 M (for the EOS H4), or if one NS's mass significantly exceeds 2 M (for EOSs MS0 and MS2). As reference points, we list in Table III also the changes in the number of GW cycles from spin-orbit and spin-spin effects. Indeed, from the leading-order spin-orbit (1.5 PN) and spinspin (2 PN) contributions to the GW phasing [18,67], one has where The (dimensionless) spins of a BNS system, χ A and χ B , are likely to be small in magnitude. The parameters β and σ are maximized when two spins are aligned with the direction of the orbital angular momentum,L. The limits on ∆N β and |∆N σ | are listed in Table III where we have used the (dimensionless) spins of the double pulsar system that is the only double NS system where two spins are precisely measured. When the double pulsar evolves to the time of its merger in 86 Myr from now, one has |χ A | 0.014 and |χ B | 0.00002 [68], assuming a canonical moment of inertia 10 38 kg m 2 for NSs. As we can see from Table III, if the spins of BNSs to be discovered by GW detectors are comparable to that of the double pulsar, the inclusion of spins only affects the number of GW cycles at percentage level at most. In addition, because ground-based detectors could observe BNSs from a population different from the one observed with pulsar timing, we also give ∆N β and |∆N σ | in Table III when the (dimensionless) spin of the fastest rotating pulsar ever observed, PSR J1748−2446ad (P = 1.4 ms) [69], is used for both NSs 10 . Even in this extreme case with |χ A | = |χ B | 0.26 (assuming a canonical mass 1.4 M ), ∆N dipole is still larger (or comparable in the case of the Advanced LIGO) than the upper limits of ∆N β and |∆N σ |.
The dephasing quantity, ∆N dipole , is nevertheless a crude indicator for realistic detectability. In reality, one has to consider various degeneracy between binary parameters, the waveform templates that are used for detection and parameter estimation, the power spectral density (PSD) of noises in GW detectors, S n ( f ), the signal-to-noise ratio (SNR) of an event ρ, and so on. In order to obtain more quantitative estimates of the constraints on dipolar radiation that can be expected from GW detectors, one would need to compute Bayes factors between two alternatives [70] or apply cutting-edge parameter-estimation techniques, for example the MCMC or nested sampling [71]. However, given the limited scope of our analysis, for simplicity, here, we adopt the Fisher-matrix approach [17,62,72,73], although we are aware of the fact that for events with mild SNR (ρ ∼ 10), the Fisher matrix can have pitfalls [74]. In Appendix A we review the main Fisher-matrix tools that we use. We summarize in Fig. 6 their dimensionless noise spectral density f S n ( f ) [23,65,66], and show in the figure also an hypothetical BNS signal. For all the studies, we fix the luminosity distance to D L = 200 Mpc for aLIGO, CE, and ET. Indeed, within such a distance, aLIGO alone is supposed to observe 0.2-200 BNS events annually at design sensitivity [75]. With the four-site network incorporating LIGO-India at design sensitivity, the number of detectable BNS events will double [75]. Therefore, it is a realistic setting to discuss BNS events for aLIGO; to be conservative, we only consider a two-detector network for aLIGO in our study. CE and ET have better sensitivities, thus will have larger SNRs for these events; besides, they will be able to detect a larger number of BNSs, including those with unfavourable orientations. Using the standard cosmological model [76], the redshift associated to D L = 200 Mpc is z 0.0438, and we take it into account in our Fisher-matrix calculation, even if it generates a small effect. Moreover, we always report masses in the rest frame of a BNS system.
The Fisher matrix is constructed as usual from the Fourierdomain waveformh( f ) [62,72,73], with ∂ ah ( f ) being the partial derivative to the parameter labeled "a" (see Appendix A for definitions and notations). We use the waveform parameters ln A, ln η, ln M, t c , Φ c , (∆α) 2 to construct the 6 × 6 Fisher matrix, Γ ab . The inverse of the Fisher matrix is the correlation matrix for these parameters, from where we can read their uncertainties and correlations [17,62,72,73]. In Fig. 7 we plot in dashed lines the uncertainties in |∆α| obtained with three GW detectors (aLIGO, CE, and ET) for an asymmetric BNS with rest-frame masses m A = 1.25 M and m B > 1.25 M , located at D L = 200 Mpc. For a BNS of masses, for example, (1.25 M , 1.63 M ) which are the most probable masses for the newly discovered asymmetric double-NS binary pulsar PSR J1913+1102 [77], we find that aLIGO, CE, and ET can detect its merger at 200 Mpc with ρ = 10.6, 450, and 153, respectively, after averaging over pattern functions and assuming two detectors in each case. The characteristic strain of such a BNS is illustrated in the frequency domain in Fig. 6. In the large SNR limit, the uncertainties in |∆α| scale with the SNR as ρ −1/2 . In Fig. 8, we give the correlations between parameters obtained from the matched-filter analysis. We find that due to its low-frequency sensitivity ET can break some degeneracy between parameters better than aLIGO and CE do.
In Fig. 7, we show with solid lines the maximum values of |∆α| at 90% CL from pulsars for 11 EOSs (calculated from Fig. 5). If for some NS's mass range a solid line (which is associated to a certain EOS) is above a dashed line (associated to a certain detector), then for NSs described by that EOS, the corresponding GW detector has potential to further improve • aLIGO has potential to further improve the current limits from binary pulsars with a discovery of a BNS of suitable masses, if the EOS of NSs is one of (or similar to) H4, MS0, MS2, SLy4, and WFF2; • CE and ET, due to their low-frequency sensitivity and better PSD curves, are able to significantly improve current limits from binary pulsars on the DEF's parameters, no matter what the real EOS of NSs is.
We stress that those conclusions are obtained with a Fishermatrix analysis, and should be made more robust in the future using more sophisticated tools, notably Bayesian analysis.

Constraints outside the spontaneous-scalarization regime
With the results above it is fairly straightforward to calculate the limits from aLIGO, CE and ET on |α 0 | when β 0 is outside the spontaneous scalarization regime, i.e., β 0 −4, and compare them to existing limits from the Solar system and pulsars [16]. For completeness, we present the relevant results here. Shibata et al. [20] have shown that for small α 0 there exists a simple relation between α A , α 0 and m A as long as spontaneous scalarization does not set in (see Eq. (44) in Ref. [20]), which in our notation reads 11 11 In principle, there is still a weak dependence on α 0 in A (A) β 0 . However, this dependence becomes very small for |α 0 | 10 −2 , as it scales with α 2 0 . Therefore the α 0 -dependence is absolutely negligible for the parameter space explored here.
With this equation at hand one can directly convert the limits from ground-based GW detectors of Fig. 7, for any given β 0 −4, into limits for |α 0 | via Figure 9 gives the results for two different mass configurations and the EOS AP4. A more stiff EOS would generally lead to less constraining limits for ground-based GW detectors and binary pulsars. As one can see, in the range β 0 −4 current Solar system and pulsar tests are already clearly more constraining than what aLIGO is expected to obtain. For CE and ET, only inspirals with a very massive component will provide constraints that are better than present limits, for a limited range of β 0 (see also aLIGO [78] and ET [78,79] limits from a NS-BH inspiral for the special case of β 0 = 0, i.e. the Jordan-Fierz-Brans-Dicke gravity). By the time CE or ET is operational, however, the expected limits from GAIA [80] and SKA [11] will have left little room for ground-based GW observatories in the regime. The space-based GW observatories LISA [17,18] and DECIGO/BBO [19] could in principle also provide limits on the DEF theory from an inspiral of a NS into an intermediate mass black hole, provided such BHs exist. However, the resulting limits on |α 0 | are not expected to be better than limits from future ground-based GW observatories [17]. It is worthy to mention that for very large (positive) β 0 , say, β 0 10 2 -10 3 , massive NSs might develop instabilities [81,82], which is beyond the scope of Fig. 9.

B. Dynamical scalarization
In addition to the nonlinear gravitational self-interaction testable with binary pulsars, GW detectors probe the nonlinear interactions between coalescing NSs. Dynamical scalarization stems from the interplay between these two regimes of For comparison we have plotted Solar system limits (grey) and the limits from PSR J1738+0333 (magenta), which currently gives the best limit for β 0 3. The limit from Cassini [56] and the limit expected from GAIA [80] are also shown. strong gravity and thus offers a promising means of complementing pulsar timing constraints on scalar-tensor theories.
Numerical relativity simulations have demonstrated that dynamical scalarization can significantly alter the late-time behavior of a BNS system. If this transition occurs before merger, the sudden growth of effective scalar couplings impacts the system's gravitational binding energy and energy flux so as to shorten the time to merger [20,24,26].
The prospective detectability of this effect was investigated in Refs. [37,83] using Bayesian model selection. The authors sought to recover injected inspiral waveforms containing dynamical scalarization with template banks constructed from similar waveforms. The injected signals and template banks used PN waveforms augmented with various non-analytic models of dynamical scalarization. To mimic the abrupt activation of the dipole emission at the onset of dynamical scalar-ization, Ref. [83] added a −1 PN correction modulated by a Heaviside function to a GR waveform, i.e., signals of the form where Ψ −1 PN ( f ) = b f −7/3 , and b and f * are parameters of the model. Injected signals were recovered with a template bank of waveforms of the same form. In Ref. [37], the authors injected waveforms constructed in Ref. [25] by integrating the 2.5 PN equations of motion combined with a semianalytic model of scalarization, then performed parameter estimation using both templates that included −1 PN and 0 PN scalar-tensor effects throughout the entire inspiral and those that modeled their sudden activation as in Ref. [83] Combined, these analyses provide a loose criterion for whether a dynamically scalarizing BNS system could be distinguished from the corresponding system in GR by aLIGO.
The key characteristic of such systems is the frequency f DS at which dynamical scalarization occurs. To be distinguishable from a GR waveform, a significant portion of the dynamically scalarized signal's SNR must occur after f DS , or equivalently, f DS must be sufficiently lower than the merger frequency. Using waveforms of the form of Eq. (20), the authors found in Ref. [83] that dynamical scalarization can only be observed with aLIGO if f DS 50-100 Hz. In only one injection considered in Ref. [37] was dynamical scalarization detectable, occurring at f DS ≈ 80 Hz. Understandably, these analyses rely on some initial assumptions that may bias these estimates away from the real detectability criteria, such as the limited range of masses and EOS considered and ignoring any degeneracies introduced by the merger and ringdown portions of the waveform or by the inclusion of spins. Ignoring these subtleties for the moment, we investigate whether the pulsartiming constraints described in Sec. III can exclude the possibility of observing dynamical scalarization with aLIGO using the conservative detectability criterion from Refs. [37,83] that scalarization must occur by f DS 50 Hz.
We consider binary systems composed of NSs with masses ranging from 1.3 M to 1.9 M . We compute f DS within the "post-Dickean" (PD) framework, a resummation of the PN expansion formulated in Ref. [27]. This model introduces new dynamical degrees of freedom that capture the nonperturbative growth of the scalar field using a semi-analytic feedback loop. This approach provides a mathematically consistent backing to previous models of dynamical scalarization [25]. The model incorporates a certain flexibility in the choice of resummed quantities; we adopt the m (RE) , F (φ) scheme outlined in Table I of Ref. [27] because it was found to give the best agreement with numerical computations of quasiequilibrium configurations [26]. For clarity, we dress quantities defined in the PD framework with tildes and leave quantities defined in the PN framework unadorned; in the limit that no resummation is performed, the PD quantities reduce to their PN analogs.
Within the PD framework, the effective scalar coupling of each NS is promoted to a function of both the asymptotic scalar field ϕ 0 and the local scalar field in which the body is immersed, i.e.α A =α A (ϕ 0 , ϕ A ). Unlike in the PN treatment, this coupling evolves as the BNS coalesces. Similarly, the inertial mass of each bodym A (ϕ 0 , ϕ A ) evolves in the PD framework. However this mass varies by no more than 0.01%, so in practice, one can simply use the PN mass m A in place ofm A .
We define the mass-averaged scalar coupling of the system asᾱ where precise definitions ofm A andα A are given in Eqs. (A3) and (A4) in Ref. [27]. Note that for equal-mass binaries, we haveα A =α B =ᾱ.
Following the work of Ref. [27], we compute the massaveraged scalar coupling as a function of frequency for binaries on quasi-circular orbits to 1 PD order. The average scalar coupling is plotted in Fig. 10 for equal-mass binaries for theories that saturate the pulsar-timing constraints at 90% CLs, as given in Table II. Scalarization occurs earlier for larger mass systems, with an ordering (by EOS) determined by the magnitude of To compute this quantity, one takes the difference in effective scalar couplings of NSs (of equal baryonic mass) with infinitesimally different asymptotic scalar fields ϕ 0 ; however, for non-spontaneously scalarized stars, β A is given approximately by provided that |α A | is sufficiently small. Binaries with spontaneously scalarized stars begin with an appreciable effective scalar coupling at large separations that continues to grow as they coalesce. In light of this remark, we note that there is no observational distinction between spontaneous (or induced) scalarization and dynamical scalarization that occurs at sufficiently low frequencies; for example, compare the scalarization of (1.7 M , 1.7 M ) systems composed of NSs with the EOSs SLy4, AP4, and WFF1 (the dark green, red, and beige curves in the lower left panel of Fig. 10, respectively). The sharp feature for the WFF1 EOS in the (1.9 M , 1.9 M ) system occurs because of the relatively low mass at which spontaneous scalarization occurs for this particular EOS. We provide a more detailed analysis of this phenomenon in Appendix B. Similarly abrupt transitions occur for other EOSs in more massive binary systems with individual masses 2 M .
We adopt the method introduced in Ref. [26] to extract f DS . The average effective scalar coupling can be closely fit by the piecewise function where Θ is the Heaviside function, a 1 and x DS are fitting parameters, and x ≡ (G * MΩ/c 3 ) 2/3 . In practice, we identify x DS with the peak in the second derivative of the lefthand side . We use the limits on (α 0 , β 0 ) at 90% CLs, given in Table II, for each EOS. The corresponding GW frequency is given along the top axis, with f GW = Ω/π. Dashed vertical lines highlight the conservative detectability criterion for aLIGO that f DS 50 Hz, derived in Refs. [37,83].
of Eq. (24) with respect to x. The gravitational wave frequency at which dynamical scalarization occurs is then given by f DS = Ω DS /π. In Ref. [27], the PD prediction was found to reproduce numerical-relativity results to within an error of 10% with this fitting procedure. The dynamical scalarization frequencies for the configurations considered in Fig. 10 are given in Table IV for theories constrained at the 68% and 90% CLs. Systems containing spontaneously scalarized stars (i.e., those with appreciable effective scalar coupling even in isolation) are demarcated as scalarizing below 1 Hz; as noted above, these systems would be indistinguishable to GW detectors from those that dynamically scalarize below 1 Hz. For clarity, we highlight the systems in Table IV that scalarize (dynamically or spontaneously) below 50 Hz. Recall that, under our definition, induced scalarization occurs in binaries comprised of one initially scalarized star and one initially unscalarized star; this asymmetry cannot be achieved in equal-mass systems like those discussed above.
We next consider the onset of dynamical scalarization in unequal-mass systems. For the sake of compactness, we show in Table V the dynamical scalarization frequencies for binaries with NS masses of 1.2 M to 1.9 M with just the MPA1 EOS. We find that the total mass plays a more important role in determining the onset of dynamical scalarization than the mass ratio. Fixing the total mass, we find that scalarization occurs earlier in more asymmetric binaries of lower mass (e.g., M 3.2M for the MPA1 EOS). None of the systems listed in Table V undergo induced scalarization. As before, we highlight the systems in Table V that scalarize below 50 Hz. To summarize, Tables IV and V demonstrate that binarypulsar constraints cannot entirely rule out the possibility of dynamical scalarization occurring at frequencies f DS 50 Hz at 90% CL. Initial detectability studies -Refs. [37,83] discussed above -suggest that this early scalarization should be observable with aLIGO (although these conclusions should be confirmed with future work in light of the limitations of these works; see above). Thus, failure to detect dynamical scalarization in future GW observations could provide tighter constraints on the parameters (α 0 , β 0 ) in DEF theory than pulsar timing. However, as can be seen from Table IV, the prospects of producing such complementary constraints depend critically on the observed NS masses and the EOS of NSs.

V. CONCLUSIONS
In this paper, we have studied the scalarization phenomena [14,24] in the massless mono-scalar-tensor theory of gravity of DEF with pulsar timing and laser-interferometer GW detectors on the Earth. We now summarize the key conclusions of our analysis.
1. The spontaneous scalarization phenomenon [14] occurs at different NS mass ranges for different EOSs [20]. Therefore, in a well-timed relativistic binary-pulsar system with a specific NS mass, the scalar-tensor gravity might be stringently constrained for some EOSs whose spontaneous-scalarization phenomenon occurs near that specific NS mass. However, in general, strong scalarization could still take place if NSs are described by an EOS whose scalarization occurs at a mass different from the one observed.
2. Combining two well-timed binary-pulsar systems with quite different NS masses, one could in principle constrain the scalar-tensor theory with whatever EOS Nature provides us. Using MCMC simulations, we showed in Sec. III that, combining five binary pulsars [32,35,[45][46][47] that best constrain gravitational dipolar radiation, we can already bound the scalarization parameter, β 0 , to be −4.28 at 68% CL and −4.38 at 90% CL, for any of the eleven EOSs that we have considered.  Table III).

5.
We performed a Fisher-matrix study of BNS inspiral signals using aLIGO, CE, and ET. We found that for BNSs at a luminosity distance D L = 200 Mpc, where we expect to observe those sources, aLIGO can still improve the limits from binary pulsars for a couple of EOSs with BNSs of suitable masses. CE (whose bandwith starts at 5 Hz) can improve the current limits for all EOSs, while ET (whose bandwith starts at 1 Hz) will provide us with even more significant improvements over current constraints for all EOSs. This is mainly due its better low-frequency sensitivity. Our conclusions for aLIGO differ from the one obtained in Refs. [16,33,34], where the authors concluded that pulsar timing would do better than aLIGO in constraining scalar-tensor theories. The main reason of this difference comes from a better understanding and larger span of the NS masses and EOSs during the past two decades [31,32], and the different PSD for aLIGO used in Refs. [16,62]. If we restricted the analysis to the same NS masses and the same EOS used in Ref. [16], we would recover the same conclusions as in Ref. [16] (see Fig. 7).
6. We investigated dynamical scalarization in equal-mass and unequal-mass BNS systems. With the criterion that the dynamical scalarization transition frequency must fall below ∼ 50 Hz [37,83] to be detectable, we found aLIGO could be able to observe this phenomenon given the constraints obtained from binary-pulsar timing, even away from the scalarization windows. We found that the prospects for observing dynamical scalarization with GW detectors depends critically on the NS EOSfor example, dynamical scalarization of NSs with the MS0 EOS could not be detected with aLIGO. Producing new constraints on scalar-tensor theories from GW searches for dynamical scalarization requires waveform models that can faithfully reproduce this nonperturbative phenomenon; ultimately, these conclusions should be revisited once such models are developed.
Our comparisons between binary pulsars and GWs made use of the current limits of the former and the expected limits of the latter. It shows that advanced and next-generation ground-based GW detectors have potential to further improve the current limits set by pulsar timing. Nevertheless, the binary-pulsar limits will also improve over time, especially if suitable systems filling the scalarization windows are discovered in future pulsar surveys. Better mass measurements of currently known pulsars will also help in narrowing down the constraints, especially with PSRs J1012+5307 [45] and J1913+1102 [77], whose observational uncertainties in masses are still large, and they might have the right masses to close the windows below 2 M . To reach this goal, the next generation of radio telescopes, such as FAST and SKA will play a particularly important role [22,84]. ultra-relativistic BNS depicted in the top panel of Fig. 11, this transition is completed after an evolution of only 0.2 Hz.