Towards a framework for testing general relativity with extreme-mass-ratio-inspiral observations
Alvin J. K. Chua Sonke Hee Will J. Handley Edward Higson Christopher J. Moore Jonathan R. Gair Michael P. Hobson Anthony N. Lasenby



firstpage–lastpage
Exact description of the boundary theory of the Kitaev Toric Code with open boundary conditions
    Stefan Kehrein
    December 30, 2023
===============================================================================================



 
Extreme-mass-ratio-inspiral observations from future space-based gravitational-wave detectors such as LISA will enable strong-field tests of general relativity with unprecedented precision, but at prohibitive computational cost if existing statistical techniques are used. In one such test that is currently employed for LIGO black-hole binary mergers, generic deviations from relativity are represented by N deformation parameters in a generalised waveform model; the Bayesian evidence for each of its 2^N combinatorial submodels is then combined into a posterior odds ratio for modified gravity over relativity in a null-hypothesis test. We adapt and apply this test to a generalised model for extreme-mass-ratio inspirals constructed on deformed black-hole spacetimes, and focus our investigation on how computational efficiency can be increased through an evidence-free method of model selection. This method is akin to the algorithm known as product-space Markov chain Monte Carlo, but uses nested sampling and improved error estimates from a rethreading technique. We perform benchmarking and robustness checks for the method, and find order-of-magnitude computational gains over regular nested sampling in the case of synthetic data generated from the null model.





gravitational waves – methods: data analysis – methods: statistical





 


§ INTRODUCTION



The spate of gravitational-wave (GW) sources found by Advanced LIGO in its first two observing runs <cit.> have opened up a new branch of multimessenger astronomy—one that extends its reach beyond electromagnetic radiation for the first time, and into the gravitational sector. Traditional electromagnetic telescopes will work together with a ground-based GW detector network <cit.>, pulsar timing arrays <cit.>, and future space-based GW detectors <cit.> to discover and study a broad spectrum of sources that are astrophysical or cosmological in origin.

Astronomy with GW observations will also improve our understanding of gravitation and fundamental physics, by granting access to unprecedented tests of general relativity (GR) and its alternatives in the dynamical strong-field regime <cit.>. Several such tests have been performed on data from Advanced LIGO's first observing run, with no evidence for deviation from GR to date <cit.>; these include various GR consistency checks for measured signals, as well as the placing of constraints on non-tensorial GW polarisations in generic metric theories, and on the graviton Compton wavelength in massive-gravity theories.

One particular test focuses on the late-stage phase evolution of GW signals from black-hole binary mergers, which encode several strong-field effects that are not observable for binary pulsars. This approach is a particular implementation of the parametrised post-Einsteinian framework <cit.>, in which model-independent deviations from GR are parametrised as deformations to the post-Newtonian (PN) and phenomenological phase parameters in a GR-based inspiral–merger–ringdown waveform model <cit.>. Any phase deformations for a given signal are then constrained through comparison with the generalised model, and a Bayesian model selection framework <cit.> is used to perform a null-hypothesis test of GR.

Merging binary systems have provided the only GW signals detected so far, and are expected to be ubiquitous across the bandwidths of current and future interferometric detectors. They will be an important type of source for ESA's space-based detector LISA <cit.>, in the form of massive-black-hole (MBH) binary mergers and the extreme-mass-ratio inspirals (EMRIs) of stellar-origin compact objects into MBHs within galactic nuclei <cit.>. EMRIs in particular have the potential to facilitate stringent tests of GR <cit.>; this is because the phase of the GW signal can be precisely tracked over the large number of observable orbits spent by the compact object in the strong field of the central black hole.

Model-independent parametric tests of GR with EMRIs may be performed using an EMRI analogue of the generalised waveform model for comparable-mass black-hole binary mergers. One such model <cit.> introduces deformations at the level of the background metric around the central MBH, by constructing analytic EMRI waveforms <cit.> on a bumpy black-hole spacetime that retains an approximately conserved energy, angular momentum and Carter constant along each geodesic orbit <cit.>. This model shows the potential of EMRIs for testing the Kerr solution in GR, as it can place much tighter constraints on the metric deformations than current X-ray observations <cit.>.

However, the prospect of doing precision science with EMRIs is accompanied by a high degree of technical difficulty. The length and complexity of EMRI signals leads to computationally expensive models and a highly multimodal likelihood surface in Bayesian inference problems, which exacerbates the already significant challenge of evaluating the evidence for model selection. Modern algorithms such as nested sampling can explore comparable-mass merger likelihoods efficiently <cit.> and are used to compute Bayes factors in the LIGO framework <cit.>, but improved techniques are required to adapt these algorithms for tests of GR with EMRIs.

In this paper, we take a first step towards developing a framework for testing GR with EMRI observations. The generalised EMRI waveform model described above is trialled in a Bayesian model selection framework based on the LIGO tests, using a nested-sampling algorithm that is tailored for high-dimensional and multimodal likelihoods <cit.>. We assess the viability of a product-space method <cit.> in accelerating the convergence of nested sampling on EMRI likelihoods, and apply a rethreading technique <cit.> that provides error estimates on the Bayes factors obtained through this method.

As the generalised EMRI model is too computationally unwieldy for method development purposes, we first present a proof of principle on a toy waveform model that mirrors several of its key qualitative features. The rethreading technique is empirically validated, and is used to compare the errors attained by the product-space method and regular nested sampling, as a function of total likelihood evaluations. We find that for the toy model likelihood, the product-space method reduces by an order of magnitude the number of likelihood calls taken to reach a standard deviation of 5% on the final GR/non-GR Bayes factor.

A similar but scaled-down analysis is then performed for the generalised model, with sampling restricted to a subset of EMRI parameters. Due to the greater complexity of the likelihood in this case, sampling bias is more likely to occur in shorter nested-sampling runs; we briefly demonstrate how a population of such runs may be combined to reduce this bias and ensure faster convergence in practice. Results for this more realistic EMRI likelihood show that the product-space method still attains better precision than regular nested sampling at the same level of computational cost, and indicate an even greater improvement in terms of cost to reach 5% error (potentially reducing the required number of likelihood calls by two orders of magnitude).

In Section <ref>, we give an overview of the generalised EMRI waveform model and define a toy surrogate that shares some of its relevant features. The statistical framework for testing GR with this model is set out in Section <ref>; we briefly summarise the key components of the LIGO test infrastructure, and describe the product-space method and rethreading technique that are used to adapt it for EMRI tests in this work. We then investigate in Section <ref> the viability of our proposed framework through a mock test of GR with the toy model, before extending our analysis to the generalised model in Section <ref>. Throughout this paper, we adopt geometrised units such that c=G=1. Greek (spacetime) indices run from 0 to 3, while Latin (space) indices run from 1 to 3. The base-10 logarithm is denoted by , while the natural logarithm is denoted by ln.



