Geodesic Models of Quasi-periodic-oscillations as Probes of Quadratic Gravity

The study of gravity near compact objects is among the last missing pieces of the grand program aimed at testing general relativity (GR) at all sub-galactic scales (for a recent review, see Berti et al. 2015). Because of their simplicity, black holes (BHs) are particularly well suited for testing gravity in the strong-field regime, which characterizes the dynamics near the horizon (for a recent review of BH-based tests of gravity, see Yagi & Stein 2016). X-rays emitted by matter accreting into stellar mass BHs provide a very promising probe of the inner region of the accretion disk, which is believed to be bounded by the Innermost Stable Circular Orbit (ISCO) of the BH. In this region, the gravitational field cannot be described by Newtonian theory or by a weak-field expansion of GR: the strong-field regime of gravity is manifest there and can be tested there.

Some X-ray spectroscopy features have been used as diagnostic tools of the inner disk region of BH accretion. These comprise the soft X-ray continuum emission, from which estimates of the ISCO location, and thus the BH rotation rate have been obtained (McClintock et al. 2011) and the broad iron Kα line and reflection spectrum, whose shape carries information on a variety of GR effects in the inner disk region (Fabian 2013). Among timing diagnostics, the multiple Quasi-periodic Oscillations (QPOs), which occur simultaneously in the X-ray flux from accreting stellar mass BHs and neutron stars (van der Klis 2006) are especially promising.⁷ Most QPO models, including the two models we adopt here, i.e., the relativistic precession model (RPM; Stella et al. 1999) and the epicyclic resonance model (ERM; Török et al. 2005), involve frequencies associated to the orbital motion of matter in the inner disk, which is directly determined by the characteristics of the strong gravitational field in this region. Both the RPM and (a recently introduced extension of) the ERM (Stuchlík & Kološ 2015) aim at interpreting three QPO signals that have been observed simultaneously in a number of accreting neutron star systems and, so far, in only one accreting BH system, GRO J1655−40. These three signals comprise (1) a low frequency (LF) QPO at ${\nu }_{\mathrm{LF}}$, which is the so-called type C QPO in BH systems and Horizontal Branch QPO in neutron star systems (Casella et al. 2005), with frequencies of up to tens of Hz and (2) twin high frequency (HF) QPOs, at ${\nu }_{\mathrm{lower}}$ and ${\nu }_{\mathrm{upper}}$, with frequencies of several hundred Hz in BHs and around ∼1 kHz NSs. Since these QPOs are detected as incoherent signals in the power spectra of high-time resolution X-ray light curves, their signal-to-noise ratio (S/N) scales linearly with the source count rate and thus with the effective area of X-ray instrumentation. Most currently available QPO measurements have been obtained with the Proportional Counter Array instrument on board the Rossi X-ray Timing Explorer (RXTE/PCA). By exploiting monolithic Silicon Drift Detector technology (Feroci et al. 2010) the next generation X-ray astronomy satellites, which are currently being studied, including LOFT (Feroci et al. 2016), eXTP (Zhang et al. 2016), and STROBE-X (Wilson-Hodge et al. 2017) will achieve an order of magnitude increase in effective area with respect to RXTE/PCA and thus obtain high-precision measurements of simultaneous QPO signals from a variety of BH systems will then become possible. In this paper, we investigate the way in which QPO as measured with the eXTP Large Area Detector (eXTP/LAD, factor of ∼6 larger area than RXTE/PCA) may afford testing the strong-field/high-curvature regime of GR against some alternative theories.

Modified gravity theories can be introduced either by using a bottom-up approach, in which one considers phenomenological parametrizations of BH spacetimes (or of other observable quantities) depending on a set of parameters, or by using a top-down approach, in which specific modifications of GR, possibly inspired by fundamental physics considerations, are adopted (Psaltis 2009). No practical and sufficiently general parametrization of deviations from GR in the strong-field regime has yet been proposed: therefore, we shall follow a top-down approach. Since we are interested in testing the strong-field/large-curvature regime of gravity, we shall consider the so-called "quadratic gravity theories," which are the simplest and most natural modifications of GR in this regime.

In quadratic gravity theories, the Einstein–Hilbert action is modified by including quadratic terms in the curvature tensor, coupled with a scalar field. These couplings can be interpreted as the first terms in an expansion taking into account all possible curvature invariants; such expansion (which is suggested by low-energy effective string theories Gross & Sloan 1987) could make the theory renormalizable (Stelle 1977).

The action (in vacuum), which includes all quadratic curvature invariants, generically coupled to a single scalar field, can be written as (see for instance Yunes & Stein 2011; Berti et al. 2015 and references therein)

