Oscillations of vertically integrated relativistic tori – I. Axisymmetric modes in a Schwarzschild space-time

Abstract

This is the first of a series of papers investigating the oscillation properties of relativistic, non-self-gravitating tori orbiting around a black hole. In this initial paper we consider the axisymmetric oscillation modes of a torus constructed in a Schwarzschild space-time. To simplify the treatment and make it as analytical as possible, we build our tori with vertically integrated and vertically averaged quantities, thus transforming the eigenvalue problem into a set of coupled ordinary differential equations. The tori are also modelled with a number of different non-Keplerian distributions of specific angular momentum, and we discuss how the oscillation properties change when different distributions of angular momentum are considered. Our investigation progresses by steps. We first consider a local analysis in Newtonian gravity and determine the properties of acoustic wave propagation within these objects, as well as the relations between acoustic and epicyclic oscillations. Next, we extend the local analysis to a general relativistic framework. Finally, we perform a global analysis and determine both the eigenfunctions and the eigenfrequencies of the axisymmetric oscillations corresponding to the p modes of relativistic tori. These behave as sound waves globally trapped in the torus and possess eigenfrequencies appearing in the simple sequence 2:3:4:…, independently of the distribution of angular momentum considered. The properties of the modes investigated here are in good agreement with those observed in recent numerical simulations, and could have a number of different applications. In X-ray binary systems containing a black hole candidate, for instance, p-mode oscillations could be used to explain the harmonic relations in the high-frequency quasi-periodic oscillations observed. In systems comprising a massive torus orbiting a black hole, on the other hand, p-mode oscillations could be used to explain the development or the suppression of the runaway instability.

accretion, accretion discs, black hole physics, hydrodynamics, relativity

1 Introduction

Waves and normal-mode oscillations in stars have long since attracted attention and a vast literature, investigating them both in Newtonian [see, for instance, Cox (1980) or Unno et al. (1989)] and in general relativistic regimes [see Stergioulas (1998) for a review], is now available. On the other hand, waves and normal-mode oscillations in geometrically thin discs around compact objects have been studied much less, both within Newtonian gravity [see Kato (2001) for a review] and within a relativistic framework (Okazaki, Kato & Fukue 1987; Perez et al. 1997; Silbergleit, Wagoner & Rodriguez 2001; Kato 2001; Rodriguez, Silbergleit & Wagoner 2002). The literature is even more scarce when one considers geometrically thick discs, which so far have been investigated mostly in connection with their stability properties both in Newtonian gravity (Papaloizou & Pringle 1984, 1985; Blaes 1985) and in general relativity (Kojima 1986). An explanation for why discoseismology has not yet reached the development and the level of sophistication that are now possible in asteroseismology is due, at least in part, to the fact that only recently have accretion discs been recognized as fundamental astrophysical objects, present in various forms at all scales. Nowadays, however, periodic and quasi-periodic variations are currently observed in different classes of astrophysical objects containing accretion discs. Although many different models have been proposed for the interpretation of the rich phenomenology associated with these quasi-periodic oscillations (QPOs), there seems not to be yet a widely accepted mechanism for most of the observed sources [see van der Klis et al. (2000) for a review]. Clearly, a systematic investigation of the oscillation properties of discs could shed some light on this.

As in stars, oscillation modes in discs are, in general, the consequence of restoring forces responding to perturbations, and these offer a way for classifying oscillations. In accretion discs, in particular, a first restoring force is given by the centrifugal force, which is responsible for the appearance of the so called inertial waves, tightly related to the orbital motion of the disc and hence to epicyclic oscillations. A second restoring force is offered by pressure gradients, and the oscillations that are produced in this way are related to p modes and have close connections with the propagation of sound waves in the perturbed fluid. A third restoring force is the gravitational field in the direction orthogonal to the orbital plane. If a portion of the disc is perturbed in the vertical direction, in fact, the vertical component of the gravitational field will produce a harmonic oscillation across the equatorial plane with oscillation frequency equal to the orbital frequency. These oscillations are related to corrugation waves[see Kato, Fukue & Mineshige (1998) for an overview on disc oscillations].

Two complementary approaches have been followed traditionally for studying the perturbations of equilibrium models of geometrically thin accretion discs. The first one is a local approach, and it has been used extensively to investigate the propagation of waves with the inclusion of many contributing physical effects, such as buoyancy, stratified atmospheres and magnetic fields (see, among others, Lubow & Ogilvie 1998). Being local, these approaches derive dispersion relations in which the frequency of the perturbation is expressed as a function of the spatial position within the object and as the linear combination of different contributions, each related to the different physical effect (Kato 2001). The second approach is a global one and it is based on the fact that, under suitable conditions, Eulerian perturbations to all physical quantities can be expressed as an eigenvalue problem, the eigenvalue being the frequency of oscillation. In practice, and for the simplest case, this amounts to solving a second-order partial differential equation once appropriate boundary conditions are provided (see Nowak & Wagoner 1991; Ipser & Lindblom 1992; Silbergleit et al. 2001).

This paper is devoted to both a local and a global perturbative analysis of axisymmetric modes of oscillation of relativistic tori in the Cowling approximation (i.e. in an approximation in which the perturbations of the space-time are neglected: Cowling 1941). Fluid tori differ from geometrically thin discs mostly in having equilibrium configurations in which the orbital motion is intrinsically non-Keplerian and in which pressure gradients play an important role, giving rise to an extended vertical structure. As a result, a consistent investigation of these configurations which would account for the coupling of the oscillations in the radial and vertical directions requires necessarily a two-dimensional treatment involving a set of partial differential equations. While this problem is solvable with presently available techniques (this has indeed already been done for relativistic rotating stars: Yoshida & Eriguchi 1997), it would require a massive use of expensive numerical calculations, leaving little room for a physical interpretation. Furthermore, in contrast with the case of relativistic stars, the eigenvalue problem for relativistic tori is for the main part unsolved, and resorting to a fully numerical solution at this stage would prevent one from appreciating most of the basic physics behind the oscillation modes of these objects. For this reason, we will here concentrate on a simpler model for the torus which can be handled for the main part analytically, and postpone the use of the two-dimensional fully numerical analysis to a subsequent work.

The main simplification in the models discussed here is that the vertical structure of the tori is accounted for by integrating the relevant quantities along the direction perpendicular to the equatorial plane. Doing so removes one spatial dimension from the problem, which can then be solved by integrating simple ordinary differential equations. The local approach, in particular, will allow us to derive the local dispersion relation obeyed by the oscillations in relativistic non-Keplerian discs, while the global approach will provide us with the eigenfunctions and eigenfrequencies of the system.

There are several motivations behind this study. First, we want to extend the relativistic discoseismology analysis carried out so far for thin discs to systems having a non-negligible contribution coming from pressure gradients. Secondly, we intend to interpret and clarify some of the numerical results found in the time evolution of ‘toroidal neutron stars’, i.e. compact and massive tori orbiting a Schwarzschild black hole (Zanotti, Rezzolla & Font 2003). Thirdly, we want to assess the possible connections between the oscillation modes of relativistic tori and the rich X-ray phenomenology observed in QPOs. Finally, we want to investigate the possibility that the axisymmetric oscillations of thick discs could provide a criterion for the development or the suppression of the runaway instability (Abramowicz, Calvani & Nobili 1983).

