# Quasi-periodic oscillations from relativistic ray-traced hydrodynamical tori

B. Mishra, F.H. Vincent, A. Manousakis, P.C. Fragile, T. Paumard, W.

Kluźniak

► To cite this version:

B. Mishra, F.H. Vincent, A. Manousakis, P.C. Fragile, T. Paumard, et al. Quasi-periodic oscillations from relativistic ray-traced hydrodynamical tori. Monthly Notices of the Royal Astronomical Society, 2017, 467 (4), pp.4036-4049. 10.1093/mnras/stx299. hal-01582852

HAL Id: hal-01582852

https://hal.science/hal-01582852v1

Submitted on 4 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Quasi-periodic oscillations from relativistic ray-traced hydrodynamical tori

B. Mishra,<fnref n="1" />* F. H. Vincent,<fnref n="1" />* A. Manousakis,<fnref n="1" />* P. C. Fragile,<fnref n="2" /> T. Paumard<fnref n="3" /> and W. Kluźniak<fnref n="1" />

Accepted 2017 January 31. Received 2017 January 31; in original form 2015 October 23

## ABSTRACT

We simulate an oscillating polytropic torus with constant specific angular momentum around a Schwarzschild black hole. The goal is to search for quasi-periodic oscillations in the light curve of the torus. The torus is subjected to an initial radial, vertical or diagonal (combination of radial and vertical) velocity perturbation. The hydrodynamical simulations are performed using the general relativistic magnetohydrodynamics code COSMOS++ and ray traced using the general relativistic ray-tracing code GYOTO. We found that a horizontal velocity perturbation triggers the radial and plus modes, whereas a vertical velocity perturbation triggers the vertical and X modes. The diagonal perturbation gives a combination of the modes triggered in the radial and vertical cases. The breathing mode is excited in all cases. We find that the radial eigenmode frequency is independent of the equation of state used. We conclude that the radial and plus modes we find in our simulations are identical to f and $o_1$ modes reported in Rezolla et al.

**Key words:** accretion, accretion disc – hydrodynamics – radiative transfer – stars: black holes.

## 1 INTRODUCTION

Black hole binaries, and especially black holes in low-mass X-ray binaries (LMXBs), have been extensively studied in the past decades (for a review, see Lewin & van der Klis 2010; Remillard & McClintock 2006). The advent of the *Rossi X-ray Timing Explorer* (RXTE) in the late 1990s, in particular, allowed the study of rapid quasi-periodic oscillations (QPOs) in these sources, as well as in accreting neutron stars. QPOs appear as peaks (of non-zero width, hence their quasi-periodic nature) in the power density spectra (PDS, the simplest expression of which is the square modulus of the Fourier transform) of some black hole binaries.

Prior to the launch of the *RXTE* mission, it had been suggested that QPOs may be observed in LMXBs in the kilohertz range as a result of inhomogeneities of flow in the inner accretion disc around neutron stars and that they could be used to test predictions of Einstein's Theory of General Relativity (GR) in the strong field regime (Kluzniak, Michelson & Wagoner 1990). In this paper, we will be only interested in the high-frequency QPOs in black hole sources, with typical frequencies of 40–450 Hz (Morgan, Remillard & Greiner 1997; Remillard et al. 1999, 2002; Belloni, Méndez &

Sánchez-Fernández 2001; Miller et al. 2001; Strohmayer 2001a,b; Altamirano & Belloni 2012; Belloni, Sanna & Méndez 2012). These QPOs exhibit characteristic frequencies of the order of the Keplerian frequency (i.e. the orbital frequency on a circular time-like geodesic in the equatorial plane) very close to the central black hole. The highest stable Keplerian frequency is reached at the marginally stable orbit, usually called the innermost stable circular orbit (ISCO) in the current literature. However, in a thick accretion disc or in a torus, in principle, the orbiting fluid can exhibit quasi-stable motion at frequencies as high as that in the marginally bound orbit. The ISCO frequency of a $10 M_{\odot}$ Schwarzschild black hole is ≈220 Hz, which is in the typical range for high-frequency QPOs of black hole binaries. Consequently, it is reasonable to investigate whether QPOs might be due to phenomena arising in the strong-field region very close to the black hole event horizon. If this picture is true, QPOs could potentially be used to constrain the black holes properties and in particular its spin parameter (Abramowicz & Kluźniak 2001) and might be used as probes of strong-field GR in the relatively distant future, when precise enough data will allow for comparison of various theories of gravity (see e.g. Johannsen & Psaltis 2011).

In order to understand the high-frequency variability, numerous models have been proposed in the last few years. However, we are still far from understanding these phenomena. In black holes, high-frequency QPOs are associated with the so-called steep power-law spectral state of X-ray binaries (Remillard & McClintock 2006).

<footnote n="1">¹Nicolaus Copernicus Astronomical Center, Bartycka 18, Warsaw, PL-00-716, Poland</footnote>

<footnote n="2">²Physics and Astronomy College of Charleston, 66 George Street, Charleston, SC 29424, USA</footnote>

<footnote n="3">³LESIA, Observatoire de Paris, CNRS, UPMC, Université Paris-Diderot, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote>* E-mail: mbhupe@camk.edu.pl (BM); frederic.vincent@obspm.fr (FHV); antonism@camk.edu.pl (AM)</footnote>
This spectral state itself is not well understood. The presence of QPOs in this state may be related to Comptonization of disc photons in a hot toroidal corona (Bursa & Adámek 2014), but we will explore other possibilities as well. Recently, Dexter & Blaes (2014) advocated a different model that accounts for both the steep power-law spectral state and approximately for the 3:2 frequency ratio of twin-peak QPOs. One goal of investigations such as ours is to try to narrow down the choice between models for high-frequency QPOs by exploring their ability to account for observations. As of today, a large diversity of models is still advocated. Blobs of matter orbiting in the inner accretion disc have been an early favourite as the cause of a quasi-periodic modulation to the X-ray flux (Cui, Zhang & Chen 1998; Kluźniak 1998; Stella & Vietri 1998), but these have been ruled out (Barret et al. 2005), at least in the neutron star LMXBs, by the high degree of coherence of the oscillation (high-quality factor *Q*). Although in black holes the observed *Q* value is low,<fnref n="1" /> so blobs are not ruled out, the presence of QPOs probably reflects modulations imprinted on the photon flux by intrinsic disc oscillations, whose normal modes have been studied extensively (Wagoner 1999; Kato 2001; Miranda, Horák & Lai 2015), or by normal modes of other accretion structures, e.g. tori.

Upon analysing certain properties of twin-kilohertz QPOs in neutron stars, Kluźniak & Abramowicz (2001) proposed that they could be accounted for by a resonance between two possible motions of the accreting fluid. The resonance model predicted that high-frequency QPOs in black holes should occur in pairs, as they do in neutron stars, and that the two QPO frequencies should be in the ratio of small integer numbers, e.g. 2:1 or 3:2 ratio. These predictions were born out with the discovery of a 450 Hz QPO in the binary black hole system GRO J1655-40 (Strohmayer 2001a), in which a 300 Hz QPO had been reported previously. Abramowicz & Kluźniak (2001) pointed out that the two frequencies, 450 Hz and 300 Hz, are in a 3:2 ratio, and constrained the spin of the black hole assuming a particular resonance model. The resonance model has been used to relate the black hole spin to the frequencies expected also for other types of resonances that can occur in nearly Keplerian discs in strong gravity (Török et al. 2005; Stuchlík, Kotrlová & Török 2013). In spite of the successful predictions of the resonance model, it is important to examine other possibilities, such as that the ratio may only be approximately 3:2, and may simply reflect the ratio of two eigenfrequencies of an accretion structure, as already noted by Blaes, Arras & Fragile (2006).

