© ESO, 2013

1. Introduction

Black hole binary systems exhibit a wide range of variability patterns. Among these, high-frequency quasi-periodic oscillations (HFQPO) appear as small narrow peaks in the power density spectrum. Their frequencies, which range from several tens to a few hundred Hertz, correspond to the rotation frequency in the inner part of the disk surrounding the black hole (see more details in the review of van der Klis 2006), making them an interesting tool to study strong-field gravity. Even though they are the fastest phenomena observed in the vicinity of black holes, HFQPO are also extremely elusive, especially when compared to lower frequency quasi-periodic oscillations (LFQPO). For this reason, we do not have a rich domain of observations from which to draw constraints for the models. This has led to a wide variety of models (see e.g. Lai & Tsang 2009, for a very brief review) that focus on one or another of the observed characteristics of the HFQPO. As more and more data are obtained, especially with next-generation telescopes such as LOFT (Bozzo et al. 2011), we will be able to gain a better picture of what a model must explain. Nevertheless, we already can make a stringent list of constraints based on the eight objects with recurrent HFQPO (see for example the review by Remillard & McClintock 2006). Among the most important of those constraints, we note:

• 1-

  The emitted flux is modulated at a level of a few percent in X-ray,and this modulation increases with energy. This first andfundamental point, often not discussed by models, is the aim ofthis paper.

• 2-

  The frequency of the modulation has a small but significant range, which must be explained.

• 3-

  The HFQPO occur sometimes, but not always, in pairs close to a 2:3 ratio. Those occurrences, and not just the close ratio, need to be elucidated.

• 4-

  Another point often neglected in models for HFQPO is the fact that they sometimes co-exist with the more ubiquitous LFQPO.

The model we are investigating here is based on the existence of the Rossby wave instability (RWI) at the inner edge of the disk as was discussed in Tagger & Varniere (2006). Indeed, it was shown that the existence of an innermost stable circular orbit (ISCO) around a black hole makes the disk prone to the RWI because the vorticity profile has a natural extremum. This dynamic instability results in the formation of large-scale spiral density waves and Rossby waves that can reach high amplitudes. Depending on the disk’s physical parameters, different modes of this instability will be selected, most often with azimuthal mode number m = 2 or 3 (Tagger & Varniere 2006). This gives a natural explanation for some of the observed characteristics of HFQPO. Nevertheless, the important step of computing the actual flux modulation when such an instability occurs in a disk was still missing. Here, we focus on this point and therefore restrict ourselves tothe case of the HFQPO alone. Indeed, we are trying to explore the capacity of one particular HFQPO model to modulate the X-ray flux, so having only one instability in the disk makes the results more conspicuous. To this end, we use the purely hydrodynamic version of the RWI (Lovelace et al. 1999), without the extra destabilizing effect provided by a poloidal magnetic field (Tagger & Varniere 2006) or the extra stabilizing effect of a toroidal magnetic field (Yu & Li 2009).

This is the continuation of the long string of studies on the impact of the RWI in disks, such as those done in the context of protoplanetary disks (Regály et al. 2012) and the galactic centre (Falanga et al. 2007). Whereas the first case is highly different from ours, the second is more comparable and we use similar tools. However, Falanga et al. (2007) restricted themselves to the 2D case and to the special context of the galactic centre source Sgr A*.

Our main interest here is thus to answer the key question of whether the RWI is able to modulate the flux of a microquasar at a level coherent with observations. In a more practical manner, we also address the necessity of full 3D simulations of the RWI and examine to what extent it is sufficient to restrict oneself to 2D simulations. Those 2D simulations are more accessible in terms of computing and storage resources, therefore allowing us to do much larger parameter studies and with a higher time resolution that could be compared with future LOFT observations. These questions will be addressed by computing the light curve of a microquasar subject to the RWI, as observed by a distant observer. This computation will be performed by a general relativistic ray-tracing code (Vincent et al. 2011). Because the main general relativistic effect that affects the RWI is the existence of an ISCO, we use a pseudo-Newtonian potential to model this in hydrodynamical simulations, instead of considering the full general relativistic case. This simplifying assumption is sufficient to derive a proof of principle of the ability of the RWI to modulate the flux of microquasars.

Section 2 briefly presents the RWI and then develops the hydrodynamical simulations that allow us to model the unstable accretion structure. The following section is focused on the ray-tracing procedure used to compute observable quantities. Section 4 presents the light curves that were obtained, and Sect. 5 gives conclusions and prospects.

2. Hydrodynamical simulation of the RWI

The RWI has been discussed in multiple astrophysical contexts, from the galactic centre to protoplanetary disks. It can be seen as the form the Kelvin-Helmholtz instability takes in differentially rotating disks and has a similar instability criterion. For 2D (vertically integrated) barotropic disks, the RWI can be triggered if an extremum exists in a quantity ℒ, which is the inverse vortensity¹$\mathrm{ℒ}\mathrm{=}\frac{\mathrm{\Sigma \Omega }}{\mathrm{2}{\mathit{\kappa }}^{\mathrm{2}}}\frac{\mathit{p}}{{\mathrm{\Sigma }}^{\mathit{\gamma }}}\mathit{,}$(1)where p is the pressure, Σ is the surface density, Ω is the rotation frequency, ${\mathit{\kappa }}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathrm{\Omega }\mathit{/}\mathit{r}\hspace{0.17em}\mathrm{d}\mathit{/}\mathrm{d}\mathit{r}{}^{\mathrm{(}}\mathit{r}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{)}}$ is the squared epicyclic frequency, and γ is the adiabatic index. An extremum of vortensity can thus typically arise from an extremum of the epicyclic frequency, as is the case close to the ISCO, or from an extremum in the density profile, which gives rise to the RWI in the disk. Once it is triggered, it leads to large-scale spiral density waves and Rossby vortices. The dominant mode (i.e. the number of vortices) is dependent on the disk conditions (Tagger & Varniere 2006).