$$\begin{eqnarray}S & = & \displaystyle \frac{1}{16\pi }\displaystyle \int \sqrt{-g}{d}^{4}x\left[R-\displaystyle \frac{1}{2}{{\rm{\nabla }}}_{a}\phi {{\rm{\nabla }}}^{a}\phi +{f}_{1}(\phi ){R}^{2}\right.\\ & & +{f}_{2}(\phi ){R}_{{ab}}{R}^{{ab}}+{f}_{3}(\phi ){R}_{{abcd}}{R}^{{abcd}}+{f}_{4}{(\phi )}^{* }{RR}],\end{eqnarray} \tag{ 1 }$$

where ${f}_{i}(\phi )$ ($i=1,2,3,4$) are generic coupling functions, ${}^{* }{RR}\equiv \tfrac{1}{2}{R}_{\mu \nu \rho \sigma }{\epsilon }^{\nu \mu \lambda \kappa }$ ${R}^{\rho \sigma }{}_{\lambda \kappa }$, with ${\epsilon }^{\mu \nu \rho \sigma }$ the Levi-Civita tensor. Two relevant cases of this class of theories are (1) Einstein–dilaton–Gauss–Bonnet (EDGB) gravity (${f}_{1}=\tfrac{\alpha }{4}{e}^{\phi }$, ${f}_{2}=-4{f}_{1}$, ${f}_{3}={f}_{1}$, ${f}_{4}=0$), in which the quadratic corrections reduce to the Gauss–Bonnet invariant, ${R}_{\mathrm{GB}}^{2}\equiv {R}^{2}-4{R}_{\mu \nu }^{2}+{R}_{\mu \nu \rho \sigma }^{2}$ (Kanti et al. 1996; Moura & Schiappa 2007); and (2) Dynamical Chern–Simons (DCS) gravity (Jackiw & Pi 2003; Alexander & Yunes 2009; Delsate et al. 2015; ${f}_{1}={f}_{2}={f}_{3}=0$ and ${f}_{4}=\tfrac{\beta }{4}\phi$).

In general, the equations of motion of the action (1) have third-(or higher-)order derivatives, and the theory is subject to Ostrogradsky's instability (Woodard 2007). To avoid this feature, the theory should be treated as an effective-field theory, valid only up to second order in the curvature, in the limit of small couplings f_i. In this way, ghosts and other pathologies disappear (for a discussion, see Berti et al. 2015). This limit also requires us to expand the functions ${f}_{i}(\phi )$ up to linear order in ϕ, i.e., ${f}_{i}(\phi )\approx {\eta }_{i}+\tfrac{{\alpha }_{i}}{4}\phi$, with ${\eta }_{i},{\alpha }_{i}$ such that the corrections are small compared to the leading Einstein–Hilbert term (Pani et al. 2011; Yunes & Stein 2011; Berti et al. 2015). The only exception is EDGB gravity, whose equations of motion are second order, avoiding Ostrogradsky's instability. Therefore, EDGB gravity can be treated as an "exact" theory, and its coupling constant α can, in principle, be a finite quantity.

BH solutions to theories based on action (1) have been found in various particular cases by Mignemi & Stewart (1993), Kanti et al. (1996), Pani & Cardoso (2009), Yunes & Pretorius (2009), Yunes & Stein (2011), Kleihaus et al. (2011), Yagi et al. (2012), and Kleihaus et al. (2014). Stationary, axisymmetric, BH solutions can be found in closed analytical form to any order in a small-spin and small-coupling expansion (Pani et al. 2011; Ayzenberg & Yunes 2014; Maselli et al. 2015b). Remarkably, to leading order in the coupling, the metric depends only on two constants, ${\alpha }_{3}\equiv {\alpha }_{\mathrm{GB}}=\alpha$ and ${\alpha }_{4}\equiv {\alpha }_{\mathrm{DCS}}=\beta$, whereas it is independent of ${\eta }_{i}$ ($i\,=1,2,3,4$) and of ${\alpha }_{\mathrm{1,2}}$. Thus, all stationary solutions to the effective-field theory introduced above reduce, to second order in the coupling parameter, to two families, namely BH solutions to EDGB gravity and to DCS gravity (or, at most, solutions to theories with different choices of ${f}_{3}(\phi )$, ${f}_{4}(\phi );$ we remark that these theories are equivalent to EDGB and DCS gravity in the small-coupling limit). The former case is the only one that can be defined beyond the small-coupling approximation; BH solutions in this case were obtained numerically for generic spin and coupling by Kleihaus et al. (2011, 2014) and in closed analytical form to fifth order in the spin and to seventh order in the coupling parameter by Maselli et al. (2015b). In the latter case, spinning DCS BHs were first obtained by Yunes & Pretorius (2009) and Konno et al. (2009) to leading order in the spin and by Yagi et al. (2012) to quadratic order in the spin.