Based on the resonance model of Abramowicz & Kluźniak (2001), it has been suggested that axisymmetric, up and down, motion at the vertical epicyclic frequency can also be directly observable if it occurs in the inner regions of the disc. In particular, a torus oscillating along its own axis may cause periodic variations in its radiative flux observed at infinity (Bursa et al. 2004; Kluźniak, Abramowicz & Lee 2004; Kluźniak et al. 2004; Lee, Abramowicz & Kluźniak 2004). Such toroidal structures may be formed spontaneously in the magnetohydrodynamic (MHD) evolution of an accretion disc, (e.g. Machida & Matsumoto 2008). Rezzolla, Yoshida & Zanotti (2003) performed a numerical study of radially oscillating tori. To produce synthetic light curves and PDS, Bursa et al. (2004) performed general relativistic ray tracing of a torus undergoing simultaneous vertical and radial eigenmode oscillations with a 3:2 frequency ratio, finding that the (higher fre-

quency) vertical oscillation appreciably modulates the light curve only in GR, and not in Newtonian gravity. Schnittman & Rezzolla (2006) performed a similar study of the emission from a numerically simulated thick torus responding to a radial perturbation. Schnittman & Bertschinger (2004) considered a hotspot radiating isotropically on nearly circular equatorial orbits around a Kerr black hole and also performed general relativistic ray tracing to produce synthetic light curves and PDS, finding that the hotspots would have to be very elongated in the azimuthal direction to suppress the (unobserved) harmonic content. The hotspot model has also been investigated by Bakala et al. (2014). The Rossby wave instability may also lead to the appearance of QPOs in black hole binaries (Tagger & Varnière 2006; Vincent et al. 2013).

Geometrically thick discs or tori (Abramowicz, Jaroszynski & Sikora 1978; Abramowicz & Fragile 2013), as well as slender tori (Blaes 1985; Abramowicz et al. 2006), have been of great interest in studies of QPOs using analytical as well as numerical simulations. Rezzolla et al. (2003) studied high-frequency QPOs in oscillating tori in the Newtonian as well as GR framework. Kluźniak (2005), Abramowicz et al. (2006) and Blaes et al. (2006) analytically calculated linear modes in relativistic slender hydrodynamical tori and also suggested that these modes can be excited in numerical simulations. Straub & Šrámková (2009) extended these results to non-slender tori. Mazur et al. (2013) and Vincent et al. (2014) studied analytic tori subjected to radial, vertical, shear and expansion perturbations using relativistic ray tracing.

Ray tracing of relativistic hydrodynamical simulations of an oscillating torus was performed for the first time by Schnittman & Rezzolla (2006) for radial velocity perturbations. They found multiple high-frequency peaks in the simulated light curve, and in this paper we are going to revisit the frequency identification with particular modes that were suggested by the authors. In order to generalize this hydrodynamical work, Schnittman, Krolik & Hawley (2006) used relativistic ray tracing to compute the light curves from 3D global relativistic MHD simulations of accreting tori around black hole. However, these authors did not see any convincing evidence of QPOs in their MHD simulations.

This paper aims to further develop the study of oscillating tori as a model for high-frequency QPOs and in particular twin-peak QPOs. To this end, we perturb the 4-velocity of equilibrium tori, use relativistic hydrodynamical simulations to compute their evolution and relativistic ray tracing to integrate the radiative transfer equation and produce PDS. Such computations have been done in the past by Schnittman & Rezzolla (2006). We perform a study of ray-traced light curves of oscillating relativistic tori in the Schwarzschild metric, extending the analysis to more general perturbations (including vertical motions of the torus), and focusing on geometrically slender tori. In this work, we were able to identify five prominent QPO peaks in the PDS of the light curves of the numerically simulated tori, corresponding to the lowest frequency, axisymmetric eigenmodes, which are well understood analytically (Kluźniak 2005; Abramowicz et al. 2006; Blaes et al. 2006). One of the main limitations of this study is the fact that no self-consistent model is given for the origin of the perturbation of the equilibrium torus. This point may be addressed in future work.

The paper is organized as follows. In Section 2, we present the numerical methods governing the general relativistic hydrodynamical (GRHD) equations to setup the torus simulation and perform relativistic ray tracing of the results. In Section 3, we discuss the simulated light curves and the numerical artefacts found within them. In Section 4, we give results inferred from PDS extracted from the light curves. In Section 5, we summarize our findings.

<footnote n="1">¹ However, the observed *Q* value may not be intrinsic, owing to the time of integration of many thousands of oscillation periods required to produce a PDS with a statistically significant high-frequency QPO.</footnote>
<table><caption>**Table 1.** Initial torus setup. The last column gives the oscillation modes that dominate the power spectra (see Section 3). The radial, vertical, plus, breathing and ‘ex’ modes are listed in short form as R, V, +, B and X.</caption><thead><tr><th>Perturbation</th><th>R<sub>in</sub>/M</th><th>R<sub>c</sub>/M</th><th>R<sub>cross</sub>/R<sub>c</sub></th><th>Detected modes</th></tr></thead><tbody><tr><td>Radial</td><td>9.5</td><td>10.0</td><td>0.05</td><td>R, +, B</td></tr><tr><td>Radial</td><td>14.5</td><td>15.0</td><td>0.03</td><td>R, +, B</td></tr><tr><td>Vertical</td><td>9.5</td><td>10.0</td><td>0.05</td><td>V, X, B</td></tr><tr><td>Vertical</td><td>14.5</td><td>15.0</td><td>0.03</td><td>V, X, B</td></tr><tr><td>Diagonal</td><td>9.5</td><td>10.0</td><td>0.05</td><td>R, V, +, X, B</td></tr><tr><td>Diagonal</td><td>14.5</td><td>15.0</td><td>0.03</td><td>R, V, +, X, B</td></tr></tbody></table>

Additionally, various tests to see effects of the background density on torus oscillations, the mode identifications and a comparison of our thick torus simulations with the existing hydro simulations are presented in Appendices A, B and C, respectively. We use $G = 1 = c$ throughout the paper.

# 2 NUMERICAL SETUP AND QPO MODES

## 2.1 GRHD simulation

We perform only hydro simulations in this work, although recently Mazur et al. (2016) performed radiative simulations to study QPOs. In this paper, we numerically solve the equations of relativistic hydrodynamics derived from the following stress-energy tensor for the torus fluid

$$T_{\mu\nu} = (\rho + \rho\epsilon + P)u_{\mu}u_{\nu} + Pg_{\mu\nu}, \quad (1)$$

where $\epsilon$ is the specific internal energy, $\rho$ is the rest-mass density, $P$ is the fluid pressure, $u^\mu = g^{\mu\nu}u_\nu$ is the contravariant component of fluid 4-velocity and $g_{\mu\nu}$ is the Schwarzschild metric tensor. We use the spherical Kerr–Schild coordinate system in the GRHD simulations. From the conservation of rest mass and energy-momentum, we arrive at the following conservation equations that we evolve numerically using the general relativistic MHD code COSMOS++ (Anninos, Fragile & Salmonson 2005),

$$\partial_t D + \partial_i (DV^i) = 0, \quad (2)$$

$$\partial_t \mathcal{E} + \partial_t (-\sqrt{-g} T_j^i) = -\sqrt{-g} T_\lambda^k \Gamma^\lambda_{0k}, \quad (3)$$

$$\partial_t S_j + \partial_i (\sqrt{-g} T_j^i) = \sqrt{-g} T_\lambda^k \Gamma^\lambda_{jk}. \quad (4)$$

where $D = W\rho$ is the generalized fluid density, $W = \sqrt{-g}u^i$ is the generalized boost factor, $V^i = u^i/u^i$ is the transport 3-velocity, $g$ is the metric determinant,

$$\mathcal{E} = -\sqrt{-g} T_0^0 = -W u_t \rho h - \sqrt{-g} P \quad (5)$$

is the total energy density, $h = 1 + \epsilon + P/\rho$ is the specific enthalpy and

$$S_v = \sqrt{-g} T_j^0 = W u_j \rho h \quad (6)$$

is the covariant momentum density. $\Gamma$ is the Christoffel symbol describing metric connection coefficients. We solve the GRHD equations using the high-resolution shock-capturing option described in Fragile et al. (2012) and Dibi et al. (2012).

We choose a stationary hydrodynamical torus of cross-sectional radius, $R_{cross} = (R_c - R_{in}) = 0.5 M$ around an $M = 10 M_\odot$ black hole. We accomplish this by setting the inner and central (density and pressure maximum) radii of the torus for each model, as given in Table 1. Initially, we chose $R_c = 10 M$. However, in

