# Impact of the gravity of a Schwarzschild black hole upon the Rossby wave instability

F. Casse, P. Varniere, Z. Meliani

► To cite this version:

F. Casse, P. Varniere, Z. Meliani. Impact of the gravity of a Schwarzschild black hole upon the Rossby wave instability. Monthly Notices of the Royal Astronomical Society, 2017, 464 (3), pp.3704-3712. 10.1093/mnras/stw2572. hal-01554351

HAL Id: hal-01554351

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

Submitted on 24 Aug 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
# Impact of the gravity of a Schwarzschild black hole upon the Rossby wave instability

F. Casse,<fnref n="1" />* P. Varniere<fnref n="1" /> and Z. Meliani<fnref n="2" />

Accepted 2016 October 6. Received 2016 October 3; in original form 2016 July 12; Editorial Decision 2016 October 3

## ABSTRACT

In an early work, the Rossby wave instability was proposed to explain the variability thought to originate in the close vicinity of black holes but this was done in the pseudo-Newtonian approach. Here we present the first general relativistic (GR) hydrodynamics simulations of this instability not only proving its theorized existence in a full GR environment but also studying the effect of the strong gravity on the instability. To that end, we performed a set of simulations increasingly closer to the black hole with our new GR version of the MPI-AMRVAC code. This allows us to study the minute changes in the behaviour of the instability. We found that a pseudo-Newtonian approach gives adequate results provided that the time-shifting induced by the black hole gravity is taken into account. Hence, to view the disc as a distant observer would, a full GR ray-tracing post-treatment of the simulations is a must.

**Key words:** accretion, accretion discs – hydrodynamics – instabilities – relativistic processes – stars: black holes.

## 1 INTRODUCTION

The Rossby wave instability (RWI) was first introduced for a galactic disc by Lovelace & Hohlfeld (1978) and dubbed, at the time, the negative-mass instability. The difficulty of finding a physical setup fulfilling its criterion made it rarely studied outside of carefully crafted configurations (Lovelace et al. 1999; Li et al. 2000, 2001). Over the past decade, the RWI has been the subject of increasing interest from a wide range of astrophysical communities ranging from planet formation (Varnière & Tagger 2006; Meheut et al. 2010; Lyra & Mac Low 2012; Lin 2012) to fast variability of microquasars (Tagger & Varnière 2006; Varniere, Tagger & Rodriguez 2011; Varnière, Tagger & Rodriguez 2012; Vincent et al. 2013) or the flares of Sgr A* (Tagger & Melia 2006; Vincent et al. 2014). That last case is because circularized matter falling towards a black hole will meet the criteria for the RWI. Such phenomena have the propensity to create a wide variety of light curves displaying time-scales coherent with the observed radiative emission coming from the galactic centre.

Here we will look in more detail at how the location of the Rossby vortices with respect to the last stable orbit of the black hole does influence its physics and what changes, if any, we should be looking for in observations. To that end, we will run a set of five general

relativistic hydrodynamic (GRHD) simulations, each fulfilling the RWI criterion at a different location in the disc.

We want to stress here that the simulations of our set are independent of one another. Indeed, considering in a single simulation an annulus of gas falling on to the black hole would prevent us from disentangling general relativistic (GR) effects from dynamical effects stemming from the motion of the gas.

In the next section, we will briefly present the basic features of the RWI and the GRHD numerical code MPI-AMRVAC used to perform the GRHD simulations presented in this paper. In Section 3, we will describe the initial setup of the various simulations as well as the diagnostics to track the changes in the behaviour of the instability. We also present in this section the overall methodology. Section 4 is devoted to the presentation of our 2D simulations and the complete display of the RWI features depending on the location of the instability with respect to the central black hole. In Section 5, we describe the 3D GRHD simulations of the RWI and compare them to previous 2D computations. Finally, Section 6 will present our conclusions.

## 2 RWI IN GRHD

As we are interested in tracking the changes in the evolution of the RWI as the instability develops closer and closer to a non-rotating black hole, we will briefly summarize its basic features in the Newtonian case as well as the code we will be using.

<footnote n="1">¹APC – Sorbonne Paris Cité, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris 10, rue Alice Domon et Léonie Duquet, F-75205 Paris Cedex 13, France</footnote>

<footnote n="2">²LUTH – Laboratoire de l'Univers et de ses Théories, CNRS UMR 8102, Observatoire de Paris, Université Paris Diderot, 5 place Jules Janssen 92195 Meudon Cedex, France</footnote>

<footnote>* E-mail: fcasse@apc.univ-paris7.fr</footnote>
Our presentation of the RWI is by no mean exhaustive and for a proper review see, for example, Lovelace & Romanova (2014).

## 2.1 Basic features of the RWI

