The Shape of an Accretion Disc in a Misaligned Black Hole Binary

We model the overall shape of an accretion disc in a semi-detached binary system in which mass is transfered on to a spinning black hole the spin axis of which is misaligned with the orbital rotation axis. We assume the disc is in a steady state. Its outer regions are subject to differential precession caused by tidal torques of the companion star. These tend to align the outer parts of the disc with the orbital plane. Its inner regions are subject to differential precession caused by the Lense-Thirring effect. These tend to align the inner parts of the disc with the spin of the black hole. We give full numerical solutions for the shape of the disc for some particular disc parameters. We then show how an analytic approximation to these solutions can be obtained for the case when the disc surface density varies as a power law with radius. These analytic solutions for the shape of the disc are reasonably accurate even for large misalignments and can be simply applied for general disc parameters. They are particularly useful when the numerical solutions would be slow.


INTRODUCTION
We consider a semi-detached binary system with an accretion disc around a black hole the spin axis of which is misaligned with the orbital plane of a companion star. This can arise in a binary system when one star has an asymmetric supernova explosion which gives a velocity kick to the collapsing compact object (Shklovskii 1970;Sutantyo 1978;Dewey & Cordes 1987;Bailes 1989). Even a small asymmetry in the supernova can give a large velocity kick to the newly formed black hole or neutron star (Brandt & Podsiadlowski 1995;Martin, Tout & Pringle 2009).
The tidal torque of the companion causes each ring in the disc to precess about an axis parallel to the binary orbital axis. The precession rate increases with increasing radius from the black hole. This torque can be balanced by the viscosity in the disc to create a steady state disc if the orbital period of the companion is short compared with the viscous timescale in the outer parts of the disc. The disc becomes warped and twisted when the inner regions of the disc remain aligned with the spin of the central black hole.
Microquasars are binary systems in which material is accreted from a normal star on to a compact object with relativistic radio jets (Mirabel & Rodríguez 1999). The inner accretion disc aligns with the spinning black hole by the Bardeen & Petterson (1975) effect. Jets produced in this region are expected to be aligned with the spin of the black hole. We can measure the inclination of the jets and the inclination of the binary orbit and find the misalignment angle. GRO J1655-40 might be misaligned by an angle of at least 15 • . Hjellming & Rupen (1995) measured a jet inclination of 85 • ± 2 • to the line-of-sight and Greene, Bailyn & Orosz (2001) measured the orbital axis inclination to be 70.2 • ± 1.9 • . For the microquasar V4641 Sgr (SAX J1819.3-2525), Orosz et al. (2001) deduced an orbital inclination angle in the range 60 • −70 • along the line-of-sight and a jet inclination of less than 10 • so this microquasar is thought to be misaligned by at least 50 • . Martin, Tout & Pringle (2008) and Martin, Reis & Pringle (2008) examined the microquasars GRO J1655-40 and V4641 Sgr. They used the steady state disc model of Martin, Pringle & Tout (2007) warped by the Bardeen & Petterson (1975) effect but included no binary torque term so their discs must effectively extend to infinity in order to flatten at the outside. Such a disc is aligned with the spinning black hole at the centre and at the outside the inclination tends to that of the incoming material. We consider here what effect a binary torque has on such a steady state disc.