The plan of the paper is as follows. In Section 2 we briefly review the local analysis of oscillation modes in vertically integrated Newtonian tori. Section 3, on the other hand, introduces the basic assumptions and equations employed in the definition of our general relativistic, vertically integrated torus. These equations will then be used to study axisymmetric oscillations both locally, in Section 4, and globally, in Section 5. We will first consider configurations with constant distributions of specific angular momentum, and subsequently distributions of specific angular momentum that are either linear or power laws in the cylindrical radial coordinate. Finally, Section 8 contains our conclusions and the prospects for further investigations. Hereafter Greek indices are taken to run from 0 to 3 and Latin indices from 1 to 3; unless stated differently, we will use units in which G=c=M_⊙= 1.

2 Newtonian Tori: A Local Analysis

We here briefly present a local analysis of the oscillation modes of vertically integrated tori within a Newtonian framework. Some of the results presented in this section have already been discussed in the literature, but will serve as a useful reference for the general relativistic treatment presented in Section 4, and will help in the physical interpretation of the relativistic results.

2.1 Assumptions and equations

Consider an extended perfect fluid configuration orbiting a central object which is the only source of the gravitational potential (i.e. the orbiting fluid is non-self-gravitating). Introduce now a cylindrical coordinate system, (ϖ, φ, z) the origin of which is at the centre of the central object, and the z-axis of which is oriented along the direction of the orbital angular momentum vector.

A simplified description of this system can be obtained by removing the dependence of the various physical quantities on the vertical coordinate z. Mathematically, this is done by integrating the relevant physical quantities along the vertical direction. Physically, this corresponds to collapsing the vertical structure of the torus on to the equatorial plane, but is quite different from just considering an equatorial slice of the vertically extended torus.

Using this approach, it is possible to define the vertically integrated pressure P:  

(1)

and the vertically integrated rest-mass density Σ:  

(2)

where H=H(ϖ) is the local ‘thickness’ of the torus. An equation of state (EOS) for the fluid in the torus needs to be specified, and we find it convenient to use here a simple barotropic EOS of polytropic type, i.e.  

(3)

where k and γ≡ d ln p/d ln ρ are the polytropic constant and the adiabatic index, respectively. The vertically integrated quantities in (1) and (2) need also to be related through an ‘effective’ EOS, and this can be done if we define an ‘effective’ adiabatic index Γ≡ d ln P/d ln Σ, so that  

(4)

with the constants and Γ playing the role of the polytropic constant and of the adiabatic index, respectively. Note that while equation (4) mimics a polytropic EOS, it does not represent a vertically integrated polytropic EOS (unless, of course, the adiabatic index is equal to 1). More importantly, the adiabatic index Γ is not constant but depends on both ϖ and z. This complication, however, can be removed if one assumes that both the pressure and the rest-mass density have a weak dependence on height, so that they can be accurately expressed in terms of their values at the equator, i.e.  

(5)

 

(6)

Using this assumption, P≈ 2Hp₀ and Σ≈ 2Hρ₀, so that Γ≈ d ln p₀/d ln ρ₀=γ and equation (4) can effectively be written as  

(7)

In all of the calculations reported here we have used Γ= 4/3, but the results do not change qualitatively when different polytropic indices are used.

Being in circular non-Keplerian motion with angular velocity Ω, the torus will have velocity components only in the ϖ- and in the φ-directions, with vertical averages given by  

(8)

and  

(9)

With these assumptions, the dynamics of the torus is fully described by a set of equations comprising the Euler equations in the ϖ- and φ-directions (Shu 1992),  

(10)

 

(11)

the equation for the conservation of mass,  

(12)

and the Poisson equation for a vertically integrated gravitational potential Ψ,  

(13)

where is the vertically integrated rest-mass density of the central object, and is the only source of the gravitational potential.

Stationary and axisymmetric models for the tori can be built after setting to zero the terms in (10)–(12) involving derivatives with respect to the t- and φ-coordinates. The solutions obtained in this way represent the background equilibrium solutions over which harmonic Eulerian perturbations of the type  

(14)

can be introduced. Note that Eulerian perturbations are indicated with a ‘δ’ to distinguish them from the corresponding Lagrangian perturbations, which we will instead indicate with a ‘Δ’. In expression (14), σ is the frequency of the perturbation (in general a complex number), k_ϖ is the wavenumber in the ϖ-direction (hereafter simply k), and δq≡δP/Σ has been introduced to describe the perturbation in the pressure.

We have here neglected the perturbations in the gravitational potential and therefore set δΨ= 0. This is referred to as the Cowling approximation (Cowling 1941), and for a fluid configuration that is non-self-gravitating (i.e. one for which the gravitational potential is such that Ψ_fluid+δΨ_fluid= 0), the Cowling approximation is actually an exact description of the pulsations (Ipser & Lindblom 1992).

Note that the harmonic spatial dependence expressed in (14) is the signature of the local approach and is valid as long as the wavelength of the perturbations considered is smaller than the length-scale of the radial variations in the equilibrium configuration [this is basically the condition for the Wenzel-Kramers-Brillouin (WKB) approximation] or, stated differently,  

(15)

Introducing the perturbations (14) in the equilibrium model given by equations (10)–(13) with time-derivative terms omitted and retaining only the first-order terms, we derive the following perturbation equations (Shu 1992):  

(16)

 

(17)

 

(18)

where c²_s≡ dP/dΣ=ΓP/Σ is the local sound speed. We can now cast equations (16)–(18) into a simple matrix form as  

(19)

where κ_r is the epicyclic frequency in the radial direction¹ and is defined as  

(20)

Note that the radial epicyclic frequency is equal to the orbital frequency for Keplerian orbital motion, i.e. κ²_r=Ω² for Ω∼ϖ^−3/2, and is zero for motions with constant specific angular momentum, i.e. for ℓ≡Ωϖ²= constant.

Being a homogeneous linear system, a non-trivial solution to the set of equations (19) exists if the determinant of the coefficients matrix is equal to zero, which then provides the dispersion relation  

(21)

This approximate form of the dispersion relation was first applied to waves in accretion discs by Okazaki et al. (1987) and then reconsidered by several authors in more general situations (see Nowak & Wagoner 1992; Ipser 1994; Silbergleit et al. 2001). The two terms in the dispersion relation (21) are most easily interpreted when considered separately. To this end, consider a torus composed of collisionless particles and with specific angular momentum increasing outwards. A fluid element that is infinitesimally displaced from its equilibrium orbit but conserves its angular momentum unchanged will start oscillating in the radial direction as a result of a restoring centrifugal force. These oscillations are called inertial oscillations and their frequency is the radial epicyclic frequency κ_r (ϖ). In compressible fluids, on the other hand, a restoring force due to pressure gradients is also present and is responsible for acoustic oscillations with frequency kc_s. Both of these terms contribute to the right-hand side of the dispersion relation (21) and are collectively referred to as ‘inertial-acoustic waves’. Following a standard convention (Kato 2001), we will here identify the high-frequency modes (σ²≳κ²_r) with the p modes, or inertial-acoustic modes (Kato & Fukue 1980) of the vertically integrated torus.

3 Relativistic Tori: Assumptions and Equations

