Numerical relativity simulations of thick accretion disks around tilted Kerr black holes

# Numerical relativity simulations of thick accretion disks around tilted Kerr black holes

Vassilios Mewes Departamento de Astronomía y Astrofísica, Universitat de València, Doctor Moliner 50, 46100, Burjassot (València), Spain    José A. Font Departamento de Astronomía y Astrofísica, Universitat de València, Doctor Moliner 50, 46100, Burjassot (València), Spain Observatori Astronòmic, Universitat de València, C/ Catedrático José Beltrán 2, 46980, Paterna (València), Spain    Filippo Galeazzi Institut für Theoretische Physik, Max-von-Laue-Straße 1, 60438 Frankfurt, Germany    Pedro J. Montero Max-Planck-Institute für Astrophysik, Karl-Schwarzschild-Str. 1, 85748, Garching bei München, Germany    Nikolaos Stergioulas Department of Physics, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece
###### Abstract

In this work we present 3D numerical relativity simulations of thick accretion disks around tilted Kerr black holes. We investigate the evolution of three different initial disk models with a range of initial black hole spin magnitudes and tilt angles. For all the disk-to-black hole mass ratios considered () we observe significant black hole precession and nutation during the evolution. This indicates that for such mass ratios, neglecting the self-gravity of the disks by evolving them in a fixed background black hole spacetime is not justified. We find that the two more massive models are unstable against the Papaloizou-Pringle (PP) instability and that those PP-unstable models remain unstable for all initial spins and tilt angles considered, showing that the development of the instability is a very robust feature of such PP-unstable disks. Our lightest model, which is the most astrophysically favorable outcome of mergers of binary compact objects, is stable. The tilt between the black hole spin and the disk is strongly modulated during the growth of the PP instability, causing a partial global realignment of black hole spin and disk angular momentum in the most massive model with constant specific angular momentum . For the model with non-constant -profile we observe a long-lived non-axisymmetric structure which shows strong oscillations of the tilt angle in the inner regions of the disk. This effect might be connected to the development of Kozai-Lidov oscillations. Our simulations also confirm earlier findings that the development of the PP instability causes the long-term emission of large amplitude gravitational waves, predominantly for the multipole mode. The imprint of the BH precession on the gravitational waves from tilted BH-torus systems remains an interesting open issue that would require significantly longer simulations than those presented in this work.

###### pacs:
04.25.D-, 04.30.Db, 95.30.Lz, 95.30.Sf, 97.60.Lf 97.10.Gz

## I Introduction

Stellar mass black hole–torus systems are believed to be the end states of binary neutron star (NS-NS) or black hole neutron star (BH-NS) mergers, as well as of the rotational gravitational collapse of massive stars. BH-torus systems emit gravitational waves (GWs), which may eventually provide the direct means to study their actual formation and evolution and prove the current hypotheses that associate them to gamma-ray burst (GRB) engines Woosley (1993); Janka et al. (1999). Such a possibility is out of reach to electromagnetic observations due to their intrinsic high density and temperature. For an overview of the event rate estimates of NS-NS and BH-NS mergers that are observable with initial and advanced LIGO see e.g. Abadie et al. (2010); Dominik et al. (2013, 2014).

Our theoretical understanding of the formation of BH-torus systems and their evolution relies strongly on numerical work. In recent years a significant number of numerical relativity simulations have shown the feasibility of the formation of such systems from generic initial data (see e.g. Rezzolla et al. (2010); Kyutoku et al. (2011); Hotokezaka et al. (2013a, b); Kastaun and Galeazzi (2015) for recent progress). In particular, the 3D simulations of Rezzolla et al. (2010) (see also references therein) have shown that unequal-mass NS-NS mergers lead to the self-consistent formation of massive accretion tori (or thick disks) around spinning black holes, thus meeting the necessary requirements of the GRB’s central engine hypothesis. However, if the energy released in a (short) GRB comes from the accretion torus, the BH-torus system has to survive for up to a few seconds Rees and Meszaros (1994). Any instability which might disrupt the system on shorter timescales, such as the runaway instability Abramowicz et al. (1983) or the Papaloizou-Pringle instability (PPI hereafter) Papaloizou and Pringle (1984), could pose a severe problem for the prevailing GRB models.

Using perturbation theory, Papaloizou and Pringle found that tori with constant specific angular momentum are unstable to non-axisymmetric global modes. Such modes have a co-rotation radius within the torus, located in a narrow region where waves cannot propagate. Waves can still tunnel through the co-rotation zone and interact with waves in the other region. The transmitted modes are amplified due to the feedback mechanism provided by the reflecting boundaries at the inner and outer edges of the torus. The manifestation of the PPI, as Hawley (1991) first showed numerically through hydrodynamical simulations in the fixed background metric of a Schwarzschild BH, is in the form of counter-rotating epicyclic vortices, or “planets’”, with planets emerging from the growth of a mode of order . Moreover, BH-torus systems are characterized by the presence of a cusp-like inner edge where mass transfer driven by the radial pressure gradient is possible. If due to accretion the cusp moves deeper inside the disk material, the mass transfer speeds up leading to a runaway process that destroys the disk on a dynamical timescale (see e.g. Font and Daigne (2002) for test-fluid simulations in general relativity of the occurrence of this instability and Montero et al. (2010) for axisymmetric simulations where self-gravity was first taken into account). In most recent numerical relativity simulations Rezzolla et al. (2010); Hotokezaka et al. (2013a); Neilsen et al. (2014); Kastaun and Galeazzi (2015) the BH-torus systems under consideration did not manifest signs of dynamical instabilities on short dynamical timescales, because the non-constant angular momentum profiles of the massive disks that formed in the simulations seem to make them stable against the development of the runaway instability. However, on longer timescales non-axisymmetric instabilities can easily develop Korobkin et al. (2011); Kiuchi et al. (2011).

Korobkin et al Korobkin et al. (2011) have studied non-axisymmetric instabilities in self-gravitating disks around BHs using three-dimensional hydrodynamical simulations in full general relativity. Their models incorporate both moderately slender and slender disks with disk-to-BH mass ratios ranging from 0.11 to 0.24. While no sign of the runaway instability was observed for these particular models, unstable non-axisymmetric modes were indeed present. More recently Korobkin et al. (2013) observed that by a suitable choice of model parameters, namely constant angular momentum tori exactly filling or slightly overflowing their Roche lobe, a rapid mass accretion episode with the characteristics of a runaway instability sets in. The astrophysical significance of such fine-tuned models is uncertain as, e.g. they do not seem to be favored as the end-product of self-consistent numerical relativity simulations of binary neutron star mergers Rezzolla et al. (2010); Hotokezaka et al. (2013a); Kastaun and Galeazzi (2015); Sekiguchi et al. (2015). Correspondingly, the 3D general relativistic numerical simulations of BH-torus systems of Kiuchi et al. (2011) showed that an non-axisymmetric instability grows for a wide range of self-gravitating tori orbiting BHs. Their models included torus-to-BH mass ratios in the range and both constant and non-constant angular momentum profiles in the disks. Kiuchi et al. (2011) found that the non-axisymmetric structure persists for a timescale much longer than the dynamical one, becoming a strong emitter of large amplitude, quasi-periodic GWs, observable by forthcoming detectors.

The vast majority of numerical relativity simulations to date that have studied the formation of BH-torus systems have led to systems in which the accretion torus is aligned with the equatorial plane of the central Kerr BH. There are, however, several ways of arriving at tilted accretion disks around Kerr BH. Possible scenarios include asymmetric supernova explosions in binary systems Fragos et al. (2010) or via BH-NS mergers in which the BH spin is misaligned with the orbital plane of the binary. See also Fragile et al. (2001); Maccarone (2002) for arguments that most BH-torus systems should be tilted. Recently, full GRHD simulations of BH-NS mergers with tilted initial BH spins have been performed by Foucart et al. (2011, 2013); Kawaguchi et al. (2015), which can result in tilted accretion tori around the BH remnant. In these simulations massive remnant disks only formed for high initial BH spins. As shown in Foucart (2012), the disk mass in BH-NS mergers increases with increasing initial BH spin and decreases with increasing initial BH mass. This is due to the size of the ISCO (innermost circular stable orbit) of the BH in the merger. The ISCO grows with BH mass and decreases with the spin magnitude of the BH. If the ISCO is large enough, the NS will be ’swallowed’ entirely by the BH before being tidally disrupted, leaving no accretion torus behind. In order to estimate disk masses of resulting from BH-NS mergers, one needs thus an estimate for the initial BH masses in these system. One such method, via population synthesis considerations, favor larger BH masses Belczynski et al. (2008, 2010); Fryer et al. (2012), with a peak around and a mass gap between the lightest expected BHs and NS masses. This means that these massive BH would need very large initial spins in order to be able form massive remnant disks after the BH-NS merger.

General relativity hydrodynamics and MHD simulations of such kind of tilted thick accretion disks around BHs have been performed by Fragile and collaborators within the test-fluid approximation (see e.g. Fragile and Anninos (2005); Fragile et al. (2007)). These authors have carried out an extensive program to study both the dynamics and observational signatures of thick accretion tori around tilted Kerr black holes.

Motivated by these previous works we present in this paper an extended set of numerical relativity simulations of tilted accretion disks around Kerr BHs in which the self-gravity of the disks is taken into account. In this first study of the effects of the disk self-gravity in tilted BH-torus systems, we have chosen low mass central black holes, and small spins that are unlikely to produce massive tilted disks in astrophysical mergers for the reasons outlined above. As a first study, we want to investigate the effects of the disk self-gravity around tilted BH with our existing initial data, where some of our initial models were known to be unstable against the PPI in the untilted case Korobkin et al. (2011). We plan on studying these tilted BH-torus systems with more astrophysically realistic data in future works. Our simulations incorporate the effects of the self-consistent evolution of BH-torus systems with misaligned spins via the solution of the full 3+1 Einstein equations coupled to the hydrodynamics equations. Due to the Lense-Thirring effect, the inner regions of the disks start precessing and might become twisted and warped affecting their dynamical behavior (for a definition of twist and warp see e.g. Nelson and Papaloizou (2000)). Those inner regions might also be forced into the equatorial plane of the BH due to the so-called Bardeen-Petterson effect Bardeen and Petterson (1975). Our simulations allow us to investigate these issues and find out how the dynamics of the tilted disks affect the central BH, assessing in turn the importance of neglecting the disk’s self-gravity on the overall dynamics of these systems.

This paper is organized as follows: In Section II we summarize the mathematical and numerical framework we employ in our numerical study, explaining in detail the methods we use for the spacetime and matter evolution. We also describe several diagnostics we use to monitor and analyze the simulations. In Section III we describe the setup of our initial data, both in the untilted and tilted cases. We also provide a table which contains the key parameters of the models considered in this study. Our results are presented and discussed in detail in Section IV. Finally, we present our conclusions in Section V. Throughout the paper, we use units in which where , and are the speed of light, the gravitational constant and the solar mass, respectively. Latin indices run from 1 to 3 while Greek indices run from 0 to 3.

