1 Introduction

The problem of binary motion in a gravitational system is one of the oldest in astronomy. Newton’s explanation of Kepler’s laws was a milestone of science. Relativistic corrections to Keplerian motion under Newton’s laws have been studied since the time of Einstein [1]. These corrections can be calculated within the post-Newtonian (PN) formalism as an expansion in the ratio v/c of the typical relative velocity v between the two components of the binary to the speed of light c. Low order corrections can be observed in the motion of binary pulsar systems where typical values of v/c are \(\mathcal {O}(10^{-4})\) [2]. The recent observations of coalescing compact binary systems (see e.g. [3,4,5,6,7,8,9]) have much higher values of \(v/c \sim \mathcal {O}(10^{-1})\). These observations therefore probe regimes where the gravitational field is strong and also dynamical.

The relation between the binary motion and the observed gravitational waves was first expounded by Peters and Mathews [10]. This work obtained the gravitational wave phase by calculating the Einstein quadrupole formula for underlying Keplerian motion. Whilst this analysis gives a very approximate fit for LIGO data [11] it has been known for a long time that more accurate fits require a large number of relativistic corrections [12]. Several different approaches have been able to model the expected signal for the inspiral and coalescence of compact objects with sufficient accuracy within standard general relativity. These include the Effective-one-body models (see e.g. [13,14,15,16,17,18,19,20,21]), the Phenomenological models (see e.g. [22,23,24,25,26,27,28]) and the surrogate models [29,30,31]. Some corrections have also been calculated in modified gravity theories [32, 33]. In this paper we restrict our attention to the binary neutron star merger GW170817 event [34]. For this event, the contribution of the merger to the signal-to-noise-ratio is not significant, whence the inspiral regime and standard PN methods dominate [35,36,37,38,39].

It has now become routine to use the observed gravitational wave observations to constrain relativistic corrections to binary motion [5, 40,41,42,43,44,45]. The TIGER framework [46] looks for a deviation from the calculated GR value by introducing a number of phenomenological deviation parameters that take the value zero in GR. In this approach the deviation parameters are assumed to be independent of the physical parameters of the binary, such as mass and spin.

These phenomenological non-GR parameters implicitly represent a certain class of GR alternatives. In different GR alternatives, these parameters may have different physical meanings. For instance, in phenomenological modifications of GR based on the massive graviton assumption, such that the effective Newtonian potential is of Yukawa type with a non-standard graviton dispersion relation, the non-GR parameter of 1PN order is proportional to the graviton’s mass squared [47]. In other alternative theories of gravity, radiative multipole moments differ from those of GR and their deviation from GR values contributes to every PN order phase term [48, 49]. In general, alternatives to GR may naturally depend on an additional set of parameters such that the non-GR parameters are in general functionally related. This motivates testing GR via the simultaneous variation of such non-GR parameters introduced as additional terms in the GW phase. In GR tests two types of variations were considered: each of these parameters was varied individually whilst keeping the others fixed, e.g. [5, 44, 45], or all the parameters were varied simultaneously [43, 50]. The results showed that simultaneous variation produced wide and non-informative posteriors on the value of each individual deviation parameter, closely consistent with the prior. Reported bounds on non-GR parameters obtained from the double pulsar J0737-3039 [51] were computed individually, one at a time [43]. This is because in binary pulsars the orbital period changes at essentially a constant rate and when the PN order is increased, the corresponding bounds quickly become rather loose and significantly less informative than those from GW’s binaries coalescence.