WARPED DISC MODEL
We model the evolution and structure of warped accretion discs according to the simplified formalism described by Pringle (1992). We consider the disc to be made up of annuli of width dR and mass 2πΣRdR at radius R from the central star of mass M1 with surface density Σ(R, t) at time t. The angular momentum at radius R is L = (GM1R) 1/2 Σl = Ll where l = (lx, ly, lz) is a unit vector describing the direction of the angular momentum of a disc annulus with |l| = 1.
We use equation (2.8) of Pringle (1992) setting ∂L/∂t = 0 and add a term to describe the Lense-Thirring precession (if the central object is a black hole) and a binary tidal torque term to give The Lense-Thirring precession term is given by where ωLT = 2GJ c 2 (3) (Kumar & Pringle 1985). The angular momentum of the black hole J = Jj, with j = (jx, jy, jz) and |j| = 1, can be expressed in terms of the dimensionless spin parameter a such that This approximation of the torque is valid if ΩLT ≪ aΩ where Ω 2 = GM/R 3 (Kumar & Pringle 1985). This corresponds to R ≫ Rg = GM/c 2 . The inner edge of an accretion disc is the innermost stable orbit around the black hole. For a nonrotating black hole this is at 6 Rg and decreases down to Rg for a maximally rotating black hole. The binary torque term is given by We discuss this torque further in the next Section. There are two viscosities. First ν1 corresponds to the azimuthal shear (the viscosity normally associated with accretion discs) and secondly ν2 corresponds to the vertical shear in the disc which smoothes out the twist. The second viscosity acts when the disc is non-planar. With the α-parameterisation the viscosities are given by ν1 = α1csH and ν2 = α2csH where cs is the sound speed and H is the scale height of the disc. The dimensionless parameters α1 1 and α2 must be determined experimentally. Papaloizou & Pringle (1983) find that the viscosities can be related by to first order. The more detailed analysis of warped discs by Ogilvie (1999) and Ogilvie (2000) supports this approach. They include, into equation (1), an extra precession term of the form where ν3 is a third viscosity parameter related to the precessional effects with ν3 = α3csH. However, Lodato & Pringle (2007) find that α3 ≪ α2 and we may neglect this term. In order to obtain analytic results we take both viscosities to have a power law form so that where ν10, ν20 and β are all constants and R0 is some fixed radius. In a steady flat (unwarped) disc the surface density far from the inner edge is then where x = 0 for a steady accretion disc. For a steady decretion disc in which matter is added at a steady rate at the inner radius x = 1/2. This value of x gives the radial dependence of the surface density at late times in the evolution of a decretion disc of fixed mass. An accretion disc reaches steady state on the viscous timescale given by As long as the mass transfer rate remains constant over this timescale then the disc can reach this steady state solution.

TIDAL TORQUE
Lubow & Ogilvie (2000) and Ogilvie & Dubus (2001) find the tidal torque in the frame of the binary orbit, of radius R b , to be approximately to first order in the angle of tilt between the disc and the orbital plane. Here ez is a unit vector in the direction of the binary orbital axis, M2 is the mass of the binary companion star and R b is the binary separation. The Laplace coefficient of celestial mechanics can be approximated for small z by We note that the tidal torque satisfies l.T tid = 0 so that the tidal torque is zero when the disc and orbit are aligned. Thus even for large tilt angles the tidal torque term is likely to have this form, although with modified coefficients. This tidal torque is averaged over the orbital period of the system which we assume to be much shorter than the viscous timescale at the outer edge of the disc.

RELEVANT DISC RADII
We expect the centre of the disc to be aligned with the spin of the central black hole by the Bardeen-Petterson effect. We expect the tidal torques to align the outer regions with the binary orbit by an analogous physical process. We note that if the outer parts of the disc are closer to counter-alignment than alignment, then the whole disc is counter-rotating. The inner parts of the disc are flat and counter-aligned with the black hole spin. The counter aligned disc solution is oppositely twisted to an aligned solution (Scheuer & Feiler 1996;Martin, Pringle & Tout 2007). The torque on the black hole is always towards alignment of the black hole with the outer parts of the disc because the angular momentum of the disc is effectively very large compared to that of the hole (King et al. 2005). This is because the angular momentum of the binary orbit is transferred to the outer parts of the disc when the tidal torque acts (Martin, Tout & Pringle 2008).
By comparing the sizes of various terms in the governing equation (1) we are able to get an initial idea of the structure that the disc is likely to have when it is in a steady state. In either case, aligned or counter aligned discs, there are three radii of physical significance.

The Lense-Thirring radius RLT
A Lense-Thirring radius RLT where the effects of Lense-Thirring precession and viscosity balance was defined by Scheuer & Feiler (1996) and Martin, Pringle & Tout (2007) by balancing the first term of equation (1) involving ν2 with the Lense-Thirring term and replacing ∂/∂R with 1/R. Well inside this radius we expect the disc to be flat and aligned with the spin of the central black hole.

The tidal radius R tid
We define in a similar way a tidal radius, R tid , to be where the tidal term balances the same viscous term in the disc so that If a disc extends beyond this radius then the viscous effects are dominated by the binary torque term. Thus we expect that well outside this radius the binary torque term, coupled with the viscosity, flattens the disc and aligns with the binary orbital plane. The precise form of R tid will become clearer in Section 6.2 when we use it to simplify equation (26).

