# Testing General Relativity using Quasi-Periodic Oscillations from X-ray
Black Holes: XTE J1550-564 and GRO J1655-40

Katherine Rink1 Ilaria Caiazzo2, Jeremy Heyl1,
1Department of Physics and Astronomy, University of British Columbia,
Vancouver, BC V6T 1Z1, Canada
2TAPIR, Walter Burke Institute for Theoretical Physics, Mail Code 350-17,
Caltech, Pasadena, CA 91125, USA
Current address: UMass Dartmouth; email<EMAIL_ADDRESS><EMAIL_ADDRESS>Sherman Fairchild Fellowemail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We use the Relativistic Precession Model (RPM) (Stella et al., 1999a) and
quasi-periodic oscillation (QPO) observations from the Rossi X-ray Timing
Explorer to derive constraints on the properties of the black holes that power
these sources and to test General Relativity (GR) in the strong field regime.
We extend the techniques outlined by Motta et al. (2014a, b) to use pairs of
simultaneously measured QPOs, rather than triplets, and extend the underlying
spacetime metric to constrain potential deviations from the predictions of GR
for astrophysical black holes. To do this, we modify the RPM model to a Kerr-
Newman-deSitter spacetime and model changes in the radial, ecliptic, and
vertical frequencies. We compare our models with X-ray data of XTE J1550-564
and GRO J1655-40 using robust statistical techniques to constrain the
parameters of the black holes and the deviations from GR. For both sources we
constrain particular deviations from GR to be less than one part per thousand.

###### keywords:

black hole physics, gravitation, relativistic processes, accretion discs,
X-rays: XTE J1550-564, X-rays: GRO J1655-40

††pubyear: 2021††pagerange: Testing General Relativity using Quasi-Periodic
Oscillations from X-ray Black Holes: XTE J1550-564 and GRO J1655-40–References

##  1 Introduction

The X-ray flux from accreting black holes presents interesting nearly periodic
fluctuations, called quasi-periodic oscillations or QPOs, that are usually
found upon inspection of the power density spectrum (PDS) (Miyamoto &
Kitamoto, 1991; Takizawa et al., 1997; van der Klis, 2005, 2006). Such
fluctuations are detected as narrow peaks in the PDS, and often appear at
three distinct frequencies. When three frequencies are observed
simultaneously, they usually consist of two higher frequencies ($\sim$100 Hz)
and one lower frequency (ranging from a few mHz to $\sim$ 30 Hz) (Franchini et
al., 2016).

The true nature of QPOs is still unknown, but the frequencies detected are
similar to the orbital and epicyclic timescales of the inner accretion disk of
a Kerr black hole, hinting to the possible origin of these phenomena. Several
models have been proposed to explain QPOs, and they have mostly focused on
either the low-frequency or the high-frequency oscillations. For the low-
frequency QPOs, the most commonly invoked mechanisms are either some
instabilities in the accretion flow or geometric oscillations (Chakrabarti &
Molteni, 1993; Tagger & Pellat, 1999; Cabanac et al., 2010), and in particular
Lense-Thirring precession (Stella & Vietri, 1998; Stella et al., 1999b;
Wagoner et al., 2001; Schnittman et al., 2006; Ingram et al., 2009). The
latter consists of a nodal precession of orbits inclined to the equatorial
plane caused by a spinning compact object dragging the surrounding spacetime
around with it (i.e. the frame dragging effect). The high-frequency
oscillations have been observed more rarely in black hole systems, and are
more obscure; the proposed models include Doppler modulation of orbiting
hotspots in the inner disk, oscillation modes of a pressure-supported torus,
nonlinear resonances, gravity and pressure modes in the accretion disk (see
Lai et al., 2013; Belloni & Stella, 2014, and references therein).

One of the models that tries to explain the different, high and low,
frequencies together is the Relativistic Precession Model (RPM) (Stella et
al., 1999a). The theory argues that the coronal disk (the innermost part of
the accretion disk) is perturbed and tilted out of the equatorial plane. The
three QPOs would then appear in the PDS if a particularly hot region in the
coronal disk formed due to random perturbations in temperature. As the node of
the orbit precesses, the flux will vary slowly at the nodal precession
frequency. The fastest variation would then be at the Kepler frequency,
corresponding simply to a hot blob of material moving through its orbit.
Finally, the intermediate frequency would correspond to the periastron
precession. Given the geometry of elliptical orbits one would expect the
material to be compressed most at the periastron, so as this location
precesses about the black hole, the flux could vary at this frequency as well.
This would result in the observed object’s PDS showing peaks in these three
frequencies.

For the X-ray binaries XTE J1550-564 and GRO J1655-40, Motta et al. (2014a, b)
published a catalogue of the QPO frequencies observed in these systems. In the
GRO J1655-40 PDS, three frequencies are observed simultaneously. Under the
assumption that all three frequencies result from an oscillation at the same
radius in the accretion disk, Motta et al. (2014a) derived the mass and spin
of the black hole, and the radius of the oscillation. For XTE J1550-564, Motta
et al. (2014b) used a pair of simultaneously observed frequencies, along with
an estimate of the black hole mass from other observations, to derive
constraints on the spin of the black hole and the radius of the oscillation.
In both studies, they used the derived values to calculate the expected
frequencies for a range of radii and compared these values with other
observations for each source. Here we will take this idea to its natural
conclusion.

In both sources, Motta et al. (2014a, b) present about 30 simultaneous
measurements of QPOs where a high frequency and a low frequency signal were
observed simultaneously. We use all these measurements to fit for the mass and
spin of each black hole as well as other parameters that we expect to vanish
for astrophysical black holes. In our calculations, we assume that the
physical properties of the black holes do not change from measurement to
measurement, that the RPM picture is the correct explanation for the QPOs
observed in these systems, and, when two frequencies are observed
simultaneously, that they come from the same region of the accretion disk and
we can successfully match them to a particular RPM frequency. Furthermore, we
use a robust fitting technique that is insensitive to outliers while
equivalent to a $\chi^{2}$-test in the limit of normally distributed data. In
this way, we can test the validity of the model and identify observations for
which the model may not fit or the frequencies were matched incorrectly.