Whilst varying each one individually captures a generic deviation from GR if it is strong enough, a generic modified gravity theory is expected to differ from GR for several of the non-GR parameters at different PN orders. This also means that one cannot directly compare the results obtained on each parameter individually with a theoretical calculation in an alternate theory of gravity (note, however, that in some specific cases of alternative theories of gravity, one can make such a comparison by using perturbative approach, see, e.g., [52,53,54,55]). This aspect was also emphasised in the work [50], where the authors used the Fisher information matrix to estimate constraints on non-GR parameters and found that in the best case of multiband observations the deviations at orders below 3PN are less than \(1\%\). They reported the deviations of \(\sim 100\%\) and higher for simultaneous measurement of four and more deformation parameters in the case of ground-based observatories or LISA observations alone. In addition, the non-GR parameters are correlated with the PN coefficients and between themselves, which can be seen for example in Figure 7 of [43]. To remedy these limitations, the method of principal component analysis was proposed in [56]. A similar method, known as the singular value decomposition approach, was proposed in [57]. These methods are applied to the Fisher information matrix that estimates errors in the parameter values. However, Fisher-matrix-based estimates are not reliable at low SNRs (see, e.g. [12] and also [58]). This limits the measurability of the higher-order PN terms, as well as their non-GR contributions.

In this work, we allow the waveform phase to vary in additive non-GR parameters sector. Using our numerical results, we construct covariance matrix of these parameters whose diagonal elements represent their variance and off-diagonal terms represent their covariance that defines their pairwise correlation. Taking only the diagonal elements and ignoring the off-diagonal ones, as it was done in many works, results in ignoring the information about the parameters correlation. In a simple example of a two-dimensional Gaussian distribution, the principal directions of the concentration ellipse define new orthogonal variables whose variances are described optimally. Accordingly, such a transformation provides more accurate measurement of the leading order (uncorrelated) non-GR parameters (see Eqs. (7) and (8) below). A similar approach based on diagonalization of the non-GR parameters covariance matrix is given in the work [59], which appeared in the arXiv after submission of our paper.

2 Posterior distributions of the non-GR parameters

When the rate of change of frequency is small relative to the squared frequency, the expansion of a GW signal h(t) in the Fourier domain can be obtained from an expansion in the time domain using the Stationary Phase Approximation [12]. For the dominant frequency, the Fourier domain waveform can be written as

$$\begin{aligned} \tilde{h}(f) = A(f)e^{i\Psi (f)} \end{aligned}$$
(1)

where \(\tilde{h}(f)\) is the Fourier transform of h(t) and

$$\begin{aligned} \Psi (f)= & {} \sum ^{8}_{n=0}\left( \frac{f}{f_{{\tiny ref}}}\right) ^{(n-5)/3}\left[ \Psi _{n}+\Psi ^{\ln }_{n}\ln \left( \frac{f}{f_{{\tiny ref}}}\right) \right] . \end{aligned}$$
(2)

Here \(\Psi _{1}\equiv 0\) and the only non-zero logarithmic terms are \(\Psi ^{\ln }_{5}\) and \(\Psi ^{\ln }_{6}\). The expansion order n corresponds to (n/2)PN term. In this expansion the 2.5PN term becomes indistinguishable from the binary coalescence phase term \((-\phi _{c}-\pi /4)\) and the 4PN term is indistinguishable from the binary coalescence time term \(2\pi ft_{c}\). The gravitational wave amplitude A can also be expanded as a post-Newtonian series, but here we will keep only the terms at Newtonian order and hence \(A(f) \sim f^{-7/6}\). Factors of a reference frequency \(f_{{\tiny ref}}\) are included as in [56] to render the expansion coefficients explicitly dimensionless. In what follows, we shall take the “natural”reference frequency \(f_{{\tiny ref}}=1/(\pi M)\). This choice matches the expansion (2) of the TaylorF2 approximant form (see, e.g., [60]), which is also implemented in the LALSuite software package [61].Footnote 1

The physical parameters that enter our waveform are the binary masses \(m_1\) and \(m_2\), the orbit-aligned dimensionless spin components of the stars \(\chi _{z1}\) and \(\chi _{z2}\) (we consider high-spin priors), the dimensionless tidal deformability parameters \(\Lambda _1\) and \(\Lambda _2\), which already appear in the 5PN term (here we present the posterior of the dimensionless combined tidal deformability \(\tilde{\Lambda }\)), and the following set of non-GR parametersFootnote 2:

$$\begin{aligned} \delta \Psi _{0}, \delta \Psi _{1}, \delta \Psi _{2}, \delta \Psi _{3}, \delta \Psi _{4}\,, \end{aligned}$$
(3)

that represent the absolute deviations in their corresponding PN coefficients,

$$\begin{aligned} \Psi _{n}\longrightarrow \Psi _{n}+\delta \Psi _{n}\,. \end{aligned}$$
(4)

The work [62] criticises parametric tests of GR, which are not based on explicit non-GR waveform models, appealing to the fact that the fractional precision of non-GR parameters determination is not parametrisation-invariant. For a particular non-GR theory in mind, these parameters are functions of the theory parameters, whose values, given the corresponding waveform model, can be estimated. Although the physical interpretation of our parameters is not in direct one-to-one relationship with the PN expansion coefficients of TaylorF2, our parameters still set observational constraints, even on theories in which more than one PN coefficient differs from its GR value. This is the case in most modified theories of GR.

Our goal is to measure values of such deviations regardless of any specific non-GR theory and these values are “absolute”in the sense that the functions values do not depend on their variables choice. To obtain the posterior density distributions of these parameters we exploit the PyCBC inference package [63] where we had fixed the reported value [34] of the luminosity distance \(40^{+8}_{-14}\) Mpc. In our runs, we used uniform prior distributions and varied all the parameters simultaneously. The supplementary file [64] contains posteriors of the binary physical parameters, assuming that GR is correct, i.e. the values of (3) are set to zero, Fig. 1, one- and two-dimensional posteriors of all the parameters, Fig. 2, and posteriors of the non-GR parameters, Fig. 3. Our results show that the dimensionless combined tidal deformability \(\tilde{\Lambda }\) is correlated with other parameters, and in particular with \(\delta \Psi _{4}\) non-GR parameter. The non-GR parameters are strongly correlated, in particular \(\delta \Psi _{n}\) and \(\delta \Psi _{n+1}\), for \(n=1,2,3\). The off-diagonal (correlation) terms corresponding to these parameters in the \((5\times 5)\) covariance matrix \(\Sigma _{nm}\) (see the supplementary file [64]) make these parameters less accurately measured.Footnote 3 Moreover, because of the correlation, we cannot infer directly the accuracy of these parameters from their 1D posteriors. Thus, we are to define a related set of uncorrelated non-GR parameters.Footnote 4 In order to statistically isolate non-GR parameters from each other, we diagonalise their covariance matrix derived from the sample data arrays from the PyCBC inference package [63]. This approach does not rely on the high SNR values that are required for the Fisher information matrix estimate. We solve the eigenvalue-eigenvectors problem and derive the transformation matrix \(T_{nm}\), which diagonalises the covariance matrix (see the supplementary file [64]). Its transposed form \(T_{mn}\) relates the correlated and new uncorrelated non-GR parameters \(\delta \hat{\Psi }_{n}\):

$$\begin{aligned} \delta \hat{\Psi }_{m}=T_{mn}\delta {\Psi }_{n}\,. \end{aligned}$$
(5)

These matrix elements \(T_{mn}\) will obviously depend on the detector noise spectrum and on the source parameters as well. The transformation matrix \(T_{mn}\) reads

$$\begin{aligned} \left[ \begin{array}{ccccc} 0.303 &{} -0.947 &{} -0.107 &{} -0.400\times 10^{-2} &{} 0.469\times 10^{-3} \\ 0.953 &{} 0.302 &{} 0.298\times 10^{-1} &{} 0.557\times 10^{-3} &{} -0.304\times 10^{-3} \\ 0.424\times 10^{-2} &{} -0.111 &{} 0.980 &{} 0.163 &{} 0.117\times 10^{-1} \\ -0.121\times 10^{-4} &{} -0.143\times 10^{-1} &{} 0.161 &{} -0.968 &{} -0.190 \\ 0.973\times 10^{-4} &{} -0.911\times 10^{-3} &{} 0.197\times 10^{-1} &{} -0.190 &{} 0.982 \\ \end{array}\right] \end{aligned}$$
(6)