The warp radius Rwarp
The overall warping of the disc is caused by the misalignment between the Lense-Thirring precession which dominates at small radii and the tidal precession which dominates at large radii. The rate of warping of the disc with radius is likely to be greatest when these two precessional effects balance. We call this the warp radius. It is It is independent of the magnitude of the viscosity because we are now balancing just the two external torques.
We may then write Thus over all we expect the disc to be warped at radii between the Lense-Thirring radius RLT and the tidal radius R tid , with the rate of warping greatest at around the warp radius Rwarp.

NUMERICAL MODELS
We are unable to solve the full non-linear disc evolution equation analytically and thus in general it is necessary to solve the equations numerically. We use the numerical scheme described by Pringle (1992), including both a binary torque precession term and a Lense-Thirring precession term. The binary torque is implemented in the form with τν 2 = R 2 /ν2, corresponding to just the first term in the expansion in equation (12). We have assumed that the warping is small so that ez.l ≈ 1 (Ogilvie & Dubus 2001). The Lense-Thirring term is We use 200 logarithmically distributed grid points between the inner radius Rin and the outer radius Rout = 10 4 Rin. We choose two different set of boundary conditions. In the first we have zero torque at the boundaries so that Σ = 0 at R = Rin and R = Rout which allows mass to flow freely through both boundaries. The surface density is chosen initially to correspond to a steady state for a flat accretion disc with the corresponding accretion rateṀ , (Pringle 1992). We add mass at R add = 9000 Rin at a rate that keeps the surface density constant at that point with angular momentum corresponding to the local Keplerian values and aligned with ez. We expect this to be in the outer flat part of the disc. Inside this radii the disc acts as a standard accretion disc. We set ∂l/∂R = 0 at R = Rin and R = Rout. These boundary conditions are realistic but to compare directly with our analytic models (Section 6) we use a second set of boundary conditions. We compute numerical models in which initially Σ =Ṁ /(3πν1) and the surface density is fixed at Rin and Rout. The disc is aligned with the black hole at the inside so at R = Rin, l = (sin θ0, 0, cos θ0) and at the outside R = Rout, l = (0, 0, 1).
We choose the viscosities ν1 = 0.1 (R/Rin) 2 and ν2 = 1.0 (R/Rin) 2 so that β = 2 and the local viscous timescale τν = R 2 /ν is independent of radius. This minimises computational run times. We have chosen ν1/ν2 = 0.1 which corresponds to α1 = 0.22 with equation (6). With SPH simulations of warped discs, Lodato & Pringle (2007) find that the ratio ν1/ν2 is small. Because ΩLT ∝ R −3 and Ω tid ∝ R 3/2 we may fix the various relevant radii by choosing the strengths of the tidal and Lense-Thirring terms. Models shown in the next sections are at a time of 500 (R 3 in /GM1) 1 2 = 50 τν 1 = 500 τν 2 so that they have reached steady state.

Misalignment of θ0 = 20 •
We first choose parameters so the Lense-Thirring radius is RLT = 10 Rin and the tidal radius is R tid = 100 Rin. We have 10 Rg < RLT < 60 Rg and so we can use the approximation to the Bardeen-Petterson effect in equation (3). For β = 2 the warp radius is then Rwarp = 21.5 Rin (equation 16). The outer disc and tidal precession is aligned with the z-axis. The Lense-Thirring precession axis is tilted at an angle of θ0 = 20 • to the z-axis. The tilt angle of the disc as a function of radius is shown in the left panel of Fig. 1. The relevant radii are also marked. Note that the disc itself is not only tilted but that the tilt varies with radius so that the overall shape of the disc is a twisted warp (right panel).
We see that the maximum rate of change of tilt does occur at about Rwarp. Once we are well within RLT or well outside R tid the disc does indeed become locally flat and aligned with the locally dominant precession axis. We see that both the numerical models with different boundary conditions and the analytical model are very similar for this small misalignment angle. In Fig. 2 we plot the surface density for the two numerical models in Fig. 1. The models are similar except in the regions close to the boundaries as expected.