When one uses a triplet of coincident frequencies as in Motta et al. (2014a),
one can only determine three values: the mass and spin of the black hole and
the radius at which the oscillation occurs. In principle, one might not be
able to find values that are consistent with a triplet, so one could compare
the values predicted by the model with separate observations; both of these
provide tests of the RPM picture and also of the spacetime within the model.
However, the conclusions drawn rest on the single observation. Here we use the
entire suite of observations to constrain the parameters of the black hole,
constrain deviations from the expectations of GR, determine whether the
observations are well fit by the RPM, and if not, which observations do not
fit the RPM.

To look at deviations from GR, we employ a black hole spacetime calculated
within the context of GR that retains all of the symmetries of the Kerr
spacetime, so that the necessary frequencies are straightforward to calculate
in closed form, and yet in general contains additional curvature due to known
sources that we expect to be absent in astrophysical black holes;
specifically, a significant electromagnetic field or a uniform energy density
with an equation of state similar to that of the cosmological constant. One
can interpret our results to constrain the energy density around the black
hole with respect to empty space within GR or to constrain the underlying
gravitational theory to be close to GR.

Although there is a large literature developing constraints on gravity (Will,
2014) from experiments in laboratories (Adelberger et al., 2003), the Solar
System (Chiba et al., 2007; Allemandi & Ruggiero, 2007), and pulsars in
binaries (Stairs, 2003) through the parameterized post-Newtonian (PPN)
formalism (Nordtvedt, 1968; Will & Nordtvedt, 1972), these results do not
translate into understanding gravity in the strong-field limit near a black
hole. Developing a fully consistent alternative theory of relativistic gravity
to GR is extremely challenging. The alternatives sometimes fail to yield even
a consistent picture of spacetime, are often only defined for small
perturbations about flat spacetimes and rarely result in stable black holes.
One technique to make progress is to alter the black hole metric in a
parameterized fashion inspired by the PPN formalism (Ghasemi-Nodehi & Bambi,
2016).

More recent studies exploring GR alternatives for fitting observed QPO
triplets from GRO J1655-40 include using quadratic gravity theories (Maselli
et al. (2017)), the Einstein-Dilaton-Gauss-Bonnet theory (Maselli et al.
(2015)), the Bardeen metric and Johannsen-Psaltis metric (Bambi (2015)), and
quadrupolar metrics such as the Kerr-Q and Hartle-Thorne metrics (Allahyari &
Shao (2021)). We take the approach to work within GR but change the sources,
which in turn yields to potential deviations in the spacetime that we can hope
to measure.

The question of whether GR applies to black holes has very profound
implications. Most of our fundamental theories in physics have a range of
applicability and break down at some particular scale. For example, the
electroweak theory separates into electromagnetism and weak theory at energies
lower than $\sim$ 200 GeV, while strong and electroweak theories are thought
to unify at very high energies. Since most of the current tests of GR are
performed in the weak gravitational fields present in our solar system, the 6
to 7 orders of magnitude higher gravitational potential and $\sim$20 orders of
magnitude greater curvature around a black hole provide a more extreme
environment to test GR in the strong field regime. Much progress has been made
in this sense in the field of gravitational waves. However, the wavelength of
gravitational waves is necessarily comparable to the size of the colliding
objects, which therefore limits the scope for probing in detail the spacetime
surrounding compact objects; therefore, electromagnetic signals provide a
crucial complementary window into strong gravity.

One of the most subtle consequences of GR is the “no-hair” theorem, for which
black holes can be fully characterized by their mass, angular momentum, and
charge (Israel, 1967, 1968; Carter, 1968, 1971; Wald, 1971). Because we expect
no charge on astrophysical black holes, and we can neglect the mass in the
immediate vicinity of the black hole, the spacetime that surrounds a black
hole can be nearly exactly described by the Kerr metric. The only way to test
this theorem is to probe the spacetime very close to the hole. Fortunately,
the X-ray emission of accreting black holes carries information about the
inner region of the accretion disk, within a few gravitational radii
($R_{g}=GM/c^{2}$) from the hole, encoded in the fast variability of its
spectrum.

##  2 Geodesics in the Kerr-Newman-de-Sitter Metric

There exists four primary metrics for describing the theoretical types of
black holes. The Schwarzschild metric describes a non-rotating, uncharged
black hole; this is the simplest case and is not expected to be common. The
Kerr metric is the most widely used for it describes a rotating and uncharged
black hole. The Reissner-Nordström metric describes the least probable
scenario, a charged and non-rotating black hole. The metric that is modified
in this paper for the derivation of the frequency equations is the Kerr-Newman
metric. This metric describes a charged and rotating black hole.

Black holes can theoretically carry charge, but an accreting astrophysical
black hole is immersed in plasma and would discharge quickly, so we expect the
observed black holes to have negligible electromagnetic charge. A measurement
of the effects of charge on the spacetime surrounding a black hole constrains
deviations from GR or, in principle, the presence of additional charged
components in the Universe. Here we add an extra parameter to the Kerr-Newman
metric that results from an uniform energy density near the black hole with an
equation of state similar to that of the cosmological constant. These
constituents yield the Kerr-Newman-deSitter Metric (Carter, 1973, 2009):