The matrix acts on the column vector \([\delta {\Psi }_{0}\,\,\delta {\Psi }_{1}\,\,\delta {\Psi }_{2}\,\,\delta {\Psi }_{3}\,\,\delta {\Psi }_{4}]^{T}\) and produces the column vector \([\delta \hat{\Psi }_{0}\,\,\delta \hat{\Psi }_{1}\,\,\delta \hat{\Psi }_{2}\,\,\delta \hat{\Psi }_{3}\,\,\delta \hat{\Psi }_{4}]^{T}\), where the superscript “T”stands for the transpose. The eigenvalue-eigenvectors problem computations were performed using the Householder reduction followed by the QL algorithm. Details of these methods can be found in the book [65].

To illustrate the result of this transformation we present one- and two-dimensional posteriors of the non-GR parameters in Fig. 1. Posteriors of the binary parameters together with the uncorrelated non-GR parameters are presented in Fig. 4 in the supplementary file [64].

Fig. 1
figure 1

Posterior density distributions of the uncorrelated non-GR parameters. The central plots are 2D marginal posteriors, where the black contours show the 50% and 90% credible regions. The upper and the right plots are the 1D marginal posteriors, where the median and 90% credible intervals are indicated by the dashed lines. Red contours are defined by vanishing non-GR parameters indicated by intersection of the green dashed lines. They define confidence regions of the 2D joint distributions. For the \(\delta \hat{\Psi }_{1}\) and \(\delta \hat{\Psi }_{2}\) joint distribution, the GR confidence region is degenerated to the single point (0, 0)

Full size image

One can see that the new non-GR parameters are indeed almost uncorrelated. Moreover, the measured values of the new non-GR parameters are overall more constrained, as compared to the correlated ones. We can present our results in terms of the relative  \(\delta {\psi }_{n}=\delta {\Psi }_{n}/\Psi _{n}\) values of the non-GR parameters considered in other works, e.g. [44, 59]. Accordingly, we compute the values of the relative correlated \(\delta {\psi }_{n}\) and uncorrelated \(\delta \hat{\psi }_{n}\) non-GR parameters (except for 0.5PN ones, which are taken to be identical to the absolute ones, for in GR 0.5PN term vanishes) as follows. First, we compute the relative correlated non-GR parameters by taking the values of the masses and spins (see Fig. 2 in [64]) and computing the first four PN terms (see, e.g. [60, 61]) and dividing by them the corresponding non-GR correlated parameters (see Fig. 2 in [64]),

$$\begin{aligned}{} & {} \vert \delta {\psi }_{0}\vert =1.05\%, \vert \delta {\psi }_{1}\vert =2.40\%, \vert \delta {\psi }_{2}\vert =57.4\%, \nonumber \\{} & {} \quad \vert \delta {\psi }_{3}\vert =40.4\%, \vert \delta {\psi }_{4}\vert =91.5\%. \end{aligned}$$
(7)

Second, we compute the relative uncorrelated non-GR parameters by taking the sample data arrays of the masses, spins, and non-GR parameters from our PyCBC inference and constructing from them the arrays of the relative non-GR parameters. Then we proceed along the same steps as above and construct their covariance matrix, diagonalise it by means of the transformation matrix, and apply its transpose to the computed relative uncorrelated non-GR parameters. As a result, we derive

$$\begin{aligned}{} & {} \vert \delta \hat{\psi }_{0}\vert =0.605\%, \vert \delta \hat{\psi }_{1}\vert =1.92\%, \vert \delta \hat{\psi }_{2}\vert =5.78\%,\nonumber \\{} & {} \quad \vert \delta \hat{\psi }_{3}\vert =48.6\%, \vert \delta \hat{\psi }_{4}\vert =104\%. \end{aligned}$$
(8)

These results show that despite of simultaneous variations, our constraints are more stringent for the leading order non-GR parameters.

