Virial tests for post-Newtonian stationary black-hole–disk systems 5 Piotr Jaranowski 1 0 Wydzial Fizyki, Uniwersytet w Bia lymstoku, Cio lkowskiego 1L, 15–245 Bia lystok, Poland 2 E-mail: [email protected] n a J Patryk Mach, Edward Malec, Micha l Pir´og 9 InstytutFizyki im. Mariana Smoluchowskiego, Uniwersytet Jagiellon´ski, L ojasiewicza 11, ] 30-348 Krak´ow, Poland c q E-mail: [email protected], [email protected], [email protected] - r g Abstract. We investigated hydrodynamical post-Newtonian models of selfgravitating [ stationary black-hole–disk systems. The post-Newtonian scheme presented here and also in our recent paper is a continuation of previous, purely Newtonian studies of selfgravitating 1 v hydrodynamical disks rotating according to the Keplerian rotation law. The post-Newtonian 4 relativisticcorrectionsaresignificantevenatthe1PNlevel. The1PNcorrection totheangular 0 velocitycanbeoftheorderof10%ofitsNewtonianvalue. Itcanbeexpressedasacombination 1 of geometric and hydrodynamical terms. Moreover, in contrast to the Newtonian Poincar´e– 2 Wavre theorem, it depends both on the distance from the rotation axis and the distance from 0 the equatorial plane. . 1 Inthetechnicalpart of thisnote wederivevirial relations valid up to1PN order. Weshow 0 that they are indeed satisfied by our numerical solutions. 5 1 : v i X 1. Introduction r Inarecent paper[1]weinvestigated post-Newtonian modelsofselfgravitating gaseous disksthat a rotate according to the Keplerian rotation law. The analysis presented there is a continuation of our previous studies, where such disk systems were investigated in Newtonian theory [2, 3]. In the Newtonian framework we posed the following problem: Suppose one observes a selfgravitating stationary gaseous disk around a central object (modeled by a point-mass) that rotates according to the Keplerian rotation law, that is with the angular velocity vφ = ω /r3/2, 0 0 where r is the distance from the rotation axis, and ω is a constant. For the disk consisting of 0 test particles we have ω = √GM , where M is the mass of the central object, and G is the 0 c c gravitational constant. What does the observed value of ω correspond to in the case where 0 the mass of the disk is comparable with M ? Is it still the central mass M , the sum of the c c two masses, or some nontrivial combination of them? It turns out that the selfgravity speeds up the rotation of the disk—it rotates faster than this would follow from the Keplerian formula involving the central mass M only. Moreover, the way in which ω depends on the central c 0 mass and the mass of the disk is prescribed by the geometry of the disk. Thus, in principle, it is possible to measure the masses of Keplerian disks, whenever their geometry is known. This procedure was applied to the accretion disk in the AGN of NGC 4258, where the Keplerian rotation curve was measured in the maser emission [3]. In [1] we extended the Newtonian analysis to the first post-Newtonian approximation (1PN). Selfgravitating stationary gaseous disks were investigated before in full relativity (cf. [4, 5]). We decided on the post-Newtonian scheme, because of its conceptual simplicity. In particular, the notion of the Keplerian rotation has a clear meaning in the post-Newtonian scheme. The main result obtained in the 1PN approximation is that the angular velocity profile is affected in two different ways—some parts of a disk can be speeded up and the others slowed down. Furthermore, the sum of the Newtonian and post-Newtonian components of the angular velocity is not anymore a function of the cylindrical radius only, but in general it depends on radial and vertical coordinates [1]. In this paper we sketch briefly the main equations that constitute the 1PN model and then derive virial-type relation that can be used to test the obtained numerical solutions. Suitable virial tests valid in the Newtonian case were presented in [6] and [2]. We discuss them here for clarity. The post-Newtonian virial identities given here are new and have not been discussed in [1]. In the last section of this paper we also show that they are satisfied by our numerical models with the accuracy similar to that of Newtonian solutions. 2. Description of the model Our 1PN black-hole–disk models are constructed assuming the metric of the form U(x,y,z) (U(x,y,z))2 A (x,y,z) ds2 = g dxµdxν = 1 2 2 (dx0)2 2 i dxidx0 µν − − c2 − c4 − c3 ! U (x,y,z) + 1 2 dx2+dy2+dz2 , (1) − c2 (cid:18) (cid:19) (cid:0) (cid:1) where we use Cartesian coordinates x0 = ct,x1 = x,x2 = y,x3 = z, and c is the speed of light. We write the energy-momentum tensor as M c2uα uβ Tαβ = c BH BHδ(x)+ρ(c2+h)uαuβ +pgαβ, (2) √g u0 BH where the first component describes the point particle (it is proportional to the Dirac delta distribution and models the central black hole) at rest, located at the origin of the coordinate system; the second one is the energy-momentum tensor of the disk matter. Here M denotes the c mass of the point particle; g is the determinant of the metric g = det(g ). The four-vectors µν − uα and uα denote the four-velocities of the central point-mass and the fluid, respectively. The BH symbolρdenotes thebaryonic rest-mass density, his the specificenthalpy, and p is thepressure. In the following sections we will also work with the three-velocity, defined as vi = cui/u0, i= 1,2,3. In the remaining part of the article we use standard cylindrical coordinates (r,z,φ). We consider a stationary, selfgravitating, axially and equatorially symmetric polytropic disk, rotatingaroundacentralpointmassM accordingtotheKeplerianrotationallaw vφ = ω r−3/2. c 0 0 We assume that the disk is geometrically bounded by the inner and outer radius r and r , in out respectively. We introduce the notation according to which any quantity ξ (if it is necessary) is separated into its Newtonian ξ and post-Newtonian ξ part according to the general pattern 0 1 ξ = ξ +ξ /c2. Following [7] we derive basic equations which, separated into their Newtonian 0 1 and post-Newtonian parts, read: ∆U = 4πG(ρ +M δ(x)), (3) 0 0 c h = U + dr(vφ)2r+C , (4) 0 − 0 0 0 Z ∂ A ∆A = 2 r φ 16πGr2ρ vφ, (5) φ r − 0 0 ∆U = 4πG M UD(0)δ(x)+ρ +2p +ρ h 2U +2r2(vφ)2 , (6) 1 c 0 1 0 0 0− 0 0 (cid:16) (cid:16) 3 (cid:17)(cid:17) h = U vφA +2h (vφ)2r2 dr(vφ)4r3 h2 4h U 2U2 C , (7) 1 − 1− 0 φ 0 0 − 0 − 2 0− 0 0− 0 − 1 Z where ∆ denotes the flat Laplacian with respect to coordinates (x1,x2,x3), C and C are 0 1 integration constants, and UD is the gravitational potential due to the disk only, i.e., U = 0 0 GM /x +UD. They follow from the conservation law, Tαβ = 0, the continuity equation − c | | 0 ∇α (ρuα)= 0, and Einstein equations α ∇ R G R g = 8π T , (8) µν − 2 µν c4 µν where R is the Ricci tensor and R denotes the Ricci scalar. µν The above system of equations is closed by assuming an equation of state. Our numerical solutions are obtained for a polytropic equation of state of the form p = Kργ, or equivalently h = Kγ/(γ 1)ργ−1 (for the 0th order solution) and h = (γ 1)h ρ /ρ (for the 1PN part), 0 − 0 1 − 0 1 0 where K and γ > 1 are constants. Equations (3–7) should be solved with respect to the Newtonian gravitational potential U (r,z), the post-Newtonian gravitational potential U (r,z), the rotational potential A (r,z) 0 1 φ and the Newtonian and post-Newtonian enthalpy h (r,z) and h (r,z). Any other quantity (the 0 1 density ρ(r,z), the pressure p(r,z), etc.) can be obtained from these five functions. Numerical solutions are obtained as follows. We use the classic Self-Consistent Field (SCF) scheme (cf. [8]) to solve the set of Eqs. (3) and Eq. (4). Given the Newtonian potential Y and 0 the enthalpy h , we can solve Eq. (5) at once. Finally, we use again the SCF scheme to solve 0 the set of Eqs. (6) and (7). For the Keplerian rotation law the above method converges for all values of the parameters that we have tested. 3. Results The main result of our post-Newtonian scheme is the correction to the Newtonian angular velocity profile. The well-known theorem by Poincar´e and Wavre states that Newtonian stationary barotropic disks (or stars) rotate with the angular velocity that can depend on the distance r from the rotation axis only [9]. It turns out that already the 1PN correction can be significant and in general it depends also on z. It can be expressed as A φ φ φ φ v (r,z) = ∂ v +2rh ∂ v . (9) 1 −2rvφ r 0 0 r 0 0 Note that the above formula involves both geometric and hydrodynamical factors. Acceptable 1PN solutions should satisfy the following (quite stringent) conditions: i) 1 U /c2 U /c4, ii) 2GM /c2 r , iii) c c , where c is a speed of sound. Sample 0 1 c in s s ≫ | | ≫ | | ≪ ≫ numerical models that do satisfy the above conditions are presented in [1]. It turns out that the 1PN correction to the velocity can be of the order of 10% of the Newtonian value. The reader interested in the details of these model may consult [1]. In the remainder of this paper we focus on the construction of virial tests that can be applied to our numerical scheme. 4. Virial identities In this section we use Cartesian and cylindrical coordinates. It is implicitly assumed that Latin indices refer to Cartesian coordinates. The virial identity that can be used to test the Newtonian solution (including the central point-mass) was obtained in [6]. It reads E +2E +2E = 0, (10) pot kin therm where E = 1 d3x ρ(U GM /x) is the total potential energy, E = 1 d3x ρv vi is pot 2 R3 0− c | | kin 2 R3 0i 0 the bulk kinetic energy and E = 3 d3x p is the internal thermal energy. We assume, as R therm 2 R3 R a virial test parameter, the value R E +2E +2E pot kin therm ǫ = . E (cid:12) pot (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) In order to obtain the post-Newt(cid:12)onian relations we rewr(cid:12)ite Eqs. (5) and (6) in a slightly different form. Equation Eqs. (5) can be written in Cartesian coordinates as ∆Ai = 16πGρ vi. − 0 0 Equation (6) is split in two parts: ∆U′ = 4πG ρ +2p +ρ h 2U +2r2(vφ)2 , (11) 1 1 0 0 0− 0 0 ∆U′′ = 4πGM(cid:16) UD(0)δ(x), (cid:16) (cid:17)(cid:17) (12) 1 c 0 ′ ′′ ′′ where U = U +U . The solution for U is 1 1 1 1 GM UD(0) U′′ = c 0 . 1 − x | | Consider a vector 1 1 ai = xl∂ A + A ∂iAk xi∂ Ak∂lA . l k k l k 2 − 2 (cid:18) (cid:19) Its divergence reads 1 1 ∂ ai = xl∂ A + A ∆Ak = 16πG xl∂ A + A ρ vk. (13) i l k 2 k − l k 2 k 0 0 (cid:18) (cid:19) (cid:18) (cid:19) For a finite disk (ρ of compact support), A tends to zero sufficiently fast, and 0 k x2ai 0, as x = x2+y2+z2 . | | → | | → ∞ Thus, by integrating Eq. (13) over R3, and makinpg use of the Gauss theorem, we see that 1 0 = d3x xl∂ A + A ρ vk. R3 l k 2 k 0 0 Z (cid:18) (cid:19) Integrating by parts, one can get rid of the term with ∂ A . This yields l k 5 0 = d3x ρ A vk xl∂ ρ vk A , R3 −2 0 k 0 − l 0 0 k Z (cid:18) (cid:16) (cid:17) (cid:19) where only 0th order terms are differentiated. The above relation can be also written in terms of the vector components in cylindrical coordinates. Because of symmetry assumptions, we have A vk = A vφ, and k 0 φ 0 xl∂ ρ vk A = xl ρ vk A = ρ A vφ +xl∂ ρ vφ A . l 0 0 k ∇l 0 0 k 0 φ 0 l 0 0 φ (cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) This yields 7 0 = d3x ρ A vφ xl∂ ρ vφ A . R3 −2 0 φ 0 − l 0 0 φ Z (cid:18) (cid:16) (cid:17) (cid:19) The virial relation following from Eq. (11) can be obtained in a similar way. It is enough to consider the divergence of the vector 1 1 bi = xl∂ U′ + U′ ∂iU′ xi∂ U′∂lU′. l 1 2 1 1− 2 l 1 1 (cid:18) (cid:19) It yields 1 ∂ bi = xl∂ U′ + U′ ∆U′, i l 1 2 1 1 (cid:18) (cid:19) and analogously 1 0= d3x xl∂ U′ + U′ ρ +2p +ρ h 2U +2r2(vφ)2 . R3 l 1 2 1 1 0 0 0− 0 0 Z (cid:18) (cid:19)(cid:16) (cid:16) (cid:17)(cid:17) Many different forms of the above relation can be obtained by ‘playing’ with Eq. (7). A helpful relation that can be used here is xl∂ U′′ = U′′. − l 1 1 ′ ′ ′ ′ In the analogy to the Newtonian case we choose as virial test parameters ǫ = (ǫ +ǫ )/ǫ ′′ ′′ ′′ ′′ | a b a| and ǫ = (ǫ +ǫ )/ǫ , where | a b a| 7 ǫ′ = d3x ρ A vφ, (14) a R3 2 0 φ 0 Z ǫ′ = d3x xl∂ ρ vφ A , (15) b − R3 l 0 0 φ Z (cid:16) (cid:17) 1 ǫ′′ = d3x U′ ρ +2p +ρ h 2U +2r2(vφ)2 , (16) a − R3 2 1 1 0 0 0− 0 0 Z (cid:16) (cid:16) (cid:17)(cid:17) ǫ′′ = d3x xl∂ U′ ρ +2p +ρ h 2U +2r2(vφ)2 . (17) b R3 l 1 1 0 0 0− 0 0 Z (cid:16) (cid:16) (cid:17)(cid:17) In Table 1 we report results of the convergence tests for a sample system. In our implementation, numerical precision is controlled by the resolution of the grid, the maximum number L of Legendre polynomials used in the angular expansion of the solutions of the scalar andvectorPoissonequations,andavalueofthemaximaldifferencebetweendensitydistributions obtained in the last two consecutive iterations ρ . (In each iteration we compute the quantity tol (k+1) (k) ρ = max ρ ρ . Here index k numbers subsequent iterations; indices i and j refer err i,j| i,j − i,j| to different grid nodes. The iteration procedure is stopped, when ρ ρ .) err tol ≤ The results of the virial test do depend on the grid resolution and almost do not depend on the numbers of Legendre polynomials and the precision in convergence procedure. Table1. Typicaldependenceoftheresultsontheresolutionofthenumericalgrid,themaximum number of the Legendre polynomials L, and the tolerance coefficient ρ . These results are tol obtained for a polytropic fluid with polytropic index γ = 5/3 and the Keplerian rotation law vφ = ω /r3/2. The mass of the disk is M = 0.327M , the inner and outer radii are r = 50R 0 0 d c in S and r = 500R respectively. out S ′ ′′ Resolution L ρ ǫ ǫ ǫ tol 100 100 100 10−6 2.53 10−5 2.29 10−4 3.51 10−4 200 × 200 100 10−6 6.36×10−6 5.87×10−5 8.87×10−5 400 × 400 100 10−6 1.56×10−6 1.93×10−5 2.29×10−5 800 × 800 100 10−6 3.66×10−7 3.65×10−6 4.62×10−6 1200× 1200 100 10−6 1.44×10−7 1.56×10−6 2.49×10−6 1600 × 1600 100 10−6 6.66×10−8 5.43×10−7 1.39×10−6 × × × × 400 400 50 10−6 1.56 10−6 1.93 10−5 2.27 10−5 400 × 400 75 10−6 1.56×10−6 1.93×10−5 2.29×10−5 400 × 400 100 10−6 1.56×10−6 1.93×10−5 2.29×10−5 400 × 400 125 10−6 1.56×10−6 1.93×10−5 2.30×10−5 400 × 400 150 10−6 1.56×10−6 1.93×10−5 2.30×10−5 × × × × 400 400 100 10−5 1.52 10−6 1.94 10−5 2.29 10−5 400 × 400 100 10−6 1.56×10−6 1.93×10−5 2.29×10−5 400 × 400 100 10−7 1.59×10−6 1.93×10−5 2.29×10−5 400 × 400 100 10−8 1.59×10−6 1.93×10−5 2.29×10−5 × × × × Acknowledgments Thisresearchwascarriedoutwiththesupercomputer‘Deszno’purchasedthankstothefinancial supportof the European Regional Development Fund in the framework of the Polish Innovation Economy Operational Program (contract no. POIG. 02.01.00-12-023/08). The work of PJ was partially supported by the Polish NCN grant Networking and R&D for the Einstein Telescope. PM and MP acknowledge the support of the Polish Ministry of Science and Higher Education grant IP2012 000172 (Iuventus Plus). References [1] P. Jaranowski, P.Mach, E. Malec, and M. Pir´og, arXiv:1410.8527v1 [gr-qc] 30 Oct 2014 [2] P. Mach, E. Malec, and M. Pir´og, ActaPhys. Pol. B44, 107 (2013) [3] P. Mach, E. Malec, and M. Pir´og, ActaPhys. Pol. B43, 2141 (2012) [4] S. Nishida and Y.Eriguchi, ApJ. 427, 429 (1994) [5] S. Nishida, Y.Eriguchi, and A. Lanza, ApJ. 401, 618 (1992) [6] P. Mach, Mon. Not.R. Astron. Soc. 422, 772 (2012) [7] T. Damour, P.Jaranowski, and G. Sch¨afer, Phys. Lett.B513, 147 (2001) [8] J.P. Ostrikerand J.W-K. Mark, ApJ.151, 1075 (1968) [9] J.L. Tassoul, Theory of Rotating Stars, Princeton Univ.Press, Princeton, NJ1978