Changing the Warp Radii
Here we demonstrate the effect of changing the outer two warp radii. We have three models each with the constant inclination and surface density boundary conditions. We choose RLT and R tid so that Rwarp = 21.5 remains unchanged. The first model has the same warp radii as that in Fig. 1 with RLT = 10 Rin and R tid = 100 Rin. The second has RLT = 3.16 Rin and R tid = 10 3 Rin and the third RLT = R tid = Rwarp = 21.5 Rin. We use Σ =Ṁ /(3πν1) at Rin and Rout.
In Fig. 3 we plot the inclinations and twists of these models. The steepest rate of change of tilt angle occurs at around Rwarp in all models. As R tid increases the outer parts of the disc become aligned with the outer precession vector (orbital plane) at larger radii. As RLT decreases the inner disc aligns with the black hole spin at smaller radii. So Rwarp tells us where the disc is warped RLT and R tid tell us how sharp the warp is.

Misalignment of θ0 = 60 •
In Fig. 4 we plot the inclination and twist of a disc which is misaligned by 60 • . The numerical model with the fixed inclination and surface density boundary conditions is similar to the analytic solution even at this high inclination. Figure 2. The surface density for numerical models with black hole spin misaligned to the binary orbital axis by θ 0 = 20 • . An accretion disc with both the Lense-Thirring effect and a binary torque term with R LT = 10 R in , R tid = 10 2 R in and Rout = 10 4 R in . All models have ν 1 , ν 2 ∝ R 2 . The solid line has the zero torque boundary conditions Σ = 0 at R in and Rout. The dotted line has Σ =Ṁ /(3πν 1 ) at R in and Rout as in our analytic models.

Misalignment of θ0 = 80 •
In Fig. 5 we plot the inclination and twist of a disc which is misaligned by θ0 = 80 • . We use only the boundary conditions Σ =Ṁ /(3πν1) at Rin and Rout. At this high inclination the non-linear terms have a bigger effect and the model with the same boundary conditions as the analytic solution is not such a good fit.

ANALYTICAL MODELS
To obtain a full solution to the warped disc problem it is necessary to resort to numerical means. However it is possible to obtain a good analytical approximation to the disc structure for θ 60 • . Here we construct such a solution and then investigate to what extent it provides a satisfactory fit to the full solution.
The full equation governing the evolution of disc tilt is non-linear but, as noted by Scheuer & Feiler (1996), when the disc tilt is small enough the non-linear term may be neglected. Once the equation has been reduced to a linear one, analytical techniques may easily be employed (see also Martin, Pringle & Tout 2007).
For small tilts to the z-axis the local unit vector l = (lx, ly, lz) that describes the direction of the angular momentum can be written to first order as l = (lx, ly, 1) with lx, ly ≪ 1. We then work in terms of the complex variable W = lx + ily.
(20) Figure 1. Models with the black hole spin misaligned to the binary orbital axis by θ 0 = 20 • . Left: The inclination of the disc. Right: The twist of the disc. The models are of an accretion disc with both the Lense-Thirring effect and a binary torque term with R LT = 10 R in , R tid = 10 2 R in and Rout = 10 4 R in . The solid lines have the boundary conditions Σ = 0 at R in and Rout. The dotted lines are a numerical model with Σ =Ṁ/(3πν 1 ) at R in and Rout which corresponds to the analytic boundary conditions. The dashed line is the matched analytical solution described in Section 6.3. All models have ν 1 , ν 2 ∝ R 2 . Figure 3. Numerical models of an accretion disc with both the Lense-Thirring effect and a binary torque term with the constant inclination and surface density boundary conditions. The black hole spin is misaligned to the binary orbital axis by 20 • and ν 1 and ν 2 ∝ R 2 . All three models have Rwarp = 21.5 R in . The solid lines have R LT = 3.16 R in and R tid = 10 3 R in . The dotted lines are the numerical models in Fig. 1 with R LT = 10 R in and R tid = 100 R in . The dashed lines have Rwarp = R LT = R tid = 21.5 R in . We use Σ =Ṁ /(3πν 1 ) at R in and Rout.  We take the scalar product of equation (1) with l and find This has the solution where Cacc and C dec are constants. The first term is the accretion disc term for which ν1Σ = const and the second term is the decretion disc term which has ν1Σ ∝ R − 1 2 (Pringle 1991).
For an accretion disc we set L = 0 at R = 0 and L = (GM1R) 1 2 Σ everywhere so C dec = 0. For a decretion disc we have Cacc = 0 and so more generally we find where x = 0 for an accretion disc and x = 1/2 for a decretion disc.