We note that the presence of the non-GR parameters in the waveform model does affect the values of the physical parameters corresponding to GR waveform model (for comparison see Fig. 1 in the supplementary file.) However, the 90% credible intervals of the “new”values of physical parameters and the “old”ones largely overlap.

3 Analysis of the results: Measures of deviations from GR

From the posteriors of the non-GR parameters one can immediately infer that the GR prediction, i.e. \(\delta \Psi _{n}=0\), \(n=0,...,4\) falls well within their 90% credible intervals. However, one may prefer to have additional quantitative estimates of GR validity. In this section, we measure deviations from GR based on our derived results.

3.1 Integral non-GR phase

According to our approach (2), (4), the waveform phase is an additive combination of the GR and non-GR components,

$$\begin{aligned} \Psi (f)=\Psi _{{\tiny GR}}(f)+\Psi _{{\tiny non-GR}}(f)\,, \end{aligned}$$
(9)

where

$$\begin{aligned} \Psi _{{\tiny non-GR}}(f) = \sum ^{4}_{n=0}\left( \frac{f}{f_{{\tiny ref}}}\right) ^{(n-5)/3}\delta \Psi _{n}\, \end{aligned}$$
(10)

is the non-GR part of the waveform phase. It indicates the phase difference due to potential non-GR contribution and thus represents another, “effective”non-GR term. To measure this term it is reasonable to introduce is cumulative value given by the integral non-GR phase

$$\begin{aligned} \Psi _{{\tiny int-non-GR}}=\int _{f_{{\tiny min}}/f_{{\tiny ref}}}^{f_{{\tiny max}}/f_{{\tiny ref}}}\Psi _{{\tiny non-GR}}(f)d\left( f/f_{{\tiny ref}}\right) \,. \end{aligned}$$
(11)

The integral non-GR phase can be computed by using the estimated values of the non-GR parameters and the expression (10) above. Note that the non-GR part of the waveform phase is a vector in the non-GR parameter subspace spanned by the basis \(\{\delta \Psi _{n},\,n=0,...,4\}\) which is invariant under the transformation (5) of the uncorrelated non-GR parameters. To see that one should also diagonalise the corresponding frequency components \(\{(f/f_{{\tiny ref}})^{(n-5)/3},\, n=0,...,4\}\) by the transposed of the matrix (6). Accordingly, the quantity \(\Psi _{{\tiny int-non-GR}}\) is invariant as well.

We may also estimate uncertainty \(\sigma \) of the non-GR phase by using the propagation of errors approach. The variance of the function (10) is

$$\begin{aligned} \sigma ^2(f) = \sum ^{4}_{n,m=0}\left( \frac{f}{f_{{\tiny ref}}}\right) ^{(m-5)/3}\Sigma _{nm}\left( \frac{f}{f_{{\tiny ref}}}\right) ^{(n-5)/3}\,. \end{aligned}$$
(12)

Integrating this expression in the frequency range, we derive the variance of the non-GR phase

$$\begin{aligned} \sigma ^2_{{\tiny int}}=\int _{f_{{\tiny min}}/f_{{\tiny ref}}}^{f_{{\tiny max}}/f_{{\tiny ref}}}\sigma ^2(f)d\left( f/f_{{\tiny ref}}\right) \,. \end{aligned}$$
(13)

For our data we find

$$\begin{aligned} \Psi _{{\tiny int-non-GR}} = 0.0447\pm 25.3000\,. \end{aligned}$$
(14)

Thus, the GR prediction \(\Psi _{{\tiny int-non-GR}}=0\) falls well within (14).

3.2 Confidence regions and deviation from GR percentile

Here we present another ways to measure GR validity. First, we construct confidence regions whose boundary corresponds to the predicted by GR zero values of the non-GR parameters. This approach is based on the kernel density estimation (kde). For a 2D kde of each pair of the non-GR parameters we used the Python function scipy.stats.gaussian_kde from the SciPy library, that is implemented by PyCBC by the default. The confidence regions are shown in Fig. 1 and also in Fig. 3 in the supplementary file.