The RWI is a global instability found in a differentially rotating disc exhibiting an extremum of the vorticity $\mathcal{L}$, whose definition for an unmagnetized 2D disc is

$$ \mathcal{L} = (\nabla \times \mathbf{v})_{\perp} / \Sigma, \quad (1) $$

where $\mathbf{v}$ is the velocity of the disc matter while $\Sigma$ is its height-integrated density ($\perp$ stands for the component perpendicular to the disc plane). The position of the extremum is called the corotation radius $r_c$ of the wave, meaning that at this radius the wave and the gas have the same velocity.

For a long time, the existence of such an extremum was a bottleneck that made this instability underused. Since then, a few scenarios creating bumps in the vorticity have been found, leading to a more in-depth study of the RWI, especially for a protoplanetary disc.

Once there is an extremum of $\mathcal{L}$ in the disc, a standing wave pattern can form, which creates vortices. As the energy and angular momentum are transferred from inside the corotation (positive energy) to outside the corotation, we get an exponential growth of the vortices corresponding to the linear growth stage of the instability. Then, either the bumps are eaten away, leading to a decrease of the instability, which then dies with the extremum, or it is replenished (by whatever caused it in the first place) and we get a saturation phase.

## 2.2 General relativistic formalism in MPI-AMRVAC

As we intend to study the evolution of the RWI in the close vicinity of a Schwarzschild black hole, it is mandatory for us to use a numerical code encompassing a full GRHD framework. From our experience with the existing MPI-AMRVAC code (van der Holst et al. 2012), we have implemented full GR support in it. Here we briefly present our approach and a test of the code.

### 2.2.1 GRHD equations of polytropic fluids

The governing equations of GRHD stem from fluid mass, momentum and energy conservation expressed in a GR framework. A perfect non-magnetized fluid is characterized by its momentum 4-vector $\rho u^\alpha$ and stress-energy tensor $T^{\mu\nu}$, namely,

$$ T^{\mu\nu} = \rho h u^\mu u^\nu + P g^{\mu\nu}, $$

where $\rho$ is the proper density of the fluid, $u^\mu$ is the contravariant 4-velocity, $P$ the thermal pressure and $h$ the specific enthalpy of the fluid. The gravitational influence of the central object appears through the presence of a non-flat metric $g^{\mu\nu}$, whose expression depends on the nature of the central object. The general expression of a line element is, using a (3 + 1) splitting of space-time,

$$ ds^2 = -\alpha^2 c^2 dt^2 + \gamma_{ij} (dx^i + \beta^i dt)(dx^j + \beta^j dt) $$

where $\gamma_{ij}$ is the spatial metric tensor while $\alpha$ is the lapse function and $\beta^i$ is the shift vector. Let us note that Latin indices stand for spatial coordinates while Greek indices are linked to all space-time coordinates. In this first study, we will consider non-rotating black holes described in Boyer-Lindquist coordinates $(r, \theta, \varphi)$ by the Schwarzschild metric, namely,

$$ \alpha^2 = \left(1 - \frac{2r_g}{r}\right) = \frac{1}{\gamma_{rr}}, \quad \gamma_{0\theta} = r^2, \quad \gamma_{\varphi\varphi} = r^2 \sin\theta, $$

where $r_g = GM/c^2$ is the gravitational radius of the black hole. The temporal coordinate $t$ used in our description is consistent with time experienced by a so-called Newtonian observer located at infinity where the metric tensor converges towards a flat spherical metric expression. Mass and momentum conservation is expressed as the vanishing covariant derivative of the fluid momentum 4-vector and stress-energy tensor:

$$ \rho u_{;\mu}^{\mu} = 0; \quad T_{;\mu}^{\mu\nu} = 0 $$

where $; \mu$ stands for the covariant derivative with respect to coordinate $\mu$. The corresponding set of GRHD equations is

$$ \begin{aligned} \partial_t(\sqrt{\gamma}D) + \partial_j(\sqrt{\gamma}D(\alpha v^j - \beta^j)) &= 0, \\ \partial_t(\sqrt{\gamma}S_i) + \partial_j(\sqrt{\gamma}[S_i(\alpha v^j - \beta^j) + \alpha P \delta_i^j]) &= \sqrt{\gamma}\{-(\xi - P)\partial_i\alpha + \frac{\alpha}{2}(S^j v^k + P\gamma^{jk})\partial_i\gamma_{jk} + S_j\partial_i\beta^j\}, \quad (2) \end{aligned} $$