Solution
Martin et al (2007) and Scheuer & Feiler (1996) have demonstrated how to obtain analytic solutions for the linear warp equation with a power law density profile and a single added precession term. Here we wish to add two precession terms, the Lense-Thirring precession at small radii and the tidal precession at large radii. We have not been able to find analytic solutions for this case. However, as we have seen from the numerical results, at most radii in the disc one of the precession terms dominates. Thus we may obtain an approximate solution to the disc by matching two solutions. At small radii we use the solution of Martin et al (2007) for a disc with Lense-Thirring precession. At large radii we obtain a solution for a disc with only tidal precession. We then match these solutions at the intermediate radius Rwarp where the precession terms are of comparable but small magnitude.

Accretion Disc with Tidal Precession
In steady state without the non-linear term and with only the tidal torque term equation (1) becomes and, substituting our expression for L (equation 23) we find We add the x-component of equation (25) to i = √ −1 times the y-component to get W = lx + ily. Then where and the linearised tidal torque is up to O(W 2 ) with equation (12). Here R b is the binary separation and we have used lz ≈ 1 (Ogilvie & Dubus 2001). We only take the first term in the power series. This introduces an error which is relatively large in the outer parts of the disc. However, our outer boundary condition means that the torque itself at the outer edge of the disc is zero. We expect the disc to remain aligned with the binary orbit until R ≪ R b so the actual error remains small.
With the binary torque warp radius (equation 14) we simplify equation (26) governing the shape of the disc to where y = R/R tid . We solve this equation for accretion discs with x = 0 to find the shape of each disc.
For an accretion disc with x = 0, equation (29) becomes We make the substitution z = y We then set This is the modified Bessel equation with solution where I and K are modified Bessel functions of the first and second kind and A and B are, possibly complex, constants of integration. Substituting for our original variables we find We expect the outer edge of the disc to align with the binary orbit. Thus we expect W → 0 as R → R b in the frame of the orbit. The outer boundary condition should be dW/dR = 0 at R = Rout. However when Rout ≫ R tid then we can use W = 0 at R = Rout.

Matching the solutions
We expect the solution W tid to be valid at large disc radii. At small disc radii we make use of the solution given by Martin et al (2007) for a disc warped by Lense-Thirring precession.
The general solution in this case was found to bẽ (Martin, Pringle & Tout 2007) in the frame in which the spin of the black hole is in the z-direction r = R/RLT. Here again I and K are modified Bessel functions of the first and second kinds and C and D are, possibly complex, constants of integration. A disc subject to both Lense-Thirring and tidal torques is dominated by this shape when R ≪ Rwarp.
In order to match solutions we need first to transform WLT to the frame of the binary orbit, so that the axes of both solutions are aligned. Note that we still formally require all misalignment angles to be small. To transform this solution to the frame of the binary we obtain WLT = cos θ0ℜ(WLT)+sin θ0 where θ0 is the angle between the spin axis of the black hole and the normal to the orbital plane. We note that we are introducing non-linear terms by this transformation but because it is only a rotation of the axes it does not affect the shape. We now take these two solutions, the inner solution WLT and the outer solution W tid , and match them at the intermediate radius Rwarp, the radius at which the two precession terms have an equally small effect. In this way we find a matched analytical solution valid at all radii throughout whole disc. We are working in the frame of the binary. We expect the disc to be aligned with the black hole spin and so take This gives one relation between the constants of integration C and D. At the outer edge we expect the disc to be aligned with the binary orbit and so take This gives one relation between the constants of integration A and B.
At R = Rwarp we require the solution for W to be continuous so that Because W satisfies a second order differential equation we also require that the derivative of W be continuous so that These two conditions give two further relations between the constants A, B, C and D so we now have enough to solve for the complex constants A, B, C and D and find the complete matched solution for the disc. The form of the constants is complicated and does not add anything particular to our understanding so we don't reproduce them here. We do however construct many such matched solutions and compare them with our numerical models.