As for the Newtonian case, we will here assume that the torus does not contribute to the space-time metric, which we will take to be that external to a Schwarzschild black hole. This is the simplest, yet non-trivial, metric to consider and, more importantly, it allows for a direct comparison with the numerical calculations of Zanotti et al. (2003). [The extension of this work to the Kerr metric will be presented in a future work (Rezzolla & Yoshida, in preparation).] Furthermore, since we are really interested in the portion of the space-time in the vicinity of the equatorial plane (i.e. for values of the spherical angular coordinate |θ−π/2|≪ 1), we will write the Schwarzschild metric in cylindrical coordinates and retain the zeroth-order terms in the ratio (z/r) (Novikov & Thorne 1973). In this case, the line element takes the form  

(22)

where  

(23)

The basic equations to be solved in order to construct equilibrium models are the continuity equation for the rest-mass density ρ,  

(24)

and the conservation of energy-momentum,  

(25)

where the symbol ∇ refers to a covariant derivative with respect to the metric (22). In equation (25), T^αβ≡ (e+p) u^αu^β+pg^αβ are the components of the stress-energy tensor of a fluid with isotropic pressure p and total energy density e. Since we are dealing with a curved background space-time, but we want a close comparison with Newtonian expressions, it is useful to introduce an orthonormal tetrad carried by the local static observers and defined by the one-form basis  

(26)

In this frame, the components of the fluid four-velocity are denoted by (with μ=t, ϖ, z, φ), and the three-velocity components are then given by  

(27)

In analogy with the Newtonian case presented in Section 2, we use the vertically integrated quantities defined in (1)–(2) and assume that they obey an effective polytropic equation of state of the type (4). Enforcing the conditions of hydrostatic equilibrium and of axisymmetry simplifies the hydrodynamical equations considerably, reducing them to Bernoulli-type equations (Kozlowski, Jaroszynski & Abramowicz 1978):  

(28)

where ℓ≡−u_φ/u_t is the specific angular momentum (i.e. the angular momentum per unit energy). Simple algebraic manipulations show that the only relevant component of equations (28), i.e. the radial one, can be written explicitly as  

(29)

which represents a generic condition for the existence of hydrostatic equilibrium of a fluid configuration orbiting around a Schwarzschild black hole. Note also that a number of cancellations have removed any dependence on the spatial derivative of the angular velocity from the right-hand side of equation (29).

Being interested in vertically integrated expressions, we need to integrate both sides of equation (29). Interestingly, in the space-time (22) the right-hand side of (29) is a function of the cylindrical radial coordinate only, and hence its vertical integration simply yields a new multiplicative factor 2H. The left-hand side, on the other hand, cannot be managed unless a simplifying assumption is made, and this amounts to considering that the total enthalpy (e+p) is only a weak function of the z-coordinate so that, effectively, it is possible to substitute it with its vertically averaged value , where  

(30)

With this assumption, the vertically integrated equation (29) can be written as  

(31)

where E is the vertically integrated energy density  

(32)

We next introduce perturbations in the velocity and pressure with a harmonic time dependence of the type  

(33)

where δQ≡δP/(E+P), and where we have defined the averaged velocity perturbations as  

(34)

After perturbing the hydrodynamical equations (24)–(25) and taking the vertically integration of the resulting equations, we derive the following set of equations:  

(35)

 

(36)

 

(37)

where A≡ 1/(u_t)². Equations (35)–(37) represent the ϖ- and φ-components of the perturbed relativistic Euler equations as well as the perturbed continuity equation.

4 Perturbations of Relativistic Tori: A Local Analysis

A local analysis of the perturbed hydrodynamical equations (35)–(37) can be performed after introducing a harmonic radial dependence of the type , where k is the radial wavenumber and, again, we will assume that λ= 2π/k≪L. Simple considerations allow us now to remove some of the terms in the continuity equation (37), thus simplifying it further. In particular, it is easy to realize that the third and fourth terms in equation (37) are very small as compared with the first and second ones, to which they are similar in nature. More specifically, both the second and the fourth terms involve the radial velocity perturbation but with coefficients proportional to k and 1/L, respectively. Assuming that the condition (15) is satisfied, it is then reasonable to drop the fourth term of equation (37). Similarly, the ratio between the third and the first terms can be approximated as  

(38)

where we have approximated the first term as ∼δP/(ΓP) ∼δQ/c²_s and the third one as . While is the ratio of two perturbed quantities and of the order of unity, the typical velocity is of the order of the orbital velocity at the marginally stable circular orbit (i.e. ) so that, overall, the ratio in (38) is much smaller than unity and the third term in equation (37) can therefore also be dropped. Finally, it should be noted that the third term on the left-hand side of equation (36) can also be discarded, as it would introduce a purely imaginary part in the dispersion relation (42). This term, which is due to a time derivative of the pressure in the Euler equations, has a purely relativistic origin and has been found to provide negligible contributions in the numerical solution of the system (35)–(37).

We can now introduce quantities that are more directly related to the ones used in the Newtonian example of Section 2, namely  

(39)

so that the linearized perturbation equations in the unknowns δU, δW and δQ can be written as a homogeneous linear system:  

(40)

where the matrix of coefficients M depends both on k and σ. Imposing the condition that the determinant of M to be zero, we obtain the dispersion relation  

(41)

A non-trivial solution of (41) is given by  

(42)

where now c²_s≡ dP/dE is the square of the relativistic sound velocity in the vertically integrated disc. The dispersion relation (42) represents the relativistic generalization of (21) in which we have defined the relativistic radial epicyclic frequency for an extended fluid object κ_r as  

(43)

A number of comments should now be made about this expression. First, while the definition (43) reduces to the corresponding expression (20) when the Newtonian limit is taken, it has a new important correction coming from pressure gradients (cf. the second term in the second round brackets, which is a function of ϖ only). Secondly, in those cases in which the (radial) pressure gradients can be neglected (e.g. for very thin discs), the fluid motion is essentially Keplerian and expression (43) coincides with the relativistic radial epicyclic frequency for a point-like particle. In this case, and only in this case, (43) reduces to the expression for the radial epicyclic frequency derived by Okazaki et al. (1987) which, for a Schwarzschild space-time in spherical coordinates, is  

(44)

An important feature of the relativistic radial epicyclic frequency (44), and which distinguishes it from the corresponding expression (43), is that such a function is not monotonic in ϖ but has a maximum at a given radius, thus indicating that in general relativity oscillation modes can be trapped near the inner edge of an accretion disc (Kato & Fukue 1980). When pressure gradients are taken into account, however, this ceases to be true and the epicyclic frequency becomes monotonically decreasing with radius (cf. Fig. 5, later)

Figure 5.

The upper panel shows the relativistic sound velocity for representative tori models having a constant [solid line; model (c2)], a linear [short-dashed line; model (l2)] or a power-law [long-dashed line; model (pl4)] distribution of specific angular momentum. The lower panel, instead, shows the dimensionless radial epicyclic frequency for the three models shown in the upper panel, with the horizontal lines showing the fundamental frequencies for models (l2) and (pl4). All quantities are expressed in normalized units.

Finally, it should be noted that, as for the Newtonian case, expression (43) predicts a zero epicyclic frequency for orbital motion with constant specific angular momentum. Note that this holds true only when the specific angular momentum is defined as ℓ≡−u_φ/u_t, in which case  