In this article, we discuss the possibility of using QPOs as tools to test GR against quadratic gravity theories, extending previous results from Vincent (2013) and Maselli et al. (2015a, 2015b). It is worth remarking that the QPO diagnostics has also been exploited in Bambi (2012, 2015), and in Franchini et al. (2017) as a probe of hairy BH solutions that emerge in standard GR in the presence of ultralight fields.

This paper is organized as follows. In Section 2, we work out the characteristic orbital frequencies—i.e., the azimuthal and epicyclic frequencies—for rotating BHs in EDGB gravity and DCS gravity. Since DCS gravity has to be considered as an effective theory, in our analysis, we assume β to be a small quantity, while we allow for finite values of the EDGB coupling parameter, though a theoretical bound does exist for α, i.e.,

$$\begin{eqnarray}&&0\lt \displaystyle \frac{\alpha }{{M}^{2}}\lesssim 0.691,\end{eqnarray} \tag{ 2 }$$

where M is the BH mass.

In this paper, we shall employ the BH solution in EDGB gravity derived in Maselli et al. (2015b) to fifth order in the spin and to seventh order in α, and the BH solution in DCS gravity to fifth order in the spin (thus extending the results of Yagi et al. 2012) and to second order in the coupling parameter β. For the latter, the explicit expressions of the metric tensor and the scalar field are quite long and are available in a Mathematica notebook provided in this tar.gz package. In Section 3, we discuss the QPO models we adopt in our analysis, in which the observed frequencies are related to the azimuthal and epicyclic frequencies of BH geodesics. In Section 4, we discuss the method we propose in order to test GR against EDGB gravity and DCS gravity by using future high-precision QPO measurements. We show that, using this approach, large-area X-ray instrumentation such as that being studied for LOFT-P, eXTP, and STROBE-X can significantly constrain the parameter space of these theories.

We use geometric units in which $G=c=1$ and consider a spin parameter ${a}^{\star }=J/{M}^{2}\lt 0.5$, such that truncation errors are expected to be of ${ \mathcal O }({a}^{\star 6})\approx 1 \%$ at the most.

In thin accretion disks around a rotating BH the stream of matter follows nearly equatorial $(\theta =\pi /2)$ and nearly circular geodesics at a specific radius r₀. For an axially symmetric spacetime described by the line element

$$\begin{eqnarray}&&{{ds}}^{2}={g}_{{tt}}{{dt}}^{2}+{g}_{{rr}}{{dr}}^{2}+{g}_{\theta \theta }d{\theta }^{2}+{g}_{t\phi }{dtd}\phi +{g}_{\phi \phi }d{\phi }^{2},\end{eqnarray} \tag{ 3 }$$

small perturbations of the particle trajectories along the radial and vertical directions $r={r}_{0}+\delta r$ and $\theta =\pi /2+\delta \theta$, lead to oscillations around the equilibrium configuration characterized by the two epicyclic frequencies

$$\begin{eqnarray}&&{\nu }_{r}^{2}=\displaystyle \frac{1}{{(2\pi )}^{2}}\displaystyle \frac{{({g}_{{tt}}+{\rm{\Omega }}{g}_{t\phi })}^{2}}{2{g}_{{rr}}}\displaystyle \frac{{\partial }^{2}{ \mathcal U }}{\partial {r}^{2}}\left({r}_{0},\displaystyle \frac{\pi }{2}\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\nu }_{\theta }^{2}=\displaystyle \frac{1}{{(2\pi )}^{2}}\displaystyle \frac{{({g}_{{tt}}+{\rm{\Omega }}{g}_{t\phi })}^{2}}{2{g}_{\theta \theta }}\displaystyle \frac{{\partial }^{2}{ \mathcal U }}{\partial {\theta }^{2}}\left({r}_{0},\displaystyle \frac{\pi }{2}\right),\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Omega }}=2\pi {\nu }_{\phi }$ is the particle angular velocity and ${\nu }_{\phi }$ its azimuthal frequency (see Maselli et al. 2015a for technical details). The effective potential ${ \mathcal U }(r,\theta )={g}^{{tt}}-2{{lg}}^{t\phi }+{l}^{2}{g}^{\phi \phi }$ depends on the metric functions and the ratio between the particle's angular momentum and energy l = L/E per unit mass. For a specific radius r, the three frequencies of a Kerr BH $({\nu }_{\phi },{\nu }_{\theta },{\nu }_{r})$ are functions of the mass and spin parameter ${a}^{\star }=J/{M}^{2}\in [-1,1]$ only, where J is the BH angular momentum:

$$\begin{eqnarray}&&{\nu }_{\phi }^{\mathrm{GR}}=\displaystyle \frac{1}{2\pi }\displaystyle \frac{{M}^{1/2}}{{r}^{3/2}+{a}^{\star }{M}^{3/2}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\nu }_{r}^{\mathrm{GR}}={\nu }_{\phi }^{\mathrm{GR}}{\left(1-\displaystyle \frac{6M}{r}+8{a}^{\star }\displaystyle \frac{{M}^{3/2}}{{r}^{3/2}}-3{{a}^{\star }}^{2}\displaystyle \frac{{M}^{2}}{{r}^{2}}\right)}^{1/2},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\nu }_{\theta }^{\mathrm{GR}}={\nu }_{\phi }^{\mathrm{GR}}{\left(1-4{a}^{\star }\displaystyle \frac{{M}^{3/2}}{{r}^{3/2}}+3{{a}^{\star }}^{2}\displaystyle \frac{{M}^{2}}{{r}^{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 8 }$$

In Newtonian gravity, the three frequencies coincide, while in GR ${\nu }_{\phi }\geqslant {\nu }_{\theta }\gt {\nu }_{r}$. In particular, for ${a}^{\star }=0$ the azimuthal and vertical components are equal, whereas the radial one vanishes at the ISCO. As an example, in Figure 1, we show these quantities for a BH with $M=10\,{M}_{\odot }$ and spin parameter ${a}^{\star }=(0,0.7)$.

Zoom In Zoom Out Reset image size

Epicyclic frequencies in the EDGB gravity have been computed in Maselli et al. (2015a) for slowly rotating BHs at the linear order in the angular momentum, as a function of the coupling constant $\alpha /{M}^{2};$ they show differences with respect to the GR case, which increase for higher BH spin, and are potentially observable by future large-area X-ray satellites. Motivated by these results, Maselli et al. (2015b) improved the templates for $({\nu }_{\phi },{\nu }_{\theta },{\nu }_{r})$ by extending them up to the fifth order in ${a}^{\star }$. Based on this higher-order expansion, more rapidly spinning BHs can be considered, thus exploring regions in parameter space that are amenable to show significant departures from GR.

To be consistent with the formalism developed for the EDGB and the DCS theories, in our analysis, we will expand Equations (6)–(8) as power series of ${a}^{\star }$ neglecting terms that are ${ \mathcal O }({a}^{\star 6})$ and higher. As an example, in Figure 2, we plot the relative percentage difference between the epicyclic frequencies ${\nu }_{\phi }$ and ${\nu }_{r}$ computed at $r=1.1{r}_{\mathrm{ISCO}}$ in GR and in EDGB or DCS theory:

$$\begin{eqnarray}&&{\epsilon }_{{\nu }_{\phi }}=\displaystyle \frac{{\nu }_{\phi }^{\mathrm{EDGB},\mathrm{DCS}}-{\nu }_{\phi }^{\mathrm{GR}}}{{\nu }_{\phi }^{\mathrm{GR}}},\quad {\epsilon }_{{\nu }_{r}}=\displaystyle \frac{{\nu }_{\phi }^{\mathrm{EDGB},\mathrm{DCS}}-{\nu }_{r}^{\mathrm{GR}}}{{\nu }_{r}^{\mathrm{GR}}},\end{eqnarray} \tag{ 9 }$$

as functions of the coupling parameters and of the BH spin. As expected, for larger values of $(\alpha /{M}^{2},\beta /{M}^{2})$ and faster rotation rates the relative difference increases; it can be as high as 3% for the equatorial frequency in EDGB gravity. For the vertical component ${\nu }_{\theta }$ the relative difference is of the same order as that of ${\nu }_{\phi }$. These values decrease when the frequencies are computed in DCS. The bottom panel of Figure 2 shows that in this case ${\epsilon }_{{\nu }_{\phi }}\ll 1$ for all coupling parameter values allowed by the theory.

Zoom In Zoom Out Reset image size

The epicyclic frequencies are the basic ingredients of the geodesic models that we shall describe in this section: the RPM (Stella et al. 1999) and the modified ERM (Török et al. 2005). Both approaches interpret the simultaneous occurrence of the LF and twin HF QPOs in terms of geodesic frequencies.

According to the RPM, the upper and lower HF QPOs coincide with the azimuthal frequency ${\nu }_{\phi }$, and the periastron precession frequency, ${\nu }_{\mathrm{per}}={\nu }_{\phi }-{\nu }_{r}$. The LF QPO mode instead, is identified with the nodal precession frequency, ${\nu }_{\mathrm{nod}}={\nu }_{\phi }-{\nu }_{\theta }$. These three QPO signals $({\nu }_{\phi },{\nu }_{\mathrm{per}},{\nu }_{\mathrm{nod}})$ are assumed to be generated at the same orbital radius. Although the first application of the RPM to BH systems can be traced back to the original paper (Stella et al. 1999), the first complete exploitation of the model has been made possible by the discovery of three simultaneous QPOs in GRO J1655−40, as measured by the RXTE/PCA with 0.5%–1.5% accuracy

$$\begin{eqnarray}&&{\nu }_{\phi }={441}_{-2}^{+2}\,\mathrm{Hz},{\nu }_{\mathrm{per}}={298}_{-4}^{+4}\,\mathrm{Hz},{\nu }_{\mathrm{nod}}={17.3}_{-0.1}^{+0.1}\,\mathrm{Hz}.\end{eqnarray} \tag{ 10 }$$

By assuming the Kerr metric and fitting the three values in terms of the RPM frequencies, a precise estimate of the emission radius $r=(5.68\,\pm 0.04)M$, and of the BH mass $M=(5.31\pm 0.07)\,{M}_{\odot }$ and angular momentum ${a}^{\star }=J/{M}^{2}\,=0.290\pm 0.003$ was obtained (Motta et al. 2014a). Independent optical observations of the same source lead to mass values in good agreement with those obtained through the RPM $[{M}_{\mathrm{opt}}=(5.4\pm 0.3)\,{M}_{\odot }].$ However, there exist systematic uncertainties in different methods for estimating BH spin, as spectral continuum measurements and analysis of the Fe spectral line profiles predict BH spin in the range $0.65\lt {a}_{\mathrm{sc}}^{\star }\lt 0.75$ and $0.94\lt {a}_{\mathrm{Fe}}^{\star }\lt 0.98$, respectively (Motta et al. 2014b), which show a discrepancy with the QPO-based value and among each other.

The ERM builds on the possibility that in BHs with twin HF QPOs, their centroid frequencies are such that ${\nu }_{\mathrm{upper}}/{\nu }_{\mathrm{lower}}=3/2$ (Török et al. 2005). This suggests an underlying mechanism based on nonlinear resonances. At the first order in the vertical and radial displacements, deviations $\delta r$ and $\delta \theta$ from geodesic circular motion can be described by the following equations:

$$\begin{eqnarray}&&\ddot{\delta r}+{\omega }_{r}^{2}{\delta }_{r}=\delta {a}_{r},\quad \ddot{\delta \theta }+{\omega }_{\theta }^{2}{\delta }_{r}=\delta {a}_{\theta },\end{eqnarray} \tag{ 11 }$$

where dots represent time derivatives, ${\omega }_{\theta ,r}=2\pi {\nu }_{r,\theta }$, and $\delta {a}_{r},\delta {a}_{\theta }$ are forcing terms. Nonlinear resonances show interesting common features with QPOs, e.g., they occur only over a finite $\delta \nu$, allow for frequency combinations, and for sub-harmonic modes. If $\delta r\gg \delta \theta$, the mixing term $\delta \theta \delta r$ cannot be neglected and it must be included in the linearized equation for the vertical component (11). For $\delta {a}_{r}=\delta {a}_{\theta }=0$, the radial component yields the solution ${\delta }_{r}=A\cos ({\omega }_{r}t)$, while the vertical oscillation obeys the Mathieu equation

$$\begin{eqnarray}&&\ddot{\delta \theta }+{\omega }_{\theta }^{2}[1+{Ah}\cos ({\omega }_{r}t)]\delta \theta =0,\end{eqnarray} \tag{ 12 }$$

where A and h are known constants. The above equation describes parametric resonances such that

$$\begin{eqnarray}&&\displaystyle \frac{{\nu }_{r}}{{\nu }_{\theta }}=\displaystyle \frac{2}{n},\quad n=1,2,3\ldots .\end{eqnarray} \tag{ 13 }$$

In GR, ${\nu }_{\theta }\gt {\nu }_{r}$ and therefore we may associate the vertical component to the larger of the HF QPOs. It is interesting to note that a nonzero forcing term along the θ-direction, i.e., $\delta {a}_{\theta }\ne 0$, allows for a combination of the frequencies

$$\begin{eqnarray}&&{\nu }_{-}={\nu }_{\theta }-{\nu }_{r},\quad {\nu }_{+}={\nu }_{\theta }+{\nu }_{r},\end{eqnarray} \tag{ 14 }$$

which still satisfies the observational evidence for a small integer ratio ${\nu }_{\mathrm{upper}}/{\nu }_{\mathrm{down}}=3/2$, as long as

$$\begin{eqnarray}&&{\nu }_{\mathrm{upper}}={\nu }_{\theta },\quad {\nu }_{\mathrm{lower}}={\nu }_{-},\end{eqnarray} \tag{ 15 }$$

oder

$$\begin{eqnarray}&&{\nu }_{\mathrm{upper}}={\nu }_{+},\quad {\nu }_{\mathrm{lower}}={\nu }_{\theta }\end{eqnarray} \tag{ 16 }$$

are considered. The ERM, originally developed to interpret only the twin HF QPOs, has been extended to interpret the LF QPO mode in terms of ${\nu }_{\mathrm{nod}}$ (as in the RPM), and thus interpret the simultaneous occurrence of the three frequencies associated with GRO J1655−40 (Stuchlík & Kološ 2016). Among all possibilities discussed by these authors, we consider here the case in which the two HF QPOs are identified with ${\nu }_{\mathrm{upper}}={\nu }_{\theta }$ and ${\nu }_{\mathrm{lower}}={\nu }_{-}$. With this choice, the BH parameters have been constrained to $M=(5.1\pm 0.1)\,{M}_{\odot }$, ${a}^{\star }=0.274\pm 0.003$, and $r=(5.67\,\pm 0.05)M$, which are reasonably close to those derived with the RPM. Similar results hold also for different combinations, since the fundamental scale of the effect is set by the ISCO frequency.

As an example, in Figure 3, we show the two sets of frequencies $({\nu }_{\phi },{\nu }_{\mathrm{per}},{\nu }_{\mathrm{nod}})$ and $({\nu }_{\theta },{\nu }_{\mathrm{nod}},{\nu }_{-})$ employed by the RPM and the ERM approaches. All values are computed in GR for a $5\,{M}_{\odot }$ BH with different spin, as a function of radius in ${r}_{\mathrm{Isco}}$ units.

Zoom In Zoom Out Reset image size

To test the ability of eXTP/LAD to constrain quadratic theories of gravity, we follow the same data analysis procedure described in Maselli et al. (2015a), based on simulations of two QPO triplets with different values of their frequencies. In the following, we briefly summarize the basic steps of this method.

We consider a prototype BH of mass $\bar{M}=5.3\,{M}_{\odot }$ and selected values of its ${\bar{a}}^{\star }$, and adopt different values of the coupling parameter for each modified gravity theory. By using the analytic relations derived in Maselli et al. (2015b), we calculate two sets of epicyclic frequencies for two different emission radii ${r}_{\mathrm{1,2}}$, in both EDGB or DCS gravity. Based on these sets, we calculate the three QPO signals expected in the case of the RPM and ERM; we then use the corresponding two QPO triplets as center values to generate $N={10}^{5}$ samples from Gaussian distributions with standard deviations $({\sigma }_{\varphi },{\sigma }_{\mathrm{nod}},{\sigma }_{\mathrm{per}})$, obtained by rescaling the error bars in Equation (9) by the ratio of the RXTE/PCA and eXTP/LAD effective areas. Finally, we use the geodesic frequencies for a Kerr BH (Equations (6)–(8)) to compute $2N$ triplets of ${(M,{a}^{\star },r)}_{j=\mathrm{1,2}}$, from which the mean values of the BH mass, spin, and emission radii, and the covariance matrices ${{\rm{\Sigma }}}_{1},{{\rm{\Sigma }}}_{2}$ associated with the two sets are derived.

In order to assess whether these distributions are consistent with ${M}_{1}={M}_{2}$ and ${a}_{1}^{\star }={a}_{2}^{\star }$, we define the three variables

$$\begin{eqnarray}&&{\rm{\Delta }}M={M}_{1}-{M}_{2},{\rm{\Delta }}{a}^{\star }={a}_{1}^{\star }-{a}_{2}^{\star },{\rm{\Delta }}r={r}_{1}-{r}_{2}\end{eqnarray} \tag{ 17 }$$

checking that the normal distribution ${ \mathcal N }({\boldsymbol{\mu }},{\rm{\Sigma }})$, with ${\boldsymbol{\mu }}=({\rm{\Delta }}M,{\rm{\Delta }}{a}^{\star },{\rm{\Delta }}r)$ and ${\rm{\Sigma }}={{\rm{\Sigma }}}_{1}+{{\rm{\Sigma }}}_{2}$, is consistent with a Gaussian of zero mean. Given the ${\chi }^{2}$ variable

$$\begin{eqnarray}&&{\chi }^{2}={({\boldsymbol{x}}-{\boldsymbol{\mu }})}^{{\rm{T}}}{{\rm{\Sigma }}}^{-1}({\boldsymbol{x}}-{\boldsymbol{\mu }}),\end{eqnarray} \tag{ 18 }$$

the values of ${\chi }^{2}=\delta$, define the confidence levels (CL) of $({\rm{\Delta }}M,{\rm{\Delta }}{a}^{\star },{\rm{\Delta }}r)$ for a specific choice of δ. We repeat this analysis assuming both the RPM and the ERM to interpret the QPOs, assuming the observed triplet as given by ${({\nu }_{\phi },{\nu }_{\mathrm{nod}},{\nu }_{\mathrm{per}})}_{i=\mathrm{1,2}}$ and ${({\nu }_{\theta },{\nu }_{-},{\nu }_{\mathrm{nod}})}_{i=\mathrm{1,2}}$, respectively. In the RPM the location of the emission radii can be, in principle, chosen freely; following Maselli et al. (2015a), we adopt ${r}_{1}=1.1{r}_{\mathrm{Isco}}$ and ${r}_{2}=1.4{r}_{\mathrm{Isco}}$, which give rise to QPO frequencies comparable to those observed in GRO J1655−40. (Somewhat different choices of the radial coordinates, as long as they are $\leqslant 2{r}_{\mathrm{Isco}}$, do not alter significantly our results.) In the ERM approach, instead, for each resonance there is only one radius that satisfies condition (13). Specifically, we choose r₁ and r₂ in order to obtain ${\nu }_{\theta }/{\nu }_{-}=3/2$ and ${\nu }_{\theta }/{\nu }_{-}=5/2$.

The top panels of Figure 4 show the CL at $1\sigma$, $2\sigma$, and $3\sigma$ in the $({\rm{\Delta }}M,{\rm{\Delta }}{a}^{\star })$ parameter space, for BHs with ${\bar{a}}^{\star }=0.5$ and $\alpha /{\bar{M}}^{2}=(0.6,0.4,0)$, in the case in which the QPO frequencies are described by the RPM. The red cross identifies the origin of the plane, and corresponds to the null hypothesis, for which the two samples of QPO triplets derive from the same distribution, i.e., also the input data are generated from the Kerr metric. The panel on the right shows the case $\alpha /{\bar{M}}^{2}=0$ as a consistency check of our method; it is apparent that the three CLs are centered on $({\rm{\Delta }}M,{\rm{\Delta }}{a}^{\star })=(0,0)$, and thus that the data are consistent with ${M}_{1}={M}_{2}$ and ${a}_{1}^{\star }={a}_{2}^{\star }$. The left and central plots refer to BH configurations with $\alpha /{\bar{M}}^{2}\gt 0$: for these models ${\rm{\Delta }}M$ and ${\rm{\Delta }}{a}^{\star }$ are both incompatible with 0. In particular, the central upper panel of Figure 4 shows that if the coupling parameter is $\alpha /{\bar{M}}^{2}=0.4$, the data are already incompatible with the Kerr metric at more than the $3\sigma$ level. Higher BH spins (${a}^{\star }\gt 0.5$) would exclude even smaller values of the coupling parameter. To clarify this point, we carry out the same analysis by varying the BH spin. Figure 5 shows the $3\sigma$ CL for $\alpha /{\bar{M}}^{2}=(0.6,0.4)$ and ${a}^{\star }=(0.3,0.4,0.5)$. For a fixed value of the EDGB coupling parameter, increasing BH spins make simulated QPO frequencies depart more and more from GR predictions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bottom panels of Figure 4 show the CLs for the same BH configurations described above and QPO frequencies calculated in accordance with the ERM. The central panel is for $\alpha /{\bar{M}}^{2}\sim 0.4$, i.e., the value excluded at slightly more than $3\sigma$ CL also in this case. In general, the two frameworks provide similar constraints, as shown in Figure 6, in which we make a direct comparison between the two QPO models. We note that the ellipses obtained for the ERM are slightly larger than those computed for the RPM, irregardless of the EDGB coupling parameter. We have also computed the ellipses for RPM and ERM choosing the same r₂, but the difference between the two approaches still persists.

Zoom In Zoom Out Reset image size

Much different results are obtained when azimuthal and the epicyclic frequencies are generated by using DCS theory. The results of our analysis are shown in Figure 7, for a BH with the same mass and spin considered before (i.e., $\bar{M}=5.3\,{M}_{\odot }$, ${a}^{\star }=0.5$), and DCS coupling parameter $\beta /{\bar{M}}^{2}=0.1$. We emphasize that the latter value is the highest compatible with a small-coupling approximation of an effective theory of gravity, free of pathologies. It is apparent that the cross corresponding to the null hypothesis (i.e., GR-Kerr frequencies apply) lies between 2σ and 3σ CL. Therefore, in this case, we would not be able to set statistically significant constraints on the theory, even for the maximum value $\beta /{\bar{M}}^{2}=0.1$.

Zoom In Zoom Out Reset image size

Finally, we consider the outcome of our study, when QPO measurements with LOFT-type S/N (as in Maselli et al. 2015a) are also considered. The two plots of Figure 8 show the $3-\sigma$ ellipses computed for the RPM in EDGB and DCS and for both eXTP-type and LOFT-type QPO data (for the effective area of LOFT/LAD we use $\times 15$ the effective area of RXTE/PCA, as in Maselli et al. 2015a). In both gravity theories, a higher S/N leads to a significant improvement of the excluded parameter space. This is especially evident for DCS, which can now be constrained even for $\beta /{\bar{M}}^{2}=0.1$.

Zoom In Zoom Out Reset image size

We note that solar system experiments based on measurements of the geodesic precession and tests of the Newton's law, have already constrained the DCS coupling to $\sqrt{\beta }\lt { \mathcal O }({10}^{8}\,\mathrm{km})$ (Ali-Haimoud & Chen 2011), which is much weaker than the requirement $\beta /{M}^{2}\ll 1$ for $M\ll {10}^{8}\,{M}_{\odot }$, i.e., for the stellar mass sources considered in this paper. For EDGB gravity, the best bound on the coupling parameter comes from observations of low mass X-ray binaries, leading to $\alpha \lesssim 47\,{M}_{\odot }^{2}$ (Yagi 2012), which is looser than the theoretical bound (2) for $M\lesssim 8.2\,{M}_{\odot }$.

X-ray QPOs emitted by accreting BHs represents a very promising tool to investigate stationary spacetimes in a genuine strong-field and high-curvature regime, even though their interpretation is still open to different possibilities. Future instrumentation such as the eXTP/LAD holds the promise not only to shed light on the origin of QPOs, but also to constrain modified gravity theories in extreme astrophysical environments.

In this paper, we have extended the analysis presented in Maselli et al. (2015a) in several directions: (1) by considering DCS gravity in addition to EDGB gravity; in an effective-field-theory approach, these are the only quadratic gravity theories admitting BH solutions other than Kerr; (2) by calculating accurate geodesic frequencies for higher values of BH spin, up to ${a}^{\star }=0.5$ (as opposed to 0.2); (3) by adopting two different geodesic QPO models, the RPM and the ERM. Most of our simulations here were carried out adopting the S/N expected for an X-ray instrument of effective area ($\sim 3.5\ {{\rm{m}}}^{2}$), comparable to that envisaged for the eXTP/LAD instrument (which is about a factor of ∼6 larger than that of the RXTE/PCA). This is at variance with the simulations in Maselli et al. (2015a), which used instead a larger effective ($\sim 9\ {{\rm{m}}}^{2}$ close to that being studied for LOFT).

Our results can be summarized as follows.

• 1.  
  Both the RPM and the ERM models provide viable frameworks to test alternative gravity theories to GR, through high-precision X-ray measurements of QPOs.
• 2.  
  Even for moderately fast-rotating BHs, with spin parameter ${a}^{\star }\sim 0.5$, an eXTP-type mission may set stringent constraints on the EDGB coupling parameter.
• 3.  
  As already noted in Maselli et al. (2015a), our ability to distinguish between GR and EDGB increases with the spin of the accreting object. By extrapolating our results here we suggest that maximally spinning BHs, with ${a}^{\star }\sim 1$, are the best probes of quadratic gravity.
• 4.  
  With an eXTP-type mission, none of the models considered here set useful bounds on DCS theory. This is mainly due to the small-coupling parameter ($\leqslant 0.1$), which is consistent with the requirement that DCS be an effective theory of gravity, and to the fact that DCS gravity introduces smaller corrections to GR BH geometries.
• 5.  
  A LOFT-like observatory, with its improved S/N, would constrain more tightly the parameter space of both modified theories we have considered. In particular, this would allow us to set new bounds on DCS gravity, which would otherwise be insensitive to the QPO diagnostic.

Our study shows that the application of geodesic models to future high S/N QPO measurements holds the potential to test alternative gravity theories in the strong-field regime, with the results appearing to be especially promising for EDGB theory.

P.P. acknowledges support from FCT-Portugal through the project IF/00293/2013. L.S. acknowledges funding from the ASI-INAF contract I/004/11/1. This work was supported by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904 and by the COST action CA1610 "GWverse."