2.1. The RWI as a model for HFQPO

Due to general relativistic effects, the epicyclic frequency profile will show a maximum in the vicinity of the black hole’s ISCO, creating a vortensity extremum. Tagger & Varniere (2006) showed that it would lead to the growth of the RWI at the inner edge of black hole binaries when the disk nears its ISCO. From these facts, it seems natural that the RWI model may be capable of accounting for points 2 to 4 in the Introduction above (which will not be investigated in detail here).

• 2-

  The frequency range: the exact location of the resultingvortensity extremum will depend on the density profile of the diskon top of the epicyclic frequency. Thus the location of thelaunching of the instability can vary depending on the diskproperties. As the Rossby vortices orbit around the central blackhole with the disk’s rotation frequency at the radius of launching,this variation of disk properties gives rise to a range of differentpossible frequencies of the signal.

• 3-

  Pairs of frequencies: different modes can dominate the signal depending on the physical state of the disk. This was adressed in Tagger & Varniere (2006) and more recently in Varniere et al. (2012b). The 2:3 frequency pairs can be naturally explained by a superposition of azimuthal modes m = 2 and m = 3.

• 4-

  Co-existence of HFQPO and LFQPO: we tackled this point by showing the ability of the RWI to co-exist with another instability, we proposed to be at the origin of the LFQPO (Varniere et al. 2012a), the accretion ejection instability (AEI). A disk giving rise to both AEI and RWI will thus naturally show both types of quasi-periodic oscillations.

It seems therefore natural to examine in more detail whether the RWI is strong enough to give rise to a modulation that could explain the observations.

2.2. Numerical setup

To perform the hydrodynamical simulation of the RWI, we use the Message Passing Interface-Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC) developed by Keppens et al. (2011). The numerical scheme is the same for all refinement levels, namely the total variation diminishing Lax-Friedrich scheme (see Tóth & Odstrčil 1996) with a third-order accurate Koren limiter (Koren 1993) on the primitive variables.

We consider a geometrically thin disk and neglect self-gravity. The Euler equations in cylindrical coordinates (r,ϕ,z) read $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{v}\mathit{\rho }\mathrm{\right)}& \mathrm{=}& \mathrm{0}\mathit{,}\\ {\mathit{\partial }}_{\mathit{t}}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{v}\mathit{\rho }{v}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathit{p}& \mathrm{=}& \mathrm{-}\mathit{\rho }\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}\mathit{,}\end{array}$where ρ is the mass density of the fluid, v its velocity, and p its pressure. We consider a barotropic flow, i.e. the entropy is constant in the entire system. The pressure is then p = Sρ^γ, with the adiabatic index γ = 5/3 and the constant S related to entropy. The sound speed is given by ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma p}\mathit{/}\mathit{\rho }\mathrm{=}\mathit{S\gamma }{\mathit{\rho }}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}$ and the temperature by T ∝ p/ρ ∝ Sρ^γ−1. All the distances are normalized by the ISCO radius r₀. The same applies for temperatures normalized at ISCO by T₀ = 10⁷ K, while all times are in units of the ISCO period t_ISCO.

The general relativistic effects have to be taken into account to study the innermost region of the black hole disk. To provide a proof of principle of the ability of the 3D RWI to modulate the flux of microquasars, we implement the minimum model necessary to mimic some of these general relativistic effects by using a pseudo-Newtonian potential Φ_G. For our simulations, we use the common Paczyńsky & Wiita (1980) gravitational potential ${\mathrm{\Phi }}_{\mathrm{G}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{r}}_{\mathrm{g}}}\mathit{,}$(4)where ${\mathit{r}}_{\mathrm{g}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{c}}^{\mathrm{2}}}$ is the gravitational radius, G the gravitational constant, M_∗ the black hole mass, c the speed of light in vacuum. With this gravitational potential, one can define the inner limit below which there is no stable circular orbit, the ISCO, located at r₀ = 6   r_g. We stress that this choice of potential allows us to mimic only some aspects of the Schwarzschild metric. The effect of the black hole spin is not taken into account. It is also important to keep in mind that this pseudo-Newtonian potential, although widely used, is far from giving an exact description of general relativity. However, the only characteristic of the Schwarzschild metric that is of primordial importance as far as RWI is concerned is the existence of an ISCO, and this is correctly modelled by Eq. (4).

We have also developed 2D RWI simulations with alternative pseudo-Newtonian potentials that take into account the rotation of the black hole. These alternative potentials were taken from Artemova et al. (1996). The RWI also develops well in these potentials, and future works will be dedicated to studying the effect of the black hole’s rotation on the observables.