$\begin{split}ds^{2}=\frac{\Delta_{r}^{2}}{\rho^{2}}&\left(dt-\frac{a}{\Xi}\sin^{2}(\theta)d\phi\right)^{2}-\frac{\rho^{2}}{\Delta_{r}^{2}}dr^{2}-\frac{\rho^{2}}{\Delta_{\theta}}d\theta^{2}-\\\
&~{}~{}~{}\frac{\Delta_{\theta}\sin^{2}(\theta)}{\rho^{2}}\left(adt-\frac{r^{2}+a^{2}}{\Xi}d\phi\right)^{2}\end{split}$
(1)

where

$\displaystyle\Delta_{r}^{2}$ $\displaystyle=$
$\displaystyle(r^{2}+a^{2})(1-H^{2}r^{2})-2Mr+Q^{2}$ (2)
$\displaystyle\Delta_{\theta}$ $\displaystyle=$ $\displaystyle
1+H^{2}a^{2}\cos^{2}(\theta)$ (3) $\displaystyle\Xi$ $\displaystyle=$
$\displaystyle 1+H^{2}a^{2}$ (4) $\displaystyle\rho^{2}$ $\displaystyle=$
$\displaystyle r^{2}+a^{2}\cos^{2}(\theta)$ (5)

In the above expressions, $t$ is a timelike coordinate, $r$ is a radial
coordinate, $a$ is the spin of the black hole, $\theta$ is the angle above or
below the plane of the accretion disk, $\phi$ is the angle around the black
hole within the plane of the accretion disk, $Q$ is the charge of the black
hole, and finally $H^{2}$ is the cosmological constant parameter. For
simplicity, the $t-$coordinate in Eq. 1 corresponds to the proper time
measured by a distant observer. We will use the value of $Q^{2}$ throughout
the rest of the paper because we interpret $Q^{2}$ as a parameter to change
the structure of the spacetime rather than to quantify the charge of the black
hole.

Furthermore, if the value of $H^{2}$ does not vanish, the spacetime as given
is not asymptotically flat; therefore, Eq. 1 should be interpreted as a chart
that is valid near the black hole in the region where the QPOs originate,
typically about ten gravitational radii. Beyond this region, Eq. 1 is to be
considered attached to another chart that is asymptotically flat.

In the RPM picture (Stella et al., 1999a), the observed frequencies of the
QPOs correspond to three different frequencies in the orbital motion, believed
to be the Kepler ($\Omega_{\phi}$), the nodal precession ($\Omega_{\rm
node}$), and the periastron precession ($\Omega_{\rm peri}$) frequencies. The
Kepler frequency represents the angular frequency at which the object
completes one full orbit about the spin axis of the black hole from the point
of view of a distance observer. It is defined as the value of $d\phi/dt$ for a
geodesic with a constant value of $r$ and $\theta=\pi/2$ (the equatorial
plane). The nodal precession frequency is the difference between the Kepler
frequency and the vertical epicyclic frequency ($\omega_{\theta}$). The latter
is the frequency of oscillations in the $\theta$-direction about the orbits
with constant $r$ in the equatorial plane. This can be thought of as the
precession of the two points at which the tilted orbit intersects the
equatorial plane. The periastron precession frequency is the difference
between the Kepler frequency and the radial epicyclic frequency ($\omega_{r}$)
and represents the precession of the point of closest approach, where we’d
expect to see peak emission.

The three observed frequencies can be represented as follows

$\displaystyle\Omega_{\phi}$ $\displaystyle=$
$\displaystyle\frac{u^{\phi}}{u^{t}}$ (6)
$\displaystyle\Omega_{\textrm{\scriptsize node}}$ $\displaystyle=$
$\displaystyle\Omega_{\phi}-\omega_{\theta}$ (7)
$\displaystyle\Omega_{\textrm{\scriptsize peri}}$ $\displaystyle=$
$\displaystyle\Omega_{\phi}-\omega_{r}\,.$ (8)

To derive the three measured QPO frequency equations for this model, their
corresponding geodesic equations are needed to find the necessary 4-velocity
vector and its perturbations. The geodesic equation is used to derive
$\Omega_{\phi}$, $\omega_{r}$, and $\omega_{\theta}$.

We begin the derivations for a general stationary, axisymmetric spacetime (see
Ghasemi-Nodehi et al., 2020, for a similar treatment) by finding the Kepler
frequency, $\Omega_{\phi}$, using the geodesic equation with respect to
radius, $r$. Since $u_{r}=g_{rr}u^{r}$, and $u^{r}=\frac{dr}{d\tau}$, we can
re-write the geodesic equation using:

$\frac{du_{r}}{d\tau}=\frac{d}{d\tau}g_{rr}u^{r}=\frac{d}{d\tau}\left(g_{rr}\frac{dr}{d\tau}\right)\,.$
(9)

However, since the object’s radial velocity is constant over time, we can
write the change in this velocity as
$\frac{d}{d\tau}\left(g_{rr}\frac{dr}{d\tau}\right)=0$, simplifying the
geodesic equation to:

$\frac{1}{2}(g_{\mu\nu,r}u^{\mu}u^{\nu})=0\,.$ (10)

For the Kepler frequency, we assume that the orbit is at a constant radius in
the equatorial plane, so the sum reduces to three terms

$\frac{1}{2}(g_{tt,r}u^{t}u^{t}+2g_{t\phi,r}u^{t}u^{\phi}+g_{\phi\phi,r}u^{\phi}u^{\phi})=0\,.$
(11)

Because the Kepler frequency is $\Omega_{\phi}=u^{\phi}/u^{t}$, dividing the
equation by $u^{t}u^{t}/2$ yields a quadratic equation for the Kepler
frequency

$(g_{tt,r}+2g_{t\phi,r}\Omega_{\phi}+g_{\phi\phi,r}\Omega_{\phi}^{2})=0\,.$
(12)

Using the coefficients from the metric and solving the quadratic equation for
the Kepler frequency gives

$\Omega_{\phi}=\frac{(Mr-H^{2}r^{4}-Q^{2})(1+H^{2}a^{2})}{r^{2}\sqrt{Mr-H^{2}r^{4}-Q^{2}}-a(Mr-H^{2}r^{4}-Q^{2})}\,.$
(13)