(45)

To the best of our knowledge, this result has not been discussed before in the literature for finite-size fluid configurations. Furthermore, it is relevant for the physical interpretation of the modes found in the numerical calculations of Zanotti et al. (2003), as it indicates that the oscillations found in those simulations cannot be associated (at least within the validity of a vertically integrated approach) with epicyclic oscillations, but should be associated with other restoring forces, most notably pressure gradients.

5 Perturbations of Relativistic Tori: A Global Analysis

A global analysis of the axisymmetric oscillation modes of a relativistic torus consists of solving the system of equations (35)–(37) as an eigenvalue problem, treating the perturbed quantities as eigenfunctions and the unknown σ as the eigenfrequency [see Rodriguez et al. (2002) for the solution of the equivalent problem for relativistic thin discs].

It is worth remarking that equation (28) basically states that a condition for the existence of stationary, geometrically thick fluid configurations orbiting around a Schwarzschild black hole is that the distribution of angular momentum is a non-Keplerian one. As long as this condition is met, in fact, the left-hand side of equation (28) will not be zero and the pressure gradients will expand the disc in the vertical direction. Of course, different distributions of angular momentum will define tori with different properties, such as different inner and outer radii, different positions for the maximum in density, etc. Among the infinite number of possible angular momentum distributions, those that produce tori with a finite radial extension are of particular interest, since they allow for the existence of globally trapped, axisymmetric modes of oscillation.

Because of the degeneracy in the functional form for ℓ(ϖ), we have constructed several different models which we have reported in Table 1 and which we will be discussing in detail in the following sections. This approach to the problem is clearly more expensive as the eigenvalue problem needs to be solved for a large number of models; however, it also allows one to investigate how the axisymmetric oscillations depend on the angular momentum distribution. When features are found that do not depend on a specific choice for the distribution of angular momentum and that are thus indication of a universal behaviour, this approach can be very rewarding, as we will discuss in the following sections.

Table 1.

Main properties of the equilibrium considered. From left to right the columns report: the type of angular momentum distribution (see Section 7 for a description of , α, β and q), the inner and the outer radii of the torus and , the radial position of the maximum density in the torus , and the value of the potential gap chosen at the inner edge of the torus.

However, before discussing the results of the numerical solution of the eigenvalue problems, it is useful to review briefly the numerical approach followed in the solution of the eigenvalue problem and the boundary conditions imposed.

5.1 Numerical method

Equations (35)–(37), supplemented with suitable boundary conditions at both edges of the numerical grid, represent a standard two-point boundary value problem, for which a ‘shooting’ method can be used. In practice, given a trial value for the eigenfrequency σ, two solutions for the unknown quantity (either δQ or δU) are found starting from the two edges of the numerical grid at the inner, ϖ_in, and outer, ϖ_out, radii of the torus. These two solutions are then matched at an arbitrary point ϖ_M in the domain, where the Wronskian of the left and right solutions is evaluated. This procedure is iterated until a zero of the Wronskian is found with the desired accuracy, thus determining both the eigenfrequency and the eigenfunctions of the mode considered.

Note that, because of the linearity of the present analysis, the solutions are invariant after a rescaling of all the quantities by a constant. In practice, this freedom allows one to choose the initial guess for one of the eigenfunctions at one edge of the numerical grid to be a freely specifiable number, which we typically set to 1. Furthermore, since the space-time is uniquely described once the mass M of the black hole has been fixed, we can use it to rescale all of the relevant quantities involved in the calculations. The new dimensionless quantities, indicated with a tilde, are then defined as  

(46)

 

(47)

 

(48)

 

(49)

 

(50)

 

(51)

Unless specified differently, all of the results presented hereafter will be given in terms of these dimensionless quantities.

5.2 Boundary conditions

Suitable boundary conditions need to be specified at the inner and outer edges of the vertically integrated torus for the solution of the set of equations (35)–(37). As for oscillations in stars, we here assume the perturbed surface of the torus to be where the Lagrangian pressure perturbation vanishes, i.e. where  

(52)

Because the Eulerian and Lagrangian descriptions are linked through the relation  

(53)

where is the Lie derivative along the Lagrangian displacement three-vector , condition (52) can also be written as  

(54)

from which it follows that Δp= 0 =δp for a polytropic fluid configuration the rest-mass density of which vanishes at the surface (Tassoul 1978). Note also that, while in the Cowling approximation the Eulerian perturbation of the metric, δg_ab, is zero, the Lagrangian perturbation is, in general, non-zero.

Using the relation between the perturbed four-velocity u^α and the Lagrangian perturbation of the metric (see Friedman 1978),  

(55)

which, in our case, reduces to the relation  

(56)

we can write the boundary condition (52) simply as  

(57)

Expression (57), taken as a linear relation between δQ and δU, provides us with the boundary conditions to be imposed at the inner and outer edges of the torus.

6 Constant Specific Angular Momentum Tori

The simplest and most studied choice for the distribution of angular momentum is the one in which ℓ=ℓ_c= constant, with ℓ_c being chosen in the interval between the specific angular momentum at the marginally stable orbit and the one at the marginally bound orbit ℓ_mb= 4, to yield a disc of finite size. Besides yielding a simpler approach, constant specific angular momentum tori benefit from having a radial epicyclic frequency which is identically zero [cf. equations (43) and (45)], thus leaving pressure gradients as the only possible restoring forces. As a result, a constant angular momentum distribution provides an interesting limit for the properties of p-mode oscillations that will not be influenced by centrifugal effects.

The procedure for the construction of the background model for the torus can be summarized as follows. First, after a value ℓ_c for the constant specific angular momentum has been chosen, the locations of the cusp ϖ_cusp and of the maximum density in the torus²ϖ_max are determined as the positions at which the Keplerian specific angular momentum is equal to ℓ_c. The inner and outer edges of the torus can then be calculated after requiring that they lie on the same equipotential surface, i.e.  

(58)

where ΔW_in is the value of the potential gap chosen at the inner edge of the torus and is therefore a measure of how much the torus ‘fills’ its Roche lobe. Clearly, the case in which ΔW_in= 0 corresponds to a torus filling its outermost closed equipotential surface (i.e. the one possessing a cusp).

Having determined ϖ_in, it is then possible to calculate the angular velocity distribution as  

(59)

where Ω_c is angular velocity at the radial position of cusp ϖ_cusp. Finally, after making a choice for the parameters of the EOS (and therefore for k and γ), the equation for the hydrostatic equilibrium (31) can be integrated, either numerically or analytically (Kozlowski et al. 1978), to build the torus model which will serve as background for the introduction of the perturbations.

When the specific angular momentum is constant within the torus, the perturbation equations also become simpler and, in particular, the second term in (36) vanishes exactly (cf. equation 45) and the resulting set of equations can then be written as  

(60)

 

(61)

 

(62)

Equations (60)–(62) have been solved numerically for a number of different models, the main properties of which, such as the inner and outer radii and , the position of the maximum in the surface density distribution , and the value of the specific angular momentum at the inner edge of the disc ℓ_c≡ℓ(ϖ_in), are summarized in Table 1. As a representative example, we present in Figs 1 and 2 the results for a torus model with ℓ_c= 3.75[i.e. model (c2) of Table 1].