2.3. Disk setup and boundary conditions

The grid is cylindrical with the unit length defined as the ISCO radius r₀. The radial coordinate r is in the range [0.8,6], the full azimuthal coordinate spans [0,2π], and the vertical coordinate z lies in the range [0,0.2] for the 3D simulations. For the fluid simulations, we considered only the upper part of the disk because the mid-plane is a symmetry plane for the RWI (Meheut et al. 2010, 2012). For the ray-tracing computations, the full disk is considered by symmetrizing the upper part. We used a fixed and homogeneous grid with a resolution n_r × n_ϕ = (384,72) in 2D and n_r × n_ϕ × n_z = (384,72,72) in 3D, which is slightly higher than the resolution used in Meheut et al. (2010).

As we describe gravitation by the pseudo-Newtonian potential in Eq. (4), the epicyclic frequency will show a maximum at a radius slightly higher than the ISCO radius r₀. This epicyclic frequency maximum gives rise to an extremum of the initial vortensity profile, which is shown in Fig. 1. The existence of this vortensity extremum does not depend on the choice of initial density profile, even though it does influence its precise location in the disk. Considering the density, we chose a typical power-law profile with a slope of −3/2. This choice, while reasonable, allows us to obtain a steeper vortensity extremum and hence a stronger instability. We note that the precise value of the rms (root mean square) of the modulation depends on this choice of slope. However, we are here only interested in proving the existence of the flux modulation and not in comparing a precise value of modulation to observations. If future instruments give better data on the rise of the HFQPO, we might be able to use these to constrain the disk model that could reproduce the observed data.

The initial density profile then reads

${\mathit{\rho }}_{\mathrm{M}}\mathrm{\equiv }\mathit{\rho }\mathrm{(}\mathit{r,\varphi ,z}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{2}}\left[\mathrm{1}\mathrm{+}\tanh \left(\frac{\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{B}}}{\mathit{\sigma }}\right)\right]{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathit{\alpha }}\mathit{,}$(5)the hydrostatic equilibrium gives the vertical profile of the density $\mathit{\rho }\mathrm{\left(}\mathit{r,\varphi ,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{M}}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{\mathit{\gamma S}{\mathit{\rho }}_{\mathrm{M}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\left[\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\mathit{r}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{r}}_{\mathrm{g}}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{r}}_{\mathrm{g}}}\right]\right)}^{\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\mathit{,}$(6)where r_B = 1.3   r₀ gives the position of the density maximum (its value was chosen to fit the simulations of Tagger & Varniere 2006), r_g =    r₀/6 is the gravitational radius, σ = 0.05   r₀ gives the width of the plunging region, α = −3/2 is the density slope, and ρ₀ is the minimum density of the simulation. These parameters, particularly the choice of density slope, can modify the structure and characteristics of the RWI. In our case, the existence of an extremum in κ limits the impact of the density power-law index, and any reasonable choice will exhibit the instability.

We also considered a 2D disk with the surface density defined as Σ ∝ r^−3/2 in the outer region. The initial conditions of the 2D and 3D simulations are then highly different. In both cases, the azimuthal velocity is determined by force balance: ${\mathit{v}}_{\mathit{\varphi }}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{r}}{\mathrm{(}\mathit{r}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{r}}_{\mathrm{g}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\mathit{\alpha S\gamma }{\mathit{\rho }}_{\mathrm{M}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\mathit{S\gamma }{\mathit{\rho }}_{\mathrm{M}}^{\mathit{\gamma }\mathrm{-}\mathrm{2}}\frac{\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\alpha }\mathrm{+}\mathrm{1}}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{2}\mathit{\sigma }{\cosh }^{\mathrm{2}}\frac{\mathit{r}\mathrm{-}{\mathit{r}}_{\mathit{B}}}{\mathit{\sigma }}}}\mathrm{·}$(7)The density and velocity profiles given in Eqs. ((6), (7)) include high-amplitude gradients and do not give an exact numerical equilibrium. For this reason, the initial conditions of our simulations are numerically computed from this pseudo-equilibrium. This means that a first simulation is done without any perturbations and is run until the disk has reached a permanent stage which is chosen similar to the initial conditions. This disk is then perturbed with small-amplitude (~10^-3   v_ϕ) random perturbations on the radial velocity, which act as seeds for the instability.

Fig. 1 Inverse vortensity profile, normalized to its extremum value. The instability is triggered at the location of the extremum, i.e. at r ≈ 1.4   r0.

The inner boundary condition is a no-inflow condition and the outer boundary condition is a null radial velocity condition. Meheut et al. (2012) have shown that the boundary conditions do not significantly change the growth rate of the RWI due to its confinement in the vortensity bump. Moreover, because the inner edge of the simulation is inside the ISCO, any reflected wave would not reach the region of interest for the instability. Therefore, the boundary conditions are not determinant for these simulations. For the 3D simulations, symmetric boundary conditions are implemented at the mid-plane, and we use a null vertical velocity at the grid upper boundary limit, which is situated outside of the disk.

2.4. Time evolution of the RWI

Figure 2 shows the density map of the z = 0 plane of a 3D RWI simulation when the instability is completely developed. This section presents the time evolution of the instability (at 2D and 3D), from its launch to its complete development.