where $D = \rho W$ is the relativistic mass density with $W = (1 - v_i^i)^{-1/2}$ being the Lorentz factor of the fluid. The relativistic inertia of the fluid is $\xi = \rho h W^2$ while $S_i = \xi v_i$ is the relativistic momentum where $v_i$ is the covariant coordinate velocity of the fluid ($u_i = Wv_i$) normalized to the speed of light $c$. For simplicity and robustness, we choose to replace the energy conservation by a simple polytropic law linking the thermal pressure to the proper gas density, namely $P = C_o \rho \bar{\gamma}$ where $C_o$ and $\bar{\gamma}$ are two positive constants. Relativistic gas thermodynamics prevents the use of a standard polytropic equation of state linking the internal energy of the gas to its density. Following Meliani et al. (2004) and Mignone & McKinney (2007), we can derive such a relation by considering the properties of the distribution function of a relativistic gas (Taub 1948; Mathews 1971). This leads to an expression for the internal energy density $U$ as

$$ U = \frac{P}{\Gamma - 1} + \sqrt{\frac{P^2}{(\Gamma - 1)^2} + \rho^2 c^4}, $$

where $\Gamma = C_p/C_V$ is the ratio of the specific heat capacities. Then, the enthalpy can be directly computed from the previous equation since

$$ \rho h c^2 = \frac{1}{2} \left\{ (\Gamma + 1)U - (\Gamma - 1) \frac{\rho^2 c^4}{U} \right\}. $$

We developed a GRHD version of the MPI-AMRVAC code that solves the aforementioned conservative equations and, thus, provides the temporal evolution of conservative variables ($D, S_i$). The structure of the code remains fully parallel by using an octree data setting and does conserve the same scaling properties as in Porth et al. (2014). Solving the set of GRHD conservation equations requires access simultaneously to both conservative ($D, S_i$) and primitive ($\rho, v_i, W, \xi$) variables. Switching from one set of variables to the other can be achieved by solving a non-linear equation relating the Lorentz factor to the gas enthalpy, namely,

$$ f(W) = D^2 h^2 c^4 (W^2 - 1) - S_i S^i = 0. $$

A Newton-Raphson algorithm can then get the zeros of such a function with quadratic efficiency. This procedure must be repeated for every cell of the computational domain at every time step, leading to a slight increase of the computational time compared to classical computations.
### 2.2.2 GRHD code validation

In this paper, we focus on non-spinning black holes as described by the Schwarzschild metric shown in equation (2). Such a metric is stationary and exhibits a spherical symmetry. The inclusion of GR components in the special relativistic MPI-AMRVAC code then requires us to test the ability of the code to handle the radial motion of matter influenced by the gravity of the black hole. One simple but efficient test is to describe the GR Bondi accretion regime (Hawley, Smarr & Wilson 1984). In this test, one considers the pressureless free-falling motion of matter along a radial geodesic. Neglecting pressure prevents any interaction between the different elements of matter composing the flow. This flow is then subject only to the gravitational force generated by the central body. Regarding our polytropic equation of state, we achieve this condition by imposing $C_o \ll 1$.

Radial geodesics in Schwarzschild space-time are described by a very simple equation, namely $du_\mu/ds = 0$. This is the expression of the conservation of the binding energy of the fluid element while free-falling on to the black hole. The normalization of the 4-velocity of matter $u_\mu u^\mu = -1$ provides the expression of the contravariant radial velocity $v^r$ while mass conservation provides the expression of the relativistic density $D$, namely (Hawley et al. 1984),

$$v^r(r) = \sqrt{\frac{2r_g}{r}} \left(1 - \frac{2r_g}{r}\right), \\ D(r) = \frac{D_o}{r^2 \sqrt{\frac{2r_g}{r} \left(1 - \frac{2r_g}{r}\right)}} \quad (3)$$

We set up a spherically symmetric GRHD simulation of this flow where we initially tune the density and the radial velocity $\sqrt{v_r v^r}$ to constant values while the orthoradial and azimuthal velocities are null. The computational domain ranges from $2.2r_g$ to $200r_g$ while we set $C_o = 10^{-10}$ to render the thermal pressure negligible compared to the mass energy. We ran the simulation up to time $t = 1410 r_g/c$, at which it reaches a nearly perfect steady state (the residual simulation drops below $10^{-11}$). In Fig. 1, we display the radial evolution of both density and radial velocity at different times, showing that they converge towards the Bondi solution represented by the diamond symbols. The relative overall variance between the final stage of the simulation and the Bondi solution is smaller than 0.1 percent, which proves that the simulation has converged towards the expected solution.

## 3 SETUP AND DIAGNOSTICS FOR THE RWI

### 3.1 Numerical setup

To trigger and follow the growth of the RWI inside a rotating disc, we have to design initial conditions exhibiting an extremum of the vortensity while achieving a perfect balance between the various forces acting on the gas. We choose the simplest setup possible for such a disc, namely, having a localized enhancement of the density or bump.<fnref n="1" /> We, therefore, choose the following radial density profile:

$$\Sigma(r) = \Sigma_o \sqrt{\frac{6r_g}{r}} \left(1 + 3 \exp\left\{-\left(\frac{r-r_c}{2\Delta}\right)^2\right\}\right), \quad (4)$$

<figure><img src="image_2.png" /><figcaption>Figure 1. Radial contravariant velocity and relativistic density of a spherically symmetric pressureless flow free-falling on to a non-rotating black hole. Both quantities reach a steady state matching the Bondi accretion regime described in Section 2.2.2 and represented by the diamond symbols.</figcaption></figure>

where $r_c$ is the corotation radius and $\Delta$ is the typical width of the density bump. In the context of GRHD, the radial equilibrium is enforced by tuning the rotational momentum $S^\varphi$ such that the radial projection of the centrifugal force balances both the radial gravitational force and the thermal gradient. On the other hand, there is no gravitational force acting on the gas in the $\theta$ direction, which requires a balance between the $\theta$-component of the centrifugal force and the pressure gradient component along the same direction. The 3D setup of the disc is then the solution of the following implicit system:

$$\begin{align}
\partial_{\theta} P &= \alpha S^{\varphi 2} \cos \theta, \notag \\
S^{\varphi}(r, \theta) &= \sqrt{\frac{\xi}{r} \left( \xi \frac{\partial_r \alpha}{\alpha} + \partial_r P \right)}, \tag{5}
\end{align}$$

<footnote n="1">¹ While mostly used for simplicity and the cleanness of the results, this could originate from some blob of matter falling on to the system and getting circularized before becoming accreted, as was proposed for Sgr A*.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Initial density distribution in the (r, θ) plane of an accretion disc harbouring a density bump at $r_c = 12r_g$.</figcaption></figure>

which we solved numerically. The resulting density distribution can be seen in Fig. 2 for $r_c = 12r_g$.

To complete the initial disc setup, we choose to use a polytropic relation linking the thermal pressure to the density. The polytropic index is set to 3 while $C_o = 1.8 \times 10^{-4}$. The values of these constants have been selected such that the ratio of the sound speed to the rotational velocity at the disc midplane remains constant throughout the disc ($c_s/\sqrt{v_\varphi v^\varphi} \sim 0.04$). This choice is consistent with a thin accretion disc structure except near the bump where the local disc scale height is a significant fraction of its spherical radius. Such a thin disc displaying a localized bloated part meets all the criteria to trigger a very rapid growth of the RWI, hence preventing any interference from the boundaries occurring within the computational domain.

We seed the instability near the bump by adding random velocity perturbations to the disc equilibrium:

$$ \delta v_r = \xi_{ijk} v_\varphi \exp \left( - \left( \frac{r-r_c}{\Delta} \right)^2 \right), $$

where $\xi_{ijk}$ is a random variable verifying $\langle\xi\rangle = 0$ and $\langle\xi^2\rangle = 10^{-10}$ when averaged over the whole computational domain.

The computational domain of all the simulations presented in this paper is $r \in [0.5r_c, 1.5r_c]$ and $\varphi \in [0, 2\pi]$. The spatial resolution for the 2D simulations is $192 \times 600$ cells while full 3D simulations will be performed using the same resolution and extension in the $(r, \varphi)$ plane and $\theta \in [0.47\pi, \pi/2]$ with 32 cells in that direction. The boundaries of the computational domain are periodic in the azimuthal direction while the radial ones are designed to absorb the velocity/density perturbations from the random perturbations spread near the density bump.

To time-advance the set of equations presented in the previous section, we used a Harten, Lax and van Leer solver linked to a Koren slope limiter (Koren 1993). A typical 2D simulation requires approximately 10 000 time steps on 20 processors and lasts less than an hour.

3D simulations take into account the full disc structure and, thus, have to deal with the presence of a low-density region above and below the disc. We assume the disc to be embedded in such a low-density region, which should not interfere with the RWI development. To design such a buffer zone, we impose that below a typical floor density (typically a hundredth of the disc midplane density), no motion is allowed. Typical 3D simulations exhibiting the aforementioned resolution ran for about 30 h on 20 processors to reach the same physical time as the 2D simulations.

## 3.2 Methodology

In this paper, we intend to study the influence of a Schwarzschild black hole on the evolution of the RWI depending on how close to the compact object this instability develops in the surrounding

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Rotational velocity $\sqrt{v_\varphi v^\varphi}$ of the disc normalized to its value at the corotation radius for the various 2D initial setups.</figcaption></figure>

disc. To this end, we will perform a set of simulations for different corotation radii, namely $r_c \in \{300, 60, 30, 18, 12\}r_g$. To compare the different simulations, we scaled distances in each simulation with respect to their corotation radius. Accordingly, time will be scaled to $T_c(r_c)$, the rotation period at the corotation radius as measured by a local observer located at a distance $r_c$ from the black hole. These various simulations represent distinct test cases where we study the growth of the RWI occurring at locations closer and closer to the black hole.