The complete expression for the Kepler frequency shows where in the spacetime
each of three parameters $M$, $Q$ and $H$ play a role. If we take $Q=H=a=0$,
we recover Kepler’s third law

$\Omega_{\phi}^{2}=\frac{M}{r^{3}}\,,$ (14)

allowing only $Q^{2}$ to have a finite value while the others vanish yields

$\Omega_{\phi}^{2}=-\frac{Q^{2}}{r^{4}}\,,$ (15)

and finally, allowing only $H^{2}$ to have a finite value yields

$\Omega_{\phi}^{2}=-H^{2}.$ (16)

Thus, relative to the mass of the black hole, the charge affects the frequency
at smaller radii and $H$ affects the frequency at larger radii.

For the remaining two frequencies, the value of $u^{t}$ is needed for the
circular orbit. We begin by noting that

$u\cdot u=u_{\alpha}u^{\beta}=g_{\alpha\beta}u^{\alpha}u^{\beta}=1\,.$ (17)

Because the nodal and periastron precession frequencies can be thought of as
perturbations to the Kepler frequency, a circular orbit is assumed in the
derivation of $u^{t}$. This assumption results in only the $t$ and $\phi$
components playing a role, such that

$1=g_{tt}(u^{t})^{2}+2g_{t\phi}u^{\phi}u^{t}+g_{\phi\phi}(u^{\phi})^{2}\,.$
(18)

Remembering that $\Omega_{\phi}=u^{\phi}/u^{t}$, and therefore
$u^{\phi}=\Omega_{\phi}u^{t}$:

$1=g_{tt}(u^{t})^{2}+2g_{t\phi}\Omega_{\phi}(u^{t})^{2}+g_{\phi\phi}(\Omega_{\phi}u^{t})^{2}\,,$
(19)

where the Kepler frequency, $\Omega_{\phi}$, was previously derived. This
leaves the above equation with only one unknown, $u^{t}$, and hence:

$u^{t}=(g_{tt}+2g_{t\phi}\Omega_{\phi}+g_{\phi\phi}\Omega_{\phi}^{2})^{-\frac{1}{2}}\,.$
(20)

To find the radial epicyclic frequency, $\omega_{r}$, we use the conserved
quantities for the orbit. The quantity $u_{\phi}$ is the angular momentum
($L$) and $u_{t}$ is the energy ($E$). With these new definitions, and a
corresponding change of indices from lowered to raised on all metric
coefficients, we can rewrite Eq. 18 as follows:

$1=g^{tt}E^{2}+2g^{t\phi}EL+g^{\phi\phi}L^{2}\,.$ (21)

To perturb the orbit, additional terms are added to account for movement in
the radial direction.:

$1+E_{2}=g^{tt}E^{2}+2g^{t\phi}EL+g^{\phi\phi}L^{2}+g_{rr}(u^{r})^{2}\,,$ (22)

where $u^{r}$ is the radial component of the four velocity, $E_{2}$ is the
additional energy of the perturbation and the values of $E$ and $L$ are held
fixed at the values for a circular orbit. We perturb the position about the
circular orbit which has $r=r_{0}$ and use Eq. 18 and the geodesic equation of
the circular orbit (Eq. 12) to yield the second-order perturbation

$E_{2}=\frac{\left(r-r_{0}\right)^{2}}{2}\left(g_{,rr}^{tt}E^{2}+2g_{,rr}^{t\phi}EL+g_{,rr}^{\phi\phi}L^{2}\right)+g_{rr}(u^{r})^{2}\,.$
(23)

By dividing the above equation by two, it becomes clear that the
$\frac{1}{2}g_{rr}(u^{r})^{2}$ term behaves as kinetic energy (classically of
the form $\frac{1}{2}mv^{2}$). The above equation is therefore a conservation
equation, with the remaining terms accounting for the effective potential
energy,

$E_{2}=\frac{1}{2}k\left(r-r_{0}\right)^{2}+\frac{1}{2}mv^{2}$ (24)

where $v=u^{r}=dr/d\tau$, $m=g_{rr}$ and

$k=\frac{1}{2}\left(g_{,rr}^{tt}E^{2}+2g_{,rr}^{t\phi}EL+g_{,rr}^{\phi\phi}L^{2}\right).$
(25)

In the frame of the orbiting particle, the oscillation frequency is given by

$\omega^{2}=\frac{k}{m}=\frac{g_{,rr}^{tt}E^{2}+2g_{,rr}^{t\phi}EL+g_{,rr}^{\phi\phi}L^{2}}{2g_{rr}}\,.$
(26)

The next step is to convert everything from proper time to time in the
observer’s reference frame. After converting reference frames, we are left
with our radial solution:

$\omega_{r}^{2}=\frac{g_{,rr}^{tt}E^{2}+2g_{,rr}^{t\phi}EL+g_{,rr}^{\phi\phi}L^{2}}{2g_{rr}(u^{t})^{2}}\,.$
(27)

This equation can be simplified further by recognizing that
$u^{\phi}=\Omega_{\phi}u^{t}$ and dividing out the $(u^{t})^{2}$ term in the
denominator

$\omega_{r}^{2}=\frac{g_{tt,rr}+2g_{t\phi,rr}\Omega_{\phi}+g_{\phi\phi,rr}\Omega_{\phi}^{2}}{2g_{rr}}\,.$
(28)

Because no assumptions were made in the derivation of $\omega_{r}$ that would
not also apply to $\omega_{\theta}$, the solution for the final frequency
$\omega_{\theta}$ can be carried out in the exact same manner but holding $r$
constant and perturbing $\theta$. For brevity, only the final result is
included here:

$\omega_{\theta}^{2}=\frac{g_{tt,\theta\theta}+2g_{t\phi,\theta\theta}\Omega_{\phi}+g_{\phi\phi,\theta\theta}\Omega_{\phi}^{2}}{2g_{\theta\theta}}\,.$
(29)

##  3 Method

It is perhaps most straightforward to explain our technique and the nature of
the dataset by looking at the frequencies for the source XTE J1550-564 (Motta
et al., 2014b). Fig. 1 depicts for this source the observed frequency pairs,
that is, a low-frequency and a high-frequency QPO observed simultaneously. The
low-frequency ones, which are interpreted to be the nodal precession
frequency, are depicted along the horizontal axis and the high-frequency ones,
which are interpreted as either the Kepler (larger) or periastron-precession
(smaller) frequency, along the vertical. Motta et al. (2014b) also observed
many single frequencies which they interpreted as the low-frequency component,
not depicted here. They fit the data points, using a prior on the black hole
mass of 9.1 $\textrm{M}_{\odot}$, through a five-step process that includes
solving for the radius of each QPO given a model of fixed mass and spin and
then finding the ensemble of spins consistent with the mass constraint and the
observational data. For Fig. 1 we calculate the three angular QPO frequencies
as a function of radius for a particular black hole mass and spin. We recast
these calculations as the periastron and Kepler frequencies as a function of
the nodal precession frequency and fit the observed pairs of frequencies to
the corresponding functions, converting from angular to spatial frequencies.
Unlike Motta et al. (2014b), we do not use the values of the low-frequency QPO
if a high-frequency QPO is not measured simultaneously.

Figure 1: The observed low-frequency QPOs from XTE J1550-564 interpreted to be
nodal precession frequencies and the simultaneously observed high-frequency
QPOS from Motta et al. (2014b). The dashed curves give the values predicted by
the RPM for the parameters determined by Motta et al. (2014b),
$M=9.1~{}\textrm{M}_{\odot}$ and $a/M=0.34$ ($\chi^{2}=1456$) . The solid
curves trace the best-fitting model (the smallest value of $\chi^{2}$, 386)
with $M=11.5~{}\textrm{M}_{\odot}$ and $a/M=0.49$. The Kepler frequency and
the periastron precession frequency ($\nu_{\phi}$ and
$\nu_{\textrm{\scriptsize peri}}$) are equal to each other at the innermost
stable orbit.

An examination of Fig. 1 shows that the values of mass and spin determined by
Motta et al. (2014b) do not fit all of the data well. In fact the value of
$\chi^{2}$ is 1456 with 31 degrees of freedom. Motta et al. (2014b) use other
observations to fix the mass of the black hole. This choice forces the values
of the Kepler frequency within the model to lie along the upper envelope of
the observations. On the other hand, the periastron precession frequency does
trace the bulk of the observations. The solid curves depict the best-fitting
model with a $\chi^{2}$ value of 386. This is a better but still poor fit to
the observations. Although this model traces the bulk of the highest frequency
QPOs, it misses the bulk of the middle frequency QPOs. In this case, it is the
QPO observations themselves that create the tension.

If any of the assumptions underlying the RPM fail, the observed QPOs need not
fit the model. In particular, if some oscillations correspond to frequencies
other than the nodal, periastron precession and Kepler frequencies, or if, on
occasion, there are perturbations at several radii, then outlying observations
will result. The underlying assumption behind $\chi^{2}$-fitting of a dataset
is that the residuals between the data and the model are normally distributed.
Where the RPM fails for one reason on another, this is not the case.

For this reason, we seek a technique that is robust against the presence of a
few outliers in our dataset and yet is efficient in determining the best-
fitting model with the fewest observations. To motivate our figure of merit to
fit the observed QPO frequencies, we start with the traditional $\chi^{2}$
estimator,

$\chi^{2}=\sum_{i}\left(X_{i}\right)^{2}=\left(\sum_{j}X_{j}\right)^{2}+\sum_{i}\left(X_{i}-\frac{1}{n}\sum_{j}X_{j}\right)^{2}$
(30)

where $n$ is the number of measurements and

$X_{i}=\frac{y_{i}-f(x_{i}|\textrm{model parameters})}{\sigma_{i}}$ (31)

is the residual of a particular measurement $y_{i}$ with respect to the model,
inversely weighted by the uncertainty in that measurement. The first term in
Eq. 30 is the product of the mean value of $X_{i}$ and the number of
measurements, squared. The second term is the product of the variance in the
value $X_{i}$ and the number of measurements, also squared. To derive a robust
figure of merit, we can replace the mean and variance with robust estimators
of location and scale. Two such estimators that may be familiar are the median
and the median absolution deviation (MAD),

$\textrm{MAD}=\textrm{median}\left\\{\left|X_{i}-\textrm{median}\left(X_{i}\right)\right|\right\\}.$
(32)

Although these two estimators are robust with a breakdown point of fifty
percent, they are not efficient, which means that, in the case of normally
distributed measurements, one would need many more measurements to obtain the
same constraints on the models as one would achieve with the traditional
$\chi^{2}$ estimator. We use two alternative statistics based on pairs of
measurements: the Hodges-Lehmann-Sen estimator $H\\!L$ (Hodges & Lehmann,
1963; Sen, 1963) and $Q_{n}$ (Rousseeuw & Croux, 1993), which are defined as

$\displaystyle H\\!L$ $\displaystyle=$
$\displaystyle\textrm{median}\left\\{\frac{X_{i}+X_{j}}{2}:i\leq j\right\\},$
(33) $\displaystyle Q_{n}$ $\displaystyle=$ $\displaystyle
c_{n}\textrm{~{}first quartile
of~{}}\left\\{\left|X_{i}-X_{j}\right|:i<j\right\\},$ (34)