Fig. 2 Mid-plane density map of the 3D simulation when the m = 1 mode dominates. The axes are in units of the ISCO radius.

At first, the disk is assumed to be 2D and the vertical structure of the disk is neglected. The growth of the instability is shown in Fig. 3, where the time evolution of the perturbation of density is plotted on a logarithmic scale. This allows the identification of both the linear phase of the instability, when the amplitude of the perturbations grows exponentially (thus linearly in logarithm), and the saturation, which is due to non-linearities. In the 2D case, the initial conditions give rise to a rapid growth of the instability on a time scale of a few periods at ISCO, as can be seen in Fig. 3. During the linear phase, the perturbations can be separated in the different azimuthal modes: $\mathit{X}\mathrm{=}\sum _{\mathit{m}\mathrm{\in }N}{\mathit{X}}_{\mathit{m}}\hspace{0.17em}\sin \mathrm{(}\mathit{m\omega t}\mathrm{-}\mathit{m\phi }\mathrm{)}\exp {\mathit{\gamma }}_{\mathit{m}}\mathit{t}\mathit{,}$(8)where X = Σ for 2D simulations and X = ρ for 3D. The quantities γ_m are the growth rates of each mode, ω the characteristic frequency of the instability depending on the position of the extremum of ℒ, m the azimuthal mode number, and X_m the amplitude of mode m. Therefore, X_{m    =    0} corresponds to the axisymmetric part of the density and the frequency of the mode m is the multiple mω of the frequency of the fundamental mode. Figure 3 also shows the time evolution of the amplitude of the strongest modes. During the linear phase, the dominant mode is m = 3, and later on, the disk is dominated by the fundamental mode m = 1, with important contributions from the modes m = 2 and 3. This evolution of the oscillation modes depends on the disk’s astrophysical properties: different modes dominate for different initial disk conditions.

Fig. 3 Amplitude of the surface density perturbation in logarithmic scale (solid line) for the 2D simulation. The amplitudes of the strongest modes are also plotted in dotted and dashed lines. The vertical dashed lines show the changes of dominant modes, indicated in italics at the top of the figure. The fundamental mode m = 1 eventually dominates.

This fact is illustrated in the 3D simulations. Indeed, the initial conditions differ largely between the 2D and 3D simulations, with a different surface density radial profile and absolute value. The characteristic velocity of the waves is then different. This modifies the time-scale for the growth of the instability as can be seen in Fig. 4. This difference is not due to the dimensionality (Meheut et al. 2012; Lin 2012). Nevertheless, the RWI is triggered, and saturation of the 3D instability is reached after a few tens of periods at ISCO. The dominant modes are also different during the linear phase, but non-linearities still tend to favour the lowest azimuthal mode.

Fig. 4 Same as Fig. 3 for 3D simulations. However, the time scale is different.

We note that, because the epicyclic frequency profile is not evolving with time, the location of the extremum of vortensity will stay the same. Moreover, the vortices grow in a density extremum where they cannot migrate due to the two density slopes of opposite sign. As a consequence, the Rossby vortices will not migrate during the simulation (for more details, see Meheut et al. 2012).

From the hydrodynamical simulations, we obtain the density and velocity at every point of our grid with a sufficient time sampling (of 40 frames in 2D and 6 frames in 3D per orbit at the ISCO) to follow the RWI in the disk from its rise to its saturation. In the next step, those will be used as input for the ray-tracing code in order to trace the impact of the RWI on the observed emission.

3. Raytracing the RWI with GYOTO

3.1. Raytracing the accretion disk

To obtain images, and especially light curves, that could be compared with observations, we used the general relativistic ray-tracing code GYOTO (Vincent et al. 2011). Null geodesics are computed backward in time in the Schwarzschild metric from a distant observer to the emitting disk. The coordinates used by GYOTO are spherical-like and denoted , where is used to differentiate the GYOTO spherical radius from the fluid simulations cylindrical radius r. Once a backward-integrated geodesic hits the disk, the density and 3-velocity at the point of emission is linearly interpolated at the time of emission from the results of the MPI-AMRVAC computations.

Fig. 5 Emission dates of photons at the footpoint of null geodesics connected to the observer’s screen, plotted with the same color scale for an inclination of 5° (left) and 85° (right). The grid of emission positions is projected onto the plane of the accretion disk. The color bar is expressed in units of the ISCO orbital time tISCO. The lensing effect of the black hole is clear for an inclination of 85° on the right panel, leading to a higher concentration of emission points on the side of the black hole opposite to the observer.

3.1.1. Radiative processes at 2D: blackbody

The 2D disk is assumed to be optically thick, so the integration is stopped at the first encounter of the disk. The emission is supposed to follow the Planck law: the only parameter needed is thus the temperature, which can be easily derived from the density according to the following computations.

The gas being assumed ideal: $\mathit{p}\mathrm{=}\frac{\mathit{\rho }}{\mathit{\mu }{\mathit{m}}_{\mathrm{u}}}\hspace{0.17em}\mathrm{k}\hspace{0.17em}\mathit{T,}$(9)where μ is the mean molecular weight, m_u is the atomic mass constant, and k is the Boltzmann constant. Assuming that the plasma is made of pure hydrogen, μ = 1. Using the relation R = N_A   k between the ideal gas constant R, the Avogadro number N_A, and the Boltzmann constant together with the expression of the sound speed, this yields: $\mathit{T}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{A}}\hspace{0.17em}{\mathit{m}}_{\mathrm{u}}}{\mathit{\gamma }}\hspace{0.17em}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathit{R}}\mathrm{·}$(10)This allows the computation of the emitted specific intensity at the surface of a 2D disk: ${\mathit{I}}_{\mathit{\nu }}\mathrm{=}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{,}$(11)where B_ν is the Planck function.