## Ii Mathematical and numerical framework

### ii.1 Spacetime evolution

The simulations reported in this paper are performed using the publicly available Einstein Toolkit (ET) code ET (); Löffler et al. (2012). The ET is a code for relativistic astrophysics simulations, which uses the modular Cactus framework Cac () (consisting of general modules commonly called ‘thorns’) and provides adaptive mesh refinement (AMR) via the Carpet driver Car (). The set of complex time-dependent partial differential equations is integrated using the method-of-lines where the time integration is performed using a 4th order Runge-Kutta algorithm. The toolkit provides numerical solvers and the necessary infrastructure to evolve the Einstein equations,

 Gμν=8πTμν, (1)

coupling the evolution of spacetime to the dynamic of the matter content through the stress-energy tensor . The methods to evolve the Einstein equations are based on the 3+1 ADM Arnowitt et al. (2008) splitting of spacetime, in which the metric becomes

 ds2=(−α2+βiβi)dt2+2βidtdxi+γijdxidxj, (2)

where is the lapse function, the shift vector, and is the spatial metric tensor. The extrinsic curvature tensor is constructed from the spatial metric as follows:

 Kij≡−12α(∂t−Lβ)γij, (3)

where stands for the Lie derivative with respect to the shift vector.

The left hand side of the Einstein field equations is evolved using the McLachlan code Brown et al. (2009); Reisswig et al. (2011), which solves a conformal-traceless “” formulation of the Einstein equations known as BSSN Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998).

The BSSN evolution equations are evolved using a fourth-order accurate, centered finite-differencing operator, while the advection terms for the shift vector are evolved with a fourth-order upwind stencil. We apply fifth-order Kreiss-Oliger Dissipation to all spacetime variables to achieve overall fourth-order for the spacetime evolution. Using the BSSN evolution system and the gauge conditions described in Appendix A, the McLachlan code can evolve stably a single puncture for several light crossing times, see e.g. Löffler et al. (2012).

We apply Sommerfeld outgoing boundary conditions for all the BSSN evolution variables with an extrapolation to include the part of the boundary conditions that does not behave like a simple outgoing wave Alcubierre et al. (2000).

### ii.2 Matter Evolution

The evolution of the hydrodynamics is performed by GRHydro Baiotti et al. (2005); Löffler et al. (2012); Mösta et al. (2014), that solves the general relativistic hydrodynamics equations in flux-conservative form, in the so-called Valencia formulation Martí et al. (1991); Banyuls et al. (1997); Font (2008), using high-resolution shock-capturing (HRSC) methods. The GRHD evolutions equations are described in Appendix B.

As a grid-based code GRHydro cannot evolve regions without matter content due to the breakdown of the EOS and evolution equations in regions of zero rest-mass density. As is customary, a low density atmosphere (typically several orders of magnitude below the maximum density of the initial data) fills those grid points that should belong to vacuum regions. In all our simulations we set the density of the atmosphere to times the initial maximum rest-mass density in the torus. Although fully evolved, the dynamical effects of the atmosphere regions can be neglected.

We use the piecewise parabolic method (PPM) Colella and Woodward (1984) to reconstruct the primitive variables at the cell interfaces and Marquina’s flux formula Donat and Marquina (1996); Aloy et al. (1999) to compute the numerical fluxes. Differently from the original implementation reported in Ref. Baiotti et al. (2007a), we reconstruct the quantities instead of the three-velocities . This guarantees that the velocities reconstructed at the cell boundaries remain subluminal even under extreme conditions like those encountered near the apparent horizon (AH) of the BH.

### ii.3 Diagnostics

We measure the BH mass during the evolution using the AHFinderDirect thorn which implements the AH finder described in Thornburg (2003). The BH spin magnitude and direction are measured using the QuasiLocalMeasures thorn Dreyer et al. (2003); Schnetter et al. (2006) (see also the review of quasi-local methods in Szabados (2009)). We measure the spin direction using the so-called flat-space rotational Killing vector method of Campanelli et al. (2007), which can be derived using Weinberg’s pseudotensor in Gaussian coordinates and is equal to the Komar angular momentum integral when the latter is expressed in a foliation adapted to the axisymmetry of spacetime Mewes et al. (2015a). To track the moving location of the puncture during the evolution, we use the PunctureTracker thorn. Gravitational waves are extracted using a multipole expansion of the Weyl scalar at different radii. The thorns Multipole and WeylScal4 compute the respective quantities. We calculate the mass accretion rate as the instantaneous flow of matter through the AH, using the Outflow thorn. To quantify the growth of non-axisymmetric modes in the disk, we use the thorn GRHydroAnalysis which computes 3D Fourier integrals of the density, described in Baiotti et al. (2007b). The amplitude of the -th mode is given by:

 Dm=∫α√γρe−imϕd3x. (4)

Since the computational domain is very large (see section III.3 below), the total mass of the atmosphere can be a non-negligible fraction of the total rest-mass of the torus. We therefore have to ignore cells corresponding to the atmosphere when integrating the torus rest-mass density and in the calculations of the non-axisymmetric modes in the disk. Despite the non-negligible mass of the atmosphere, this mass is not dynamically important because it is spatially homogeneous (and therefore its gravitational pull on the central BH-torus system is zero) and its coordinate velocity is set to zero.

### ii.4 Twist and tilt angles

In order to keep track of the response of the tilted accretion disk and to check for Lense-Thirring precession and the occurrence of the Bardeen-Patterson effect, we measure two angles, the twist (precession) and tilt (inclination) of the disk Nelson and Papaloizou (2000). We closely follow the description of the angles given in Fragile and Anninos (2005), with a modification in the calculation of the twist, and furthermore adapt the way they are calculated.

The main idea is to split the disk up into a series of annuli and calculate the angular momentum vector of the matter, , for each individual annulus. Using this vector, we define two angles, the twist

 (5)

and the tilt

 ν(r)=∠(SBH,JDisk(r)), (6)

where

 P(a,b)=a−a⋅b|b|2b, (7)

is the projection of vector onto the plane with normal .

The vector is constructed in the following way: We project the BH spin onto the -plane, and then rotate the resulting vector by about the -axis. The cross product of and then lies in the equatorial plane of the BH and the twist at is throughout the disk. A sketch showing the construction of these vectors and twist and tilt angles is provided in Fig. 1.

The twist is then a measure of how much the angular momentum vector of each annulus has precessed in the (dynamically changing) equatorial plane of the BH, while the tilt gives the angle between the BH spin and the angular momentum vectors for each annulus. The disk is said to be twisted if the twist becomes a function of radius, , and said to be warped if the tilt becomes a function of radius, .

The angular momentum vector for each annulus is calculated following the procedure outlined in the section on spin in Weinberg (1972) (pages 46f.). We calculate from the anti-symmetric angular momentum tensor

 Lμν=∫d3x(xμTν0−xνTμ0), (8)

in the following way:

 Jx=L23,Jy=J31andJz=L12. (9)

The relevant components of the stress-energy tensor are

 Ti0=ρh(W2αvi−Wαβi)+Pβiα2, (10)

where Latin indices run from to .

We note that our procedure differs slightly from that given in Fragile and Anninos (2005), where the components of the disk shell angular momentum vectors were computed using the intrinsic spin as follows:

 Sα=12ϵαβγδJβγUδ, (11)

where is the 4-dimensional Levi-Civita symbol and is the total 4-velocity of the system, with being the 4-momentum.

Using the extrinsic spin, the contribution to the total angular momentum stemming from choosing a reference point is accounted for in the linear momentum of the center of mass of the system. We choose to use the angular momentum tensor , calculated about the center of the BH instead, using formula (8) in the calculation of and . We therefore pick a reference frame centered on the BH, as it is moving around the grid and we are interested in the total angular momentum of the disk shell about the BH center. We note that the calculation of the twist and tilt angles is in fact gauge dependent. We have checked the impact of the gauge choice on the twist and tilt angles by measuring the angles for the initial data with different lapse and shift profiles to find that the angles are calculated correctly (as the initial tilt and twist profiles are parameters of our initial data), which also serves as an important consistency check for the correct measure of the initial angles. In the actual numerical implementation, we do not use the inverse cosine to calculate the angles, as the function becomes very inaccurate when the vectors become close to being parallel or anti-parallel, instead we use the atan2 function which is free of this deficiency. The tilt is then calculated using the following formula:

 ν(r)=tan−1[|^SBH×^JDisk(r)|,^SBH⋅^JDisk(r)] (12)

where the hat indicates unit vectors.

For the twist , we need to be more careful. In the formula above, we implicitly assume that , because Eq. (12) will always return angles . Although it does not generally make much sense to consider angles between two 3D vectors, because there is an up and down direction, we are nevertheless interested in directional (signed) angles for the twist. The reason is that we want to be able to track cumulative twists larger than . In order to calculate the directional angle from the reference vector to the target vector in a fixed sense of rotation, we need a third reference vector that always lies above the plane spanned by the two original vectors. As the plane in which both vectors live is constructed to be the equatorial plane of the BH, we can therefore choose this vector to be . This information allows us (because the cross product contains this information in the sign it returns) to calculate twist angles in the range with the following formula:

 σ(r)=tan−1[|JBH×JDisk(r)|,JBH⋅JDisk(r)] (13)