order to distinguish between the nearly degenerate vertical and plus-mode frequency at $R_c = 10 M$, we also considered a second central radius at $R_c = 15 M$. In both cases, the inner radius of the torus is set to $R_{in} = R_c - 0.5$. All the simulations in our study assume a constant specific angular momentum distribution, which is defined as $l = -u_\phi/u_t$, where $u_\phi$ and $u_t$ are the covariant azimuthal and time components of the 4-velocity, respectively. We solve for the initial internal energy distribution of the torus $\epsilon(r, \theta)$. To initialize the setup, we assume a polytropic equation of state $P = \rho\epsilon(\gamma - 1) = \kappa\rho^\gamma$, so the density is given by $\rho = [\epsilon(\gamma - 1)/\kappa]^{1/(\gamma - 1)}$, with $\gamma = 5/3$. We also use $\gamma = 4/3$ to reproduce model c of Zanotti, Rezzolla & Font (2003) (Appendix C). We can choose $\kappa$ to set the density (and mass) normalization of the initial torus. Once the initial setup is fixed, the gas temperature of the torus is fixed by the gravity. The torus is evolved with time using ideal gas equation of state $P = (\gamma - 1)\rho\epsilon$.

The initial torus setup is subjected to a vertical or radial velocity perturbation using velocity fields corresponding to the vertical and radial eigenmodes, respectively, in the slender torus limit (Bursa et al. 2004; Blaes et al. 2006). However, we note that our tori do not strictly satisfy this limit, as they are extended solutions with eigenmodes not expected to correspond to uniform velocity fields (Blaes et al. 2007). The diagonal perturbation is obtained by adding these two uniform velocity fields. The equilibrium torus is defined by the covariant 4-velocity,

$$u_\mu = (u_t, 0, 0, u_\phi). \quad (7)$$

After applying the velocity perturbation, the covariant 4-velocity is given by

$$\mathbf{u} = \begin{cases} (u_t, \eta, 0, u_\phi), & \text{radial} \\ (u_t, 0, K\eta, u_\phi), & \text{vertical}, \\ (u_t, \eta, K\eta, u_\phi), & \text{diagonal} \end{cases} \quad (8)$$

where

$$K = -\sqrt{\frac{g_{\theta\theta}(R_c)}{g_{RR}(R_c)}}. \quad (9)$$

This parameter is chosen such that radial and vertical 4-velocity perturbations will have the same magnitude. The parameter $\eta$ normalizes the perturbation. We choose $\eta = 0.003 V_{\text{Kepler}}$, where $V_{\text{Kepler}} = 2\pi R_c v_K$ is the Keplerian speed at the centre of torus, with $v_K$ being the orbital frequency. Illustrations of the various velocity perturbations are shown in Fig. 1.

The various initial setups we study are given in Table 1. The computational domain along radial direction is defined by $R_{\text{min}} = (R_c - 0.75) M$ and $R_{\text{max}} = (R_c + 0.75) M$. The domain along the polar direction is defined by $\theta_{\text{min}} = 86.4^\circ$ and $\theta_{\text{max}} = 93.6^\circ$ for the torus at $R_c = 10 M$ and $\theta_{\text{min}} = 87.66^\circ$ and $\theta_{\text{max}} = 92.34^\circ$ for the torus with $R_c = 15 M$. We performed the hydrodynamics simulations with a linearly spaced grid resolution, $n_r \times n_\theta = 256 \times 256$ along both radial and polar directions. We also checked one simulation with $512 \times 512$ to confirm that the results we obtained are not sensitive to the chosen resolution. All the models presented in Table 1 are computed for a total coordinate time of $\Delta t = 20 T_K$ (orbital period at the centre of the torus), where $T_K = 198.6 M$ for $R_c = 10 M$ and $T_K = 365.0 M$ for $R_c = 15 M$.

We also performed various test runs to check the robustness of our results, and carried out several additional simulations to clarify the results and relate them to previous work. These are summarized in Figs C1–C3 in Appendix C. Different background densities and their effects are explored in Appendix A.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Three initial velocity perturbations, equation (8) to the torus. Plots only present the initial torus with pressure maximum at $R_c = 15 M$. The torus setup at $R_c = 10 M$ looks qualitatively similar at corresponding locations within the torus.</figcaption></figure>

<table><caption>**Table 2.** Analytically predicted frequencies (using $\gamma = 5/3$ in Section 2.2) of the five lowest order oscillation modes of slender tori for the two central radii considered in this study, in units of the Keplerian frequency at the centre of the torus, $v_K$.</caption><thead><tr><th>R<sub>c</sub>/M</th><th>Radial (R)</th><th>Vertical (V)</th><th>X</th><th>Plus (+)</th><th>Breathing (B)</th></tr></thead><tbody><tr><td>10</td><td>0.63</td><td>1</td><td>1.18</td><td>0.98</td><td>1.66</td></tr><tr><td>15</td><td>0.77</td><td>1</td><td>1.26</td><td>1.18</td><td>1.69</td></tr></tbody></table>

## 2.2 Analytic eigenfrequencies of QPO modes

The analytical expressions of the various torus eigenmodes that we are interested have been calculated in Kluźniak (2005), Abramowicz et al. (2006), Blaes et al. (2006) and Straub & Šrámková (2009). These are the radial mode (R), the vertical (θ or V) mode, the X mode, the plus mode (+) and the breathing mode (B). We shall follow the notations of Blaes et al. (2006) in defining the frequencies of the eigenmodes in the slender torus limit. These expressions are valid for a polytropic equation of state; in this work, we assume a polytropic index of $\gamma = 5/3$, except for the simulations reported in Appendices A, B and C, which use $\gamma = 4/3$.

The lowest order eigenmode frequencies of generalized polytropic index tori with pressure maxima at $R_c$ are

$$ \tilde{\sigma}_R^2 = \tilde{\omega}_R^2 = \left(1 - \frac{6}{R_c}\right), \qquad (10) $$

$$ \tilde{\sigma}_{\theta}^2 = \tilde{\omega}_{\theta}^2 = 1, \qquad (11) $$

$$ \tilde{\sigma}_{X}^{2} = \sqrt{\tilde{\omega}_{R}^{2} + \tilde{\omega}_{\theta}^{2}}, \qquad (12) $$

$$ \tilde{\sigma}_{+}^{2} = \frac{1}{2n} [a - b^{1/2}], \qquad (13) $$

$$ \tilde{\sigma}_{B}^{2} = \frac{1}{2n} [a + b^{1/2}], \qquad (14) $$

where, $\tilde{\omega}_R$, $\tilde{\omega}_\theta$ are radial and vertical epicyclic frequencies scaled by the Keplerian orbital frequency $2\pi/T_K$ at the centre of the torus. The a and b are defined by

$$ a = (2n + 1) (\tilde{\omega}_{\theta}^2 + \tilde{\omega}_{R}^2), \qquad (15) $$

$$ b = [4n(n + 1)(\tilde{\omega}_{\theta}^2 - \tilde{\omega}_{R}^2)^2 + (\tilde{\omega}_{\theta}^2 + \tilde{\omega}_{R}^2)^2], \qquad (16) $$

where, $n = 1/(\gamma - 1)$. Table 2 shows the predicted values of these frequencies for $\gamma = 5/3$ in above expressions. The numerical values of mode frequencies identified in our numerical simulations