Fig. 3 shows the behaviour of the rotational velocity scaled to its value at $r_c$ for each simulation. The corresponding velocity profiles are not exactly identical as the gravity is not an exact power law. However, they do not significantly differ one from another. Such a similar setup will then allow us to compare the evolution of the RWI in the various simulations.

To compare our computations with previous works, we start our analysis with the simulation having the corotation radius located farthest away from the last stable orbit of the black hole, namely $r_c = 300r_g$. For several reasons, we have selected a disc configuration where we expect strong and fast growth of the instability. First, this ensures that any parasitic instability will not have time to develop before the RWI reaches its non-linear stage and, therefore, this warrants that we are studying only the RWI. Secondly, it prevents any interaction with the edges of the simulation box as the spiral waves do not reach them before the end of the simulation, as can be seen in Fig. 4.

This is not an issue when running a simulation with $r_c = 300r_g$ because there is no hard limit for the simulation box. Nevertheless, as we get closer to the last stable orbit of the black hole, we need to ensure that we are still confronted with the same boundary conditions. Hence, we set the radial extent of the smallest simulation to $[6r_g, 18r_g]$ so that it will match the relative radial size of the farthest simulation while ensuring that the simulation is not influenced by the unstable part of the flow located inside the last stable orbit of the disc.
<figure><img src="image_5.png" /><figcaption>Figure 4. Modified logarithmic colour map of the density perturbations as a function of the radius and time for $\varphi_0 = \pi$ at $r_c = 300r_g$. Thick black lines represent the limit of the box on which we did the Fourier transform. Time is normalized with respect to the rotation period $T(300r_g)$ as measured by an observer located at $r_c = 300r_g$.</figcaption></figure>

### 3.3 Set of RWI diagnostics: near-Newtonian case at $r_c = 300r_g$

The first simulation considered in this paper, at $r_c = 300r_g$, is consistent with a near-Newtonian case. Indeed, that far away from the last stable orbit of the black hole, the gravitational field generated by the compact object is very similar to a Newtonian one. We can then compare the different features of this first computation with previous simulations of the RWI performed in a pure Newtonian case. To be certain that we are actually dealing with the RWI, the first thing to check is the criterion in equation (1). The aforementioned disc setup has been designed solely for that purpose so the criterion is actually fulfilled, as can be seen in Fig. 5.

On top of that, we see that the behaviour of $\mathcal{L}$ is relatively similar for all the simulations in the set with only a difference of about 10 per cent near the extremum. This will allow us to compare the results and growth rate of the RWI between the different runs. Fig. 5 also shows that the vortensity profiles get flatter as the inner edge of the simulation is getting closer and closer to the last stable orbit around the black hole. Such flattening denotes the impact of a full GR description of the black hole's gravity upon the velocity.

Once we know the criterion is fulfilled, we can seed some random perturbations around the corotation and follow how they evolve with time. This is what we see in Fig. 4, which shows, stacked in time, the relative density perturbations:

$$ \frac{\rho(r, \varphi_0, t) - \rho_0}{\langle \rho(r, \varphi, t) \rangle_{\varphi}}. \qquad (6) $$

Such a plot makes it easy to track propagating waves as time evolves. Indeed, we can follow the evolution of the density perturbations and the spiral waves that develop from the corotation zone where seed perturbations were added. This ensures that the evolution is as expected for the RWI and also that no problem occurs from the interaction with the boundaries of the simulation box.

Another way to ensure the identity of the instability is to decompose it into its compressional ($u_{ij}^l$) and torsional ($u_{ij}^r - u_{ij}^\varphi$) components, which can be seen in Fig. 6. On those plots, we can track the evolution of the spiral waves as with the density, though much earlier. The torsional component is also the best way to look at the formation and evolution of the vortices. On those plots, the

<figure><img src="image_6.png" /><figcaption>Figure 5. RWI criterion $\mathcal{L}$ for the entire set of 2D simulations normalized to $\mathcal{L}^*(r_c)$ where $\mathcal{L}^*$ corresponds to the criterion of a disc without any density bump. The subplot represents a zoom of the minima where we can actually see the small differences in our setups.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Logarithmic colour map of the divergence ($u_{ij}^l$) and linear colour map of the curl ($u_{ij}^r - u_{ij}^\varphi$) of the velocity field showing the spiral waves and vortices after one orbit at the corotation radius $r_c = 300r_g$.</figcaption></figure>

size of the inner disc radius can be inferred from the size of the central dark spot, which represents the event horizon of the black hole.