and then select the directional angle by:

 σ(r)={σ(r):SBH⋅(SBH×JDisk(r))≥02π−σ(r):SBH⋅(SBH×JDisk(r))<0

We also compute the instantaneous precession and nutation rates for both the BH spin vector and the total disk angular momentum vector (which we calculate by summing the components of all shells). For the BH, we compute its precession and nutation about the -axis, while the total precession and nutation of the disk angular momentum vector is calculated about the BH spin axis.

## Iii Initial data and setup

Our initial setup is a thick, self-gravitating axisymmetric accretion disk in equilibrium around a rotating BH. Such a system is built following the approach laid out in Stergioulas (2011) that we briefly describe below. The solver to build the initial data (ID) first computes models of self-gravitating, massive tori around non-rotating BHs. Then, for our simulations of tilted disks around rotating BHs we retain the hydrodynamical content of a model and replace the spacetime by a tilted Kerr spacetime in quasi-isotropic coordinates. The generated initial data is then interpolated onto Cartesian coordinates and evolved with the Einstein Toolkit. We are not aware of a method to generate self-consistent data for self-gravitating tilted disks around Kerr black holes and therefore resort to this method. In order to keep the constraint violations as low as possible, we use disk models with rest masses of only a few percent of the central BH mass and with small to moderate BH spins. When replacing the original spacetime of the initial data with the tilted Kerr spacetime in the computational domain, both the torus rest-mass and the torus gravitational mass will change somewhat, due to the fact that the volume element of the spacetime , where is the determinant of the four-dimensional metric , as well as the lapse, have changed. In order to assess the difference in and due to the replacement of the spacetime we calculate both quantities as in Stergioulas (2011), as follows:

 M0=∫ρut√−gd3x=∫√γWρd3x, (14)

and

 MT =∫(−2Ttt+Tμμ)α√γd3x. (15)

We find the largest difference in and in model C1Ba03b30, where we obtain the following fractional differences: and .

### iii.1 Self-gravitating accretion disks

Our ID are self-gravitating disks around a Schwarzschild BH in quasi-isotropic (QI) coordinates Stergioulas (2011). The ID is fully described by the four metric potentials , , and that characterize the metric of a stationary, axisymmetric spacetime in spherical polar QI coordinates,

 ds2=−e2νdt2+e2α(d¯r2+¯r2dϕ2)+B2e2ν¯r2sin2θ(dϕ−ωdt)2, (16)

where the isotropic radius is given in terms of the areal radius by

 ¯r=12(r−M+√r2−2Mr), (17)

where is the mass of the BH. The ID also provides the pressure and orbital angular velocity , measured by an observer at infinity at rest. The evolution time is measured in terms of the initial orbital period at the radius of the initial maximum of the rest-mass density, . These two quantities, together with the EOS, are sufficient to obtain the remaining hydrodynamical quantities. The ID is constructed using the polytropic EOS described in Eq. (58). The components of the 3-velocity in Cartesian coordinates are given by

 vi=(y(ω−Ωα),−x(ω−Ωα),0) (18)

where is the lapse function. The Lorentz factor is calculated directly from the 3-velocity

 W=1√1−vivi=1√1−(ω−Ω)2B2e4ν(x2+y2). (19)

The matter content of the ID only fills the computational domain up to the event horizon of the central BH. In our simulations we do not excise the BH but rather treat it as a puncture. We noticed that when evolving the BH as a puncture with matter content without excision, the simulations were not long-term stable and failed at early stages of the evolution due to errors in the matter evolution at the location of the puncture. The reason is a very rapid pile up of matter at this location, which leads to non-physical values of the hydrodynamical variables at the puncture. These lead to failures in the GRHydro evolution code during the conversion from conserved to primitive variables. In Faber et al. (2007) a successful method to evolve hydrodynamics in the presence of punctures was presented. This method checks for unphysical values of the hydrodynamical variables at the immediate vicinity of the puncture and resets them to physical values. In practice it checks for the positivity of the conserved energy () and an upper bound for :

 γijSiSj≤τ(τ+2ρ). (20)

In our simulations, we use a similar treatment inside the BH horizon, namely we specify a fraction of the volume of the AH in which we reset the matter fields to atmosphere values and the stress-energy tensor to zero, which is essentially the same technique as described in Reisswig et al. (2013). In this way, we avoid using a moving excision boundary for the hydrodynamical variables as in Löffler et al. (2006). From the point of view of the evolution it is safe to do this, as the evolution of the hydrodynamics inside the horizon cannot influence the evolution outside the horizon. In fact, we checked that the fraction of the AH in which we apply such atmosphere resetting has no effect on the evolution, as expected. We checked this by applying the atmosphere reset for different fractions of the AH. In all cases, the evolution outside the AH was unaffected, confirming that the region inside the AH is causally disconnected. This has been employed for instance in the so-called “turduckening“ of initial BH data Brown et al. (2009). For the practical implementation of this approach we use a spherical surface that contains the shape of the AH and apply the atmosphere method in a fraction of the minimum radius of that surface. Our procedure is illustrated in Fig. 2 where we plot the initial density profile along the -axis for model C1Ba0b0 of our sample (see Table 1 below). The AH is initially located at and there exists a smooth density profile across the horizon. The remaining hydrodynamical quantities are treated in the same way. By using this approach we are able to evolve matter fields in the presence of punctures for very long times without the need for moving hydrodynamical excision zones.

### iii.2 Tilted Kerr spacetime in improved quasi-isotropic coordinates

We set up a Kerr black hole tilted about the -axis by an angle in the improved quasi-isotropic coordinates proposed by Liu et al. (2009). In those coordinates, the radius of the horizon does not shrink to zero in the extreme Kerr limit, but approaches . The initial mass of the Kerr BH is chosen to be equal to the value used in the ID calculation of the self-gravitating torus. The spin parameter varies for different runs. The list of 21 models we use in our investigation is summarized in Table 1.

We perform the rotation of the BH about the -axis in the following way: we first rotate by an angle the coordinates to the tilted coordinates,

 ⎛⎜⎝^x^y^z⎞⎟⎠=⎛⎜⎝1000cosβ0−sinβ00sinβ0cosβ0⎞⎟⎠⎛⎜⎝xyz⎞⎟⎠. (21)

We then calculate the spatial metric and extrinsic curvature by performing coordinate transformations of the respective expressions given in Liu et al. (2009) (in spherical polar coordinates) to our tilted Cartesian coordinates . We set the initial shift to and choose the following initial lapse profile:

 α=21+(1+M2r)4. (22)

We find that this initial lapse profile greatly helps with the conservation of irreducible mass and spin during the initial gauge transition the system undergoes. The bigger the difference between the initially chosen lapse profile and the profile the system attains after the transition to the puncture coordinates, the larger the errors in irreducible mass and spin. As all quantities have been set up in the tilted coordinates , we also need to apply rotations to vectors and tensors, in order to obtain the correct components of the shift, spatial metric and extrinsic curvature in the original coordinate system that is going to be used in the simulations. We obtain the components of the shift vector with

 β=R^β, (23)

and the components of the spatial metric and extrinsic curvature with

 γ = R^γRT, (24) K = R^KRT, (25)

where is the rotation matrix rotating from to . Our choice of the initial BH spin magnitude is currently limited by our computational method and is therefore somewhat smaller than spins expected in certain astrophysical scenarios. For example, as the simulations of misaligned BH-NS mergers of Foucart et al. (2011, 2013); Kawaguchi et al. (2015) have shown, accretion disks form in these systems only for relatively high () initial BH spins; see also Foucart (2012) for disk mass predictions from BH-NS mergers. These initial spins in astrophysical merger scenarios are larger than the spins investigated in our study . The reasons for our choice of smaller spins are the very demanding resolution requirements when evolving highly spinning BHs, see, e.g. Lousto et al. (2012). Furthermore, spin conservation seems to be affected non-linearly with higher spin magnitude Marronetti et al. (2008), which we have also observed in Mewes et al. (2015b). The choice of astrophysically more realistic BH spin magnitudes would have been prohibitively expensive (see, e.g. Lousto et al. (2012) for details) and is beyond the scope of our study of self-gravitating accretion disks around tilted Kerr BHs. Nevertheless, our results with dimensionless spin values of 0.5 should provide a first qualitative understanding of what may happen at dimensionless spin values of 0.7. We also note that the existence of the ’mass gap’ in BH masses in compact mergers has been questioned in Kreidberg et al. (2012). While the peak of BH masses in these new considerations is still around , very small BH masses in mergers are at least not ruled out as in previous studies. Furthermore, the BH mass predictions population synthesis models are based on assumptions about the exact details of the supernova explosion mechanism, see for instance Belczynski et al. (2012), where no mass gap is found using a delayed supernova explosion model assuming a 100-150 ms instability growth time. As we explained in the introduction, the smaller the BH mass, the smaller the BH spin needs to be in order to tidally disrupt the NS during the merger in order to leave a thick accretion torus around the central BH. It therefore seems that some of the combination of BH spins and mass-ratios studied in this work (while having been primarily motivated on computational grounds) are a possible (although maybe not the most likely) outcome from BH-NS mergers. Based on these considerations, model D2, the lightest in our study, is the most astrophysically relevant model of the three models studied. We choose the initial tilt angle . In the simulations of tilted BH-NS mergers of Foucart et al. (2011), the authors found that for initial inclinations , the mass of the resulting accretion disk would be around of the initial neutron star mass, whereas higher initial inclinations produced a sharp drop in the mass of the resulting disk. For those those pre-merger inclinations , the resulting tilt of the post-merger disk was around . In the recent simulations of tilted BH-NS mergers of Kawaguchi et al. (2015), the authors report a post-merger disk with a tilt , for an initial tilt of . These findings seem to be in accordance with the probability distributions of post-merger tilt angles for BH-NS mergers obtained in Stone et al. (2013), which show a sharp cut-off in the probability distribution for post-merger tilt angles .

### iii.3 Numerical setup

The use of mesh-refinement techniques is of fundamental importance in our simulations where different physical scales need to be accurately resolved. For this reason, we use the Carpet driver that implements a vertex-centered AMR scheme adopting nested grids Schnetter et al. (2004).

Figure 3 shows isocontours of the initial rest-mass density profile for a representative disk model (C1Ba0b0) in a vertical cut (-plane) together with the initial innermost AMR grid structure. In our simulations, we use 13 levels of mesh refinement, with the outer boundary placed at . We perform a large series of simulations using two different resolutions, high resolution runs with a highest resolution of and simulations with half this resolution, see Table 1. Every refinement level except the innermost one is half the size of the preceding level. The innermost refinement level, as can be seen in Fig. 3, is more than half the size of the next refinement level. All simulations have been performed using a CFL factor of . In the outermost refinement levels, we use the same time step as in level , in order to prevent instabilities in the spatially varying damping parameter . Using the same (smaller) time step in the outer levels is necessary, even though the spatially varying parameter of Alic et al. (2010) takes into account the time step limitation of the parameter described in Schnetter (2010). Details of the accuracy and convergence properties of our simulations are provided in Appendix C.

## Iv Results

### iv.1 Surface plots

We start describing the morphology and dynamics of a representative model of our sample. The evolution of model C1Ba01b30 is shown in Fig. 4. This figure displays surface plots of the logarithm of the rest-mass density at six different stages of the evolution. The final panel shows the morphology of the system after the disk has completed 20 orbits around the central BH. The interpolation of the ID from QI spherical polar coordinates to Cartesian coordinates, along with the inclusion of a moderate BH spin in the originally non-rotating BH, results in a significant perturbation of the equilibrium model that triggers a phase of oscillations of the torus around its equilibrium. These oscillations are present throughout the simulation and, as the disk is initially filling its Roche lobe, they induce a small accretion process of matter through the cusp towards the BH. This does not reduce significantly the total rest-mass of the torus as we show below.

The oscillations of the disk are damped due to numerical viscosity and also by the formation of shock waves that propagate inside the disk and convert kinetic energy into thermal energy. These outward-propagating shocks, also found in the simulations of Korobkin et al. (2011), are produced by in-falling material along the funnel walls which reaches inner high-density regions and bounces back. One such shock wave can be clearly seen on the right portion of the disk in the 12 orbits panel of Fig. 4. The initial tilt of the BH spin () twists and warps the innermost parts of the disk closer to the BH. This gradual effect becomes evident on the time series shown in Fig. 4 and by the final 20 orbits the disk is significantly misaligned with the (horizontal) -axis.

In Fig. 5 we show surface plots of the rest-mass density for the models C1B, NC1 and D2 for the lower spin simulations () at the final time of the evolution . The effects of the initial BH tilt on the morphology of the disk at the final time are more clearly pronounced for the heavier models C1B (top row) and NC1 (middle row). Although the spin of the BH is very low and the spacetime is therefore nearly spherically symmetric, the initial tilt angle nevertheless causes a different evolution. In contrast, the significantly lighter model D2 is seen to be hardly affected by the initial tilt by the end of the evolution, at least for these simulations with low initial spin. Similarly, in Fig. 6, we plot the same rest-mass density surface plots as in Fig. 5 for all our models with higher initial spin. Comparing the evolution of the disk morphology with the lower initial spin figure, we see that there are now significant changes in the morphology when increasing the initial tilt angle for all models. The effects of evolving the disk in the tilted Kerr spacetime are now more pronounced, as the higher spin causes a significant deviation from spherical symmetry in the spacetime. The growth of non-axisymmetric modes associated with the PPI (see below) in the two more massive models C1B and NC1 causes an alignment of the overall disk angular momentum vector with the BH spin. This interesting effect will be discussed in Section IV.6 below, where we comment on the BH precession and nutation properties.

In Figures 7 and 8 we plot linearly spaced isocontours and the logarithm of the rest-mass density in the plane perpendicular to the total angular momentum vector of the disk. This kind of “equatorial” plots allows for a better visualization of the possible growth of non-axisymmetric structures in the disk. We note that due to the twisting and warping of the disk as a result of the initial BH tilt, the choice of the plane on to which project the isocontours for an adequate visualization is not straightforward and careless choices (as e.g. the simple choice of the equatorial -plane) may hide important pieces of information on the dynamics and morphology. The top row of Fig. 7 corresponds to models C1Ba01b5 (left), C1Ba01b15 (middle), and C1Ba01b30 (right), and all of them are displayed at the time at which the growth of the fastest growing nonaxisymmetric structure in the disk saturates (see Section IV.4 below where we discuss the mode growth of the PPI for the different models of our sample). The middle row of Fig. 7 corresponds to models NC1a01b5 (left), NC1a01b15 (middle), and NC1a01b30 (right), also shown at the time the fastest growing mode has maximum amplitude, while the bottom row shows models D2a01b5 (left), D2a01b15 (middle), and D2a01b30 (right), which are displayed at the final time of the evolution () when the amplitudes of the corresponding PPI modes (if present) are largest. Correspondingly Fig. 8 shows the same type of isocontour plots as Fig. 7 but for the models with higher initial spin. The top row of Fig. 8 corresponds to models C1Ba03b5 (left), C1Ba03b15 (middle), and C1Ba03b30 (right), the middle row of Fig. 7 corresponds to models NC1a03b5 (left), NC1a03b15 (middle), and NC1a03b30 (right) and the bottom row shows models D2a05b5 (left), D2a05b15 (middle), and D2a05b30 (right). As in Fig. 7, the times shown are at either saturation of the fastest growing non-axisymmetric structure or when the maximum mode amplitudes are largest.

The inspection of the different panels displayed in Figs. 7 and 8 shows noticeable differences in the flow morphology. On the one hand all six models C1B (top rows in both figures), corresponding to initial BH spins of and , respectively, exhibit a very prominent spiral arm once the stationary accretion phase is reached. This structure, which is visible for all three tilt angles considered, and , is associated with the growth of the non-axisymmetric PPI mode, which is the fastest-growing mode and has the largest amplitude (see also Fig. 11). The middle rows, corresponding to the six models NC1, show the development of a spiral arm as well, however not as clearly pronounced as in the case of C1B. On the other hand the isodensity contours of all three models D2 displayed in the bottom row of Fig. 7 for a BH spin barely show any morphological deviation from the initial axisymmetry. Only the innermost regions closest to the BH show some slight non-axisymmetric structure, particularly for the model with the largest tilt angle () shown on the right panel. Model D2 is considerable lighter than model C1B, and model NC1 has an initial non-constant specific angular momentum profile and according to Kiuchi et al. (2011) more massive tori with const profiles favor the appearance of the PPI with respect to less massive tori (see in particular Figs. 2 and 3(a) of Kiuchi et al. (2011) for their model C1 which has very similar properties regarding size and BH-to-disk mass ratio to our model C1B). Our results seem to only partially confirm those previous findings as there seems to be some weak dependence of a moderate PPI growth with increasing tilt angles.

The dependence on the BH spin seems however to be more significant for the light model D2, as shown in the panels at the bottom row of Figures 7 and 8. For the higher spin runs shown in Fig. 8, all three panels show non-axisymmetric morphological features. All three snapshots correspond to the same final time of the evolution () and the dynamical differences between the three models are only due to the initial BH tilt angle. For (left panel) the spiral structure is visible and dominant but, at the same time, the mode seems to have reached an almost similar amplitude (see the bottom panel of Fig. 11). Such feature becomes clearly dominant the larger the BH inclination becomes, as depicted in the central regions of the middle and right panels.

### iv.2 Maximum rest-mass density evolution

We next show the evolution of the maximum rest-mass density in the disks in Fig. 9, normalized to the corresponding initial central density of each model (see Table 1). For model D2 (bottom panel), the effect of tilting the BH has very little effect on the evolution of for a BH spin of . However, for , the evolution of differs for the different initial tilt angles. Model D2a05b5 exhibits an upward drift from orbits onwards to bring to the level of the lower spin runs towards the end of the simulation. This might be connected to the growth of the m=1 non-axisymmetric mode in this model, as discussed in the previous section. In the initial stages of the evolution, is higher for larger tilt angles, a dependence that is also observed in the accretion rate (see Section IV.3). The maximum rest-mass density stays below its initial value for almost the entire evolution for both models D2 and NC1 (middle panel). This is very different in model C1B (top panel) where the occurrence and saturation of the PPI has a drastic effect on the evolution of . When the PPI starts growing, a situation which is accompanied by an exponential growth of the accretion rate, grows to up to twice its initial value, reaching the highest value for the untilted model C1Ba00. It subsequently settles to a value that is very similar for all initial spins and tilt angles, but the time of growth and the shape of the “lump” in the density evolution are different for the different spins and initial tilt angles. Model NC1 shows features from both models D2 and C1B described above. Similarly to model D2, stays below its initial value for almost the entire evolution, with model NC1a00 being an exception where it grows above the initial value briefly. There is, as in models C1B, a rapid growth in the central density, albeit somewhat smaller for the non-constant angular momentum disks, associated with the growth of the PPI in all the models. The magnitudes at which the mode growth saturates are similar, while the time at which the growth sets in depends on the initial spin and tilt angle.

### iv.3 Evolution of the accretion rate

Fig. 10 shows the time evolution of the mass accretion rate for all models of our sample, in units of . The colors used for the different models as well as the line style follow the same convention as in the rest-mass density evolution figure (Fig. 9). The top panel displays the mass flux for models C1B, the middle panel for models NC1, and the bottom panel for models D2. The mass flux is computed as the instantaneous flux of matter across the AH surface which is parameterized in spherical coordinates by the polar angle and azimuthal angle and is given by:

 ˙M=2πr2∫2π0∫π0Dvrsinθdϕdθ, (26)

where is the radius of the AH and is the radial velocity of the fluid crossing the sphere. All three panels of Fig. 10 show that after a transient initial phase the mass accretion rate is seen to tend asymptotically to a fairly constant value.

In the previous fixed spacetime simulations of Fragile and Anninos (2005) it was found that the tilt and the twist of the disks are not strongly responsive to . Our results for a dynamical spacetime validate those earlier findings, as the final dispersion found in once this quantity reaches steady-state values is relatively small irrespective of the BH spin and inclination angles ( for models C1B and NC1, and for model D2). However, the transition to the steady-state does seem to show a clear dependence on the BH spin. We find that the initial values are consistently smaller the larger the BH spin, about 2 orders of magnitude for models D2 after 1 orbit and about half that value for the more massive disks C1B and NC1. Another trend, which is most clearly seen in model D2, is the initial dependence of the accretion rate on the tilt angle as well. We find that the higher the initial tilt angle, the higher the accretion rate in the early stages of the evolution. This is because the innermost stable circular orbit (ISCO) of a rotating Kerr BH has an ellipsoidal shape, attaining the largest size at the equatorial plane of the BH. For prograde orbits such as those in our work, the higher the BHs spin, the smaller the size of the ISCO and the closer can the matter orbit around the BH before falling in. When tilted, this centrifugal barrier becomes however less and less effective and the accretion rate becomes therefore higher for larger inclinations.

A common feature present in the mass-flux plots of the simulations of Fragile and Anninos (2005) was a late-time “bump” which the authors attributed to the reflection of an outgoing wave within the disk at the outer edge. That feature does not seem to be so apparent in our models apart from maybe model C1B which shows a bump around somewhat similar to that reported by Fragile and Anninos (2005). However, in our interpretation, the phase of rapid growth in in less than about 5 orbits seems to have a different origin. This conclusion is reached by comparing the time at which the growth of the PPI saturates (discussed in Figs. 11 and 12 in the following section) with the saturation time of the accretion rate for model C1B, which yields a significant close agreement. We note that this peak in the accretion rate is the highest of all our models, , lasting for less than about 1 ms. Therefore, we identify the exponential growth and subsequent saturation of the mass accretion rate with the growth and saturation of the PPI in our models.

The total rest mass accreted during the evolution ranges within , and for models C1B, NC1 and D2, respectively. The accretion is facilitated by the outward transport of angular momentum, which we describe in section IV.8 below. The growth of the PPI, to which we attribute the high peak accretion rates, is very effective in transporting angular momentum outwards. While we have performed our study without considering magnetic fields, it is known that magnetic fields are very effective at transporting angular momentum in accretion disks Balbus and Hawley (1991), via the development of the so-called magnetorotational instability Velikhov (1959); Chandrasekhar (1960). Not taking into account the effects of magnetic fields in accretion disks seems to underestimate the accretion rate Gold et al. (2014), see also recent GRMHD simulations of BH-NS mergers Paschalidis et al. (2015); Kiuchi et al. (2015) and the accretion rates reported therein.

### iv.4 PPI growth in the disk and BH movement

In their seminal paper on the PPI Papaloizou and Pringle (1984), the authors showed that vertically thick, radially slender tori with constant specific angular momentum profiles are unstable against the development of global non-axisymmetric modes. Two of our models, C1B and D2, have constant specific angular momentum profiles, while model NC1 has a radial dependence of , where . In Zurek and Benz (1986), the authors showed that tori with non-constant angular momentum profiles might become PP-unstable as well. In order to check for the onset and the growth of non-axisymmetric instabilities, we monitor the mode amplitudes computed using Eq. (36), which is basically a Fourier decomposition of the azimuthal distribution of the density in the disk Zurek and Benz (1986). In simulations of accretion tori susceptible to the development of the PPI performed on spherical grids, it is customary to induce the instability with small initial non-axisymmetric density or velocity perturbations. As outlined any non-axisymmetric perturbation will trigger the growth of the instability provided the torus is not stable against it. We do not actively seed the PPI in any of our models. Instead, there are two main sources of perturbations in our simulations, both connected to the fact that we evolve the BH-disk system in a Cartesian grid. First, the exact axisymmetry of the initial data is lost when we interpolate the initial data onto the Cartesian grid of the evolution, introducing small, non-axisymmetric perturbations in all interpolated quantities. This is due to the fact that in Cartesian grids neighboring cells lie on concentric shells only in the limit of infinite resolution. The second source of perturbation is the so-called junk-radiation leaving the BH while the gauge evolves from the values specified by the initial data to the puncture gauge attained during the evolution. Here, we need to distinguish between models evolved around Schwarzschild BHs and the models evolved around tilted Kerr BHs. For the former, the junk-radiation leaving the BH initially should be axisymmetric. Indeed, when analyzing this initial burst of radiation with the Weyl scalar , we find that the amplitudes of the multipoles (which are axisymmetric) are the largest. Nevertheless, the initial junk-radiation has also non-vanishing amplitudes in the other multipoles. Again, the reason for the presence of these non-axisymmetric multipoles seems to be connected to the evolution in the Cartesian grid, and specifically to the existence of mesh refinement boundaries, which are not axisymmetric even in the continuum limit (see Zlochower et al. (2012) for interference patterns produced by the radiation crossing mesh refinement boundaries). Finally, in the case of the disk being around a tilted Kerr BH, the junk-radiation reaching the disk is now truly non-axisymmetric even without taking into account the effects of the Cartesian grid and the mesh refinement boundaries, as the axis of rotation is tilted with respect to the plane.

We show in Fig. 11 the evolution of the amplitude of the first four non-axisymmetric modes in the disk for models C1Ba01b30 (top), NC1a01b30 (middle) and D2a01b30 (bottom). Those amplitudes are computed using Eq. (4). In addition, in the inset of all three plots we also show as a black solid curve the evolution of the polar distance of the center of the BH from the origin of the computational grid, which we label .

Model C1Ba01b30 is the one that most clearly is found to develop the PPI, as signaled by the exponential growth of the mode. Modes also show a phase of exponential growth lasting for about 2 orbital periods (between ) but the saturation amplitudes reached are 1 to 2 orders of magnitude smaller than for . As mentioned in the previous section, the growth of the PPI is accompanied by an exponential growth of the mass accretion rate (see Fig. 10). For this model the PPI clearly saturates at around , when the amplitude of the mode quickly drops by an order of magnitude. After that, the amplitude of all the modes we have analyzed remains rather constant. The strength of the modes becomes similar at late times while the mode still retains a significantly larger amplitude.

The polar distance of the BH from the origin for model C1Ba01b30 closely follows the behavior of the fastest growing PPI mode. This distance is found to grow at the same actual rate as the mode, because the circular motion of the overdensity hump (well visible in the top-rightmost panel of Fig. 7) causes the BH to start moving in a spiral trajectory. After mode saturation the polar distance of the BH from the origin remains almost constant, reflecting a similar trend in the evolution. This spiral motion of the BH as a result of the development of the PPI has also been observed in the simulations of Korobkin et al. (2011).

For the non-constant angular momentum model NC1a01b30 (middle panel) the growth rate of the mode is smaller, although it is still the dominant mode, as in model C1Ba01b30. The growth of the mode saturates somewhat later at around but, contrary to the evolution of C1Ba01b30, it does not drop significantly after saturation, remaining fairly constant thereafter. Furthermore, it remains the dominant mode at late times while the three remaining modes attain similar amplitudes. The evolution of the polar distance of the BH from the origin shows a small secular drift after the saturation of the PPI.

On the other hand, model D2a01b30 (bottom panel) shows a very distinct behavior. The amplitude of all four modes and the polar distance of the BH remain below the values attained in models C1Ba01b30 and NC1a01b30, typically in the range between 1 to 2 orders of magnitude. Furthermore, the and modes grow at approximately the same rate and to similar final values during the evolution. The growth rate of is very similar to the growth rate of the modes and . We also note that at early times the mode shows the largest amplitude (purple line; not so clearly noticeable in the previous two models), perhaps an artifact of the Cartesian grid we use in our simulations. Nevertheless, the amplitude of this mode at late times is sufficiently smaller than the dominant amplitudes (mostly that of the mode) to safely consider its effect on the dynamics negligible. From the mode analysis, we can therefore say that model D2 is essentially PP-stable. We will return to this point in section IV.8, where we analyze the angular momentum transport in the models.

The results we have just described connect the exponential growth and saturation of the mass accretion rate with the growth and saturation of the PPI in the models where it has unambiguously taken place. To further justify this, we plot in Fig. 12 the evolution of the mode for all models considered in our sample. Focusing on the evolution of the models C1B displayed in this figure (top panel) we can clearly see that the temporal order of the saturation of the mode growth closely coincides with the various peaks present in the mass accretion rate evolutions shown in Fig. 10 (compare curves of the same color). The comparison among the three types of models also shows that the saturation mode growth for models C1B is only slightly larger than for models NC1, their late-time values being fairly similar, but about 2 orders of magnitude larger than in models D2. The slope of the mode growth is also steeper in models C1B, while the behavior of the mode evolution in the case of the light tori D2 almost seems to indicate that the development of non-axisymmetric instabilities is only marginal for such models. A striking feature worth highlighting from the middle panel of Fig. 12 is the non-existing dependence of the late time mode amplitude on the BH spin and tilt angle for the non-constant angular momentum tori NC1.

To demonstrate the connection between the growth of non-axisymmetric modes in the disk with the movement of the central BH, we plot in Fig. 13 the position of the BH center on the plane for models C1Ba01b30 (top), NC1a01b30 (middle) and D2a01b30 (bottom). This figure shows several interesting features connected to the growth of the modes in the tori. The motion of the BH is caused by the formation of a “quasi-binary system” between the BH and the overdensity lump. It is this binary motion that causes the long-term emission of GWs reported in Kiuchi et al. (2011). We will return to this in the discussion of the GW emission connected to the PPI in section IV.9. For model C1Ba01b30, the development of the dominant structure in the disk exerts a small kick in the BH. As a result the BH starts moving in a spiral trajectory until the PPI saturates and the mode amplitude drops, at which point the BH is roughly located at . After saturation, the position of the BH is a combination of linear and circular motion. In model NC1a01b30, the linear motion is stronger from the beginning, causing a linear shift in the growing spiral motion. Upon saturation of the PPI, the circular motion continues with a smaller radius, in accordance to the evolution of the mode which remains at a rather constant value after saturation and up until the end of the evolution. Finally, for model D2a01b30 we have a superposition of circular and linear motion again, as the modes and are of almost equal magnitude during the evolution. Here, the overall distance covered by the BH is significantly smaller than for the heavier models C1Ba01b30 and NC1a01b30, as the mode amplitudes are much smaller for D2 and the disk is also lighter.

We argue that the non-axisymmetric instability we observe in our models is indeed the PPI based on the azimuthal mode evolution in the disk and the development of the overdensity lump. However, we do not claim to have presented a direct, undeniable proof that the instability observed in our simulations is the analytical PPI studied in Papaloizou and Pringle (1984). We also note that there exist other instabilities in self-gravitating disks which display spiral arm formation in the disk or even disk fragmentation, such as the axisymmetric local Toomre instability Toomre (1964), but as Christodoulou and Narayan (1992) have shown, the Toomre instability is unlikely to affect radially slender accretion tori such as the ones considered in our study.

### iv.5 PPI saturation and black hole kick

The complete 3D trajectory of the BH and its spin vector is plotted for model C1Ba01b30 in Fig. 14. In this plot each small sphere with an attached arrow corresponds to the position of the BH and the magnitude (scaled for better visualization) and direction of its spin vector. The positions are plotted at equal time intervals, so we see that the BH is moving much faster during the final stages of the PPI than before its occurrence or after its saturation. We can also see that there is a movement in the vertical direction, which is caused by the fact that the “binary” motion of the BH and overdensity lump does not lie in the equatorial plane anymore due to the initial tilt. While the PPI is growing, the BH-torus system therefore moves up and down in an oscillatory fashion. As we have seen in the mode plots for models C1B (see the top panel of fig. 12), the saturation of the PPI is very fast. The rapid saturation of the PPI (and the corresponding destruction of the overdensity lump) result in a mild vertical kick of the BH+torus system. To see that this is indeed connected to the saturation of the PPI, we plot in Fig. 15 the linear momentum radiated away by gravitational waves in the -direction for models C1Ba01. From the plot, we can clearly see that the time of maximum emission of linear momentum corresponds exactly to the time the PPI saturates for each model. The radiated linear momentum has been calculated with the pyGWAnalysis package Reisswig and Pollney (2011).

### iv.6 BH precession and nutation

Contrary to the test-fluid simulations of Fragile and Anninos (2005); Fragile et al. (2007), as we are evolving the full spacetime, we can monitor the response of the BH to the evolution of the disk in the initially tilted configuration. In particular we can measure the total precession of the BH spin about the -axis and its rate, as well as the inclination of the spin with respect to the -axis and the corresponding nutation.

In Fig. 16 we plot the evolution of the precession rate of the central BH about the -axis for models C1B (top), NC1 (middle) and D2 (bottom panel) in radians per orbital period. The first thing to mention is that the constant angular momentum disks D2 and C1B show a clear dependence of the precession rate with the initial BH spin magnitude. Namely, for models D2, the larger the spin the smaller the precession rate, and for models C1B smaller spins lead to higher precession rates. However, for the non-constant angular momentum torus NC1 the precession rate is very similar irrespective of the spin magnitude, particularly at the early stages (). In models D2 the evolution of is very smooth, and changes sign for the higher spin runs after about 17 orbits, whereas it remains fairly constant for the lower spin simulations.

A fact worth emphasizing is that there is essentially no dependence of the precession rate on the tilt angle for constant angular momentum tori. This is particularly true for models D2 but also models C1B show this lack of dependence up until the growth of the PPI. In all models C1B we observe a prominent modulation of the precession rate when the PPI enters the final stages of its growth and then saturates, as well as the change in sign for the higher spin runs. The magnitude of the precession rate is about an order of magnitude higher compared to models D2, which we attribute to the higher disk-to-BH mass ratio models C1B possess. The evolution of the BH precession rate for models NC1 seems to fit in between the two other models, the magnitude of the precession rate is somewhere in between but not very different evolutions are observed for the two BH spins considered. We also note from Fig. 16 that models C1Ba01 and D2a01 attain a fairly constant precession rate by the end of our simulations, namely and for models C1Ba01 and D2a01, respectively.

The total precession of the BH spin about the -axis is shown in Fig. 17. The lower spin models C1Ba01, that show the highest BH precession rate, have completed a quarter of a full precession cycle by the end of our simulations at . Given that the precession rate for those models appears to settle to a rather constant value (see Fig. 16) we would need a prohibitively expensive evolution of about orbital periods for a full precession cycle of the central BH for those models. If we assume that the precessing Kerr BH radiates GW as a freely precessing rigid body, it will radiate at and in the , multipole modes Zimmermann and Szedenits (1979). This means we would have to evolve the BH-torus system for at least one full BH precession cycle in order to see the slow modulation of the GW signal caused by the BH precession. It remains an interesting open issue to see if the upcoming advanced GW detectors will be able to measure these GWs in the low tens of Hz.

The imprint of the growth of non-axisymmetric instabilities in the disk on the BH response can be seen even clearer in Fig. 18. This figure shows the evolution of the nutation rate for all models, that is, the evolution of the temporal variation of the tilt angle of the BH with the -axis, . For models C1B (top panel) there is a clear period of alignment of the BH spin with the -axis around the time when the PPI starts to grow. This is signaled by the rapid fall and rise of the nutation rate around the mark for the slow spin () BHs, but it is also (less) visible in the more rapidly rotating BHs (dashed lines). Upon saturation of the instability the nutation rate becomes almost constant and fairly close to being again zero for the BH models (see solid lines). This behavior is consistent with the evolution of the BH tilt angle itself displayed in Fig. 19, which, for the case of models C1B with small spins, shows two distinct zones of constant spin inclination, particularly for models with and (model with still shows a negative slope at late times). The transition to zero nutation rate in the C1B models after PPI saturation seems to take a longer time than that displayed in Fig. 18.

In addition Fig. 18 also shows that the maximum nutation rate attained increases with increasing initial tilt angle. This is very clearly seen in models C1B but it is also visible in most of the models D2 and NC1. In particular the evolution of for the non-constant angular momentum models NC1 with (dashed curves) does not seem to reach a local maximum, at least in the timescales we can afford in our simulations, despite these models are also affected by the PPI. This might be explained by the fact that the non-axisymmetric mode does not drop as sharply after the saturation of the PPI in the models NC1 as it does in the C1B models (see Fig. 12). The growth of the mode in the disk seems to cause an alignment of the BH spin with the -axis. Furthermore, the nutation rate is much larger for the C1B models. Finally, models D2, in the absence of any strong growth of PPI modes, do not show strong modulation of the nutation rate. The modulation of this rate is therefore closely connected to the existence of non-axisymmetric modes in the disk.

The corresponding evolution of the inclination of the BH spin from the -axis is seen in Fig. 19. The panels show that a significant realignment of the BH spin and the -axis takes place by the end of the simulations in models C1B, i.e. in those that develop the PPI most significantly. For all models we see that the final inclination is smaller the bigger the initial tilt angle is, for equal spin magnitudes. This trend also seems to be more pronounced the higher the initial spin magnitude is.

### iv.7 Disk twist and tilt

We turn now to describe in this section the response of the disks as they evolve in the tilted spacetime, using the diagnostics we have introduced in Section II.4. In Figures 20 and 21, we plot the evolution of the tilt and precession of the total angular momentum vector of the disk . This vector is the sum of the angular momentum vectors of each shell, as described in Section II.4. We will refer to these as the global disk tilt and global disk precession around the BH spin, in order to distinguish them from the analysis of the twist and tilt angles in the individual disk shells.

The evolution of the global disk tilt, shown in Fig. 20, shows resemblances with the evolution of the tilt of the BH spin (see Fig. 19). For models C1B (top panel) we see that the development of the PPI not only causes a partial realignment of the BH spin with the -axis, but also a realignment of the BH spin and disk angular momentum. This alignment seems to become more constant towards the end of the simulation than the alignment of BH spin and the -axis. We also note that the amount of realignment is reversed for the global disk tilt for models C1B. Here the higher initial spin magnitude leads to a lower amount of alignment, in contrast with the evolution of the BH tilt. On the other hand, the alignment between BH spin and total disk angular momentum for models NC1 (middle panel) is not as pronounced as the realignment of the BH spin with the -axis. As in the case of models C1B, the trend of more alignment with higher spin for the BH tilt is inverted for the evolution of the global disk tilt. There are also oscillations in the evolution of the global disk tilt for models NC1, which are not visible in the other two set of models. These oscillations seem to be caused by the persistent non-axisymmetric mode in models NC1. Finally, for models D2 (bottom panel), neither the initial spin magnitude nor the initial tilt angle seem to affect the evolution of the global disk tilt, which is furthermore hardly changed during the entire evolution for all initial spins. The development of the PPI therefore seems to cause an alignment of the BH spin with the total disk angular momentum, which is larger for smaller initial spins. As both and decrease during the growth of the PPI, it is clear that the BH spin is tilted into the orbital plane of the disk.

In Fig. 21 we plot the global precession of the total disk angular momentum vector about the BH spin axis for models C1B (top), NC1 (middle) and D2 (bottom panel). For all models, the higher the spin (dashed curves), the larger the global precession, as expected. The initial tilt angle does not seem to influence strongly the evolution for models C1B and D2, while the evolution for NC1 shows oscillations towards the end of the evolution. Note that these oscillations are superimposed on the slow growth of the precession, which always has a positive slope, as expected. The oscillations are caused by the wobbling and smaller precession cycles about the direction of the global vector. These oscillations are then visible in the projection of the disk angular momentum vector onto the equatorial plane of the BH. The amplitude of the oscillations is larger the smaller the initial tilt angle. The non-axisymmetric modes that survive much longer in the case of models NC1 could be causing these oscillations in the global disk tilt and precession.

In Fig. 22 we show spacetime () diagrams of the radial disk tilt profile for the entire time evolution for all tilted models in order to see if and which disk models become warped during their evolution. The initial radial tilt profile is flat and the panels cover the initial radial extent of each model. For models C1B (top two rows) we see that by the end of the simulations the radial tilt profile for larger radii is significantly lowered compared to the initial tilt angle. Furthermore, there is a clear distinction in the radial tilt profile before and after the saturation of the PPI for the lower spin models (), in accordance with what we have observed for the global disk angular momentum tilt (see Fig. 20). For the higher spin models () the tilt profile evolution is broken up into two regions as well. In these cases, the peak close to the origin reemerges after the saturation of the PPI. Likewise, models NC1 (two middle rows) also show a peak close to the origin, but contrary to the previous set of models, the late time tilt profiles are oscillating in the regions closest to the central BH. The oscillations are stronger for both, radii closer to the origin and smaller initial BH spin. As described before, we believe that these oscillations are a consequence of the long-lasting non-axisymmetric mode in the disk.

The oscillations of the tilt angle close to the origin in models NC1 might be connected to the so called Kozai-Lidov (KL) mechanism Kozai (1962); Lidov (1962), where the eccentricity and inclination of particle orbits around an object that is itself in a binary are interchanged in an oscillatory fashion. Recently Martin et al. (2014) performed simulations of the KL mechanism in hydrodynamical disks around a component of a binary system, finding oscillations of the inclination angle. The persistent structure in the disks of models NC1 can be thought of as an eccentric distribution of matter. While our disks are not evolved around a BH in a binary system, we have nevertheless shown in Section IV.4 that the growth of the non-axisymmetric modes causes the BH to start moving in a “quasi-binary” with the overdensity lump in the disks. As seen in Fig. 13, for models NC1, the long lasting “planet” causes the BH to move in a continuous elliptic fashion in the projection onto the -plane. We believe this movement of the BH induced by the PPI serves as the correct boundary conditions for the KL mechanism to occur in models NC1. In Martin et al. (2014) the authors estimate the timescale of the KL oscillations in their disks as being , where is the orbital period of the central object in the binary. In our case, the periods of the tilt angle oscillation and the orbital motion of the Kerr BH are roughly equal for all models. We plan to further investigate the correctness of this interpretation of the tilt angle oscillations in a future work.

To complete the analysis of Fig. 22, we note that models D2 (two bottom rows) become clearly warped, with all models showing a peak near the origin. The profiles of the low spin models D2a01 are very similar in shape for the three initial tilt angles. The peak in models D2a05 is seen to oscillate stronger the smaller the initial tilt angle. Note that for none of the models we see an alignment of the inner region of the disk with the equatorial plane of the BH (this would mean there). We therefore find no occurrence of the Bardeen-Petterson effect in our simulations, at least in the timescales considered. This is in agreement with the results found in the fixed-spacetime inviscid simulations with no angular momentum transport of Fragile and Anninos (2005) and also with the GRMHD simulations of Fragile et al. (2007) that included angular momentum transport driven by the magneto-rotational instability. The observed profile and evolution of the disk warp seems to be in qualitative agreement with the analytic work on warp propagation as bending-waves in tilted thick accretion disks; see, for instance, the analytic tilt profiles in Lubow et al. (2002), the analytic model of linear warp evolution of Foucart and Lai (2014) and the application of the linear warp evolution model in Franchini et al. (2015). We note that while the Bardeen-Petterson effect is not manifest in our models, we nevertheless do observe a global partial realignment of the BH spin with the disk angular momentum, caused by the growth of the PPI. This alignment of the disk angular momentum with the BH spin has also been observed in the post-merger evolution of a tilted accretion disk in Kawaguchi et al. (2015). The authors conclude that the transport of angular momentum by non-axisymmetric shock waves in the disk might be responsible for such a Bardeen-Petterson-like effect. Elucidating this issue deserves further studies.

Finally, to check for Lense-Thirring precession, we plot the evolution of the inner region disk twist from orbits onwards for all our tilted models in Fig. 23. We see that all models become twisted (as varies with ) in the regions for radii up to . The plots further indicate solid precession of the disks, as the radial profile of is increasing smoothly in time almost independently of the radius in the outer regions of the disk and does not remain for large radii. This is consistent with the solid-body precession found in the fixed spacetime simulations of Fragile and Anninos (2005). Models D2 actually show a slightly growing twist for larger radii. In this trend, the higher initial spin models show flatter profiles, while there is some radial modulation of the profiles for the lower spin simulations. A striking feature of the plots is the much smaller difference in accumulated twist for large radii between different initial spin magnitude for models NC1. In models C1B and D2 the higher spin models show a much higher twist at large radii than their low spin counterparts. As the long survival of the mode in the NC1 models seems to have an effect on the evolution of the twist, we plot this evolution only up to the saturation of the PPI for these models in order to properly visualize the solid body precession of the disk.

### iv.8 Angular momentum transport

In this section we investigate the transport of angular momentum in the disks. To illustrate the transport of angular momentum during the evolution of our models, we plot spacetime diagrams showing the magnitude of the angular momentum in radial shells, computed using Eq. (9). Such diagrams are displayed in Fig. 24. The magnitude of the angular momentum in each shell, , has been rescaled to the global maximum value attained in any shell during the evolution. Fig. 24 shows the evolution of for the three untilted models C1Ba0b0 (top panel), NC1a0b0 (middle panel) and D2a01b0 (bottom panel). This figure reveals that compared to the latter two models, model C1Ba0b0 shows a radical redistribution of angular momentum in the radial direction. This drastic difference is due to the development of the PPI in model C1Ba0b0. Following the initial perturbation we see that the contours of start growing in a wave-like manner, transporting angular momentum outwards. This happens up until the saturation of the PPI at about 10 orbital periods. As noted by Balbus (2003) during the growth phase the PPI is thought to be an effective mechanism for angular momentum transport in thick accretion disks, where a right combination of rotation and pressure at the boundaries of the disk allows the growing non-axisymmetric modes to transport angular momentum outwards. In Zurek and Benz (1986), the authors showed as well that the development of nonaxisymmetric instabilities is closely connected to the transport of angular momentum. The authors observed that disks with a larger than the critical value given in Papaloizou and Pringle (1984), , were still unstable, although the angular momentum transport was slower in those cases. After saturation, the new distribution of angular momentum seems PP-stable and does not exhibit drastic changes for the rest of the evolution.

We can compare the evolution of to the development of the non-axisymmetric mode in the disk, which drops to a value of around of its maximum value after saturation (see top panel of Fig. 12). As described in Balbus (2003), the development of the PPI crucially relies on the properties of both the inner and outer boundaries of the disk. In Blaes (1987) it was shown that a small amount of accretion (which modifies the inner boundary of the disk) was sufficient to saturate the growing PPI. Furthermore, wider tori are known to be not as violently unstable as slender tori Hawley (1991). In subsections IV.3 and IV.4 we already showed that the accretion rate and the growth of the non-axisymmetric mode saturate at the same time, but the diagram of Fig. 24 also shows that the outer angular momentum boundary of the disk is pushed outwards with growing amplitude as well. It therefore seems that both mechanisms are saturating the PPI simultaneously in model C1Ba00.

Model NC1a0b0, displayed in the middle panel of Fig. 24, has a non-constant specific angular momentum distribution, and shows a very different evolution of . The inner region of the disk shows a gentle reduction of angular momentum during the entire evolution and no sudden and drastic redistribution is found as occurred in model C1Ba0b0. The outer contour grows to encompass the entire outer domain shown in the plot. Another difference is the almost complete absence of variability of the inner region of the disk, which is in contrast to model C1Ba0b0, where oscillatory behavior is seen in the inner boundary of the disk too. The specific angular momentum profile of model NC1a0b0 is between constant and Keplerian, with Keplerian disks being essentially PP-stable. Model NC1a0b0 still develops the PPI, albeit in a more mild manner than model C1Ba0b0, with the torus of model NC1a0b0 being initially wider.

Finally, model D2a01b0 (bottom panel), the lightest of our models, is PP-stable as already shown in subsection IV.4. The corresponding spacetime diagram of shows no pronounced outward transport of angular momentum as the evolution proceeds, the profile remaining very similar to its initial configuration. The inner region with a high angular momentum magnitude is thinned out slightly, and the profile spreads a bit overall, but the changes are nowhere near as pronounced as for models C1Ba0b0 and NC1a0b0.

As we have seen in the analysis of the modes, the amplitude of the mode drops sharply upon the saturation of the PPI in the C1Ba00 model. The subsequent evolution of the angular momentum profile remains more constant, similar to what we observe for the entire evolution of the PP-stable model D2a01b0. There is a clear split between the stages of the evolution: the first during which the PPI develops and saturates, and the subsequent evolution of an essentially PP-stable torus.

Model NC1a00, on the other hand, shows a very different behavior: as already seen in the mode analysis, the mode amplitude remains roughly at its saturation value for the rest of the evolution, which is lower than the values attained by models C1B. The persistent mode is seen to continuously transport angular momentum outwards till the end of the simulation. The restructuring of the disk by the transport of angular momentum is more long-lived, compared to the drastic change the saturation of the PPI brings about in model C1Ba00. These findings seem to confirm those found in Zurek and Benz (1986). Using the same kind of spacetime diagrams we have also checked the dependence of the angular momentum transport on the BH spin and on the initial tilt angle for all models of our sample. We find that the angular momentum transport shows weak dependence with the tilt angle (the higher the tilt angle the slowest the transport) and shows essentially no dependence with the spin, at least for the moderate values of the BH spin we could afford in our study. This is a remarkable result, as it demonstrates that models C1B and NC1 are PP-unstable for a wide range of initial BH spin magnitudes and tilt angles.

### iv.9 Gravitational waves

The recent numerical simulations of Kiuchi et al. (2011) have shown that the growth and saturation of the PPI makes BH-torus systems strong emitters of large amplitude, quasi-periodic gravitational waves, potentially detectable by forthcoming ground-based and spacecraft detectors. It was found in particular that the non-axisymmetric structure survives with an appreciable amplitude after the saturation of the PPI nonlinear growth, which leads to the emission of quasi-periodic GWs with a large amplitude. (First estimates through full GRHD simulations of the GW detectability by the PPI instability in self-gravitating BH-torus systems were obtained in Korobkin et al. (2011) where the nonlinear saturation phase was almost reached - see also van Putten (2001) for earlier proposals and alternative estimates of the gravitational waves emitted by non-axisymmetric instabilities in such systems.) In this section we turn our attention to analyze the implications of the PPI in our tilted BH-torus systems regarding the emission of GWs, namely the dependence of the resulting gravitational waveforms and spectra on the BH spin magnitude and initial tilt angle.

In Fig. 25 we plot the real part of the mode of the outgoing part of the complex Weyl scalar for models C1Ba01 and NC1a01 and for all initial tilt angles . This quantity has been computed at an extraction radius . Fig. 25 shows that, in agreement with Kiuchi et al. (2011), all models display strong emission of gravitational waves and that the emission persists for many dynamical timescales well after the saturation of the PPI. The peak amplitude for models C1Ba01 is about an order of magnitude higher than for models NC1a01. We recall that model C1B has a slightly higher initial disk-to-BH mass ratio than model NC1 ( and , respectively) and that the development of the PPI in models C1B is more pronounced than in the non-constant specific angular momentum models NC1 (see Fig. 7). These two features are responsible for the difference in the peak amplitudes between the two initial models. We also note that the peak amplitude in the C1Ba01 models, reached at PPI saturation, is much higher than the amplitude of the remaining signal, whereas in models NC1a01 the variations are not so pronounced. As we showed in Fig. 12, upon PPI saturation the dominant PPI mode drops to about of its peak value in models C1B, while it remains at a similar strength for models NC1.

Fig. 25 also shows that the GW signal has a weak dependence with the initial tilt angle, particularly for the non-constant specific angular momentum models NC1. The smallest peak amplitudes are attained for the most tilted BH spacetime (). In any event, the differences found in the values spanned by the peak amplitudes with regard to are not too significant.

We use the fixed frequency integration (FFI) described in Reisswig and Pollney (2011) in the integration of the data to obtain the GW strain. In Figs. 26-28 we plot the effective strain as a function of frequency for all models, placing the sources at 1 MPc, in order to see any possible imprint of the initial spin magnitude and tilt angle on the GW spectrum. As described above, we use the same time-step in the 7 outermost refinement levels ( and 0.64 for higher and lower resolution runs, respectively), which means that the time-step at the extraction radius of is sufficiently small to sample the GW multipole signal with a of 1.28 (2.56) for the higher (lower) resolution runs, respectively. Due to the Nyquist criterion, we can therefore resolve the GW spectra up to (40) kHz. The top panels of Figs. 26 and 27 correspond to the mode and the bottom panels to the mode. The highest power is usually found at a frequency of about kHz. For the effective strain of the mode for models C1B shown in the top panel of Fig. 26, we can clearly identify two trends. On the one hand there is a difference at high frequencies between the models with either no or low () initial BH spin (models C1Ba00 and C1Ba01) and the high spin models (). The spectra for these two groups show a different slope from onwards. On the other hand, for models C1Ba00 and C1Ba01, an increase in the initial tilt angle leads to an increase in the power for frequencies . This also seems to be present at the beginning of the spectrum although it is somewhat less clear.

Correspondingly, the power spectra of the mode shown in the bottom panel of Fig. 26 neatly splits the models in three groups, according to the initial spin magnitude of the BH. We find that the slope of the spectra becomes steeper from onwards the higher the BH spin. In all cases, the slope of the effective strain with the frequency shows essentially no dependence with the initial tilt angle for the various groups of models characterized by the same BH spin. The power in model C1Ba00 is at least an order of magnitude larger for all frequencies than that of the tilted models. This is because the spiraling movement in the equatorial plane of both the BH and the overdensity PPI “planet” emit GWs predominantly in the mode which are emitted in a direction perpendicular to the equatorial plane. This is an optimal situation for the model with aligned BH and disk spins. However for the tilted models the “orbit” of the BH and of the overdensity lump in the disk is not confined to the equatorial plane anymore because of the reaction of the disk to the tilted BH spin. This results in a decrease in the GW power.

The spectra of the non-constant specific angular momentum models NC1 displayed in Fig. 27 is markedly different to those of constant specific angular momentum tori. Namely, their dependence on the frequency is similar for all modes (same slope) irrespective of the initial tilt angle and of the BH spin. No split in two or three families is found according to the BH spin. The effective strain of the mode for models NC1 shows a plateau in the low frequency part (up to ) for models NC1a03b15 and NC1a03b30. The spectra shows a prominent peak around for all models but these two. For the mode we find that the curve for the untilted model NC1a00 lies more than an order of magnitude above the curves corresponding to the initially tilted disks, as seen also for the C1Ba00 model in the previous figure.

Finally, in Fig. 28 we show the effective strains of the mode (top panel), (2,2) mode (middle panel) and (3,1) mode (bottom panel) for models D2. In the spectra of the mode there is a clear variation with the initial tilt and BH spin magnitude, with both increasing the power of the mode for all frequencies. Contrary to model C1B and similar to model NC1, we find now that model D2 does not show the grouping of models depending on the BH spin, and the associated frequency dependence (i.e. steeper slopes for higher spins). As this model is PP-stable and is built with an initially constant distribution of specific angular momentum, its dynamics is somewhat between that of the other two models. The spectra of the mode (middle panel) show in particular a very different behavior to those shown for models C1B and NC1, namely that the model with zero BH spin does not show significantly larger power than the models with and . Since models D2 do not develop the PPI, there is no strong emission of GW in the mode due to the absence of the spiral motion of both BH and overdensity ‘planet’ that were observed for the other two initial models. Therefore, the spectrum of the untilted model D2a01b0 does not show a higher power density than the rest of the tilted D2 models. The two models D2a05b15 and D2a05b30 both show two peaks in the low frequency part of the spectrum, as well as a flatter slope between and .

The absence of the PPI in the D2 models results in a GW emission where the mode is not the dominant mode. To show this, we plot the effective strain of the mode in the bottom panel of Fig. 28. The power is comparable with that of the and modes, which is not the case for models C1B and NC1 that both develop the PPI. The plot of the effective strain of the mode shows that same trend as that of the mode, i.e. the power density increases with spin magnitude and initial tilt angle. Furthermore, we can see clear quasi-periodic oscillations in the low frequency part of the spectrum for models D2a05, as well as a developing plateau which becomes flatter with higher initial tilt angle from .

In order to estimate the detectability of the GWs emitted by these systems, we follow the analysis in Kiuchi et al. (2011). We assume that the BH-torus system of models NC1 will radiate GWs for the entire accretion timescale (which is the time needed for the entire disk to accrete). This is because in these models, the overdensity lump and the BH form a long-lasting quasi-binary system. Assuming that the accretion rate will remain at the levels it attained by the end of the simulations (see Fig. 10), the lowest estimate for the accretion timescale for models NC1 is . From this timescale, we can estimate the number of GW cycles at the peak frequency, which yields cycles as a lowest estimate for models NC1. The peak amplitude of the GW strain will then be amplified by , while we assume no such amplification for neither models C1B (where the PPI is rapidly damped after saturation) nor for models D2 (which are PP-stable). We note that these estimates are in very good agreement with the corresponding findings reported in Kiuchi et al. (2011). The corresponding peak amplitudes for the BH-torus systems located at , together with the advanced LIGO (aLIGO) sensitivity curve, are plotted in Fig. 29, showing that the GWs emitted by models NC1 could be detectable due to the long-lasting emission of GWs due to the persistent mode. Note that the planned factor-of-3 upgrade of aLIGO Hild (2012); Adhikari (2014); Miller et al. (2015) will further improve the detectability.

## V Discussion

In this paper we have presented the results of an extended set of numerical relativity simulations of massive accretion disks around tilted Kerr black holes. We have considered three different thick accretion disks of varying mass and specific angular momentum magnitude and distribution, along with different black hole configurations with two spin magnitudes and four initial tilt angles . On the one hand the motivation for this work has been to extend the investigations of Fragile and Anninos (2005); Fragile et al. (2007) in the test-fluid approximation to full general relativity, analyzing the effects of the self-gravity of the disk on the BH-torus dynamics. On the other hand, our work has also served to enlarge the parameter space of the existing numerical relativity simulations Korobkin et al. (2011); Kiuchi et al. (2011) by accounting for tilted BH-torus systems for the first time.

For the models with disk-to-BH mass ratios of we have studied, we have found that the assumption of using a fixed background spacetime is not complete as we have observed significant precession and nutation of the tilted black hole as a result of the disk evolution that cannot be accounted for in fixed spacetime simulations. For some of our models the precession rate attains a fairly constant value by the end of the evolution. The BH nutation rate is also seen to be drastically modulated for those models that develop the PPI, showing that the development of non-axisymmetric modes in the disk exerts a torque on the central BH.

Our simulations have revealed that the development of the PPI seems to be universal for the range of initial spin magnitudes and tilt angles investigated. Models C1B and NC1, both PP-unstable when the central BH is non-rotating or untilted, remain PP-unstable when the BH is rotating or tilted. The reverse also seems to be true: model D2, PP-stable for an untilted, non-rotating BH remains PP-stable for all spins and initial tilt angles considered. As we mentioned before, model D2, having the smallest mass ratio in our study, is the most likely outcome of BH-NS mergers based on astrophysical considerations. Thus, our findings are relevant to gauge the practical role of the PPI in thick post-merger accretion disks around black holes. In agreement with the numerical relativity simulations of Kiuchi et al. (2011) the growth of the PPI mode manifests itself in the formation of a counter-rotating overdensity lump (or “planet”) that forms a “quasi-binary” with the central BH during its existence. This causes the BH to start moving in a spiral trajectory for as long as the “planet” exists. For models C1B the “planet” disperses quickly upon the saturation of the PPI, while in models NC1 the overdensity structure persists much longer, with the mode amplitude remaining at its level of saturation. When the PPI develops around a tilted Kerr BH, the binary motion of BH and “planet” is not restricted to the -plane but also shows some vertical motion. For models C1B in particular, the rapid destruction of the non-axisymmetric structure causes a mild kick to the BH-torus system in the vertical direction which is also imprinted in the time evolution of the linear momentum radiated away by gravitational waves in the same direction.

The evolution of the disk around the tilted BH can cause a significant twisting (differential precession) and warping (spatially varying tilt) in the disk. We have monitored the evolution of the twist and tilt of radial shells during the simulations. For models C1B (NC1) we have found a phase of rapid (mild) realignment of the total angular momentum vector of the disk, , and the BH spin, , during the growth of the PPI. We attribute this alignment to the development of the structure in the disk. Our simulations have also confirmed the presence of significant differential twisting of the disks due to Lense-Thirring precession. For all models, the cumulative twist is higher for higher initial BH spins, as expected, and the outer regions of the disks precess as a solid body. The evolution of the tilt profiles of models D2 is similar to what Fragile and Anninos (2005) observed, with the development of a noticeable peak in the innermost region of the disk. For models C1B the PPI causes a drastic change in the tilt profile upon saturation. Correspondingly, models NC1 develop strong tilt oscillations in the regions close to the BH after the PPI saturates and the structure persists till the end of the simulations. By interpreting the long-lived “planet” as an eccentric distribution of matter, we suggest that the so-called Kozai-Lidov mechanism, in which eccentricity and inclination are exchanged in an oscillatory manner, might explain the tilt angle oscillations we observe, in analogy with the recent findings of Martin et al. (2014). In none of our models we observe the Bardeen-Petterson effect, as in the inner region of the disks.

The PPI is thought to activate the outward transport of angular momentum. This is indeed the case in our simulations, as the evolution of the angular momentum profiles shows. For model C1Ba0b0, the angular momentum profile contours start growing in a wave-like manner up until the sudden saturation of the PPI. After saturation there is no more transport of angular momentum and the evolution resembles that of the PP-stable model D2a01b0. For model NC1a0b0 the persistence of the structure in the disk continuously transports angular momentum outwards for the entire evolution of the disk.

As Kiuchi et al. (2011) showed, the development of the PPI and the corresponding overdensity lump in the disk cause the long-term radiation of gravitational waves, predominantly in the multipole mode, as expected from the radiation emitted by a binary system. Our simulations have confirmed these results, showing that models C1B and NC1 do indeed radiate mainly in that particular mode. The rapid destruction of the structure in models C1B causes the amplitude of the emitted GW signal to drop significantly after the PPI saturates, while the peak amplitude is closely correlated with the time of saturation. Models NC1 on the other hand, emit at an amplitude similar to that attained at saturation, as the structure survives until the end of the simulations. We have also calculated the effective strains of the GW signals. While models C1B show different spectra for different initial spin magnitudes and tilt angles, there seems to be no such trend for models NC1. For those two set of models, the GW emission is predominantly in the mode while in models D2 there is also significant power in the , mode. In addition, the strain spectra of models D2 show a clear trend to higher values with increasing initial spin and tilt angle. Finally, we have briefly touched on the issue of the possible modulation of the GW signal caused by the BH precession. For the fastest precessing models of our sample, models C1Ba01, we have shown that we would need about 80 orbits for a complete BH precession cycle in order to see such an effect, assuming the precession rate remains constant beyond the 20 orbit mark we could afford in our simulations. We plan on exploring such a long-term simulation in the future in order to obtain the complex GW signal from the kind of precessing BH-torus systems we have explored in this work.

###### Acknowledgements.
We would like to thank the anonymous referees for their very helpful comments and suggestions to improve the manuscript. It is a pleasure to thank Chris Fragile for his careful reading of the manuscript and for providing useful comments. VM would like to thank the developers of the Einstein Toolkit for their invaluable support during the duration of the project. He would also like to thank Ewald Müller and the MPA for hospitality and support during his two research visits. FG thanks Luciano Rezzolla for useful discussions on the gravitational wave signals. FG gratefully acknowledges financial support from the “NewCompStar” COST Action MP1304 during his visit to the University of Valencia. This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) through grants AYA2010-21097-C03-01 and AYA2013-40979-P, by the Generalitat Valenciana (PROMETEOII-2014-069), by the Deutsche Forschungsgemeinschaft (DFG) through its Transregional Center SFB/TR7 “Gravitational Wave Astronomy”, and by the Helmholtz International Center for FAIR within the framework of the LOEWE program launched by the State of Hesse. The simulations were performed on the Hydra HPC system at the Rechenzentrum Garching.