<table><caption>**Table 3.** Summary of mode frequencies observed in Figs 7 and 8. The ‘-’ denotes absent modes. In the model names, R,V and D correspond to radial, vertical and diagonal perturbations, respectively. The ‘Thick’ and ‘Thin’ refer to the optical depth. The values with * correspond to degenerate frequencies for the vertical and plus modes at $R_c = 10$.</caption><thead><tr><th>Model</th><th>R</th><th>V</th><th>X</th><th>+</th><th>B</th><th>+/R (Table 2)</th><th>+/R</th><th>B/V</th></tr></thead><tbody><tr><td>R10Thick</td><td>0.628</td><td>-</td><td>-</td><td>0.966</td><td>1.64</td><td>1.55</td><td>1.54</td><td>-</td></tr><tr><td>R10Thin</td><td>0.628</td><td>-</td><td>-</td><td>0.960</td><td>1.64</td><td>1.55</td><td>1.53</td><td>-</td></tr><tr><td>V10Thick</td><td>-</td><td>0.989</td><td>1.17</td><td>-</td><td>1.63</td><td>-</td><td>-</td><td>1.65</td></tr><tr><td>V10Thin</td><td>-</td><td>0.987</td><td>1.20</td><td>-</td><td>1.64</td><td>-</td><td>-</td><td>1.66</td></tr><tr><td>D10Thick</td><td>0.628</td><td>0.989*</td><td>1.17</td><td>0.989*</td><td>1.62</td><td>1.55</td><td>1.57</td><td>1.64</td></tr><tr><td>D10Thin</td><td>0.628</td><td>0.994*</td><td>1.20</td><td>0.994*</td><td>1.64</td><td>1.55</td><td>1.58</td><td>1.65</td></tr><tr><td>R15Thick</td><td>0.772</td><td>-</td><td>-</td><td>1.18</td><td>1.68</td><td>1.52</td><td>1.53</td><td>-</td></tr><tr><td>R15Thin</td><td>0.773</td><td>-</td><td>-</td><td>1.16</td><td>1.68</td><td>1.52</td><td>1.50</td><td>-</td></tr><tr><td>V15Thick</td><td>-</td><td>0.997</td><td>1.26</td><td>-</td><td>1.69</td><td>-</td><td>-</td><td>1.69</td></tr><tr><td>V15Thin</td><td>-</td><td>0.996</td><td>-</td><td>-</td><td>1.70</td><td>-</td><td>-</td><td>1.70</td></tr><tr><td>D15Thick</td><td>0.771</td><td>0.996</td><td>1.25</td><td>1.16</td><td>1.68</td><td>1.52</td><td>1.50</td><td>1.69</td></tr><tr><td>D15Thin</td><td>0.773</td><td>0.998</td><td>-</td><td>1.16</td><td>1.68</td><td>1.52</td><td>1.50</td><td>1.68</td></tr></tbody></table>

for $\gamma = 5/3$ are reported in Table 3. We found that both sets of frequencies are in good agreement.

## 2.3 Ray tracing

The numerical data obtained from the `cosmos++` code are ray traced using the open-source code,<fnref n="2" /> `gyoto` (Vincent et al. 2011). An observer is placed at a distance, $R_{obs} = 10^4 M$, from the black hole. Photons of energy $E_{obs}$ are traced backwards in time by integrating the geodesic equation in Schwarzschild geometry. When the photon reaches the GRHD simulation box, two different radiation mechanisms are considered. The optically thick case assumes that blackbody radiation is emitted at the last scattering surface of the torus. Here scattering refers to Thomson scattering. We thus integrate the Thomson optical depth

$$ d\tau_e = 0.4 \rho \, ds, \qquad (17) $$

where the rest-mass density $\rho$, and the proper length in the emitter's frame, $s$, are expressed in cgs units. The locus of points where $\tau_e$ reaches unity is defined as the torus last scattering surface. Photons emitted from this surface will experience no Thomson scattering, and their trajectories are thus simply defined by gravitation

<footnote n="2">² Freely available at gyoto.obspm.fr</footnote>
(they follow the null geodesics of the metric). The blackbody radiation is emitted from this surface at the local temperature (assuming a perfect-gas law, so that $T \propto P/\rho$) and photon energy $E_{\text{em}}$ (which is of course different from $E_{\text{obs}}$ due to relativistic effects). The optically thin case assumes that bremsstrahlung radiation is emitted by the torus. Here, the radiative transfer equation is integrated using the well-known emission and absorption coefficients for bremsstrahlung radiation (see e.g. Rybicki & Lightman 1986).

The GRHD simulations are scale free in density, meaning that the mass density obtained from one given solution can be multiplied by an arbitrary constant factor to get another well-defined solution. As a consequence, we are free to choose the mass-density normalization. We choose it by imposing that the temperature at the last scattering surface be of the order of $10^7$ K, which fixes the density (typically, $\rho_c = 10^{-5}-10^{-6}$ g cm$^{-3}$). Because of the way these two emission models are constructed, the optically thick case mostly traces the evolution of the torus surface (of last scattering), while the optically thin case traces the variation of its full volume but weighted towards the densest regions.

For both the optically thick and thin cases, the torus 3-velocity and density are interpolated trilinearly in 2D space (R, θ) and time. When the density is evolving too abruptly between adjacent cells to allow interpolation in this simple way (which is the case close to the surface of the torus where the density gradient becomes quite steep), the closest neighbouring grid zone is used.

## 3 LIGHT CURVES

By creating and combining ray-traced images that span the entire duration of the simulations, we can create synthetic light curves. In doing so, we have noticed two numerical artefacts, that we mainly attribute to the hydro simulations. In this section, we investigate these artefacts in more detail. First, we describe increasing and decreasing secular trends in the computed light curves. Secondly, we address noisy behaviour in the light curves. We also perform numerical experiments within our post-processing routines to explore ways to remove the noise.

### 3.1 Trends in light curves

The raw, computed light curves in Fig. 2 show examples of increasing (left-hand panel) and decreasing (right-hand pane) secular trends. These trends amount to about a 1 per cent increase or decrease in the flux over the duration of our simulations. In order to investigate these trends further, we consider the $R_c = 15 M$ case, as the trends are more pronounced in it. The left-hand panel of Fig. 2 is for the vertically perturbed, optically thin case, while the right-hand panel shows the optically thick light curve for the same simulation. These trends seem to be related to the growth of the total mass in the central regions of the torus in this simulation. In Fig. 3, we show the scaled total mass in the central parts of the tori by defining two central regions, $\rho > 99$ per cent $\rho_{\text{max}}$ and $\rho > 85$ per cent $\rho_{\text{max}}$. The dotted lines show fits to the underlying secular trends. The slope for the $R_c = 15$ case in the region with $\rho > 0.99 \rho_{\text{max}}$ is 0.009, while for $R_c = 10$, this slope is 0.001. These secular trends in the total mass explain both trends in the light curves, as well as why the light-curve trends are more pronounced in the $R_c = 15$ case. As the central density increases, the temperature of the gas will decrease (as $T \propto P/\rho$). This explains why the optically thick light curve decreases with time, as the assumed blackbody radiation will weaken with lower temperatures. However, the optically thin light curve will show the opposite trend, because the bremsstrahlung emission

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Raw light curves for the vertically perturbed torus at $R_c = 15 M$ and $\gamma = 5/3$ at an inclination angle of $i = 10^\circ$. The horizontal axis shows time in units of Keplerian orbital period, $T_K$, at pressure maximum of torus. The left-hand panel shows the optically thin case, while the right one shows the optically thick.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Secular trend of the total mass in the very central region of the torus, defined by $\rho > 0.99 \rho_{\text{max}}$ (solid) and $\rho > 0.85 \rho_{\text{max}}$ (dashed) for $R_c = 10 M$ (red, bottom) and $R_c = 15 M$ (black, top) using $\gamma = 5/3$. The normalization of this plot is arbitrary and is chosen simply for viewing convenience. The dotted lines are fits to the slopes of the underlying trends. The slope for the $R_c = 15$ torus is nine times larger than the slope for the $R_c = 10$ torus.</figcaption></figure>

coefficient in our range of parameters varies as $\propto \rho^2$. The secular trends of our light curves are thus consistent with the underlying GRHD simulations.

For the rest of this paper, all secular trends are removed from our synthetic light curves. We have confirmed that this does not change the dominant peaks in the PDS, and thus, does not affect our conclusions. However, by removing these trends, some of the low-frequency noise goes away. The trend-free light curves from each of our simulations and each of our radiation models are shown in Figs 4 and 5, for an observer inclination of $i = 10^\circ$ (nearly face-on view). For all cases, the maximum relative flux variation is of the order of 1 per cent. This is in reasonable agreement with X-ray observations of QPOs in black hole systems.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Trend-free light curves of the oscillating torus with $R_c = 10 M$ and $\gamma = 5/3$, assuming an optically thin (left) or optically thick (right) emission, for three initial velocity perturbations, radial (top), vertical (middle) and diagonal (bottom). The percentage values in each panel correspond to the maximum relative variation of the flux.</figcaption></figure>

## 3.2 Noise in light curves