§ WAVEFORM MODELS



The generalised model considered in this work is based on the analytic kludge (AK) formalism of <cit.>, which we summarise in Section <ref>. This GR-based EMRI model is the most computationally efficient one available, and as such has been widely used in scoping out data analysis for space-based GW detectors; it is qualitatively resemblant to more accurate EMRI models, and its quantitative fidelity may also be improved through frequency corrections <cit.>. In Section <ref>, we describe the construction by <cit.> of AK waveforms on a family of generic modified-gravity black-hole spacetimes. These spacetimes are parametrised by metric deformations (or “bumps”) of different sizes, which show up in the resultant “bumpy AK” (bAK) model as perturbations to the phase evolution at different orders. The sinusoidal toy model we introduce in Section <ref> is designed to mimic the parameter dependence of the bAK phase evolution, at a fraction of the computational cost.



 §.§ Analytic kludge model



In the AK model, the instantaneous orbit of the compact object around the central black hole is approximated as Newtonian; the orbital parameters of this Keplerian ellipse are evolved over time with mixed-order PN expressions to simulate relativistic effects such as radiation reaction (inspiralling and circularising) and orbital precession (apsidal and Lense–Thirring). The waveform is then generated along the orbital trajectory using the <cit.> mode-sum approximation.

The inertia tensor for an EMRI with masses (μ,M≫μ) is given by I^ij(t)=μ x^i(t)x^j(t), where 𝐱 is the position of the compact object relative to the central black hole. For an instantaneous Newtonian orbit with orbital frequency ν, the second time derivative of I^ij may be decomposed into n-harmonics of ν as Ï^ij=∑_nÏ_n^ij. The three independent components of Ï_n^ij are given by

    Ï_n^11=a_n+c_n,  Ï_n^12=b_n,  Ï_n^22=c_n-a_n,

with <cit.>

    a_n=   -n𝒜[J_n-2(ne)-2eJ_n-1(ne)+2/nJ_n(ne)+2eJ_n+1(ne)
       -J_n+2(ne)]cos(nΦ),
    
    b_n=   -nA(1-e^2)^1/2[J_n-2(ne)-2J_n(ne)+J_n+2(ne)]
       ×sin(nΦ),
    
    c_n=    2AJ_n(ne)cos(nΦ),
    𝒜=    μν̃^2/3,

where the J_n are Bessel functions of the first kind, e is the orbital eccentricity, Φ(t) is the mean anomaly (i.e. Φ̇=2πν), and ν̃:=2πν M is the dimensionless orbital angular frequency.

At the detector location, it is convenient to work in an orthonormal coordinate frame {𝐩̂,𝐪̂,𝐫̂} with

    𝐩̂=𝐫̂×𝐋̂/|𝐫̂×𝐋̂|,  𝐪̂=𝐩̂×𝐫̂,

where 𝐫̂ points from detector to source and 𝐋 is the orbital angular momentum of the binary. In the transverse–traceless gauge, the leading-order gravitational radiation at the detector due to a source at luminosity distance D is given by <cit.>

    h_ij=2/D(P_ikP_jl-1/2P_ijP_kl)Ï^kl,   h^+,×=1/2h^ijH_ij^+,×,

with the polarisation and projection tensors

    H_ij^+=p̂_ip̂_j-q̂_iq̂_j,   H_ij^×=p̂_iq̂_j+q̂_ip̂_j,   P_ij=δ_ij-r̂_ir̂_j,

where δ_ij is the Kronecker delta.

From (<ref>)–(<ref>), the two GW polarisation amplitudes h^+,× for an EMRI may be written in terms of the Peters–Mathews harmonic decomposition as <cit.>

    h^+,×=1/D∑_n A^+,×_n,
       A^+_n=C^+_aa_n+C^+_bb_n+C^+_cc_n,
       A^×_n=C^×_aa_n+C^×_bb_n,

where the coefficients C^+,×_a,b,c depend on the angular configuration (λ,γ̃,α,θ_K,ϕ_K,θ_S,ϕ_S) of the EMRI: the inclination λ between 𝐋 and the black-hole spin 𝐒; the azimuth γ̃ of periapsis in the orbital plane (relative to 𝐋×𝐒); the azimuth α of 𝐋 projected onto the spin-equatorial plane; the orientation (θ_K,ϕ_K) of 𝐒; and the sky location (θ_S,ϕ_S). The first two angles are intrinsic to the source, while the rest are defined relative to a fixed ecliptic-based coordinate system <cit.>. Explicit expressions for the coefficients are given by Eqs. (10) and (18)–(25) in <cit.>.

In a relativistic EMRI, the frequency ν, eccentricity e and inclination λ change over time due to radiation reaction, while the angles (γ̃,α) precess as well. The AK model evolves (ν,e) with 3.5PN expressions, while approximating λ as constant due to its slow evolution over the full inspiral <cit.>. The fluxes (ν̇,ė) are given by <cit.>

    ν̇=96/10πην̃^11/3/M^2(1-e^2)^7/2(1+73e^2/24+37e^4/96)+𝒪(ν̃^13/3),
       ė=-304/15ην̃^8/3e/M(1-e^2)^5/2(1+121e^2/304)+𝒪(ν̃^10/3),

where η:=μ/M is the mass ratio, and only the leading 2.5PN terms are explicitly presented here. (See Eqs. (28) and (30) in <cit.> for the full expressions.)

The angular rates (γ̇̃̇,α̇) determine the precession of periapsis in the orbital plane, γ̇̃̇+α̇, and the Lense–Thirring precession of the orbital plane around the black-hole spin axis, α̇. They are given by <cit.>

    γ̇̃̇=3ν̃^5/3/M(1-e^2)+𝒪(ν̃^2),
       α̇=2ãν̃^2/M(1-e^2)^3/2,

where ã:=|𝐒|/M^2 is the dimensionless spin angular momentum. (See Eq. (29) in <cit.> for the full version of (<ref>).)