We see that for this near-Newtonian case ($r_c = 300r_g$), the compressional and torsional components behave as expected with the spiral waves propagating from the corotation where vortices are also formed. From the number of vortices, we can find which instability modes are present. While this is easy when there is only a very limited number of modes or when one of them is strongly dominant, it generally requires using the Fourier space to dissociate the various instability modes fully. From Fig. 9, we can see that the strength distribution of instability modes is not uniform, hinting that one or more low $m$ modes are getting stronger than the high $m$ modes we get from the start [see Tagger & Varnière (2006) for a more detailed mode study].
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Logarithmic colour map of the divergence ($u_{\phi}^i$) and linear colour map of the curl ($u_{\phi\phi}^i - u_{\phi\varphi}^i$) of the velocity field showing the spiral waves and vortices after one orbit at the corotation radius $r_c = 12r_g$.</figcaption></figure>

# 4 IMPACT OF THE GRAVITY OF A SCHWARZSCHILD BLACK HOLE UPON THE RWI

Now that we have recovered the near-Newtonian case at $r_c = 300r_g$, we can perform a series of simulations where a similar bloated disc is considered with a corotation radius $r_c$ closer and closer to the black hole. We choose to run the set with $r_c \in \{300, 60, 30, 18, 12\}r_g$. For reasons of space, and since they are very similar, we present the diagnostics only for $r_c = 12r_g$.

## 4.1 Comparison in real space

When looking at the simulation in real space, namely the $(r, \varphi)$ domain, we obtain very similar diagnostics as we did for the near-Newtonian case. For example, we see in Fig. 7 the compressional and rotational components of the velocity for the simulation at $r_c = 12r_g$, which is to be compared with Fig. 6 where the corotation was 25 times further away from the black hole (the black circle at the centre is, hence, 25 times smaller). We see that the main behaviour of the RWI is unaffected by its closeness to the black hole, its overall behaviour being indistinguishable from the near-Newtonian case.

Nevertheless, we detect some slight differences when looking more closely at the density perturbations and especially their time evolution as seen in Fig. 8. Indeed, according to our simulation design, the inwardly propagating sound wave (generated by the initial perturbations) should reach the inner border after one local orbital time $T_c(r_c)$ while the outwardly propagating sound wave will reach the outer border after approximately 1.27 $T_c(r_c)$. Our simulation setup makes these values independent of $r_c$. Looking at Fig. 4, we recover this feature but not in Fig. 8 where both waves seem to reach the borders approximately at the same time. This discrepancy stems from a GR local time dilatation occurring as one gets closer and closer to the black hole. Indeed, general relativity shows that the clock rate $d\tau/dt$ of a referential located at distance $r$ from a Schwarzschild black hole is provided by the lapse function, namely $d\tau/dt = \alpha(r)$. In contrast to the near-Newtonian simulation where the lapse function $\alpha$ is very close to unity over the whole computational domain, the simulation having $r_c = 12r_g$ exhibits a radially increasing $\alpha$ that leads to a stretching (along the time axis) of the left part of the plot as time dilatation is more significant on the inner side of the corotation radius than on its outer side. One has to keep in mind at this stage that the velocity and density maps presented in the previous figures are *not* the image that a Newtonian observer would perceive, since light emitted from the disc is also subject to a relativistic distortion induced by the black hole gravity.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Modified logarithmic colour map of the density perturbations as a function of the radius and time for $\varphi_0 = \pi$ at $r_c = 12r_g$. Thick black lines represent the limit of the box on which we did the Fourier transform. Time is normalized to the rotation period $T(12r_g)$ as measured by an observer located at $r_c = 12r_g$.</figcaption></figure>

## 4.2 Comparison in the Fourier domain

Since all the plots presenting physical quantities in the real space domain are so similar, we have decided to study the evolution of the growth rate of the first few instability modes in the Fourier domain. Indeed, it is an interesting observable diagnostic as the presence of a perturbation mode in the disc would be translated as a flux modulation, which could be observable.<fnref n="2" /> Here we will look at the growth rate of those modes to see if being closer to the last stable orbits influences what we would observe.

For each simulation of our set, we have computed the Fourier transform of the rest density $\rho$ inside the box that encompasses the entire span of the spiral wave up to one orbit at the corotation. The boundary of the box is visible in Figs 4 and 8 and is similar for all the simulations in between.

In Fig. 9, we display the measured growth rate of the first 19 azimuthal modes of the instability for all the simulations of our set, meaning $r_c \in \{300, 60, 30, 18, 12\}r_g$. At the top, we have the early behaviour of the RWI mode, namely before $T_c/2$. This confirms what we can see directly on the plots of the divergence and curl of the velocity, Figs 6 and 7, namely that high $m$ modes are dominant in the disc with the maximum growth rate in this early rise of the instability being for the mode $m=12$ closely followed by the modes between 10 and 14.