3.1.2. Radiative processes at 3D: Bremsstrahlung

For the 3D computations, the only radiative process considered is Bremsstrahlung. As the disk is purely hydrodynamic, there is no synchrotron emission. Moreover, Compton scattering is neglected because we are only interested here in a proof of principle, not in a detailed study of emission processes. The Bremsstrahlung emission is assumed to be thermal, so that the emission coefficient j_ν and absorption coefficient α_ν are related via Kirchhoff’s law: ${\mathit{j}}_{\mathit{\nu }}\mathrm{=}{\mathit{\alpha }}_{\mathit{\nu }}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathit{.}$(12)The emission coefficient for thermal Bremsstrahlung is given by (Rybicki & Lightman 1979) ${\mathit{j}}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\frac{{\mathrm{2}}^{\mathrm{5}}\mathit{\pi }{\mathit{e}}^{\mathrm{6}}}{\mathrm{3}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{3}}}{\left(\frac{\mathrm{2}\mathit{\pi }}{\mathrm{3}\mathit{k}{\mathit{m}}_{\mathrm{e}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathit{\rho }}{{\mathit{m}}_{\mathrm{u}}}\right)}^{\mathrm{2}}\exp \left(\mathrm{-}\frac{\mathit{h\nu }}{\mathit{kT}}\right)\mathit{,}$(13)where e is the electron charge, m_e is the electron mass, and h is the Planck constant. Here, we assume that the disk is made of pure hydrogen and that the emission is isotropic in the emitter’s frame (hence the 1/4π initial factor). Moreover, the Gaunt factor is neglected as most of the radiation arises from locations in the disk where hν ≈ kT.

Once the emission coefficient is computed, the absorption coefficient is derived by using Eq. (12).

3.1.3. Dependency on temperature at 2D and 3D

We investigate the dependency on temperature of the 3D emission process (Bremsstrahlung) as compared to the 2D case (blackbody). For temperatures around 10⁷ K, from where most of the flux arises, the 3D emission coefficient follows ${\mathit{j}}_{\mathit{\nu }}^{\mathrm{Br}}\mathrm{\propto }{\mathit{T}}^{\mathrm{2.5}}\exp {}^{\mathrm{(}}\mathrm{-}\mathit{h\nu }\mathit{/}\mathit{k}{\mathit{T}}^{\mathrm{)}}$, whereas the 2D specific intensity follows ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{BB}}\mathrm{\propto }\exp {}^{\mathrm{(}}\mathrm{-}\mathit{h\nu }\mathit{/}\mathit{k}{\mathit{T}}^{\mathrm{)}}$. The dependency on temperature is thus much more important in 3D, and the emission arises only from the hottest parts of the disk. In the 2D case, however, the emission is more spread out.

All these computations allow us to determine the specific intensity in the emitter’s frame, I_νem. To compute the specific intensity in the observer’s frame I_νobs, the frame invariant quantity I_ν/ν³ is used: ${\mathit{I}}_{{\mathit{\nu }}_{\mathrm{obs}}}\mathrm{=}{\mathit{\nu }}_{\mathrm{obs}}^{\mathrm{3}}\mathit{/}{\mathit{\nu }}_{\mathrm{em}}^{\mathrm{3}}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }{}_{\mathrm{em}}}$. The quantity ν_em can be related to the emitter’s 4-velocity. The 3-velocity computed by the fluid simulations is used to determine the emitter’s 4-velocity, and the observed specific intensity is thus at hand.

3.2. Computing the disk image and the light curve

Before producing an image of the disk, one needs to compute the relativistic time delay for each point of the disk in order to take into account the multiple time steps of the fluid simulations that will correspond to the same observed time. Figure 5 shows the emission date of photons for two extreme inclinations (5° and 85°). It also shows that the emission points are spread differently along the disk, depending on the inclination. The distribution is homogeneous and isotropic for low inclinations, but is very inhomogeneous and anisotropic for high inclinations. This anisotropy is due to the projection effect of the mapping of the observer’s screen onto a very inclined disk, resulting in the emission points being aligned preferentially along the direction perpendicular to the line of sight. The inhomogeneity is due to the lensing effect of the black hole that concentrates emission points on the side of the black hole opposite to the observer. Figure 5 also clearly shows that the higher the inclination, the more time slices of data will be required. The respective difference between the maximal and minimal emission dates on the primary image are approximately 0.1   t_ISCO (left-hand panel) and 0.9   t_ISCO (right-hand panel). However, due to the strong beaming effect at high inclination, the bulk of the total specific intensity in one image comes from a small part of the disk. Thus, the emission dates of the photons reaching this small part of the disk are close to each other, and the final effect of time delay is only marginal on the light curve. We have checked that the difference between the exact light curve and a light curve computed without taking into account the time delay is at the level of one to a few percent only.