ACCURACY OF THE ANALYTIC SOLUTIONS
We plot these analytical solution in the dashed lines for the parameters of the numerical models in Fig. 1, 3, 4 and 5 alongside the numerical calculations. At an inclination θ0 = 20 • the solutions match very well. Even at θ = 60 • they match well when the same boundary conditions are used. Though deviations are larger at 80 • much of the shape is still reproduced. We note that the analytical solutions are easily found for any β. The numerical solution takes much longer to evaluate when β is small and the timestep, controlled by the centre of the disc, becomes very small.

Error in neglecting the non-linear terms
We have assumed that the warp is gradual enough that we can neglect the non-linear term l.∂ 2 l/∂R 2 = − |∂l/∂R| 2 . The relative error in neglecting this term is given by (41) Figure 6. The effect of neglecting the non-linear terms in the matched analytical solution for a disc with the spin axis of the black hole misaligned with the binary orbital axis by θ 0 = 20 • . We choose R LT = 10 R in , R tid = 100 R in and ν 2 /ν 1 = 10 as in Fig. 1.
We can justify neglecting this term if E1 ≪ 1. With similar power-law viscosities We plot this for the matched solution for a misalignment of θ0 = 20 • in Fig. 6. We see that the error is at a maximum of just over 1 per cent at R = RLT. In Fig. 7 we plot the maximum value of E1 for different misalignment angles between the black hole spin and the binary orbit. It begins to dominate after θ0 = 60 • and this accounts when most of the error in the model for θ = 80 • .

SHAPE OF THE DISC
Given a position in cylindrical polar coordinates (R, φi) in the x − y plane of the binary orbit, we can find the height, z, of the disc at that point. In Cartesian coordinates the disc can be described by where θ(R) = cos −1 (lz) is the angle of the direction vector of the disc, l, at radius R, to the z-axis. The azimuthal angle of l is φ = tan −1 (ly/lx). In Fig. 8 we picture this disc solution, with a logarithmic scale in the radial direction and Rin mapped to the origin, using the plotting package Povray. We can see the warp and twist of the disc for a disc with a misalignment of θ0 = 45 • with the warp radii RLT = 10 Rin and R tid = 100 Rin. These models can be used to investigate disc properties Figure 7. The maximum error in neglecting the non-linear terms for a model with R LT = 10 R in and R tid = 100 R in and with ν 2 /ν 1 = 10 for a range of misalignments θ 0 between the black hole spin and the orbital angular momentum of the binary system. and compare with observations. We illustrate this by examining the obscuration of the central black hole by the disc at various inclinations. In Fig. 9 we plot the visibility of the central black hole on the sky for θ0 = 20 • and 60 • for the numerical model with zero torque boundary conditions and for the analytically matched solution. The binary orbital axis is inclined at i to the line of sight. The angle φ is the phase around the disc. It changes on the precession timescale of the black hole (Martin, Tout & Pringle 2008). This is long compared with the orbital period so it is essentially fixed relative to the observer. As expected, the higher the misalignment, the more obscured the black hole is. In each plot the region between the two lines would be obscured. Known systems would therefore tend to be observationally selected to be aligned.

CONCLUSIONS
We have made both numerical and analytical models of accretion discs warped by the Lense-Thirring effect and a binary companion. The models have ν ∝ R β with β = 2 so that the viscous timescale is constant throughout the disc. This means that the models converge at their fastest rate. More realistically β = 3/4 and the numerical disc evolves slowly because the small viscous timescale at the inner edge dominates. On the other hand our analytical matched models can be constructed just as quickly for any β. Up to a misalignment of about θ = 60 • they provide a very good fit to the disc structure and so may be used where many calculations are required in a short time. This is necessary when attempting a fit to observations with unknown misalignment. With a relatively simple analytical fit it is easy to compute all sorts of disc properties such as the extent of irradiation by the central engine at any point or the integrated spectrum without the need to compute a numerical model for every set of parameters. We have found that known systems are observationally selected to be aligned. Figure 9. The visibility of the black hole with a disc with the numerical (solid lines) and matched analytical disc solutions (dashed lines) with R LT = 10 R in , R tid = 100 R in and a misalignment angle of θ 0 = 20 • (left) and 60 • (right). The binary orbit is inclined at i to the line of sight and φ is the phase around the disc. The black hole itself and the central parts of the disc are obscured between the two lines (the shaded region).