As time evolves, we see on the bottom panel of Fig. 9 that, while the high $m$ modes stay at the same level as they were originally, the lower $m$ modes become dominant. Ultimately, this will lead to the merging of the vortices, as has been observed in all the simulations of the RWI that were performed for a longer period (see for example Varnière & Tagger 2006; Vincent et al. 2013). Here we cannot push the simulation that far, as the bump was put in *by hand* and will eventually be destroyed by the instability. The study in the Fourier space allows us to recover the start of this merging with the low $m$ mode having the strongest growth rate as early as one orbital period.

<footnote n="2">² Indeed, that several low modes can have a high enough growth rate to be detectable has been proposed as an explanation for the apparent multiplicity of frequency-related High Frequency Quasi-Periodic Oscillations (HFQPO) in microquasars (Tagger & Varnière 2006).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Growth rate $\gamma$ of the first 19 azimuthal modes for every location of the corotation radius $r_c$ studied in 2D. The top panel stands for growth rates measured during the initial linear phase of the instability while the bottom one displays growth rates measured during the second half of the simulation.</figcaption></figure>

The most interesting result of Fig. 9 is the overall similarity of the behaviour as the corotation gets closer to the last stable orbit of the black hole. While this shows that, even in a full GR simulation, the basic features of the RWI stay the same, we do see a decrease in the growth rate. While we try to keep the disc structure, Fig. 3, and instability criterion, Fig. 5, as similar as possible between all the simulations, there is a net decrease as the corotation radius gets closer to the black hole. It is difficult to untangle this change in the initial conditions from a potential GR effect.

## 5 RWI IN FULLY 3D ACCRETION DISCS NEAR A NON-ROTATING BLACK HOLE

It was shown previously that the 2D behaviour and 3D behaviour of the RWI were similar enough to warrant the use of a 2D simulation when comparing with observations until we get better observational data (Vincent et al. 2013). However, this work was done in the pseudo-Newtonian potential approach. Here we have the opportunity to test this result without an approximation.

The setup we took is very similar to the one we took for the 2D case, as we aim to compare the results of both the 2D and 3D simulations. This setup is presented in Section 3, and Fig. 2 depicts the mass density of the disc with $r_c = 12r_g$. As the 3D simulations are both more time- and space-consuming, we choose to focus only on the near-Newtonian case at $r_c = 300r_g$ and the case where the inner edge of the disc is at the last stable orbit $r_c = 12r_g$. Especially after showing how similar the 2D runs are, both extreme cases will give us a good view of the 3D behaviour.

Again, we start comparing the simulations in the real spatial domain, looking, for example in Fig. 10, at the behaviour of the divergence for the cases at both $r_c = 300r_g$ and $r_c = 12r_g$. Note that this time the colour map is linear to identify the contours of the various vortices (at the disc midplane and along the bump height) and, therefore, it cannot be compared directly to Figs 6 and 7. We see once again that the behaviour of the two extreme simulations of our set is very similar.

We also identify the behaviour already noted in Vincent et al. (2013) corresponding to the 3D simulations having dominant low $m$ modes earlier than in the 2D simulations. Indeed, we did recover this behaviour of merging vortices towards low $m$ modes in a few of the 2D simulations that we ran longer than one orbit but choose not to exploit the data here as the bump was starting to be modified by the instability and the results could not be easily compared.

To look at this behaviour in more detail, we performed the same time-domain analysis as in the 2D case, this time doing the Fourier transform for $\theta \in [0.49\pi, 0.5\pi]$, which represents most of the height of the disc without the bump. We choose to cut out part of

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Linear colour map of the divergence ($u_{;i}^i$) at $r_c = 300r_g$ and $r_c = 12r_g$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Growth rate $\gamma$ of the first 19 azimuthal modes for $r_c = 12r_g$ and $r_c = 300r_g$ computed from both the 2D and 3D simulations. The left-hand panel stands for growth rates measured during the initial linear phase of the instability while the right one displays growth rates measured during the second half of the simulation.</figcaption></figure>

the bump in our analysis, since having the entire bump would also add a large part of the low-density region. Not having the entire bump means that some modes are underestimated but this way we do not have any interference from the low-density region in the Fourier transform. Fig. 11 compares the growth rate $\gamma$ of the first 19 azimuthal modes for $r_c = 12r_g$ and $r_c = 300r_g$ computed from both the 2D and 3D simulations. We see that the overall behaviour between the 2D and 3D runs is in good agreement. In the early phase of the instability, the growth rates of the 3D simulation computed in the disc only roughly match the 2D ones. As previously stated, the slightly lower growth rate for a few modes is understandable as we have artificially cut part of the bump out of our Fourier transform, so a mode there does not have the full growth rate. For the growth rate computed between $0.5T_c$ and $T_c$, we see, as was already noted by Vincent et al. (2013), that the low $m$ modes in the 3D run are stronger than in the 2D cases. Indeed, those 3D modes are much closer to the non-linear evolution we see in 2D shortly after one orbit when the vortices have had a little more time to merge. Despite that the 3D simulations are closer to the full development of the instability, the resulting differences between the 2D and 3D computations remain small. The 2D simulations are then likely to be a good approximation of the RWI evolution, especially in the context of comparing with observations from available instruments, whose spatial resolution prevents us from distinguishing small-scale details in the disc structure.