Fig. 6 Image of a 2D optically thick disk emitting blackbody radiation with an inclination of 85° and in the presence of the RWI. The first-order image (the distorted complete ring) clearly shows the spiral shape of the emitting region. The beaming effect makes the emission brighter on the approaching side of the disk (here, on the left of the image). The second-order image is the portion of ring at the bottom of the image. It is due to photon swirling around the black hole before reaching the observer. The third-order image is the thin ring of illuminated pixels at the center of the image, which is due to photons orbiting around the black hole very close to the event horizon before reaching the observer. This ring is the so-called black hole silhouette; it is truncated here due to optical thickness.

Figure 6 shows the image (i.e. the map of specific intensity) of the accretion disk at an inclination² of 85°, approximately at the time when the fundamental mode of the RWI dominates. Each pixel of the image is obtained by interpolating between the set of simulated data at the time of emission of the photon, as stressed in Sect. 3. Here, around 40 different data time slices are used for computing the image. The instability has started to grow, and one can identify the spiral density waves of the RWI. On top of this modulation due to the instability, there is a clear beaming effect: matter moving toward the observer is brighter (here, on the left side of the image). The secondary image of the disk is visible as a semicircle at the center of the primary image.

Figure 7 depicts the image of a 3D disk subject to the RWI, approximately at the time when the fundamental mode of the RWI dominates. Here, around six different data time slices are used for computing the image: the time sampling has been reduced as 3D data are much heavier than 2D data. We checked at 2D that this reduced sampling still allows us to retrieve similar results. As compared to Fig. 6, the 3D disk’s emission is much more concentrated on the inner parts of the disk. This is due to the much stronger dependency of the 3D emission on temperature as compared to the 2D case, as explained in Sect. 3.1.3. As the hottest parts of the disk are concentrated close to the radius of launching of the instability (see Fig. 2), only these inner regions are seen.

Fig. 7 Image of a 3D disk emitting Bremsstrahlung radiation with an inclination of 85° and in the presence of the RWI. The overall emission is much more concentrated in the vicinity of the Rossby vortex as compared to Fig. 6.

4. Modulation of the light curve by the RWI

This section presents and analyzes the light curves obtained by raytracing the RWI hydrodynamical simulations. These are computed by summing the specific intensities over all solid angles on the observer’s sky. This boils down to summing the disk images over all pixels, for all observation times.

4.1. 2D case: high time resolution

Fig. 8 Light curves of a microquasar subject to the 2D RWI at 1 keV (left) and 2 keV (right) at different inclinations. The flux is normalized to its initial value. The inclination is 5° (dash-dotted black), 45° (dashed blue), or 85° (solid red).

We first analyze the main characteristics of the light curve modulation in the 2D case, where we can have a much better time resolution between the hydrodynamical snapshots³. To catch the details of the light curve evolution, 40 frames of fluid simulations are computed during each period of rotation at ISCO. The ray-tracing code then interpolates between these fluid simulation results in order to compute the specific intensity map at different times of observation. We compared the results obtained with the one from a higher (80 frames) and lower (20 frames and 6 frames) time resolution. While the overall behaviour was similar, 40 frames allow us to get a more detailed light curve, similar to that with 80 frames. We therefore decided to do all the runs in 2D at 40 frames per period of rotation at ISCO. We computed the light curve from the moment the RWI started to grow in the disk until it reaches its non-linear states, as shown by the amplitude evolution in Fig. 3.

At zero inclination, the light curve will slowly drift toward smaller values due to the black hole’s accretion during the simulation. To determine the flux fluctuation that is only due to the RWI, it is important to subtract this continuous drift from the light curve. This is done by simply subtracting the light curve obtained at 1° of inclination from every other light curve. As the GYOTO code uses Boyer-Lindquist spherical-like coordinates, the z-axis is singular, and it is not possible to derive a light curve at exactly zero inclination. However, the light curve is only marginally impacted by beaming at 1° as compared to exactly 0°. The error introduced by subtracting the light curve at 1° of inclination is a small underestimation of the amplitude of the modulation that does not impact our result here. Indeed, we are only looking for a proof of principle that the RWI could modulate the flux at a level compatible with the observed HFQPO.

Figure 8 shows the light curves obtained for different values of the inclination parameter at two different values of the energy of the observed photon: 1 and 2 keV. Initially when the instability has a very low amplitude, the modulation of the flux is very low. After a few ISCO rotations, the instability is growing exponentially and modulates the flux at a detectable level. The period T_mod of the modulation equals the period of rotation of the Rossby vortices, i.e. the period of rotation at r ≈ 1.4   r₀ (see Fig. 1), that is to say T_mod ≈ 1.6   t_ISCO.

Due to the general relativistic beaming effect, when the Rossby vortex is on the approaching side of the disk, the resulting flux is boosted at high inclination. The opposite is true for the receding side of the disk. This beaming effect is also visible in the light curve substructures that can be observed at high inclination. For instance, an m = 2 mode will lead to two sharp peaks in the light curve at high inclination, related to the passage of the Rossby vortices at the approaching side of the disk. If the inclination is low, this effect is much fainter, resulting in smaller substructures in the light curve. The oscillation frequency of the light curve evolves during the simulation with a high frequency after a few rotations. After five rotations, a mixture of modes between one frequency and twice this frequency, can be identified. At the end of the simulation, however, the mode with the lowest frequency dominates. This evolution corresponds to the evolution of modes seen in Fig. 3: the dominant mode initially has a frequency 3ω, then 2ω, and eventually the fundamental mode dominates.