In the original AK model, the two precession rates and the orbital frequency can differ significantly from their actual values in a relativistic EMRI. Following <cit.>, we correct the starting angular rates (Φ̇_0,γ̇̃̇_0,α̇_0) through a parameter-space map (M,ã,ν)↦(M',ã',ν') such that

    Φ̇_0(M',ã',ν')=ω_r(M,ã,ν),
       γ̇̃̇_0(M',ã',ν')=ω_θ(M,ã,ν)-ω_r(M,ã,ν),
       α̇_0(M',ã',ν')=ω_ϕ(M,ã,ν)-ω_θ(M,ã,ν),

where ω_r,θ,ϕ are the fundamental frequencies of radial, polar and azimuthal motion <cit.> on the starting Kerr geodesic. This map does not eradicate the accumulated phase error of AK waveforms over the full inspiral, but has negligible computational cost and greatly improves the quantitative accuracy over short times.

An EMRI has 14 degrees of freedom in the AK model, which neglects the spin of the compact object; the parameters of the model are chosen to decouple the seven source-intrinsic degrees of freedom from the observer-dependent extrinsic ones. We define the set of dimensionless AK parameters as

    Θ_AK=Θ_int∪Θ_ext,
       Θ_int={(μ/M_⊙),(M/M_⊙),ã,e_0,cosλ,Φ_0,γ̃_0},
       Θ_ext={ν̃_0,α_0,cosθ_K,ϕ_K,cosθ_S,ϕ_S,(D/Gpc)},

where the subscript zero denotes the value taken by a quantity at the arbitrary starting time t=0.



 §.§ Bumpy analytic kludge model



The family of generically deformed black-hole spacetimes derived separately by <cit.> and <cit.> provides a model-independent setting in which to construct modified-gravity EMRI waveforms for testing the Kerr metric solution in GR. These bumpy black holes are not required to satisfy the Einstein equations (as done in <cit.>), but are Kerr-like through their stationarity, axisymmetry, and admission of an approximate second-rank Killing tensor. They are also more general than the modified Kerr spacetimes considered in other proposed EMRI tests <cit.>, as they allow for mass-moment deformations beyond quadrupole order.

In Boyer–Lindquist coordinates, the components of the Kerr metric around a black hole with mass M and spin angular momentum a=ãM are given by <cit.>

    g^K_tt=-(1-2Mr/Σ),   g^K_tϕ=-2Marsin^2θ/Σ,   g^K_rr=Σ/Δ,
       g^K_θθ=Σ,   g^K_ϕϕ=(r^2+a^2+2Ma^2rsin^2θ/Σ)sin^2θ,

where Σ:=r^2+a^2cos^2θ and Δ:=r^2-2Mr+a^2. The Kerr metric is stationary and axisymmetric, and hence may be cast in Lewis–Papapetrou form via the partial coordinate transformation

    (ρ,z)=(√(Δ)sinθ,(r-M)sinθ),

with the temporal and azimuthal coordinates (t,ϕ) left unchanged.

<cit.> consider a linear deformation of the Kerr metric in Lewis–Papapetrou form, and apply the inverse of the transformation (<ref>) to obtain its components in Boyer–Lindquist-like coordinates (r,θ); this ensures that the deformed metric remains stationary and axisymmetric, i.e. it admits temporal and azimuthal Killing vectors t^μ and l^μ such that

    ∇_(μt_ν )=∇_(μl_ν )=0.

The deformed metric components may then be written as

    g_μν=g^K_μν+ϵ h_μν,

where ϵ≪1 is a bookkeeping parameter for the metric deformation h_μν. Other Kerr-like properties are included by requiring the second-rank tensor ξ_μν=Δ t_(μl_ν)+r^2g_μν to approximately satisfy the Killing tensor equation, i.e.

    ∇_(λξ_μν )=𝒪(ϵ^2),

and by requiring |h_μν|=𝒪(1/r^2) as M/r→0 (such that the deformed metric retains asymptotic flatness, along with its original mass and spin angular momentum).

With the above constraints, the only nonzero components of h_μν are h_tt, h_tϕ, h_rr and h_ϕϕ, which depend on the black-hole parameters (M,ã) and three arbitrary radial functions γ_i,i∈{1,3,4}. By representing the nonzero components and radial functions as power series in M/r, <cit.> obtain expressions for h_μν,n,2≤ n≤5 in terms of γ_i,n, where these quantities are the coefficients of the (M/r)^n terms in the corresponding series. If the inclination angle of a geodesic orbit is approximated as constant (as done in the AK model), a metric deformation that is purely[For simplicity, the allowed deformations in the bAK model are restricted to “pure” bumps and their linear combinations, e.g. a ℬ_4 bump is not the most general 4th-order deformation, but defined as one for which the ℬ_2,3,5 coefficients are all zero. While the model might not fully represent deviations from GR at sub-leading order, it is an adequate surrogate in the present work.] nth-order in M/r turns out to be fully specified by a set of three coefficients, ℬ_n:={γ_1,n,γ_4,n,γ_3,n+1}; we refer to such a deformation as a ℬ_n bump.

The AK formalism is then used to construct EMRI waveforms on the bumpy black-hole spacetime, which is parametrised by the set of coefficients ⋃_nℬ_n,2≤ n≤5 in addition to (M,ã). Both the long-timescale radiation reaction fluxes and the short-timescale precession rates of the EMRI are altered by the ℬ_n coefficients, since they perturb the three first integrals of motion along a timelike geodesic with four-velocity u^μ: the energy E=t_μ u^μ and orbital angular momentum L_z=l_μ u^μ, which remain conserved, and the analogue of the Carter constant Q=ξ_μνu^μ u^ν, which is conserved at linear order in the metric deformation.

By matching the turning points of motion for an instantaneous geodesic orbit (E,L_z,Q) with those of a precessing Keplerian orbit (ν,e,λ) in the AK model, <cit.> compute the leading-order corrections caused by the ℬ_n deformations to the fluxes and angular rates (<ref>)–(<ref>). Each set of corrections at nth-order in M/r depends only on a single linear combination ϵ_n:=γ_1,n+2γ_4,n of the ℬ_n coefficients (γ_3,n+1 is at sub-leading order), which effectively reduces the additional degrees of freedom in the extended model to four. The corrections are given by

    δν̇_n=8/5πην̃^(2n+9)/3/M^2(1-e^2)^n+5/2g_ν,nϵ_n,
       δė_n=-16/5ην̃^(2n+6)/3/M(1-e^2)^n+3/2g_e,nϵ_n,
       δγ̇̃̇_n=ν̃^(2n+1)/3/M(1-e^2)^n-1g_γ̃,nϵ_n,
       δα̇_n=-α̃ν̃^(2n+2)/3/M(1-e^2)^n-1/2g_α,nϵ_n,

where explicit expressions for the eccentricity-dependent factors g_·,n are given by Eqs. (301), (302) and (335)–(346) in <cit.> (with g_γ̃,2=1/2 and g_α,2=1).



In summary, the bAK model comprises i) adding the corrections (<ref>)–(<ref>) to the corresponding evolution equations (<ref>)–(<ref>) in the AK model, and ii) extending the set of model parameters to

    Θ_bAK=Θ_AK∪Θ_ℬ,
       Θ_ℬ={ϵ_n | 2≤ n ≤5},

where we assume ϵ_n>0 for simplicity. The deformation parameters ϵ_n then determine the magnitudes of the ℬ_n bumps, which manifest as phase drifts in the corresponding bAK waveforms. Figure <ref> illustrates how these waveforms dephase over time relative to the AK waveform (where ϵ_n=0 for all n). Full EMRI waveforms can have up to ∼10^5 observable cycles, and hence might be able to constrain the leading-order bumps down to ϵ_n∼10^-7 <cit.>. For the sake of comparison, combined results for the LIGO events GW150914 and GW151226 in a less conservative single-parameter analysis placed no upper bound below ∼10^-1 on the dimensionless deformation parameters of the generalised black-hole binary merger model <cit.>.



 §.§ Sinusoidal toy model



Motivated by the phase behaviour of the bAK waveforms, we define a deformed sinusoidal model in which the phase evolution has similar dependence on a set of deformation parameters analogous to Θ_ℬ. Waveforms from this toy model are several orders of magnitude faster to compute than bAK waveforms, and as such are useful for building intuition during method development. They are qualitatively equivalent to the special case of bAK waveforms from circular, equatorial EMRIs over short timescales (or alternatively, in the geodesic limit η=0).

Our toy model consists simply of a sine wave parametrised by amplitude A and angular frequency Ω, along with four parameters that deform the phase at different orders in time; these deformation parameters are denoted ϵ_n by way of analogy to the bAK model. The toy waveform is given by[Although full GW models require two independent waveform modes for parameter estimation, the likelihood obtained with the single toy model mode is qualitatively sufficient for the purposes of this work.]

    h=Asin[Ωt̃(1+∑_n=2^5ϵ_n(t̃/τ)^n-1)],

where Ω, t̃ and τ are all dimensionless. The timescale-like quantity τ has been introduced to control the overall strength of the phase drift, and is not treated as a parameter in the model. Henceforth we fix τ=2T_obs, where T_obs is some specified waveform duration, and define the set of toy model parameters as

    Θ_toy={A,Ω}∪Θ_ℬ,

where Θ_ℬ is defined as in (<ref>). 



§ STATISTICAL FRAMEWORK



The framework presented in this paper is based on the model-independent and parametric test infrastructure introduced by <cit.>, in which the combinatorial submodels of a generalised waveform model are used in a null-hypothesis test of GR, and a nested-sampling algorithm is employed to compute the Bayesian evidence for each submodel. An overview of the pertinent concepts is given in Sections <ref>–<ref>. In Section <ref>, we describe the product-space nested-sampling method considered by <cit.>; it provides an evidence-free computation of the Bayes factor for each submodel (with respect to the null submodel), thus mitigating the prohibitive computational difficulty of performing the test with a generalised EMRI waveform model. The estimation of Bayes factor errors in the method also requires a new rethreading technique <cit.>, which we outline in Section <ref>.



 §.§ Gravitational-wave likelihood



In the basic matched-filtering framework for GW data analysis, data from a single interferometric detector is written as the time series x=h+n, where h is the detector response to a passing GW and n is the detector noise (typically approximated as a Gaussian and stationary random process). For a LISA-like detector with three arms, two independent signals h_I,II may be obtained; these are related on average to the two GW polarisations h^+,× by

    h_I,II=√(3)/2(F_I,II^+h^++F_I,II^×h^×),

where the antenna pattern functions F_I,II^+,× <cit.> depend on the sky location and polarisation angle of the source in a detector-based coordinate system. (See Eqs. (15)–(17) in <cit.> for the explicit expressions in ecliptic coordinates.) Doppler modulation of the waveform phase may also be included in h^+,× to account for the orbital motion of LISA.

For compactness, we write the GW signal as a complex time series h=h_I+ih_II.[This notation is also compatible with the single real waveform mode h of the sinusoidal toy model.] The signal-to-noise ratio (SNR) of h is given by ρ:=√(⟨ h|h⟩)=√(⟨ h_I|h_I⟩+⟨ h_II|h_II⟩), with the noise-weighted inner product ⟨·|·⟩ defined as <cit.>

    ⟨ a|b⟩=2∫_0^∞ df ã^*(f)b̃(f)+ã(f)b̃^*(f)/S_n(f),

where S_n is the power spectral density of n. If the detector noise is assumed to be white (as done in Section <ref>), (<ref>) simplifies to

    ⟨ a|b⟩=1/S_n∫_0^T_obsdt a^*(t)b(t),

which may be computed directly from the time-domain waveforms. In Section <ref>, S_n is given instead by an analytic approximation <cit.> to the LISA noise power spectral density.[The noise model used corresponds to a down-scoped version of LISA, and thus gives more conservative results; the mission design has now been restored to an earlier configuration with higher sensitivity.] Throughout this paper, we work with waveforms that are renormalised with respect to some reference waveform h' and specified amplitude ρ', i.e. h→ρ'h/√(⟨ h'|h'⟩); hence it is more useful to think of “SNR” here as an amplitude relative to the defined norm in (<ref>), and not to the actual noise in the data.

Given GW data x that contains a signal h(θ_*) corresponding to the model parameter values θ_*, the Bayesian likelihood =(x|θ) is defined as <cit.>

    (θ)∝exp(-1/2⟨ x-h(θ)|x-h(θ)⟩).

We consider only EMRI waveforms with ρ≳10 throughout this paper, which is a conservative choice that is consistent with typical values of the threshold SNR for reliable detection <cit.>. In this regime, the noise term in the log-likelihood is suppressed by a factor of ≳100 relative to the leading term, and may be neglected for the purposes of model selection. Hence, to simplify analysis in Sections <ref>–<ref>, we assume a particular noise realisation of n=0 in the synthetic data (i.e. x=h(θ_*)).



 §.§ Null-hypothesis test



From Bayes' theorem, the posterior probability of a model hypothesis  given the data x may be written as

    (|x)=Π/(x),

where =(x|) and Π=() are, respectively, the evidence and prior probability for the model. The model evidence is typically obtained by marginalising the likelihood =(x|θ,) over the model parameters, i.e.

    =∫ dθ (θ)π(θ),

where π=(θ|) is the parameter prior.

Two model hypotheses _i,j may be compared through the ratio of their posterior probabilities, which quantifies the degree of belief in one model over the other. The logarithm of this posterior odds ratio is defined as

    P^i_j:=ln[(_i|x)/(_j|x)]=B^i_j+ln(Π_i/Π_j),

where B^i_j:=ln(_i/_j). The quantity B^i_j is the logarithm of the Bayes factor, which is commonly used as an equivalent substitute for the posterior odds ratio through the implicit assumption Π_i=Π_j. A general scale for the interpretation of posterior odds ratios <cit.> is given in Table <ref>.



For a generalised waveform model that extends a GR-based model _GR through a set Θ_ℬ of N deformation parameters, we may define 2^N submodels whose deformation parameter sets are given by the 2^N subsets of Θ_ℬ. The submodel corresponding to the null set is _GR itself; the remaining 2^N-1 submodels are modified-GR models, which we denote collectively by the hypothesis _modGR=⋁_m≠GR_m. Even though _modGR is a collection of nested submodels, <cit.> observe that the individual submodel evidences are logically disjoint, due to the distinct integration measure on each submodel parameter space. Hence the posterior odds ratio for _modGR over _GR simplifies to

    (_modGR|x)/(_GR|x)   =∑_m≠GR(_m|x)/(_GR|x)
       =1/2^N-1∑_m≠GR_m/_GR,

where the second equality follows from the assumptions that Π_modGR=Π_GR and Π_m=Π_m' for all m,m'≠GR.

This approach provides the basic framework for a test of GR with the bAK model, where the modified-GR hypothesis is compared against the null hypothesis (the AK model) through evaluation of the submodel Bayes factors in (<ref>). A single EMRI source is considered for the assessment of our methods, but it is straightforward to generalise (<ref>) for a population of sources <cit.>. In the bAK model (where N=4), the 16 submodels may be indexed by 0≤ m≤15, whose nonzero digits in binary representation specify the deformation parameters included in each submodel. Hence the AK model is denoted _0≡_0000, while the full bAK model is _15≡_1111. We use the convention that the binary digits of m from right to left correspond to ϵ_n with n=2,3,4,5 respectively; for example, _1010 has only two deformation parameters {ϵ_3,ϵ_5}, and is equivalent to the ϵ_2=ϵ_4=0 slice of the full model.

Although the Bayesian evidence (<ref>) has a built-in Occam penalty on model complexity, it does not fully account for the multiplicity effect <cit.>. As a larger number of deformation parameters is considered, it becomes more likely that one particular parameter will cause the inclusive submodels to fit the data well by chance. In other words, a uniform submodel prior (i.e. Π_m=Π_m' for all m,m'≠0) might not be appropriate for the nested submodels in this framework, and more thorough prescriptions for assigning model prior probability <cit.> should be investigated. However, the number of deformation parameters in the bAK model is small; furthermore, the focus of this work is the evaluation of (<ref>) with improved precision, and not the validation of its accuracy. Hence we retain a uniform submodel prior for simplicity, and it follows from (<ref>) and (<ref>) that

    =B^modGR_GR=ln[∑_m=1^15exp(B^m_0)]-ln15.




 §.§ Nested sampling



For most statistical problems, (<ref>) admits no analytic solution and is impractical to evaluate through direct numerical integration, even over a parameter space of modest dimensionality. A wealth of alternative techniques for computing the evidence have thus been developed; these range from simple estimates that use Laplace's approximation <cit.> or Chib's method <cit.>, to more sophisticated sampling strategies based on concepts such as thermodynamic integration and simulated annealing <cit.>. We employ in this work the nested sampling strategy introduced by <cit.>, which provides an accurate and computationally efficient means of simultaneously exploring the posterior surface and evaluating the evidence. Nested sampling has been shown to be suitable for GW likelihoods <cit.>, and is one of the two algorithms used to compute model evidences in LIGO tests of GR <cit.>.

In nested sampling, the multidimensional integral (<ref>) is written in the one-dimensional form <cit.>

    =∫_0^1dX (X),

where the prior mass X is given by

    X(λ)=∫_(θ)>λdθ π(θ),

i.e. the integral of the parameter prior over the interior of the likelihood contour =λ. A set of  initial points is first sampled from the prior; this set of “live” points is then evolved across parameter space by discarding at the i-th iteration the point θ_i with the lowest likelihood value λ_i, and replacing it with one drawn from the prior but within the contour _i=λ_i.

<cit.> shows that the prior mass corresponding to each λ_i may be approximated as

    X_i≈exp(-i/),

which shrinks exponentially from X_0:=1 to zero as the live points converge on and navigate the bulk of the posterior surface. Upon the satisfaction of suitable convergence criteria, the algorithm is truncated and the evidence (<ref>) may be approximated as

    ≈∑_i>0w_i_i,

where the weights w_i are given by a Riemann sum rule, e.g. w_i=X_i-1-X_i. The set of discarded live points θ_i (often termed dead points) also serves as a set of posterior samples, through the assignment of posterior probability p_i=w_i_i/ to each point.

The prior-mass approximation (<ref>) introduces statistical error into the posterior probabilities p_i and the evidence estimate (<ref>), via the weights w_i. The uncertainty in the evidence estimate depends on the absolute error of each w_i, and is dominated by the Poisson variability in the number of iterations taken to reach the posterior bulk. Hence (<ref>) is log-normally distributed, and the standard deviation of ln in nested sampling scales with the number of live points as <cit.>

    σ_ln∝1/√().


It is often difficult to sample from the prior under the constraint >λ, since the likelihood contours =λ might in general be multimodal or degenerate. In this work, we make use of the nested-sampling implementation PolyChord <cit.>, which mitigates these difficulties through the incorporation of clustering and slice sampling algorithms. It also exhibits good scaling with the dimensionality of the parameter space, and in that sense is an improved successor to the widely adopted nested-sampling tool MultiNest <cit.>. The two main runtime parameters in PolyChord that determine sampling resolution and reliability are, respectively,  and ; the latter is the number of randomly oriented one-dimensional slices sampled at each iteration in order to decorrelate the new live point from the discarded point.



 §.§ Product-space nested sampling



Methods exist for obtaining Bayes factors without explicitly evaluating evidences, and these are especially useful when the number of competing models is large. The most well known is the Savage–Dickey density ratio for nested models (and generalisations thereof <cit.>), where the Bayes factor between a null model and an encompassing one is given by the ratio between the posterior and prior probabilities at the null point in the larger model space. Another class of methods uses the product-space representation <cit.> for a prespecified and indexed collection of competing models; this approach involves exploring the set of model indices and each model space simultaneously, and has been investigated in the context of Markov chain Monte Carlo (MCMC) algorithms <cit.> and nested sampling <cit.>.

Product-space sampling and reversible-jump MCMC <cit.> are modern examples of transdimensional frameworks for Bayesian model selection <cit.>. One disadvantage of product-space methods (as compared to reversible-jump MCMC) is the required prior specification of all the competing models. This is not an issue in our framework since the bAK submodels are fully defined by the power set of Θ_ℬ, and the ℬ_n bumps are truncated at n=5 (their effects are exponentially suppressed as n increases). In the product-space representation, the 16 submodels _m in Section <ref> may be thought of as distinct “slices” of some hypermodel ; the parameter space of  is the combination[More precisely, the hypermodel parameter space is the vector bundle of submodel spaces over the discrete set of model indices.] of all the submodel spaces, and is parametrised by

    Θ=Θ_gen∪{m},

where Θ_gen is the parameter set for the generalised waveform model (i.e. (<ref>) or (<ref>)) and m is the submodel index.

The posterior probability for m is given by <cit.>

    (m|x,)   =∫ dθ (θ,m|x,)
       =1/_∫ dθ (θ,m)π(θ,m)
       =π(m)/_∫ dθ_m (θ_m)π(θ_m|m)
       =π(m)_m/_,

where _ and _m are the evidences for the hypermodel and mth submodel respectively. To obtain the third equality, we have decomposed θ into the parameters θ_m that are included in the submodel _m, and the parameters ϕ_m that are excluded; the integral ∫ dϕ_m π(ϕ_m)=1 then factors out of the expression.

When sampling in the hypermodel space, it is convenient to assign a uniform prior π(m)=1/16 on the submodel index, since it is straightforward to restore the assumption ∑_m≠0Π_m=Π_0 in post-processing (i.e. through the final term in (<ref>)). We then have

    B^m'_0=ln(_m'/_0)=ln[(m=m'|x,)/(m=0|x,)],

such that the Bayes factors in (<ref>) may be obtained directly from the sampled posterior distribution for m. A nested-sampling implementation of the product-space method was found to be viable by <cit.> through application to a cosmological problem. In this paper, we use a similar implementation to compare a larger number of models in our GW test of GR, and perform an improved investigation of the errors on the obtained Bayes factors.



 §.§ Error estimation from rethreading



If the Bayes factors in (<ref>) are obtained by evaluating the individual submodel evidences with regular nested sampling, the standard deviation of  scales as 1/√() (from propagation of the log-evidence error (<ref>)). In the product-space method, however, the Bayes factors are ratios of posterior probabilities p_m:=(m|x,); computing them is then a parameter estimation problem, with a different associated uncertainty that depends on the posterior errors σ_p_m. These errors arise from the prior-mass approximation (<ref>), as in the case of σ_ln, but also from the dimensional reduction in the likelihood reparametrisation (θ)→(X) <cit.>. Nevertheless, <cit.> observe that σ_p_m are determined by the relative errors of the nested-sampling weights, and hence are typically <σ_ln. This is the key factor that underpins the improved efficiency of product-space nested sampling.

An algorithm for estimating the errors of p_m (or of any quantity derived from p_m) over a single nested-sampling run has been proposed by <cit.>. The technique is motivated by the fact that any number of runs r with N_r live points may be merged into a single run with =∑_rN_r <cit.>, by combining and reordering all of their dead points (and adjusting  in (<ref>) accordingly). Conversely, it is also possible to unravel a single nested-sampling run into its constituent “threads”, i.e. a set of  independent runs r with N_r=1. This is achieved by tracking the birth and death order of the sampling points; each thread is then formed from the sequence of replacements for an original live point.

A distribution of runs with the original number of live points may be obtained rapidly from the set of threads, by using bootstrap resampling on the set and recombining (or “rethreading”) the sample threads. These new runs contain only points that are present in the original run, but yield statistical variance in the estimated weights and p_m. In this work, we use the rethreading technique to generate 10^3 realisations of  from a single nested-sampling run, and thus to evaluate its mean and standard deviation directly. The rethreading results are also validated against those from repeated runs for the toy model in Section <ref>.



§ RESULTS: SINUSOIDAL TOY MODEL



To demonstrate that product-space nested sampling with rethreading can explore qualitative EMRI likelihoods with improved efficiency, we first apply it to synthetic data generated from the toy EMRI model defined in Section <ref>. Two data sets x are considered: one from the “GR” submodel _0000, and one from the deformed submodel _0010 with a “ℬ_3 bump”. The signal parameters for x_0000 and x_0010 are (A,Ω)=(1,1) and (A,Ω,ϵ_3)=(1,1,-1.9) respectively. Both waveforms are generated with duration T_obs=10^4 (such that they contain ∼10^3 cycles) and a dimensionless sampling rate of unity, then renormalised to SNR ρ=10.





As the toy waveforms are near-sinusoidal, the GW likelihood (<ref>) contains damped oscillations in Ω; to a lesser extent, such oscillations are also present for parameters that are highly correlated with frequency (e.g. M) in more realistic EMRI models. Figure <ref>(a) shows the ϵ_3–Ω slice of ln for the submodel _0010 with the data x_0000, where Ω is seen to be relatively well localised in the moderate-SNR regime. However, the likelihood is highly degenerate in the deformation parameters. In the case of Figure <ref>(a), any value for ϵ_3 that falls below some threshold ≈-2 ceases to deform the waveform significantly, and hence has negligible effect on the value of ln. This results in an extended region of high likelihood along the ϵ_3 axis, even though the true GR signal lies on the ϵ_3=0 boundary of the _0010 space.

The deformation parameters are also degenerate among themselves, which further complicates the likelihood surface if the true signal is deformed. For the data x_0010, the signal is not even contained in submodels without ϵ_3; however, Figure <ref>(b) shows that regions of high likelihood exist in the ϵ_4–ϵ_5 slice of ln for _1100, and that a ϵ_3=-1.9 deformation can be well approximated by a larger value of ϵ_4 or ϵ_5 (or a combination of the two). Although such degeneracies are partially broken by the additional degrees of freedom in more realistic generalised GW models, they are still expected to hamper parameter estimation for the individual deformations <cit.>. Nevertheless, they have a less severe impact on the null-hypothesis test in Section <ref>, where the aim is not to distinguish among the deformed submodels.

For the two data sets, both regular and product-space nested sampling are used to obtain  with an associated standard deviation σ_P. In the regular case, (<ref>) is evaluated piecewise by computing the individual submodel evidences, and the standard nested-sampling estimates for each log-evidence error (<ref>) are propagated to yield σ_P. For the product-space method,  is obtained through a single nested-sampling run in the hypermodel space, and its approximate error is calculated in two ways: using the single-run rethreading technique, and from 50 repetitions of the same product-space run. The latter procedure is only to demonstrate that the methods in Sections <ref> and <ref> are working as expected, and is not performed in Section <ref>.

The product-space method manifestly provides computational savings over regular nested sampling in the null-hypothesis test, as it samples in just one parameter space instead of 16. However, the hypermodel space is more complex than a single submodel space and requires additional runtime to explore effectively, such that differences in efficiency cannot simply be estimated from the number of spaces sampled. To assess the gains when using the product-space method, we vary the PolyChord runtime parameter  to obtain different degrees of precision on each corresponding evaluation of , and compare the number of likelihood calls taken by the two methods to achieve the same σ_P. The other PolyChord runtime parameter is set as =30, which has been chosen empirically to ensure the convergence of regular evidence estimates.

It is also instructive to study how regular and product-space nested sampling perform on the penultimate step in the evaluation of , i.e. the individual submodel Bayes factors B^m_0. Figure <ref> shows the B^m_0 and associated errors that are obtained from the two methods with =2500, for the GR data x_0000. The Occam penalty on model complexity is clearly observed in both sets of results, as the relative evidence for each submodel decreases with the number of parameters it contains. However, the product-space method appears to systematically give Bayes factors that are more pronounced (negative), and in tension with the regular results. This is likely because the entire hypermodel space is explored with the same number of live points allocated to each submodel space in the regular method, leading to a slight degree of sampling bias. The Bayes factor errors for both methods are nevertheless comparable, since the errors on a posterior over m are smaller than those on submodel evidence evaluations (as discussed in Section <ref>).





In Figure <ref>, the error σ_P of  for the two data sets x_0000 and x_0010 is plotted against the number of likelihood calls for a sequence of regular and product-space nested-sampling runs with 100≤≤2500. Both the rethreading and repetition error estimates for the product-space method are included; they are seen to agree well, with the latter showing more scatter (since they are computed from only 50 evaluations of , as opposed to 10^3 realisations in the rethreading technique). The rethreading error estimates are further validated through a reduced chi-squared test against the sample mean μ of  in Figure <ref>, where χ^2=∑_i(P_i-μ)^2/(12σ_P^2)≈1 for both data sets. We also find μ_0000∼-1 in the GR case and μ_0010∼1 in the ℬ_3 case, which is by design from our choices of SNR and ϵ_3 for the synthetic data.

For both data sets, it is clear that product-space nested sampling is effective at reducing the computational cost required to reach a given level of precision (or alternatively, at providing greater precision with a given number of likelihood calls). In the GR case, the average gain in efficiency (i.e. the mean horizontal distance between the blue and green curves in Figure <ref>) is a factor of around 24. Furthermore, as the likelihood surface over the hypermodel space is less complex for x_0000, nested sampling explores it nearly as efficiently as each of the 16 submodel spaces. This is seen by comparing regular and product-space runs of equal , where the number of likelihood calls taken by the latter is almost exactly 16 times smaller, and its associated error is slightly lower as well.

In the ℬ_3 case, larger overall errors σ_P are obtained for both methods, and the average reduction in computational cost with product-space nested sampling is reduced to a factor of around 9. As expected, this is largely caused by the increased complexity of the likelihood surface over the hypermodel space for x_0010: additional modes are present in the parameter spaces of the eight submodels containing ϵ_3, as well as in other submodel spaces due to deformation parameter degeneracies such as that from Figure <ref>(b). When comparing regular and product-space runs of equal , the number of likelihood calls taken by the latter is only 13 times smaller, and its associated error is now slightly higher. It follows from these results that computational savings will likely be smallest for data generated from _1111; we discuss this in Section <ref>.

Another contribution to the overall difference in σ_P between the two data sets arises from the construction of  itself. For the ℬ_3 data, the parameter space of the null submodel _0000 is a region of lower posterior probability in the product-space approach, and hence the relative posterior error on (m=0) is higher. This error propagates into every submodel Bayes factor B^m_0 via (<ref>), which increases the final error σ_P. The effect is also present for regular nested sampling, since _0000 has a higher relative evidence error as well. However, the increase in σ_P only becomes significant in the strong-deformation regime, which is unlikely for actual tests of GR; furthermore, it can also be mitigated in practice (e.g. by using the actual submodel index prior ∑_m≠0Π_m=Π_0 when sampling).



§ RESULTS: BUMPY ANALYTIC KLUDGE MODEL



For a more realistic generalised EMRI waveform model in the product-space approach, the dimensionality of the hypermodel parameter space is considerably higher. In the case of the bAK model, Θ=Θ_AK∪Θ_ℬ∪{m} is 19-dimensional and parametrised by 14 GR parameters, four deformation parameters, and the submodel index. The likelihood surface over the GR parameter space is known to be highly multimodal with a large information content <cit.>. A full exploration of this space is hampered by the significant computational cost of waveform generation (even though the AK formalism already provides the cheapest EMRI waveforms available), and is beyond the scope of the present work. We instead fix all but seven of the hypermodel parameters, allowing only the component masses, deformation parameters and submodel index to vary. For our synthetic data, the intrinsic GR parameters (<ref>) of the signal are chosen as Θ_int=(1,6,0.1,0.1,0.9,0,0); the waveform is two months long with an initial frequency of 2 mHz (such that it contains ∼10^4 cycles), and is sampled at 0.2 Hz.



Instead of considering synthetic data from a deformed submodel as in Section <ref>, we restrict analysis here to the more realistic case of a GR signal, and investigate the effect of SNR on sampling performance. The two data sets studied in this section are both generated from _0000 with the same GR parameters, but are renormalised to SNRs ρ=10 and ρ=100 respectively. Figure <ref> shows the B^m_0 and associated errors that are obtained from regular and product-space nested sampling with =500, for the moderate-SNR case ρ=10. The Occam penalty is again apparent, but the submodel Bayes factors are generally lower than in Figure <ref> even though the SNR is the same; this is likely due to the EMRI's orbital evolution reducing degeneracy in the deformed parameters. As in Figure <ref>, there also appears to be a slight systematic difference between the Bayes factors from the two methods. Both the regular and product-space methods correctly favour the null submodel _0000, although the latter does so with smaller errors and fewer likelihood calls for the same number of live points.



In Figure <ref>, the error σ_P of  for the ρ=10 and ρ=100 data sets is plotted against the number of likelihood calls for several product-space nested-sampling runs with varying PolyChord runtime parameters  and . These are compared against a single regular nested-sampling run for each SNR value (due to the considerably higher computational cost of the regular method). It is clear that the efficiency gains obtained for the sinusoidal toy model are still present for the bAK model, in that every product-space run shown has both better precision and lower computational cost than the two regular runs. As indicated by the results in Section <ref>, these gains might be diminished in the case of data that is generated from more complex submodels such as _1111. However, since there is strong prior expectation for an actual EMRI signal to be well described by GR, the order-of-magnitude savings observed here will likely be close to what is obtained in practice.

Furthermore, since the log-evidence error (<ref>) scales as 1/√() and the number of likelihood calls increases approximately linearly with , an extrapolation of both regular nested-sampling runs to 5% error (∼10^8 calls) indicates that for the bAK likelihood, the reduction of computational cost from the product-space method is boosted to around two orders of magnitude. A direct verification of this statement is impractical for the present work, since the cost of each ∼10^7-call run in Figure <ref> is ≈4000 core hours (a single call to the bAK likelihood takes ≈1.5 s). Nevertheless, the scaling of errors on nested-sampling evidences is both well understood and reliable (as indicated by the smoothness of the blue curves in Figure <ref>), which lends credence to the validity of such an extrapolation.



Significant sampling bias is introduced into the results for the high-SNR case ρ=100; this is because optima in the likelihood surface become more localised, and hence require higher sampling resolution/reliability to map out accurately. The values of  and their associated errors from the 10 high-SNR product-space runs (indicated by the light green points in Figure <ref>) are plotted in Figure <ref>. Computing the reduced chi-squared statistic for these runs gives χ^2≈5, such that σ_P for any given run appears to underrepresent the observed scatter within the set of runs. Since the rethreading technique has been validated in Section <ref>, the large χ^2 suggests that the hypermodel space is not being sampled consistently between runs. If any particular run misses a region of high posterior probability, then its threads will not contain enough information to provide a rethreading estimate that represents its true error; this could occur if  is not large enough to find sharply defined modes, or if replacement live points are correlated with discarded ones due to inadequate . The sampling bias observed here (and to a lesser extent in Figures <ref> and <ref>) might be better characterised with new diagnostic tests <cit.> in future work.

In Figure <ref>, we illustrate through a clustered chi-squared test that the obtained  and σ_P actually show a closer fit to two distinct sampling distributions with differing mean (since χ^2≈1 for each cluster). This analysis is not to be taken at face value, since  clearly has a unique underlying value. It does however highlight the possibility of drawing an erroneous conclusion from the null-hypothesis test due to sampling bias, especially if the runtime parameters  and  are set too low to adequately explore high-SNR likelihoods. Furthermore, with the added complexity of the EMRI likelihood surface, it is computationally unfeasible to systematically increase these parameters without bound. We thus propose a strategy for making the product-space method robust to sampling bias; the idea is simple, and made possible by the same premise that facilitates the rethreading technique.

As discussed in Section <ref>, nested sampling permits both the unravelling and interweaving of independent runs. Although threads from a single run might have correlations that reflect the sampling bias present in that run, such correlations are more diffuse in a set of threads from a large population of runs, and sampling bias will be reduced in runs that are randomly reconstructed from these threads. Hence the best way to utilise all of the information in the 10 high-SNR product-space runs is to unravel them into their constituent threads (5900 of them in this case) and rethread these into realisations of a run with =5900, which produces the value =-1.16±0.02. The error is far smaller than that of any run in Figure <ref>, which is perhaps unsurprising given the combination of information from all the runs. To demonstrate that this procedure does actually serve to decorrelate the individual threads, 10 new runs with =590 are constructed through rethreading; their results (the grey points in Figure <ref>) are seen to exhibit none of the sampling bias present in the original runs.



§ CONCLUSION



In this paper, we have adapted, combined and assessed a variety of recent modelling/statistical techniques to devise a preliminary framework for testing GR with EMRI observations from future space-based GW detectors. A generalised EMRI waveform model <cit.> and its toy surrogate are trialled in a null-hypothesis test developed for LIGO sources <cit.>; the method of product-space nested sampling <cit.> with rethreading error estimates <cit.> is shown to systematically increase computational efficiency by an order of magnitude over regular evidence-based sampling.

The results and observations presented here are quite general; they are relevant not just for the outlined EMRI test of GR, but indeed any similar parametrised test that uses a generalised waveform model to describe GW sources in modified gravity (although the need to reduce computational cost is most strongly motivated for EMRIs). Product-space nested sampling with rethreading is furthermore shown to be efficient, robust and hence potentially useful for a broad range of model selection problems beyond the null-hypothesis test in this work. While much of the present analysis is exclusive to nested-sampling theory, some results (e.g. the characterisation of posterior versus evidence errors) might also be applicable to other algorithms such as product-space MCMC.

Although the computational savings afforded by our proposed methods are promising, this work is only the first step in developing a practical infrastructure for testing GR with future EMRI observations. The framework should eventually incorporate other techniques for increased efficiency, such as reduced-order quadratures <cit.> to accelerate individual likelihood evaluations, or dynamic nested sampling <cit.> to improve sampling convergence. Finally, the actual accuracy and instructiveness of the null-hypothesis test must also be validated on data sets containing realistic source signals and detector noise, e.g. as performed by <cit.> for the constraints on deformation parameters in the LIGO test, but with additional focus on the final posterior odds ratio.



§ ACKNOWLEDGEMENTS


We thank Nicolas Yunes, Michele Vallisneri and Stephen Taylor for useful comments on the manuscript. AJKC acknowledges support from the Jet Propulsion Laboratory (JPL) Research and Technology Development programme. SH thanks the Science and Technology Facilities Council (STFC) for financial support. CJM acknowledges financial support provided under the European Union's H2020 ERC Consolidator Grant “Matter and strong-field gravity: New frontiers in Einstein's theory” grant agreement no. MaGRaTh646597, and networking support by the COST Action CA16104. Parts of this work were performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England and funding from STFC. Parts of this work were also undertaken on the COSMOS Shared Memory system at DAMTP, University of Cambridge operated on behalf of the STFC DiRAC HPC Facility; this equipment is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. Parts of this work were also carried out at JPL, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. ©2018. All rights reserved.

mnras