Some of the light curves presented in Figs 4 and 5 exhibit rather noisy behaviour (the clearest cases are the middle right panel of Fig. 4 and the top right panel of Fig. 5). In this section, we explore the source of this noise via three numerical tests. We focus primarily on the radially perturbed, $R_c = 15 M$, optically thick case, as it exhibited the greatest amount of noise. In the first test (upper left panel of Fig. 6), we contrast the light curve from Fig. 5 with a new light curve resulting from increasing the torus density normalization by 33 per cent. The resulting light curve is considerably less noisy. The main effect of increasing the density normalization in the optically thick case is to enlarge the surface of last scattering. Making this surface larger effectively improves our ability to resolve it. There are also some subtle differences in the nature of the torus motion further away from its centre. However, an important point is that the PDS remains largely unchanged between the two cases (upper right panel of Fig. 6).

The middle panel of Fig. 6 investigates how the light curve varies when only a small portion of the torus is considered. Here, the reference light curve is compared to a light curve computed in such a way that only a small area of the torus, close to the maximum relativistic beaming region, contributes to the observed flux. It is interesting to note that this 'zoomed in' light curve is much smoother than the reference light curve. This fact can be easily explained: light from different parts of the torus reach the observer at different times

due to the finite propagation speed of light and light-bending effects. Each small region of the torus emits a rather smooth signal, but the varying arrival times of these smooth signals of different intensities can lead to a noisy light curve. Interestingly, the power in the plus mode is reduced in this case. This is because the plus mode involves relative motion of the inner and outer parts of the torus, unlike the dominant motion in the R (or V) mode.

Finally, the lower panel of Fig. 6 shows a light curve generated using the same mass-density scaling as in the reference simulation, but dividing the value of the temperature at the last scattering surface by a factor of 100. Note that when we change the density normalization or the temperature, the PDS displays the same dominant peaks with the possible exception of the harmonics. The amplitude of the variation is linked to how much the Planck function varies when the temperature is varied by some amount. But this effect is strongly dependent on the value of the initial temperature. For instance, at an observer energy of 1 keV, if $T$ varies by 1 per cent around $T = 10^6$ K, the Wien tail varies by $\approx 10$ per cent, whereas, if $T = 10^8$ K, the variation of the blackbody emission is reduced to $\approx 1$ per cent because we are then in the Rayleigh-Jeans part of the spectra. This experiment thus shows the importance of the absolute value of the temperature at the last scattering surface.

These three numerical tests highlight the high level of complexity embedded in the light curves of oscillating tori. However, it is important to note that the PDS in all our tests always display the
<figure><img src="image_7.png" /><figcaption>Figure 5. Same as Fig. 4 with $R_c = 15 M$.</figcaption></figure>

same dominant peaks. This assures us that such subtleties do not affect the conclusions and findings of this paper.

# 4 RESULTS

## 4.1 Power density spectra

The Lomb–Scargle PDS (Press & Rybicki 1989) of all trend-free light curves (two different radial positions with two different inclinations) are shown in Figs 7 and 8. The first obvious feature appearing in the PDS is that all of them are dominated by peaks corresponding to the mode(s) of initial perturbation. More interesting is the presence of additional peaks in all PDS. Given that the diagonal perturbation PDS appear simply as the sum of the radial and vertical PDS, we shall mainly focus on these latter ones in our discussions.

In the PDS of the radially perturbed torus in Figs 7 and 8, the prominent peaks correspond to the radial, plus and breathing modes, as well as the sum of the radial and plus modes and higher harmonics of the radial mode. We checked carefully that the presence of the plus and breathing modes is not a numerical artefact due to interaction of the torus with its surrounding background. We did this by exploring different background densities, as described in Appendix A.

In the PDS of the vertically perturbed $R_c = 10 M$ torus in Fig. 7, the vertical and breathing modes are present in both the optically thin and thick cases. In addition, the optically thick case shows a

prominent peak at the X-mode frequency, though it is much less prominent in the optically thin case. This is due to the fact that in the optically thick case we are tracing the surface, and in the optically thin case the volume. The X mode, because of its nature, causes more fluctuations on the surface than in the volume. Fig. 8 ($R_c = 10 M$) has similar features as Fig. 7, except that we see almost no evidence for the X mode in left (optically thin) panel of Fig. 8.

Three higher frequency harmonics and combinations of frequencies also appear with strong power in both tori. They correspond to twice the vertical frequency, twice the X-mode frequency and the sum of the vertical and X-mode frequencies. We note that multiplying two sinusoidal functions gives rise to a signal varying at the sum of their frequencies. The ‘vertical + X’ peak could thus be interpreted as a non-linear coupling between these two modes. Pseudo-Newtonian hydrodynamic studies have also found a similar eigenmode content in corresponding simulations (Parthasarathy, Manousakis & Kluźniak 2016).

In all our PDS except for the optically thick, vertically perturbed case, there are clearly dominant modes corresponding to the initial velocity perturbations, with all other oscillation modes appearing significantly weaker. Only in the optically thick, vertically perturbed cases are the ‘secondary’ peaks of comparable strength to the vertical epicyclic peaks. We note that the velocity eigenvector field for a vertically oscillating, non-slender torus, contains within it the velocity field reminiscent of an X mode (see fig. 5 of Blaes et al. 2007). As illustrated in Fig. 9, in a non-slender torus, the vertical and radial modes do not correspond to uniform motion. Therefore, it seems
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Three ways of smoothing the light curves (the reference case, in black) of the radially perturbed, $R_c = 15M$, optically thick simulation. The top panel shows the effect of increasing the mass-density normalization by a factor of 33 per cent. The middle panel shows the light curve emitted by a small region of the torus, located close to the maximum relativistic beaming region. The bottom panel shows the effect of changing the absolute value of temperature at the last scattering surface. See the text for more details. The right-hand panel illustrates the corresponding changes to the PDS. Note that with the exception of the plus-mode peak in the middle panel, the only major effect is on the harmonics at $v > 2v_K$.</figcaption></figure>

most likely that the prominent X mode is a result of the initial perturbations and not evidence of non-linear coupling.

Another observation that comes from our data is that the breathing mode is present in all our reported simulations. A high contrast between torus central density and background density can enhance the breathing mode. However, the breathing mode is always present even if we choose the lowest possible background density, and its amplitude converges to some value for low background densities (see Appendix A). This strongly suggests that the excitation of the breathing mode in all our simulations is a consequence of the adopted velocity field perturbation, and not a numerical artefact.

Finally, we do not note a strong inclination dependence in our optically thin PDS, unlike in Bursa et al. (2004), whose findings were based on an incompressible torus. However, the inclination dependence is a very strong function of the torus equation of state (Bursa, private communication). The inclination dependence reported in Schnittman & Rezzolla (2006) is also obtained for a different equation of state ($\gamma = 4/3$), and for a much thicker torus.

## 4.2 Observational features

To recapitulate, the most important conclusion of our analysis is that the dominant modes appearing in the PDS are the radial, plus

and breathing modes for a radial perturbation, and the vertical, X and breathing modes for a vertical perturbation (Table 3). Note that the numerical values for each mode frequency are taken as the centroid of a Gaussian fit of width $\Delta v = 0.2$ around the peak in the PDS. We also find that certain mode combinations occur in a nearly 3:2 frequency ratios. For example, in all the radial and diagonal perturbation cases, the ratio of the plus to radial mode frequencies is $\approx 1.5$ (Table 3, column 7), in agreement with earlier analytic (Fig. 3, Blaes et al. 2006) and numerical (Schnittman & Rezzolla 2006) studies. For the $R_c = 15M$ diagonal perturbation case, the mode frequency ratio 1.50 is closer to 3:2 compared to the one inferred from analytic calculations of non-slender tori (Straub & Šrámková 2009). For the vertical perturbations, we find the ratio of the breathing to vertical mode frequencies to be $\approx 1.7$ (Table 3, column 8), consistent with a $\gamma = 5/3$ gas (Blaes et al. 2006).

## 5 CONCLUSIONS

In total, we report 22 hydrodynamic simulations of oscillating tori in the Schwarzchild metric. We ray trace the results (only the ones reported in Figs 4 and 5) to obtain PDS. We concentrated on three types of velocity perturbations: radial, vertical and diagonal. The important conclusions of the paper concern the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Power density spectra (PDS), computed from the trend-free light curves (Fig. 4) of the oscillating, $R_c = 10$ M torus with radial, vertical and diagonal velocity perturbations (top to bottom). The horizontal axis corresponds to the frequency normalized by the Keplerian frequency at the torus centre $v_K$. The solid black and dashed red curves show the PDS for two inclinations, $i = 10^\circ$ and $i = 70^\circ$, respectively. The left- and right-hand panel correspond to optically thin and thick cases, respectively. The vertical grey lines are located at the values of frequencies of the five lowest order oscillation modes, as given in Table 2.</figcaption></figure>

