Testing the post-Newtonian expansion with GW170817

Observations of gravitational waves from compact binary mergers have enabled unique tests of general relativity in the dynamical and non-linear regimes. One of the most important such tests is constraints on the post-Newtonian (PN) corrections to the phase of the gravitational wave signal. The values of the PN coefficients can be calculated within standard general relativity, and these values are different in many alternate theories of gravity. It is clearly of great interest to constrain the deviations based on gravitational wave observations. In the majority of such tests which have been carried out, and which yield by far the most stringent constraints, it is common to vary these PN coefficients individually. While this might in principle be useful for detecting certain deviations from standard general relativity, it is a serious limitation. For example, we would expect alternate theories of gravity to generically have additional parameters. The corrections to the PN coefficients would be expected to depend on these additional non-GR parameters, whence, we expect that the various PN coefficients to be highly correlated. We present an alternate analysis here using data from the binary neutron star coalescence GW170817. Our analysis uses an appropriate linear combination of non-GR parameters that represent absolute deviations from the corresponding post-Newtonian inspiral coefficients in the TaylorF2 approximant phase. These combinations represent uncorrelated non-GR parameters which correspond to principal directions of their covariance matrix in the parameter subspace. Our results illustrate good agreement with GR. In particular, the integral non-GR phase is Ψint-non-GR=0.0447±25.3000\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Psi _{{\tiny int-non-GR}} = 0.0447\pm 25.3000$$\end{document} and the deviation from GR percentile is pnDev-GR=25.85%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{{\tiny Dev-GR}}_{n}=25.85\%$$\end{document}.


I. 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 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 ∼ 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.
It has now become routine to use the observed gravitational wave observations to constrain relativistic corrections to binary motion [5,[35][36][37][38][39][40].The TIGER framework [41] 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 [42].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 [43,44].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 the tests of [5,39,40] two types of variations were considered: each of these parameters was varied individually whilst keeping the others fixed, or all the parameters were varied simultaneously.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 [38].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.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 [38].To remedy these limitations, the method of principal component analysis was proposed in [45].A similar method, known as the singular value decomposition approach, was proposed in [46].These methods are applied to the Fisher information matrix that estimates parameter values.However, Fishermatrix-based estimates are not reliable at low SNRs (see, e.g.[12] and also [47]).This limits the measurability of the higher-order PN terms, as well as their non-GR contributions.
In this work, in order to statistically isolate non-GR parameters we diagonalise their covariance matrix derived from the sample data arrays from the PyCBC inference package [48].This approach does not rely on the high SNR values that are required for the Fisher information matrix estimate.Using the eigenvectors of the covariance matrix we construct new uncorrelated non-GR parameters (that are linear combinations of the original ones) corresponding to the principal directions of the covariance matrix.In this work, we consider simultaneous variation of non-GR parameters introduced additively to the PN phase expansion parameters from 0PN to 2PN order, although it could be extended to other orders.We demonstrate our method on the signal GW170817 that was identified by the LIGO-Virgo detector network as the first gravitational wave observation of a binary neutron star inspiral [34].
Although, the combined signal-to-noise ratio of 32.4 for GW170817 was the largest of any gravitational wave event observed to date, our primary interest in this system is because of its low mass and long inspiral phase in the detector band.The chirp mass was found to be M = 1.188 +0.004 −0.002 M and with the restriction of the component spins to the range inferred in binary neutron stars, the component masses were found to be in the range 1.17 − 1.60 M , with the total system mass of 2.74 +0.04 −0.01 M .The signal source sky position was localised to an area of 28 deg 2 with 90% probability and had a luminosity distance of 40 +8 −14 Mpc.Our paper is organised as follows.In Sec.II we introduce the waveform model and its TaylorF2 approximation up to the 4PN order.We also introduce a set of non-GR additive parameters to the corresponding PN coefficients.All these parameters are given uniform prior distributions and they are varied simultaneously.We derive the corresponding posteriors of all the waveform parameters.Next, diagonalizing of the non-GR parameters covariance matrix, we construct new uncorrelated non-GR parameters and repeat the procedure.In Sec.III we analyse posteriors of the non-GR parameters and define two different methods to measure deviations from GR.We summarise our results in the Conclusion.

II. 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 where h(f ) is the Fourier transform of h(t) and Here Ψ 1 ≡ 0 and the only non-zero logarithmic terms are Ψ ln 5 and Ψ 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 (−φ c − π/4) and the 4PN term is indistinguishable from the binary coalescence time 2πf t 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 ) ∼ f −7/6 .Factors of a reference frequency f 0 are included as in [45] to render the expansion coefficients explicitly dimensionless.In what follows, we shall take the "natural" reference frequency f ref = 1/(πM ).This choice matches the expansion (2) of the TaylorF2 approximant form (see, e.g., [49]), which is also implemented in the LALSuite software package [50].
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 χ z1 and χ z2 (we consider high-spin priors), the dimensionless deformability parameters Λ 1 and Λ 2 , which already appear in the 5PN term (here we present the posterior of the dimensionless combined tidal deformability Λ), and the following set of non-GR parameters 1 : that represent the absolute deviations in their corresponding PN coefficients, To obtain the posterior density distributions of these parameters we exploit the PyCBC inference package [48].
In particular, we use the marginalized phase model and the emcee pt sampler with 2000 walkers, 12000 effective samples, with a maximum of 1000 samples per chain, and four temperatures.We used (uniform) prior distributions and varied all the parameters simultaneously.Oneand two-dimensional posteriors of all the parameters are shown in Fig. 1 and posteriors of the non-GR parameters are shown separately in Fig. 2. One can notice that the dimensionless combined tidal deformability Λ is correlated with other parameters, and in particular with δΨ 4 non-GR parameter.
One can infer from Fig. 2 that the non-GR parameters are strongly correlated, in particular δΨ n and δΨ n+1 , for n = 1, 2, 3.The off-diagonal (correlation) terms corresponding to these parameters in the (5 × 5) covariance matrix Σ nm make these parameters less accurately measured.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. 2 To construct such a set we need to diagonalize the covariance matrix Σ nm .We solve the eigenvalue-eigenvectors problem and derive the transformation matrix T nm , which diagonalises the covariance matrix.Its transposed form T mn relates the correlated and new uncorrelated non-GR parameters δ Ψn : These matrix elements T mn will obviously depend on the detector noise spectrum and on the source parameters as well.Explicitly, the matrix elements which define the transformation are: 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 [51].
To illustrate the result of this transformation we present one-and two-dimensional posteriors of the binary parameters together with the uncorrelated non-GR parameters in Fig. 3 and separately the posteriors of the non-GR parameters in Fig. 4. 1 Here we restricted the non-GR parameters up to the 2PN order.
Measurements of the higher orders are less accurate and their posteriors do not converge well.We include point-particle PN terms from GR up to 4PN and tidal deformation terms up to 6PN. 2 Ideally, we could construct a complete set of all 10 uncorrelated parameters.However, such parameters would represent a mixture of the physical binary parameters and the non-GR parameters and this would alter their original physical meaning.Thus, to avoid this situation, we keep the non-GR sector in the parameter space analytically (but not statistically) isolated.
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, i.e., more accurate.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.However, the 90% credible intervals of the "new" values of physical parameters and the "old" ones largely overlap.
From the posteriors of the non-GR parameters one can immediately infer that the GR prediction, i.e. δΨ 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.