## 6 CONCLUSIONS AND OBSERVATIONAL CONSEQUENCES

The RWI has previously been proposed to explain several phenomena occurring in the vicinity of black holes (Vincent et al. 2013, 2014) but its existence in a full-GR environment had never been demonstrated up to now. Here we remedy this by performing 2D and 3D GRHD simulations of the RWI that enable us to prove at the same time:

(i) that the instability actually exists when taking into account the full effects of a non-spinning black hole

(ii) that its local behaviour is only mildly affected by the strong gravity generated from a non-spinning black hole.

Nevertheless, while the GR effects are mild on the growth and saturation of the instability, this is not what an observer at infinity would perceive, first because the local clock rate depends on the location of the emitting region but also because of the way the emitted photons would travel to the observer. Indeed, a remote observer would perceive distorted instability patterns as a consequence of the time dilatation occurring near the black hole. For a Schwarzschild

black hole as presented here, this effect is relatively weak as the typical variation length scale of the lapse function $\alpha$ is large compared to the size of the region where the instability develops. This implies that to depict the RWI near a non-spinning black hole, we can keep a mostly pseudo-Newtonian approach for the fluid dynamics as long as we take into account the time dilation. The second effect needed to compute observables is to take into account the way the emitted light would arrive at the observer through GR ray-tracing (see, e.g., Falanga et al. 2007). For a Schwarzschild black hole, this will be limited to high-inclination systems and still it would be hard to constrain as it is entangled with other effects for most of the observables we have access to; see, for example, Varniere & Vincent (2016) for the impact on an instability of the same family. To see some potential effect in the observations, one would need enough precision in time to follow the changes in the shape of the pulse profile within a few orbits. Indeed, while the mode structure is not widely influenced by GR, the way it is perceived by an observer changes as seen in Varniere & Vincent (2016) and can lead to detectable differences in the pulse profile. Addressing such differences, in the context of a more realistic setup for observables like the flares of Sgr A*, is beyond the scope of this paper as it would also require us to take into account the potential spin of the black hole.

## ACKNOWLEDGEMENTS

We acknowledge the financial support from the UnivEarthS Labex programme of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02).

## REFERENCES

* Falanga M., Melia F., Tagger M., Goldwurm A., Bélanger G., 2007, ApJ, 662, L15

* Hawley J. F., Smarr L. L., Wilson J. R., 1984, ApJ, 277, 296

* Koren B., 1993, Numerical Methods for Advection Diffusion Problems. Vieweg, Braunschweig, p. 117

* Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023

* Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874

* Lin M.-K., 2012, MNRAS, 426, 3211

* Lovelace R. V. E., Hohlfeld R. G., 1978, ApJ, 221, 51

* Lovelace R. V. E., Romanova M. M., 2014, Fluid Dynamics Res., 46, 041401

* Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805

* Lyra W., Mac Low M.-M., 2012, ApJ, 756, 62

* Mathews W. G., 1971, ApJ, 165, 147

* Meheut H., Casse F., Varniere P., Tagger M., 2010, A&A, 516, A31

* Meliani Z., Sauty C., Tsinganos K., Vlahakis N., 2004, A&A, 425, 773

* Mignone A., McKinney J. C., 2007, MNRAS, 378, 1118
Porth O., Xia C., Hendrix T., Moschou S. P., Keppens R., 2014, ApJS,
214, 4

Tagger M., Melia F., 2006, ApJ, 636, L33

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

Taub A. H., 1948, Phys. Rev., 74, 328

van der Holst B., Keppens R., Meliani Z., Porth O., van Marle
A. J., Delmont P., Xia C., 2012, MPI-AMRVAC: MPI-Adaptive Mesh
Refinement-Versatile Advection Code, Astrophysics Source Code Li-
brary (ascl:1208.014)

Varnière P., Tagger M., 2006, A&A, 446, L13

Varniere P., Vincent F. H., 2016, A&A, 591, A36

Varniere P., Tagger M., Rodriguez J., 2011, A&A, 525, A87

Varnière P., Tagger M., Rodriguez J., 2012, A&A, 545, A40

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

Vincent F. H., Paumard T., Perrin G., Varniere P., Casse F., Eisenhauer F., Gillessen S., Armitage P. J., 2014, MNRAS, 441, 3477

This paper has been typeset from a *TEX/LATEX* file prepared by the author.