where the value of $c_{n}$ depends on the size of the sample and scales the
value of $Q_{n}$ to be the standard deviation in the case of normally
distributed data. Although both of these estimators appear to require sorting
the complete list of pairs and ${\cal O}(n^{2}\log n)$ operations, there are
efficient techniques111A python library is available at
https://github.com/UBC-Astrophysics/qn_stat., ${\cal O}(n\log n)$, to
determine both (Rousseeuw & Croux, 1993; Monahan, 1984).

Because both of these estimators rely on ranks like the median and MAD, they
are robust against outliers with breakdown points of 29 and 50 percent,
respectively. A 50 percent breakdown point means that one could change the
values of half of the data without changing the estimator. However, they are
both efficient. The expected value of the standard deviation of the mean is
97% of that $H\\!L$ for normally distributed data. The corresponding value for
the median is 80%. The comparison for $Q_{n}$ and MAD are 88% and 62%. We
propose

$R^{2}=n\left[H\\!L\left(X_{i}\right)\right]^{2}+n\left[Q_{n}\left(X_{i}\right)\right]^{2}$
(35)

to fit the observed QPO frequencies. The $R^{2}$ estimator is an unbiased
estimator for $\chi^{2}$ for normally distributed data. Furthermore, the
$H\\!L$ and $Q_{n}$ estimators trade a modest amount of efficiency for
robustness. They require slightly larger samples (about 40% larger) to achieve
the same precision on model parameters as the traditional $\chi^{2}$ estimator
for normally distributed data, but are much less sensitive to outliers.

Fig. 2 shows the result of the fitting process. The most obvious difference is
in the model values for $\nu_{\textrm{\scriptsize peri}}$ that now avoid the
two lower points that were apparent outliers in Fig. 1. This model has a lower
spin of $a/M=0.31$, similar to the value of Motta et al. (2014b) but a larger
mass of $M=19.0~{}\textrm{M}_{\odot}$. The value of $R^{2}$ is 57 compared to
$\chi^{2}$ of 386 for the best fit in Fig. 1. Using the values derived by
Motta et al. (2014b) yields a value of $R^{2}$ of 185.

Here we note a key qualitative difference between our fit and that of Motta et
al. (2014b). The value of the mass is significantly larger in our fit which
means that the maximum value of the nodal precession frequency is less than 10
Hz. Motta et al. (2014b) presented measurements of single QPOs with
frequencies up to nearly 20 Hz which they interpreted as the nodal precession
frequency. This interpretation requires a smaller mass than the one we
obtained here. The value of $R^{2}$ is similar to the number of degrees of
freedom, 31. This result is consistent with the fact that the bulk of the
measurements are well characterized by the RPM and that the uncertainties have
been well estimated. We achieve this without explicitly identifying outliers.

Figure 2: Same as Fig.1 but using the robust $R^{2}$ fitting technique. The
solid curves trace the best-fitting model using the robust technique with
$M=19.0~{}\textrm{M}_{\odot}$ and $a/M=0.31$. The value of $R^{2}$ is 57
compared to $\chi^{2}$ of 386 for the best fit in Fig. 1.

A crucial final part of the analysis is the estimation of the uncertainties in
the results. We estimate these uncertainties by a bootstrap resampling of the
dataset, in which we select the same number of measurements as observed with
replacement, so that the same measurement could be included many times or not
at all. Fig. 3 and 4 depict results of the bootstrap resample for the
$\chi^{2}$ and $R^{2}$ analysis. The best-fitting values for the original
dataset are depicted in blue for both techniques. In particular, the
$\chi^{2}$ analysis yields a smaller mass and larger spin than the $R^{2}$
analysis. However, the best-fit $\chi^{2}$ model lies within the 99%
confidence region of the $R^{2}$ analysis and vice versa. The uncertainties on
the mass from the $R^{2}$ analysis are more than twice as large as for the
$\chi^{2}$ analysis, but the uncertainties on the spin are comparable.
Although the constraints are weaker with the $R^{2}$ analysis, we choose to
use it to constrain and quantify potential deviations from GR so that we can
obtain the most conservative estimates of these quantities.

Figure 3: The results of the bootstrap resampling and $\chi^{2}$ analysis of
the data for XTE J1550-564. The ranges depicted each contain 99% of the values
that results from the fits. The blue lines indicate the best-fitting
parameters for the original sample. Figure 4: The results of the bootstrap
resampling and robust $R^{2}$ analysis of the data for XTE J1550-564. The
ranges depicted each contain 99% of the values that results from the fits. The
blue lines indicate the best-fitting parameters for the original sample.

##  4 Results

Our focus so far has been to demonstrate and explain our technique for
analyzing the QPO data with XTE J1550-564 as an example and by focusing on the
standard GR parameters of mass and spin. However, in addition to XTE
J1550-564, Motta et al. (2014a) published an analysis of RXTE observations of
GRO J1655-40. For this object, they found a simultaneous triplet of QPOs from
which they could determine the mass and spin of the black hole (and the radius
of the oscillations) from the RPM equations, yielding
$M=5.3~{}\textrm{M}_{\odot}$ and $a/M=0.29$.

Fig. 5 depicts the best-fitting model for GRO J1655-40 using the $R^{2}$
estimator compared with the model of Motta et al. (2014a). Again the robust
$R^{2}$ estimator finds a larger mass resulting in smaller frequencies at the
innermost stable orbit. Motta et al. (2014a) present observations of single
QPOs up to 25 Hz which would be difficult to reconcile with the best-fitting
model if these single QPOs are manifestations of the nodal precession
frequency. Turning to the constraining power of the robust fitting technique,
we find in Fig. 6 that the distribution of bootstrapped samples is centered
about $6.7~{}\textrm{M}_{\odot}$ and the values determined by Motta et al.
(2014a) lie within the range of the fits.