A. Integral non-GR phase
According to our approach (2), (4), the waveform phase is an additive combination of the GR and non-GR components, where is the non-GR part of the waveform phase.In order to measure the non-GR component value we introduce the integral non-GR phase The integral non-GR phase can be computed by using the estimated values of the non-GR parameters [see Fig. 2] and the expression (8) above.
We may also estimate uncertainty σ of the non-GR phase by using the propagation of errors approach.The variance of the function ( 8) is Integrating this expression in the frequency range, we derive the variance of the non-GR phase For our data we find Thus, the GR prediction Ψ non-GR = 0 falls well within (12).

B. 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.gaussiankde from the SciPy library, that is implemented by PyCBC by the default.The confidence regions are shown in Fig. 2 and Fig. 4.
Second, we calculate the deviation from GR percentile p 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 δΨ n 's as an array using the standard reading hdf file command.In order to get array for each new non-GR parameter δ Ψn , we use the expression (6).This gives us array for each new non-GR parameter δ Ψn .To get the probability density p(δ Ψn ) we use the Python function scipy.stats.gaussiankde.
To get the value of the deviation from GR percentile, we performed the integration numerically by using the vegas package, which exploits the adaptive multidimensional Monte Carlo integration.

IV. 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 × 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 separately.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,39,40].In comparison to the principal component analysis [45] or the singular value decomposition approach [46] applied to the Fisher matrix, our computations provide a more efficient technique that achieves more stringent constraints on parameters estimation.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.For the δ Ψ1 and δ Ψ2 joint distribution, the GR confidence region is degenerated to the single point (0, 0).

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

FIG. 2 .
FIG. 2. Posterior density distributions of the non-GR parameters from Fig. 1.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.