Figure 1.

Eigenfunctions for δQ=δP/(E+P) as a function of the radial coordinate for a constant angular momentum disc. Only the fundamental mode f and the first three overtones, denoted by o₁, o₂ and o₃, have been reported. The data refer to model (c2) of Table 1 and the units on the vertical axis are arbitrary.

Figure 2.

Eigenfunctions for as a function of the radial coordinate for a constant angular momentum disc. Only the fundamental mode f and the first three overtones, denoted by o₁, o₂ and o₃, have been reported. The data refer to model (c2) of Table 1 and the units on the vertical axis are arbitrary.

Fig. 1, in particular, shows the solution for the eigenfunction δQ both for the fundamental mode f and for the first three overtones o₁, o₂ and o₃. Similarly, Fig. 2 shows the fundamental mode and the first three overtones for the eigenfunction δU, which has always one node less than δQ. Note that the choice of a very small potential barrier introduces complications in the numerical solution of the eigenvalue problem. For a model with ΔW_in= 0, in fact, the fluid elements at the inner edge are just marginally stable to accretion on to the black hole. This behaviour is reflected in the eigenfunctions δQ and δU, which tend to diverge at the inner and outer edges of the torus. To avoid this problem (still partially visible in Figs 1 and 2) we have used a small but non-zero value, ΔW_in=−10⁻⁵, for the potential barrier.

A careful look at Figs 1 and 2 reveals that the modes f and o₂ have similar properties at the inner edge; furthermore, during each half-period, they have signs which are opposite to those of the modes o₁ and o₃. This is an important feature, indicating that not all modes will behave in the same manner at the inner edge where mass loss and accretion on to the black hole take place. While a conclusion on the this behaviour cannot be drawn on the basis of the present linear analysis, in which the total mass is conserved exactly ( at the inner and outer edges for all modes), the functional form of the eigenfunctions for δQ and δU seems to suggest that mass accretion through the cusp will be possible preferably at the frequencies corresponding to the f and o_2n modes [see also the discussion below for a comparison with the numerical calculations presented in Zanotti et al. (2003)].

The eigenfrequencies corresponding to model (c2), as well as for all of the other models in Table 1, are listed in Table 2, where we have included the frequencies up to the fourth overtone. A careful look at Table 2 reveals that the computed eigenfrequencies are, at least for the first few modes, in a sequence 2:3:4:…, to a good precision. As the order of the mode increases, this simple sequence is no longer followed and a different one appears. The fact that the lowest order modes appear at frequencies that are in small integer ratios is not particularly surprising if p modes are to be interpreted as sound waves trapped within the cavity represented by the confined torus. In this case, in fact, one would indeed expect that modes should be trapped and with wavelengths that are multiples of the trapping length-scale, i.e. λ_n=[2/(2 +n)]L, where n= 0, 1, …. Stated differently, the results reported in Table 2 are consistent with the idea that p modes can be associated with sound waves trapped in the torus and having frequencies that are multiples of the fundamental one and given by σ_n=c_s/λ_n=[(2 +n)/2]f.

Table 2.

Eigenfrequencies of the fundamental f mode and of the first four overtones o_n of the vertically integrated models described in Table 1. All the frequencies are given in normalized units.

Note that while these p modes follow the same definition and are similar in nature to the ones discussed by Nowak & Wagoner (1991, 1992) for thin discs, their trapping is not produced by a turning of the relativistic radial epicyclic frequency at smaller radii. Rather, the relativistic radial epicyclic frequency in discs in which pressure gradients play an important role has been found to be monotonically decreasing for all of the disc models considered here [cf. equation (43), Fig. 5, later]. An important consequence of this is that the p modes for thick discs are not restricted to be trapped in special parts of the disc, as happens for the corresponding p modes for thin discs (Nowak & Wagoner 1991, 1992). In particular, they need to be present only in the (rather small) inner regions of the disc, where they can produce only a slight modulation in the disc emission. Similarly, they are not constrained to be trapped in the outer regions of the disc, where the location of the outer radius remains uncertain and the corresponding frequencies very small. On the contrary, the modes discussed here are present over the whole torus, which then behaves as a single trapping cavity. This is a substantial difference, giving these modes the possibility of possessing frequencies comparable to observations and of causing non-negligible modulations in the disc emission.

Having solved the eigenvalue problem for p modes in vertically integrated relativistic tori with constant specific angular momentum, we can now proceed to a comparison with the fully non-linear two-dimensional simulations discussed by Zanotti et al. (2003). This is presented in Fig. 3, where we have plotted the power spectra of the L₂ norm of the rest-mass density for the torus model (c2) in Table 1. The solid line, in particular, refers to initial data in which a global perturbation has been introduced in terms of the radial velocity field of a relativistic spherical accretion solution [cf. equation (15) of Zanotti et al. (2003)]. Clearly, the power spectrum for the continuous line shows the presence of peaks appearing in a simple small integer sequence 2:3:4:…, providing convincing evidence that the oscillations triggered in the simulations of Zanotti et al. (2003) correspond indeed to p-mode oscillations of perturbed relativistic tori. It should also be noted that this sequence is not the same as shown by the peaks in the power spectrum of the mass accretion rate on to the black hole, which instead are only in the sequence 1:2:3:…[cf. fig. 7 of Zanotti et al. (2003)]. This difference remains puzzling, but could be explained on the basis of the behaviour of the eigenfunctions for the f and o_2n modes at the inner edge discussed above. In this case, in fact, all of the o_2n+1 modes would have vanishingly small perturbations at the inner edge of the disc and would not produce a mass accretion at those frequencies, as shown by fig. 7 of Zanotti et al. (2003).

Figure 3.

Power spectra of the L₂ norm of the rest-mass density for the torus model (c2) in table 1 of Zanotti et al. (2003). Different lines refer to models with different initial perturbations. The solid line, in particular, refers to initial data in which a global perturbation has been introduced in the radial velocity field only (cf. equation 15). The dashed line, on the other hand, refers to initial data in which a perturbation based on the computed eigenfunctions of the o₁ mode has been introduced both in the radial velocity and in the density. The two spectra have been rescaled in order to match in the power of the fundamental frequency.

The knowledge of the eigenfunctions of p modes in relativistic tori can also be used to go a step beyond a simple comparison with numerical simulations and actually use the latter to perform investigations of ‘numerical’ discoseismology. In other words, having a complete picture of the oscillation properties and an accurate two-dimensional numerical code to simulate them, it is possible to investigate the dynamical response of a relativistic thick disc as a result of the introduction of suitably selected perturbations. As a concrete example, we report with the dashed line in Fig. 3 the power spectrum of the L₂ norm of the rest-mass density for the torus model (c2) in Table 1. The important difference from the corresponding spectrum indicated with a solid line is that the dashed line refers to a simulation having as initial data perturbations based on the computed eigenfunctions for δU and δQ of the o₁ mode. This represents a selective excitation of the o₁ mode and, as a result, the corresponding power in the first overtone is increased by almost a factor of 10 (the power in the o₂ mode is instead decreased by a similar amount). This behaviour further confirms that the p modes derived in this perturbative analysis correspond to the modes simulated numerically by Zanotti et al. (2003). Finally, it should be noted that the power spectrum corresponding to initial data containing only an o₁-mode perturbation shows that the other modes have also been excited, most notably the fundamental one. This is due partly to a mode-mode coupling which transfers energy from one mode to the other ones, but also to the error introduced using eigenfunctions derived for a vertically integrated model, and which cannot reproduce the corresponding eigenfunctions of a two-dimensional disc exactly.