## Appendix A The BSSN evolution equations

The BSSN formulation introduces a set of new auxiliary variables for the evolution, defined in terms of the ADM quantities: the conformal factor

 ϕ≡ln(112detγij), (27)

the conformal spatial metric

 ~γij≡e−4ϕγij, (28)

the trace of the extrinsic curvature

 K≡γijKij, (29)

the conformal trace-free extrinsic curvature

 ~Aij≡e−4ϕ(Kij−13γijK), (30)

and the conformal connections

 ~Γi≡~γjk~Γijk. (31)

Using these auxiliary variables together with appropriate gauge choices leads to a strongly hyperbolic system Beyer and Sarbach (2004) that allows for long-term stable numerical evolutions.

The evolution equations used in McLachlan are the following:

 ∂0α = −α2f(α,ϕ,xμ)(K−K0(xμ)), (32) ∂0K = −γij~Di~Djα+α(~Aij~Aij+13K2) (33) +4π(E+γijSij), ∂0βi = α2G(α,ϕ,xμ)Bi, (34) ∂0Bi = e−4ϕH(α,ϕ,xμ)∂0~Γi−ηi(Bi,α,xμ), (35) ∂0ϕ = −16(αK−∂kβk), (36) ∂0~γij = −2α~Aij+2~γk(i∂j)βk−23~γij∂kβk, (37) ∂0~Aij = e−4ϕ(α~Rij+αRϕij−~Di~Djα)TF (38) +αK~Aij−2α~Aik~Akj+2~Ak(i∂j)βk −23~Aij∂kβk−8παe−4ϕSTFij, ∂0~Γi = −2~Aij∂jα (39) +2α(~Γikl~Akl+6~Aij∂jϕ−23~γijK,j) −~Γj∂jβi+23~Γi∂jβj+13~γikβj,jk+~γjkβi,jk −16πα~γikSk,

where we have used the following shorthand notation, .

The coupling to the stress-energy tensor is done via the usual projections with the spatial metric and normal vector ,

 E ≡ nμnνTμν (40) = 1α2(T00−βiT0i+βiβjTij), (41) Sij ≡ γiμγjνTμν, (42) S ≡ Sii=γijSij, (43) Si ≡ −γiμnνTμν=−1α(T0i−βjTij). (44)

This system constitutes the so-called -variant of the BSSN formulation (see Pollney et al. (2011) for other possible variations).

McLachlan uses the 1+log slicing for the evolution of the lapse Bona et al. (1995):

 f(α,ϕ,xμ) ≡ 2α, (45) K0(xμ) ≡ 0, (46)

and the following -driver shift condition Alcubierre et al. (2003) for the evolution of the shift:

 G(α,ϕ,xμ) ≡ 34α2, (47) H(α,ϕ,xμ) ≡ e−4ϕ. (48) ηi(Bi,α,xμ) ≡ ηBiq(r)