Figure 5: Same as Fig.1 but using the robust $R^{2}$ fitting technique for GRO
J1655-40. The dashed curves give the values predicted by the RPM for
parameters determined by Motta et al. (2014a) $M=5.3~{}\textrm{M}_{\odot}$ and
$a/M=0.29$. The solid curves trace the best-fitting model using the robust
technique with $M=9.1~{}\textrm{M}_{\odot}$ and $a/M=0.27$. Figure 6: The
results of the bootstrap resampling and robust $R^{2}$ analysis of the data
for GRO J1655-40. The ranges depicted each contain 99% of the values that
results from the fits. The blue lines indicate the best-fitting parameters for
the original sample.

We have been conservative in restricting the data to QPO doublets because the
identification of the QPO frequencies as RPM frequencies is most certain in
those cases and our ultimate goal is to derive constraints on GR, not the
properties of the black hole for which the usefulness of the singlet QPOs
trumps the uncertainties in their identification. To place constraints on the
additional metric parameters, we fit the mass, spin, charge and cosmological
constant to the data and to the bootstrapped samples to obtain confidence
regions on all of the parameters simultaneously, as depicted in Fig. 7. In
this broader class of models, the masses and spins are again typically larger
than those obtained by Motta et al. (2014a) and Motta et al. (2014b), but the
masses are smaller than obtained in the restricted GR fits.

The key results here are the constraints on the parameters $Q^{2}$ and $H^{2}$
depicted in the figure in natural units of the black hole mass. The value of
$Q^{2}/M^{2}$ is constrained at the $10^{-3}$ level, and $H^{2}M^{2}$ at the
$10^{-6}$ level. To understand the significance of these constraints on the
components of the metric, we have to include the typical radius where the QPOs
are excited within the model. Using the observed Kepler frequencies and the
best-fitting masses, we find that the QPOs are typically excited around a
radius of $10~{}GM/c^{2}$. Comparing either the terms in the metric itself
(Eq. 1) or the square of the angular frequencies (Eq. 14 to 16) yields the
following fractional constraints on the deviation of the metric coefficients
or frequencies from the GR values:

$\left|\frac{Q^{2}}{Mr}\right|\lesssim
10^{-4}~{}\textrm{and}~{}\left|\frac{H^{2}r^{3}}{M}\right|\lesssim 10^{-3}$
(36)

where we have used $r=10M$. These constraints are the key result of this work.

GRO J1655-40

XTE J1550-564

Figure 7: The results of the bootstrap resampling and robust $R^{2}$ analysis
of the data for XTE J1550-564 (left) and GRO J1655-40. The ranges depicted for
$M$ and $a/M$ each contain 99% of the values that results from the fits. The
ranges for $Q^{2}/M^{2}$ and $H^{2}M^{2}$ contain 95% and 90% of the values
respectively. The blue lines indicate the best-fitting parameters for the
original sample.

##  5 Summary and Conclusions

Three types of QPOs have been observed coming from the accretion disks around
black holes. Though the main mechanism behind these frequencies is still
unknown, the Relativistic Precession Model (Stella et al., 1999a) appears to
match the data well when tested. This model was therefore used to test GR in
the strong field regime using the Kerr-Newman-deSitter metric, which describes
a spinning and charged black hole with the addition of a cosmological constant
like parameter.

Thus far, the cosmological constant has been assumed to be negligible close to
the black-hole horizon. However, if we are to find a theory that connects
general relativity and quantum mechanics, it is essential that we understand
each theory to their fullest extent and make no approximations in the coming
age of enhanced instrumentation. We therefore used a nonzero cosmological
constant parameter to constrain any type of excess energy density that could
be in the vicinity of the black hole.

After calculating three frequency equations to describe the three observed
QPOs in a Kerr-Newman-deSitter spacetime, our models were tested against RXTE
data. We compared our fit using four black hole parameters to a fit using the
standard Kerr metric as well as the fits derived in Motta et al. (2014a).
Their fits relied on a single instance of simultaneously observed QPOs. This
study fit for black hole parameters using models optimized to all observed
QPOs, instead of using calculated values based on simultaneous QPO
observations or values taken from the literature to test RPM.

Finally, statistical tests were conducted to compare the three models more
precisely. We were able to make fractional constraints on deviations from GR
at a radius of $10M$ to be

$\left|\frac{Q^{2}}{Mr}\right|\lesssim
10^{-4}~{}\textrm{and}~{}\left|\frac{H^{2}r^{3}}{M}\right|\lesssim 10^{-3}.$

Interpreting the results of our bootstrap resampling analysis, we see that the
charge parameter heavily peaked around $Q^{2}=0$ for both sources, as expected
(Fig.7). However, our cosmological constant parameter $H^{2}$ was rarely zero,
which suggests that there could be excess energy density or some similarly
behaving parameter yet unaccounted for in GR. It is also possible the nonzero
nature of the $H^{2}$ parameter fit is a result of one or more QPOs being
incorrectly classified.

We conducted our analysis using robust methods for testing GR, prioritizing
reduction in QPO identification uncertainty over the inclusion of more data
and using fitting techniques insensitive to outliers. Our methods were
therefore less reliable in recovering values for mass and spin that closely
agreed with those derived from X-ray observations. This does not diminish the
accuracy of our GR constraints as our best fit values fall within the
estimated uncertainties of past studies that relied on a priori assumptions
not considered here (Orosz et al., 2011; Beer & Podsiadlowski, 2002).

This study demonstrates that robust tests of GR are possible in the strong
field regime using QPO observations. Our methods, when applied to higher
resolution data, will help confirm the RPM theory and further constrain the
spacetime around black hole regimes. More observations of coincident QPOs with
precise timing would also help the accuracy of our $R^{2}$ analysis to further
constrain our fits.

Building on the growing literature of black hole X-ray timing analysis
techniques, future studies will be able to compare observational results
across models and further our understanding of the spacetime around black
holes. These proof-of-concept studies are essential in attaining this goal if
we wish to probe unknown parameter space and conduct precision tests of GR.