dominant modes, the possibility of non-linear coupling, an absence of inclination dependence in optically thin cases and the distribution of mode power for different radiative transfer models. To summarize:

(i) We found that a uniform radial perturbation of the velocity field can trigger both the radial and plus modes. Likewise, a uniform vertical perturbation can trigger the vertical and X modes. We also found that the response appears to be independent of where we place the pressure maximum of the torus. We argue that the presence of both the plus mode for the radial perturbation and the X mode for the vertical perturbation is owing to the non-uniformity of the velocity associated with the vertical and radial eigenmodes of non-slender tori. This point is illustrated in Fig. 9.

(ii) We did not find any clear dependence on inclination in our ray-traced light curves for the optically thin cases, in contrast to previous results in the literature. We suspect that this is due to the fact that our tori are very thin and compressible.

(iii) We found that the two different radiative models we considered can give significantly different light curves and hence different PDS. In particular, the optically thin and thick cases give different powers in the plus and X modes.

(iv) We have also found that our hydrodynamical simulations of non-slender tori show a frequency ratio close to 3:2 for the plus to radial mode. The PDS exhibit ratios slightly closer to 3:2 than predicted by analytic calculations.

(v) We also assess two numerical artefacts found in our light curves. The first concerns secular (increasing/decreasing) trends in the overall flux, while the second deals with noise in the light curves. We checked that these artefacts do not affect our conclusions.

(vi) In addition, owing to the small cross-sectional radius of the torus, we also see evidence of the breathing mode in all our simulations. We confirmed its appearance by performing various low- and high-density background setups. In all such setups, the breathing mode is present but is more prominent for larger background densities (Appendix A). At lower background densities ($b_0 = 10^{-6}$, $10^{-7}$, $10^{-8}$), the corresponding PDS look similar. This ensures that breathing mode is not a numerical artefact.

(vii) We also shed some light on the issue of mode identification. Our simulations suggest that the $f$ and $o_1$ eigenmodes in vertically integrated models (fig. 2 of Rezzolla et al. 2003) correspond to our radial and plus modes. A detailed discussion is presented in Appendix B.

(viii) We also closely reproduced model c of Zanotti et al. (2003). We found that the *fundamental* (f) eigenfrequency, $v \approx 52$ Hz, is independent of the equation of state of the gas (Fig. C2). Conversely, the plus-mode eigenfrequency shows a clear dependence on the chosen equation of state. The plus mode shows an approximately 13 per cent increase in eigenfrequency for $\gamma = 4/3$ compared to $\gamma = 5/3$ (Fig. C2).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Same as Fig. 7, but for the $R_c = 15 M$ torus.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Qualitative illustration of eigenmodes of a non-slender torus. See Parthasarathy, Manousakis & Kluźniak (2015) for snapshots of a simulation.</figcaption></figure>

## ACKNOWLEDGEMENTS

We thank Jiří Horak, Omer Blaes, Michal Bursa and Odele Straub for discussions. It is a pleasure to thank Olindo Zanotti for help in the comparison of our results. MB also thank to College of Charleston SC for hosting during final stage of project. We also thank the referee, Luciano Rezzolla for useful suggestions. Research was supported by the Polish NCN grant 2013/08/A/ST9/00795. Part of the work was supported by the Polish NCN grants 2013/09/B/ST9/00060, 2014/15/N/ST9/04633. The simulations were performed using the PROMETHEUS supercomputer in the PL-Grid infrastructure.

## REFERENCES

* Abramowicz M. A., Fragile P. C., 2013, Living Rev. Relativ., 16, 1

* Abramowicz M. A., Kluźniak W., 2001, A&A, 374, L19

* Abramowicz M., Jaroszynski M., Sikora M., 1978, A&A, 63, 221

* Abramowicz M. A., Blaes O. M., Horák J., Kluzniak W., Rebusco P., 2006, Class. Quantum Gravity, 23, 1689

* Altamirano D., Belloni T., 2012, ApJ, 747, L4

* Anninos P., Fragile P. C., Salmonson J. D., 2005, ApJ, 635, 723

* Bakala P. et al., 2014, MNRAS, 439, 1933

* Barret D., Kluźniak W., Olive J. F., Paltani S., Skinner G. K., 2005, MNRAS, 357, 1288

* Belloni T., Méndez M., Sánchez-Fernández C., 2001, A&A, 372, 551

* Belloni T. M., Sanna A., Méndez M., 2012, MNRAS, 426, 1701

* Blaes O. M., 1985, MNRAS, 216, 553

* Blaes O. M., Arras P., Fragile P. C., 2006, MNRAS, 369, 1235

* Blaes O. M., Šrámková E., Abramowicz M. A., Kluźniak W., Torkelsson U., 2007, ApJ, 665, 642

* Bursa M., Adámek K., 2014, in Stuchlík Z., Török G., Pechaček T., eds, Proceedings of Ragtime 14–16: Workshops on Black Holes and Neutron Stars. Silesian University, Opava, p. 31

* Bursa M., Abramowicz M. A., Karas V., Kluźniak W., 2004, ApJ, 617, L45

* Cui W., Zhang S. N., Chen W., 1998, ApJ, 492, L53

* Dexter J., Blaes O., 2014, MNRAS, 438, 3352

* Dibi S., Drappeau S., Fragile P. C., Markoff S., Dexter J., 2012, MNRAS, 426, 1928

* Fragile P. C., Gillespie A., Monahan T., Rodriguez M., Anninos P., 2012, ApJS, 201, 9

* Fragile P. C., Straub O., Blaes O., 2016, MNRAS, 461, 1356

* Johannsen T., Psaltis D., 2011, ApJ, 726, 11

* Kato S., 2001, PASJ, 53, 1
Kluźniak W., 1998, ApJ, 509, L37

Kluźniak W., 2005, Astron. Nachr., 326, 820

Kluźniak W., Abramowicz M. A., 2001, Acta Phys. Pol. B, 32, 3605

Kluzniak W., Michelson P., Wagoner R. V., 1990, ApJ, 358, 538

Kluźniak W., Abramowicz M. A., Lee W. H., 2004, in Kaaret P., Lamb F. K., Swank J. H., eds, AIP Conf. Ser., Vol. 714, X-ray Timing 2003: Rossi and Beyond. Am. Inst. Phys., New York, p. 379

Kluźniak W., Abramowicz M. A., Kato S., Lee W. H., Stergioulas N., 2004, ApJ, 603, L89

Lee W. H., Abramowicz M. A., Kluźniak W., 2004, ApJ, 603, L93

Machida M., Matsumoto R., 2008, PASJ, 60, 613

Mazur G. P., Vincent F. H., Johansson M., Šramková E., Török G., Bakala P., Abramowicz M. A., 2013, A&A, 554, A57

Mazur G. P., Zanotti O., Sądowski A., Mishra B., Kluźniak W., 2016, MNRAS, 456, 3245

Miller J. M. et al., 2001, ApJ, 563, 928

Miranda R., Horák J., Lai D., 2015, MNRAS, 446, 240

Montero P. J., Rezzolla L., Yoshida S., 2004, MNRAS, 354, 1040

Morgan E. H., Remillard R. A., Greiner J., 1997, ApJ, 482, 993

Parthasarathy V., Manousakis A., Kluźniak W., 2015, Proc. Polish Astron. Soc., Vol. 3, Numerical Simulations of Oscillating Tori in a Pseudo-Newtonian Potential

Parthasarathy V., Manousakis A., Kluźniak W., 2016, MNRAS, 458, 666

Press W. H., Rybicki G. B., 1989, ApJ, 338, 277

Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49

Remillard R. A., Morgan E. H., McClintock J. E., Bailyn C. D., Orosz J. A., 1999, ApJ, 522, 397

Remillard R. A., Sobczak G. J., Muno M. P., McClintock J. E., 2002, ApJ, 564, 962

Rezzolla L., Yoshida S., Zanotti O., 2003, MNRAS, 344, 978

Rybicki G. B., Lightman A. P., 1986, Radiative Processes in Astrophysics. Wiley-VCH, Hoboken, NJ