There is a final aspect of the axisymmetric p modes discussed so far that is worth underlining. This is illustrated in Fig. 4, where we have plotted the values of the eigenfrequencies for the fundamental mode f (solid line) and for the first overtone o₁ (dashed line) for a large number of tori. All of the points on the solid and dashed curves in Fig. 4 represent the solution of the eigenvalue problem and the two sequences have been calculated for tori with fixed ‘centre’, but having different radial extents . As one would expect for modes behaving effectively as sound waves, the eigenfrequencies decrease as the radial extent of the torus increases. Furthermore, the frequencies for f and o₁ maintain a harmonic ratio 2:3 for all the values of . Most importantly, however, the filled dot shown in Fig. 4 for represents the radial epicyclic frequency for a circular orbit at , i.e. .

Figure 4.

Eigenfrequencies for the fundamental f (solid line) and for the first overtone o₁ (dashed line) of axisymmetric p modes for ℓ= constant tori. The two lines refer to sequences of tori having the same radial position for the ‘centre’, but different radial extents . The filled dot represents the radial epicyclic frequency for a circular orbit at and thus the value at which the fundamental p-mode frequency tends in the limit of a vanishing torus size.

This is an important result which provides a simple interpretation of the nature of the axisymmetric p modes of thick discs: the radial epicyclic frequency at the position of maximum density represents the value at which the fundamental p-mode frequency tends to the limit of a vanishing torus size. The importance of this result is that, while different models for the tori, using, for instance, different equations of state or different angular momentum distributions, will produce different slopes for the solid and dashed lines in Fig. 4, all of the sequences will terminate on the filled dot when L→ 0, that is, when finite-size effects will cease to be relevant and each torus will effectively behave as a ring of particles in circular orbits. Because this occurs only in the limit L→ 0, it should not be surprising that it is valid also for sequences of tori with constant distributions of specific angular momentum for which, as discussed in Section 4, the epicyclic frequency is zero. Indeed, this result is totally independent of the distribution of specific angular momentum and will be confirmed in the following section, where tori with non-constant specific angular momentum distributions are considered. An analytic proof of this conclusion is provided in Appendix A.

7 Non-Constant Specific Angular Momentum Tori

While the study of constant angular momentum discs offers several advantages and simplifies the equations, realistic discs are likely to have angular momentum distributions that are not constant. For this reason, and in order to assess the validity of the results derived so far with more generic distributions of angular momentum, we have extended our mode analysis also to discs with non-constant angular momentum. The first step in this direction consists of determining which distribution of specific angular momentum ℓ=ℓ(ϖ) should be specified in the construction of an equilibrium model. This choice is, to some extent, arbitrary with the only constraint being given by Rayleigh's criterion for the dynamical stability against axisymmetric perturbations (Tassoul 1978). This condition, however, is not very strong and simply requires that dℓ/dϖ≥ 0, so that even a constant angular momentum distribution is stable, although only marginally.

Some guidance in this choice can come from numerical simulations, and indeed calculations of torus formation performed by Davies et al. (1994) with smooth particle hydrodynamics (SPH) techniques have shown that the final configuration consists of a core object surrounded by a torus the angular momentum distribution of which is close to a power law in radius, with index 0.2. In our computations we have therefore adopted two simple and different distributions for the specific angular momentum that can be expressed analytically. We refer to the first one as the ‘linear’ distribution, in which  

(63)

where α and β are adjustable constant coefficients. Similarly, we refer to the second one as the ‘power-law’ distribution, in which instead  

(64)

where q is also an arbitrary parameter. The values chosen for and q for the representative models discussed here have been summarized in Table 1. Note that we have deliberately chosen small values for both α and q, and this is to avoid the construction of unperturbed models that differ significantly in radial extent from the ones built with a constant specific angular momentum distribution.³

The construction of the equilibrium models proceeds in this case as discussed in the previous section, namely, by first fixing an initial value for the specific angular momentum at the cusp, then by calculating the radial extent of the torus in terms of the potential barrier ΔW_in, and finally by integrating numerically equation (29) for the chosen distribution of Ω(ϖ).

The full system of the perturbed equations (35)–(37) can now be rewritten in a more compact form as  

(65)

 

(66)

 

(67)

and can be solved by adopting the same numerical technique as used for the constant angular momentum models.

The upper panel of Fig. 5 shows the relativistic sound velocity for representative unperturbed tori models having a constant (solid line), a linear (short-dashed line) or a power-law (long-dashed line) distribution of specific angular momentum [the data refer to models (c2), (l2) and (pl4) of Table 1, respectively]. Note that, while all of the curves refer to models with a polytropic EOS, each curve does not depend on the value chosen for the polytropic constant k and hence on the mass of the disc (see Appendix B for a proof of this). The lower panel of Fig. 5, on the other hand, shows the radial epicyclic frequency κ_r as calculated from expression (43) for the three models shown in the upper panel (we recall that κ_r= 0 when ℓ= constant). The horizontal lines in the lower panel refer to the fundamental frequencies for models (l2) and (pl4), and their relative position with respect to the curves for κ_r will be important for the appearance of an evanescent-wave region in the inner parts of the torus (see Fig. 6 and the discussion below).

Figure 6.

Schematic propagation diagram for a thick disc. The thick solid line shows the values of the relativistic radial epicyclic frequency in the inner regions of the disc, while the horizontal dashed line shows the value of the fundamental frequency f. The data are presented in normalized units and refer to model (pl2) of Table 1. See the main text for a discussion.

We have summarized in the different panels of Fig. 7 the eigenfunctions for δQ, δU and δW for tori with a linear [panels (a)–(c)] and a power-law [panels (d)–(f)] distribution of specific angular momentum. The eigenfunctions refer, in particular, to models (l1) and (pl3) of Table 1 and are shown in the fundamental as well in the first three overtones. Similarly, in Table 2 are also summarized the eigenfrequencies computed for a number of tori with non-constant specific angular momentum.

Figure 7.

Eigenfunctions for and as a function of the radial coordinate for tori with a ‘linear’[panels (a)–(c); model (l1)] and a ‘power-law’[panels (d)–(f); model (pl3)] distribution of the specific angular momentum. While the units used on the vertical axes are arbitrary, it should be noted that the values in panel (e) are one order of magnitude smaller than the corresponding ones in panel (b). In all panels only the fundamental mode f and the first three overtones have been reported.

A rapid analysis of the eigenfrequencies in Table 2 as well as of the eigenfunctions in the various panels of Fig. 7 indicates that, overall, the results for non-constant distributions of specific angular momentum do not differ, at least qualitatively, from those computed for constant distributions. Most importantly, the eigenfrequencies continue to appear in the simple sequence 2:3:4:…, at least for the lowest-order modes considered here.