## Acknowledgements

This work has been supported by the Natural Sciences and Engineering Research
Council of Canada through the Discovery Grants program and Compute Canada.

Data Availability

The measurements used in this paper can be found in Tab. 1 of Motta et al.
(2014a) and Tab. 3 of Motta et al. (2014b).

## References

  * Adelberger et al. (2003) Adelberger E. G., Heckel B. R., Nelson A. E., 2003, Annual Review of Nuclear and Particle Science, 53, 77
  * Allahyari & Shao (2021) Allahyari A., Shao L., 2021, Testing No-Hair Theorem by Quasi-Periodic Oscillations: the quadrupole of GRO J1655$-$40 (arXiv:2102.02232)
  * Allemandi & Ruggiero (2007) Allemandi G., Ruggiero M. L., 2007, General Relativity and Gravitation, 39, 1381
  * Bambi (2015) Bambi C., 2015, The European Physical Journal C, 75
  * Beer & Podsiadlowski (2002) Beer M. E., Podsiadlowski P., 2002, MNRAS, 331, 351
  * Belloni & Stella (2014) Belloni T. M., Stella L., 2014, Space Sci. Rev., 183, 43
  * Cabanac et al. (2010) Cabanac C., Henri G., Petrucci P.-O., Malzac J., Ferreira J., Belloni T. M., 2010, MNRAS, 404, 738
  * Carter (1968) Carter B., 1968, Physical Review, 174, 1559
  * Carter (1971) Carter B., 1971, Phys. Rev. Lett., 26, 331
  * Carter (1973) Carter B., 1973, in Black Holes (Les Astres Occlus). pp 57–214
  * Carter (2009) Carter B., 2009, General Relativity and Gravitation, 41, 2873
  * Chakrabarti & Molteni (1993) Chakrabarti S. K., Molteni D., 1993, ApJ, 417, 671
  * Chiba et al. (2007) Chiba T., Smith T. L., Erickcek A. L., 2007, Phys. Rev. D, 75, 124014
  * Franchini et al. (2016) Franchini A., Motta S. E., Lodato G., 2016, Monthly Notices of the Royal Astronomical Society
  * Ghasemi-Nodehi & Bambi (2016) Ghasemi-Nodehi M., Bambi C., 2016, European Physical Journal C, 76, 290
  * Ghasemi-Nodehi et al. (2020) Ghasemi-Nodehi M., Lu Y., Chen J., Yang C., 2020, European Physical Journal C, 80, 504
  * Hodges & Lehmann (1963) Hodges J. L., Lehmann E. L., 1963, Ann. Math. Statist., 34, 598
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Israel (1967) Israel W., 1967, Physical Review, 164, 1776
  * Israel (1968) Israel W., 1968, Communications in Mathematical Physics, 8, 245
  * Lai et al. (2013) Lai D., et al., 2013, in Feeding Compact Objects: Accretion on All Scales. pp 57–61 (arXiv:1212.5323), doi:10.1017/S1743921312019205
  * Maselli et al. (2015) Maselli A., Gualtieri L., Pani P., Stella L., Ferrari V., 2015, The Astrophysical Journal, 801, 115
  * Maselli et al. (2017) Maselli A., Pani P., Cotesta R., Gualtieri L., Ferrari V., Stella L., 2017, The Astrophysical Journal, 843, 25
  * Miyamoto & Kitamoto (1991) Miyamoto S., Kitamoto S., 1991, ApJ, 374, 741
  * Monahan (1984) Monahan J. F., 1984, ACM Trans. Math. Softw., 10, 265–270
  * Motta et al. (2014b) Motta S. E., et al., 2014b, Monthly Notices of the Royal Astronomical Society
  * Motta et al. (2014a) Motta S. E., et al., 2014a, Monthly Notices of the Royal Astronomical Society
  * Nordtvedt (1968) Nordtvedt K., 1968, Physical Review, 169, 1017
  * Orosz et al. (2011) Orosz J. A., Steiner J. F., McClintock J. E., Torres M. A. P., Remillard R. A., Bailyn C. D., Miller J. M., 2011, ApJ, 730, 75
  * Rousseeuw & Croux (1993) Rousseeuw P. J., Croux C., 1993, Journal of the American Statistical Association, 88, 1273
  * Schnittman et al. (2006) Schnittman J. D., Homan J., Miller J. M., 2006, ApJ, 642, 420
  * Sen (1963) Sen P. K., 1963, Biometrics, 19, 532
  * Stairs (2003) Stairs I. H., 2003, Living Reviews in Relativity, 6, 5
  * Stella & Vietri (1998) Stella L., Vietri M., 1998, ApJ, 492, L59
  * Stella et al. (1999a) Stella L., Vietri M., Morsink S. M., 1999a, The Astrophysical Journal Letters
  * Stella et al. (1999b) Stella L., Vietri M., Morsink S. M., 1999b, ApJ, 524, L63
  * Tagger & Pellat (1999) Tagger M., Pellat R., 1999, A&A, 349, 1003
  * Takizawa et al. (1997) Takizawa M., et al., 1997, ApJ, 489, 272
  * Wagoner et al. (2001) Wagoner R. V., Silbergleit A. S., Ortega-Rodríguez M., 2001, ApJ, 559, L25
  * Wald (1971) Wald R. M., 1971, Phys. Rev. Lett., 26, 1653
  * Will (2014) Will C. M., 2014, Living Reviews in Relativity, 17, 4
  * Will & Nordtvedt (1972) Will C. M., Nordtvedt Kenneth J., 1972, ApJ, 177, 757
  * van der Klis (2005) van der Klis M., 2005, Ap&SS, 300, 149
  * van der Klis (2006) van der Klis M., 2006, Rapid X-ray Variability. Cambridge University Press, pp 39–112