Schnittman J. D., Bertschinger E., 2004, ApJ, 606, 1098

Schnittman J. D., Rezzolla L., 2006, ApJ, 637, L113

Schnittman J. D., Krolik J. H., Hawley J. F., 2006, ApJ, 651, 1031

Srámková E., Torkelsson U., Abramowicz M. A., 2007, A&A, 467, 641

Stella L., Vietri M., 1998, ApJ, 492, L59

Straub O., Šrámková E., 2009, Class. Quantum Gravity, 26, 055011

Strohmayer T. E., 2001a, ApJ, 552, L49

Strohmayer T. E., 2001b, ApJ, 554, L169

Stuchlík Z., Kotrlová A., Török G., 2013, A&A, 552, A10

Tagger M., Varnière P., 2006, ApJ, 652, 1457

Török G., Abramowicz M. A., Kluźniak W., Stuchlík Z., 2005, A&A, 436, 1

Lewin W., van der Klis M., eds, 2010, Compact Stellar X-ray Sources. Cambridge Univ. Press, Cambridge

Vincent F. H., Paumard T., Gourgoulhon E., Perrin G., 2011, Class. Quantum Gravity, 28, 225011

Vincent F. H., Meheut H., Varniere P., Paumard T., 2013, A&A, 551, A54

Vincent F. H., Mazur G. P., Straub O., Abramowicz M. A., Kluźniak W., Török G., Bakala P., 2014, A&A, 563, A109

Wagoner R. V., 1999, Phys. Rep., 311, 259

Zanotti O., Rezzolla L., Font J. A., 2003, MNRAS, 341, 832

Zanotti O., Font J. A., Rezzolla L., Montero P. J., 2005, MNRAS, 356, 1371

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Figure B1.** Pseudo-colour plot of mass density for the thinnest torus model with $R_c = 8.352$, $R_c - R_{in} = 0.2$ and $M_{BH} = 10 M_\odot$.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

<figure><img src="image_12.png" /><figcaption>**Figure A1.** PDS for various background densities. The torus is located at $R_c = 15 M_\odot$ and subjected to a radial perturbation. The value $b_0 = 10^{-7}$ corresponds to the background chosen in the simulations underlying Figs 7 and 8.</figcaption></figure>

## APPENDIX A: BACKGROUND EFFECT

Our numerical torus simulations require the presence of a low-density background fluid. However, it should be chosen in such a way that the modes of oscillations in the torus are neither suppressed nor enhanced due to the interaction between the torus surface and the atmosphere around it. In our main simulations listed in Table 3, we identify the plus modes in all radially perturbed tori. The breathing mode is identified in all our simulations presented in Figs 7 and 8.

In this appendix, we present a set of simulations with different background densities to explore whether or not the source of the plus and the breathing modes is an interaction of the torus surface with the background. In order to answer this question, we define the background density, $\rho_b = b_0\rho_{max}$, where $\rho_{max}$ is the central density of the torus, and vary $b_0$, which parameterizes the ratio of the background to central densities of the torus. Fig. A1 shows the PDS for five different background densities. Note that these PDS are computed from the $L_2$ norm of mass density, which is defined by $||\rho||^2 = \sum_{i=1}^{n_f} \sum_{j=1}^{n_\theta} \rho_{ij}^2$, where $n_f$ and $n_\theta$ are the number of shells along the radial and polar directions, instead of from light curves. It appears that the higher density backgrounds induce larger amplitudes in both the plus and breathing modes, although these modes are both present (in the radial perturbation case) at all background levels. The highest background density, which is $b_0 = 10^{-4}$, also slightly reduces the power in the radial mode. This may be due to an increase in the resistance of the background to the radial motion of the torus.

The fact that the PDS for $b_0 = 10^{-8}$, $10^{-7}$, $10^{-6}$ are almost identical gives a compelling argument that the plus and the breathing modes seen in our simulations are not triggered by the interaction of the torus surface and the background, but rather are present in the decomposition of the initial perturbation into normal modes. The breathing mode is present in all three types of perturbation considered here, while the plus mode is not excited by the vertical perturbation. The larger power in the plus and breathing mode for higher background density may be caused by the initial mismatch of pressure and gravity on the surface of the torus – the initial torus corresponds to an equilibrium solution in vacuum and its structure does not take into account the pressure exerted by the background
fluid present in the simulation.<fnref n="3" /> It may also be due to the motion of the torus through the background medium. The latter effect seems likely to be responsible for the enhancement of the plus mode, which is only present in radially perturbed tori. An object (torus in our case) moving in a high-density medium (background around torus) will be subject to resistance in its motion and hence, the torus will get squeezed along the direction of its motion. This squeezing during its motion may enhance the plus-mode (and possibly also breathing mode) eigenvector.

## APPENDIX B: MODE IDENTIFICATION

We found that the frequencies present in the PDS of the ray-traced light curves of perturbed 'nearly slender' tori (i.e. ones of a cross-section semimajor axis that is small compared to the radial distance of the torus from the black hole) match well to the corresponding analytic expressions (Blaes et al. 2006), which also have counterparts in slightly thicker tori (Straub & Šrámková 2009; Fragile, Straub & Blaes 2016). However, previous numerical work, beginning with Zanotti et al. (2003), studied tori with very large cross-sections, for which there are no corresponding analytic expression. Here, we would like to address the nature of the eigenmodes in such simulations.

In a series of studies, Zanotti et al. (2003), Montero, Rezzolla & Yoshida (2004) and Schnittman & Rezzolla (2006) suggested that the oscillation frequencies found in their simulations correspond to the eigenmodes of vertically integrated<fnref n="4" /> tori, as discussed in Rezzolla et al. (2003). The authors call this set of modes p modes. The two lowest-frequency eigenmodes were denoted by f (for fundamental) and o₁ (first overtone). Initially, it was thought that each of the detected frequencies was an eigenfrequency of a particular mode, but Zanotti et al. (2005) revisited the simulations and found that exciting 'overtones of the fundamental mode above o₁ is increasingly difficult,' and identified all frequencies in the simulation higher than f and o₁ as harmonics of f or o₁ or their combination. Thus, only the two lowest frequency eigenmodes of Zanotti et al. (2003) or Schnittman & Rezzolla (2006) are of interest here.

We note that the radial velocity eigenfunctions of the f and o₁ of the vertically integrated tori (fig. 2 of Rezzolla et al. 2003) correspond to the radial and plus mode of slender tori (Blaes et al. 2006), the former mode corresponding to a uniform radial velocity field, the latter to a velocity field antisymmetric with respect to a reflection in the vertical mid-plane of the torus. Further, as expressly noted by Zanotti et al. (2005) at the end of their section 5.1, 'the fundamental modes are closely related to the epicyclic oscillations' and 'the eigenfrequency of the fundamental p mode tends in the limit of vanishing size to the radial epicyclic frequency at the disc rest-mass density maximum.' We posit that (at least for the nearly slender tori considered in this work) the fundamental mode is in fact the radial epicyclic mode (Kluźniak 2005; Blaes et al. 2006). Regarding the o₁ mode, Blaes et al. (2006) stated that the plus mode in a compressible torus may indeed be characterized as a p mode. We therefore expect

that the radial (R) and plus (+) modes discussed by us correspond to the f and o₁ modes. See Appendix C below for a discussion of the dependence of the frequencies on the size of the torus and on the equation of state, which further confirms this point of view.

So far the identification of modes in our simulations was based on the coincidence of the observed frequencies with the eigenfrequencies of slender tori. In order to make the mode identification unambiguous, we tried two additional thinner torus ($R_c - R_{in} = 0.2$) setups using $\gamma = 4/3$ with $R_c = 8.352$ and $R_c = 15.0$ to distinguish the plus mode. This setup is perturbed using a Bondi kick to its radial velocity (as in Zanotti et al. 2003) with an amplitude of $\eta_1 = 10^{-4}$, which corresponds to strength of perturbation (Zanotti et al. 2003). This is a smaller amplitude perturbation than in any other of the simulations reported in this work. An examination of the motion of the radial position of the point of maximum density reveals that, correspondingly, the radial motion of the torus is also much smaller than in the other simulations. We were led to choosing such a small amplitude perturbation to ensure that this nearly slender torus is neither blown apart by the velocity kick, nor removed from the simulation domain. In an approach similar to that of the other simulations reported in this paper, we tested different computational domains to make sure that the inferred results from the simulations of these thinnest tori are not numerical artefacts. The tori with pressure maximum at $R_c = 8.352$ and $R_c = 15.0$ show similar oscillating patterns, and the corresponding eigenfrequencies scale as expected with $R_c$. In this description, we shall focus on the torus with $R_c = 8.352$.