Fig. 9 Dispersion of the light curve points in Fig. 8 at 1 keV (left) and 2 keV (right). Each point is obtained by computing the dispersion of the light curve points over two orbital periods. The inclination is 5° (black triangles), 45° (blue diamonds), or 85° (red crosses).

The comparison between the two values of energy of the observed photon shows that similar shapes are obtained for the two cases, but a higher modulation amplitude is reached for higher energy as expected due to the Planck law. Indeed, because the ISCO temperature is 10⁷ K, the maximum frequency of the Planck law is closer to 2 keV that to 1 keV.

Figure 9 shows the light curve rms, computed with a sliding window with a width of two orbital periods. The amplitude of the modulation is growing linearly with a change of slope at around 5   t_ISCO, corresponding to the saturation of the instability (see Fig. 3). The maximum level of modulation is around 4% at 1 keV and 8% at 2 keV. This is somewhat higher than most observations of HFQPO. However, the maximum’s exact value is influenced by the initial ℒ profile (see Eq. (1) and Fig. 1). Here, we are only interested in seeing if a reasonable setup can modulate the flux to a level similar to the one observed. We are keeping the detailed study of the growth rate and saturation level for a forthcoming paper (Meheut et al. 2013).

4.2. Full 3D case

We then turned to 3D simulation to see if there was any difference in the observables. There, we could not get a similar time resolution as in 2D because of the hydrodynamical data size. We therefore used the 40 frames per orbit resolution for only one period for confirmation, while we made a longer term light curve at the much lower time resolution of six frames per orbit at the ISCO. As we see below, this lower resolution is still able to catch the broad aspects of the light curve, although it will not allow us to obtain the finest details.

Figure 10 depicts the light curve and its rms for a 3D disk subject to the RWI, for photons of observed energy equal to 2 keV, and for three different values of inclination. The first important result is that the RWI is capable of modulating the light curve, similar to what had already been shown in the 2D case, with a modulation of a few percent. This is a strong argument that makes the RWI a reliable model of microquasar HFQPO.

The domination of the different modes is also easily seen in Fig. 10, with the mode m = 2 dominating at the beginning of the simulation, and eventually the mode m = 1. Here, the mode m = 1 dominates after around 15 ISCO periods, whereas it dominates after around 30 periods in Fig. 4. This difference is due to the fact that the ray-tracing simulations are initiated when the instability is already strong enough to give a non-negligible rms. As a consequence, the 15 first ISCO periods of the 3D hydrodynamics computations were not used as their rms is extremely low.

The right-hand panel of Fig. 10 shows the same characteristics as in the equivalent 2D Fig. 9. The rms shows a linear profile with a clear change of slope around t ≈ 7   t_ISCO. This is linked to mode m = 1 starting to dominate over mode m = 2.

This being given, there are two main differences between the 2D and 3D results:

• 1-

  The 3D modulation grows more slowly than its 2D counterpart.This comes from the differences in initial disk conditions,particularly the choice of the density profile: different densityprofiles give different shapes for the extremum of ℒ, which in turngives different growths of the instability. Figures 3and 4 already showed that the3D growth rate is slower than its 2D counterpartwhich translates directly to the light curve. Moreover, as statedabove, the 3D light curve initial point is 15 ISCO periods after thelaunch of the instability. We stress that the index of the density profile is not constrained by observations. As the difference of growth time depends on the choice of the power-law index, it cannot be taken as an intrinsic difference between 2D and 3D RWI.

• 2-

  The light curve exhibits a slight asymmetry in the 3D case (the oscillation is not exactly centered on the initial flux value). This could be explained by the strong dependency on the temperature of the 3D emission, as explained in Sect. 3.1.3. During the simulation, the temperature of the Rossby vortex increases by a few percent, due to accretion. Its emission thus becomes greater and greater. This translates to a slight drift of the light curve maxima toward higher values.

Fig. 10 Left: light curves of a microquasar subject to the 3D RWI at 2 keV at an inclination of 5° (dash-dotted black), 45° (dashed blue), or 85° (solid red). Right: dispersion of the light curve points of the left panel. Each point is obtained by computing the dispersion of the light curve points over two orbital periods. The inclination is 5° (black triangles), 45° (blue diamonds), or 85° (red crosses). For both panels, the time-scale is different compared to Figs. 8 and 9.

4.3. Are 3D simulations necessary to account for observations?

Comparing Figs. 8 to 10, it first appears that the light curve is modulated both in 2D and 3D, which is the main result of this paper. The dependency of the light curve as a function of the inclination parameter and frequency of the radiation is also very similar. As these general features of the 2D and 3D computations are alike, an interesting conclusion is that 2D results are sufficient to analyze the general observable characteristics of the RWI, at least until we get a better constraint on the profiles in the disk.

This is particularly interesting when considering the difference in terms of computing time and memory resources between a 2D and 3D simulation. The hydrodynamical data at a given time is around 100 times heavier at 3D than 2D (typically respectively 50 Mo and 0.5 Mo).The time needed to compute an image from 3D data is generally 15 times what is needed for a 2D computation.