This result may appear surprising in light of the fact that the the radial epicyclic frequency defined by (43) does not vanish in the case of non-constant specific angular momentum distributions, so that the computed eigenfrequencies contain also a contribution coming from the inertial oscillations [cf. the dispersion relation (42)]. However, it is important to bear in mind that the p modes discussed here behave essentially as sound waves trapped inside the torus. As a result, while their absolute frequencies will be different when tori with different radial extensions are considered, the sequence in which they appear will basically follow the relation σ_n=[(2 +n)/2]f.

The ‘universality’ in which the lowest order p-mode seem to appear has an interest of its own, but may also have important implications for the high-frequency QPOs observed in binaries containing a black hole candidate (Strohmayer 2001; Remillard et al. 2002). In at least three systems, in fact, the QPOs have been detected with relatively strong peaks that are in a harmonic ratio of small integers 1:2, 2:3 or 1:2:3 (Abramowicz & Kluzniak 2001), and a possible explanation of this phenomenology in terms of the oscillations of a geometrically thick, non-Keplerian disc of small radial extent orbits around the black hole will be presented in a separate paper (Rezzolla et al. 2003).

Before concluding this section, we should underline an important qualitative difference that emerges when non-constant distributions of specific angular momentum are considered. This is related to the behaviour of the eigenfunctions for the radial velocity perturbation δU at the inner edge of the disc. A rapid comparison of Fig. 2 with the corresponding panels (b) and (e) of Fig. 7 shows that the eigenfunctions become vanishingly small at the inner edge of the disc in the case of a ‘power-law’ distribution of specific angular momentum. This does not seem to be case for the behaviour at the outer edge, where the eigenfunctions maintain very large values. Similarly, this behaviour at the inner edge is not encountered with the other two distributions considered [note that the values the on vertical axis of panel (e) are one order of magnitude smaller than the corresponding ones on panel (b) and in Fig. 2].

A physical explanation for this behaviour can be found by comparing the behaviours of the radial epicyclic frequency curves for the different distributions of angular momentum. The lower panel of Fig. 5, in particular, shows that the radial epicyclic frequencies for both the linear and the ‘power-law’ distributions of specific angular momentum are monotonically decreasing for increasing radii. (This has been found to hold for all of the disc models considered here.) However, while the two frequencies tend to be comparable in the outer regions of the torus, they differ considerably in the inner regions, where the epicyclic frequency for a power-law distribution increases inward much more rapidly. More importantly, the epicyclic frequency is always smaller than the fundamental frequency for a torus with a linear distribution of angular momentum, while this is not the case for a power-law distribution (cf. Table 2).

In this latter case, then, there will be regions of the disc in which σ²−κ²_r∼k²c²_s < 0 and thus where the amplitude of any acoustic wave is forced to be zero. [Other modes such as the g modes can, however, exist in this region, in which they are actually trapped (Kato & Fukue 1980).] Such a region is referred to as an evanescent-wave region and effectively represents the part of the disc in which the centrifugal barrier produced by rotation makes the radial epicyclic frequency so high that all inward-moving acoustic waves cannot propagate but are reflected back to larger radial positions. In those parts of the disc where the fundamental frequency is larger than the radial epicyclic one, on the other hand, σ²−κ²_r > 0, and sound waves can propagate undamped at least in a perfect fluid. Stated differently, the results found here indicate that, when acoustic waves are considered, the inner regions of the torus are of evanescent type for discs with a power-law distribution of angular momentum, while no evanescent region exists in a constant or in a linear distribution of specific angular momentum.

Finally, it is worth pointing out that the behaviour discussed above can be used to interpret some of the most recent results on the occurrence of the runaway instability (Font & Daigne 2002a, b; Zanotti et al. 2003). In a series of relativistic hydrodynamics calculations of a torus orbiting and accreting on to a rotating black hole, in fact, Font & Daigne (2002b) have found that, while tori with constant specific angular momentum distributions lead to a runaway instability, a slight outward increase as a power law of the specific angular momentum has a dramatic stabilizing effect, suppressing the instability completely. While this behaviour reflects a process which is fundamentally non-linear, it is indeed consistent with the interpretation that the inner regions of thick discs prevent the inward propagation of perturbations and hence tend to suppress the accretion of mass on to the black hole. Because the runaway instability is fed by the changes in the space-time produced by the mass accretion, a considerable reduction of the latter can be the cause of the suppression of the former. A more detailed investigation is necessary to validate this hypothesis.

8 Conclusions

We have presented the first investigations of the oscillation properties of relativistic, non-self-gravitating tori orbiting around a black hole. More precisely, we have here considered the axisymmetric oscillation modes of a geometrically thick disc constructed in a Schwarzschild space-time. Like relativistic stars, relativistic tori are extended objects with non-trivial equilibrium configurations governed by the balance between gravitational, centrifugal and pressure forces. Unlike relativistic stars, however, the contributions coming from centrifugal forces are not small in relativistic tori, and this amplifies the role of oscillations restored by centrifugal forces. This leads to oscillation properties that can be considerably different from those encountered in stars.

Because thick discs have angular momentum distributions that are intrinsically non-Keplerian, we have modelled them with a number of different distributions of specific angular momentum, which have been chosen either to be constant within the torus, or to have a radial dependence that is linear or a power law. Furthermore, in order to keep the treatment as simple as possible and handle equations analytically whenever possible, we have built the models for the tori using vertically integrated and vertically averaged quantities.

Rather little is still known about the oscillation modes of relativistic tori, and our approach to the problem has therefore progressed by steps. In particular, our first step has been that of considering the dispersion relation for acoustic waves propagating within these objects. We have done this first in a Newtonian framework in which equations are simpler and so is their physical interpretation. We have then extended the local analysis to a general relativistic framework and in particular to a Schwarzschild black hole space-time. While a local analysis and the resulting dispersion relation are valid only for oscillations with wavelengths small when compared with the typical length-scale in the tori, the ones derived here have been important to clarify the relation between the acoustic waves and the other waves that play a fundamental role in fluids orbiting in a central potential, i.e. inertial (or epicyclic) waves. In particular, it has been possible to show that, in general, both acoustic and inertial oscillations are present in the perturbations of thick discs, and that the latter can be removed only in the special case of a constant distribution of specific angular momentum. This result, which was already known in Newtonian physics, has here been shown to hold also for extended relativistic fluid configurations. To the best of our knowledge this result has never been discussed in the literature before.

Going a step further and beyond a local analysis, we have used the same mathematical setup to perform a global analysis and determine both the eigenfunctions and the eigenfrequencies of the axisymmetric oscillations. Also in this case, the assumption of a vertically integrated equilibrium has simplified the eigenvalue problem considerably, translating it into a set of coupled ordinary differential equations which have been solved using standard techniques and for a number of different models for the tori. The modes found in this way correspond to the p modes of relativistic tori and are characterized by eigenfrequencies appearing in a simple sequence of small integers 2:3:4:…, at least for the first lower order modes. This important feature does not depend on the distribution of angular momentum used to build the tori, and basically reflects the fact that p modes are to be interpreted as sound waves trapped within the cavity represented by the confined torus.

The properties found here with a linear global analysis are in very good agreement with the numerical results found in the time evolution of perturbed ‘toroidal neutron stars’ having constant distributions of specific angular momentum (Zanotti et al. 2003). This agreement with fully non-linear two-dimensional simulations provides convincing evidence that the assumption of a vertically integrated equilibrium has not subtracted important information from the present analysis.