Second, we calculate the deviation from GR percentile \(p^{{\tiny Dev-GR}}_{n}\) in the following way. The PyCBC inference hdf file contains the posterior for original non-GR PN parameters. We extract the effective samples from the inference hdf file by using pycbc_inference_extract_samples command to get a new hdf file. From this new hdf file, we read the posterior for each original non-GR PN parameters \(\delta \Psi _n\)’s as an array using the standard reading hdf file command. In order to get array for each new non-GR parameter \(\delta \hat{\Psi }_{n}\), we use the expression (5). This gives us array for each new non-GR parameter \(\delta \hat{\Psi }_{n}\). To get the probability density \(p(\delta \hat{\Psi }_{n})\) we use the Python function scipy.stats.gaussian_kde.

The original non-GR PN parameters are correlated, but the new non-GR parameters are almost uncorrelated. Exploiting this fact, we construct the probability density P in 5 dimensions for the new non-GR parameters by multiplying the individual probability densities \(p(\delta \hat{\Psi }_{n})\)’s of the new non-GR parameters,

$$\begin{aligned} P(\delta \hat{\Psi }_{n};n=0,...,4) = \prod _{n=0}^{4} p(\delta \hat{\Psi }_{n})\,. \end{aligned}$$
(15)

Next, we define the level hypersurface \(\Sigma :\,P=\text{ const }\) corresponding to GR prediction, i.e. to zero values of the new non-GR parameters,

$$\begin{aligned} \Sigma :\,P(\delta \hat{\Psi }_{n};n=0,...,4) = P(\delta \hat{\Psi }_n=0;n=0,...,4)\,. \end{aligned}$$
(16)

Finally, to get the deviation from GR percentile, we integrate P over the domain \(D(\Sigma )\) enclosed by this hypersurface,

$$\begin{aligned} p^{\text{ Dev-GR }}_{n} = \int _{D(\Sigma )} P \,d^5(\delta \hat{\Psi }_{n})\,. \end{aligned}$$
(17)

To get the value

$$\begin{aligned} p^{{\tiny Dev-GR}}_{n}=25.85\% \end{aligned}$$
(18)

of the deviation from GR percentile, we performed the integration numerically by using the Vegas package, which exploits the adaptive multidimensional Monte Carlo integration.

4 Conclusion

We have analysed deviations from GR with the example of the binary neutron star signal GW170817. A set of 5 non-GR parameters up to and including the 2PN order was introduced into the GW waveform in the TaylorF2 phase approximant. These parameters represent absolute deviations of the post-Newtonian coefficients from their GR values. By using the PyCBC package we computed their posterior density distributions for the simultaneous variation of all these non-GR parameters using uniform prior distributions. Next, diagonalizing their \((5\times 5)\) covariance matrix, we derived linear combinations of these parameters that represent a set of uncorrelated non-GR parameters. This approach allows the first study of each of the new non-GR parameters independently. Each of these uncorrelated parameters corresponds to a covariance matrix principal direction in the 5D subspace of the parameter space.

Our results provide more stringent constraints on GR than those presented in [5, 44, 45]. In comparison to the principal component analysis [56] or the singular value decomposition approach [57] applied to the Fisher matrix, our computations provide a more efficient technique that achieves more stringent constraints on parameters estimation. The reason is that diagonalization of the Fisher matrix followed by the parameters estimation run gives less stringent posteriors as compared to diagonalization of the covariance matrix constructed from the computed data. This is due to requirement of the high SNR value for a good Fisher information matrix estimate. A detailed analysis and comparison of both methods on other GW events may be more informative and requires additional investigation, which is beyond the scope of our paper.

Although the approach presented here, based on orthogonalisation of the covariance matrix, looks quite promising, one can also analyse a “hybrid method”, that combines both the principal component analysis of the Fisher matrix and the orthogonalisation of the covariance matrix constructed from the computed data. A more detailed study of this possibility is left for future work.

5 Supplementary information

Our manuscript has accompanying supplementary file [64].