However, we stress that the present work does not allow us to compare the detailed time evolution of the 3D RWI with the 2D case, due to our choice of time sampling in the ray-tracing computations (this choice is dictated by the computing time and memory resources needed for the 3D simulations). As the 3D RWI displays specific features in the z direction (Meheut et al. 2010), it may be that specific observable characteristics could be obtained only by resorting to high time-resolution 3D simulations. Nevertheless, these features would be small corrections to the general trend that stays close to the 2D results and would not be within reach of current instruments. Indeed, the 2D time sampling of 40 frames per ISCO period implies a time resolution of around 10^-4 s for the observed light curve. This is far beyond current instrumental capabilities.

Moreover, we stress that the radiative processes and radiative transfer are treated in a very simplified way in this study for the 3D case. The 2D simulations are thus sufficient only when one is not interested in studying precisely the radiative properties of the disk.

5. Conclusion

The RWI was previously proposed as a model for HFQPO, and we have now demonstrated its ability to modulate the flux coming from the disk. Using 2D and 3D hydrodynamical simulations, we have also been able to study how the amplitude of this modulation evolves as a function of the source inclination and of the radiation frequency. The 2D simulations have been shown to be sufficient to recover the broad characteristics of the light curve.

By using hydrodynamical simulations, we were able to focus only on the RWI and its effects on the light curve. If this can be considered for the case where HFQPO occur alone in a softer state, as is the case for the 67 Hz modulation of GRS 1915+105, for example (Morgan et al. 1997), HFQPO are more commonly observed in the steep power-law or hard-intermediate state simultaneously with a LFQPO. In Varniere et al. (2012a), we have demonstrated in 2D the ability of the RWI to co-exist with another instability that could give rise to the LFQPO. Our future work will be devoted to that particular case, above all, how it

influences the observables. Indeed, this state is much more frequent during microquasar outbursts than the softer state we studied here.

Acknowledgments

This work has been financially supported by the French GdR PCHE and Campus Spatial Paris Diderot. Some of the simulations were performed using HPC resources from GENCI-CINES (Grant 2012046810).

References

All Figures

	Fig. 1 Inverse vortensity profile, normalized to its extremum value. The instability is triggered at the location of the extremum, i.e. at r ≈ 1.4   r0.
In the text

	Fig. 2 Mid-plane density map of the 3D simulation when the m = 1 mode dominates. The axes are in units of the ISCO radius.
In the text

	Fig. 3 Amplitude of the surface density perturbation in logarithmic scale (solid line) for the 2D simulation. The amplitudes of the strongest modes are also plotted in dotted and dashed lines. The vertical dashed lines show the changes of dominant modes, indicated in italics at the top of the figure. The fundamental mode m = 1 eventually dominates.
In the text

	Fig. 4 Same as Fig. 3 for 3D simulations. However, the time scale is different.
In the text

Fig. 5 Emission dates of photons at the footpoint of null geodesics connected to the observer’s screen, plotted with the same color scale for an inclination of 5° (left) and 85° (right). The grid of emission positions is projected onto the plane of the accretion disk. The color bar is expressed in units of the ISCO orbital time tISCO. The lensing effect of the black hole is clear for an inclination of 85° on the right panel, leading to a higher concentration of emission points on the side of the black hole opposite to the observer.

In the text

Fig. 6 Image of a 2D optically thick disk emitting blackbody radiation with an inclination of 85° and in the presence of the RWI. The first-order image (the distorted complete ring) clearly shows the spiral shape of the emitting region. The beaming effect makes the emission brighter on the approaching side of the disk (here, on the left of the image). The second-order image is the portion of ring at the bottom of the image. It is due to photon swirling around the black hole before reaching the observer. The third-order image is the thin ring of illuminated pixels at the center of the image, which is due to photons orbiting around the black hole very close to the event horizon before reaching the observer. This ring is the so-called black hole silhouette; it is truncated here due to optical thickness.

In the text

	Fig. 7 Image of a 3D disk emitting Bremsstrahlung radiation with an inclination of 85° and in the presence of the RWI. The overall emission is much more concentrated in the vicinity of the Rossby vortex as compared to Fig. 6.
In the text

	Fig. 8 Light curves of a microquasar subject to the 2D RWI at 1 keV (left) and 2 keV (right) at different inclinations. The flux is normalized to its initial value. The inclination is 5° (dash-dotted black), 45° (dashed blue), or 85° (solid red).
In the text

	Fig. 9 Dispersion of the light curve points in Fig. 8 at 1 keV (left) and 2 keV (right). Each point is obtained by computing the dispersion of the light curve points over two orbital periods. The inclination is 5° (black triangles), 45° (blue diamonds), or 85° (red crosses).
In the text

Fig. 10 Left: light curves of a microquasar subject to the 3D RWI at 2 keV at an inclination of 5° (dash-dotted black), 45° (dashed blue), or 85° (solid red). Right: dispersion of the light curve points of the left panel. Each point is obtained by computing the dispersion of the light curve points over two orbital periods. The inclination is 5° (black triangles), 45° (blue diamonds), or 85° (red crosses). For both panels, the time-scale is different compared to Figs. 8 and 9.

In the text