A number of motivations are behind the investigations carried out in this paper. First, and as mentioned above, relativistic tori possess a rich variety of oscillation modes which are still essentially unexplored, in stark contrast with what is known (both in Newtonian and in relativistic regimes) about geometrically thin discs. In this respect, the present work was intended as a first investigation of the discoseismology of geometrically thick discs. Secondly, our analysis was aimed at interpreting and clarifying some of the results found by Zanotti et al. (2003), who had investigated numerically the dynamical response of massive tori to perturbations but had not been able to characterize all of the mode properties. Finally, by investigating the response of tori to perturbations, we intended to highlight possible connections between the oscillation modes of these objects and all those astrophysical scenarios in which geometrically thick discs may play an important role.

Some of the results discussed here could find application in more general contexts. One such application is offered by the quasi-periodic X-ray phenomenology observed in X-ray binary systems containing a black hole candidate (Strohmayer 2001; Remillard et al. 2002). In these systems, in fact, the X-ray emission is modulated quasi-periodically with power spectra having peaks in a harmonic ratio of small integers 1:2, 2:3 or 1:2:3 (Abramowicz & Kluzniak 2001). While there are a number of possible explanations for this behaviour, it could be interpreted simply in terms of the p-mode oscillations of a geometrically thick, non-Keplerian disc of small radial extent orbiting around the black hole (Rezzolla et al. 2003). Another such application is offered by the runaway instability which has recently been considered by a number of authors (Font & Daigne 2002a, b; Zanotti et al. 2003) through fully non-linear numerical calculations. In particular, the numerical evidence that the instability can be suppressed if the torus has a distribution of angular momentum that is slightly increasing outwards can be interpreted simply in terms of the behaviour of the eigenfunctions for the radial velocity perturbations, which become vanishingly small in the inner regions of the disc. This, in turn, is produced by the appearance of a region of evanescent-wave propagation in the inner regions of the torus that prevents the inward propagation of perturbations, and reduces the accretion of mass on to the black hole, thus suppressing the instability.

As mentioned in the Abstract, the present work represents the first of a series of papers aimed at a systematic investigation of the oscillation properties of relativistic tori. We are presently investigating extensions of the present approach in which the p modes are investigated both when the background space-time is that of a Kerr black hole (Rezzolla & Yoshida, in preparation), and when deviations from axisymmetry are present. In addition to this, work is also in progress to extend the solution of the eigenvalue problem to a fully two-dimensional model for the torus. The results of these investigations will be presented in future papers.

Acknowledgments

It is a pleasure to thank N. Andersson, O. Blaes, T. Font, W. Kluzniak and J. Miller for many useful discussions. Financial support for this research has been provided by the MIUR and by the EU Network Programme (Research Training Network Contract HPRN-CT-2000-00137). SY is grateful for support from the grant distributed to Professor Y. Eriguchi at the University of Tokyo from the 21st Century COE program of JSPS. The computations were performed on the Beowulf Cluster for numerical relativity ‘Albert100’, at the University of Parma.

¹

In principle, an epicyclic oscillation takes places also for motions away from the background orbital plane. The corresponding frequencies are the vertical epicyclic frequencies but these will not be considered here.

²

Note that the position of the maximum density in the torus is traditionally referred to as the ‘centre’ of the torus, although its coordinate centre is, of course, at ϖ=z= 0.

³

It is not difficult to realize that a constant specific angular momentum distribution yields the most compact tori among those having the same angular momentum at the cusp.

Appendix

Appendix A: On the Values of the Eigenfrequencies when L → 0

In the limit of vanishing torus sizes, L→ 0, all of the eigenfunctions can be assumed to have only a simple linear dependence on ϖ, i.e.  

(A1)

 

(A2)

 

(A3)

where a, b, …, f are just constant coefficients. Using this ansatz in the system of perturbed equations (35)–(37), we obtain the following system of equations:  

(A4)

 

(A5)

 

(A6)

where  

(A7)

 

(A8)

and where it is easy to realize that the radial epicyclic frequency is then simply given by [cf. equation (43)]  

(A9)

Each of the equations in the system (A4)–(A6) can be written symbolically as α_iϖ+β_i= 0, i= 1, 2, 3 with α_i and β_i being just shorthands for the longer expressions in (A4)–(A6). Furthermore, because each equation in (A4)–(A6) must hold for any ϖ, we need to impose that all of the coefficients α_i and β_i vanish simultaneously. This generates a system of six equations in the six unknowns a, b, …, f:  

(A10)

 

(A11)

 

(A12)

 

(A13)

 

(A14)

 

(A15)

which, in matrix form, can be written as  

(A16)

A non-trivial solution to the system (A16) is possible if the determinant of the matrix A vanishes and, after lengthy but straightforward calculations, this condition yields  

(A17)

Since ΓP/E+P∼c²_s, all of the terms proportional to this quantity can be neglected in (A17), which then reduces to  

(A18)

A non-trivial solution of equation (A18) is then given by  

(A19)

thus proving that the eigenfrequencies tend to the value of the epicyclic frequency in the limit ot L→ 0.

Appendix B: On the Sound Speed in Rotating Polytropic Tori

We here show that a non-self-gravitating, polytropic fluid configuration orbiting in hydrostatic equilibrium around a Schwarzschild black hole has a sound velocity which is invariant under changes of the polytropic constant k. To prove this, consider the equation of hydrostatic equilibrium (28) which, in the simplified metric (22), can be rewritten as  

(B1)

The right-hand-side of equation (B1) depends only on the kinematic and geometric properties of the disc, namely on the angular momentum distribution (through Ω) and on the external gravitational field of the central black hole (through ν). Both of these quantities depend on ϖ only.

Introducing the relation γ= 1 + 1/n, the polytropic EOS p=kρ^γ can also be written in terms of the Emden function Θⁿ≡ρ as  

(B2)

so that equation (B1) effectively becomes  

(B3)

where we have defined ψ≡k(n+ 1)Θ.

Equation (B3) can be integrated analytically to give  

(B4)

where is an arbitrary constant.

The important result contained in equation (B4) is that ψ is a function of ϖ only and will not depend, therefore, on the specific value chosen for k. As a result, once γ (and thus n) is fixed, any transformation K→Kα, where α is an arbitrary constant, must be accompanied by a corresponding transformation Θ→Θ/α. An important consequence of this is that all of the quantities given by p/ρ, dp/dρ, p/e and dp/de are invariant under changes in k. Equally invariant is the sound speed defined as .

To appreciate the physical implications of this result, it is useful to recall that the polytropic constant plays here the role of determining the total rest mass of the torus. Once the other parameters of the background model have been fixed (i.e. the distribution of angular momentum and the potential gap), in fact, the rest mass of the torus is calculated from the integral  

(B5)

As a consequence, different choices for k will yield tori with different density distributions Σ and hence rest masses, while maintaining the same radial extent. The result contained in equation (B4) states, therefore, that it is possible to build tori with largely different rest-mass densities and yet have them all share the same sound velocity.

Two additional remarks are worth making. The first one is that while this result has been derived for a disc with a non-zero vertical extent, it applies also for a vertically integrated model. The second remark is that, while we have been here concentrating on relativistic models, the same result can be shown to apply in Newtonian physics when the central black hole is replaced by a generic spherically symmetric gravitational potential.