For a torus this close to slender limit, a movie of the simulation clearly shows motions corresponding to the plus mode and the breathing mode. It is also possible to identify the frequencies directly, by noting the time of repetition of the patterns. In Fig. B1, we show snapshots from pseudo-colour plot of mass density for the torus with $R_c = 8.352$. The snapshots show a torus cross-section shape evolving in the expected way for a velocity field corresponding to the plus and breathing modes (Blaes et al. 2006). The left-hand column in Fig. B1 shows torus evolution for $t = 100.00, 194.12, 288.15 GM/c^3$ from top to bottom. Similarly, the second column shows snapshots for $t = 45.5, 94.02, 145.59 GM/c^3$. The apparent periodicity of the radial elongation of the torus is $\Delta t \approx 189 GM/c^3$, which translates to a frequency of $\sim 106$ Hz, comparing favourably with the plus-mode frequency of 106.11 Hz in the slender torus limit for a $10 M_{\odot}$ black hole and $\gamma = 4/3$. Vertical motions corresponding to expansions and contractions of the torus are also apparent with thickness maxima, e.g. at $t \approx 45$ and $t \approx 145 GM/c^3$, suggesting a period of $\Delta t \approx 98 GM/c^3$. This period corresponds to a frequency of 203 Hz, which is in good agreement with the corresponding slender torus limit for the breathing mode of $\sim 203.64$ Hz.

A PDS extracted from the $L_2$ norm of mass density (Fig. B2) for this simulation yields exactly same frequency values as those obtained from a visual inspection of the pattern at the beginning of the simulation, as described in the previous paragraph. In the PDS, a peak at $\sim 106$ Hz and at $\sim 203$ Hz corresponds to plus and breathing mode eigenfrequencies, respectively. These frequency values are in perfect agreement with the analytic estimates, equation (13) and equation (14). Owing to the very small amplitude of radial motion of the density maximum, or perhaps to the incompressible nature of the fundamental mode in the slender torus limit, the radial mode is not apparent in Fig. B2.

<footnote n="3">³ We note that a low-amplitude breathing mode was even observed in simulations of an *unperturbed* torus (Srámková, Torkelsson & Abramowicz 2007).</footnote>

<footnote n="4">⁴ We note that the fluid motions in the breathing mode, as well as in the vertical mode are primarily in the vertical direction and hence would not have been revealed by an analysis of a vertically integrated torus. The X mode also has a significant vertical velocity component.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure B1.** Pseudo-colour plot of mass density for the thinnest torus model with $R_c = 8.352$, $R_c - R_{in} = 0.2$ and $M_{BH} = 10 M_{\odot}$. The polytropic index $\gamma = 4/3$. The time in each snapshot is shown in $GM/c^3$ units. The colour bar shows mass density in black hole mass units for a $10 M_{\odot}$ black hole. The left- and right-hand panels show one complete cycle for plus and breathing modes, respectively. An online movie has been provided in the supplementary materials.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure B2.** PDS computed using $L_2$ norm of mass density for torus with pressure maximum at $R_c = 8.352$ and $R_c - R_{in} = 0.2$. The polytropic index is $\gamma = 4/3$. The peaks at 106 and 203 Hz correspond to the plus and breathing modes, respectively.</figcaption></figure>

## APPENDIX C: COMPARISON WITH PREVIOUS SIMULATIONS

The simulations listed in Table 3, correspond to oscillations of a nearly slender torus ($R_c - R_{in} = 0.5$). However, if we extend these simulation setups to thick torus models, we would expect to reproduce the results reported for previous thick torus simulations (e.g. Zanotti et al. 2003; Schnittman & Rezzolla 2006). Specifically, we attempt to reproduce model c in Zanotti et al. (2003). In this setup, the torus is centred at $R_c = 8.352$. Its inner edge is at $R_{in} = R_{cusp} = 4.576$ and the outer edge at $R_{out} = 15.889$. The black hole mass in our setup is $M_{BH} = 10 M_{\odot}$, consistent with Schnittman & Rezzolla (2006), whereas Zanotti et al. (2003) chose $2.5 M_{\odot}$. However, since the frequencies scale with the black hole mass, one can easily compare all of these results. Unlike our other simulations reported in Figs 7 and 8, this setup is perturbed using a Bondi radial velocity kick with amplitude $\eta_1 = 0.06$ (Zanotti et al. 2003). Also, these simulations have not been ray traced to compute light curves. Instead, the PDS are computed using the $L_2$ norm of mass density ($||\rho||^2$).
<figure><img src="image_15.png" /><figcaption>**Figure C1.** PDS computed using the $L_2$ norm of mass density for linear high (dashed red), logarithmic low (green dot-dashed) and high (solid black) resolution grids. The torus pressure maximum is initially at $R_c = 8.352$, and the inner edge at $R_{in} = 4.57$. The HR and LR correspond to resolution 256 × 256 and 250 × 84 (the same as model c of Zanotti et al. 2003), respectively. The peaks at 51 and 81 Hz correspond to the radial and plus mode, respectively.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure C2.** PDS computed using the $L_2$ norm of mass density for torus with $R_c = 8.352$, $R_{in} = 4.57$, for $\gamma = 4/3$ (solid line) and $\gamma = 5/3$ (dashed line). The peaks at 51 and 106 Hz correspond to radial and breathing mode. The 71 and 81 Hz peaks correspond to plus-mode eigenfrequency.</figcaption></figure>

In Fig. C1, PDS for three different simulations with low- and high-resolution logarithmic grids (green dot-dashed and solid black) and a linearly spaced radial grid (red dashed) are shown. There are three dominant peaks at 51, 81 and 106 Hz, consistent with the lowest three peaks inferred from model c of Zanotti et al. (2003), and Schnittman & Rezzolla (2006), who obtained 57, 85 and 114 Hz, respectively. Schnittman & Rezzolla (2006) and Zanotti et al. (2003)

<figure><img src="image_17.png" /><figcaption>**Figure C3.** Radial (black, dots) and plus (red, plus) mode frequency plotted against the cross-section size of the torus. The torus pressure maximum is initially at $R_c = 8.352$, the polytropic index is $\gamma = 4/3$. The largest torus corresponds to model c of Zanotti et al. (2003).</figcaption></figure>

identified the two lowest modes as both being pressure-driven oscillations. In Fig. C2, we compare PDS from simulations using two different polytropic indices. The overlap of the lowest-frequency peak for $\gamma = 4/3$ and $5/3$ suggest that this is not, in fact, a p mode. Instead, this result strongly suggests that the mode is driven by gravity, consistent with our identification as the radial epicyclic mode (Kluźniak 2005; Blaes et al. 2006). The second peak, which we associate with the plus mode, does show a dependence on the polytropic index $\gamma$, consistent with the plus mode being identified as a p mode.

We have further tested the radial- and plus-mode frequency dependence on the cross-sectional radius of the initial torus. We employed the aforementioned setup (centred at $R_c = 8.352$) while the cross-sectional size is varied from 1.0 to 11.3. As usual, we chose a background density of $b_0 = 10^{-7}$. Fig. C3 shows the identified frequencies of radial (black, dots) and plus (red, + signs) modes as a function of the size of the torus ($R_{out} - R_{in}$) using $\gamma = 4/3$. For small cross-sections, the radial- and plus-mode frequencies tend towards the nearly slender torus limit, i.e. ~71 and ~106 Hz, respectively. Note that our nearly slender torus setups reported in Appendix B also show a peak at 106 Hz. As the torus width increases, both the radial- and plus-mode eigenfrequencies decrease in a similar way and maintain a ratio close to 3:2. This trend is comparable to fig. 4 of Rezzolla et al. (2003) for their f and $\omega_1$ modes. We conclude that these are, in fact, the radial and plus modes identified in our thin-torus simulations.

<footnote>This paper has been typeset from a TeX/LA TeX file prepared by the author.</footnote>